# Kinetic theory of stellar systems: A tutorial

Chris Hamilton, Jean-Baptiste Fouvry

► To cite this version:

Chris Hamilton, Jean-Baptiste Fouvry. Kinetic theory of stellar systems: A tutorial. Physics of Plasmas, 2024, 31 (12), 10.1063/5.0204214. hal-05343133

HAL Id: hal-05343133

https://hal.science/hal-05343133v1

Submitted on 4 Nov 2025

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
# Kinetic theory of stellar systems: A tutorial

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

Special Collection: **2024 Early Career Collection**, **Reviews and Tutorials in Basic Plasma Phenomena, Waves, and Instabilities**, **Papers from the 65th Annual Meeting of the APS Division of Plasma Physics**

Chris Hamilton

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

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

; Jean-Baptiste Fouvry

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

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

Check for updates

Phys. Plasmas 31, 120901 (2024)

https://doi.org/10.1063/5.0204214

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

View
Online

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

Export
Citation

## Articles You May Be Interested In

### Low temperature plasma biomedicine: A tutorial review

*Phys. Plasmas* (August 2014)

### Hamiltonian description for magnetic field lines in fusion plasmas: A tutorial

*Phys. Plasmas* (September 2023)

### Whistlers, helicons, and lower hybrid waves: The physics of radio frequency wave propagation and absorption for current drive via Landau damping

*Phys. Plasmas* (September 2015)

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

AIP
Publishing

Physics of Plasmas

Special Topics Open
for Submissions

Learn More

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

AIP
Publishing
# Kinetic theory of stellar systems: A tutorial EP

Cite as: Phys. Plasmas **31**, 120901 (2024); doi: 10.1063/5.0204214

Submitted: 20 February 2024 · Accepted: 17 October 2024 ·

Published Online: 6 December 2024

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

View Online

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

Export Citation

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

CrossMark

Chris Hamilton<fnref n="1" />,ª ID and Jean-Baptiste Fouvy<fnref n="2" /> ID

## AFFILIATIONS

**Note:** This paper is part of the Special Collection: Papers from the 65th Annual Meeting of the APS Division of Plasma Physics.

**Note:** Paper PT2 1, Bull. Am. Phys. Soc. **68** (2023).

## ABSTRACT

Stellar systems—globular and nuclear star clusters, elliptical and spiral galaxies and their surrounding dark matter haloes, and so on—are ubiquitous characters in the evolutionary tale of our Universe. This tutorial article is an introduction to the collective dynamical evolution of the very large numbers of stars and/or other self-gravitating objects that comprise such systems, i.e., their *kinetic theory*. We begin by introducing the basic phenomenology of stellar systems, and explaining why and when we must develop a kinetic theory that transcends the traditional two-body relaxation picture of Chandrasekhar. We then study the individual orbits that comprise stellar systems, how those orbits are modified by linear and nonlinear perturbations, how a system responds self-consistently to fluctuations in its own gravitational potential, and how one can predict the long-term evolutionary fate of a stellar system in both quasilinear and nonlinear regimes. Though our treatment is necessarily mathematical, we develop the formalism only to the extent that it facilitates real calculations. Each section is bolstered with intuitive illustrations, and we give many examples throughout the text of the equations being applied to topics of major astrophysical importance, such as radial migration, spiral instabilities, and dynamical friction on galactic bars. Furthermore, in the 1960s and 1970s, the kinetic theory of stellar systems was a fledgling subject which developed in tandem with the kinetic theory of *plasmas*. However, the two fields have long since diverged as their practitioners have focused on ever more specialized and technical issues. This tendency, coupled with the famous obscurity of astronomical jargon, means that today relatively few plasma physicists are aware that their knowledge is directly applicable in the beautiful arena of galaxy evolution, and relatively few galactic astronomers know of the plasma-theoretic foundations upon which a portion of their subject is built. Yet, once one has become fluent in both Plasmaish and Galacticese, and has a dictionary relating the two, one can pull ideas directly from one field to solve a problem in the other. Therefore, another aim of this tutorial article is to provide our plasma colleagues with a jargon-light understanding of the key properties of stellar systems, to offer them the theoretical minimum necessary to engage with the modern stellar dynamics literature, to point out the many direct analogies between stellar- and plasma-kinetic calculations, and ultimately to convince them that stellar dynamics and plasma kinetics are, in a deep, beautiful and *useful* sense, the same thing.

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0204214

## I. INTRODUCTION

The evolution of our Universe is a great unfolding drama in which galaxies play a fundamental role. These galaxies grow and die; they pulse and throb; they swallow gas and they belch it back out again; and they merge together and they tear one another apart. The study of the collective dynamical behavior that occurs inside galaxies—and of the violent lives they lead, the tumultuous dark matter environments in which they sit, and the star clusters they cradle—is best termed the *kinetic theory of stellar systems*. This is the subject with which we are concerned in this tutorial article.

The kinetic theory of stellar systems rests upon two main theoretical pillars, and our tutorial article aims to introduce colleagues to both of them. The first pillar is the theory of *orbits*. By this, we mean the orbital motion of stars and/or dark matter particles in gravitational fields, or, more broadly, the dynamics of point particles subject to some Hamiltonian flow. The basic language of this topic is angle-action coordinates. These are not only the natural variables with which to describe smooth orbits in quiescent gravitational fields but also to understand how these orbits are modified by *perturbations* of various kinds. These perturbations tend to have the most dramatic effect on

<footnote n="1">¹Institute for Advanced Study, Einstein Drive, Princeton, New Jersey 08540, USA</footnote>

<footnote n="2">²CNRS and Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote>ª) Invited speaker. **Author to whom correspondence should be addressed:** chamilton@ias.edu</footnote>
orbits near resonances, where the frequency of the perturbation matches the natural frequency of the orbit. The first aim of our tutorial article is to introduce colleagues to these concepts. We employ as little heavy formalism as possible, assuming that readers already have a passing familiarity with Hamiltonian mechanics, our aim being to help them understand how they apply in the stellar dynamical context.

The second pillar upon which our subject is built is the collective dynamics of many bodies, and here our theoretical foundation is plasma kinetics. Why is this the case in a treatise on gravity? Well, every high school student knows that electromagnetism and gravity are intellectually adjacent subjects: up to a sign, the Coulombic electrical force and the Newtonian gravitational force are equivalent. What is less well appreciated is that this adjacency extends much further than simple $\pm 1/r^2$ force laws. Indeed, the mathematical methods developed to probe the collective behavior of protons and electrons in astrophysical and fusion plasmas are in many cases precisely the same as those needed to describe the dynamics of stars, spiral arms, and dark matter in galaxies like the Milky Way. Concepts like two-body relaxation, linear response theory, Landau damping, resonant wave-particle interactions, and phase mixing all had origins in plasma theory but are now part of the standard lexicon of galactic dynamics. However, these two fields long ago went their separate ways, and in much of today's literature are separated by a near-impenetrable wall of jargon. In the second part of this tutorial article, we develop the tools needed to describe collective behavior in galaxies. We aim to provide colleagues with the theoretical minimum needed to appreciate galactic-dynamical problems. In so doing, we aim to convey—to those on both the stellar-dynamical and plasma-theoretic sides of the aisle—just some of the deep connections between the two fields, and to encourage much greater crosstalk between these two estranged, yet intimately related, subjects.

In contrast with most other treatises on stellar (or “galactic”) dynamics and galaxy evolution, we approach our subject from the point of view of theoretical physics rather than practical astronomy. Thus, we spill almost no ink over such crucial concepts as mass-to-light ratios and chemical evolution and dust extinction (though we happily refer the reader to the many excellent texts that do).¹⁻³ Our models are necessarily rather simple, in some cases even crude, since our goal here is to provide physical insight rather than to accurately

reproduce observational data. On the other hand, we do not simply develop mathematical formalism for its own sake. Almost every formal calculation we present is accompanied by a concrete example of its application. We hope this approach underscores our key message, which is that the kinetic theory of stellar systems is not merely the mathematical fancy of a few specialists. Instead, it can and should aid our understanding of galaxy evolution; and it passes the basic test of all good physics, by telling us things about the real world that we did not already know.

## A. A (very) brief astronomy lesson

To set the context for this tutorial article, and to get an idea of the characteristic scales involved, we invite you to consider Fig. 1. Panel (a) shows an image of a *disk galaxy* called the Southern Pinwheel, captured by the la Silla Observatory in the mountains of Chile. We have chosen this image not because it is visually stunning—though who could deny that?—but because it exhibits several of the paradigmatic features of disk galaxies, and is in fact thought to resemble our own Galaxy, the Milky Way. Like the Milky Way, the Southern Pinwheel has an elongated or *barred* structure at its center and exhibits beautiful *spiral arms*. The red knots strung along the arms are $H\alpha$ emission from regions of recent star formation, meaning stars have been born there in the last few tens of Myr (where $1 \text{ Myr} \equiv 10^6$ years), while the bluer light comes from much older stars, up to ages of $\sim 10$ Gyr (where $1 \text{ Gyr} \equiv 10^9$ years; the Big Bang happened 13.7 Gyr ago, and galaxies formed in the first few Gyr that followed).³ The image is hap-hazardly obscured by thin, dark *dust lanes*, regions wherein visible starlight has been absorbed. The whole galaxy is a few tens of kiloparsecs (kpc) across, where $1 \text{ pc} \sim 3 \times 10^{16} \text{ m}$. (The distance from the Sun to the next-closest star, Proxima Centauri, is of order $\sim 1$ pc.) From our perspective, we are seeing the Southern Pinwheel galaxy more or less face-on, but if we were to see it from the side we would find that it was very thin, with a central bulge—like two fried eggs glued back to back. It consists of around $N \sim 10^{11}$ stars, most of which are orbiting in this very thin, nearly two-dimensional disk plane, and they mostly move on (quasi-)circular orbits around the center of the galaxy, all rotating in the same direction. This means that the motion of stars in the galaxy is highly ordered, a.k.a. *dynamically cold*. Another

<figure><img src="image_14.png" /><figcaption>FIG. 1. (a) The spiral galaxy M83 (the Southern Pinwheel). Credit: ESO. (b) The globular cluster M80. Credit: NASA, The Hubble Heritage Team, STScI, AURA.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIG. 2. Schematic illustration of typical orbits in (a) cold disks and (b) hot spheres (idealized versions of the archetypical systems shown in Fig. 1). In a cold disk, most orbits are nearly circular and confined to low inclinations relative to the disk plane. By contrast, in a hot sphere the orbits are near-isotropically distributed and have a broad range of eccentricities.</figcaption></figure>

way to say this is that the stellar disk, like a plasma beam, has low-velocity dispersion. The key dynamical timescale in this disk is the azimuthal orbital period of one of the stars, which is typically $t_{dyn} \approx 100$ Myr, and this timescale increases with orbital radius. Since galaxies are ~10 Gyr old, we deduce that most stars have completed fewer than 100 orbits in the lifetime of the galaxy.

Now let us contrast this with panel (b) of Fig. 1, which is an image of the globular cluster M80. This is another optical image, taken with the Hubble Space Telescope. Globular clusters are much smaller than galaxies, with typical radii of ~10 pc—in fact, a galaxy like the Southern Pinwheel typically contains many hundreds or even thousands of globular clusters. Despite their relative smallness, globular clusters have a much higher density of stars than do galaxies. The cluster M80 contains something like $N \sim 10^5$ stars, so the typical distance between stars is $\sim(10^5/\text{pc}^3)^{-1/3} \sim 0.02$ pc, of order 50 times smaller than the distance from the Sun to Proxima Centauri. Globular clusters typically contain very little gas and hence exhibit almost no ongoing star formation. This is why there is a much narrower range of colors in this image compared to the one beside it. Moreover, unlike disk galaxies, globular clusters are (nearly) spherically symmetric, so the image of M80 would not look very different if we were able to view it from another angle. Also unlike disk galaxies, they are disordered, dynamically hot systems, meaning that the motion of their stars has no strongly preferred direction (like the isotropic Maxwellian in an electrostatic plasma). More mathematically, at any given point $\mathbf{r}$ in the cluster, the local phase space distribution function $f(\mathbf{r}, \mathbf{v})$ depends primarily on the magnitude, rather than the direction, of the velocity $\mathbf{v}$. Though they are of a comparable age to galaxies—having formed ~10 Gyr ago—globular clusters are much more “dynamically old” than galaxies, since the typical orbital period of their stars is only $t_{dyn} \sim 10^5$ yr. Thus, each star in a globular cluster has undergone $\sim 10^5$ orbits in its lifetime.

To summarize: disk galaxies are dynamically young and cold, globular clusters are dynamically old and hot.

## B. Two mental models

The brief astronomy lesson above motivates us to introduce two simple mental models which we will refer to throughout this tutorial article: the cold, rotating axisymmetric disk and the hot, nonrotating spherical cluster—see Fig. 2 for an illustration.

These are really just (over-)simplified versions of disk galaxies and globular clusters, respectively, and constitute clean examples of stellar systems that exhibit completely different dynamical behaviors. Some of the characteristic properties of these two prototypical systems are summarized in Table I.

We will find that the cold disk is analogous to a nonthermal plasma (with collective oscillations/instabilities and wave-particle interactions playing a key role), while the hot sphere is more analogous to a thermal plasma or a classical gas (with a kinetic evolution primarily driven by discrete, local two-body encounters). However, there are many subtleties and differences to be encountered along the way, and gravitational systems are capable of behaviors that have no analog in gas or plasma kinetics.

It is also worth mentioning that throughout this article, for concreteness, we often refer to the motion of “stars.” However, we are also very much interested in the motion of dark matter particles, which obey the same set of equations as stars, and are just as important a component of galactic dynamics.

## C. Plan for the rest of this article

The rest of this tutorial article is organized as follows. In Sec. II, we introduce some of the basic properties of N-body self-gravitating systems, discuss the extent to which self-gravitating systems play by the same thermodynamic rules as gases and plasmas, and give a qualitative idea of how these systems evolve on long timescales. Then, we begin with the theory of orbits. In Sec. III, we approximate a galaxy’s gravitational potential by a smooth “mean field” potential and study the orbits of individual stars in this mean field, with the help of angle-

<table><caption>TABLE I. Characteristic scales for the two prototypical systems considered in this article, the disk galaxy and the globular cluster.</caption><thead><tr><th></th><th>Disk galaxy</th><th>Globular cluster</th></tr></thead><tbody><tr><td>Mental model</td><td>Cold, rotating<br>axisymmetric disk</td><td>Hot, nonrotating<br>sphere</td></tr><tr><td>Size</td><td>10s of kpc</td><td>10 pc</td></tr><tr><td>Age</td><td>10 Gyr</td><td>10 Gyr</td></tr><tr><td>Number of stars <i>N</i></td><td>10<sup>11</sup></td><td>10<sup>5</sup></td></tr><tr><td>Orbital period <i>t</i><sub>dyn</sub></td><td>100 Myr</td><td>0.1 Myr</td></tr></tbody></table>
action variables. In Sec. IV, we show how individual orbits respond to both linear and nonlinear perturbations and describe the theory of resonant trapping. Our proper study of kinetic theory begins in Sec. V, where we introduce the Klimontovich equation for an N-body system's phase space distribution function (DF) and derive a generic set of equations governing its evolution. In Sec. VI, we study the self-consistent linear response of stellar systems to weak perturbations, and then, in Sec. VII, we use those linear solutions to construct a quasilinear theory of statistically averaged evolution, culminating in the Balescu–Lenard equation. Our final topic, in Sec. VIII, is nonlinear kinetic theory, particularly in the presence of orbit trapping and with applications to bar-halo friction and the saturation of spiral instabilities.

Throughout these sections, we continually draw direct analogies between gravitational and plasma systems, as well as giving many concrete examples of the equations working in practice. Finally, in Sec. IX, we touch on some other areas of galactic dynamics that have plasma-theoretic analogues, and give references for further reading.

We mention here that the standard reference in this subject is the textbook Galactic Dynamics⁵ by Binney and Tremaine, which we will hereafter refer to as BT08. While our article is far narrower in scope than BT08, much of what we discuss cannot be found there or in any other textbook or publication. The main exceptions to this rule are the introductory sections Secs. II and III, where by necessity there is significant overlap with BT08. Where helpful, we point the reader to specific sections of BT08 in which more details can be found.

Finally, this tutorial grew out of notes that we both compiled as lecturers of the Kinetic Theory of Stellar Systems course at Oxford. Lecture materials, problem sets, etc., can be found online.⁵

## II. BASIC PROPERTIES OF N-BODY SELF-GRAVITATING SYSTEMS

The simplest theoretical approach to describing a stellar system is to imagine it as a gas of stars, and to apply the traditional methods of equilibrium thermodynamics and two-body collision theory. In this section, we will attempt to do just that, and we will find that while these classical notions can be useful in characterizing stellar systems, they also break down completely in some circumstances—and that in the strict sense, in self-gravitating systems there is no such thing as equilibrium. The insights we gain from this failure will point us toward the more complete kinetic description that we will develop in later sections.

### A. Virial equilibrium

The first thing to mention is that the speeds of stars in galaxies and star clusters are nearly always of order tens or hundreds of km s⁻¹, much less than the speed of light. This means that relativistic corrections to Newtonian gravity can almost always be ignored. There are exceptions, particularly for stars orbiting around supermassive black holes in galactic centers,⁶ but for our purpose, Newtonian gravity is an excellent approximation. Also, while stars do exhibit a spectrum of different masses, for simplicity we will always assume that our system consists of a single mass species, unless we explicitly state otherwise.

We therefore consider a system of N ≫ 1 stars of equal mass *m* (usually, but not necessarily, taken to be of order a solar mass *M*<sub>☉</sub>), interacting only via Newtonian gravity. The gravitational potential energy due to the interaction between two stars with positions **r** and **r'** is

$$U(\mathbf{r}, \mathbf{r}') = -\frac{Gm^2}{|\mathbf{r} - \mathbf{r}'|}, \quad (1)$$

where $G = 4.3 \times 10^{-3}$ pc $M_{☉}^{-1}$ km s⁻¹ is Newton's gravitational constant (here $M_{☉} = 2.0 \times 10^{30}$ kg is the mass of the Sun). Newton's second law then tells us that the force on particle *i* is

$$m \frac{d\mathbf{v}_i}{dt} = -\frac{\partial}{\partial \mathbf{r}_i} \sum_{j \neq i}^{N} U(\mathbf{r}_i, \mathbf{r}_j) = -Gm^2 \sum_{j \neq i}^{N} \frac{\mathbf{r}_i - \mathbf{r}_j}{|\mathbf{r}_i - \mathbf{r}_j|^3}, \quad (2)$$

where $\mathbf{v}_i \equiv \mathbf{dr}_i/dt$ is the *i*th star's velocity.

Gravity is an exclusively attractive force—two masses always attract and never repel—which means that unlike in plasma, there is no shielding, nor is there any sensible notion of a “Debye length” over which forces are screened. Since *N* is large (~10¹¹ in the cold disk, ~10⁵ in the hot sphere), and since every star is attracting every other star, what stops the entire system from collapsing under its own weight and forming a huge black hole? The answer depends on the context. In the hot sphere, the random motions of the stars counterbalance the gravitational attraction (and can be thought of as an effective “pressure”). In the case of a cold disk, it is the rotational support (i.e., the fact that the disk is spinning) that counterbalances gravity.

To put these notions on a more mathematical basis, we now derive the virial theorem. For an isolated N-body system, we can choose the origin of our coordinates to be the fixed center of mass. Then, the moment of inertia of the N-body system with respect to the origin is

$$I = \sum_{i=1}^{N} m \mathbf{r}_i^2. \quad (3)$$

Let us now take the second derivative of this quantity with respect to time

$$\frac{d^2 I}{dt^2} = 2 \left( \sum_{i=1}^{N} m \mathbf{v}_i^2 + \sum_{i=1}^{N} m \mathbf{r}_i \cdot \frac{d\mathbf{v}_i}{dt} \right). \quad (4)$$

Using the equation of motion (2), we can write the second term on the right-hand side as

$$\begin{aligned} \sum_{i=1}^{N} m \mathbf{r}_i \cdot \frac{d\mathbf{v}_i}{dt} &= -Gm^2 \sum_{i=1}^{N} \sum_{j \neq i}^{N} \frac{\mathbf{r}_i - \mathbf{r}_j}{|\mathbf{r}_i - \mathbf{r}_j|^3} \\ &= -\frac{1}{2} Gm^2 \sum_{i=1}^{N} \sum_{j \neq i}^{N} \frac{(\mathbf{r}_i - \mathbf{r}_j) \cdot (\mathbf{r}_i - \mathbf{r}_j)}{|\mathbf{r}_i - \mathbf{r}_j|^3} \\ &= \sum_{i=1}^{N} \sum_{j>i}^{N} U(\mathbf{r}_i, \mathbf{r}_j). \end{aligned} \quad (5)$$

For an equilibrium system, the moment of inertia can fluctuate around some average value taken over a few dynamical times, but that average value must not itself change. Denoting such a time average by $\langle \cdot \rangle$, we have that for a system in virial equilibrium

$$2K + W = 0, \quad (6)$$

where
$$K = \left\langle \frac{1}{2} m \sum_{i=1}^{N} \mathbf{v}_i^2 \right\rangle, \quad W = \left\langle \sum_{i<j}^{N} U(\mathbf{r}_i, \mathbf{r}_j) \right\rangle \quad (7)$$

are the averaged total kinetic and total potential energies of the system, respectively. In other words, in virial equilibrium, the kinetic and potential parts of the energy of the system balance each other to within an order-unity factor.

A globular cluster (right panel of Fig. 1) is an excellent example of a system in virial equilibrium. Its smooth, near-spherical shape reflects the fact that it is not, on average, undergoing any violent bulk motions, despite the stars inside it whizzing around like a swarm of bees on the timescale $t_{dyn}$. Note that even though a cluster is in virial equilibrium, that equilibrium *does* still evolve; but the evolution of this “quasi-stationary state” occurs over a timescale much longer than $t_{dyn}$, as we will see (Sec. II C).

The virial theorem allows us to make quick order-of-magnitude estimates about typical random velocities of stars, $\sigma$. For instance, consider that a globular cluster has a mass $M \approx Nm \sim 10^5 M_\odot$ and a radius $R$ of a few parsecs. Using the virial theorem, we can calculate the typical velocity dispersion of stars in this cluster. Roughly, we have $K \approx M\sigma^2/2$ and $W \approx -GM^2/R$, so using (6), we get

$$\sigma \sim \sqrt{\frac{GM}{R}} \approx 10 \text{ km s}^{-1} \left( \frac{M}{10^5 M_\odot} \right)^{1/2} \left( \frac{R}{4\text{pc}} \right)^{-1/2}, \quad (8)$$

in agreement with what is observed (and much less than $c$, justifying our use of Newtonian gravity).

Finally, let us contrast globular clusters with an example of a system, which is very much *not* in virial equilibrium. Figure 3 shows the Antennae system, which consists of two galaxies that have collided and are in the process of merging.

For this system, Eq. (6) does not hold; instead, energy is rapidly being transferred between kinetic and potential forms on the timescale $t_{dyn}$ as the two galaxies perform their cosmic dance. However, we know from observations and numerical simulations of such galaxies

<figure><img src="image_16.png" /><figcaption>FIG. 3. The Antennae galaxies NGC4038 and NGC4039. These two colliding galaxies are in the process of merging into a single structure. It is easy to see by eye that this system is going to undergo significant morphological changes on the typical timescale $\sim t_{dyn}$. Hence, it is not in virial equilibrium. Credit: ESA/Hubble.</figcaption></figure>

that they ultimately *do* virialize and that this virialization is nearly complete after only a few $t_{dyn}$ (which in this case still means a few $10^9$ years or more). This rather rapid evolution toward a virialized state is called *violent relaxation*.⁷ The precise mechanics of violent relaxation, and the question of which particular virial equilibrium state the system will fall into, are a matter of ongoing research not only in stellar dynamics but also in plasma physics.⁸

For the rest of this article, we will only consider virialized stellar systems whose bulk properties are evolving on timescales $\gg t_{dyn}$. Isolated disk galaxies and globular clusters both fall into this category.

