Kinetic theory of stellar systems: A tutorial

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 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 $±1/r2$ 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).^1–3 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.

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 $α$ emission from regions of recent star formation, meaning stars have been born there in the last few tens of Myr (where 1 Myr $≡106$ years), while the bluer light comes from much older stars, up to ages of $∼10$ Gyr (where 1 Gyr $≡109$ years; the Big Bang happened 13.7 Gyr ago, and galaxies formed in the first few Gyr that followed).³ The image is haphazardly 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 pc $∼3×1016$ m. (The distance from the Sun to the next-closest star, Proxima Centauri, is of order $∼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∼1011$ 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 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 $tdyn≳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 $∼10pc$—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∼105$ stars, so the typical distance between stars is $∼(105/pc3)−1/3∼0.02pc$⁠, 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 $r$ in the cluster, the local phase space distribution function $f(r,v)$ depends primarily on the magnitude, rather than the direction, of the velocity $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 $tdyn∼105$ yr. Thus, each star in a globular cluster has undergone $∼105$ orbits in its lifetime.

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

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.

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

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.

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.

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 (⁠ $∼1011$ in the cold disk, $∼105$ 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.

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 $tdyn$⁠. 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 $tdyn$⁠, as we will see (Sec. II C).

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 $tdyn$ as the two galaxies perform their cosmic dance. However, we know from observations and numerical simulations of such galaxies that they ultimately do virialize and that this virialization is nearly complete after only a few $tdyn$ (which in this case still means a few 10⁹ 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 $≫tdyn$⁠. Isolated disk galaxies and globular clusters both fall into this category.

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 $tdyn$⁠. 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.

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, $−|δE|$⁠. By Eq. (11), this will cause the cluster to achieve a higher temperature by $δT=(−|δ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→∞$⁠.

Evaporation is a real and important effect in cluster dynamics; indeed, there are certain star clusters called open clusters with $N∼103$ (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^4,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).¹⁰ However, in general, one has to be very careful about naive application of “standard” statistical mechanical results to gravitational systems.

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, $trelax$⁠.

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 $|Δv|≈|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.¹¹ 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.

Here, we give the simplest version of Chandrasekhar's calculation of the relaxation time,¹² based on his theory of two-body 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, $δ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=−∞$ to $t=∞$⁠. 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 $δv⊥$⁠.

We seem to have a problem here, since taking either $bmin→0$ or $bmax→∞$⁠, or both, causes the total velocity impulse to diverge (albeit in the most undramatic way, i.e., logarithmically). So, what values of $bmin$ and $bmax$ should we take? Usually, one takes $bmax$ to be comparable to the scale of the system $∼R$⁠, since there cannot be any encounters more distant than this (this is different from plasma, in which one can safely take $bmax$ of order the Debye screening length¹³). Moreover, one typically takes $bmin≈b90$ 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 $Λ≈R/b90$ is good enough for obtaining order-of-magnitude estimates, which is all we care about here.

Plugging in characteristic numbers from Table I, we find that for a globular cluster $t2BR∼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 $t2BR∼1010$ 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 $≳$ 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, $tcoll∼(Λp/ log Λp)ωp−1$⁠, where $Λp$ is the plasma parameter (the number of particles per Debye sphere) and $ωp$ is the plasma frequency.¹⁶ Thus, there is a direct analogy between the plasma parameter $Λ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/Λp≪1$ and $1/N≪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).

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

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 $trelax$ and various other evolutionary properties,¹⁴ although the value of the Coulomb logarithm has to be tweaked.^22,23 Moreover, one can use Chandrasekhar's ideas to spin up a full kinetic theory—i.e., a kinetic equation of the form $∂f/∂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 $t2BR≫tdyn$⁠), 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 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 $trelax≈(2000/619)2×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.

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

As we mentioned in Sec. II B, on the timescale $∼tdyn$ 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 $tdyn$⁠. 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.

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: $r(t)=vt+r(0)$ and $v=v(0)=$ 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.^26–28 

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 $Φ(r)$⁠, where $r≡|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, $L=r×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 $(ϕ,R)$ polar coordinate system centered at $r=0$⁠.

In Fig. 9, we show three orbits integrated in the potential (28) in the $(ϕ,R)$ plane, each with the same $Ra=b$ but different $Rp$⁠, as labeled above each panel.

Orbit (a) is nearly radial, with $Rp≪Ra$⁠, whereas orbit (c) is nearly circular, with $Rp≈Ra$⁠, and orbit (b) is somewhere in-between.

The orbit labels $(Rp,Ra)$⁠, 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 $r$⁠. However, it turns out that the phase of the orbits also exhibit rather simple behavior, as we investigate next.

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 $Φ(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, $LZ$⁠).^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 guaranteed³⁴).

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 $Φ(r)$ are of this nature.

The upshot of quasiperiodicity is that stellar orbits can be thought of as a superposition of oscillators with fixed frequencies $Ω$⁠. 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.

A perfectly valid choice for the canonical coordinates $(q,p)$ is the position and velocity $(r,v)$⁠. However, the lack of straight line orbits in stellar systems means that $(r,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(r,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.

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³⁹ and Arnold.⁴⁰ 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(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 $(θ,J)$ corresponding to the corresponding perfectly periodic system and then describe the real (perturbed) system by expressing its Hamiltonian in these coordinates, $H(θ,J)$⁠, deriving the corresponding equation of motion $dθ/dt=∂H/∂J$ and $dJ/dt=−∂H/∂θ$⁠, etc. We will see how this construction works in practice starting in Sec. V B.

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

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.

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 $Ω(J)$⁠. However, a disadvantage of angle–action coordinates is that we can rarely make explicit the mapping $(r,v)↔(θ,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 $Rg$⁠, 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.

For the remainder of this tutorial article, we assume that we can always construct angle–action variables $(θ,J)$⁠, such that the mean field Hamiltonian depends only on $J$ (which is itself a nontrivial technical problem⁴⁴). 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 $v→J$ and $r→θ$ (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., $Ω(J)$ is typically a complicated (nonlinear) function of $J$⁠, whereas in plasma we have $Ω∝J∝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., $δΦ(r,t)=δΦ(θ,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.

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

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.

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

We now fix the pattern speed to be $Ωp=0.154 V0/R0$⁠, so that the corotation radius of circular orbits is $RCR=6.5R0$⁠. We take the spiral strength to be $ε=0.01$⁠, the pitch angle $p=60°$⁠, and the initial phase $ϕp=0$⁠. We initialized the orbits of nine test stars in this potential, putting each star on a circular orbit with speed $V0$ and azimuthal location $ϕ(t=0)=0$⁠. The initial conditions of the stars therefore differed only in their initial orbital radii $R(0)$⁠, which we sampled linearly between $RILR$ and $ROLR$ inclusive. We integrated each orbit for a total time $20Tp$⁠.