but much shorter than the relaxation time. Assuming the system's
equilibrium state is invariant under time reversal, *f* must satisfy the
master equation

$$
\begin{equation}
\begin{split}
\frac{\partial f(\mathbf{J})}{\partial t} = \frac{(2\pi)^d}{m} \int d\mathbf{J}' d\Delta\mathbf{J} d\Delta\mathbf{J}' \frac{1}{2} & \left[ w(\Delta\mathbf{J}, \Delta\mathbf{J}' | \mathbf{J}, \mathbf{J}') [-f(\mathbf{J})f(\mathbf{J}') \right. \\
& \left. + f(\mathbf{J} + \Delta\mathbf{J})f(\mathbf{J}' + \Delta\mathbf{J}') \right] + w(\Delta\mathbf{J}, \Delta\mathbf{J}' | \mathbf{J} - \Delta\mathbf{J}, \mathbf{J}') \\
& \times \left[ f(\mathbf{J} - \Delta\mathbf{J})f(\mathbf{J}') - f(\mathbf{J})f(\mathbf{J}' + \Delta\mathbf{J}') \right].
\end{split}
\tag{130}
\end{equation}
$$

(The factor of 1/2 accounts for the fact that, when we integrate over all
possible **J', ΔJ** and **ΔJ'**, the contribution from each diagram in Fig. 26
accounts for all possible interactions, and both diagrams are equiva-
lent. The "accounting trick" of using both diagrams and then halving
the answer leads very simply to the Balescu–Lenard equation.)

By expanding the integrand in Eq. (130) for weak interactions,
i.e., for Δ**J**, Δ**J'** ≪ **J**, **J'**, up to second order in small quantities, we can
show that

$$
\frac{\partial f(\mathbf{J})}{\partial t} = \frac{\partial}{\partial \mathbf{J}} \cdot \int d\mathbf{J}' \left[ \mathbf{A}(\mathbf{J}, \mathbf{J}') \cdot f(\mathbf{J}') \frac{\partial f}{\partial \mathbf{J}} + \mathbf{B}(\mathbf{J}, \mathbf{J}') \cdot f(\mathbf{J}) \frac{\partial f}{\partial \mathbf{J}'} \right], \quad (131)
$$

where A is the 3 × 3 matrix

$$
\begin{align}
\mathbf{A}(\mathbf{J}, \mathbf{J}') &= \frac{(2\pi)^d}{2m} \int d\Delta\mathbf{J} d\Delta\mathbf{J}' w(\Delta\mathbf{J}, \Delta\mathbf{J}' | \mathbf{J}, \mathbf{J}') \Delta\mathbf{J} \otimes \Delta\mathbf{J} \nonumber \\
&\equiv \frac{(2\pi)^d}{2m} \frac{\langle \Delta\mathbf{J} \otimes \Delta\mathbf{J} \rangle_{\tau}}{\tau}, \tag{132}
\end{align}
$$

where 〈 Δ**J** ⊗ Δ**J** 〉<sub>τ</sub> is the expectation value of Δ**J** ⊗ Δ**J** after a time interval τ for a given test star action **J** and field star action **J'**. The matrix **B** is identical to **A** except we replace 〈 Δ**J** ⊗ Δ**J** 〉<sub>τ</sub> → 〈 Δ**J** ⊗ Δ**J'** 〉<sub>τ</sub>. In a homogeneous system, **J** and **J'** can be taken to be the (suitably normalized) linear momenta of the test and field star; then, momentum conservation implies Δ**J** = −Δ**J'** so that **A** = −**B**. With this, one can easily recover equation (2.36) of the analogous calculation for a homogeneous plasma by Diamond *et al.*<sup>16</sup>

Now we put the above results together. Since by assumption stars
at a given action are randomly distributed in angles, we can calculate
the expectation value 〈 Δ**J** ⊗ Δ**J** 〉<sub>τ</sub> by averaging over the initial phases
θ<sub>0</sub>, θ'<sub>0</sub>. [This is analogous to the integral over impact parameters *b* in
Chandrasekhar's theory, as in Eq. (17).] Thus, we have

$$
\langle \Delta \mathbf{J} \otimes \Delta \mathbf{J} \rangle_{\tau} = \int \frac{d\theta_0}{(2\pi)^d} \frac{d\theta'_0}{(2\pi)^d} \delta(\mathbf{J}(\theta_0, \mathbf{J}, \theta'_0, \mathbf{J}', \tau) \otimes \delta(\mathbf{J}(\theta_0, \mathbf{J}, \theta'_0, \mathbf{J}', \tau), \quad (133)
$$

where δ**J** is given in Eq. (129). Plugging (133) and (129) in to (132)
and taking the limit τ → ∞, and making use of the following two
identities:

$$
\psi_{nn'}^d(\mathbf{J}, \mathbf{J}', \omega_R) = [\psi_{-n, -n'}^d(\mathbf{J}, \mathbf{J}', -\omega_R)]^* \quad \text{for } \omega_R \text{ real,} \quad (134\text{a})
$$

$$
\lim_{\tau \to \infty} \left[ \frac{e^{i\kappa\tau} - 1}{x^2\tau} \right] = 2\pi\delta(x), \quad (134\text{b})
$$

we find the following expression for the matrix A:

$$
\mathbf{A}(\mathbf{J}, \mathbf{J}') = \pi (2\pi)^d m \sum_{\mathbf{n}, \mathbf{n}'} \mathbf{n} \otimes \mathbf{n} \delta(\mathbf{n} \cdot \mathbf{\Omega} - \mathbf{n}' \cdot \mathbf{\Omega}') |\psi_{\mathbf{n}\mathbf{n}'}^d(\mathbf{J}, \mathbf{J}', \mathbf{n}' \cdot \mathbf{\Omega}')|^2. \quad (135)
$$

The result for **B** is identical to Eq. (135) except we replace the factor **n** ⊗ **n** with −**n** ⊗ **n'**. Putting the explicit formulas for **A**(J, J') and **B**(J, J') in to the kinetic equation (131), we recover the Balescu-Lenard equation (125).

Thus, just as in plasma physics,¹⁶ the Balescu–Lenard collision
operator can be interpreted as arising from the uncorrelated pairwise
scattering of *dressed* stars as they traverse their (otherwise unper-
turbed) mean field orbits.

C. Examples

Implementing the Balescu–Lenard equation (125) in practice
is hard, and has only been attempted a handful of times. Here, we
mention some of these implementations, starting with an artificial
one-dimensional gravitational system (Sec. VII C 1), then moving
on to a hot sphere (Sec. VII C 2) and finally a cold disk (Sec.
VII C 3). Other applications include the inhomogeneous
Hamiltonian mean-field model,¹²³ as well as the relaxation of the
eccentricities¹²⁴ and orientations¹²⁵ of stellar orbits around super-
massive black holes.

1. Example: Balescu-Lenard in one-dimensional
"gravity"

Recently, Roule et al.¹²⁶ considered the evolution of a one-
dimensional “gravitational” N-body system. Their system is gravita-
tional in the sense that density and potential are linked through
Poisson’s equation. Yet, because it is one-dimensional the resulting
bare interaction between particles is ψ(x, x') ∝ |x - x'|, and the force
between any two particles is proportional to sgn(x - x'). In this simple
setup, there is a one-to-one mapping between action *J* and energy *E*,
and the authors tend to use *E* to characterize orbits. The initial DF the
authors choose for the particles is closely analogous to the isotropic DF
that generates the three-dimensional Plummer sphere whose orbits we
studied in Sec. III A 1. Crucially for the application of Balescu–Lenard
theory, the DF is linearly stable.

Roule et al.<sup>126</sup> simulated this one-dimensional N-body system directly many times over with random initial conditions, in order to produce an ensemble of systems over which to average. This was to be compared with the prediction from the Balescu–Lenard equation (125). They also performed simulations “without collective effects.” To do this they first integrated the mean field motion of random ensembles of massive but noninteracting particles in their chosen model. They then used the fluctuating gravitational potential generated by these particles as the input to a test particle simulation. Since test particles, by definition, do not produce potential fluctuations of their own, they experience no friction [so F<sub>1</sub> = 0 in Eq. (119)], but they do experience diffusion (F<sub>2</sub> ≠ 0). These test particle simulations were to be compared to the predictions of “Landau” theory, i.e., diffusion along the lines of Balescu–Lenard but without collective amplification (Sec. VII A 1).

Figure 27 shows Roule et al.’s¹²⁶ results. The red crosses in this figure show measurements of the diffusion coefficient in energy of particles over a time long compared to the dynamical time t<sub>dyn</sub> but short compared to the relaxation time t<sub>relax</sub>. The red solid line shows the theoretical “Balescu–Lenard” prediction. It demonstrates an excellent fit, confirming that Eq. (125) is capturing the physics of this system. Meanwhile, the blue crosses and blue solid line shows that the
“Landau” theory is an excellent predictor of the diffusion coefficient in the case without collective effects.

Let us emphasize that the red and blue data in Fig. 27 are plotted on different y axes. In particular, the Landau diffusion coefficient is ≈ 1 in these units, which would correspond to a relaxation time $t_{relax} \sim Nt_{dyn}$, similar to Chandrasekhar's theory [Eq. (19)]. What is very notable is that over a large range of energies $E$, the Balescu–Lenard diffusion coefficient is ten times smaller than the Landau one, corresponding to a relaxation time $t_{relax} \sim 10Nt_{dyn}$. In other words, this is a self-gravitating system (to the extent that this is “gravity” at all) in which collective effects tend to quench, rather than amplify, the rate of evolution.

While this result is unusual, it is not unique. A similar result was obtained by Weinberg,⁴⁵ who studied the dynamical friction drag on a satellite falling into a spherical galaxy (a similar calculation to the one we will present in Sec. VII B 1). Weinberg found that when collective effects were included, the wake that the satellite induced in the galaxy was highly symmetric and that this led to a drag which was weaker than in the bare case without collective effects. The lesson we draw from these examples is that collective effects can have a major impact upon the evolution of inhomogeneous stellar systems, but it may be difficult to guess in advance what that impact will be.

## 2. Example: Balescu–Lenard in a hot sphere

On the surface, Balescu–Lenard theory is much more general than Chandrasekhar's theory of two-body relaxation (Sec. II C 1). One might therefore ask (i) whether it contains Chandrasekhar's theory as a limiting case, and (ii) given that Chandrasekhar's theory works very well in hot spheres (Sec. II C 3), does a direct application of the results of Balescu–Lenard theory to those systems give an answer that differs from Chandrasekhar's theory?

The answer to question (i) is yes: one can recover Chandrasekhar's theory from the Balescu–Lenard formalism by ignoring collective effects and confining oneself to interactions on small scales over which the curvature of mean field orbits can be neglected (replacing them with straight lines¹¹⁶). To answer question (ii) requires a lot more work. Fouvry et al.²⁴ investigated it by evaluating the right

hand side of Eq. (125) for an isochrone sphere, which is a toy model of a spherical star cluster. Using a known (stable) DF, they calculated the Balescu–Lenard flux F as a sum of contributions $F^\ell$ from spherical harmonic quantum numbers $\ell = 0, 1, 2, \dots$. As illustrated in Fig. 28, they found that on small scales $\ell \gg 1$, the modulus $|F_\ell|$ scales as $\propto \ell^{-1}$. Hence, a sum over $\ell$ recovers the small-scale logarithmic divergence familiar from Chandrasekhar's theory, a result already obtained by Weinberg¹²⁷ when considering dynamical friction (see Fig. 5 therein). Such a divergence is always present in a kinetic theory based on weak $\pm 1/r^2$ interactions. This strongly suggests that Chandrasekhar's theory—which is much simpler to implement than the Balescu–Lenard theory—is being recovered by Eq. (125) on small scales. This motivated Fouvry et al.²⁴ to calculate the Balescu–Lenard flux as accurately as possible on large scales ($\ell \le \ell_{crit}$, with $\ell_{crit} \approx 12$, say), and to calculate the Chandrasekhar flux on small scales ($\ell > \ell_{crit}$), and then add the two contributions together. The advantage of this over just using Chandrasekhar's theory is that one is not plagued by any large-scale divergence (which is normally cured by imposing a maximum impact parameter $b_{max}$), and in principle, one accounts for any resonant and/or collective phenomena that are ignored by Chandrasekhar. The result of this exercise, which is shown in Fig. 29, is a slightly improved agreement between the contour map of $\partial f_0 / \partial t$ in action space computed from theory as compared to N-body simulations. As we discussed in Sec. II C 3, the agreement based on Chandrasekhar's theory alone is already good, and as such the large-scale Balescu–Lenard “correction” is small.

In summary, perhaps the most intellectually satisfying approach to spherical cluster kinetics is to couple a large-scale Balescu–Lenard flux with a small-scale Chandrasekhar flux. However, in practice, Chandrasekhar's theory applied to all scales (with an additional cutoff at large scales $b_{max}$) is almost as good a description and is both much simpler to understand and easier to implement.

## 3. Example: Balescu–Lenard in a cold disk

As a final example, let us discuss the dynamics of an initially stable, razor thin stellar disk, following the N-body simulations of

<figure><img src="image_40.png" /><figcaption>FIG. 27. Diffusion coefficients of particles in a one-dimensional self-gravitating system (Sec. VII C 1). Crosses show the results of N-body simulations, while solid curves show the theoretical results based on the Balescu–Lenard equation (125) (when collective interactions are switched on) or the corresponding Landau equation (when they are switched off). Note that the vertical axis scale is different for red and blue points. Figure reproduced from Roule et al., Phys. Rev. E 106, 044118 (2022). Copyright 2022 American Physical Society.¹²⁶</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>FIG. 28. Balescu–Lenard flux, $F_\ell$ [see Eqs. (119)–(121)], for fluctuations with spherical harmonic quantum number $\ell$ (an angular scale $\sim \ell^{-1}$). The magnitude of this flux scales as $\ell^{-1}$ for large $\ell \gg 1$. This confirms that the Balescu–Lenard equation recovers the Coulomb logarithm behavior from Chandrasekhar's theory in the limit of local encounters (Sec. II C 1). Figure reproduced from Fouvry et al., Mon. Not. R. Astron. Soc. 508, 2210 (2021). Copyright 2021 Royal Astronomical Society.²⁴</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>FIG. 29. Contour plot of the evolution rate $\partial f_0 / \partial t$ in action space, for a simple model of a spherical cluster as measured from (a) direct N-body simulations and (b) a theory which combines the Balescu–Lenard equation (125) on large scales with the Chandrasekhar theory on small scales. Figure reproduced from Fouvry et al., Mon. Not. R. Astron. Soc. 508, 2210 (2021). Copyright 2021 Royal Astronomical Society.²⁴</figcaption></figure>

Sellwood¹²⁸ and the kinetic predictions from Fouvry et al.¹²⁹ The model employed by these authors is the stable Mestel disk, with a DF $f_0$ whose $J_R = 0$ slice we plotted with a dotted black curve back in Fig. 22(a). Notably, this DF is very smooth in action space. However, Sellwood’s¹²⁸ simulations of this disk shows that after many dynamical times $t_{dyn}$, a grooved/ridged feature always develops in the action space DF [see Fig. 30(a)]. The time taken for this feature to develop increases linearly with the number of particles $N$ employed, suggesting it is driven by the finite-$N$ noise. Remarkably, this newly modified DF is found in the simulations to be unstable to the formation of exponentially growing spiral waves.

In fact, this modified DF is precisely the one studied by De Rijcke et al.⁹⁴ in their linear stability analysis (Sec. VI B 2), the $J_R = 0$ slice through which is shown with a solid black curve in Fig. 22(a). The corresponding spiral Landau mode they calculated from linear theory [which is shown in Fig. 22(b)] matches that extracted by Sellwood¹²⁸ from his simulations, confirming that his spiral waves are true linear instabilities of the modified DF. Thus, our stable, isolated stellar disk has driven itself unstable!

We can usefully contrast this behavior with that of a plasma. If we initiated, e.g., a finite-$N$ electrostatic plasma with a Maxwellian DF, then the plasma would exhibit thermal fluctuations $\delta\Phi$, $\delta f \propto N^{-1/2}$,

but its mean-field Maxwellian $f_0$ would simply never evolve. That is because the Maxwellian is the ultimate maximum entropy state in a homogeneous plasma, a perfect thermodynamic equilibrium. Yet, for self-gravitating systems, there is no such thermodynamic equilibrium (Sec. II B), and so even a linearly a stable stellar disk is never at rest: instead, it gradually carves a groove in itself and drives itself unstable, even without any external influences. It is as if a stable plasma were to spontaneously grow a bump on its own tail!

Before the system becomes unstable, its evolution can be predicted using the Balescu–Lenard equation (125), which after all is supposed to describe the long-term evolution of stable systems driven by amplified finite-$N$ noise. Fouvry et al.¹²⁹ evaluated the right hand side of equation (125) for the same Mestel disk DF used by Sellwood.¹²⁸ The resulting prediction is shown in Fig. 30(b). While there are some important differences between this panel and the result in panel (a), the key physics is quite clear: the grooved/ridged DF is formed because of dressed, resonant coupling between finite-$N$ fluctuations in the disk, which happen to be strongly enhanced along a particular direction in action space (in this case $\mathbf{n} = (2, -1)$, shown with light blue lines in the figure). Unfortunately, the Balescu–Lenard equation is so cumbersome to implement in practice that nobody has managed to push the theory any further than this.

<figure><img src="image_43.png" /><figcaption>FIG. 30. Change in the DF $f_0(L, J_R)$ for an initially smooth, razor thin stellar disk, (a) as measured in an N-body simulation, (b) as predicted by the Balescu–Lenard equation (125) [Figure adapted from Fouvry et al., Astron. Astrophys. 584, A129 (2015). Copyright 2015 European Southern Observatory.¹²⁵] The DF develops a ridge/groove feature which renders it unstable to spiral Landau modes—see Fig. 22.</figcaption></figure>
## D. Breakdown of Balescu-Lenard theory: Fundamental differences between galaxies and plasmas

As discussed at the beginning of this section and in Sec. V E, any linear or quasilinear kinetic theory will break down if strongly nonlinear effects such as orbit trapping are important to the evolution. We will develop a method for dealing with such systems in Sec. VIII. However, there is another way in which the Balescu-Lenard theory can break down even in the absence of strong nonlinearities. Namely, the Balescu-Lenard equation (125) is derived under a key assumption which is that we may ignore the direct influence of the Landau modes. This assumption is often a poor one when applied to real galaxies. This important point is often neglected in the literature, and a discussion of it allows us to highlight some fundamental differences between stellar systems and plasmas.

First, the Balescu-Lenard theory is designed to describe the evolution of an N-body system in which dressed star-star interactions alone drive the evolution. Mathematically, this assumption enters the Balescu-Lenard theory when one uses Eq. (117) for the potential fluctuations. This equation is valid on timescales long compared to the damping time of the system's most weakly damped Landau mode $\omega_m$, which requires

$$t_{\text{relax}} \gg 1/|\text{Im } \omega_m|. \quad (136)$$

This requirement should already be making us worried. Recall that for real galaxies $t_{\text{relax}}$ may be of order only a few Gyr, and many stellar systems are capable of harboring weakly damped Landau modes that live for many $t_{\text{dyn}}$, which is $\approx 250$ Myr in the Milky Way. Moreover, the right-hand side of (136) by definition diverges during the transition from a stable to unstable DF described in Sec. VII C 3, so the Balescu-Lenard theory is clearly invalid in that regime. To make matters worse, in this regime the Balescu-Lenard equation actually gives nonsensical answers. To see this, note that a Landau mode with $\text{Im } \omega_m \to 0$ corresponds to a divergence in $|\psi|^d$ [Eq. (115)] at some real frequency, and the Balescu-Lenard equation (125) involves an integral of $|\psi|^d|^2$ over all real frequencies. Thus, a naive application of the Balescu-Lenard theory to this marginally stable regime will necessarily return an unphysical, divergent result.

How might we do better? Note that in plasma kinetics, the classic quasilinear (CQL) theory of wave-particle interactions⁹⁹ does just the opposite of the Balescu-Lenard theory: it includes the Landau modes while ignoring the dressed discreteness noise. (Throughout this discussion, we use the abbreviation “CQL” to distinguish the classic theory of particles interacting with Landau modes from the generic quasilinear scheme of constructing kinetic equations by multiplying two linear quantities together, see Sec. V D). This assumption is justified in certain plasma regimes. In fact, in near-Maxwellian electrostatic plasmas, the Balescu-Lenard and CQL regimes are physically separated by the Debye length $\lambda_D$.

Let us be more precise. In an electron plasma, the dressing of finite-N noise amounts to Debye shielding (Fig. 23). This means that on scales larger than $\lambda_D$ the amplitude of “dressed noise” drops to zero and one can ignore particle-particle interactions altogether, i.e., the contribution to the Balescu-Lenard equation from these scales is negligible. Meanwhile, there exist Landau modes on all scales: these are Langmuir waves, which involve longitudinal oscillations of electrons as they try to maintain the plasma’s quasi-neutrality. Langmuir waves oscillate very rapidly—typically their associated timescale, the inverse

plasma frequency $\omega_p^{-1}$, is the shortest physical timescale in the system. Now, we know from Landau's theory of wave-particle damping that these modes will damp only if they can resonate with electrons with velocities $v \approx \omega_p/k$, where $k$ is the wavenumber of the Langmuir oscillation in question. Thus, for small $k$ (large scales), the only electrons capable of resonating with the wave are those moving extremely fast, and since these are few in number, small $k$ Langmuir waves tend to be only very weakly damped, and so wave-particle interactions must be taken into account in the kinetic theory via a CQL collision operator. Conversely, for very large $k$, the waves are able to resonate with relatively slowly moving particles (the “thermal bulk”), of which there are many—hence, large-$k$ Langmuir oscillations are typically very strongly damped. Again, it turns out that the separation between these two regimes occurs at the Debye scale $k_D = 2\pi/\lambda_D$. In short, in an electron plasma, the electric potential fluctuations on scales $k \gg k_D$ are dominated by collectively dressed noise and therefore well-described by Balescu-Lenard theory (a “noise-only” regime), whereas on scales $k \ll k_D$ Landau modes dominate and one may use CQL theory (“Landau-modes-only”).

The reader has already anticipated where this argument fails when applied to self-gravitating systems. In this case, there is no Debye shielding, so there is nothing to suppress discreteness noise on any scale. On small scales (say, the interparticle distance) collective effects are very weak, so the power spectrum just reflects bare Poisson noise. This is why Chandrasekhar's theory works so well there (Sec. VII C 2). On large scales (comparable to the size of the system), there typically exists a combination of weakly damped Landau modes and collectively amplified finite-N noise, and crucially, it is near the frequencies of weakly damped Landau modes that the collective amplification of the finite-N noise tends to be most pronounced. In other words, the power spectrum of the continuum resulting from dressed noise ($\propto |\psi_{nn'}^d|^2$) is sharply peaked at frequencies $n' \cdot \Omega'$, which are very close to the (real part of the) frequency $\omega_m$ of weakly damped Landau modes. This argument shows the logical gap that can arise when naively applying the Balescu-Lenard equation to stellar systems. For in many cases, it is only worth putting in the effort to implement the Balescu-Lenard equation if its predictions are significantly different from those of the much simpler Chandrasekhar two-body relaxation theory, i.e., if there is significant collective dressing of finite-N noise on large scales. However, this dressing in turn is typically only significant if the system harbors weakly damped Landau modes. It is therefore somewhat self-contradictory to derive Balescu-Lenard theory by arguing these modes fade away rapidly—if they really did so, there might be no interesting collective behavior to speak of, and therefore no need for the Balescu-Lenard equation! (The 1D system discussed in Sec. VII C 1 is an exception to this rule, since in that case the collective effects tended to suppress interactions rather than amplify them, and there are no weakly damped Landau modes.)

It follows that typically in stellar systems, amplified noise and Landau modes cannot naively be decoupled. A comprehensive kinetic theory must therefore account for both types of potential fluctuation self-consistently. This is also a challenge that arises in plasma kinetics since neither the Balescu-Lenard or CQL theory is strictly valid on scales $k \sim k_D$, and a general theory must incorporate both effects. Several plasma theorists attempted to formulate such a general theory in the 1960s, notably Rogister and Oberman,¹³⁰,¹³¹ Oberman¹³² and Hatori.¹³³,¹³⁴ These authors considered a homogeneous electrostatic
plasma and derived a set of self-consistent coupled equations for the particle distribution function and the energy in the electric field fluctuations, accounting for both Balescu–Lenard and CQL-type effects. It turns out that the divergence of the Balescu–Lenard part of the theory when Im $\omega_m \to \pm 0$ is accompanied by a corresponding divergence in the opposite sense in the CQL contribution. This renders the generalized theory divergence-free for any Im $\omega_m$. (In fact, this cancellation of divergences arises already in linear theory—see Hamilton and Heinemann.¹³⁵) A complete, unified kinetic theory of stellar systems would account for both types of potential fluctuations (dressed particle noise and Landau modes), while also taking account of the inhomogeneity of the system using angle-action variables. This is a topic of ongoing research.¹³⁵,¹³⁶

## VIII. NONLINEAR KINETIC THEORY

The evolution of a system's mean DF $f_0$ is always “nonlinear” in the sense that it depends on the product of two fluctuations $\delta\Phi$ and $\delta f$ [Eq. (118)]. In Sec. VII, we discussed a quasilinear approximation to this evolution in which those fluctuations were each calculated using linear theory. Here, we discuss nonlinearities that cannot be captured with a quasilinear approach. We begin in Sec. VIII A by discussing some of the ways in which these nonlinearities develop, and we identify one that is typically most important for stellar systems, namely, when orbits are strongly distorted and/or trapped by a single coherent perturbation. Then, in Sec. VIII B, we write down and solve a kinetic equation which describes this regime, using the pendulum machinery that we developed for individual orbits in Sec. IV F. Finally, we give two concrete examples of this nonlinear kinetic theory in practice: first in the context of galactic bar–dark matter halo interactions and dynamical friction (Sec. VIII B 1), and then regarding the saturation of spiral instabilities (Sec. VIII B 2).

### A. Which nonlinear regimes are relevant for stellar systems?

Consider a stellar system that is initially completely quiet ($\delta f = \delta\Phi = 0$), and subject it to some perturbation. Then, we know that at least for a short time, the system will be well described by linear theory (Sec. VI), and so, e.g., its potential/density response will consist of a combination of amplified wakes and Landau modes. If the system is unstable (and in some cases even if it is weakly stable), the response will grow large enough that the system will eventually cease to be well described by linear theory.

Why, physically, does linear theory break down? Let us list a few of the main possibilities.

* One reason can be that the DF $f_0$ relaxes quasilinearly, in the manner already discussed in Sec. VII. Since linear theory is based on a constant $f_0$ [e.g., Eq. (94)], such an evolution will eventually break the initial linear prediction (effectively “resetting the initial value problem” with a new DF and fluctuation spectrum). This will be the dominant source of nonlinearity if the evolution time of $f_0$ according to the relevant quasilinear theory—which is naively on the order of $t_{relax} \sim \epsilon^{-2}t_{dyn}$, where $\epsilon \sim |\delta\Phi/H_0|$ is the typical fluctuation strength, see Eq. (118)—is shorter than the timescale over which any alternative nonlinearities can develop.

* One such alternative, which is especially important if the system is unstable, is that nonlinear evolution can arise in the form of direct coupling between Landau modes. That is, the (initially independent) Landau mode solutions attain such large amplitudes that they no longer respond solely to the mean Hamiltonian $H_0$ and the action space gradient $\partial f_0 / \partial J$ (as would be dictated by linear theory), but instead start to notice each other's gravitational fields. Landau mode coupling can lead to modifications both of the potential fluctuation spectrum and the mean stellar distribution function $f_0$. For this effect to be significant, the system usually has to harbor at least two Landau modes that reach amplitudes $\epsilon \sim 1$ on a short enough timescale that they can interact before being destroyed by some other process (such as the quasilinear evolution referred to above). The specific nonlinear physics involved here depends on whether (i) only a handful of Landau modes are excited to large amplitudes, or (ii) a broad spectrum of many such modes is excited.

Case (i), in which only a handful of modes are important, has been studied to some extent in galactic dynamics, mostly in the context of global spiral–spiral interactions¹³⁷,¹³⁸ (in what plasma physicists would refer to as “three wave coupling”). However, the timescale required for even a rather vigorously growing global spiral instability to reach a nonlinear amplitude is usually on the order of several Gyr (see Sec. VIII B 2 for a quantitative estimate). This is often comparable to the timescale over which the bulk of the galaxy tends to evolve (Sec. II C 3). This leaves little time for such modes to interact nonlinearly (recall from Sec. I that a typical disk galaxy is only ≈ 50 orbits old!) As a result, spiral–spiral coupling and similar effects seem to be present but not dominant in galactic disk evolution.

Case (ii) on the other hand, in which many modes interact, is usually referred to by plasma physicists as the regime of “wave turbulence.”⁶⁶ That is, plasma physicists tend to think of such a system as consisting of “particles” (described by $f_0$) and “waves” (Landau modes, which dominate the fluctuation spectrum). Then, they construct kinetic theories in which wave–wave coupling plays a dominant role, and in which the evolution of $f_0$ is deemed secondary or just ignored altogether. Some studies of this flavor have been conducted in galactic dynamics, mainly with reference to tightly wound spiral waves.¹³⁹⁻¹⁴¹ Some authors in the 1970s¹⁴²,¹⁴³ took the plasma analogy so far as to construct versions of spiral structure theory governed by a nonlinear Schrödinger equation, and attempted to describe spirals as solitons/modulational instabilities in the way that Zakharov¹⁴⁴ described Langmuir collapse. However, no serious numerical evidence was ever offered to bolster this formalism. Moreover, as in case (i), it is only approximately correct to freeze the stellar distribution function $f_0$ in order to isolate the evolution of nonlinearly interacting modes. Even if it were possible, the dynamical youth of galaxies means that the number of wave cycles over which such modes can possibly have interacted is much less than is normally assumed in plasma theory.

* A third possibility is for the system to be “dynamically degenerate.” For example, with $\partial f_0 / \partial J = 0$ and $\partial H_0 / \partial J = 0$, Eq. (93) exactly vanishes and the linearized dynamics of fluctuations becomes trivial. The system’s evolution is then intrinsically nonlinear, and one must resort to statistical closure theory to make
progress.⁶⁸ Astrophysically such regimes occur when describing the diffusion of stellar orbital orientations around supermassive black holes¹⁴⁵,¹⁴⁶ or dynamical friction in harmonic cores in the center of galaxies.⁴⁹,¹⁴⁷ This regime bears deep connexions with (anomalous) two-dimensional diffusion in magnetized plasmas.¹⁴⁸

* A fourth possibility is that stellar orbits can be trapped, or otherwise strongly distorted, by a single coherent perturbation. For instance, stars can be trapped at the corotation resonance of a rotating bar or spiral [Fig. 14(e)]. This effect can be important if the libration time $t_{\text{lib}}$ of trapped orbits is shorter than or comparable to the evolution time of the distribution function $f_0$ in the relevant region of phase space, and shorter than the timescale over which some other effect (like incoherent collisions or the sweeping of resonances through phase space) is able to interrupt the libration process, as discussed in Sec. IV F 2.

The last of these nonlinear effects, orbit trapping by a single perturbation, seems to be important in sculpting the evolution of both real and simulated galaxies,⁵¹ and there is good evidence that it has occurred in our own Galaxy.¹⁴⁹,¹⁵⁰ Therefore, for the rest of this section, we will focus only on this type of nonlinearity. We already showed in Sec. IV F how one can describe individual orbits in this regime using the pendulum approximation. Below, we show how a collection of pendulum orbits can be described in a kinetic formalism, although we will ignore self-consistency (i.e., Poisson's equation). We will, however, allow for the possibility that the nonlinearity is interrupted due to stochastic collisions, by appending an *ad hoc* diffusion operator to our (otherwise collisionless) kinetic equation.

We emphasize, though, that in certain contexts more than one of the above nonlinear effects can be important in different parts of the system and at different times. For instance, in Sellwood and Carlberg's¹⁵¹ simulations of grooved unstable stellar disks a spiral Landau mode initially grows exponentially in accordance with linear theory. Once it reaches a large enough amplitude, it begins to trap stars nonlinearly at its corotation resonance, which causes its amplitude to saturate (Sec. VIII 2 2). At the same time, the spiral interacts quasilinearly with stars at the Lindblad resonances, draining its angular momentum and carving a new groove in the DF at the location of the ILR. This newly grooved disk is unstable to a new set of linear spiral eigenmodes, and so on in a recurrent cycle.¹⁵²,¹⁵³

## B. Kinetic theory in the pendulum approximation

In Sec. IV F, we described the nonlinear dynamics of individual stars orbiting near resonances with a long-lived potential perturbation rotating rigidly with some fixed pattern speed. We did this by isolating a particular resonance, defining an associated set of slow and fast angle-action variables, and averaging out the fast angles. We were left with dynamics in the slow angle-action space $(\phi, I)$ governed by the pendulum Hamiltonian $h$ [Eq. (65)], all at a given (fixed) fast action $J_f$. We now study how an ensemble of stars behaves near such a resonance, following Hamilton *et al*.¹⁵⁴ Let us therefore consider a subset of stars that all share the same fast action $J_f$. We will describe the density of these stars in the $(\phi, I)$ plane using a smooth DF which we call $f(\phi, I, t)$, defined such that the number of stars in the phase space area element $d\phi dI$ surrounding $(\phi, I)$ at time $t$ is proportional to

$d\phi dI f(\phi, I, t)$. Since $(\phi, I)$ are canonical variables, the kinetic equation governing $f$ is

$$ \frac{\partial f}{\partial t} + [f, h] = C[f], \qquad (137) $$

where $[\cdot, \cdot]$ is the Poisson bracket encoding the smooth advection in the $(\phi, I)$ plane

$$ [f, h] = \frac{\partial f}{\partial \phi} \frac{\partial h}{\partial I} - \frac{\partial f}{\partial I} \frac{\partial h}{\partial \phi}, \qquad (138) $$

and we have included a possible “collision operator” $C[f]$, which encodes the effect of any stochastic fluctuations in the potential.

For concreteness, let us consider the interaction of a collection of stars (described by $f$) interacting with a rigidly rotating stellar “bar” (a collection of millions of stellar orbits at the center of a galaxy that rotates like a dumbbell). Then, for $h$, we can use the pendulum form (65) with $k = 1$. Meanwhile, for $C[f]$, we choose a very simple diffusive form with constant diffusion coefficient $D$,

$$ C[f] = D \frac{\partial^2 f}{\partial I^2}. \qquad (139) $$

This corresponds to every star in our ensemble being kicked with a stochastic white noise forcing in the slow action, $\dot{I} = \eta(t)$, with statistics $\langle\eta(t)\rangle = 0$ and $\langle\eta(t)\eta(t')\rangle = 2D\delta(t-t')$. (For simplicity, we will ignore any such forcing in $\phi$ and $J_f$.) The diffusion coefficient $D$ can be calculated from a given theoretical model of, e.g., scattering by passing stars and molecular clouds,¹⁰⁸,¹⁵⁵,¹⁵⁶ dark matter substructure,¹⁵⁷⁻¹⁵⁹ spurious numerical heating in simulations,¹⁶⁰,¹⁶¹ or can potentially be calibrated from data²¹,¹⁶² (see below). We emphasize that we make no attempt at self-consistency here. We simply impose a diffusion coefficient by hand, and we do not account at any stage for the self-gravity of the perturbed DF (the stars interact only with the bar, not with each other¹⁶³,¹⁶⁴).

Before attempting to solve the kinetic equation (137), one more simplification is in order. Naively, Eq. (137) seems to depend on three parameters: $G, F$ [Eq. (65)], and $D$ [Eq. (139)]. However, we can reduce this to one effective parameter by introducing certain dimensionless variables. First, we note that the typical timescale for a star to diffuse all the way across the resonance (a distance $\approx 2I_h$ in action space) is the diffusion time

$$ t_{\text{diff}} \equiv \frac{(2I_h)^2}{2D} = \frac{8F}{GD}. \qquad (140) $$

Now let us define the dimensionless variables

$$ \tau \equiv \sqrt{GFt} = \frac{2\pi t}{t_{\text{lib}}}, \qquad (141a) $$

$$ j \equiv \sqrt{\frac{G}{F}} I = \frac{2I}{I_h}, \qquad (141b) $$

$$ \Delta \equiv \sqrt{\frac{G}{F^3}} D = \frac{4 t_{\text{lib}}}{\pi t_{\text{diff}}}, \qquad (141c) $$

where the libration time $t_{\text{lib}}$ is given in Eq. (67) and the island half-width $I_h$ is defined in Eq. (68). Clearly, $\tau$ is a dimensionless measure of time normalized by the libration time, $j$ is a dimensionless measure of the slow action variable relative to the resonant island width, and $\Delta$ is
the diffusion strength, i.e., the ratio of the libration timescale to the diffusion timescale. Treating f as a function of these variables, i.e., writing¹⁶⁵ $f(\varphi, j, \tau)$, Eq. (131) becomes

$$ \frac{\partial f}{\partial \tau} + j \frac{\partial f}{\partial \varphi} - \sin \varphi \frac{\partial f}{\partial j} = \Delta \frac{\partial^2 f}{\partial j^2}. \quad (142) $$

We see that in these variables, the kinetic equation depends on a single parameter, the dimensionless diffusion strength $\Delta \approx t_{\text{lib}}/t_{\text{diff}}$. The regime $\Delta \gg 1$ corresponds to very strong diffusion, whereas $\Delta \ll 1$ corresponds to very weak diffusion; the “collisionless limit” in which stars only interact with the bar and nothing else is $\Delta = 0$. We can estimate $\Delta$ heuristically by noting that $t_{\text{diff}}$ is related to the relaxation time $t_{\text{relax}}$, which is the time required for any stochastic process to change a star’s action by order of itself. Then, from Eq. (140), we have

$$ \frac{t_{\text{diff}}}{t_{\text{relax}}} \sim \left( \frac{I_h}{J_{\text{s, res}}} \right)^2 \sim \left| \frac{\Psi_1}{H_0} \right| \sim \varepsilon \ll 1. \quad (143) $$

Here, $\varepsilon \sim |\Psi_1/H_0|$ is the dimensionless strength of the bar; for the Milky Way at corotation, $\varepsilon \approx 0.02$.¹⁴⁹ Putting Eqs. (141c) and (143) together, we estimate

$$ \Delta \approx \frac{t_{\text{lib}}}{\varepsilon t_{\text{relax}}} \approx 1 \times \left( \frac{\varepsilon}{0.02} \right)^{-1} \left( \frac{t_{\text{lib}}}{2 \text{Gyr}} \right) \left( \frac{t_{\text{relax}}}{100 \text{Gyr}} \right)^{-1}. \quad (144) $$

Typically, important Galactic bar resonances like the corotation resonance have libration times of $t_{\text{lib}} \approx 2$ Gyr.¹⁵⁰ The relaxation time $t_{\text{relax}}$ depends on what we consider to be driving the diffusion, i.e., whether it is two-body relaxation (Sec. II C 1), transient spiral arms (Sec. IV F 2), finite-N granularity noise (Sec. VII), or whatever. In Sec. II C 3, we estimated that the relaxation time of angular momenta in the Milky Way disk is $t_{\text{relax}} \approx 60$ Gyr, which would give $\Delta \approx 2$. More generally, the estimate (144) tells us that $\Delta$ can be of order unity even for systems with relaxation times much longer than the age of the Universe. The important point is that it is much easier to diffuse across the width of the resonance than it is to diffuse across the whole of action space, and it is the diffusion across resonances that governs the dynamics.⁵⁹,⁶⁰,⁶²

We can solve equation (142) analytically in the limits $\Delta = 0$ and $\Delta \gg 1$, but it is easier to just solve it numerically for a given $\Delta$. To do this, let us take as our initial condition a simple linear DF

$$ f(\varphi, j, \tau = 0) = f_0(j = 0) + \alpha j \equiv f_{\text{init}}(j), \quad (145) $$

where $\alpha \equiv \partial f_0 / \partial j|_{j=0}$ is the gradient of the unperturbed DF at resonance. In Fig. 31, we show with colored contours the resulting evolution of the dimensionless auxiliary DF

$$ g(\varphi, j, \tau) \equiv \frac{f(\varphi, j, \tau) - f_0(j = 0)}{\alpha}, \quad (146) $$

whose initial value is $g_{\text{init}}(j) \equiv g(\varphi, j, \tau = 0) = j$. Each row corresponds to a different value of $\Delta$, and within each row, from left to right we plot the solution at different “times” $\Delta^{1/3} \tau$. First consider row (a), which is in the limit of very weak diffusion ($\Delta = 0.001$). Physically, we expect that in this limit the DF will phase mix (c.f. Fig. 19) around and within the resonant island. The reason for phase mixing is that in the

absence of diffusion ($\Delta = 0$), the trajectories of individual stars trace contours of constant $h$ in the $(\varphi, j)$ plane. And, since adjacent contours correspond to slightly different libration/circulation periods, the initial DF gets sheared out along these contours, producing ever finer-scale structures in the phase space. If $\Delta$ were exactly zero, and in the absence of any coarse-graining, this process would continue indefinitely. However, for our very small but finite value of $\Delta = 0.001$, diffusion is able to wipe out the very smallest-scale features, without needing to invoke any coarse-graining. We see that by the third column, $g$ is approximately zero everywhere inside the separatrix, and is smeared almost evenly on contours of constant $h$ outside of the separatrix. The DF has reached a steady state; crucially, in this state the DF is symmetric around $\varphi = 0$ at any fixed $j$.

Next consider row (c) in Fig. 31, which corresponds to rather strong [but not unrealistic!—see Eq. (144)] diffusion, $\Delta = 1$. In this case, the resonance has a much less dramatic effect on $g$ at any time. This again is as expected since the initial linear DF (145) is annihilated by the collision operator $C[f]$ [Eq. (139)]. In other words, wherever the bar perturbation induces some curvature in the DF, strong diffusion immediately tries to remove it and hence to restore the linear initial condition. Another way to say this is that the effect of diffusion is to counteract the phase mixing of orbits, in particular within the resonant island. Importantly, this means that in steady state the DF is strongly asymmetric in $\varphi$ at a given $j$, in contrast to row (a). Finally, the case $\Delta = 0.1$ in row (b) is intermediate between these two regimes.

It is also instructive to average the solution $f(\varphi, j, \tau)$ over the slow angle $\varphi$ and investigate the resulting DF

$$ \langle f \rangle_{\varphi} \equiv \frac{1}{2\pi} \int_{-\pi}^{\pi} d\varphi f(\varphi, j, \tau). \quad (147) $$

In Fig. 32, we plot the corresponding auxiliary DF [see Eq. (146)]

$$ \langle g \rangle_{\varphi} = \frac{\langle f \rangle_{\varphi} - f_0(j = 0)}{\alpha}, \quad (148) $$

as a function of $j$ for various $\Delta$ values and at different times (measured in units of the libration time $t_{\text{lib}} = 2\pi/\omega_{\text{lib}}$).

The key point is that for $\Delta \to 0$, the phase mixing of orbits in and around the resonant island causes the angle-averaged DF to be flattened there (see the blue line in panel c). Diffusion ($\Delta > 0$) suppresses this flattening effect in an attempt to preserve the initial linear profile.

We now show how these results can be applied to two real astrophysical problems.

## 1. Example: Dynamical friction on a galactic bar

As a galactic bar ages, it transfers angular momentum to its host galaxy and consequently its rotation rate slows. The mechanism responsible is dynamical friction: the bar produces a perturbation in the DF of the stars and dark matter particles that surround it (c.f. Fig. 24), and that perturbation back-reacts to produce a torque on the bar, draining its angular momentum. Problems of this kind have been the focus of many classic studies in galactic dynamics⁵¹,¹⁶⁶ and are strongly analogous to wave-particle interaction problems in plasma (e.g., Landau damping). Here, we illustrate the utility of our nonlinear pendulum kinetics, and the important role of the diffusion strength $\Delta$, by calculating the dynamical friction torque on a galactic bar in a simple model.
<figure><img src="image_44.png" /><figcaption>FIG. 31. Contours of g(φ, j, τ) [related to f through Eq. (146)], computed by solving the kinetic equation (142) numerically with the initial condition g(τ = 0) = j. Rows (a)–(c) show the time evolution of the solution for Δ = 0.001, 0.1, 1 respectively. The right column represents the steady-state solution. Adapted from Hamilton et al., Astrophys. J. 954, 12 (2023). Copyright 2023 IOP.¹⁵⁴</figcaption></figure>

We assume that the center of the bar is at the origin and that its long axis sits in the plane of an initially axisymmetric disk of stars. We let the bar rotate around its center in the φ direction with pattern speed Ωp and consider its gravitational interaction with the stars. One can write down a rather general formula for the torque T felt by the bar as follows. Let the bar have potential δΦ(r, t). Then, from Hamilton’s equation (32b), the specific torque felt by an individual star due to the bar is

$$ \frac{dL}{dt} = - \frac{\partial \delta \Phi}{\partial \theta_{\phi}}. \qquad (149) $$

Let the DF of stars be f (normalized such that ∫ dθφdθRdLdJf is the total mass in stars). Then, by Newton’s third law, the total torque on the bar is equal to

$$ T(t) = \int d\theta d\mathbf{J} f(\theta, \mathbf{J}, t) \frac{\partial \delta\Phi(\theta, \mathbf{J}, t)}{\partial \theta_{\phi}}. \qquad (150) $$

The challenge is to compute f for a given perturbation δΦ and then perform the integral (150).

The simplest way to calculate f is to resort to linear theory (and ignore diffusion) and use (94). In this case, putting δΦN(J, t) = δΦN(J) e⁻iNφΩp't and performing the integral over t' in Eq. (94), one
<figure><img src="image_45.png" /><figcaption>FIG. 32. Time evolution of the slow-angle-averaged DF, $\langle g \rangle_{\phi}$ (see Eq. (148)), for the same values of $\Delta$ as used in Fig. 31. Different panels correspond to different times, and the right-hand panel is the steady-state solution. The vertical dotted lines show the maximum extent of the separatrix, $j = \pm 2$ [see Eq. (141b)]. Adapted from Hamilton et al., Astrophys. J. 954, 12 (2023). Copyright 2023 IOP. $^{154}$</figcaption></figure>

finds a “linear torque,” $T^{\text{lin}} = \sum_{N} T_{N}^{\text{lin}}$, where the contribution from resonance N is

$$T_{N}^{\text{lin}}(t) = (2\pi)^{2}N_{\phi} \int d\mathbf{J} |\delta\Phi_{N}(\mathbf{J})|^{2} \mathbf{N} \cdot \frac{\partial f_{0}}{\partial \mathbf{J}} \frac{\sin[(\mathbf{N} \cdot \mathbf{\Omega} - N_{\phi}\mathbf{\Omega}_{\text{p}})t]}{\mathbf{N} \cdot \mathbf{\Omega} - N_{\phi}\mathbf{\Omega}_{\text{p}}}. \quad (151)$$

(We emphasize that in deriving this formula, we have ignored collective effects, because we have not accounted for the contribution to $\delta\Phi$ that arises self-consistently from the perturbed stellar distribution. While collective effects are typically very important for galaxy evolution, in the present case it turns out that they do not significantly alter the character of the response,$^{163}$ and we will continue to neglect them for simplicity). Taking the time-asymptotic limit $t \to \infty$, one arrives at the classic “LBK torque formula”$^{166}$

$$T_{N}^{\text{LBK}} = (2\pi)^{2}N_{\phi} \int d\mathbf{J} |\delta\Phi_{N}(\mathbf{J})|^{2} \mathbf{N} \cdot \frac{\partial f_{0}}{\partial \mathbf{J}} \pi\delta(\mathbf{N} \cdot \mathbf{\Omega} - N_{\phi}\mathbf{\Omega}_{\text{p}}). \quad (152)$$

The LBK formula (152) predicts that angular momentum is transferred to and from the stars exclusively at resonances. It is also directly analogous to the classic formula for the Landau damping rate of a Langmuir wave in an electrostatic plasma$^{167,168}$ or Landau damping in more general geometry.$^{100,169}$ Importantly, for many realistic DFs $f_0(\mathbf{J})$, the LBK torque is finite and negative, implying a long-term transfer of angular momentum away from the bar and hence a decay in its pattern speed. (As an example, at corotation the sign of the LBK torque is the sign of $\partial f_0 / \partial L$, which is typically negative whenever the density profile of the galaxy decreases with radius. This is similar to a plasma velocity distribution whose tail satisfies $\partial f_0 / \partial v < 0$). In practice, the time-asymptotic limit may not be valid since the torque can take several dynamical timescales to converge, by which time the pattern speed may have changed significantly.$^{170}$ Nevertheless, the LBK formula (152) is a good benchmark against which we can compare the magnitude of the torque arising from more sophisticated calculations.

Since we used the linearized solution (94) to compute $f$, we did not account for the nonlinear orbit trapping that can occur sufficiently close to each resonance, nor did we include any finite $\Delta$. We can

include these effects by redoing the calculation using our pendulum formalism (Sec. IV F). Converting to slow-fast angle-action variables $(\theta', J')$ around each resonance N, Fourier expanding the potential $\delta\Phi$ in slow angle space via $\delta\Phi(\mathbf{r}, t) = \sum_k \Psi_k(J')e^{ik\theta'}$ as in Eq. (63), and expanding the DF as $f = \sum_k f_k(J', t)e^{ik\theta'}$, one can show that the total “nonlinear torque” on the bar $T^{\text{non-lin}}$ can be written as a sum over N of contributions (see Hamilton et al. $^{154}$ for details$^{171}$)

$$T_{N}^{\text{non-lin}}(t) = 2(2\pi)^{2}N_{\phi} \sum_{k>0} k \int d\mathbf{J}' \operatorname{Im}[f_{k}(J', t) \Psi_{k}^{*}(J')]. \quad (153)$$

Here, we have used the shorthand $f_{(0,k)} = f_k$, and similarly for $\Psi_k$. In using expressions like (153), one must be careful to properly divide phase space into nonoverlapping sub-volumes that are dominated by individual resonances N (see Sec. 4.5 of Chiba and Schönrich$^{52}$ for a discussion).

In reality, the torque on the bar involves a sum over all possible resonances N and an integral over all of action space, and each such contribution comes with a slightly different value of $\Delta$. In order to isolate a representative quantity that only depends on a single value of $\Delta$, we will focus henceforth on the corotation torque density $S(t)$, which is the contribution to the torque only arising from the corotation resonance $\Omega_{\phi} = \Omega_{\text{p}}$ (so $N_{\phi} = 2$ and $k = 1$, see Sec. IV F 1) and only considering stars with a particular $J_R$, namely,

$$S(t)J_R = 4(2\pi)^2 \int d\mathbf{L} \Psi_1 \operatorname{Im} f_1. \quad (154)$$

For a given $\Delta$, we can compute $\operatorname{Im} f_1$ using the numerical solution to the kinetic equation (142) that we showed in Fig. 31. We will compare this to the benchmark LBK value [see Eq. (152)], which is the result of computing the torque density using linear theory with $\Delta = 0$:

$$S_{\text{LBK}}(J_R) = 4(2\pi)^2 \int d\mathbf{L} |\Psi_1|^2 \frac{\partial f_0}{\partial \mathbf{L}} \pi \delta[2(\Omega_{\phi} - \Omega_{\text{p}})]. \quad (155)$$

In Fig. 33(a), we show the bar’s dimensionless slowing rate $S(t)/S_{\text{LBK}}$ as a function of time for various values of $\Delta \ge 0$, fixing $J_R = 50 \text{ kpc}^{-2} \text{ Gyr}^{-1}$ (which is much less than $L_{\text{CR}}$, meaning the stars are initially on roughly circular orbits—see Fig. 9).
Note that $S_{LBK}$ is negative, so that a positive value of $S(t)/S_{LBK}$ means the bar feels a negative torque (slowing it down). In the collisionless case ($\Delta = 0$, black dotted line), we see that the slowing rate oscillates on the timescale $\approx t_{lib} \approx 2$ Gyr, and its envelope also decays over time. In the time-asymptotic limit for $\Delta = 0$, we find

$$S(t \to \infty) = 0, \quad (\Delta = 0). \qquad (156)$$

This zero torque result is a consequence of the angle-symmetry of the stellar density distribution ($\partial_{\phi}f = 0$ at fixed $j$) that arises when the DF completely phase mixes within and around the resonant island [Fig. 31(a)]. Simply put, in the fully phase-mixed state, there are the same number of stars “pushing” on the bar as “pulling” on it. This is the direct analog of O’Neil⁹⁶’s famous result in plasma theory: when nonlinear trapping is included (and diffusion neglected), the Landau damping rate of an electrostatic wave asymptotes to zero, because of the symmetrization of the “bouncing” particle distribution in the trough of the wave. However, as one deviates from $\Delta = 0$ to finite $\Delta > 0$, the behavior changes. The slowing rate in the steady state is manifestly positive for all $\Delta > 0$. This is a consequence of the fact that while phase-mixing attempts to abolish the asymmetry of the slow-angle distribution, diffusion replenishes it [e.g., Fig. 31(b)] leading to a finite $Imf_1$ in Eq. (154). Moreover, for $\Delta \gtrsim 1$, the slowing rate no longer oscillates on the libration timescale, because no trapped star is able to complete a full libration before being de-trapped by diffusion.

To make these statements more quantitative, we turn to panel (b) of Fig. 33, in which we show the time-asymptotic slowing rate $S(t \to \infty)$ as a function of $\Delta$. We see that $S(t \to \infty)$ grows almost linearly with $\Delta$ until $\Delta \approx 1$, after which it saturates at a level similar to that given by the naive LBK theory, even though LBK theory ignores both diffusion and nonlinearity! This is as expected since, roughly speaking, for any $\Delta \gtrsim 1$ there is no nonlinear trapping and the effect of “diffusion” is just to broaden the delta-function in Eq. (152). That is, diffusion relinearizes the dynamics. Again, this effect was also appreciated in the plasma theory of electric wave damping: if small-scale collisions are able to knock particles out of the wave trough before they can

complete a full bounce orbit, then O’Neil’s nonlinear trapping effects are interrupted, and the standard (quasi)-linear Landau damping/growth rate is recovered.⁵⁹,⁶⁰

These results suggest that real galactic bars will always feel some nonzero torque, since finite diffusion will always replenish some asymmetry in the angular distribution of the stars as viewed in the bar frame. They also suggest that for systems with moderate to strong diffusion ($\Delta \gtrsim 1$), the LBK formula (152) provides a robust order-of-magnitude estimate of the frictional torque. Moreover, we have treated the bar’s pattern speed, $\Omega_p$, as constant here, but very similar arguments apply to the case where the “interruption” of trapping is not due to collisions (i.e., finite $\Delta$) but rather due to the fact that a changing $\Omega_p$ causes resonance locations to sweep through phase space. In that case, the LBK formula holds in what Tremaine and Weinberg⁵¹ called the “fast regime,” namely, when the bar slows down rapidly enough that the resonances move by their own width $I_h$ on a timescale shorter than $t_{lib}$. In the opposite “slow regime,” trapped stars are dragged to larger radii as the corotation resonance moves outwards, a process that seems to have occurred in the Milky Way.¹⁵⁰ A general theory that connects smoothly the fast (LBK or Landau-like) regime and the slow (nonlinear, O’Neil-like) regime was recently completed by Chiba.⁵⁸

## 2. Example: Saturation of spiral instabilities

Finally, let us discuss how pendulum kinetics can be used to understand the saturation mechanism of spiral instabilities and predict the resulting saturation amplitude.

We already saw (Sec. VI B 2) that spiral waves can arise in stellar disks as exponentially growing, linear instabilities of the disk’s underlying DF $f_0$, particularly in cases where the disk is sufficiently cold/dense (Fig. 21) or grooved (Fig. 22). This exponential growth phase cannot continue forever. Simulations by Sellwood and Carlberg¹⁵¹ and theoretical work by Hamilton¹⁷² suggest that the saturation occurs because of orbit trapping at the corotation resonance. To understand their argument, we present Fig. 34 (adapted from Hamilton¹⁷²).

<figure><img src="image_46.png" /><figcaption>FIG. 33. Adapted from C. Hamilton, E. A. Tolman, L. Arzamasskiy, and V. N. Duarte, Astrophys. J. 954, 12 (2023). Copyright 2023 IOP.¹⁵⁴ (a) Dimensionless torque on the bar (i.e., its “slowing rate”) $S(t)/S_{LBK}$ due to stars at the corotation resonance, with $J_R = 50$ kpc$^2$ Gyr$^{-1}$ (i.e., roughly circular orbits). The libration time is $t_{lib} \approx 2$ Gyr. Colored lines show results for different $\Delta > 0$, while the dotted line shows the collisionless ($\Delta = 0$) result. (b) Steady-state value of $S(t)/S_{LBK}$ as a function of $\Delta$ (note that both axes are on a logarithmic scale, and the different lines just correspond to different choices of integration limits).</figcaption></figure>
<figure><img src="image_47.png" /><figcaption>FIG. 34. Reproduced from Hamilton, Mon. Not. R. Astron. Soc. 528, 5286 (2024). Copyright 2024 Royal Astronomical Society.¹⁷² (a) Growth and saturation of a spiral mode in an N-body simulation of a razor thin stellar disk (adapted from Ref. 151 green line). The gray line illustrates the simplification used in this paper: we separate the spiral's life into a linear phase of perfectly exponential growth, and a nonlinearly saturated phase of constant amplitude. (b) Illustration of the change to the axisymmetric DF $f_0(L, t)$ in the vicinity of a corotation resonance. The initial instability is driven by some sharp feature in the DF (here a groove, centered on $L = L_*$ with width $w_L$, see Sellwood and Carlberg¹⁵¹). In the saturated state, the DF is flattened around the resonance $L = L_{CR}$.</figcaption></figure>

The green line in panel (a) shows the logarithmic amplitude of a spiral density wave with constant pattern speed $\Omega_p$, from an N-body simulation of an initially grooved, razor thin disk. Clearly, once some noisy transients have faded away the mode undergoes a very well-defined exponential growth, confirming that it is a linear mode of the initial DF. Eventually, though, the exponential growth stops and the mode settles at a roughly constant amplitude. The gray line illustrates schematically the simplification we will use here, which is that the transition from the linear phase (exponentially growing amplitude) to the nonlinear phase (constant amplitude) happens instantaneously. In panel (b), we sketch the corresponding DF of angular momenta in the disk, $f_0(L)$, for circular orbits ($J_R = 0$). The initial unstable distribution is shown with a black line. The idea is that once the mode has grown to large enough amplitude that it starts to resonantly trap stellar orbits at corotation, the DF in the vicinity of that resonance will be flattened (as shown with the red line), just as we saw in Sec. VIII B 1. This flattening erases the sharp feature in $f_0$ that originally drove the instability.

Hamilton¹⁷² showed that by calculating the angular momentum content of the spiral in the linear (i.e., exponentially growing) and nonlinear (i.e., fully saturated) regimes, and equating the two expressions, one can read off an approximate formula for the saturation amplitude of the spiral potential. This is

$$|\delta\Phi_m| = \left(\frac{128}{9\pi^2}\right)^2 \frac{2}{m} \frac{\beta^2}{|\partial\Omega_\phi/\partial L|_{L_{CR}}}, \quad (157)$$

where $\beta = \text{Im}\,\omega_m$ is the mode's growth rate. This result can be made more physically enlightening if we use (66) and (67) to express it as

$$\omega_{\text{lib}} = \frac{256}{9\pi^2} m^{1/2} \beta \approx 2.88 m^{1/2} \beta. \quad (158)$$

In other words, saturation occurs when the nonlinear libration frequency $\omega_{\text{lib}}$ is a few times larger than the linear growth rate $\beta$ (so nonlinear orbit distortions occur before the mode amplitude can grow significantly). Note that for one-armed spirals, Eq. (158) is identical to

Eq. (12) of Dewar,¹⁷³ who studied the saturation of linear momentum transfer between electrons and a single longitudinal plasma wave.

One can also relate (157) to the surface density perturbation of the spiral $\delta\Sigma_m$ relative to the background stellar disk $\Sigma_0$,

$$\begin{aligned} \left|\frac{\delta\Sigma_m}{\Sigma_0}\right| &= \frac{32768}{81\pi^4} \xi^{-1} \cot p \left(\frac{\beta}{\Omega_p}\right)^2 \\ &\approx 0.5 \times \left(\frac{\xi}{0.5}\right)^{-1} \left(\frac{\cot p}{1.5}\right) \left(\frac{\beta/\Omega_p}{0.2}\right)^2, \end{aligned} \quad (159)$$

where $\xi$ is the active fraction of the disk (Sec. VI B 2) and $p$ is the spiral's pitch angle. Thus, spiral instabilities tend to saturate once their amplitudes reach several tens of percent of the background density.

Last of all, it is enlightening to estimate the timescale $t_{\text{sat}}$ over which saturation of realistic spirals might be achieved. To do this, first recall [Eq. (158)] that saturation occurs when $\omega_{\text{lib}}(t_{\text{sat}}) \sim \beta$. Second, a basic scaling of Eq. (67) tells us that at corotation, $\omega_{\text{lib}}(t) \sim \varepsilon(t)^{1/2}\Omega_p$, where $\varepsilon(t) \sim |\Phi_m(t)/H_0|$ is the dimensionless amplitude of the mode. Third, in the linear regime we have $\varepsilon(t) = \varepsilon(0)e^{\beta t}$. Putting these three ingredients together, we get a saturation time

$$t_{\text{sat}} \sim \frac{2}{\beta} \left[ \ln \left( \frac{\beta}{\Omega_p} \right) + \frac{1}{2} \ln \left( \frac{1}{\varepsilon(0)} \right) \right]. \quad (160)$$

In particular, if the initial fluctuation level is set by Poisson noise alone, then $\varepsilon(0) \sim N^{-1/2}$, and we have (very roughly)

$$\begin{aligned} t_{\text{sat}} &\sim \frac{1}{2\beta} \ln N, \\ &\sim 2.5 \text{ Gyr} \times \left( \frac{T_p}{250 \text{ Myr}} \right) \left( \frac{\beta/\Omega_p}{0.2} \right)^{-1} \left( \frac{\ln N}{25} \right), \end{aligned} \quad (161)$$

where $T_p = 2\pi/\Omega_p$ is the orbital period at corotation. Thus, if it has to exponentiate out of Poisson noise, a rather vigorously growing spiral with $\beta/\Omega_p \gtrsim 0.1$ will saturate after a few Gyr. More realistic galactic
noise tends to be of higher amplitude than is implied by Poisson statistics (Sec. VII C 3), though this will only change the timescale $t_{sat}$ by a logarithmic factor.

## IX. FURTHER READING

There are many more connections between plasma theory and stellar dynamics that we have not touched on in this tutorial article. Let us mention just a few of them here, alongside references for further reading.

First, the most basic connection between the two fields is the theory of orbits in spatially smooth and constant or time-varying potentials, the regularity/chaoticity of those orbits, and the associated theories of canonical transformations, Hamiltonian perturbation theory, and guiding center motion. These are classic problems in the general theory of dynamics,³⁹ and Hamiltonian perturbation theory in angle-action variables has long been a tool of theorists (see Chap. 3 of BT08), but it is only relatively recently that angle-action variables have been adopted across mainstream galactic astronomy.¹⁷⁴ They are now ubiquitous tools with which astronomers study phase space data of stellar systems, particularly from the Gaia satellite.⁴³ Chaotic orbits are known to exist in triaxial (i.e., nonaxisymmetric) galaxies and star clusters (see, e.g., the books by Contopoulos,¹⁷⁵ Merritt¹⁷⁶) but are mostly ignored in kinetic theories, which usually focus on the evolution of near-axisymmetric or near-spherical systems. Guiding center and gyrokinetic theory have close analogues in stellar dynamics, particularly for describing stellar motion near supermassive black holes, where the orbits are closed Keplerian ellipses, which can be averaged over to form Keplerian “rings” or “wires”.¹²⁴,¹⁷⁷ Despite this, the Hamiltonian theories of guiding center motion, averaged Lagrangians, ponderomotive effects, and so on developed by the likes of Cary and Brizard,¹⁷⁸ Dodin and Fisch¹⁷⁹ in plasma theory have not (yet) received much attention in galactic dynamics.

Second, another common thread between plasma kinetics and stellar dynamics lies in the construction of collisionless equilibria and their stability properties. Though we have barely touched this subject here, developing stable candidate distribution functions is a crucial part of stellar-dynamical modeling,

completely analogous to the construction of Vlasov-Poisson or Vlasov-Maxwell¹⁸⁰ equilibria in plasma. Similarly, Nyquist diagrams¹⁸¹ and Penrose’s criterion¹⁸² are just as applicable to stellar systems as they are to plasmas. On the gravitational side, perhaps the best introductions to this subject (written in a language plasma physicists can appreciate) are to be found in Secs. 4 and 5 of BT08 and the 1994 book by Palmer.¹⁸³ The theory of *nonlinear* stability, formation of BGK-like modes,¹⁸⁴ and so on have been developed to some extent in stellar dynamics,¹⁸⁵⁻¹⁸⁷ usually with direct reference to previous plasma results.

Third, we mention a major area of stellar-dynamical kinetics that we have only touched on here, the theory of spiral structure. This was one of the original areas in which the galactic plasma overlap was found to be most fruitful,⁷⁵,⁷⁷,¹¹⁵,¹⁶⁶,¹⁸⁸ since it concerns problems familiar to plasma practitioners such as wave propagation, WKB theory, wave-particle interaction, phase mixing, and so on. Despite this, the theory is far from complete: opinions still differ over the mechanisms that provoke spiral responses, the lifetimes of individual spiral patterns, the importance of gas flows and star formation, and so on—for an array of perspectives, see, e.g., reviews by Athanassoula,¹⁸⁹ Dobbs and Baba,¹⁹⁰ and Sellwood and Masters.⁹³

Finally, while in this article we have focused on isolated and rather quiescent stellar systems, long gone are the days when astronomers thought of galaxies as “island universes.” In the modern view, galaxies are continually bombarded by infalling dark matter subhaloes and dwarf galaxies¹⁹¹ and are connected to other galaxies via the cosmic web.¹⁹² Nor can one always think of “stellar” (or even stellar + dark matter) components of galaxies in isolation without also considering their gaseous components. These give rise to star formation and provide significant dynamical feedback via active galactic nuclei¹⁹³ and via turbulence, shocks, supernova bubbles, etc., in the interstellar gaseous medium (see Fig. 35 for illustration). Galaxies are therefore inherently multi-scale (both spatially and temporally), multi-component, open, dissipative systems that are invariably found out of—and often far from—equilibrium. Kinetic theories that are able to deal with such complexity are in their infancy.⁶⁵

<figure><img src="image_48.png" /><figcaption>**FIG. 35.** The Phantom Galaxy. Credit: ESA/ Webb, NASA & CSA, J. Lee and the PHANGS-JWST Team.</figcaption></figure>
We hope that our plasma colleagues consider this tutorial article an open invitation to participate in a subject that is still in its theoretical adolescence; in which simulations with realistic numbers of particles are routine;¹⁹⁴ in which observations are astonishing and improving all the time; and, most of all, in which the Universe constantly impresses upon us its complex, kaleidoscopic beauty.

## ACKNOWLEDGMENTS

This article grew out of the notes that we each compiled as lecturers of the Kinetic Theory of Stellar Systems MMathPhys course at the University of Oxford. We thank Professor James Binney, who designed the course and was its first lecturer, as well as Professor Paul Dellar and Professor Alexander Schekochihin who taught the sister courses on kinetic theory of gases and plasmas,¹⁹⁵ respectively, and we gratefully acknowledge the hospitality of Merton College, Oxford. We also thank Uddipan Banik, James Binney, Thomas Foster, Shaunak Modak, Michael Nastac, Christophe Pichon, Frank van den Bosch, Martin Weinberg, George Wong, and Tomer Yavetz for feedback on the manuscript. C. H. is supported by the John N. Bahcall Fellowship Fund at the Institute for Advanced Study. This research was supported in part by Grant No. NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Chris Hamilton:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal).

**Jean-Baptiste Fouvry:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## APPENDIX A: THE EPICYCLIC APPROXIMATION

In this Appendix, we derive the angle-action variables for near-circular orbits in a central potential $\Phi(R)$, using the so-called epicyclic approximation. To do this, let us first write down the Hamiltonian for radial motion in this potential

$$H(R, v_R) = \frac{1}{2} v_R^2 + \Phi_{\text{eff}}(R), \qquad (A1)$$

where $\Phi_{\text{eff}}(R) \equiv \Phi(R) + L^2/(2R^2)$ is the effective potential and $L = R^2 \phi$ is a constant. The Hamiltonian equations of motion then give

$$\frac{d^2 R}{dt^2} = - \frac{\partial \Phi_{\text{eff}}(R)}{\partial R}. \qquad (A2)$$

Perfectly, circular motion corresponds to the right hand side of this equation being set to zero, i.e., a minimum of the effective potential. To work out the lowest-order correction to this circular motion, we can expand $\Phi_{\text{eff}}(R)$ around $R_g$ and find

$$\frac{d^2(R - R_g)}{dt^2} = -\Omega_R^2(R - R_g) + \dots, \qquad (A3)$$

where $\Omega_R^2 = (\partial^2\Phi_{\text{eff}}/\partial R^2)|_{R_g}$. This just corresponds to a harmonic oscillator with angular frequency $\Omega_R$. This radial motion may be alternatively described using a pair of radial angle-action coordinates $(\theta_R, J_R)$ satisfying

$$J_R = \frac{\dot{R}^2 + \Omega_R^2(R - R_g)^2}{2\Omega_R} \equiv \frac{E_R}{\Omega_R}, \qquad (A4a)$$

$$R = R_g - a_R \cos \theta_R, \qquad (A4b)$$

where the amplitude of the radial epicycle is

$$a_R \equiv \frac{\sqrt{2E_R}}{\Omega_R}. \qquad (A5)$$

(Note that we have chosen our angle coordinate $\theta_R$ such that $\theta_R = 0$ corresponds to the minimum $R$, i.e., the pericenter $R = R_a$. Other authors often choose different conventions.) Similarly, the azimuthal part of the problem can also be described by angle-action variables

$$L = R_g^2 \Omega_{\phi}, \qquad (A6a)$$

$$\phi = \theta_{\phi} + \frac{2\Omega_{\phi}a_R}{\Omega_R R_g} \sin \theta_R, \qquad (A6b)$$

and the frequencies can be written in the following form in terms of $\Phi$, all evaluated at $R = R_g$:

$$\Omega_{\phi}^2 = \frac{1}{R} \frac{\partial \Phi}{\partial R} \equiv \Omega_c^2, \qquad (A7a)$$

$$\Omega_R^2 = \frac{2\Omega_{\phi} \partial(\Omega_{\phi} R^2)}{R \partial R} \equiv \kappa^2. \qquad (A7b)$$

One can also extend the epicyclic approximation to include vertical harmonic oscillations or higher order corrections, but this will not concern us here. See Sec. 3.2.3 of BT08 and Dehnen¹⁹⁶ for more details.

## APPENDIX B: THE BASIS METHOD

To solve Eq. (114), we assume that we have at our disposal a complete set of biorthogonal basis elements $(\Phi^{(p)}(\mathbf{r}), \rho^{(p)}(\mathbf{r}))$, which satisfy

$$\Phi^{(p)}(\mathbf{r}) = \int d\mathbf{r}' \psi(\mathbf{r}, \mathbf{r}') \rho^{(p)}(\mathbf{r}'), \qquad (B1a)$$

$$\int d\mathbf{r} \Phi^{(p)*}(\mathbf{r}) \rho^{(q)}(\mathbf{r}) = -\mathcal{E} \delta_p^q, \qquad (B1b)$$

where $\psi(\mathbf{r}, \mathbf{r}') = -G/|\mathbf{r} - \mathbf{r}'|$ and $\mathcal{E}$ is an arbitrary constant with units of energy. The “completeness” requirement is that we can
expand an arbitrary function $\Psi(\mathbf{r}) = \sum_p A_p \Phi^{(p)}(\mathbf{r})$ for some set of constant coefficients $A_p$. Various such basis sets do exist, and the particular set you should choose depends on the geometry of the system you are interested in. For instance, in a spherical geometry, one might take the basis to consist of products of radial Bessel functions and angular spherical harmonics (see, e.g., Kalnajs,²⁷ Petersen et al.⁹¹ for more details). The idea is to project (114) onto this basis and then read off the coefficients to get an explicit expression for $\psi^d$.

Thus, let us begin by expanding the pairwise interaction potential $\psi(\mathbf{r}, \mathbf{r}') = -G/|\mathbf{r} - \mathbf{r}'|$ for a fixed value of $\mathbf{r}'$ as

$$ \psi(\mathbf{r}, \mathbf{r}') = \sum_p u_p(\mathbf{r}') \Phi^{(p)}(\mathbf{r}), \quad (B2) $$

where the coefficient $u_p(\mathbf{r}')$ has units of (mass)$^{-1}$. To find this coefficient, we multiply both sides of this equation by $\rho^{(q)*}(\mathbf{r})$, integrate over $\mathbf{r}$, and use both properties (B1). The result is

$$ \begin{aligned} u_p(\mathbf{r}') &= -\frac{1}{\mathcal{E}} \int d\mathbf{r} \psi(\mathbf{r}, \mathbf{r}') \rho^{(p)*}(\mathbf{r}) \\ &= -\frac{1}{\mathcal{E}} \Phi^{(p)*}(\mathbf{r}'). \end{aligned} \quad (B3) $$

Plugging this back in to Eq. (B2), we get

$$ \psi(\mathbf{r}, \mathbf{r}') = -\frac{1}{\mathcal{E}} \sum_p \Phi^{(p)}(\mathbf{r}) \Phi^{(p)*}(\mathbf{r}'). \quad (B4) $$

The fact that the bare interaction can be expanded in this “quasi-separable” form is one of the main advantages of the basis method.

What we really need is the Fourier transform of this function

$$ \psi_{nn'}(\mathbf{J}, \mathbf{J}', \omega) = -\frac{1}{\mathcal{E}} \sum_p \Phi_n^{(p)}(\mathbf{J}) \Phi_{n'}^{(q)*}(\mathbf{J}'), \quad (B5) $$

where $\Phi_n^{(p)}(\mathbf{J})$ stands for the Fourier transformed basis elements naturally defined following Eq. (37). Thus, we now have an explicit expression for the functions $\psi_{nn'}(\mathbf{J}, \mathbf{J}')$, which enter the right-hand side of equation (99).

We now perform a very similar expansion of the dressed interaction $\psi^d$. By analogy with (B5), we assume that, at fixed $\omega$, this function can be expanded as

$$ \psi_{nn'}^d(\mathbf{J}, \mathbf{J}', \omega) = -\frac{1}{\mathcal{E}} \sum_{p,q} \Phi_n^{(p)}(\mathbf{J}) E_{pq}^{-1}(\omega) \Phi_{n'}^{(q)*}(\mathbf{J}'), \quad (B6) $$

where $E(\omega)$ is a dimensionless, as-yet-unknown, matrix which depends only on $\omega$ (and we have assumed its inverse exists). Injecting this decomposition into Eq. (114), we obtain

$$ \sum_{p,q} \Phi_n^{(p)}(\mathbf{J}) E_{pq}^{-1}(\omega) \Phi_{n'}^{(q)*}(\mathbf{J}') = \sum_{p,q} \Phi_n^{(p)}(\mathbf{J}) [\mathbf{M}(\omega) \mathbf{E}^{-1}(\omega)]_{pq} \Phi_{n'}^{(q)*}(\mathbf{J}') \\ + \sum_{p,q} \Phi_n^{(p)}(\mathbf{J}) l_{pq} \Phi_{n'}^{(q)*}(\mathbf{J}'), \quad (B7) $$

where we introduced the identity matrix $I$, as well as the system's dimensionless response matrix $M(\omega)$ as in Eq. (101). Since the basis is biorthogonal, we can identify the elements in Eq. (B7), and we obtain the matrix relation

$$ E^{-1}(\omega) = M(\omega) E^{-1}(\omega) + I, \quad (B8) $$

which immediately tells us that

$$ E(\omega) = I - M(\omega). \quad (B9) $$

To conclude, from Eq. (B6), we have

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

Thus, we have done what we anticipated in Sec. VI A, which is essentially to invert an operator $I - M$, with $M$ encoding collective effects. We can easily recover the bare interaction, Eq. (B5), by ignoring collective effects, i.e., by setting $M(\omega) = 0$ in Eq. (B10).

## REFERENCES

¹J. Binney and M. Merrifield, Galactic Astronomy (Princeton University Press, 1998).

²J. Falcón-Barroso and J. H. Knapen, Secular Evolution of Galaxies (Cambridge University Press, 2012).

³H. Mo, F. Van den Bosch, and S. White, Galaxy Formation and Evolution (Cambridge University Press, 2010).

⁴J. Binney and S. Tremaine, Galactic Dynamics: Second Edition (Princeton University Press, 2008).

⁵See http://www.ch-astro.com/teaching for the associated lecture materials.

⁶T. Alexander, Annu. Rev. Astron. Astrophys. 55, 17 (2017).

⁷D. Lynden-Bell, Mon. Not. R. Astron. Soc. 136, 101 (1967).

⁸R. J. Ewart, A. Brown, T. Adkins, and A. A. Schekochihin, J. Plasma Phys. 88, 925880501 (2022).

⁹A. Campa, T. Dauxois, D. Fanelli, and S. Ruffo, Physics of Long-Range Interacting Systems (OUP Oxford, 2014).

¹⁰S. Tremaine, Astrophys. J. 807, 157 (2015).

¹¹J. G. Hills and C. A. Day, Astrophys. Lett. 17, 87 (1976), see https://ui.adsabs. harvard.edu/abs/1976ApL....17...87H/abstract.

¹²S. Chandrasekhar, Astrophys. J. 97, 255 (1943).

¹³E. M. Lifshitz and L. P. Pitaevskii, Physical Kinetics (Elsevier, 1981).

¹⁴D. C. Heggie and P. Hut, The Gravitational Million-Body Problem (Cambridge University Press, 2003).

¹⁵Actually, the argument we have made here does not really apply to the geometry of a thin disk. Accounting for the disk geometry, the estimate on the right hand side of (19) is reduced by a factor like $\sim 10^4$ [see Eq. (5) of Sellwood⁹⁵ and the surrounding discussion]. Nevertheless, even this shortened two-body relaxation time is much, much longer than the age of the Universe.

¹⁶P. H. Diamond, S.-I. Itoh, and K. Itoh, Modern Plasma Physics (Cambridge University Press, Cambridge, 2010).

¹⁷P.-H. Chavanis, Eur. Phys. J. Plus 128, 128 (2013).

¹⁸B. Moore, S. Gelato, A. Jenkins, F. R. Pearce, and V. Quilis, Astrophys. J. 535, L21 (2000).

¹⁹K. J. Joshi, F. A. Rasio, and S. Portegies Zwart, Astrophys. J. 540, 969 (2000).

²⁰M. Hénon, Astrophys. Space Sci. 14, 151 (1971).

²¹N. Frankel, J. Sanders, Y.-S. Ting, and H.-W. Rix, Astrophys. J. 896, 15 (2020).

²²M. Hénon, in Dynamics of the Solar Systems, edited by A. Hayli (IAU, 1975), Vol. 69, p. 133; see https://ui.adsabs.harvard.edu/abs/1975IAUS...69..133H/abstract.

²³M. Giersz and D. C. Heggie, Mon. Not. R. Astron. Soc. 268, 257 (1994).

²⁴J.-B. Fouvry, C. Hamilton, S. Rozier, and C. Pichon, Mon. Not. R. Astron. Soc. 508, 2210 (2021).

²⁵D. Lynden-Bell and R. Wood, Mon. Not. R. Astron. Soc. 138, 495 (1968).

²⁶A. H. Boozer, Rev. Mod. Phys. 76, 1071 (2004).

²⁷F. Zonca, L. Chen, S. Briguglio, G. Fogaccia, G. Vlad, and X. Wang, New J. Phys. 17, 013052 (2015).

²⁸L. Chen and F. Zonca, Rev. Mod. Phys. 88, 015008 (2016).

²⁹J. Bovy, Astrophys. J., Suppl. Ser. 216, 29 (2015).

³⁰A. M. Price-Whelan, J. Open Source Software 2, 388 (2017).
³¹E. Vasiliev, *Mon. Not. R. Astron. Soc.* **482**, 1525 (2019).

³²J. Binney and D. Spergel, *Astrophys. J.* **252**, 308 (1982).

³³J. Binney and D. Spergel, *Mon. Not. R. Astron. Soc.* **206**, 159 (1984).

³⁴M. Hénon and C. Heiles, *Astrophys. J.* **69**, 73 (1964).

³⁵H. Goldstein, *Classical Mechanics* (Addison-Wesley, 1950).

³⁶M. D. Weinberg, *Astrophys. J.* **108**, 1398 (1994).

³⁷M. D. Weinberg, *Astrophys. J.* **108**, 1403 (1994).

³⁸M. D. Weinberg, *Astrophys. J.* **108**, 1414 (1994).

³⁹A. J. Lichtenberg and M. A. Lieberman, *Regular and Chaotic Dynamics* (Springer, 2013), Vol. 38.

⁴⁰V. I. Arnold, *Mathematical Methods of Classical Mechanics* (Springer, 1989).

⁴¹J. A. Bittencourt, *Fundamentals of Plasma Physics* (Springer, 2013).

⁴²J. A. S. Hunt, M. W. Bub, J. Bovy, J. T. Mackereth, W. H. Trick, and D. Kawata, *Mon. Not. R. Astron. Soc.* **490**, 1026 (2019).

⁴³W. H. Trick, J. Coronado, and H.-W. Rix, *Mon. Not. R. Astron. Soc.* **484**, 3291 (2019).

⁴⁴J. L. Sanders and J. Binney, *Mon. Not. R. Astron. Soc.* **457**, 2107 (2016).

⁴⁵M. D. Weinberg, *Mon. Not. R. Astron. Soc.* **239**, 549 (1989).

⁴⁶K. Kaur and S. Sridhar, *Astrophys. J.* **868**, 134 (2018).

⁴⁷U. Banik and F. C. van den Bosch, *Astrophys. J.* **912**, 43 (2021).

⁴⁸U. Banik and F. C. van den Bosch, *Astrophys. J.* **926**, 215 (2022).

⁴⁹K. Kaur and N. C. Stone, *Mon. Not. R. Astron. Soc.* **515**, 407 (2022).

⁵⁰J. Binney, *Mon. Not. R. Astron. Soc.* **495**, 895 (2020).

⁵¹S. Tremaine and M. D. Weinberg, *Mon. Not. R. Astron. Soc.* **209**, 729 (1984).

⁵²R. Chiba and R. Schönrich, *Mon. Not. R. Astron. Soc.* **513**, 768 (2022).

⁵³A.-C. Eilers, D. W. Hogg, H.-W. Rix, N. Frankel, J. A. S. Hunt, J.-B. Fouvry, and T. Buck, *Astrophys. J.* **900**, 186 (2020).

⁵⁴J. A. Sellwood and J. Binney, *Mon. Not. R. Astron. Soc.* **336**, 785 (2002).

⁵⁵J. M. Petit and M. Henon, *Icarus* **66**, 536 (1986).

⁵⁶M. Hasegawa and K. Nakazawa, *Astron. Astrophys.* **227**, 619 (1990), see https://ui.adsabs.harvard.edu/abs/1990A%26A...227.619H/abstract.

⁵⁷P. A. Sturrock, *Phys. Rev.* **141**, 186 (1966).

⁵⁸R. Chiba, *Mon. Not. R. Astron. Soc.* **525**, 3576 (2023).

⁵⁹G. L. Johnston, *Phys. Fluids* **14**, 2719 (1971).

⁶⁰S. P. Auerbach, *Phys. Fluids* **20**, 1836 (1977).

⁶¹V. N. Duarte, N. N. Gorelenkov, R. B. White, and H. L. Berk, *Phys. Plasmas* **26**, 120701 (2019).

⁶²E. A. Tolman and P. J. Catto, *J. Plasma Phys.* **87**, 855807201 (2021).

⁶³Note the subtlety associated with self-interactions, i.e. contributions from $\psi(r_i, r_j)$. As long as one assumes that there are no self-forces, this is fine.

⁶⁴Equation (85) is usually called the Collisionless Boltzmann Equation in stellar dynamics, or the Vaslov equation in plasma kinetics. But in research papers one is equally likely to see either of those names given to the Klimontovich equation (80), or to something else entirely, so caution is required.

⁶⁵C. Pichon and D. Aubert, *Mon. Not. R. Astron. Soc.* **368**, 1657 (2006).

⁶⁶B. B. Kadomtsev, *Plasma Turbulence* (Academic Press, New York, 1965).

⁶⁷T. H. Dupree, *Phys. Fluids* **9**, 1773 (1966).

⁶⁸J. A. Krommes, *Phys. Rep.* **360**, 1–352 (2002).

⁶⁹J. A. Krommes, *J. Plasma Phys.* **81**, 205810601 (2015).

⁷⁰R. W. Gould, T. M. O’Neil, and J. H. Malmberg, *Phys. Rev. Lett.* **19**, 219 (1967).

⁷¹K. F. Riley, M. P. Hobson, and S. J. Bence, *Mathematical Methods for Physics and Engineering* (Cambridge University Press, 2006).

⁷²L. Al’tshul and V. Karpman, *JETP* **20**, 1043 (1965), see http://www.jetp.ras.ru/cgi-bin/e/index/e/20/4/p1043?a=list.

⁷³A. J. Layzer, *Phys. Rev.* **129**, 897 (1963).

⁷⁴R. Bälescu, *J. Math. Phys.* **4**, 1009 (1963).

⁷⁵A. J. Kalnajs, *Astrophys. J.* **166**, 275 (1971).

⁷⁶T. Kato, *Perturbation Theory for Linear Operators* (Springer, 2013).

⁷⁷A. J. Kalnajs, *Astrophys. J.* **205**, 745 (1976).

⁷⁸J. Barré, A. Olivetti, and Y. Y. Yamaguchi, *J. Phys. A* **44**, 405502 (2011).

⁷⁹J. Barré and Y. Y. Yamaguchi, *J. Phys. A* **46**, 225501 (2013).

⁸⁰D. B. Melrose, *Instabilities in Space and Laboratory Plasmas* (Cambridge University Press, 1986).

⁸¹J. Y. Kim, Y. Kishimoto, W. Horton, and T. Tajima, *Phys. Plasmas* **1**, 927 (1994).

⁸²T. Kuroda, H. Sugama, R. Kanno, M. Okamoto, and W. Horton, *J. Phys. Soc. Jpn.* **67**, 3787 (1998).

⁸³S. D. Baalrud, *Phys. Plasmas* **20**, 012118 (2013).

⁸⁴J.-B. Fouvry and S. Prunet, *Mon. Not. R. Astron. Soc.* **509**, 2443 (2022).

⁸⁵A. Toomre and A. J. Kalnajs, in *Dynamics of Disc Galaxies* (Göteborgs Univ. and Chalmers Univ. of Technology 1991), p. 341; see https://ui.adsabs.harvard.edu/abs/1991dodg.conf..341T/abstract.

⁸⁶Note that in this article we use the nomenclature “Landau mode” regardless of whether the growth rate of the mode is positive, negative or zero.

⁸⁷J. Y. Lau and J. Binney, *Mon. Not. R. Astron. Soc.* **506**, 4007 (2021).

⁸⁸Again, the Maxwellian DF is not physically realizable in real astrophysical systems—see Sec. II B—but is well-defined in the periodic cube we are discussing here.

⁸⁹B. D. Fried and S. D. Conte, *The Plasma Dispersion Function* (Academic Press, New-York, 1961).

⁹⁰T. A. Zang, “The stability of a model galaxy,” Ph.D. thesis (MIT, 1976).

⁹¹M. S. Petersen, M. Roule, J.-B. Fouvry, C. Pichon, and K. Tep, *Mon. Not. R. Astron. Soc.* **530**, 4378 (2024).

⁹²J. Kormendy, J. Kennicutt, and C. Robert, *Annu. Rev. Astron. Astrophys.* **42**, 603 (2004).

⁹³J. A. Sellwood and K. L. Masters, *Annu. Rev. Astron. Astrophys.* **60**, 73 (2022).

⁹⁴S. De Rijcke, J.-B. Fouvry, and C. Pichon, *Mon. Not. R. Astron. Soc.* **484**, 3198 (2019).

⁹⁵J. A. Sellwood, *Rev. Mod. Phys.* **86**, 1 (2014).

⁹⁶T. O’Neil, *Phys. Fluids* **8**, 2255 (1965).

⁹⁷H. L. Berk and B. N. Breizman, *Phys. Fluids* **2**, 2235 (1990).

⁹⁸H. L. Berk and B. N. Breizman, “Enhancement of particle-wave energy exchange by resonance sweeping,” Technical Report No. DOE/ET/53088-722-Rev (UNT, 1996).

⁹⁹W. Drummond and D. Pines, “Non-linear stability of plasma oscillation,” Technical Report No. NSA-17-042107 (General Atomic Division, General Dynamics Corporation, 1961).

¹⁰⁰A. N. Kaufman, *Phys. Fluids* **15**, 1063 (1972).

¹⁰¹M. D. Weinberg, *Astrophys. J.* **421**, 481 (1994d).

¹⁰²D. C. Heggie, P. G. Breen, and A. L. Varri, *Mon. Not. R. Astron. Soc.* **492**, 6019 (2020).

¹⁰³J.-B. Fouvry and B. Bar-Or, *Mon. Not. R. Astron. Soc.* **481**, 4566 (2018).

¹⁰⁴P. Gibbon, *arXiv:2007.04783* (2020).

¹⁰⁵J. Magorrian, *Mon. Not. R. Astron. Soc.* **507**, 4840 (2021).

¹⁰⁶E. D’Onghia, M. Vogelsberger, and L. Hernquist, *Astrophys. J.* **766**, 34 (2013).

¹⁰⁷W. H. Julian and A. Toomre, *Astrophys. J.* **146**, 810 (1966).

¹⁰⁸J. Binney and C. Lacey, *Mon. Not. R. Astron. Soc.* **230**, 597 (1988).

¹⁰⁹P.-H. Chavanis, *Physica A* **391**, 3680 (2012).

¹¹⁰P.-H. Chavanis, *Astron. Astrophys.* **556**, A93 (2013).

¹¹¹R. Bälescu, *Phys. Fluids* **3**, 52 (1960).

¹¹²A. Lenard, *Ann. Phys.* **10**, 390 (1960).

¹¹³A. N. Kaufman and T. Nakayama, *Phys. Fluids* **13**, 956 (1970).

¹¹⁴A. N. Kaufman, *Phys. Fluids* **14**, 387 (1971).

¹¹⁵E. Dekker, *Phys. Rep.* **24**, 315 (1976).

¹¹⁶P.-H. Chavanis, *Eur. Phys. J. Plus* **139**, 51 (2024).

¹¹⁷J. F. Luciani and R. Pellat, *Astrophys. J.* **317**, 241 (1987).

¹¹⁸J. Heyvaerts, *Mon. Not. R. Astron. Soc.* **407**, 355 (2010).

¹¹⁹J. Heyvaerts, J.-B. Fouvry, P.-H. Chavanis, and C. Pichon, *Mon. Not. R. Astron. Soc.* **469**, 4193 (2017).

¹²⁰C. Hamilton, *Mon. Not. R. Astron. Soc.* **501**, 3371 (2021).

¹²¹N. Rostoker, *Phys. Fluids* **7**, 479 (1964).

¹²²N. Rostoker, *Phys. Fluids* **7**, 491 (1964).

¹²³F. P. C. Benetti and B. Marcos, *Phys. Rev. E* **95**, 022111 (2017).

¹²⁴B. Bar-Or and J.-B. Fouvry, *Astrophys. J. Lett.* **860**, L23 (2018).

¹²⁵J.-B. Fouvry, B. Bar-Or, and P.-H. Chavanis, *Phys. Rev. E* **99**, 032101 (2019).

¹²⁶M. Roule, J.-B. Fouvry, C. Pichon, and P.-H. Chavanis, *Phys. Rev. E* **106**, 044118 (2022).

¹²⁷M. D. Weinberg, *Astrophys. J.* **300**, 93 (1986).

¹²⁸J. A. Sellwood, *Astrophys. J.* **751**, 44 (2012).

¹²⁹J.-B. Fouvry, C. Pichon, J. Magorrian, and P.-H. Chavanis, *Astron. Astrophys.* **584**, A129 (2015).
¹³⁰A. Rogister and C. Oberman, *J. Plasma Phys.* **2**, 33 (1968).

¹³¹A. Rogister and C. Oberman, *J. Plasma Phys.* **3**, 119 (1969).

¹³²C. Oberman, “Advanced kinetic theory,” in *Physics of Hot Plasmas* (Springer, 1970) pp. 42–102.

¹³³T. Hatori, *Phys. Fluids* **11**, 1584 (1968).

¹³⁴T. Hatori, *Phys. Fluids* **12**, 1652 (1969).

¹³⁵C. Hamilton and T. Heinemann, *Mon. Not. R. Astron. Soc.* **525**, 4161 (2023).

¹³⁶C. Hamilton and T. Heinemann, arXiv:2011.14812 (2020).

¹³⁷M. Tagger, J.-F. Sygnet, E. Athanassoula, and R. Pellat, *Astrophys. J.* **318**, L43 (1987).

¹³⁸J.-F. Sygnet, M. Tagger, E. Athanassoula, and R. Pellat, *Mon. Not. R. Astron. Soc.* **232**, 733 (1988).

¹³⁹A. Toomre, *Astrophys. J.* **158**, 899 (1969).

¹⁴⁰A. Toomre, *Annu. Rev. Astron. Astrophys.* **15**, 437 (1977).

¹⁴¹J. W. K. Mark, *Astrophys. J.* **212**, 645 (1977).

¹⁴²C. Norman, *Mon. Not. R. Astron. Soc.* **182**, 457 (1978).

¹⁴³D. ter Haar, *Phys. Scr.* **20**, 291 (1979).

¹⁴⁴V. E. Zakharov, Zh. Eksp. Teor. Fiz. **35**, 908 (1972), see https://ui.adsabs.harvard.edu/abs/1972JETP...35.908Z/abstract.

¹⁴⁵B. Kocsis and S. Tremaine, *Mon. Not. R. Astron. Soc.* **448**, 3265 (2015).

¹⁴⁶J.-B. Fouvry, B. Bar-Or, and P.-H. Chavanis, *Astrophys. J.* **883**, 161 (2019).

¹⁴⁷J. I. Read, T. Goerdt, B. Moore, A. P. Pontzen, J. Stadel, and G. Lake, *Mon. Not. R. Astron. Soc.* **373**, 1451 (2006).

¹⁴⁸J. B. Taylor and B. McNamara, *Phys. Fluids* **14**, 1492 (1971).

¹⁴⁹R. Chiba, J. K. S. Friske, and R. Schönrich, *Mon. Not. R. Astron. Soc.* **500**, 4710 (2020).

¹⁵⁰R. Chiba and R. Schönrich, *Mon. Not. R. Astron. Soc.* **505**, 2412 (2021).

¹⁵¹J. A. Sellwood and R. G. Carlberg, *Mon. Not. R. Astron. Soc.* **517**, 2610 (2022).

¹⁵²J. A. Sellwood and F. D. Kahn, *Mon. Not. R. Astron. Soc.* **250**, 278 (1991).

¹⁵³J. A. Sellwood and R. G. Carlberg, *Mon. Not. R. Astron. Soc.* **489**, 116 (2019).

¹⁵⁴A. A. Tolman, E. A. Tolman, L. Arzamasskiy, and V. N. Duarte, *Astrophys. J.* **954**, 12 (2023).

¹⁵⁵A. Jenkins and J. Binney, *Mon. Not. R. Astron. Soc.* **245**, 305 (1990), see https://ui.adsabs.harvard.edu/abs/1990MNRAS.245.305J/abstract.

¹⁵⁶R. S. De Simone, X. Wu, and S. Tremaine, *Mon. Not. R. Astron. Soc.* **350**, 627 (2004).

¹⁵⁷J. Peñarrubia, *Mon. Not. R. Astron. Soc.* **490**, 1044 (2019).

¹⁵⁸A. A. El-Zant, J. Freundlich, F. Combes, and A. Halle, *Mon. Not. R. Astron. Soc.* **492**, 877 (2020).

¹⁵⁹B. Bar-Or, J.-B. Fouvry, and S. Tremaine, *Astrophys. J.* **871**, 28 (2019).

¹⁶⁰M. D. Weinberg and N. Katz, *Mon. Not. R. Astron. Soc.* **375**, 425 (2007).

¹⁶¹A. D. Ludlow, S. M. Fall, J. Schaye, and D. Obreschkow, *Mon. Not. R. Astron. Soc.* **508**, 5114 (2021).

¹⁶²Y.-S. Ting and H.-W. Rix, *Astrophys. J.* **878**, 21 (2019).

¹⁶³M. D. Weinberg, *Mon. Not. R. Astron. Soc.* **213**, 451 (1985).

¹⁶⁴D. Dootson and J. Magorrian, arXiv:2205.15725 (2022).

¹⁶⁵Note that *f* still depends parametrically on *N* and *f*<sub>*f*</sub>, though we suppress the explicit dependence here to keep the notation clean.

¹⁶⁶D. Lynden-Bell and A. J. Kalnajs, *Mon. Not. R. Astron. Soc.* **157**, 1–30 (1972).

¹⁶⁷L. D. Landau, Zh. Eksp. Teor. Fiz. **16**, 574 (1946), see https://inspirehep.net/literature/45141.

¹⁶⁸S. Ichimaru, *Phys. Rev.* **140**, B226 (1965).

¹⁶⁹R. W. Nelson and S. Tremaine, *Mon. Not. R. Astron. Soc.* **306**, 1 (1999).

¹⁷⁰M. D. Weinberg, arXiv:astro-ph/0404169 (2004).

¹⁷¹Note, however, that Hamilton¹⁵³ were working in three dimensions, since they were considering the interaction of a bar with a spherical halo. Here we are working in two dimensions, so our DF *f* is normalized differently.

¹⁷²C. Hamilton, *Mon. Not. R. Astron. Soc.* **528**, 5286 (2024).

¹⁷³R. L. Dewar, *Phys. Fluids* **16**, 431 (1973).

¹⁷⁴J. Binney, *New Astron. Rev.* **57**, 29 (2013).

¹⁷⁵G. Contopoulos, *Order and Chaos in Dynamical Astronomy* (Springer, 2002).

¹⁷⁶D. Merritt, *Dynamics and Evolution of Galactic Nuclei* (Princeton University Press, 2013).

¹⁷⁷S. Sridhar and J. Touma, *Mon. Not. R. Astron. Soc.* **303**, 483 (1999).

¹⁷⁸J. R. Cary and A. J. Brizard, *Rev. Mod. Phys.* **81**, 693 (2009).

¹⁷⁹I. Dodin and N. Fisch, *Phys. Rev. E* **77**, 036402 (2008).

¹⁸⁰P. J. Channell, *Phys. Fluids* **19**, 1541 (1976).

¹⁸¹H. Nyquist, *Bell Syst. Tech. J.* **11**, 126 (1932).

¹⁸²O. Penrose, *Phys. Fluids* **3**, 258 (1960).

¹⁸³P. Palmer, *Stability of Collisionless Stellar Systems* (Springer, 1994).

¹⁸⁴I. B. Bernstein, J. M. Greene, and M. D. Kruskal, *Phys. Rev.* **108**, 546 (1957).

¹⁸⁵A. M. Fridman and V. L. Polyachenko, *Physics of Gravitating Systems I: Equilibrium and Stability* (Springer, 2012).

¹⁸⁶J. Goodman, *Astrophys. J.* **329**, 612 (1988).

¹⁸⁷P. O. Vandervoort, *Mon. Not. R. Astron. Soc.* **339**, 537 (2003).

¹⁸⁸R. L. Dewar, *Phys. Fluids* **15**, 712 (1972).

¹⁸⁹E. Athanassoula, *Phys. Rep.* **114**, 319 (1984).

¹⁹⁰C. Dobbs and J. Baba, *Publ. Astron. Soc. Aust.* **31**, e035 (2014).

¹⁹¹J. C. Bird, S. Kazantzidis, and D. H. Weinberg, *Mon. Not. R. Astron. Soc.* **420**, 913 (2012).

¹⁹²N. I. Libeskind, R. van de Weygaert, M. Cautun *et al.*, *Mon. Not. R. Astron. Soc.* **473**, 1195 (2018).

¹⁹³P. F. Hopkins, P. Torrey, C.-A. Faucher-Giguère, E. Quataert, and N. Murray, *Mon. Not. R. Astron. Soc.* **458**, 816 (2016).

¹⁹⁴L. Wang, R. Spurzem, S. Aarseth, K. Nitadori, P. Berczik, M. B. N. Kouwenhoven, and T. Naab, *Mon. Not. R. Astron. Soc.* **450**, 4070 (2015).

¹⁹⁵A. A. Schekochihin, *Lecture Notes for the Oxford MMathPhys programme* (2022).

¹⁹⁶W. Dehnen, *Astrophys. J.* **118**, 1190 (1999).