## B. Thermal equilibrium?

We mentioned above that virial equilibria are not true equilibria of stellar systems—rather, they are just quasi-equilibria that evolve on timescales much longer than the dynamical time $t_{dyn}$. Yet why is this the case? After all, this is not what happens in either a gas or an electrostatic plasma: there, the natural thermal equilibrium distribution of particles’ positions and velocities is a homogeneous Maxwellian. Unfortunately, gravity does not play by those thermodynamic rules.

To see why, suppose we wished to describe a globular cluster as a “stellar gas” in thermal equilibrium. We could make a natural definition of the cluster’s temperature $T$ by relating it to the total kinetic energy

$$K = \frac{3}{2} N k_B T, \quad (9)$$

where $N$ is the number of stars. Then, the total energy of the cluster $E_{tot}$ is the sum of the kinetic and potential energies

$$E_{tot} = K + W. \quad (10)$$

However, from the virial theorem (6), we know that $W = -2K$, so $E_{tot} = -K < 0$. The fact that $E_{tot}$ is negative is not unusual: it just means that the system is gravitationally bound. The heat capacity is then

$$C = \frac{\partial E_{tot}}{\partial T} = -\frac{\partial K}{\partial T} = -\frac{3}{2} N k_B < 0. \quad (11)$$

Notice that $C$ is *negative*. This is strange: if you add heat to the system (increase its total energy $E_{tot}$), it cools down! In other words, the more energy is given to the cluster, the lower the amplitude of the random velocities inside it.

The negativity of the heat capacity implies that it is impossible for a globular cluster ever to reach a thermal equilibrium. For instance, consider a cluster at temperature $T$ immersed in a heat bath of the same temperature. Suppose there is some infinitesimal negative fluctuation in the cluster’s energy, $-|\delta E|$. By Eq. (11), this will cause the cluster to achieve a higher temperature by $\delta T = (-|\delta E|)/C > 0$, i.e., with higher velocity dispersion. Since the temperature is now higher, the system loses more heat to the bath, which then makes its temperature increase even further, and so on without limit.

Nor can one achieve thermal equilibrium in an isolated stellar system (i.e., in the absence of any heat bath). The reason is that any star with a sufficiently large velocity is unbound from the system and hence will escape to infinity, a process known as *evaporation*. This precludes the system from ever attaining a Maxwellian velocity
distribution function (DF), since the cluster would otherwise have a finite phase space density out to $v \to \infty$.

One might argue at this point that the Maxwellian is only an idealization, and not being able to *exactly* realize a Maxwellian DF is a rather pedantic argument for the nonexistence of thermal equilibrium. (Indeed, all speeds are ultimately bounded by the speed of light, but this does not stop us from employing a Maxwellian DF when describing plasmas!) To justify our claim further, let us calculate the rate at which stars *would* evaporate from a hypothetical Maxwellian stellar system. As we will argue in Sec. III, the motion of individual stars on the timescale $\sim t_{dyn}$ is very accurately described by a *mean field approximation*, in which we ignore the complexities of the N-body system and instead consider the stars to be moving in a *smooth gravitational potential* $\Phi(\mathbf{r})$, in precisely the same way as one imagines the unperturbed motion in a homogeneous electrostatic plasma to be a straight line. Then, the energy per unit mass of any star with velocity $v$ and position $\mathbf{r}$ is $v^2/2 + \Phi(\mathbf{r})$, meaning that any star with $v > v_{esc}(\mathbf{r}) \equiv \sqrt{-2\Phi(\mathbf{r})}$ will evaporate. Averaging this over the entire system gives us the typical escape speed

$$ \sigma_{esc}^2 \approx \frac{1}{M} \int d\mathbf{r} \rho(\mathbf{r}) v_{esc}^2(\mathbf{r}) = -\frac{2}{M} \int d\mathbf{r} \rho(\mathbf{r}) \Phi(\mathbf{r}), \quad (12) $$

where $\rho_0(\mathbf{r})$ is the mass density of stars at position $\mathbf{r}$ and $M$ is the total cluster mass. This is obviously related to the total potential energy of the system W. Being careful not to over-count the pairwise interactions, and using the virial theorem (6), we get

$$ \sigma_{esc}^2 \approx -\frac{4W}{M} = \frac{8K}{M} = 4\sigma^2, \quad (13) $$

where $\sigma^2 = N^{-1} \langle \sum_{i=1}^N v_i^2 \rangle$ is the cluster's velocity dispersion. Supposing the cluster did have a Maxwellian DF with the same velocity dispersion, we can estimate the fraction of unbound stars by

$$ f_{esc} \approx \frac{\int_{2\sigma}^{\infty} dv \, v^2 e^{-3v^2/2\sigma^2}}{\int_{0}^{\infty} dv \, v^2 e^{-3v^2/2\sigma^2}} \approx \frac{1}{140}. \quad (14) $$

In other words, we would expect about 1% of an isolated stellar system to evaporate away on the timescale of refilling of the Maxwellian tail, which is roughly a relaxation time (see Sec. II C for the definition of this term). The equivalent situation in a collisional plasma would be that 1% of the electrons disappeared every mean-free time! Thus, even if the cluster “wanted” to reach a Maxwellian steady-state, it could not do so because stars are constantly evaporating from the high-velocity tail.

Evaporation is a real and important effect in cluster dynamics; indeed, there are certain star clusters called *open clusters* with $N \sim 10^3$ (some of which you can see with the naked eye, like the Pleiades) whose lifetime is limited by evaporation to a few hundred Myr. In other words, there were star clusters visible to the dinosaurs that have since evaporated into the vastness of space and are no longer visible to us. Yet at no point was their DF ever Maxwellian.

With minor modifications, the above arguments also hold for cold, differentially rotating disks, as well as other virialized self-gravitating structures. The bottom line is that *strict thermal equilibrium is impossible in a self-gravitating system*, and there is no universal requirement for stellar systems to “tend toward the Maxwellian” as

gases and plasmas do. Fundamentally, this stems from the purely attractive nature of gravity: since the potential energy of a very loosely bound particle decays like $1/r$, there are infinitely many ways to rearrange the phase space coordinates of very distant particles with no energy cost. In the language of statistical mechanics, this means that in isolated self-gravitating systems there is no microcanonical ensemble, no canonical ensemble, energy is not extensive, and entropy is not bounded (so there exists no maximum entropy state$^{1,9}$). Concepts like entropy are still useful in stellar dynamics, particularly if the phase space is bounded by, e.g., tidal truncation (meaning the system is not truly isolated) or by the approximate conservation of certain invariants (meaning there is a quasi-maximum entropy state on short enough timescales).$^{10}$ However, in general, one has to be very careful about naive application of “standard” statistical mechanical results to gravitational systems.

### C. So what does happen to a virialized stellar system?

We have argued that there exists no unevolving thermal equilibrium state for stellar systems. What, then, does happen to stellar systems after they reach virialization? Roughly speaking, the answer is that they move slowly from one virialized quasi-equilibrium state to another on a timescale that we call the *relaxation time*, $t_{relax}$.

How do we calculate the relaxation time? Roughly speaking, it is the timescale over which a typical star “forgets” its initial orbit. A proper theory of relaxation requires an understanding of stellar orbits (Sec. III), but for now the intuition from gas or plasma kinetics is good enough. Hence, we simply ask, supposing an unperturbed star was moving in a straight line, how long would it take for its velocity $v$ to de-correlate from its initial value, i.e., to accumulate an impulse $|\Delta v| \approx |v|$?

In a gas, these impulses occur because a given molecule collides with some other molecule roughly once per mean free time. Hence, after a few such collisions, this molecule has completely forgotten its initial state [see Fig. 4(a)].

This is not so in a stellar system. Even in a dense system like a globular cluster, the distance between single stars is such that direct physical collisions are very rare.$^{11}$ Instead, in a cluster, the evolution of a star's orbit is mostly driven by a very large number of weak nudges from all the other stars it is interacting with simultaneously [Fig. 4(b)]. In the coming subsections, we will follow Chandrasekhar and estimate the relaxation timescale due to a series of such nudges (Sec. II C 1), a calculation completely analogous to the Spitzer two-body collision calculation from plasma kinetics. We will also pause to comment on the confusing words “collisional” and “collisionless” as applied to stellar systems (Sec. II C 2). Then (Sec. II C 3), we will show that Chandrasekhar's theory is very successful at describing the evolution of hot spheres (in the same way that Spitzer's theory works for describing collisional plasmas) but fails completely to describe the evolution of cold disks (as does Spitzer's theory when dealing with collisionless plasma out of equilibrium). This will be the jumping-off point where stellar dynamics ceases to consist of a series of rough estimates; thereafter, things get quantitative.

#### 1. Chandrasekhar's theory of two-body relaxation

Here, we give the simplest version of Chandrasekhar's calculation of the relaxation time,$^{12}$ based on his theory of two-body
<figure><img src="image_17.png" /><figcaption>FIG. 4. (a) Relaxation in a gas involves physical contact collisions (or at least, very short-range interactions) between molecules, meaning a molecule's velocity can be changed significantly after a single collision. (b) In a stellar system, by contrast, two-body encounters are mostly weak, long-range nudges. Later in the article, we will deal with collective effects, whereby one must think beyond bare two-body processes and account for the coherent oscillations of very large numbers of stars at once.</figcaption></figure>

relaxation, and primarily following Sec. 1.2.1 of BT08. (For a more complete treatment of Chandrasekhar theory, see Chap. 7 of BT08.) What we will present is a gravitational version of Rutherford scattering.

To begin, we will ignore all the complications associated with real stellar systems and instead consider a single “test” star moving in a statistically homogeneous box of static “field” stars of mass m. Let the test star travel in a perfectly straight line until it approaches a field star at impact parameter b. The goal now is to calculate the nudge to the test star’s velocity, $\delta v(b)$.

Since we are going to assume that the interaction between the two stars is weak, to lowest order we can approximate the trajectory of the test star as a straight line throughout the encounter (see Fig. 5). Suppose without loss of generality that the closest approach of the test star to the field star occurs at $t = 0$. We may then calculate the total impulse it receives by integrating from $t = -\infty$ to $t = \infty$. Then, the acceleration the test star feels in the direction parallel to v while flying toward the field star ($t < 0$) is exactly canceled by the retardation it feels as it flies away ($t > 0$), meaning that there is no overall velocity change in the parallel direction. Hence, we need only to calculate the total perpendicular velocity nudge $\delta v_{\perp}$.

The perpendicular acceleration field felt by the test star is (see Fig. 5 for the definition of r and θ)

$$F_{\perp} = \frac{Gm}{r^2} \cos \theta = \frac{Gm}{b^2} \left[ 1 + \left( \frac{vt}{b} \right)^2 \right]^{-3/2} . \quad (15)$$

Integrating over all time, we get

<figure><img src="image_18.png" /><figcaption>FIG. 5. A test star encounters a field star of mass m with impact parameter b. Assuming the encounter to be weak, the lowest-order nudge to the test star's velocity, $\delta v$, can be calculated by integrating along the unperturbed trajectory (straight dotted line).</figcaption></figure>

$$\delta v_{\perp} = \int_{-\infty}^{\infty} dt F_{\perp} = \frac{Gm}{bv} \int_{-\infty}^{\infty} \frac{ds}{(1 + s^2)^{3/2}} = \frac{2Gm}{bv} . \quad (16)$$

As one might expect, the more massive the field star, and the smaller the impact parameter b, the larger the velocity impulse the test star experiences. Moreover, it is sensible that the impulse decays with v, since the larger is v the shorter is the duration of the encounter. On the other hand, the formula (16) also predicts its own failure: if $\delta v_{\perp}$ is sufficiently large that it is comparable to the original speed v, then the assumed straight-line trajectory is obviously a poor approximation. Note that for a fixed m and v, this breakdown occurs at a critical impact parameter $b_{90} = 2Gm/v^2$, so-called because this in fact produces a 90-degree deflection. For impact parameters $b \le b_{90}$, one ought really to solve the two-body encounter problem without making any approximations (which is not difficult—see Sec. 3.2 of BT08). Luckily, such strong encounters are rare and make little difference to the bulk evolution of the system. We will ignore them from now on.

Chandrasekhar’s next step is to assume that relaxation in a stellar system consists of nothing more than an uncorrelated series of these weak two-body encounters. Under this assumption, we can calculate the relaxation time as follows. Let the number density of stars in the box be n. Then, in a time T, the field star passes a total of $\approx nv T 2\pi b db$ field stars at impact parameter $\in [b, b + db]$. The mean square kick to the velocity that the test star receives in time T is therefore

$$\sum \delta v_{\perp}^2 \approx \int_{b_{\max}}^{b_{\min}} nv T 2\pi b db \times \delta v_{\perp}^2 = \frac{8\pi G^2 m^2 n}{v} T \ln \Lambda, \quad (17)$$

where $\ln \Lambda \equiv \ln (b_{\max}/b_{\min})$ is the Coulomb logarithm, which is familiar from the corresponding plasma calculation of charged particle scattering.¹³

We seem to have a problem here, since taking either $b_{\min} \to 0$ or $b_{\max} \to \infty$, or both, causes the total velocity impulse to diverge (albeit in the most undramatic way, i.e., logarithmically). So, what values of $b_{\min}$ and $b_{\max}$ should we take? Usually, one takes $b_{\max}$ to be comparable to the scale of the system $\sim R$, since there cannot be any encounters more distant than this (this is different from plasma, in which one can safely take $b_{\max}$ of order the Debye screening length¹³). Moreover, one typically takes $b_{\min} \approx b_{90}$ because, as we mentioned above, strong encounters are usually rare and unimportant. In fact, while setting the Coulomb logarithm accurately is an important problem for those modelling globular cluster evolution,¹⁴ from a theorist’s standpoint the answer is it does not matter too much because the model is only heuristic anyway. For instance, encounters with impact parameters b on the
order of the size of the system certainly cannot be thought of as straight lines, since they must take into account the curved geometry of the true orbits of the stars (Sec. II C 3). The standard choice $\Lambda \approx R/b_{90}$ is good enough for obtaining order-of-magnitude estimates, which is all we care about here.

If we now set $\sum \delta v_{\perp}^2$ equal to the original squared speed $v^2$, we can identify $T = t_{2BR}$ as the two-body relaxation time

$$t_{2BR} \approx \frac{v^3}{8\pi G^2 m^2 n \ln \Lambda} \qquad (18)$$

Let us relate this more clearly to the number of stars in the system N. Let the side-length of the box be $\approx R$—which is a proxy for the radius of the star cluster or galaxy in which we are interested. Then, we know from the virial theorem that typically $v^2 \approx GNm/R$ [see Eq. (8)]. We may also set the typical crossing time to be $t_{dyn} \approx R/v$ and the typical number density to be $n \approx N/(4\pi R^3/3)$. Then, taking $b_{max} = R$ and $b_{min} = b_{90} = 2Gm/v^2$, we get $\Lambda \approx Rv^2/(Gm) \approx N$, so that

$$t_{2BR} \approx \frac{0.1 N}{\ln N} t_{dyn} \qquad (19)$$

To be even more heuristic, we will often drop the factors of 0.1 and $\ln N$ and simply state that the two-body relaxation time is on the order of $t_{2BR} \sim Nt_{dyn}$.

Plugging in characteristic numbers from Table I, we find that for a globular cluster $t_{2BR} \sim 10$ Gyr, comparable to the age of the cluster. This suggests that the simple physics of uncorrelated two-body encounters is of major importance for determining the evolution of globular clusters. However, for a disk galaxy, this theory suggests evolution should only occur on the excruciatingly long timescale $t_{2BR} \sim 10^{10}$ Gyr, nine orders of magnitude longer than the age of the Universe.¹⁵ Yet, we know from observations and simulations that disk galaxies can evolve significantly on a timescale of $\ge$ a few Gyr. The estimate (19) tells us that star-star encounters are much too inefficient to have produced this evolution, so a more sophisticated kinetic theory is required (Sec. II C 3).

Note that Eq. (19) is directly analogous to the inverse two-body collision time in an electrostatic plasma, $t_{coll} \sim (\Lambda_p/\log \Lambda_p)\omega_p^{-1}$, where $\Lambda_p$ is the plasma parameter (the number of particles per Debye sphere) and $\omega_p$ is the plasma frequency.¹⁶ Thus, there is a direct analogy between the plasma parameter $\Lambda_p$ and the number of stars N (at least for stellar systems whose evolution is driven purely by finite-N fluctuations, with little collective amplification). In particular, the inverse of these quantities, $1/\Lambda_p \ll 1$ and $1/N \ll 1$, are the small numbers with which one can expand the BBGKY hierarchy in weakly coupled plasmas and stellar systems, respectively.¹⁷ These numbers being small are equivalent to the statement that the kinetic energy of each particle is large compared to the fluctuations in the potential energy (recall that the mean potential energy in a stellar system is of order the kinetic energy, see Sec. II A).

## 2. Collisional and collisionless

It is worth commenting here on the (often extremely confusing) use of the words collisional and collisionless when describing stellar systems. These words were partly inherited from plasma theory, where they have a fairly concrete meaning: roughly, interactions that occur on scales smaller than the Debye scale are collisions, and if these can

be ignored, then the system is collisionless. In a stellar system, however, there is no shielding, no Debye length, and no natural way to decouple small-scale and large-scale gravitational potential fluctuations—every star is interacting gravitationally with every other star, always.

Sometimes in the gravitational context, the word collisionless is used as a shorthand for “a system whose naive two-body relaxation time (19) is much longer than its age” (e.g., a disk or an elliptical galaxy), while collisional systems are those with $t_{2BR}$ less than or comparable to their age (basically, star clusters). Alternatively, the word “collisional” is sometimes used to refer to self-interacting dark matter (SIDM),¹⁸ a model of dark matter involving a finite cross section for two dark matter particles to annihilate, by contrast with standard collisionless dark matter, which has no such cross section. We could give many other examples of these words being used for slightly different purposes. It is therefore no surprise that this semantic issue frequently causes confusion among nonexperts.

There are two sensible ways to deal with the issue. The first way is to call a system collisional if and only if the finite number of particles N matter to its evolution, and to call it collisionless otherwise. The second way is to avoid using the words collisional and collisionless altogether. We advocate the second way, although these words are so ingrained in the literature that sometimes this is impossible.

## 3. Where do we go from here?

It turns out that when applied to the evolution of hot spherical clusters, Chandrasekhar’s picture does a remarkably good job of estimating the true relaxation time $t_{relax}$ and various other evolutionary properties,¹⁴ although the value of the Coulomb logarithm has to be tweaked.²²,²³ Moreover, one can use Chandrasekhar’s ideas to spin up a full kinetic theory—i.e., a kinetic equation of the form $\partial f/\partial t = C[f]$, with C calculated by assuming relaxation consists of two-body encounters as above—and this works well in hot spheres too.²⁴ To demonstrate this, we present Fig. 6 (adapted from Joshi et al.¹⁹), which shows the evolution of “Lagrange radii” of a spherical globular cluster in a numerical simulation. These Lagrange radii are the radii enclosing a fixed percentage of the cluster’s mass at any given time. First of all, we see that the inner Lagrange radii gradually decrease with time (corresponding to the contraction of the cluster core, leading ultimately to the “gravothermal catastrophe”²⁵) and the outer radii slowly increase (corresponding to the expansion of the cluster’s envelope). This evolution happens very slowly (the time unit on the horizontal axis is the two-body relaxation time $t_{2BR} \gg t_{dyn}$), corresponding to the gradual evolution from one virialized state to another. Looking more closely, one can see that for each Lagrange radius in the panel there are actually two lines plotted. One of these was calculated from a direct N-body simulation of the “exact” equations of motion of all the stars, whereas the other employed a Monte-Carlo calculation based upon Chandrasekhar’s theory of two-body encounters. The fact that across much of the figure, these two lines are indistinguishable is testament to the success of Chandrasekhar theory in describing hot spherical clusters.

On the other hand, as mentioned at the end of Sec. II C 1, Chandrasekhar’s theory fails completely to describe the evolution of cold disk galaxies. First of all, disk galaxies develop long-range collective structures like the spiral arms in Fig. 1(a), and these have no place in Chandrasekhar’s picture. Also, Chandrasekhar’s theory drastically
<figure><img src="image_19.png" /><figcaption>**FIG. 6.** Long-term evolution of a globular cluster. Adapted from Joshi et al., Astrophys. J. 540, 969 (2000). Copyright 2000 IOP.¹³ The lines show the evolution of the “Lagrange radii” which enclose 0.35%, 1%, 3.5%, 5%, 7%, 10%, 14%, 20%, 30%, 40%, 50%, 60%, 70%, and 80% of the cluster mass, respectively. In each case, both the results of a direct N-body simulation and of a Monte-Carlo calculation²⁰ based on Chandrasekhar’s theory are shown. For most times and Lagrange radii, the two sets of results are indistinguishable.</figcaption></figure>

overestimates the timescale on which these systems relax. As an example of this, consider Fig. 7 (adapted from Frankel et al.²¹). This figure was generated by fitting a simple diffusion model to the orbital actions (to be introduced in Sec. III C) of stars in the Milky Way using data from Gaia, a space observatory, which has been mapping the kinematics of stars in our Galaxy since 2014. Panel (a) shows the diffusion of

the stellar angular momenta (essentially the guiding radius of near-circular orbits around the galaxy), while panel (b) shows the evolution of those same stars' mean radial actions (essentially their orbital “eccentricities”). In both panels, we see that significant evolution can occur on the timescale of a few Gyr. To get an estimate for the true relaxation time in this system, note that the orbital angular momenta of stars on circular orbits (speed $v_φ ≈ 250$ km s⁻¹) near the Solar radius in the Milky Way ($R ≈ 8$ kpc) is $L = Rv_φ ≈ 2000$ kpc km s⁻¹. The relaxation time we infer from Fig. 7(a) is therefore of order $t_{relax} ≈ (2000/619)² × 6$ Gyr ≈ 60 Gyr, and many orders of magnitude shorter than would be suggested by the Chandrasekhar estimate (19). Even though this relaxation time is still several times longer than the age of the Galaxy (~10 Gyr), it can still have a major effect on galaxy evolution, as we will see in Sec. VIII B 1.

Chandrasekhar’s theory fails because it is based on a number of false assumptions, for instance:

* Stellar systems are inherently *inhomogeneous* (see Fig. 1) but the theory assumed homogeneity.

* Because of the inhomogeneity, unperturbed stellar orbits are not endless straight lines, but instead are *quasiperiodic* (Sec. III). This leads to *resonant* interactions.

* Stars do not behave in an independent, uncorrelated way but instead tend to oscillate *collectively*. Thus, one cannot replace all interactions between stars with a superposition of uncorrelated two-body encounters.

These statements apply to both globular clusters and disk galaxies, which raises the question: Why does Chandrasekhar’s theory work to describe the former but not the latter? The answer lies in the fact that globular clusters are dynamically hot, whereas disk galaxies are dynamically cold. The hotness of clusters means that they act like thermal plasmas: collective effects are weak and successive fluctuations are mostly uncorrelated, as different parts of the system have no way to efficiently communicate with one another. On the other hand, the highly ordered orbital structure of disk galaxies means that great numbers of stars can have long-lived coherent dynamical interactions. The result is that disk galaxies, like nonequilibrium plasmas, are replete with large-scale collective oscillations and possible instabilities. Such effects are able to redistribute energy and angular momentum much more efficiently—and in a fundamentally different way—than the mechanism of incoherent two-body relaxation.

<figure><img src="image_20.png" /><figcaption>**FIG. 7.** Adapted from Frankel et al., Astrophys. J. 896, 15 (2020). Copyright 2020 IOP.²¹ Each panel shows the evolution of an orbital integral of motion [namely, (a) the component of angular momentum perpendicular to the galactic disk and (b) the radial action, which will be introduced properly in Sec. III C] of stars in the Milky Way over 6 Gyr, as inferred from Gaia data. In the right panel, different colors correspond to groups of stars with different initial Galactocentric radii.</figcaption></figure>
It follows that to describe the behavior of disk galaxies, we need to leave Chandrasekhar's theory behind and develop a new kinetic theory, one capable of dealing with the complications listed above. That is one of the aims for the rest of this tutorial article. However, before we can do kinetics, we must learn how to describe the orbits of stars in inhomogeneous, but still fixed and smooth, mean field gravitational potentials (Sec. III), and then figure out what happens to those orbits when they are perturbed (Sec. IV).

### III. ORBITS IN MEAN FIELD POTENTIALS

As we mentioned in Sec. II B, on the timescale ~ $t_{dyn}$ the orbits of stars in a virialized N-body system can be very well described by ignoring the discrete nature of the system and instead treating individual stars as test particles orbiting in a smooth potential. This is a good approximation for two reasons: (i) the gravitational force between two objects is long range, decreasing as the square of the distance between them, and (ii) the number of stars in the system, N, is very large. The consequence of (i) is that the force acting on any individual star is typically not dominated by its nearest neighbors, but rather by the collective gravity of simultaneous distant interactions. The consequence of (ii) is that the Poisson fluctuations in the number of stars in a given sub-volume tends to be small, so the collective gravitational force due to stars in this sub-volume fluctuates only weakly on the timescale $t_{dyn}$. The average value of this collective force, combining all the sub-volumes, is then typically very smooth (see also Sec. 1.2 of BT08).

Another key feature of virialized stellar systems is that they have rather simple bulk geometry: they are normally approximately spherical (in the case of globular clusters, some dwarf galaxies, etc.), or thin and axisymmetric (like a disk galaxy), or possibly ellipsoidal/triaxial (like some massive elliptical galaxies). This simple geometry leads to (mostly) simple, regular orbits, and not many chaotic orbits. In this section, we will describe how stars orbit within a simple class of smooth mean field potentials.

#### A. Mean field orbits

In elementary courses on plasma physics, there is no great heed paid to “mean field dynamics” because those dynamics are trivial: in a homogeneous electrostatic plasma, mean field “orbits” are just straight lines: $\mathbf{r}(t) = \mathbf{v}t + \mathbf{r}(0)$ and $\mathbf{v} = \mathbf{v}(0) = \text{cst}$. On the contrary, stellar systems are inherently inhomogeneous, which means that the unperturbed mean field orbits of stars are not straight lines. This is something plasma physicists also have to deal with when considering orbits in systems with nontrivial geometry and especially with magnetic fields, such as in toroidal tokamak fusion plasmas or the Earth’s magnetosphere.²⁶⁻²⁸

If the gravitational potential is $\Phi(\mathbf{r})$, the equations of motion of a test star are

$$ \frac{d\mathbf{r}}{dt} = \mathbf{v}, \qquad (20a) $$

$$ \frac{dv}{dt} = -\frac{\partial \Phi}{\partial \mathbf{r}}. \qquad (20b) $$

In practice, we cannot solve these equations of motion explicitly except in a few very simple potentials (e.g., the Kepler potential $\Phi(\mathbf{r}) \propto -1/|\mathbf{r}|$, and the harmonic potential $\Phi(\mathbf{r}) \propto |\mathbf{r}|^2$). However, we can solve the equations via numerical integration—this is called “integrating the orbit.” Modern numerical packages like Galpy,²⁹ Gala,³⁰ and

AGAMA³¹ allow one to integrate and analyze mean field orbits in pre-defined potentials with ease.

In general, numerical orbit integration reveals a wide array of possible orbit families—see Chap. 3 of BT08 for an overview. However, for simplicity, throughout this article, we will only consider one simple orbit family, which is the planar rosette orbit shown in Fig. 8.

This is (nearly) the only type of orbit that is possible in central potentials $\Phi(r)$, where $r \equiv |\mathbf{r}|$. It is therefore the overwhelmingly dominant orbit family in both our mental models, the hot sphere and the cold axisymmetric disk (Fig. 2). By symmetry, motion in central potentials conserves each component of the specific orbital angular momentum vector, $\mathbf{L} = \mathbf{r} \times \mathbf{p}$, and so each such orbit is confined to a two-dimensional plane (which we can take to be the plane $Z = 0$ without loss of generality). From now on, we describe orbits in this plane with a standard $(\phi, R)$ polar coordinate system centered at $r = 0$.

#### 1. Orbits in central potentials

We will now show how to describe the rosette-type orbit of Fig. 8 quantitatively. We see from Fig. 8 that the orbit consists of a combination of an azimuthal circulation in the positive $\phi$ direction, and a radial oscillation in $R$ between the fixed pericenter $R_p$ and apocenter $R_a$. We can calculate $R_p$ and $R_a$ in terms of the two conserved quantities of the star’s motion, namely, its specific energy and specific angular momentum (we will drop the word “specific” from now on)

$$ E = \frac{1}{2} \left( \frac{dR}{dt} \right)^2 + \frac{1}{2} L^2 + \Phi(R), \qquad (21a) $$

<figure><img src="image_21.png" /><figcaption>FIG. 8. Nearly all orbits in central potentials $\Phi(r)$ look like this: a two-dimensional planar “rosette.” Using in-plane polar coordinates $(\phi, R)$, each orbit is a combination of an azimuthal oscillation in the $\phi$ direction and a radial oscillation between $R_p$ and $R_a$. The thick black line shows the trajectory over one radial period, from $R_a$ to $R_p$ and back again.</figcaption></figure>
$$L = R^2 \frac{d\phi}{dt}. \qquad (21b)$$

The peri/apocentre are the radial turning points of the orbit; that is, they are solutions to $v_R(R_p) = v_R(R_a) = 0$, where $v_R = dR/dt$ is the radial velocity

$$v_R(R) \equiv \sqrt{2[E - \Phi(R)] - L^2/R^2}. \qquad (22)$$

In central potentials, there is then a one-to-one correspondence between the integrals of motion (E, L) and the peri/apocentric distances ($R_p, R_a$):

$$E(R_p, R_a) = \frac{R_a^2 \Phi(R_a) - R_p^2 \Phi(R_p)}{R_a^2 - R_p^2}, \qquad (23a)$$

$$L(R_p, R_a) = \sqrt{\frac{2[\Phi(R_a) - \Phi(R_p)]}{R_p^{-2} - R_a^{-2}}}. \qquad (23b)$$

We can also calculate the frequencies of the radial and azimuthal oscillations as follows. First, the radial frequency is

$$\Omega_R = \frac{2\pi}{T_R}, \qquad (24)$$

where $T_R$ is the radial period

$$T_R(R_p, R_a) = 2 \int_{R_p}^{R_a} \frac{dR}{v_R(R)}. \qquad (25)$$

The azimuthal frequency is

$$\Omega_{\phi} = \Omega_{R} \frac{|\Delta\phi|}{2\pi}, \qquad (26)$$

where $|\Delta\phi|$ is the increment in azimuthal angle over one radial period reading

$$|\Delta\phi(R_p, R_a)| = 2L \int_{R_p}^{R_a} \frac{dR}{R^2 v_R(R)}. \qquad (27)$$

(In all sensible spherical potentials, the increment $|\Delta\phi(R_p, R_a)|$ lies between $\pi$ and $2\pi$). In other words, all the characteristic quantities describing the orbit can be calculated given the two numbers ($R_p, R_a$), or equivalently the two constants of motion (E, L). This is in contrast to a homogeneous plasma, where one can label an orbit just by its velocity $v$.

To further illustrate what these rosette orbits look like, consider motion in the spherical Plummer potential, which is often used as a simple model for the gravitational potentials in globular clusters:

$$\Phi(r) = - \frac{GM}{\sqrt{r^2 + b^2}}, \qquad (\text{Plummer potential}). \qquad (28)$$

Here, $b$ is the scale radius of the potential (for a globular cluster $b \sim 1$ pc—see Table 1). The Plummer potential has very simple asymptotic forms. For $r \gg b$, it approximates the Keplerian potential of a point mass, $\Phi \approx -GM/r$. In this potential, all orbits are closed ellipses with focus at the origin, i.e., they have $|\Delta\phi(R_p, R_a)| = 2\pi$. For $r \ll b$, the (nonconstant part of the) Plummer potential is quadratic $\Phi \propto r^2$, which corresponds to motion in a spherical “harmonic core” of constant density (see BT08, Sec. 3.1). All orbits in the harmonic core are also closed ellipses, but this time, they are centered on the origin, and have $|\Delta\phi(R_p, R_a)| = \pi$.

In Fig. 9, we show three orbits integrated in the potential (28) in the $(\phi, R)$ plane, each with the same $R_a = b$ but different $R_p$, as labeled above each panel.

Orbit (a) is nearly radial, with $R_p \ll R_a$, whereas orbit (c) is nearly circular, with $R_p \approx R_a$, and orbit (b) is somewhere in-between.

The orbit labels ($R_p, R_a$), or equivalently (E, L), contain no information about the orbit’s instantaneous radial or azimuthal phase, just as labeling an orbit by $v$ in a plasma does not tell us its instantaneous position in the box $\mathbf{r}$. However, it turns out that the phase of the orbits also exhibit rather simple behavior, as we investigate next.

## B. Quasiperiodicity

A crucial feature of orbits in many smooth, regular galactic potentials revealed by numerical orbit integration is that they are

<figure><img src="image_22.png" /><figcaption>FIG. 9. Three orbits in the Plummer potential (28). In each case, the orbit starts at its apocentric distance at the location $(X, Y) = (b, 0)$. Each orbit has the same apocenter $R_a = b$ (red dashed circles) but different pericenters $R_p$ (blue dashed circles), as indicated above each panel. As a result, orbit (a) is almost radial ($J_R \gg L$), orbit (c) is almost circular ($J_R \ll L$), and orbit (b) is somewhere in-between ($J_R \sim L$)—see Sec. III C 3 for the definition of the radial action $J_R$.</figcaption></figure>
quasiperiodic. By this we mean that, very often, one can Fourier-analyze a numerically integrated orbit and write it as

$$ \mathbf{r}(t) = \sum_{n} \mathbf{r}_n e^{in\Omega t}. \qquad (29) $$

Here, the sum is over all *d*-dimensional vectors comprised of integers, where *d* is the number of dimensions, e.g., (0, 0, 1), (0, 1, 1), ... etc., in the case of *d* = 3. The vector $\Omega$ is made up of *d* frequencies that are characteristic of the orbit in question. Quasiperiodicity is also a common feature of particle orbits in tokamaks, magnetospheres, and so on.

That the orbits in Fig. 9 are quasiperiodic might have been guessed by eye, but this fact also extends to a much broader range of orbits in realistic galactic potentials.32,33 Mathematically, that an orbit is quasiperiodic can be proven in cases where the number of known independent isolating integrals of motion is (at least) equal to the number of dimensions (which is always true for orbits in spherically symmetric potentials). However, it is a nontrivial empirical fact that quasiperiodicity also extends to orbits in many realistic axisymmetric potentials $\Phi(R, Z)$ for which the orbit is fully three-dimensional and there are only two obvious independent constants of motion (E and the component of angular momentum along the symmetry axis, $L_Z$).32,33 Normally in a system with more dimensions than there are independent integrals of motion, one would expect some degree of chaos; instead, it is believed that orbits in axisymmetric potentials almost always possess a third independent integral of motion, whose analytic form we do not know but which makes the orbits regular and hence quasiperiodic (though this is not generically guaranteed34).

While triaxial potentials (like ellipsoids) can also harbor families of quasiperiodic orbits, they also tend to have many more chaotic orbits too, which arise in the parts of phase space between the regions of distinct quasiperiodic families. As such, in the rest of this tutorial article, we will suppose that the exact potential of the system is close to, or exactly, axisymmetric, and that in the limit of pure axisymmetry, there is only one regular (quasiperiodic) family of orbits we need to worry about, and no chaos. The central potentials $\Phi(r)$ are of this nature.

### C. Angle-action variables

The upshot of quasiperiodicity is that stellar orbits can be thought of as a superposition of oscillators with fixed frequencies $\Omega$. Here, we will show how this fact allows us to construct a coordinate system in which this mean field orbital motion is trivial, namely, angle-action coordinates. An alternative and more detailed introduction to these coordinates and how they are used can be found in Sec. 3.5 of BT08.

Let $(\mathbf{q}, \mathbf{p})$ be a set of canonical coordinates and momenta on phase space (a reader who is not familiar with canonical coordinates can consult Sec. D.4 of BT08). Now consider a star moving according to the generic Hamiltonian $H(\mathbf{q}, \mathbf{p}, t)$. The equations of motion for that star are

$$ \frac{d\mathbf{q}}{dt} = \frac{\partial H}{\partial \mathbf{p}}, \qquad (30a) $$

$$ \frac{d\mathbf{p}}{dt} = -\frac{\partial H}{\partial \mathbf{q}}. \qquad (30b) $$

We can also use these canonical coordinates to introduce a bilinear antisymmetric functional on phase space called the Poisson bracket. For any two functions $g(\mathbf{q}, \mathbf{p})$ and $h(\mathbf{q}, \mathbf{p})$, we define

$$ [g, h] = \frac{\partial g}{\partial \mathbf{q}} \cdot \frac{\partial h}{\partial \mathbf{p}} - \frac{\partial g}{\partial \mathbf{p}} \cdot \frac{\partial h}{\partial \mathbf{q}}. \qquad (31) $$

This will be useful later.

A perfectly valid choice for the canonical coordinates $(\mathbf{q}, \mathbf{p})$ is the position and velocity $(\mathbf{r}, \mathbf{v})$. However, the lack of straight line orbits in stellar systems means that $(\mathbf{r}, \mathbf{v})$ are awkward variables with which to formulate a kinetic theory. We would like instead to work as much as possible with quantities that are invariant along orbits. An integral of motion $I(\mathbf{r}, \mathbf{v})$ is any function of phase space coordinates that is constant along an orbit. As argued in Sec. III B, for all the mean field orbits we care about in this tutorial article (that is, orbits that are quasiperiodic) and there exist at least as many independent integrals of motion as there are degrees of freedom in the system (normally 2 or 3). One can construct as many mutually dependent integrals of motion as one wishes, since any function of integrals of motion is itself an integral of motion.

The obvious next question is: given that these integrals of motion exist, can we use them as coordinates for our phase space? The answer is yes, but if we want to preserve the Hamiltonian structure of our system—i.e., if we want the new coordinates to be canonical—then we have to be careful about which integrals to use. It turns out that there is a special set of integrals called actions $J$, which are canonical momenta on phase space,35 and can be complemented by canonical coordinates. Specifically, these actions have the special property that the Hamiltonian of the underlying quasiperiodic system $H(\mathbf{q}, \mathbf{p})$ can be written exclusively in terms of them, $H = H(J)$, and they can be complemented by canonical coordinates or angles $\theta$, which are $2\pi$-periodic. In these special angle-action coordinates $(\mathbf{q}, \mathbf{p}) = (\theta, J)$, Hamilton's equations (30) take the simple form

$$ \frac{d\theta}{dt} = \frac{\partial H}{\partial J} = \Omega(J), \qquad (32a) $$

$$ \frac{dJ}{dt} = -\frac{\partial H}{\partial \theta} = 0, \qquad (32b) $$

where we introduced the orbital frequencies $\Omega(J) = \partial H/\partial J$. In these coordinates, the orbits are straight, horizontal lines, as one has

$$ \theta(t) = \theta(0) + \Omega(J)t \quad (\text{mod } 2\pi), \qquad (33a) $$

$$ J(t) = J(0) = \text{cst.} \qquad (33b) $$

Some further properties of the angle-action coordinates are as follows. Since they are canonical, they leave the Poisson bracket (31) invariant, so that one generically has

$$ [g, h] = \frac{\partial g}{\partial \theta} \cdot \frac{\partial h}{\partial J} - \frac{\partial g}{\partial J} \cdot \frac{\partial h}{\partial \theta}. \qquad (34) $$

As a result, angle-action coordinates also leave the infinitesimal phase space volumes invariant so that

$$ d\mathbf{q} \wedge d\mathbf{p} = d\mathbf{r} \wedge d\mathbf{v} = d\theta \wedge dJ. \qquad (35) $$

Furthermore, the periodicity of the angles in *d* dimensions means that any phase space function can be Fourier expanded as
$$G(\theta, \mathbf{J}) = \sum_{\mathbf{n}} G_{\mathbf{n}}(\mathbf{J}) e^{i\mathbf{n}\cdot\theta}, \qquad (36)$$

with inverse

$$G_{\mathbf{n}}(\mathbf{J}) = \int \frac{d\theta}{(2\pi)^d} G(\theta, \mathbf{J}) e^{-i\mathbf{n}\cdot\theta}. \qquad (37)$$

Actions are also adiabatic invariants: roughly speaking, if $H$ evolves on a timescale longer than the dynamical time, $t_{dyn} \sim 1/|\Omega|$, an orbit governed by $H$ will evolve in such a way that $\mathbf{J} = \text{cst.}$ (but see Weinberg$^{36-38}$ for a more precise statement).

Note that we have not proven the existence of the angle-action variables, nor have we explained how to generate them from scratch (via the process of writing down a suitable canonical generating function). All this is standard material which is explained rather neatly (if heuristically) in Sec. 3.5 of BT08, or much more rigorously in the books by Lichtenberg and Lieberman$^{39}$ and Arnold.$^{40}$ Throughout the rest of this tutorial article, we will simply assume that the angle-action coordinates can be constructed for the problem at hand.

Also, it is important to note that angles and actions are just coordinates, and that once we have constructed them, we can use them to describe the motion of particles in Hamiltonians other than $H(\mathbf{J})$. In particular, we often care about systems in which the motion is nearly, but not perfectly, quasiperiodic (such as a spherical stellar system perturbed by an infalling satellite). In such a case, we can construct the angle-action variables $(\theta, \mathbf{J})$ corresponding to the corresponding perfectly periodic system and then describe the real (perturbed) system by expressing its Hamiltonian in these coordinates, $H(\theta, \mathbf{J})$, deriving the corresponding equation of motion $d\theta/dt = \partial H/\partial \mathbf{J}$ and $d\mathbf{J}/dt = -\partial H/\partial \theta$, etc. We will see how this construction works in practice starting in Sec. V B.

## 1. Example: Angle-action variables for a multi-periodic homogeneous system

As a first example, we consider a fictitious homogeneous stellar system, which provides a direct analogy with homogeneous plasmas. To this end, let us assume that our stellar system is placed within a

periodic three-dimensional box of side-length $L$. We also assume that the mean potential vanishes, i.e., $\Phi_0 = 0$, so that unperturbed trajectories are straight lines. In this case, the system's angle-action coordinates and the associated orbital frequencies are given by

$$\theta = \frac{2\pi}{L} \mathbf{r}, \quad \mathbf{J} = \frac{L}{2\pi} \mathbf{v}, \quad \mathbf{\Omega} = \frac{2\pi}{L} \mathbf{v}. \qquad (38)$$

This makes clear a crucial correspondence between (fictitious) homogeneous stellar systems and homogeneous electrostatic plasmas. In both systems, mean field orbits are labeled by their velocity $\mathbf{v}$, and, up to a prefactor, these velocities are also the orbit's intrinsic frequencies $\mathbf{\Omega}$. This means that (i) the gravitational/electrostatic potential $\Phi(\mathbf{r})$ is a function only of angles $\theta$ and not actions $\mathbf{J}$, and (ii) the orbital frequencies $\mathbf{\Omega}(\mathbf{J})$ are trivial (linear) functions of the actions. These two features are special to homogeneous systems, and greatly simplify the kinetic description (see, e.g., Sec. VIB 1).

## 2. Example: Angle-action variables for the 1D harmonic oscillator

Let us now consider the case of a one-dimensional harmonic oscillator whose Hamiltonian is $H = (v^2 + \omega^2 x^2)/2$. As illustrated in Fig. 10, in this case, orbits take the form of concentric circles in the phase space $(x, v/\omega)$.

We define the action via $H = J\omega$, along with the angle $\theta$ so that

$$x = \sqrt{\frac{2J}{\omega}} \sin \theta, \quad v = \sqrt{2J\omega} \cos \theta. \qquad (39)$$

Following Eq. (34), one can check that $[\theta, J] = 1$; thus, $(\theta, J)$ are canonical coordinates. Finally, orbits which followed circles in $(x, v/\omega)$ space now traverse straight horizontal lines in $(\theta, J)$ space.

## 3. Example: Angle-action variables for orbits in a central potential

For central potentials, $\Phi(R)$, the two-dimensional rosette motion can be described using action variables $J_1 = L$ and $J_2 = J_R$, where $L$ is the usual angular momentum (see Eq. 21) and $J_R$ is the radial action

<figure><img src="image_23.png" /><figcaption>FIG. 10. Phase space diagram of a harmonic oscillator, $H = (v^2 + \omega^2 x^2)/2$. Panel (a) shows particle trajectories in the physical phase space $(x, v/\omega)$, which take the form of concentric circles. Here, the action $J$ should therefore be seen as a label for the circle associated with the orbit, and the angle $\theta$ should be seen as the position along the circle. Panel (b) shows the same trajectories in angle-action space $(\theta, J)$. In these coordinates, the mean field motions are straight lines. The action $J$ is conserved, while the angle $\theta$ evolves linearly in time.</figcaption></figure>
$$J_R = \frac{1}{\pi} \int_{R_p}^{R_a} dR v_R(R), \qquad (40)$$

with the radial velocity, $v_R$, given by Eq. (22). This is directly analogous to the longitudinal adiabatic invariant $J = m \oint v_\parallel ds$ for particle motion in a magnetic mirror.41 The conjugate angles are the azimuthal and radial angles

$$\theta_1 = \theta_\phi = \phi + \int_{R_p}^{R} \frac{dR' \frac{\Omega_\phi - L/R'^2}{v_R(R')}}{v_R(R')}, \quad (41a)$$

$$\theta_2 = \theta_R = \Omega_R \int_{R_p}^{R} \frac{dR'}{v_R(R')}. \qquad (41b)$$

These angles evolve at the rate of the azimuthal frequency $\Omega_\phi$ and radial frequency $\Omega_R$, respectively.

Figure 11 (taken from Ref. 42) presents recent data from the spacecraft Gaia, showing the distribution of millions of stellar orbits in our local patch of the Milky Way disk, plotted in various slices through angle-action space.

Each panel is full of structure, some of which is due to selection effects (since we are unable to access data across the whole Galaxy) but much of which is due to dynamical interactions between stars, spiral arms, the Galaxy's central bar, and so on. Though we will not pursue any details here, understanding structures like those in Fig. 11 is a major goal of modern Galactic astronomy.

<figure><img src="image_24.png" /><figcaption>FIG. 11. The distribution of stars' orbits in our local patch of the Milky Way disk (a.k.a. the Solar neighborhood), plotted in angle-action space, using data from the Gaia satellite. This figure was reproduced from Hunt et al., Mon. Not. R. Astron. Soc. 490, 1026 (2019). Copyright 2019 Royal Astronomical Society.42 Another detailed guide to features in the action space (top left panel) can be found in Ref. 43.</figcaption></figure>
### 4. Example: Angle-action variables for near-circular orbits

We have introduced angle-action coordinates as the natural variables with which to describe quasiperiodic orbits, since those orbits are superpositions of (generally anharmonic) oscillations at discrete frequencies $\Omega(\mathbf{J})$. However, a disadvantage of angle-action coordinates is that we can rarely make explicit the mapping $(\mathbf{r}, \mathbf{v}) \leftrightarrow (\theta, \mathbf{J})$. On the other hand, we can write down an explicit mapping if we restrict ourselves to near-circular orbits (which is a highly relevant regime for describing cold disks). There, we can employ the epicyclic approximation. In this limit, each orbit consists of circular motion at the guiding center radius $R_g$, plus a superposition of harmonic oscillations in the radial and azimuthal coordinates. In Appendix A, we show how to construct explicit angle-action coordinates in the epicyclic approximation.

### D. Summary and outlook

For the remainder of this tutorial article, we assume that we can always construct angle-action variables $(\theta, \mathbf{J})$, such that the mean field Hamiltonian depends only on $\mathbf{J}$ (which is itself a nontrivial technical problem$^{44}$). It is then tempting to think of actions as the “velocity” variables in our new phase space and the angles as our new “positions.” Indeed, many formulas from plasma kinetics will more-or-less hold in stellar dynamics just by making the substitutions $\mathbf{v} \to \mathbf{J}$ and $\mathbf{r} \to \theta$ (see especially Table II). However, there are two key mathematical differences that make stellar kinetics much harder than plasma kinetics. The first is that unlike in a homogeneous plasma, the canonical coordinates change at rates that are in general, nonlinear functions of the canonical momenta, i.e., $\Omega(\mathbf{J})$ is typically a complicated (nonlinear) function of $\mathbf{J}$, whereas in plasma we have $\Omega \propto \mathbf{J} \propto \mathbf{v}$. The other key difference is that the fluctuating gravitational potential in these variables (which we will introduce in Sec. IV) depends on both the canonical coordinates and momenta, i.e., $\delta\Phi(\mathbf{r}, t) = \delta\Phi(\theta, \mathbf{J}, t)$. As we will see, these differences make the physics of stellar systems in some ways richer dynamically than that of electrostatic plasmas, but also make explicit calculation more difficult.

### IV. PERTURBED ORBITS

Now that we know how to describe stellar orbits in smooth central potentials, we are ready to consider what happens if we perturb those orbits. In particular, in this section we discuss how stellar orbits in an axisymmetric, razor-thin (i.e., two-dimensional) disk are modified when we introduce a rigidly rotating nonaxisymmetric perturbation (e.g., a spiral arm pattern). The discussion will be familiar to plasma physicists who have studied, for example, the interaction of ions with toroidal Alfvén waves in tokamaks.$^{26,28}$

#### A. Mean field and perturbations

In order to construct angle-action variables, we need a system that is integrable. Since we know that orbits in axisymmetric potentials are nearly always integrable, we will first consider an axisymmetrized version of our (in general, nonaxisymmetric) system, and construct the angle-action variables associated with motion in that axisymmetrized system. We then treat nonaxisymmetries as perturbations.

To begin our construction, at a given time let the exact, in general, nonaxisymmetric, gravitational potential of the whole two-dimensional system (including external perturbations, if they are present) be $\Phi(\mathbf{r}) = \Phi(\phi, R)$. Then, the exact Hamiltonian of a test star orbiting in this system is

$$H(\mathbf{r}, \mathbf{v}) = \frac{1}{2}\mathbf{v}^2 + \Phi(\mathbf{r}). \qquad (42)$$

The azimuthal average of $H$ is

$$\bar{H}(R, \mathbf{v}) = \frac{1}{2}\mathbf{v}^2 + \bar{\Phi}(R), \qquad (43)$$

where for any function $q(\mathbf{r})$, we define

$$\bar{q}(R) = \frac{1}{2\pi} \int d\phi q(\mathbf{r}). \qquad (44)$$

We know that we can construct angle-action variables $(\theta, \mathbf{J})$ such that $\bar{H}$ only depends on $\mathbf{J}$. Given these angle-action variables, we can express any function on phase space in terms of them via $\mathbf{r} = \mathbf{r}(\theta, \mathbf{J})$ and $\mathbf{v} = \mathbf{v}(\theta, \mathbf{J})$. In particular, we can return to our original exact Hamiltonian $H$ and split it into a $\theta$-independent “mean-field” part and a $\theta$-dependent “fluctuation”:

$$H[\mathbf{r}(\theta, \mathbf{J}), \mathbf{v}(\theta, \mathbf{J})] = H_0(\mathbf{J}) + \delta H(\theta, \mathbf{J}), \qquad (45)$$

where the $\theta$ average of some quantity $q(\theta, \mathbf{J})$ is

$$q_0(\mathbf{J}) = \int \frac{d\theta}{(2\pi)^2} q(\theta, \mathbf{J}). \qquad (46)$$

Here, the “0” subscript reminds us that $q_0$ is equivalent to the zeroth Fourier mode of $q$ when expanded in angles [see Eq. (36) above].

One can show that the “mean-field Hamiltonian” $H_0$ and the “axisymmetric Hamiltonian” $\bar{H}$ are equivalent (even though this is not generally true of the potential alone, i.e., $\Phi_0 \neq \bar{\Phi}$). Furthermore, we have that

$$\delta H(\theta, \mathbf{J}) = \Phi(\mathbf{r}) - \bar{\Phi}(R) \equiv \delta\Phi(\theta, \mathbf{J}), \qquad (47)$$

i.e., the angle-dependent part of the Hamiltonian is just the nonaxisymmetric part of the potential. Hence, from now on, we write

$$H[\mathbf{r}(\theta, \mathbf{J}), \mathbf{v}(\theta, \mathbf{J})] = H_0(\mathbf{J}) + \delta\Phi(\theta, \mathbf{J}) \qquad (48)$$

and refer to $H_0$ and $\delta\Phi$ as the mean-field Hamiltonian and the potential perturbation, respectively.

#### B. Rotating potential perturbation

We now assume that the nonaxisymmetric potential perturbation $\delta\Phi$ rotates in the $\phi$ direction with fixed pattern speed $\Omega_p$,

$$\delta\Phi = \delta\Phi(\phi - \Omega_p t, R). \qquad (49)$$

Given that the perturbation is periodic in $\phi - \Omega_p t$, using Eq. (41a) and the fact that $R$ is a function only of $\theta_R$ and $\mathbf{J}$ (through the implicit relation 41b), we can always expand $\delta\Phi$ as a Fourier series [c.f. Eq. (36)]

$$\delta\Phi(\theta, \mathbf{J}, t) = \sum_{\mathbf{n} \neq \mathbf{0}} \delta\Phi_{\mathbf{n}}(\mathbf{J}) e^{i(\mathbf{n}\cdot\theta - n_{\phi}\Omega_p t)}, \qquad (50)$$
where $\mathbf{n} = (n_\phi, n_R)$. For a perturbation with fixed azimuthal harmonic number $m > 0$ (e.g., an $m$-armed spiral with a sinusoidal azimuthal profile), the only terms that survive in the expansion (50) are those with $n_\phi = \pm m$. More generally, bar/spiral perturbations can usually be built as a superposition of a small number of $m$-harmonic components.

### C. Resonances

A star's orbit evolves according to Hamilton's equations (32). Thus, using the Fourier expansion (50), its action evolves as $\mathbf{J} = \mathbf{J}(0) + \delta\mathbf{J}(t)$, where

$$ \frac{d\delta\mathbf{J}}{dt} = -i \sum_{\mathbf{n}} \mathbf{n} \delta\Phi_{\mathbf{n}}(\mathbf{J}) e^{i(\mathbf{n} \cdot \theta - n_{\phi} \Omega_{p} t)}. \quad (51) $$

Given this equation, it will not surprise the reader that the majority of interesting orbital modifications by rigidly rotating potential perturbations occur at resonances, i.e., locations $\mathbf{J}$ in phase space where

$$ \mathbf{N} \cdot \mathbf{\Omega}(\mathbf{J}) = N_{\phi} \Omega_{p}, \quad (52) $$

for some vector of integers $\mathbf{N} = (N_\phi, N_R)$. That is because these are locations where, substituting the unperturbed trajectories (33a) and (33b), the argument of the exponential on the right-hand side of Eq. (51) is constant.

For concreteness, suppose the perturbation in question has azimuthal harmonic number $m$. Then, the most important values of $\mathbf{N}$—which are sometimes referred to as the perturbation's major resonances—are

$$ \mathbf{N}_{\text{CR}} = (m, 0) \quad \text{corotation resonance (CR)}, \quad (53a) $$

$$ \mathbf{N}_{\text{ILR}} = (m, -1) \quad \text{inner Lindblad resonance (ILR)}, \quad (53b) $$

$$ \mathbf{N}_{\text{OLR}} = (m, 1) \quad \text{outer Lindblad resonance (OLR)}. \quad (53c) $$

These resonances tend to dominate the in-plane evolution of disk galaxies driven by large-scale perturbations like spirals and bars. Resonances with $|N_R| \ge 2$ also exist, but are usually subdominant because bars and spirals have Fourier spectra dominated by the lowest few $n_R$. On the other hand, resonances with $|N_R| \ge 2$ are very important for understanding, e.g., the dynamical friction felt by an orbiting point mass.45–49

The resonances (53) have a simple physical interpretation in the important case of near-circular orbits (see Sec. III C 4 and Appendix A). In that case, the corotation resonance is located at the orbital radius $R_{CR}$ such that an orbit with guiding radius $R_g = R_{CR}$ advances in azimuth at exactly the same rate as the perturbation, i.e., $\Omega_c(R_{CR}) = \Omega_p$, with $\Omega_c$ the azimuthal frequency for circular orbits. If we now imagine keeping that star's orbit near-circular but increasing $R_g$ slightly, it will begin to circulate at a frequency slightly smaller than $\Omega_p$. As we continue this process, eventually $R_g$ will be sufficiently large that the difference $\Omega_p - \Omega_c$ is equal to the frequency of radial epicycles $\kappa$ [Eq. (A7b)]; this is the outer Lindblad resonance. Had we performed a similar experiment but instead gradually made the circular orbit smaller than $R_g$, the star would begin to orbit faster than $\Omega_p$, until $\Omega_p - \Omega_c$ coincided with $-\kappa$; this is the inner Lindblad resonance. In summary, for epicyclic orbits, we can write the resonance condition as

$$ m(\Omega_p - \Omega_c) = n\kappa, \quad (54) $$

with $n = 0, 1, -1$ for CR, OLR, and ILR, respectively. For an illustration of these resonant orbits, see, e.g., Fig. 6.10 of BT08.

### 1. Example: Orbits perturbed by a spiral potential

To illustrate the morphology of perturbed orbits and the importance of resonances, consider the motion of test stars in the potential $\Phi = \bar{\Phi} + \delta\Phi$, where $\bar{\Phi}$ is the logarithmic halo potential

$$ \bar{\Phi}(R) = V_0^2 \ln(R/R_0), \quad (55) $$

and $R_0$ is an arbitrary scale length. This potential gives rise to a flat rotation curve, meaning that all stars on circular orbits ($J_R = 0$) have the same orbital speed: $v_\phi^2 = R \partial\bar{\Phi}/\partial R = V_0^2 = \text{cst}$. This is in fact a reasonably realistic model for many galactic disks.3 For $\delta\Phi$, we choose the potential of a steadily rotating logarithmic spiral

$$ \delta\Phi(\phi, R, t) = -\epsilon V_0^2 \cos[\alpha \ln(R/R_0) - m(\phi - \Omega_p t - \phi_p)], \quad (56a) $$

$$ \alpha \equiv m \cot p. \quad (56b) $$

This corresponds to a spiral with $m$ arms, pattern speed $\Omega_p$, and pitch angle $p$ (so $p \to 0$ is a very tightly wrapped trailing spiral, which gets less and less tightly wound as $p \to 90^\circ$—for much more detail, see Sec. 6.1 of BT08). The rotation period of this spiral is simply, $T_p = 2\pi/\Omega_p$. In Fig. 12, we plot contours of $\delta\Phi/(\epsilon V_0^2)$ for $m = 2$ and $p = 60^\circ$ (a rather loosely wrapped arm).

We now fix the pattern speed to be $\Omega_p = 0.154 V_0/R_0$, so that the corotation radius of circular orbits is $R_{CR} = 6.5 R_0$. We take the spiral strength to be $\epsilon = 0.01$, the pitch angle $p = 60^\circ$, and the initial phase $\phi_p = 0$. We initialized the orbits of nine test stars in this potential, putting each star on a circular orbit with speed $V_0$ and azimuthal location $\phi(t=0) = 0$. The initial conditions of the stars therefore differed only in their initial orbital radii $R(0)$, which we sampled linearly between $R_{ILR}$ and $R_{OLR}$ inclusive. We integrated each orbit for a total time $20T_p$.

In Figs. 13 and 14, we show the resulting orbits. In each panel, the black line shows the trajectory of a star whose initial radius $R(0)$ is given at the top of the panel. we also show the corotation radius of near-circular orbits $R_{CR}$ with a dashed red line, and the corresponding Lindblad radii $R_{ILR}/OLR$ with dotted red lines. In Fig. 13, we plot the trajectories in the inertial frame of the galaxy, while in Fig. 14, we plot them in a frame that rotates azimuthally at the rate $\Omega_p$ [so the spiral potential (56) is fixed].

These figures reveal that initially circular orbits are not dramatically distorted if the star sits far from a resonance. In this case, the orbits are well described with linear perturbation theory which we introduce in Sec. IV E. On the other hand, the orbits are changed qualitatively if they are located near a resonance. In particular, panel (a) shows an orbit that is “trapped” at the ILR, which is the key resonance responsible for the formation of stellar bars. Meanwhile, panel (e) shows an orbit that is trapped at corotation. When viewed in the rotating frame, these orbits do not complete any azimuthal loops, but are rather “stuck” on one side of the spiral, and librate back and forth in the trough of the spiral potential. Trapped orbits at corotation are especially important for determining the slowdown rate of spinning stellar bars (Sec. VIII B 1), the saturation of spiral instabilities (Sec. VIII B 2), and so on. (They are also very familiar to plasma physicists, who are used to particles “bouncing” in the troughs of waves if the
<figure><img src="image_25.png" /><figcaption>FIG. 12. Contours of the spiral potential, Eq. (56), for m = 2 and p = 60° (the azimuthal phase is arbitrary). The dotted and dashed white curves show the radial locations of the ILR, CR, and OLR, respectively, assuming Ωp = V0/(6.5R0), but this is only meant to guide the eye—the logarithmic spiral (56) and the logarithmic halo (55) are both scale-free, so nothing about the figure would change if we modified Ωp apart from the radii of these white circles.</figcaption></figure>

wave's phase velocity matches the particle velocity.) Finally, panel (i) shows an orbit trapped at the OLR.

The rest of this section is devoted to a mathematical description of how stellar orbits respond to a perturbation of the form (50). We start by discussing a quantity that is rigorously conserved for such orbits (Sec. IV D). We then describe how one can approximate the evolution of linearly perturbed orbits (Sec. IV E). Finally, we show how the nonlinearly distorted orbits that can develop near resonances can be treated analytically in the pendulum approximation (Sec. IV F).

### D. Jacobi integral

One thing we can immediately say about test star motion in the Hamiltonian (48) is that it preserves neither angular momentum L nor energy E (since δΦ is both nonaxisymmetric and time-dependent). Yet, it does conserve the Jacobi integral (see Sec. 3.3.2 of BT08):

$$E_J \equiv E - \Omega_p L. \qquad (57)$$

This follows because $H_J \equiv H - \Omega_p L$ is just the Hamiltonian in the frame co-rotating with the perturbation, and in that frame, the Hamiltonian is time-independent and hence conserved (see Sec. 3.3.2 of BT08 for proof). Since $E_J = \text{cst.}$, any interaction of a test star with a constant amplitude, rigidly rotating potential perturbation must produce changes to the star's energy $\Delta E$ and angular momentum $\Delta L$ that are related by

$$\Delta E = \Omega_p \Delta L. \qquad (58)$$

If $\Delta E$ is sufficiently small, we can approximate it by $\Delta E \approx (\partial E / \partial L) \Delta L + (\partial E / \partial J_R) \Delta J_R = \Omega_\phi \Delta L + \Omega_R \Delta J_R$, following the generic definition of orbital frequencies in Eq. (32a). This allows us to

relate the change in radial action $\Delta J_R$ to the change in angular momentum $\Delta L$,

$$\Delta J_R \approx \frac{\Omega_p - \Omega_\phi}{\Omega_R} \Delta L. \qquad (59)$$

For instance, stars that corotate with the perturbation, $\Omega_\phi = \Omega_p$, do not experience changes to their radial actions to first order.

### E. Linear perturbations

Suppose the perturbation strength $|\delta\Phi|/|H_0| \sim \epsilon \ll 1$ [just as in Eq. (56)]. Then, to zeroth order in $\epsilon$, we have $\delta\mathbf{J}^{(0)} = \mathbf{0}$, i.e., the star's trajectory is just a straight line through angle-action space, $\mathbf{J} = \mathbf{J}_0$ and $\theta = \theta_0 + \Omega t$. Plugging this back in to the right-hand side of Eq. (51) and integrating forward in time, we get

$$\delta \mathbf{J}^{(1)}(t) = -i \sum_{\mathbf{n}} \mathbf{n} \delta\Phi_{\mathbf{n}}(\mathbf{J}_0) e^{i \mathbf{n} \cdot \theta_0} \frac{e^{i(\mathbf{n} \cdot \Omega \mathbf{J}_0) - n_{\phi} \Omega_p} t - 1}{i(\mathbf{n} \cdot \Omega \mathbf{J}_0) - n_{\phi} \Omega_p}. \qquad (60)$$

Thus, to first order in $\epsilon$, the perturbation (50) nudges the star's orbital action by an amount $\delta\mathbf{J}^{(1)}$, which is proportional to $\epsilon$ in magnitude and oscillatory in time, except in the case where the star's mean field orbit is in resonance with the perturber, namely, when there exists some pair of integers $\mathbf{n} = \mathbf{N} = (N_\phi, N_R)$ such that (52) holds for $\mathbf{J} = \mathbf{J}_0$. In the resonant case, (60) predicts that $\delta\mathbf{J}$ will grow linearly in time, since the Nth Fourier component of the forcing (50) is constant in the frame moving with the unperturbed stellar orbit.

In reality, it is precisely near these resonant locations that linear perturbation theory will eventually break down. In that case, one should instead turn to an alternative, nonlinear description of the dynamics.

### F. The pendulum approximation

To describe the near-resonant orbits analytically, we turn to the “pendulum approximation.” The main idea of the pendulum approximation is that near resonances, the linear approximations used in Sec. IV E only truly break down along certain special “resonant directions” in angle-action space, namely, those directions where the argument of the exponential in Eq. (50) varies slowly for some resonance $\mathbf{n} = \mathbf{N}$. This allows us to average out the oscillatory motion in the “nonresonant” directions and thereby reduce the complicated resonant dynamics to the integrable problem of a one-dimensional pendulum.

For definiteness, let us set the number of dimensions to $d = 2$. Then, $\theta = (\theta_\phi, \theta_R)$ and $\mathbf{J} = (J_\phi, J_R)$. We now make a canonical transformation (see Sec. D.4 of BT08) to a new set of coordinates, which are again angle-action coordinates of the unperturbed problem.39,50 Precisely, we map $(\theta, \mathbf{J}) \to (\theta', \mathbf{J}')$, where $\theta' = (\theta_f, \theta_s)$ consists of the “fast” and “slow” angles

$$\theta_f \equiv \theta_R, \qquad (61a)$$

$$\theta_s \equiv \mathbf{N} \cdot \mathbf{\theta} - N_\phi \Omega_p t, \qquad (61b)$$

and $\mathbf{J}' = (J_f, J_s)$ consists of the corresponding fast and slow actions

$$J_f \equiv J_R - \frac{N_R}{N_\phi} L, \qquad (62a)$$

$$J_s \equiv \frac{L}{N_\phi}. \qquad (62b)$$
<figure><img src="image_26.png" /><figcaption>FIG. 13. Orbits of stars in the combined logarithmic halo (55) and spiral (56) potentials. All stars are initially on circular orbits with speed $V_0$ but differ in their initial orbital radii $R(0)$. The dotted and dashed lines show the Lindblad and corotation radii, as in Fig. 12. In particular, in panels (a) (e), and (i) the orbit is trapped at the ILR, CR, and OLR, respectively (see Fig. 14).</figcaption></figure>

(In plasma, the analog of, e.g., Eq. (61) would be a coordinate transform $x' = x - v_{ph}t$, where $v_{ph}$ is the phase velocity of a wave moving along the $x$ axis). Having made this transformation, we may rewrite the Hamiltonian $H$ [Eqs. (48) and (49)] in terms of the new coordinates:

$$H(\theta', \mathbf{J}') = H_0(\mathbf{J}') - N_{\phi} \Omega_{\mathrm{p}} J_s + \sum_{\mathbf{k} \neq 0} \Psi_{\mathbf{k}}(\mathbf{J}') e^{i\mathbf{k} \cdot \theta'}, \quad (63)$$

where $\mathbf{k} = (k_f, k_s)$ is a vector of integers and we have expanded $\delta\Phi$ as a Fourier series in the new angles $\theta'$, i.e., written $\delta\Phi(\mathbf{r}, t) = \sum_{\mathbf{k}} \Psi_{\mathbf{k}}(\mathbf{J}') e^{i\mathbf{k} \cdot \theta'}$. The coefficients $\Psi_{\mathbf{k}}$ are easily computed for a
<figure><img src="image_27.png" /><figcaption>FIG. 14. As in Fig. 13 except in a frame that rotates in the $\phi$ direction at rate $\Omega_p$ (so the spiral pattern is stationary).</figcaption></figure>

simple bar or spiral model.⁵¹,⁵² The special thing about the form (63) of the Hamiltonian is that it has no explicit time dependence [or rather, the time dependence has been absorbed into the definition of the angle $\theta_s$; see Eq. (61)].

The fast angles $\theta_f$ evolve on the orbital timescale, whereas $\theta_s$ evolves on the much longer timescale $\sim (\mathbf{N} \cdot \mathbf{\Omega} - N_{\phi}\Omega_p)^{-1}$. Thus, we choose to average $H$ over the unimportant fast motion and work instead with the simpler Hamiltonian $h \equiv (2\pi)^{-1} \int d\theta_f H$ at fixed $J_f$.

We also expand $h$ around the resonance. To do this, we let the resonant action be $J_{\text{res}}$ [meaning that Eq. (52) is satisfied for $J = J_{\text{res}}$]. In the transformed coordinates, this is $J'_{\text{res}} = (J_f, J_{s,\text{res}})$, and so we can define

$$I \equiv J_s - J_{s,\text{res}}, \qquad \varphi_k = \theta_s + \arg(\Psi_k)/k, \qquad (64)$$

with $\varphi_k \in [-\pi, \pi]$ (and we used the shorthand $\Psi_{(0,k)} = \Psi_k$). Expanding $h$ for small $I$, and assuming that a single Fourier component $k$ dominates the perturbation, we find that at each fixed $J_f$ the
dynamics reduces to motion in the “slow plane” ($\theta_s, J_s$) dictated by the pendulum Hamiltonian

$$h(\varphi, I) = \frac{1}{2}GI^2 - F \cos k\varphi, \qquad (65)$$

where $\varphi = \varphi_k$ and

$$G(J_f) = \left. \frac{\partial^2 H_0}{\partial J_s^2} \right|_{J_{s,\text{res}}}, \quad F(J_f) = -2|\Psi_k|_{J_{s,\text{res}}}. \qquad (66)$$

Here, $J_{s,\text{res}}(J_f)$ is the resonant value of the slow action at fixed fast action, $|\Psi_k|_{J_{s,\text{res}}}$ is the strength of the bar perturbation on resonance, and $\varphi \in [-\pi, \pi]$ is just a phase-shifted slow angle. Note that in contrast with linear perturbation theory (Sec. IV E), in this case if the perturbation strength is $|\delta\Phi|/|H_0| \sim \varepsilon \ll 1$, then the “perturbing” part of the Hamiltonian (65) has relative strength $|F/h| \sim \mathcal{O}(\varepsilon^{1/2})$, and the timescale for evolution due to this perturbation is therefore $t_{\text{lib}} \sim \varepsilon^{-1/2}t_{\text{dyn}}$. Appreciating this difference is crucial when it comes to developing an appropriate kinetic theory (see Sec. V E).

The variables $(\varphi, I)$ are canonical variables for the pendulum Hamiltonian (65). The pendulum moves at constant “energy” $h$, either on an untrapped “circulating” orbit with $h > F$ or a trapped “librating” orbit with $h < F$ (the separatrix between these two families is $h = F$). The libration period for oscillations around $(\varphi, I) = (0, 0)$ is

$$t_{\text{lib}} \equiv 2\pi/\omega_{\text{lib}} \quad \text{where} \quad \omega_{\text{lib}}(J_f) \equiv \sqrt{kFG}. \qquad (67)$$

The maximum width in $I$ of the librating “island” is at $\varphi = 0$, where it spans $I \in [-I_h, I_h]$ and $I_h$ is the island half-width:

$$I_h(J_f) \equiv 2\sqrt{F/G}, \qquad (68)$$

which is independent of $k$. The contours of this Hamiltonian are plotted in Fig. 15 for the case with $k = 1$.

The separatrix between librating and circulating orbit families is shown with a black dashed line.

## 1. Example: Corotation resonance, flat rotation curve

Let us illustrate the pendulum dynamics with a concrete example. This requires stipulation of the resonance $N$ under consideration. Let us consider the corotation resonance (CR), as defined in Eq. (53a), with the resonance location at the angular momentum value $L = L_{CR}(J_R)$ for a given radial action $J_R$. Then, we have

$$N_{\phi} = m, \quad N_{R} = 0, \quad k = 1, \qquad (69a)$$

$$J_s = L/m, \quad J_f = J_R, \quad I = (L - L_{CR})/m, \qquad (69b)$$

$$\varphi = m(\theta_{\phi} - \Omega_{\text{p}}t), \quad |\Psi_1| = |\delta\Phi_{(m,0)}|. \qquad (69c)$$

We must also stipulate a background potential; the simplest choice is the logarithmic halo potential (55), which gives rise to a flat rotation curve with circular speed $V_0$. With this choice, for initially circular orbits ($J_f = J_R = 0$), we get

$$L = Rv_{\phi} = RV_0, \qquad (70a)$$

$$\Omega_c(L) = V_0/R = V_0^2/L, \qquad (70b)$$

$$G = -N_{\phi}^2 V_0^2 / L_{\text{CR}}^2 = -N_{\phi}^2 / R_{\text{CR}}^2. \qquad (70c)$$

We now combine these two choices, i.e., we consider circular, corotating orbits in a logarithmic halo. In this case, the slow action $J_s$ is just proportional to the guiding radius $R$, with the resonant value $J_{s,\text{res}}$ occurring at radius $R_{CR}$, and the slow angle is (up to a phase shift) equivalent to $m[\varphi - \int_0^t dt' \Omega_p(t')]$. Thus, in this case, one can think of the transformation into slow and fast angle-action variables simply as “moving into the rotating frame” (or in the plasma context, shifting to a frame in which the wave’s phase velocity is zero).

With these choices, the libration time (67) can be estimated as

$$t_{\text{lib}} \approx \frac{t_{\text{dyn}}}{m \varepsilon^{1/2}} = \frac{2\pi}{m \varepsilon^{1/2} \Omega_{\text{p}}}, \qquad (71a)$$

$$\approx 1 \text{ Gyr} \times \left(\frac{m}{2}\right)^{-1} \left(\frac{\varepsilon}{0.01}\right)^{-1/2} \left(\frac{t_{\text{dyn}}}{200 \text{ Myr}}\right), \qquad (71b)$$

where $t_{\text{dyn}} = 2\pi R/V_0$ is the circular orbital period at radius $R$, and all quantities must be evaluated at the resonance location $R = R_{CR}$. Clearly, $t_{\text{lib}}$ is always several times longer than the orbital period $t_{\text{dyn}}$. Moreover, Eq. (68) lets us estimate the radial width of the resonant island for these orbits

$$R_h = \frac{mI_h}{V_0} \approx 2\varepsilon^{1/2} R_{CR}, \qquad (72a)$$

$$\approx 1.6 \text{ kpc} \times \left( \frac{\varepsilon}{0.01} \right)^{1/2} \left( \frac{R_{CR}}{6.5 \text{ kpc}} \right). \qquad (72b)$$

Realistic spirals like those in the Milky Way today have $\varepsilon \sim$ several percent or more.⁵³ Thus, the resonant island at corotation can span one or more kiloparsecs in radius. This must be borne in mind whenever one makes a “narrow resonance” approximation, since $R_h$ can

<figure><img src="image_28.png" /><figcaption>FIG. 15. Pendulum phase space. Colors show contours of constant $h(\varphi, I)$ [see Eq. (65)] for $k = 1$. The black dashed line shows the separatrix between librating (i.e., trapped) and circulating (i.e., untrapped) orbit families.</figcaption></figure>
also be comparable to the radial scale over which the galaxy's underlying density profile changes significantly (≈ 3 kpc for the stellar disk of the Milky Way).

As a more quantitative example, let us apply our pendulum formalism to the nine orbits shown in Fig. 14. In Fig. 16, we show these nine trajectories through the space of azimuthal angle $\phi - \Omega_p t + \ln(R_{CR}/R_0)\cot p$ and angular momentum $L$.

Up to constants, these are equivalent to the slow angle-action pair $(\phi, I)$ near corotation. With the colored lines, we plot contours of the pendulum Hamiltonian (c.f. Fig. 15) using the parameters for the corotation resonance [Eqs. (69) and (70)]. The dashed and dotted red lines again show the ILR, CR, and OLR from bottom to top, respectively. We see that the pendulum approximation does a good job of predicting the orbit shapes near corotation. The slight mismatch in the orbit shapes is mostly because the pendulum theory ignores all but one Fourier harmonic $k$ (as such, the agreement with theory gets better (respectively worse) if we use more open (respectively tight) spirals, i.e., take $p > 60^\circ$ [respectively $p < 60^\circ$]). The pendulum approximation can even work quite well *far away* from the CR. Mostly, this is because away from corotation, orbits are not perturbed much in angular momentum so the approximation $L = \text{const.}$ is a good one. However, we see that the pendulum approximation begins to break down close to the ILR and OLR, where resonant effects not associated with corotation start to play an important role in the dynamics. Near those resonances, angular momentum is clearly not conserved, and to describe the dynamics accurately, one must instead shift to the proper slow angle-action coordinates appropriate to each resonance.

## 2. Example: Radial migration

Radial migration is the name given to any dynamical process that takes a star on a circular or nearly circular orbit in a galactic disk and alters its angular momentum (or equivalently its guiding radius) while producing little or no change in its radial action. The ubiquity of radial migration in galactic disks, and the importance of this especially for the evolution of chemical gradients, was first appreciated by Sellwood and Binney.⁵⁴ Here, we perform two numerical experiments similar to the ones shown in their paper, in order to demonstrate the difference between linear (Sec. IV E) and nonlinear (Sec. IV F) perturbations.

We consider a similar setup to the one we described in Sec. IV C 1 by integrating an array of test particles in the logarithmic halo potential $\Phi$ with an accompanying spiral potential perturbation $\delta\Phi$, except this time we modify the potential perturbation (56) by multiplying it by a Gaussian in time

$$ \delta\Phi(\phi, R, t) = -\epsilon V_0^2 \cos[\alpha \ln(R/R_0) - m(\phi - \Omega_p t - \phi_p)] \times e^{-(t-t_{\text{peak}})^2/(2\tau^2)}. \quad (73) $$

This corresponds to a transient spiral which grows and decays on the timescale $\approx \tau$, and whose amplitude reaches its peak at $t = t_{\text{peak}}$. All other parameters are set to the same values as in Fig. 13. With these parameters, we can use Eq. (71a) to estimate the libration time at peak spiral amplitude to be $t_{\text{lib}} \approx 5T_p$. For the initial conditions of the stars, we use circular orbits with radii randomly sampled from the range $R(0)/R_0 \in [4, 9]$, which is centered on the corotation radius but does not reach as far as the ILR/OLR. We choose the initial azimuthal angles

<figure><img src="image_29.png" /><figcaption>FIG. 16. The same orbits as in Fig. 14, but plotted in the space of azimuthal angle (in the rotating frame) and angular momentum (relative to the corotation resonance value). Up to constants, these are equivalent to the slow angle and slow action at corotation [Eqs. (61) and (62)]. The background colored lines are contours of the associated pendulum Hamiltonian (c.f. Fig. 15). The dashed and dotted red lines have the same meaning as in Fig. 14.</figcaption></figure>

$\phi(0)$ randomly in $[0, 2\pi]$. We integrate for a total time of $10T_p$, and always set the time of peak amplitude to be halfway through the integration, $t_{\text{peak}} = 5T_p$.

First, we integrated the equations of motion of 5000 test stars in this spiral potential, choosing the growth/decay time to be $\tau = 3T_p$. In panels (a) and (b) of Fig. 17, we show the resulting evolution of the angular momenta of a randomly selected subset of these stars, shown with different colors, as a function of azimuthal angle and time. In panel (c), we show a scatterplot of the change in angular momentum, $\Delta L \equiv L(10T_p) - L(0)$, for every star, as well as the mean $\langle \Delta L \rangle$ at each fixed $L(0)$. The red dashed line in these panels shows the corotation resonance, while the solid pink lines are displaced from corotation by $\approx L_h$ [Eq. (68)], the approximate resonance width at peak spiral amplitude. We see that close to the resonance, the spiral traps stars and causes them to undergo large shifts in their angular momenta (up to the resonance width $\approx L_h$). As the spiral's amplitude decays, it “drops off” these stars at angular momenta, which are often rather different from their initial values. The blue line in panel (c) makes it clear that there is a preference for stars initially inside corotation ($L(0) < L_{CR}$) to receive positive $\Delta L$, and vice versa. To understand this, note that in the limit where the libration time at peak amplitude is very short, a star which was originally located inside corotation will librate back and forth across $L = L_{CR}$ several times before being deposited at some new $L = L(0) + \Delta L$. This new $L$ could be either inside or outside $L_{CR}$ with roughly equal probability, so on average the star will have moved out ($\Delta L > 0$). Meanwhile, far from the resonance, the stars respond roughly adiabatically to the perturbation and so $\Delta L$ is clustered around zero. We also mention that in panel (c),
<figure><img src="image_30.png" /><figcaption>FIG. 17. Radial migration in a galactic disk. We integrate the orbits of 5000 stars that are initially on circular orbits with random azimuths $\phi(0)$ and radii chosen randomly in the range $R(0)/R_0 \in [4, 9]$. The stars are perturbed by a transient spiral potential (73) with lifetime $\approx \tau = 3T_p$. Panels (a) and (b) show the evolution of the angular momenta of a handful of these stars as functions of azimuthal angle in the rotating frame and time, respectively. Panel (c) shows the change in angular momentum experienced by each particle by the end of the simulation as a function of its initial angular momentum $L(0)$. In each panel, the red dashed line corresponds to corotation, and the solid pink lines are offset from corotation by approximately the resonance half-width $L_h = ml_h$ [Eq. (68)] at peak spiral amplitude.</figcaption></figure>

many stars gather on the line with slope −2 that passes through $\Delta L = 0$ at corotation. This line corresponds to stars that are initially located at $L(0) = L_{CR} - \ell$ being deposited at $L_{CR} + \ell$.55,56

Second, in Fig. 18, we rerun the same experiment except with a much shorter-lived spiral, $\tau = 0.3 T_p$. We see from panel (a) that in this case, no orbits manage to undergo a full libration inside the trapping region. As a result, there is no preferential direction of the transport of stellar angular momenta across the resonance, so the mean $\Delta L$ in panel (c) is approximately zero for all initial $L(0)$. In fact, since the perturbation lasts for so short a time, the resonance location is not so significant, meaning the distribution of $\Delta L$ values does not depend particularly strongly on $L(0)$. In fact, the envelope of $\Delta L$ values is significantly broader far away from corotation than it was in Fig. 17 (so the variance $\langle (\Delta L)^2 \rangle$ is relatively large even though $\langle \Delta L \rangle$ is close to zero).

From the results of these numerical experiments, we can extract two key physical principles, which apply ubiquitously in the kinetic theory of both stellar systems and plasmas:

* The efficiency with which different perturbations shuffle the integrals of motion of a distribution of particles depends on the spatio-temporal correlation spectrum of those perturbations (just as in, e.g., the classic stochastic acceleration study of Sturrock$^{57}$). It is therefore imperative to develop a theory, which can account for the statistical evolution of a distribution of particles when exposed to such a perturbation spectrum, be it externally or self-consistently generated. This is what we will do in the rest of this tutorial article, starting in Sec. V.

* If a single coherent perturbation dominates the potential $\delta\Phi$ over a timescale $\gtrsim t_{lib}$, then it will trap resonant orbits nonlinearly, and this may drive the system down a substantially different evolutionary path compared to the case with no trapping. But if some process is able to interrupt the trapping process on a timescale $\ll t_{lib}$, then nonlinear trapping can be ignored and one can apply linear perturbation theory (Sec. IV E) with confidence. This “interruption” can be due to the rapid decay of the perturbation as in the example above, or because of a time-dependent $\Omega_p$ and

<figure><img src="image_31.png" /><figcaption>FIG. 18. As in Fig. 17, except the spiral grows and decays much more quickly, $\tau = 0.3 T_p$.</figcaption></figure>
the accompanying movement of the resonance location through phase space,⁵¹,⁵⁸ or it can be due to some collisional/diffusive process, which kicks particles in and out of the resonance region. We address the last of these possibilities in Sec. VIII B. We also mention that the interruption of nonlinear trapping by time-dependent wave amplitudes, collisionality, and so on is very commonly studied in plasma theory.⁵⁹⁻⁶²

## V. FUNDAMENTALS OF KINETIC THEORY

In Sec. II, we introduced the idea that an arbitrary gravitationally bound ensemble of stars reaches virial equilibrium after a few dynamical timescales $t_{dyn}$. It then moves through a series of such equilibria on the relaxation timescale $t_{relax} \gg t_{dyn}$. Then, in Secs. III and IV, we saw how to describe the quasiperiodic orbits of stars in smooth, regular, time-independent potentials $\Phi$, became familiar with the machinery of angle-action variables, and discussed the modification of orbits by some perturbation $\delta\Phi$. Here, we will set up the remainder of the article by bringing these two sets of ideas together. Our aim is to describe a quasi-equilibrium stellar system as a collection of quasiperiodic orbits and then to develop a quantitative theory of how fluctuations in the true gravitational potential cause that quasi-equilibrium to change slowly over time.

### A. Klimontovich description

We start with the microscopic Klimontovich (DF), or empirical DF,

$$f_m(\mathbf{r}, \mathbf{v}, t) = m \sum_{i=1}^{N} \delta(\mathbf{r} - \mathbf{r}_i(t)) \delta(\mathbf{v} - \mathbf{v}_i(t)). \quad (74)$$

(Here and in what follows, the subscript “M” stands for “microscopic.”) This function encodes the exact position $\mathbf{r}_i(t)$ and velocity $\mathbf{v}_i(t)$ of every star $i$ at time $t$, which are themselves the exact solutions to the Newtonian equations of motion (2). Therefore, for indistinguishable particles, the Klimontovich function (74) contains complete information about the system. Yet, at this stage, it is totally useless because if we knew all this information, we would not need kinetic theory. Nevertheless, describing the system in this way is what will allow us to most conveniently formulate that kinetic theory.

Note that $f$ is normalized such that $\int d\mathbf{r}d\mathbf{v}f = Nm = M$, the total mass of the stars. Thus, if there is no external forcing, the exact microscopic gravitational potential is

$$\Phi_m(\mathbf{r}, t) = \int d\mathbf{r}' d\mathbf{v}' \psi(\mathbf{r}, \mathbf{r}') f_m(\mathbf{r}', \mathbf{v}', t), \quad (75)$$

where $\psi(\mathbf{r}, \mathbf{r}') = -G/|\mathbf{r} - \mathbf{r}'|$. This is simply the formal solution of the Poisson equation

$$\nabla^2 \Phi_m = 4\pi G \rho_m, \quad (76)$$

where the microscopic mass density is $\rho_m = \int d\mathbf{v}f_m$. It follows that every particle $i$ satisfies Hamilton equations (30), namely,⁶³

$$\frac{d\mathbf{r}_i}{dt} = \frac{\partial H_m(\mathbf{r}_i, \mathbf{v}_i, t)}{\partial \mathbf{v}_i}, \quad (77a)$$

$$\frac{d\mathbf{v}_i}{dt} = - \frac{\partial H_m(\mathbf{r}_i, \mathbf{v}_i, t)}{\partial \mathbf{r}_i}, \quad (77b)$$

where $H_m$ is the microscopic one-particle Hamiltonian

$$H_m(\mathbf{r}, \mathbf{v}, t) = \frac{1}{2} \mathbf{v}^2 + \Phi_m(\mathbf{r}, t). \quad (78)$$

Using Eqs. (74)–(78), it is now easy to check that $f_m(\mathbf{r}, \mathbf{v}, t)$ exactly satisfies the Klimontovich equation

$$\frac{\mathrm{d}f_m}{\mathrm{d}t} = \frac{\partial f_m}{\partial t} + \mathbf{v} \cdot \frac{\partial f_m}{\partial \mathbf{r}} - \frac{\partial \Phi_m}{\partial \mathbf{r}} \cdot \frac{\partial f_m}{\partial \mathbf{v}} = 0, \quad (79)$$

which we notice can also be written in a pleasingly coordinate-free form as

$$\frac{\partial f_m}{\partial t} + [f_m, H_m] = 0. \quad (80)$$

This form of the equation is particularly convenient when we wish to perform a canonical coordinate transform, e.g., from position and velocity to angle-action variables. Physically, Eq. (80) just says that each infinitesimal “pixel” of phase space fluid is invariant if one follows it along the trajectory prescribed by the exact Hamiltonian $H_m$ (Liouville’s theorem).

### B. Angle-averaged and angle-dependent quantities

We now take the formal step of splitting all quantities into an angle-averaged part and an angle-dependent fluctuation (which at this stage is not necessarily small). Thus, we write the microscopic DF as

$$f_m = f_0 + \delta f, \quad (81)$$

where

$$f_0(\mathbf{J}, t) = \int \frac{d\mathbf{\theta}}{(2\pi)^d} f_m(\mathbf{\theta}, \mathbf{J}, t). \quad (82)$$

Obviously, if we were to convert back to real space coordinates, then $f_0$, for a galactic disk, would be axisymmetric (independent of $\phi$), while all nonaxisymmetric fluctuations would be contained in $\delta f$. The equivalent in a homogeneous plasma would be a DF $f_0$ that depended only on velocity, with all position-dependent fluctuations carried in $\delta f$. Similarly, the scalar potential $\Phi_m$ as defined in Eq. (75) is a linear functional of the DF $f_m$, as is the microscopic Hamiltonian $H_m$ defined in Eq. (78). Thus, in correspondence with (81), we may write

$$H_m = H_0 + \delta\Phi. \quad (83)$$

By construction, $H_0$ is the part of the Hamiltonian that only depends on actions $\mathbf{J}$, while any $\theta$ dependence (and hence any nonaxisymmetry of the potential) is contained in $\delta\Phi$ (see Sec. IV A). In a homogeneous system, $H_0 = v^2/2$. Since the transformation to angle-action coordinates (Sec. III C) is canonical, it conserves the infinitesimal volumes $d\mathbf{r}d\mathbf{v} = d\mathbf{\theta}d\mathbf{J}$ and so the Klimontovich DF (74) can be expressed as

$$f_m(\mathbf{\theta}, \mathbf{J}, t) = m \sum_{i=1}^{N} \delta(\mathbf{\theta} - \mathbf{\theta}_i(t)) \delta(\mathbf{J} - \mathbf{J}_i(t)). \quad (84)$$

We note the subtlety that even though all nonaxisymmetry of the DF $f_m$ is contained in $\delta f$, the quantity $\delta f$ may also contain an axisymmetric piece. Related to this is the fact that $f_0$ is not equivalent to the axisymmetric part of $f_m$. This means that, for example, one cannot in general construct the axisymmetric part of the exact surface density
distribution $\bar{\Sigma}_n(R) = (2\pi)^{-1} \int d\phi \int dv f_m$ of a disk from $f_0(\mathbf{J})$ alone: one also needs to know the $\theta_R$-dependent part of the DF.

Let us briefly consider what happens if there are no nonaxisymmetric fluctuations (which in reality would require $N \to \infty$). Setting $\delta f = \delta\Phi = 0$ in Eqs. (81)–(83) and plugging the results into (80), we get⁶⁴

$$ \frac{\partial f_0}{\partial t} + [f_0, H_0] = 0. \qquad (85) $$

The left-hand side of this equation can be thought of as a “mean field operator” which advects particles along their mean field trajectories; the equation therefore tells us that in the absence of fluctuations, the phase space density along mean field trajectories is conserved, as we would expect. In fact, since $f_0$ and $H_0$ only depend on $\mathbf{J}$, the Poisson bracket is zero [see Eq. (34)]:

$$ [f_0, H_0] = 0. \qquad (86) $$

This means that $\partial f_0 / \partial t = 0$, i.e., $f_0$ cannot change with time even at a fixed phase space location. Put differently, any $f_0$ that is only a function of $\mathbf{J}$ is an equilibrium distribution in the absence of fluctuations, or what is called a collisionless equilibrium, or equivalently a quasi-stationary distribution—a result known as Jeans’ theorem. In a homogeneous stellar system or electrostatic plasma, Eq. (86) is equivalent to the requirement that $[f_0, \frac{1}{2}\mathbf{v}^2] = 0$, which is trivially true of any $f_0$ which only depends on velocity. In this case and in the case of inhomogeneous stellar systems, the fact that $f_0$ satisfies (86) does not mean that it is a stable equilibrium; in general, determining whether some $f_0$ is linearly stable requires the linear response theory that we will develop in Sec. VI.

Let us now reinstate the fluctuations, in particular allowing $f_0$ to be time-dependent. Plugging the expansions (81), (83) into (80) and writing the Poisson bracket out explicitly as in Eq. (34), we get, without any approximations,

$$ \frac{\partial f_0}{\partial t} = \frac{\partial}{\partial \mathbf{J}} \cdot \left( \int \frac{d\theta}{(2\pi)^d} \delta f \frac{\partial \delta\Phi}{\partial \theta} \right), \qquad (87a) $$

$$ \frac{\partial \delta f}{\partial t} + \mathbf{\Omega} \cdot \frac{\partial \delta f}{\partial \theta} - \frac{\partial \delta \Phi}{\partial \theta} \frac{\partial f_0}{\partial \mathbf{J}} = - \frac{\partial}{\partial \theta} \left( \delta f \frac{\partial \delta \Phi}{\partial \mathbf{J}} \right) + \frac{\partial}{\partial \mathbf{J}} \left( \delta f \frac{\partial \delta \Phi}{\partial \theta} - \int \frac{d\theta}{(2\pi)^d} \delta f \frac{\partial \delta \Phi}{\partial \theta} \right), \qquad (87b) $$

where $\mathbf{\Omega}(\mathbf{J}) \equiv \partial H_0 / \partial \mathbf{J}$. Equation (87a) tells us how the (nonlinear) coupling of nonaxisymmetric fluctuations $\delta f$ and $\delta\Phi$ drives evolution of the angle-averaged part of the DF $f_0$. Note that this equation takes the form of a continuity equation in action space, as it must, since we are not allowing stars to be created or destroyed. Meanwhile, Eq. (87b) tells us how the fluctuating part of the DF evolves over time.

If we are doing the problem self-consistently, then we must also couple these two equations to the Poisson equation for the fluctuations, namely,

$$ \delta\Phi(\theta, \mathbf{J}, t) = \int d\theta' d\mathbf{J}' \psi(\theta, \mathbf{J}, \theta', \mathbf{J}') \delta f(\theta', \mathbf{J}', t). \qquad (88) $$

In plasma, the interaction kernel $\psi$ has no dependence on the “actions” (velocities), which makes the Poisson equation much easier to solve in practice.

### C. Fourier expansion

At this point, it is convenient to introduce a Fourier representation for our angle-dependent quantities. That is, we use the periodicity of the angle variables $\theta$ to express $\delta f$ and $\delta\Phi$ as Fourier series, like in Eq. (36), namely,

$$ \delta f(\theta, \mathbf{J}, t) = \sum_{\mathbf{n}} \delta f_{\mathbf{n}}(\mathbf{J}, t) e^{i\mathbf{n} \cdot \theta}, \qquad (89a) $$

$$ \delta\Phi(\theta, \mathbf{J}, t) = \sum_{\mathbf{n}} \delta\Phi_{\mathbf{n}}(\mathbf{J}, t) e^{i\mathbf{n}\cdot\theta}. \qquad (89b) $$

Let us now rewrite (87a) and (87b) in Fourier form so that

$$ \frac{\partial f_0}{\partial t} = - \frac{\partial}{\partial \mathbf{J}} \cdot \sum_{\mathbf{n}} \mathrm{in} \delta\Phi_{\mathbf{n}}^*(\mathbf{J}, t) \delta f_{\mathbf{n}}(\mathbf{J}, t), \qquad (90a) $$

$$ \frac{\partial \delta f_{\mathbf{n}}}{\partial t} + \mathrm{in} \cdot \mathbf{\Omega} \delta f_{\mathbf{n}} - \mathrm{in} \cdot \frac{\partial f_0}{\partial \mathbf{J}} \delta \Phi_{\mathbf{n}} = i \sum_{\mathbf{n}'} \mathbf{n}' \cdot \left[ \frac{\partial \delta f_{\mathbf{n}-\mathbf{n}'}}{\partial \mathbf{J}} \delta \Phi_{\mathbf{n}'} - \frac{\partial \delta \Phi_{\mathbf{n}-\mathbf{n}'}}{\partial \mathbf{J}} \delta f_{\mathbf{n}'} \right]. \qquad (90b) $$

Meanwhile, the Fourier-transformed version of Eq. (88) reads

$$ \delta\Phi_{\mathbf{n}}(\mathbf{J}, t) = (2\pi)^d \sum_{\mathbf{n}'} \int d\mathbf{J}' \psi_{nn'}(\mathbf{J}, \mathbf{J}') \delta f_{\mathbf{n}'}(\mathbf{J}', t). \qquad (91) $$

The (bare) coupling coefficients $\psi_{nn'}(\mathbf{J}, \mathbf{J}')$ encode, in essence, Poisson’s equation in angle-action space, through the decomposition

$$ \psi(\theta, \mathbf{J}, \theta', \mathbf{J}') = \sum_{\mathbf{n}, \mathbf{n}'} \psi_{nn'}(\mathbf{J}, \mathbf{J}') e^{i(\mathbf{n}\cdot\theta - \mathbf{n}'\cdot\theta')}. \qquad (92) $$

In the limit of an homogeneous multi-periodic plasma (Sec. II C 1), these coefficients are greatly simplified and become $\psi_{nn'} \propto \delta_{\mathbf{n}'}^n / |\mathbf{n}|^2$. In that limit, note that (i) the only “resonances” permitted are for $\mathbf{n} = \mathbf{n}'$ and (ii) the coupling between particles does not depend on the actions $(\mathbf{J}, \mathbf{J}')$. The fact that neither (i) or (ii) hold in general for realistic stellar systems makes linear response calculations much harder than in homogeneous plasmas (Sec. VI).

### D. Linear and quasilinear approximations

The only assumption we have made in deriving Eqs. (90a) and (91) is that $H_0$ is integrable, i.e., that there exists a set of global angle-action coordinates $(\theta, \mathbf{J})$ such that $H_0 = H_0(\mathbf{J})$. If this is true, then Eqs. (90a) and (91) are exact. Sadly, they are also very difficult to solve.

To render them solvable, we will often make an additional assumption that angle-dependent fluctuations are small, i.e., $\delta f \ll f_0$ and $\delta\Phi \ll H_0$. With this, it follows naturally from Eqs. (90a) and (90b) that the fluctuation $\delta f$ evolves on the dynamical timescale $t_{\text{dyn}} \sim 1/|\mathbf{\Omega}|$, whereas the angle-averaged DF $f_0$—whose rate of change is quadratic in the small fluctuations—evolves on some much longer timescale $t_{\text{relax}} \gg t_{\text{dyn}}$. This timescale separation is what allows the system to move through a series of virial equilibria, each of which approximately satisfies Eq. (85), as we described qualitatively in Sec. II C 3. For self-consistent problems, $H_0[f_0]$ will also be changing on the timescale $t_{\text{relax}}$, and in practice every time $H_0$ changes significantly, we have to recompute the mapping $(\mathbf{r}, \mathbf{v}) \to (\theta, \mathbf{J})$.

Our assumption of weak fluctuations means that we now have a well-defined program for solving the system (90a) and (91). On the
timescale ~ t$_{dyn}$, we may linearize equation (90b) by throwing away the small terms on the right-hand side

$$ \frac{\partial \delta f_n}{\partial t} + \mathbf{in} \cdot \mathbf{\Omega} \delta f_n - \mathbf{in} \cdot \frac{\partial f_0}{\partial \mathbf{J}} \delta \Phi_n = 0 \quad (93) $$

and then solve it under the assumption that $f_0$ is fixed. If we are doing a self-consistent problem, we must couple this linearized equation with the Poisson equation (91). The simultaneous solution of these two linear equations is a nontrivial problem and is the subject of Sec. VI. Supposing we can do this successfully, we can then plug the solutions $\delta f$, $\delta \Phi$ into the right-hand side of Eq. (90a)—or, if appropriate, its ensemble-averaged cousin (118)—and solve for the evolution of $f_0$ on the much longer timescale $t_{relax}$. Such evolution of $f_0$ is called quasilinear, since it involves the multiplication of two small quantities, each of which is determined using linear theory.

We emphasize that the linear equation (93) is just the kinetic version of the linear perturbation theory that we wrote down for the orbits of individual stars in Sec. IV E. To see this, note that the formal solution to Eq. (93) is just

$$ \delta f_n(\mathbf{J}, t) = \delta f_n(\mathbf{J}, 0) e^{-i \cdot \mathbf{\Omega} t} + \mathbf{in} \cdot \frac{\partial f_0}{\partial \mathbf{J}} \int_0^t dt' e^{-i \cdot \mathbf{\Omega} (t-t')} \delta \Phi_n(\mathbf{J}, t'). \quad (94) $$

The first term on the right-hand side represents the ballistic evolution (i.e., the mean-field motion, along straight lines in angle-action space) of some initial fluctuation $\delta f(\theta, \mathbf{J}, 0)$ as governed by the mean field Hamiltonian $H_0$, i.e., the phase mixing of unperturbed orbits $\mathbf{J} = \text{cst}$, $\theta = \mathbf{\Omega}t + \text{cst}$. Phase mixing causes the initial fluctuation to oscillate rapidly in action (or, more properly, frequency) space for times $t \gg 1/|\mathbf{\Omega}|$—see Fig. 19 for illustration. [Analogously in plasma, phase mixing consists of particles moving along $\mathbf{v} = \text{cst}$, $\mathbf{r} = \mathbf{v}t + \text{cst}$ and therefore $\delta f(\mathbf{r}, \mathbf{v}, 0)$ winds up to ever smaller scales in velocity space].

The second term on the right-hand side of Eq. (94) tells us how the DF is modified by the potential perturbation $\delta \Phi$ as evaluated along the unperturbed trajectory; in other words, it is the result of applying an equation like (60) to every star. Put differently, in any linear or quasilinear kinetic theory, the only particle trajectories that are accounted for are the ballistic unperturbed trajectories set by $H_0$ (straight lines in angle-action space) plus the small (sinusoidal) wiggles instigated by $\delta \Phi$. Quasilinear transport is then a consequence of resonances between

the oscillatory forcing $\delta \Phi$ and these first-order wiggles—hence enters perturbation theory at second order.

## E. Nonlinear and nonperturbative effects

One can imagine iterating the above procedure to higher orders in perturbation theory by plugging (94) back into the right hand side of (90b), calculating $\delta f_n$ to second order in the perturbation strength, and so on.⁶⁵ At the level of particle orbits, this would correspond to adding “sinusoidal disturbances on top of sinusoidal disturbances” over and over with ever smaller amplitude. In the case where the perturbations $\delta \Phi$ are stochastic, this approach is usually given the name weak turbulence theory. Weak turbulence is the subject of many a plasma monograph such as the excellent one by Kadomtsev⁶⁶ and underlies the resonance broadening theories pioneered by Dupree⁶⁷ and others.¹⁶,⁶⁸,⁶⁹ On the other hand, in the case where $\delta \Phi$ is deterministic, the second-order calculation leads, e.g., to the classic plasma echo phenomenon first pointed out by Gould et al.⁷⁰

However, no matter how many iterations of this perturbative scheme one performs, it cannot account for fundamentally nonperturbative effects such as particle trapping by strong potential perturbations (Sec. IV F). One way to see this is to note that the perturbative iteration scheme accounts for effects that are of order $\epsilon$, $\epsilon^2$, $\epsilon^3$, etc., where $\epsilon \sim |\delta \Phi/H_0|$ measures the strength of the potential perturbations. Nonlinear trapping, though, enters the dynamics at order $\epsilon^{1/2}$ (Sec. IV F), and so cannot be recovered by summing any finite number of perturbative terms. Another way to arrive at the same conclusion is to note that trapped orbits [like the one in Fig. 14(e)] are topologically distinct from untrapped ones [c.f. Fig. 14(f)], and no amount of iterative sinusoidal deformations of the latter can ever reproduce the former.

These fundamentally nonperturbative effects can be very important for the evolution of the system as a whole, but a kinetic theory that accounts for them cannot start from Eq. (93). In the case where $\delta \Phi$ is stochastic, this regime is termed strong turbulence, and formulating an accurate kinetic theory turns out to be very difficult even in the simplest plasmas.⁶⁸ Luckily, in galaxies, $\delta \Phi$ is often dominated by a single coherent perturbation, such as a strong bar or spiral arm. In this case, a simple kinetic theory incorporating trapping can be developed using the pendulum approximation (Sec. IV F), and this will be the subject of Sec. VIII. Moreover, just as we discussed at the end of Sec. IV F 2, if nonlinear effects are somehow quenched or interrupted on a short enough timescale, then the problem is relinearized and (94) is

<figure><img src="image_32.png" /><figcaption>FIG. 19. Schematic illustration of phase mixing in angle-action space as a function of time. The gray blob represents some phase space DF fluctuation $\delta f$ evolving in the absence of external forces [first term on the right-hand side of equation (94)]. Unperturbed trajectories follow horizontal lines in angle-action space [Fig. 10(b)]. Since the orbital frequencies $\mathbf{\Omega}(\mathbf{J})$ depend on $\mathbf{J}$, particles of different actions dephase. This leads to the appearance of ever finer structures in phase space.</figcaption></figure>
rendered (approximately) valid. We will see a concrete example of this in Sec. VIII B.

## VI. SELF-CONSISTENT LINEAR RESPONSE THEORY

In this section, we study the self-consistent response of a system to weak perturbations, using the linear approximation (Sec. V D). We assume throughout that the angle-averaged DF $f_0$ is time-independent, which means that we need to solve Eqs. (91) and (93) simultaneously, but we can ignore Eq. (90a).

For simplicity, we will develop the linear theory in the absence of any possible external perturbations, i.e., we will assume the only potential fluctuations are derived self-consistently internally from the angle-dependent part of the DF $\delta f$. However, the case where external perturbations are present is not essentially different (see, e.g., Chap. 5 of BT08), and we will show a relevant example in Sec. VI C 1.

### A. Laplace transforming the linearized Vlasov-Poisson system

Equations (91) and (93) constitute an initial value problem and are naturally solved via Laplace transform.⁷¹ Let us therefore define the Laplace transform of a function $h(t)$, and its inverse via

$$ \tilde{h}(\omega) = \int_{0}^{\infty} dt h(t) e^{i\omega t}, \qquad (95a) $$

$$ h(t) = \frac{1}{2\pi} \int_{i\sigma-\infty}^{i\sigma+\infty} d\omega \, \tilde{h}(\omega) e^{-i\omega t}. \qquad (95b) $$

Recall that the Laplace transform (95a) is defined for complex $\omega$ in the upper-half plane, or more precisely for $\mathrm{Im}\,\omega > \sigma$ where $\sigma$ is chosen large enough for the integral to converge at $t \to \infty$. In the inverse transform (95b), the “Bromwich contour” runs along the line $\mathrm{Im}\,\omega = \sigma$ from $\mathrm{Re}\,\omega = -\infty$ to $\mathrm{Re}\,\omega = +\infty$, but can be deformed at will, provided it never crosses over any singularities of $\tilde{h}(\omega)$.

Laplace transforming Eq. (93), we get

$$ -\delta f_n(\mathbf{J}, 0) - i(\omega - \mathbf{n} \cdot \mathbf{\Omega}) \delta \tilde{f}_n(\mathbf{J}, \omega) = \mathrm{in} \cdot \frac{\partial f_0}{\partial \mathbf{J}} \delta \tilde{\Phi}_n(\mathbf{J}, \omega), \qquad (96) $$

with $\delta f_n(\mathbf{J}, 0)$ is the initial fluctuation in the DF. The solution to this equation is

$$ \delta \tilde{f}_{\mathbf{n}}(\mathbf{J}, \omega) = - \frac{\mathbf{n} \cdot \frac{\partial f_0}{\partial \mathbf{J}}}{\omega - \mathbf{n} \cdot \mathbf{\Omega}} \delta \tilde{\Phi}_{\mathbf{n}}(\mathbf{J}, \omega) - \frac{\delta f_{\mathbf{n}}(\mathbf{J}, 0)}{i(\omega - \mathbf{n} \cdot \mathbf{\Omega})}. \qquad (97) $$

[Note that by taking the inverse Laplace transform of this, we can recover the formal solution (94)]. Meanwhile, the Laplace-transformed version of Eq. (91) is

$$ \delta \tilde{\Phi}_{\mathbf{n}}(\mathbf{J}, \omega) = (2\pi)^d \sum_{\mathbf{n}'} \int d\mathbf{J}' \psi_{\mathbf{n}\mathbf{n}'}(\mathbf{J}, \mathbf{J}') \delta \tilde{f}_{\mathbf{n}'}(\mathbf{J}', \omega). \qquad (98) $$

Let us now put Eqs. (97) and (98) together in order to eliminate $\delta \tilde{f}_{\mathbf{n}}(\mathbf{J}, \omega)$. The result is the following integral equation for the potential fluctuation $\delta \tilde{\Phi}_{\mathbf{n}}(\mathbf{J}, \omega)$:

$$ \begin{aligned} \delta \tilde{\Phi}_{\mathbf{n}}(\mathbf{J}, \omega) = &-(2\pi)^d \sum_{\mathbf{n}'} \int d\mathbf{J}' \frac{\mathbf{n}' \cdot \frac{\partial f_0}{\partial \mathbf{J}'}}{\omega - \mathbf{n}' \cdot \mathbf{\Omega}'} \psi_{\mathbf{n}\mathbf{n}'}(\mathbf{J}, \mathbf{J}') \delta \tilde{\Phi}_{\mathbf{n}'}(\mathbf{J}', \omega) \\ &- (2\pi)^d \sum_{\mathbf{n}'} \int d\mathbf{J}' \frac{\delta f_{\mathbf{n}}(\mathbf{J}', 0)}{i(\omega - \mathbf{n}' \cdot \mathbf{\Omega}')} \psi_{\mathbf{n}\mathbf{n}'}(\mathbf{J}, \mathbf{J}'), \end{aligned} \qquad (99) $$

where $\mathbf{\Omega}' \equiv \partial H_0(\mathbf{J}')/\partial \mathbf{J}'$.

### B. Dressed wakes and Landau modes

Equation (99) is a linear operator equation of the form

$$ |\delta \tilde{\Phi}(\omega)\rangle = \hat{M}(\omega)|\delta \tilde{\Phi}(\omega)\rangle + |\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle, \qquad (100) $$

where we have written the “state” of the potential fluctuations at frequency $\omega$ as $|\delta \tilde{\Phi}(\omega)\rangle$ and introduced the abstract linear operator $\hat{M}(\omega)$. Note that (100) is not merely a schematic representation of Eq. (99); instead, it is what is known in quantum statistical mechanics as a Dyson equation.²⁸,⁷² One can in fact develop the entire response theory of inhomogeneous systems on the basis of the abstract Dyson equation and operator theory⁷³⁻⁷⁶ without reference to the explicit angle-action representation. We will not go into the details here, but instead just discuss the physics that can be gleaned from Eq. (100).

On the right-hand side of Eq. (100), we have the source term $|\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle$, which corresponds to the second term on the right-hand side of Eq. (99). This is the part of the potential fluctuations that arises purely from the inevitable, finite-N discreteness noise $\delta f_n(\mathbf{J}, 0)$, propagated ballistically along mean field trajectories. (Had we included an external potential perturbation in our calculation, it would just be added to this source term.) The other term involves the operator $\hat{M}(\omega)$ acting on the potential fluctuation $|\delta \tilde{\Phi}(\omega)\rangle$. It tells us the extent to which potential fluctuations arise because of perturbations in the DF which were themselves induced by potential fluctuations. Thus we are discussing collective effects: perturbations to the DF that seed their own gravitational fields. Glancing back at (99), we expect these collective effects to be significant if there are fluctuations at frequencies $\omega$ which are close to some of the available stellar frequencies $\mathbf{n}' \cdot \mathbf{\Omega}'$. These fluctuations will be the all more efficient provided they come from stars at locations in action space $\mathbf{J}'$ where the mean field DF has a large gradient $\partial f_0 / \partial \mathbf{J}'$.

Another way to understand (100) is as follows. We start with an initial fluctuation $|\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle$ and then ask when we ask what is the self-consistent response of the system to itself. First, the fluctuation $|\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle$ provokes a response $\hat{M}(\omega)|\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle$. The response to this response is $\hat{M}(\omega)\hat{M}(\omega)|\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle$, and so on. Equating this string of self-consistent responses to the overall potential, we have $|\delta \tilde{\Phi}(\omega)\rangle = [1 + \hat{M}(\omega) + \hat{M}(\omega)\hat{M}(\omega) + \dots]|\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle = [1 - \hat{M}(\omega)]^{-1}|\delta \tilde{\Phi}_{\text{source}}(\omega)\rangle$, which we can multiply by $[1 - \hat{M}(\omega)]$ to arrive at (100).

How do we solve equation (100)? In homogeneous plasma theory, it is easy: in that case, the operator notation we have employed here is unnecessary because $\hat{M}(\omega)$ is a scalar (up to proportionality constants, it is equal to 1 minus the dielectric function $\epsilon_k(\omega)$), which “acts” on the potential fluctuation by simple multiplication. In inhomogeneous systems, we are not so lucky. Ultimately, this is because, as we warned in Sec. III D, the potential fluctuation $\delta \tilde{\Phi}(\theta, \mathbf{J})$ depends on both phase space coordinates and momenta, and so, $\hat{M}(\omega)$ acts upon these fluctuations in a nontrivial long-range fashion. In practice, this means we must solve the operator equation by casting $[1 - \hat{M}(\omega)]$ in a matrix representation.

The usual technique for doing this is the so-called biorthogonal basis method originally developed by Kalnajs.⁷⁷ The idea is that we can expand an arbitrary potential fluctuation in real space in terms of some set of basis functions $\Phi^{(p)}(\mathbf{r})$, and similarly for a density
fluctuation using some $\rho^{(p)}(\mathbf{r})$. For self-consistent problems, it is intel-
ligent to choose these functions in pairs such that they solve the
Poisson equation, $\nabla^2\Phi^{(p)} = 4\pi G\rho^{(p)}$, and are orthogonal, i.e.,
$$ \int d\mathbf{r} \, \Phi^{(p)*}(\mathbf{r}) \, \rho^{(q)}(\mathbf{r}) = 0 \text{ unless } p=q. $$
Many such sets are known, and the precise set one should use depends on the problem at hand (but they usually consist of familiar objects like spherical harmonics, Bessel functions, etc.). Physically, each discrete pair $(\Phi^{(p)}, \rho^{(p)})$ corresponds to a way in which the stellar system can oscillate self-consistently, and the full response of the system can be built up by superposing these discrete oscillations. This is reminiscent of elementary quantum mechanics, in which one builds the wavefunction of an electron in the hydrogen atom as a superposition of orthogonal eigenfunctions. Mathematically, once we have constructed such a basis, we can project our linear response equations onto it, thereby translating a problem of coupled differential equations into one of linear algebra.

The details of the biorthogonal basis method are rather technical,
and unimportant for most of what follows, so we relegate them to
Appendix B, but the explicit matrix form of the operator $\hat{M}(\omega)$ that
results is

$$
M_{pq}(\omega) = \frac{(2\pi)^d}{\mathcal{E}} \sum_{\mathbf{n}} \int d\mathbf{J} \frac{\mathbf{n} \cdot \partial f_0 / \partial \mathbf{J}}{\omega - \mathbf{n} \cdot \mathbf{\Omega}} \Phi_{\mathbf{n}}^{(p)*}(\mathbf{J}) \Phi_{\mathbf{n}}^{(q)}(\mathbf{J}), \quad (101)
$$

where $\Phi_{\mathbf{n}}^{(p)}(\mathbf{J})$ is the Fourier transform of $\Phi^{(p)}(\mathbf{r})$ in angle space [see Eq. (36)], and $\mathcal{E}$ is an arbitrary constant with units of energy. In Eq. (101), $M(\omega)$ is usually called the *response matrix*. There is a striking similarity between Eq. (101) and the corresponding quantity in plasma theory, and we explore this connection in Sec. VI B 1. Moreover, just as in plasma theory, although we derived (101) for $\mathrm{Im}\,\omega > 0$, we can analytically continue this function to the lower-half $\omega$ plane by taking the action space integral along the Landau contour. In practice, for any given $\mathbf{n}$, stellar systems typically sustain only a finite range of resonance frequencies, $\mathbf{n} \cdot \mathbf{\Omega}$. This leads to branch cuts in the response matrix,⁷⁸,⁷⁹ as is also the case, for example, in magnetized relativistic plasmas,⁸⁰ in the toroidal ion temperature gradient mode,⁸¹,⁸² or truncated Maxwellians.⁸³ We will not pursue any details here, but refer the reader to Chap. 5 of BT08 and to Fouvry and Prunet⁸⁴ for more information.

Having expressed $\hat{M}(\omega)$ as a matrix $M(\omega)$, the nature of the solu-
tion to Eq. (100) depends crucially on whether $\det[1 - M(\omega)] = 0$ or
$\det[1 - M(\omega)] \neq 0$:

* If $\omega$ is such that $\det[1 - M(\omega)] \neq 0$ then the inverse matrix $[1 - M(\omega)]^{-1}$ exists, and the solution to Eq. (100) is simply

$$
|\delta\bar{\Phi}(\omega)\rangle = [1 - M(\omega)]^{-1} |\delta\bar{\Phi}_{\text{source}}(\omega)\rangle. \quad (102)
$$

This is true for the great majority of $\omega$ values. Equation (102) equation tells us that at these frequencies, the noise source $|\delta\Phi_{\text{source}}(\omega)\rangle$ gets amplified by collective effects. Sometimes this is referred to as the noise being *dressed* by a gravitational *wake*. In the equivalent electrostatic plasma calculation, $1 - M$ is just the dielectric function $\epsilon_k(\omega)$; for near-equilibrium plasmas, the solution (102) tells us that each “bare” charge is locally screened by the Debye cloud of opposite charges that it gathers round itself. But in the gravitational context, there is no screening, and collective effects will often make gravitational perturbations on large scales *stronger*, not weaker. Loosely speaking, whereas Debye clouds tend to *reduce the effective charge*, gravitational

wakes tend to *increase the effective mass* hence typically shorten-
ing the relaxation time (see Sec. VI C 1 for an example). This
behavior was illustrated beautifully in the context of a local shear-
ing portion of a galactic disk by Toomre and Kalnajs.⁸⁵

* Alternatively, for $\omega$ such that $\det[1 - M(\omega)] = 0$, the inverse of $[1 - M(\omega)]$ does not exist, which tells us that $|\delta\Phi(\omega)|$ can be nonzero even in the absence of a source term on the right-hand side of Eq. (100). These are the so-called Landau modes of the system,⁸⁶ and they tend to occur for a discrete set of (generally) complex frequencies $\{\omega_m\}$. The system is then able to support a discrete set of collective oscillations which, just as electromagnetic waves in vacuum, do not need to be continually driven by some noise source. Moreover, under fairly weak assumptions one can show that the solution of Eq. (100) at this frequency will have time-dependence $\propto e^{i\omega_m t}$ (see Kalnajs⁷⁵ for a precise statement). Thus, a system is linearly unstable if there exists a Landau mode whose imaginary part $\mathrm{Im}\,\omega_m > 0$. Otherwise, the system is stable.

This confluence of a “dressed continuum” of fluctuations super-
posed with a discrete set of “Landau modes” is characteristic of both
plasmas and stellar systems, although the typical responses of these
systems to perturbations can be very different, as we will illustrate with
some examples below.

Note that although we use the term Landau mode above, these
are in general not the normal modes of the linearized system, because
$\omega$ is in general not real. An alternative approach to linear response the-
ory is to use purely real frequencies, in which case the fundamental
objects in the theory are (singular) van Kampen modes.⁸⁷ One can
show that Landau modes are equivalent to a superposition of a contin-
uum of these van Kampen modes. We will not refer to van Kampen
modes again in this article—we only note that they offer an alternative
way to formulate everything we present in this section.

1. Example: Jeans instability in a homogeneous periodic cube

First, let us address directly the analogy between stellar systems and homogeneous plasmas, by calculating the response matrix $M(\omega)$ for a fictitious three-dimensional homogeneous stellar system placed inside a periodic cube. As noted in Sec. III C 1, in that case, orbits are simple (multi-periodic) straight lines, and we can take $\theta \propto \mathbf{r}$ and $\mathbf{J} \propto \mathbf{\Omega} \propto \mathbf{v}$ [see Eq. (38)]. The system’s instantaneous potential and DF are linked by

$$
\Phi(\mathbf{r}) = \int d\mathbf{r}' d\mathbf{v}' f(\mathbf{r}', \mathbf{v}') \psi(\mathbf{r}, \mathbf{r}'), \quad (103)
$$

where $\psi(\mathbf{r}, \mathbf{r}') = -G/|\mathbf{r} - \mathbf{r}'|$ is the gravitational pairwise interaction.
Making use of the periodicity of this system one has

$$
\psi(\theta, \theta') = -\frac{G}{L\pi} \sum_{\mathbf{p}\neq0} \frac{e^{i\mathbf{p}\cdot(\theta-\theta')}}{|\mathbf{p}|^2}. \quad (104)
$$

Glancing at (B4), we find that, in homogeneous systems, the natural basis elements are $\Phi^{(p)} \propto e^{i\mathbf{p}\cdot\theta'}/|\mathbf{p}|$, namely, plane waves. This is unsurprising, since the system is translationally invariant. In Eq. (104), the factor $1/|\mathbf{p}|^2$ is naturally reminiscent of the $1/|\mathbf{k}|^2$ factor that arises when Fourier transforming the electrostatic Poisson equation for plasmas.
With this, the response matrix from Eq. (101) becomes

$$M_{pq}(\omega) = \frac{GL^2}{\pi} \delta_{\mathbf{p}}^{\mathbf{q}} \frac{1}{|\mathbf{p}|^2} \int d\mathbf{v} \frac{\mathbf{p} \cdot \partial f_0}{\bar{\omega} - \mathbf{p} \cdot \mathbf{v}}, \quad (105)$$

where $\bar{\omega} \equiv (2\pi)^{-1}L\omega$. Assuming⁸⁸ that the system's mean DF follows the Maxwellian distribution

$$f_0(\mathbf{v}) = \frac{\rho_0}{(2\pi\sigma^2)^{3/2}} e^{-|\mathbf{v}|^2/(2\sigma^2)}, \quad (106)$$

where $\rho_0$ is the system's mean density, and $\sigma$ is the velocity dispersion, we find from Eq. (105) that

$$M_{pq}(\omega) = \delta_{pq} \frac{1}{|\mathbf{p}|^2} \left(\frac{L}{L_I}\right)^2 [1 + \zeta Z(\zeta)], \quad (107)$$

where $\zeta = \bar{\omega}/(\sqrt{2}|\mathbf{p}|\sigma)$ is a dimensionless frequency and $Z$ is the usual plasma dispersion function⁸⁹

$$Z(\zeta) = \frac{1}{\sqrt{\pi}} \int_{-\infty}^{\infty} du \frac{e^{-u^2}}{u - \zeta}, \quad (108)$$

and $L_I$ is the so-called *Jeans length*

$$L_I = \sqrt{\frac{\pi\sigma^2}{G\rho_0}} \approx 10 \text{ kpc} \left( \frac{\sigma}{100 \text{ kms}^{-1}} \right) \left( \frac{\rho_0}{10^8 M_{\odot} \text{ kpc}^{-3}} \right)^{-1/2} . \quad (109)$$

The main difference the expression (107) for a self-gravitating system, and the analogous one for an electrostatic plasma, is an overall minus sign.

The dispersion relation for Landau modes is then $\delta_{\mathbf{p}}^{\mathbf{q}} - M_{pq}(\omega_m) = 0$. By looking for solutions with $\mathrm{Re}\,\omega_m = 0$, one can show that instability occurs ($\mathrm{Im}\,\omega_m > 0$) if the system is larger than a critical scale, $L > L_I$ (see Sec. 5.2.4 of BT08 for details). Thus, unlike in a plasma, a homogeneous Maxwellian distribution of stars is not necessarily stable to linear perturbations, but instead can be made unstable, e.g., by making the system dynamically colder (reducing $\sigma$) or denser (increasing $\rho_0$). The physics underlying this *Jeans instability* is the competition between gravity and velocity dispersion, which we sketch in Fig. 20.

In a stellar system without any velocity dispersion, a density fluctuation $\rho_0$ over a length scale $L$ would collapse under its own gravity on the timescale $\sim 1/\sqrt{G\bar{\rho}_0}$. However, if the stars have some velocity dispersion $\sigma$, then they can escape the region of size $L$ on the timescale $\sim L/\sigma$. If the latter timescale is longer than the former, the system cannot avoid collapse: this happens if $L \gtrsim \sigma/\sqrt{G\bar{\rho}_0}$, which, up to a prefactor, is the *Jeans length* (109). The *Jeans instability* is the basic mechanism underlying large-scale structure formation in the Universe.³

## 1. Example: Spiral Landau modes in cold stellar disks

Another important example of a Landau mode is a spiral mode in a disk galaxy. In truth, it is still not completely clear what gives rise to the spiral structure that we see in real galaxies like that in the top panel of Fig. 1.

Rather than attempt a review of this fascinating subject here,⁹²,⁹³ let us restrict ourselves to considering a more artificial system, namely,

an isolated, two-dimensional stellar disk consisting of equal-mass stars and no gas. In this case, spirals can arise manifestly as linear instabilities (i.e., Landau modes with $\mathrm{Im}\,\omega_m > 0$) of an underlying DF $f_0(L, J_R)$.

A classic example of such a spiral mode was discovered by Zang,⁹⁰ which we reproduce in Fig. 21 (adapted from Petersen et al.⁹¹) The model Zang employed was a “Mestel disk”—a simple, scale-free, two-dimensional disk embedded in a rigid dark matter halo—to which he added inner and outer radial truncations. Just as with the Jeans instability (Sec. VI B 1), one can always make a disk unstable by making it sufficiently dynamically cold and/or dense. In Zang’s case, this was achieved by modifying the *active fraction* of the disk, essentially the ratio between the self-gravitating mass in the disk and that in the unresponsive halo. Panel (a) of Fig. 21 shows the contours of $\mathrm{det}[l - M(\omega)]$ [Eq. (101)] in the complex frequency plane. The yellow cross shows the location where this determinant equals zero, which gives the frequency $\omega_m$ of an unstable spiral Landau mode. Precisely, the mode rotates with pattern speed $0.88 \, \Omega_0$ and growth rate $0.12 \, \Omega_0$, where $\Omega_0$ is a typical dynamical frequency comparable to $\sim 2\pi/t_{dyn}$. The spatial structure of the spiral mode is shown in panel (b); this is calculated by finding the eigenvector $\mathbf{X}_m$ with components $X_m^p$, such that $M(\omega_m) \cdot \mathbf{X}_m = 0$ and then writing

$$\delta\Phi_m(\mathbf{r}, t) \propto \mathrm{Re} \left[ \sum_{\rho} X_m^\rho \Phi^{(\rho)}(\mathbf{r}) e^{-i\omega_m t} \right]. \quad (110)$$

(Since this is a linear eigenmode, the overall normalization is arbitrary.)

The two examples we have given so far both involved smooth mass distributions that were made sufficiently cold and/or dense that they became unstable. However, there are other ways of producing instabilities which are perhaps more familiar to plasma physicists used to the bump-on-tail paradigm.¹⁶ For instance, De Rijcke et al.⁹⁴ recently performed a linear stability analysis, like the one we have outlined in this section, for a *grooved* Mestel disk. That is, they took a Mestel disk with a stable DF—the $J_R = 0$ slice through which is shown with a dotted line in Fig. 22(a)—and then carved a “groove” in it by removing a portion of the stars in a given range of angular momentum [solid line in Fig. 22(a)].

This grooved disk is also unstable to spiral Landau modes, and Fig. 22(b) shows the shape of the dominant mode.

<figure><img src="image_33.png" /><figcaption>FIG. 20. *Jeans instability* as a competition between velocity dispersion and gravitational attraction. A density fluctuation $\rho_0$ over a length scale $L$ will collapse under its own gravity if the collapse timescale $\sim 1/\sqrt{G\bar{\rho}_0}$ is shorter than the crossing time $L/\sigma$ where $\sigma$ is the velocity dispersion.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>FIG. 21. Determination of an unstable spiral Landau mode in a truncated Mestel disk following Ref. 90. Panel (a) shows the determinant of $I - M(\omega)$ in the complex frequency plane; the yellow cross shows where this function passes through zero. Panel (b) shows contours of the corresponding spiral Landau mode. Figure adapted from Petersen et al., Mon. Not. R. Astron. Soc. 530, 4378 (2024). Copyright 2024 Royal Astronomical Society.⁹¹</figcaption></figure>

Simulations show that the Landau mode spectrum of stellar disks is often dominated by one, or at most a few, isolated frequencies.⁹⁵ This means that spirals in unstable disks are perhaps more analogous to, e.g., Alfvénic eigenmodes in tokamak plasmas than to classic bump-on-tail instabilities in homogeneous plasmas, since the latter tend to harbor a whole spectrum of unstable modes (at least one at every unstable wavenumber $k$). Because of this, the evolution of unstable stellar disks might be more aptly described by something like the single-wave kinetics of O’Neil,⁹⁶ Berk and Breizman,⁹⁷,⁹⁸ etc., than by the typical many-wave quasilinear theory of Drummond and Pines,⁹⁹ and Kaufman.¹⁰⁰ See Sec. VIII for more.

### C. Explicit wake solution

We now suppose that our system is stable, meaning that all Landau modes are damped. An astrophysical example of such a (weakly) damped mode is the (dipole) $\ell = 1$ mode of globular clusters,¹⁰¹ which represents long-lasting sloshing oscillations.¹⁰² We further assume that we may wait long enough that these modes no longer contribute to the response (i.e., we wait for $t \gg |\mathrm{Im}\,\omega_m|^{-1}$, where $\omega_m$ is the system’s least strongly damped Landau mode). In this regime, we can solve (99) directly under the assumption that $\omega$ is not associated with a Landau mode. Having constructed the solution, we will

then perform an inverse Laplace transform to get the dressed wake part of the response in the time domain.

The easiest—and most physically enlightening—way to solve (99) is to guess what form the solution should have, and then work out the exact details post hoc. So, let us look at Eq. (99). Supposing for a moment that we were to ignore collective effects ($\bar{M} = 0$), we would find

$$[\delta\tilde{\Phi}_n(\mathbf{J}, \omega)]_{\text{bare}} = -(2\pi)^d \sum_{\mathbf{n}'} \int d\mathbf{J}' \frac{\delta f_{\mathbf{n}'}(\mathbf{J}', 0)}{i(\omega - \mathbf{n}' \cdot \mathbf{\Omega}')^2} \psi_{nn'}(\mathbf{J}, \mathbf{J}'), \quad (111)$$

which we call the “bare” solution. Inverse Laplace transforming this according to Eq. (95b), deforming the Bromwich contour so that it encircles the pole at $\omega = \mathbf{n}' \cdot \mathbf{\Omega}'$ (see Fig. 1.1 of Fouvry and Bar-Or¹⁰³ for illustration) and closing the contour with a large semicircle in the lower-half plane, we would find

$$[\delta\Phi_n(\mathbf{J}, t)]_{\text{bare}} = (2\pi)^d \sum_{\mathbf{n}'} \int d\mathbf{J}' \delta f_{\mathbf{n}'}(\mathbf{J}', 0) e^{-i\mathbf{n}' \cdot \mathbf{\Omega}' t} \psi_{nn'}(\mathbf{J}, \mathbf{J}'). \quad (112)$$

The physical interpretation of this result is very clear: it is the (Fourier transform of the) potential at location $\mathbf{r}(\theta, \mathbf{J})$ that is created by summing the contributions from every star with initial locations $\mathbf{r}'(0) = \mathbf{r}'(\theta'(0), \mathbf{J}'(0))$, assuming those stars simply travel along mean field trajectories $\theta' = \mathbf{\Omega}'t + \theta'(0)$ and $\mathbf{J}' = \mathbf{J}'(0)$.

<figure><img src="image_35.png" /><figcaption>FIG. 22. Groove-induced instability in a Mestel disk. Panel (a) shows the $J_R = 0$ cut through the DF in the stable (ungrooved) and unstable (grooved) case. Panel (b) shows the unstable spiral mode as calculated from linear stability analysis of the grooved disk. The solid circle shows the location of the groove for circular orbits, while the three dashed circles show the ILR, CR, and OLR resonance locations of the resulting spiral mode. Figure adapted from De Rijcke et al., Mon. Not. R. Astron. Soc. 484, 3198 (2019). Copyright 2019 Royal Astronomical Society.⁹⁴</figcaption></figure>
The ansatz we now make is that the solution to the full Eq. (99) is exactly of this functional form, except with the bare interaction potential $\psi_{nn'}(\mathbf{J}, \mathbf{J}')$ replaced with some dressed (and frequency-dependent) interaction potential $\psi_{nn'}^d(\mathbf{J}, \mathbf{J}', \omega)$, whose form we do not yet know. That is, we assume the solution has the form

$$ \delta\tilde{\Phi}_{\mathbf{n}}(\mathbf{J}, \omega) = -(2\pi)^d \sum_{\mathbf{n}'} \int d\mathbf{J}' \frac{\delta f_{\mathbf{n}'}(\mathbf{J}', 0)}{i(\omega - \mathbf{n}' \cdot \mathbf{\Omega}')'} \psi_{\mathbf{n}\mathbf{n}'}^d(\mathbf{J}, \mathbf{J}', \omega). \quad (113) $$

This ansatz is actually very natural if one is acquainted with Rostoker's principle in plasma kinetics, which essentially says that one can consider a hot plasma as a set of completely uncorrelated particles, provided one replaces the bare Coulombic potential of each particle with its dressed, or screened, version.¹⁶ (We will use this idea to derive the Balescu–Lenard kinetic equation for stable stellar systems in Sec. VII B.) Plugging the ansatz (113) into (99), we can read off an implicit equation for $\psi^d$:

$$ \psi_{nn'}^d(\mathbf{J}, \mathbf{J}', \omega) = \psi_{nn'}(\mathbf{J}, \mathbf{J}') - (2\pi)^d \sum_{\mathbf{n}''} \int d\mathbf{J}'' \frac{\mathbf{n}'' \cdot \partial f_0 / \partial \mathbf{J}''}{\omega - \mathbf{n}'' \cdot \mathbf{\Omega}''} \\ \times \psi_{nn''}(\mathbf{J}, \mathbf{J}'') \psi_{\mathbf{n}''\mathbf{n}'}^d(\mathbf{J}'', \mathbf{J}', \omega), \quad (114) $$

where $\mathbf{\Omega}'' \equiv \partial H_0(\mathbf{J}'') / \partial \mathbf{J}''$. Thus, we will have constructed a valid solution to Eq. (99) if we can solve Eq. (114) for $\psi^d$. While at first sight it may not seem like we have made much progress with this rewriting, it turns out that (114) is actually quite easy to solve if we project the various terms onto the biorthogonal basis elements $\Phi_{\mathbf{n}}^{(p)}(\mathbf{J})$, just as in Eq. (101). We leave the details to Appendix B, and here just write down the final answer:

$$ \psi_{nn'}^d(\mathbf{J}, \mathbf{J}', \omega) = -\frac{1}{\mathcal{E}} \sum_{p,q} \Phi_{\mathbf{n}}^{(p)}(\mathbf{J}) [1 - \mathbf{M}(\omega)]_{pq}^{-1} \Phi_{\mathbf{n}'}^{(q)*}(\mathbf{J}'). \quad (115) $$

We can usefully compare this to the bare interaction, which as we show in Appendix B can be expanded as

$$ \psi_{nn'}(\mathbf{J}, \mathbf{J}') = -\frac{1}{\mathcal{E}} \sum_{p} \Phi_{\mathbf{n}}^{(p)}(\mathbf{J}) \Phi_{\mathbf{n}'}^{(p)*}(\mathbf{J}'). \quad (116) $$

Thus, we see that in passing from $\psi \to \psi^d$, we have simply “divided” the bare interaction by $[1 - \mathbf{M}(\omega)]$. This is analogous to dividing by the dielectric function $\epsilon_k(\omega)$ in electrostatic plasma theory.

We can now inverse Laplace transform Eq. (113). The assumption that Landau modes can be ignored (i.e., $t \gg |\mathrm{Im}\,\omega_m|^{-1}$) means that when we deform the Bromwich contour we need only pick up contributions from the poles at $\omega = \mathbf{n}' \cdot \mathbf{\Omega}'$, and can ignore any singularities of $\psi^d$ [i.e., zeroes of $[1 - \mathbf{M}(\omega)]$, see Eq. (115)]. The result is

$$ \delta\Phi_{\mathbf{n}}(\mathbf{J}, t) = (2\pi)^d \sum_{\mathbf{n}'} d\mathbf{J}' \delta f_{\mathbf{n}'}(\mathbf{J}', 0) e^{-i\mathbf{n}' \cdot \mathbf{\Omega}'t} \psi_{\mathbf{n}\mathbf{n}'}^d(\mathbf{J}, \mathbf{J}', \mathbf{n}' \cdot \mathbf{\Omega}'), \quad (117) $$

which is the same as (112) with $\psi \to \psi^d$. Thus, we have realized Rostoker's idea: the dressed potential fluctuation is equivalent to that carried by the initial fluctuation ($\delta f_{\mathbf{n}'}(\mathbf{J}', 0)$), transported along the mean field ($\propto e^{-i\mathbf{n}' \cdot \mathbf{\Omega}'t}$), but with the Newtonian potential of each star replaced with the dressed potential ($\psi^d$) evaluated at the mean field frequency ($\mathbf{n}' \cdot \mathbf{\Omega}'$), with contributions summed over all stars ($\sum_{\mathbf{n}'} \int d\mathbf{J}'$). Furthermore, (117) can be plugged into the right hand

side of Eq. (94) to give an explicit solution for the fluctuations in the DF $\delta f_{\mathbf{n}}(\mathbf{J}, t)$. Note that the basis functions $\Phi^{(p)}$ play explicit no role in Eq. (117); they are merely a calculational tool which is only used to construct the dressed interaction $\psi_{nn'}^d(\mathbf{J}, \mathbf{J}', \omega)$ for a given mean field $f_0(\mathbf{J})$.

## 1. Example: Response of a stellar system to a perturbing mass

So far in this section we have concentrated on the response of a stellar system to internal noise $\delta f(t = 0)$. However, an entirely analogous calculation shows that even externally imposed potential fluctuations will be “dressed” by the system’s self-gravity—just as an externally imposed point charge will gather its own Debye cloud in an electrostatic plasma. Here, we will not pursue any quantitative details, but instead show a few figures that capture the characteristic features of such a response.

First, we present Fig. 23, which is adapted from Gibbon.¹⁰⁴ This is a cartoon of a quasineutral electrostatic plasma made of blue positive charges and red negative charges. When a pair of electrical probes is inserted into the plasma, they gather around them Debye clouds consisting of particles of the opposite charge. In this way, the electrostatic potential of the probes is screened over a Debye length. This cloud is typically much smaller than the scale of the whole system.

Next, we turn to Fig. 24, which is adapted from Magorrian.¹⁰⁵ This is an analogous stellar-dynamical scenario, in which a heavy point mass is driven on a straight line orbit at constant speed through a homogeneous stellar system with Maxwellian $f_0$ (Sec. VI B 1). The size of the box in this case has size $L \approx 0.9L_0$, so the system is (weakly) stable (recall that instability occurs if $L > L_0$). The contours in panel (a) show the density response of the system according to linear theory in the case that stars only respond to the point-mass perturber but not to each other. This we call the bare response. By contrast, panel (b) shows the response density in the dressed case in which stars interact both with the perturber and with each other [the contour levels are the same as in panel (a)]. Let us make two simple observations which contrast the behavior of the stellar system shown in these panels with the plasma system sketched in Fig. 23. (i) The overdense wake in the gravitational system acts to reinforce the perturbation rather than shield it, since there are no opposite charges. Heuristically therefore, one can think of gravitational polarization as a kind of “anti-Debye-shielding.” (ii) The scale of the induced wake, especially in the dressed case, can be comparable to the scale of the system itself, in contrast with an electrostatic plasma where the Debye length is usually much smaller than the system.

Since a homogeneous box is not a particularly realistic model for a stellar system, we finally present Fig. 25, which is adapted from D’Onghia et al.¹⁰⁶ This figure shows the response density of a stellar disk to an embedded point-mass perturber, driven on a circular orbit with constant speed for one azimuthal period (250 Myr). The disk is chosen to be fairly dynamically cold, but not so cold as to be unstable. The response was calculated using a full N-body simulation in which stars interact both with the perturber and with each other. The perturbing mass was much less than 1% of the whole system’s mass. Despite its relatively small mass, the perturber provokes a strongly amplified response in the disk surface density which extends over a large azimuthal range. This phenomenon, which was first discovered
<figure><img src="image_36.png" /><figcaption>FIG. 23. Cartoon [adapted from P. Gibbon, arXiv, 2007.04783 (2020). Copyright 2020 Author.¹⁰⁴] of the Debye shielding that occurs when electrical probes are placed inside a quasineutral plasma.</figcaption></figure>

by Julian and Toomre¹⁰⁷ using linear theory, reflects the fact that cold stellar disks are very effective amplifiers of perturbations, as we discussed in Sec. II C 3. As a result, it is typical of weakly stable stellar systems to develop collective responses on a global scale.

## VII. QUASILINEAR EVOLUTION AND THE BALESCU-LENARD EQUATION

In many cases, we do not know, nor do we really care, about the details of the microscopic DF $f_m$ of a particular stellar system. What we often do care about, however, is the structure and evolution of some suitably *averaged* DF, $\langle f_m \rangle$. If we are interested in the evolution of a system driven by random fluctuations, be they internally or externally generated, then the appropriate average is an *ensemble average* over all possible microscopic realizations of those fluctuations. The central idea of kinetic theory is to use the statistics of the realizations of the fluctuations to derive an equation for the evolution of $\langle f_m \rangle$.

We expect for an ensemble average of $f_m$ over microscopic realizations to be equivalent to an average over angles and over time (on timescales $\ll t_{relax}$). Thus, $\langle f_m \rangle = f_0$ and the relevant equation describing the evolution is [c.f. Eq. (90a)]

$$ \frac{\partial f_0}{\partial t} = - \frac{\partial}{\partial \mathbf{J}} \cdot \sum_{\mathbf{n}} \text{in} \langle \delta\Phi_{\mathbf{n}}^*(\mathbf{J}, t) \delta f_{\mathbf{n}}(\mathbf{J}, t) \rangle. \quad (118) $$

In Sec. V D, we spelled out the basic idea of *quasilinear* kinetic theory, which is that if the fluctuations $\delta\Phi$ and $\delta f$ are sufficiently small, then we can predict the evolution of $f_0$ on the timescale $t_{relax} \gg t_{dyn}$ by plugging the solutions from linear response theory (Sec. VI) into the right-hand side of Eq. (118). Quasilinear theory works provided such a timescale separation exists, and provided we can ignore higher order nonlinear and nonperturbative effects like trapping, as discussed in Sec. V E.

In this section, we will develop quasilinear kinetic theory more quantitatively. Our study of internally driven evolution will culminate in the Balescu–Lenard equation which describes the slow evolution of an isolated, stable N-star system. We first describe the basic physics underlying the Balescu–Lenard equation, then simply write it down without a proper derivation and discuss some of its properties (Sec. VII A). Then, we offer a simple derivation of the Balescu–Lenard equation based on Rostoker’s principle (Sec. VII B). We give several examples of this equation being used in practice (Sec. VII C), before discussing where the theory breaks down and how it might be improved (Sec. VII D).

As in Sec. VI, we consider systems in which the fluctuations are internally generated by the finite-N noise, but we note that a very similar formalism can be developed to account for externally driven quasilinear evolution.¹⁰³,¹⁰⁸

### A. A tale of two fluxes

As a first step, we take the formal solution for the perturbation to the DF in linear theory, Eq. (94), and plug it into the right hand side of Eq. (118). The result is

$$ \frac{\partial f_0}{\partial t} = - \frac{\partial}{\partial \mathbf{J}} \cdot \mathbf{F}, \quad (119) $$

where the two contributions to the “action space flux” $\mathbf{F} = \mathbf{F}_1 + \mathbf{F}_2$ on the right-hand side are

<figure><img src="image_37.png" /><figcaption>FIG. 24. Response of a Maxwellian box of stars to a point mass perturber, adapted from J. Magorrian, Mon. Not. R. Astron. Soc. 507, 4840 (2021). Copyright 2020 Royal Astronomical Society.¹⁰⁵ The point mass (black dot) is driven through the box of size $\approx 0.9L_\odot$ at constant speed, and the wake response is computed according to linear theory. Solid blue contours denote positive overdensity, and red dotted contours denote negative overdensity, while the blue dashed curve corresponds to zero overdensity. In panel (a), the stars only feel the gravitational potential of the perturber but not of each other, while in panel (b), they interact with both the perturber and each other. The same linearly spaced contour levels are used in both panels.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>FIG. 25. Response of a stellar disk to a massive perturber driven on a circular orbit [adapted from D'Onghia et al., Astrophys. J. 766, 34 (2013). Copyright 2013 IOP.¹⁰⁶] The perturber mass is much less than 1% of the mass of the system. Still, it manages to produce a significant response over a large swathe of the disk.</figcaption></figure>

$$F_1(\mathbf{J}, t) = \sum_{\mathbf{n}} \mathrm{in} e^{-\mathrm{in} \cdot \Omega t} \langle \delta f_{\mathbf{n}}(\mathbf{J}, 0) \delta \Phi_{\mathbf{n}}^*(\mathbf{J}, t) \rangle \quad (120)$$

and

$$F_2(\mathbf{J}, t) = -\sum_{\mathbf{n}} \mathbf{n}\mathbf{n} \cdot \frac{\partial f_0}{\partial \mathbf{J}} \int_0^t dt' e^{-\mathrm{in} \cdot \Omega(t-t')} \langle \delta \Phi_{\mathbf{n}}^*(\mathbf{J}, t) \delta \Phi_{\mathbf{n}}(\mathbf{J}, t') \rangle. \quad (121)$$

These equations are not yet written in a form that is useful for calculation, but their physical meaning is already apparent. In order to make this discussion as clear as possible, let us deal with them in reverse order.

* The second flux $F_2$ is a diffusion term. This is obvious since it is proportional to the gradient of the DF $f_0(\mathbf{J})$ in action space. That is to say, this term answers the following question: “I am a star with action $\mathbf{J}$, moving along my orbit within a bath of potential fluctuations $\delta\Phi$. If I think of myself as a massless test particle, then on average, how are those potential fluctuations going to nudge me around?” It is natural that the answer involves the correlation of fluctuations sampled at action space location $\mathbf{J}$. Nothing about the expression (121) would formally change had we insisted that $\delta\Phi$ was externally imposed, but here we are going to use the solution for $\delta\Phi$ derived in Sec. VIC assuming the fluctuations are internally generated and that the system is stable. Thus, we say that the flux $F_2$ describes orbital diffusion due to dressed discreteness noise.

* On the other hand, the first flux $F_1$ provides the answer to the obvious rebuttal: “Wait! I’m not really a test particle, so diffusion cannot be the whole story. I have finite mass, and so I drive potential fluctuations of my own. To what extent is my motion altered by the part of $\delta\Phi$ that I have myself produced?” Physically, this term represents dynamical friction, which Chandrasekhar intuited back in the 1940s. The key idea is that when any massive body moves through a sea of other bodies, it creates a “wake” behind it in response (just as in Figs. 24 and 25). This wake then backreacts, or drags, on the massive body, slowing it down. It is therefore natural that the flux $F_1$ involves the correlation between $\delta\Phi$ at action space location $\mathbf{J}$ and the initial discreteness noise at the same location, $\delta f_n(\mathbf{J}, 0)$. The existence of such a term is also a necessary consequence of energy (and momentum) conservation: since $F_2$ scatters stars to larger energies on

average, in a self-consistent system, there must be another term dragging them back to lower energies.

Although we have taken all stars to be of equal mass for simplicity, one can show that in a multi-mass system, the more massive a body is the larger the wake it induces, and so the stronger the drag force it feels. This leads for instance to massive bodies losing energy and spiraling in toward the center of their host galaxy. Mathematically, for a subsystem of very massive bodies, the frictional term in Eq. (119) is typically much more important than the diffusive term. By contrast, a very light subsystem will contribute only weakly to the total gravitational potential fluctuation, and hence will experience little dynamical friction; for this subsystem, the diffusive term in Eq. (119) can dominate over the frictional one. This is also an astrophysically relevant scenario: for instance, a key problem of galactic evolution is to understand the diffusion of stellar orbits (Fig. 7) due to interactions with transient spiral waves and/or giant molecular clouds.¹⁰⁸ For the present single-mass system, though, $F_1$ and $F_2$ are typically comparable in magnitude.

The next step in deriving a proper quasilinear kinetic equation is to take the solution for the dressed potential from linear theory $\delta\Phi[\delta f(0)]$, Eq. (117), and plug it into the right hand sides of Eq. (120) and (121). One then further assumes that stellar orbits are initially uncorrelated and takes $t \gg t_{dyn}$. The details of this procedure are rather lengthy and do not convey much physical insight, so we refer the reader to, e.g., the works by Chavanis,¹⁰⁹ Fouvry and Bar-Or,¹⁰³ and here just write down the final answer. Equation (119) takes the Fokker-Planck form

$$\frac{\partial f_0}{\partial t} = \frac{\partial}{\partial \mathbf{J}} \cdot \left[ -\mathbf{A}(\mathbf{J})f_0 + \mathbf{D}(\mathbf{J}) \cdot \frac{\partial f_0}{\partial \mathbf{J}} \right], \quad (122)$$

with the friction vector

$$\mathbf{A}(\mathbf{J}) = \pi(2\pi)^d m \sum_{\mathbf{n},\mathbf{n}'} \mathbf{n} \int d\mathbf{J}' \delta(\mathbf{n} \cdot \mathbf{\Omega} - \mathbf{n}' \cdot \mathbf{\Omega}') \\ \times |\psi_{\mathbf{n}\mathbf{n}'}^d(\mathbf{J}, \mathbf{J}', \mathbf{n}', \mathbf{\Omega}')|^2 \mathbf{n}' \cdot \frac{\partial f_0}{\partial \mathbf{J}'} \quad (123)$$

and diffusion tensor
$$
\begin{equation}
\begin{split}
D(\mathbf{J}) = \pi (2\pi)^d m \sum_{\mathbf{n},\mathbf{n}'} \mathbf{n} \otimes \mathbf{n} & \int d\mathbf{J}' \delta(\mathbf{n} \cdot \mathbf{\Omega} - \mathbf{n}' \cdot \mathbf{\Omega}') \\
& \times |\psi_{\mathbf{n}\mathbf{n}'}^d(\mathbf{J}, \mathbf{J}', \mathbf{n}' \cdot \mathbf{\Omega}')|^2 f_0(\mathbf{J}')
\end{split}
\tag{124}
\end{equation}
$$

We can put these expression together to arrive at the *Balescu–Lenard equation*

$$
\begin{equation}
\begin{split}
\frac{\partial f_0(\mathbf{J})}{\partial t} ={}& \pi (2\pi)^d m \frac{\partial}{\partial \mathbf{J}} \cdot \sum_{\mathbf{n},\mathbf{n}'} \mathbf{n} \int d\mathbf{J}' \delta(\mathbf{n} \cdot \mathbf{\Omega} - \mathbf{n}' \cdot \mathbf{\Omega}') \\
& \times |\psi_{\mathbf{n}\mathbf{n}'}^d(\mathbf{J}, \mathbf{J}', \mathbf{n}' \cdot \mathbf{\Omega}')|^2 \left( \mathbf{n} \cdot \frac{\partial}{\partial \mathbf{J}} - \mathbf{n}' \cdot \frac{\partial}{\partial \mathbf{J}'} \right) f_0(\mathbf{J}) f_0(\mathbf{J}').
\end{split}
\tag{125}
\end{equation}
$$

**1. Properties of the Balescu-Lenard equation**

At first glance, the Balescu–Lenard equation might look implausi-
bly complicated, but the physics underlying it is simple. Indeed, with
hindsight, one could argue that this is the simplest equation one could
write down with all the desired properties (particularly that it must
overcome the shortcomings of Chandrasekhar theory that we listed in
Sec. II C 3). For instance:

* The equation is expressed in angle-action variables and therefore naturally accounts for the inhomogeneity of the system and its associated nontrivial orbital structure.

* Related to this, unlike in Chandrasekhar's theory, in which we had to stipulate a maximum impact parameter $b_{max}$ to fix the Coulomb logarithm [e.g., Eq. (17)], the Balescu–Lenard equation suffers from no large scale divergence and no such cutoff is needed. The largest scale of fluctuations in the system is accounted for naturally through the lowest-order basis elements $\Phi^{(p)}$ and wavenumbers $\mathbf{n}$.

* It takes the form of a Fokker-Planck equation in action space, with both friction and diffusion terms present, as it must following the physical arguments below equation (121).

* The right-hand side is an antisymmetric functional of $f_0(\mathbf{J})$ and $f_0(\mathbf{J}')$. This tells us that the system is driven by pairwise interactions. In other words, the Balescu–Lenard equation describes a collection of orbits at action space location $\mathbf{J}$ interacting with orbits at all other locations $\mathbf{J}'$. It is remarkable that in this approximation (truncating at second order in small fluctuations), the collective evolution of an amazingly complicated $N$-star system acts as if it were driven by uncorrelated two-body encounters—though they interact not through the usual bare coupling $\psi$

= $-G/|\mathbf{r} - \mathbf{r}'|$ but through the dressed pairwise interaction potential $\psi^d$. This is a direct manifestation of Rostoker's principle, which serves as the basis for the alternative derivation shown in Sec. VII B.

* The Dirac delta function on the right-hand side encodes the fact that interactions are not only pairwise but also resonant. In other words, the only meaningful interactions occur if some linear combination of frequencies of the star at **J** matches a linear combination of the frequencies of the star at **J'**. This makes sense, since pairs of stars that do not satisfy this condition will undergo an interaction that vanishes on average—each star will feel an equal “pull” and “push” from the other as it traverses many mean field orbits, with a net result of zero.

* The equation conserves mass, energy, and satisfies a H-theorem for Boltzmann entropy. In addition, the right-hand side exactly vanishes for inhomogeneous Boltzmann DFs (i.e., thermal distributions) of the form $f_0(\mathbf{J}) \propto e^{-\beta H_0(\mathbf{J})}$ (although these DFs are not realizable in actual galaxies since they contain infinite mass).

* In the limit where collective amplification is ignorable [M = 0, Eq. (101)], we can replace $\psi^d \to \psi$ and recover the Landau equation which describes the evolution of an N-body system driven by “bare” resonant interactions.¹¹⁰

Finally, in the limit of a homogeneous system, one can use the
angle-action mapping (38) and thereby recover (up to constants) the
classic Balescu-Lenard collision operator from electrostatic plasma
kinetics

$$
\begin{equation}
\begin{split}
\frac{\partial f_0(\mathbf{v})}{\partial t} \propto \frac{\partial}{\partial \mathbf{v}} \cdot \sum_{\mathbf{k}} \mathbf{k} \int d\mathbf{v}' \delta[\mathbf{k} \cdot (\mathbf{v} - \mathbf{v}')] \frac{1}{k^4 |\epsilon_{\mathbf{k}}(\mathbf{k}' \cdot \mathbf{v}')|^2} \\
\times \mathbf{k} \cdot \left( \frac{\partial}{\partial \mathbf{v}} - \frac{\partial}{\partial \mathbf{v}'} \right) f_0(\mathbf{v}) f_0(\mathbf{v}').
\end{split}
\tag{126}
\end{equation}
$$

In Table II, we provide a direct comparison between the various terms
in the homogeneous equation (126) and the corresponding inhomoge-
neous equation (125). In particular, for a homogeneous system, we
have $|\psi_{kk'}^d|^2 \propto \delta_k^k \delta_{kk'}^{-4} |\epsilon_k|^2$, where $\epsilon$ is the dielectric function. This is
why the homogeneous equation (126) involves only a single sum over
wavenumbers $k$, rather than the double sum over $\mathbf{n}, \mathbf{n}'$ that we had in
the inhomogeneous equation (125). The factor $\propto |\epsilon_k(k' \cdot v')|^2$ in Eq.
(126) is the classic dielectric shielding factor. In near-Maxwellian plas-
mas, this factor tends to suppress two-body interactions on large (big-
ger than Debye) scales, but in stellar systems, it can strongly amplify
such interactions—see Sec. VII D for a detailed discussion.

<table><caption>TABLE II. Comparison of various pieces entering the Balescu–Lenard equation for homogeneous [Eq. (126)] and inhomogeneous [Eq. (125)] stellar systems.</caption>
  <thead>
    <tr>
      <th></th>
      <th>Homogeneous</th>
      <th>Inhomogeneous</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Orbital distribution</td>
      <td>f<sub>0</sub>(v, t)</td>
      <td>f<sub>0</sub>(J, t)</td>
    </tr>
    <tr>
      <td>Basis decomposition</td>
      <td>&psi;(r, r') &propto &sum;<sub>k</sub> e<sup>ik&middot;(r&minus;r')</sup>/k<sup>2</sup></td>
      <td>&psi;(r, r') &propto &sum;<sub>p</sub> &Phi;<sup>(p)</sup>(r) &Phi;<sup>(p)*</sup>(r')</td>
    </tr>
    <tr>
      <td>Response matrix</td>
      <td>&propto 1/k<sup>2</sup> &int; dv &nabla; &middot; f<sub>0</sub>/&part;v</td>
      <td>&propto &sum;<sub>n</sub> &int; dJ &nabla; &middot; f<sub>0</sub>/&part;J &Phi;<sup>(p)*</sup>(J) &Phi;<sup>(q)</sup>(J)</td>
    </tr>
    <tr>
      <td>Dressing of perturbations</td>
      <td>&propto 1/&epsilon;<sub>k</sub> (k&prime; &middot; v&prime;)<sup>2</sup></td>
      <td>|&psi;<sup>d</sup><sub>nn'</sub>(J, J', n' &middot; &Omega;')|^<sup>2</sup></td>
    </tr>
    <tr>
      <td>Resonance condition</td>
      <td>k &middot; (v &minus; v') = 0</td>
      <td>n &middot; &Omega; = n' &middot; &Omega;'</td>
    </tr>
  </tbody>
</table>
## 2. An historical aside

The Balescu–Lenard equation was first derived in the context of homogeneous electrostatic plasmas.¹¹¹,¹¹² The main mathematical difference in the stellar dynamical case is that one must formulate the equation in angle-action variables.

The basic idea of formulating a quasilinear diffusion equation in action space seems to have been first arrived at by Kaufman and Nakayama,¹¹³ Kaufman,¹⁰⁰,¹¹⁴ in the context of inhomogeneous plasmas. The idea of doing gravitational kinetic theory in angle-action variables was introduced around the same time by Kalnajs,⁷⁵ who was predominantly interested in the linear properties of disk galaxies and their associated spiral modes (Sec. VI). A proper diffusion equation of the form $\partial_t f_0 = \partial_t \cdot [D \cdot \partial_t f_0]$ was not written down in the stellar-dynamical context until a few years later, by Dekker.¹¹⁵ For a thorough overview of the history and of linear and quasilinear kinetics in angle-action variables (including their plasma precedents), see Chavanis.¹¹⁶

The first derivation of the full Balescu–Lenard kinetic equation for stellar systems was in a somewhat opaque 1987 treatise by Luciani and Pellat.¹¹⁷ The same equation was then arrived at in 2010 by Heyvaerts¹¹⁸ using the BBGKY hierarchy; then in 2012 by Chavanis¹⁰⁹ using the Klimontovich description; again in 2017 by Heyvaerts et al.¹¹⁹ using a Fokker–Planck description; then in 2018 by Fouvry and Bar–Or¹⁰³ using stochastic equations and Novikov’s theorem; and finally in 2021 by Hamilton¹²⁰ using Rostoker’s principle, as we now detail.

## B. Alternative derivation of the Balescu–Lenard equation from Rostoker’s principle

Rostoker’s principle¹⁶,¹²¹,¹²² tells us that we can think of a secularly evolving plasma or stellar system as a collection of N uncorrelated particles undergoing two-body encounters (à la Chandrasekhar), provided we replace the bare Coulombic/Newtonian interaction $\psi(\mathbf{r}, \mathbf{r}')$ of each two-body encounter with its dressed counterpart $\psi^d$. Roughly speaking, this is because particles are weakly coupled, and because the main effect of collective interactions upon the potential fluctuation spectrum is to dress each particle with its collective wake [Eq. (117)]. Here, we use this principle to derive the Balescu–Lenard equation (125) in a short and simple way following Hamilton.¹²⁰ (Note we do not actually need to deal explicitly with fluctuations here, so we will drop the “0” subscript on the ensemble-averaged DF).

First, consider a “test” star with coordinates $(\theta, \mathbf{J})$ and a “field” star with coordinates $(\theta', \mathbf{J}')$. Mathematically, Rostoker’s principle tells us to forget about all the other stars and treat the interaction of these two stars as if they were an isolated system with specific two-body Hamiltonian [units of (velocity)²]

$$h = H_0(\mathbf{J}) + H_0(\mathbf{J}') + m\psi^d(\theta, \mathbf{J}, \theta', \mathbf{J}'), \quad (127)$$

where $H_0(\mathbf{J})$ is the mean field Hamiltonian. Here, $m\psi^d(\theta, \mathbf{J}, \theta', \mathbf{J}')$ is the dressed specific potential energy between a star at phase space location $(\theta, \mathbf{J})$ and a star at $(\theta', \mathbf{J}')$. (It consists of the usual Newtonian attraction plus collective effects; if we ignore these, then $\psi^d \to -G/|\mathbf{r} - \mathbf{r}'|$). Let us expand $\psi^d$ as a Fourier series in the angle variables:

$$\psi^d(\theta, \mathbf{J}, \theta', \mathbf{J}') = \sum_{\mathbf{n},\mathbf{n}'} e^{i(\mathbf{n}\cdot\theta - \mathbf{n}'\cdot\theta')} \psi^d_{\mathbf{n}\mathbf{n}'}(\mathbf{J}, \mathbf{J}', \mathbf{n}'\cdot\mathbf{\Omega}'). \quad (128)$$

Here, $\psi^d_{\mathbf{n}\mathbf{n}'}(\mathbf{J}, \mathbf{J}', \mathbf{n}'\cdot\mathbf{\Omega}')$ is the dressed potential interaction we derived in Sec. VI C, and we have put in the correct frequency dependence $(\omega = \mathbf{n}'\cdot\mathbf{\Omega}')$, since it is physically clear that these are the only frequencies available to the system. (Alternatively, one can cheat: just think of $\psi^d/m$ as some function of $\theta, \theta', \mathbf{J}, \mathbf{J}',$ and then identify the correct functional form at the end by comparison with 125.)

Let us treat the two-body interaction as a perturbation to each star’s motion. Using Hamilton’s equation $d\mathbf{J}/dt = -\partial\mathbf{h}/\partial\theta$, we find that to zeroth order in this perturbation the test and field stars just follow their mean field trajectories $\theta = \theta_0 + \mathbf{\Omega}t$ and $\theta' = \theta'_0 + \mathbf{\Omega}'t$ indefinitely. To first order, similarly to Eq. (60), the result of their interaction is to nudge each other to new actions $\mathbf{J} + \delta\mathbf{J}$ and $\mathbf{J}' + \delta\mathbf{J}'$, respectively, where

$$\delta\mathbf{J}(\theta_0, \mathbf{J}, \theta'_0, \mathbf{J}', \tau) = -m \sum_{\mathbf{n},\mathbf{n}'} \mathrm{i}\mathbf{n} \psi^d_{\mathbf{n}\mathbf{n}'}(\mathbf{J}, \mathbf{J}', \mathbf{n}'\cdot\mathbf{\Omega}') e^{i(\mathbf{n}\cdot\theta_0 - \mathbf{n}'\cdot\theta'_0)} \\ \times \frac{e^{i(\mathbf{n}\cdot\mathbf{\Omega} - \mathbf{n}'\cdot\mathbf{\Omega}')\tau} - 1}{\mathrm{i}(\mathbf{n}\cdot\mathbf{\Omega} - \mathbf{n}'\cdot\mathbf{\Omega}')}. \quad (129)$$

The result for $\delta\mathbf{J}'$ is the same as (129) except we replace the first factor $\mathrm{i}\mathbf{n} \to -\mathrm{i}\mathbf{n}'$. Note that because the two stars nudge each other, this pairwise interaction conserves the total energy of the pair. This is different from considering an interaction between a test star and a background star (as in Sec. II C 1).

Next, following Rostoker, we consider the relaxation of our entire system to consist of nothing more than an uncorrelated set of dressed two-body encounters. Then, it is easy to write down a master equation for the DF $f(\mathbf{J})$. To do so, we account for (a) test stars being nudged out of the state $\mathbf{J}$ and in to some new state $\mathbf{J} + \Delta\mathbf{J}$, as illustrated in Fig. 26(a), and (b) test stars being nudged in to the state $\mathbf{J}$ from $\mathbf{J} - \Delta\mathbf{J}$, as in Fig. 26(b). Processes (a) and (b) are both deterministic and time-reversible, so we must also account for their inverses, which would be represented by the same diagrams but with the arrows pointing in the opposite direction. Let the transition rate density be $w(\Delta\mathbf{J}, \Delta\mathbf{J}'|\mathbf{J}, \mathbf{J}')$. This quantity is defined such that $w(\Delta\mathbf{J}, \Delta\mathbf{J}'|\mathbf{J}, \mathbf{J}') d\Delta\mathbf{J}d\Delta\mathbf{J}'\tau$ is the probability that a given test star with initial action $\mathbf{J}$ is scattered to the volume of phase space $d\Delta\mathbf{J}$ around $\mathbf{J} + \Delta\mathbf{J}$, by a given field star with action $\mathbf{J}'$ that is itself scattered to the volume element $d\Delta\mathbf{J}'$ around $\mathbf{J}' + \Delta\mathbf{J}'$, in a time interval $\tau$ that is much longer than an orbital period

<figure><img src="image_39.png" /><figcaption>FIG. 26. Two possibilities for “scattering” in action coordinates. A test star’s action changes by $\Delta\mathbf{J}$ during an “encounter” with a field star whose action changes from $\mathbf{J}'$ to $\mathbf{J}' + \Delta\mathbf{J}'$. The DF $f(\mathbf{J})$ will be decremented if the test star is kicked out of state $\mathbf{J}$ [process (a)], but will be incremented if it is kicked in to state $\mathbf{J}$ from $\mathbf{J} - \Delta\mathbf{J}$ [process (b)]. To write down a master equation, we must also account for the inverse processes (by reversing the directions of the arrows).</figcaption></figure>