5.3. The collisionless Boltzmann equation¶
To go beyond the virial theorem, we now study the evolution of gravitational systems consisting of large numbers of particles in more depth. The fundamental equation describing collisionless systems is the collisionless Boltzmann equation. This equation describes the behavior of the distribution function \(f(\vec{x},\vec{v},t)\) of bodies moving under the influence of gravity. The distribution function gives the number of bodies in a small volume of phase space at a given time \(t\): \(\mathrm{d}N(\vec{x},\vec{v},t) \propto f(\vec{x},\vec{v},t)\,\mathrm{d}\vec{x}\,\mathrm{d}\vec{v}\), with the proportionality set by how one chooses to normalize \(f\). We will again use \(\vec{w} \equiv (\vec{x},\vec{v})\) to denote phase-space coordinates. The “collisionless” assumption implies that this “one-body” distribution function fully describes the system, that is, that there are no correlations between the number of bodies at different phase-space positions. Said another way, the assumption is that the presence of a body at one position does not affect the probability of seeing a body at another position. Thus, the \(N\)-body distribution function \(f^{(N)}(\vec{w}_1,\vec{w}_2,\ldots,\vec{w}_N,t)\) is simply given by \begin{equation} f^{(N)}(\vec{w}_1,\vec{w}_2,\ldots,\vec{w}_N,t) = \prod_{i=1}^N\,f(\vec{w}_i,t)\,, \end{equation} for all \(N\). We further normalize \(f(\vec{w},t)\) such that \(\int \mathrm{d}\vec{w}\,f(\vec{w},t) = 1\) for all times \(t\). We can then consider what happens to the total number of stars \(N(\mathcal{V},t) = \int_{\mathcal{V}} \mathrm{d}\vec{w}\, f(\vec{w},t)\) in a volume \(\mathcal{V}\) of phase space with a surface \(S\). As time passes, the number of bodies in \(\mathcal{V}\) will change as \(\mathrm{d} N(\mathcal{V},t) / \mathrm{d} t = \int_{\mathcal{V}} \mathrm{d}\vec{w}\,\partial f(\vec{w},t) / \partial t\). Assuming that no bodies are created or removed inside the volume, this change happens because bodies leave and enter \(\mathcal{V}\) through the boundary \(S\) at a rate of \(f(\vec{w},t)\,\dot{\vec{w}}_\perp\) (the subscript \(\perp\) denotes that the flow is perpendicular to the surface, with a positive rate defined as pointing outward). The total flow across the surface is the integral over the entire surface of this flow: \(\oint_S \mathrm{d} S\, f(\vec{w},t)\,\dot{\vec{w}}_\perp\), and because the number of stars must be conserved, this must be equal to the change in the number of stars in \(\mathcal{V}\): \begin{equation} \int_{\mathcal{V}} \mathrm{d}\vec{w}\,\frac{\partial f(\vec{w},t)}{\partial t} = -\oint_S \mathrm{d} S\, f(\vec{w},t)\,\dot{\vec{w}}_\perp\,. \end{equation} Using the divergence theorem from Equation (B.5), we can write the integral on the right-hand side as a volume integral as well \begin{equation} \int_{\mathcal{V}} \mathrm{d}\vec{w}\,\frac{\partial f(\vec{w},t)}{\partial t} = -\int_{\mathcal{V}} \mathrm{d} \vec{w}\,\nabla \cdot [f(\vec{w},t)\,\dot{\vec{w}}]\,. \end{equation} Because this equation must hold for any volume \(\mathcal{V}\), it must be the case that \begin{equation} \frac{\partial f(\vec{w},t)}{\partial t} +\nabla \cdot [f(\vec{w},t)\,\dot{\vec{w}}]=0\,, \end{equation} which is the continuity equation for the phase-space density. Now let’s work out the gradient in terms of the Cartesian position and momentum coordinates \((\vec{x},\vec{v})\) \begin{equation}\label{eq-coll-boltzmann-penultstep} \frac{\partial f(\vec{w},t)}{\partial t} +\frac{\partial}{\partial \vec{x}} \cdot [f(\vec{x},\vec{v},t)\,\dot{\vec{x}}] +\frac{\partial}{\partial \vec{v}} \cdot [f(\vec{x},\vec{v},t)\,\dot{\vec{v}}]=0\,. \end{equation} Because \(\partial \dot{\vec{x}} / \partial \vec{x} = \partial \dot{\vec{v}} / \partial \vec{v} = 0\) this becomes \begin{equation}\label{eq-collisionless-boltzmann-cartesian} \frac{\partial f(\vec{x},\vec{v},t)}{\partial t} +\dot{\vec{x}}\,\frac{\partial f(\vec{x},\vec{v},t)}{\partial \vec{x}} +\dot{\vec{v}}\,\frac{\partial f(\vec{x},\vec{v},t)}{\partial \vec{v}}=0\,, \end{equation} which is the collisionless Boltzmann equation. This equation in fact holds for any set of canonical coordinates \((\vec{q},\vec{p})\). The derivation up to Equation (5.27) is the same, because \(f(\vec{q},\vec{p},t) = f(\vec{x},\vec{v},t)\) (which in turn follows from the fact that \(\mathrm{d}\vec{q}\,\mathrm{d}\vec{p} = \mathrm{d}\vec{x}\,\mathrm{d}\vec{v}\) for any set of canonical coordinates as we showed in Section 3.4.2) . The last step leads to \begin{equation}\label{eq-collisionless-boltzmann} \frac{\partial f(\vec{q},\vec{p},t)}{\partial t} +\dot{\vec{q}}\,\frac{\partial f(\vec{q},\vec{p},t)}{\partial \vec{q}} +\dot{\vec{p}}\,\frac{\partial f(\vec{q},\vec{p},t)}{\partial \vec{p}}=0\,, \end{equation} in this case because \(\partial \dot{\vec{q}} / \partial \vec{q} = -\partial \dot{\vec{p}} / \partial \vec{p}\) from Hamilton’s equations. Equation (5.29) makes it much easier to write down the collisionless Boltzmann equation in coordinate systems other than Cartesian coordinates. In Cartesian coordinates, we can re-write Equation (5.28) as \begin{equation}\label{eq-collisionless-boltzmann-cartesian-wpot} \frac{\partial f(\vec{x},\vec{v},t)}{\partial t} +\dot{\vec{x}}\,\frac{\partial f(\vec{x},\vec{v},t)}{\partial \vec{x}} -\frac{\partial \Phi}{\partial \vec{x}}\,\frac{\partial f(\vec{x},\vec{v},t)}{\partial \vec{v}}=0\,, \end{equation} by using Newton’s second law. For general canonical coordinates \((\vec{q},\vec{p})\), we have a similar equation using Hamilton’s equations \begin{equation}\label{eq-collisionless-boltzmann-general} \frac{\partial f(\vec{q},\vec{p},t)}{\partial t} +\frac{\partial H}{\partial \vec{p}}\,\frac{\partial f(\vec{q},\vec{p},t)}{\partial \vec{q}} -\frac{\partial H}{\partial \vec{q}}\,\frac{\partial f(\vec{q},\vec{p},t)}{\partial \vec{p}}=0\,, \end{equation}
The collisionless Boltzmann equation describes the evolution of any collisionless system, not just those that are in equilibrium. It is the fundamental equation of collisionless dynamics and the equation that is solved by \(N\)-body simulations. If the system under consideration is self-gravitating, the density of stars sources the potential, and the Poisson equation becomes an additional evolution equation that relates the gravitational potential and the distribution function \begin{align} \nabla^2 \Phi(\vec{x},t) & = 4\pi\,G\,\rho(\vec{x},t) = 4\pi\,G\,M\,\int \mathrm{d} \vec{v}\,f(\vec{x},\vec{v},t)\,,\label{eq-spherequil-selfconsist} \end{align} where \(M\) is the total mass of the system (because of our convention that \(\int \mathrm{d}\vec{x}\mathrm{d}\vec{v}\,f(\vec{x},\vec{v},t) = 1\), we have that \(\rho(\vec{x},t) = M\,\int \mathrm{d} \vec{v}\,f(\vec{x},\vec{v},t)\) such that \(\int \mathrm{d}\vec{x}\, \rho(\vec{x},t) = M\,\int \mathrm{d}\vec{x}\,\mathrm{d} \vec{v}\,f(\vec{x},\vec{v},t) = M\)). Such a system is a self-consistent system. But the collisionless Boltzmann equation does not only apply to self-consistent systems, it more generally describes the evolution of a subset of bodies orbiting in a general smooth mass distribution that they do not necessarily contribute much mass to (for example, the globular cluster system in the Milky Way; see above).
For a system in equilibrium, the number of stars in a given region of phase-space does not change, so \(\partial f(\vec{x},\vec{v},t) / \partial t = 0\). For such systems, the collisionless Boltzmann equation reduces to \begin{equation}\label{eq-collisionless-boltzmann-equil} \dot{\vec{q}}\,\frac{\partial f(\vec{q},\vec{p})}{\partial \vec{q}} +\dot{\vec{p}}\,\frac{\partial f(\vec{q},\vec{p})}{\partial \vec{p}}=0\,. \end{equation} This is the equilibrium collisionless Boltzmann equation. In the next sections we study solutions and properties of this equilibrium equation.
Before continuing, we will consider one more corollary of the collisionless Boltzmann equation. The collisionless Boltzmann equation describes how the phase-space density \(f(\vec{x},\vec{v})\) behaves in general. But we can also consider how \(f(\vec{x},\vec{v})\) evolves along the orbit of a body. The time evolution of \(f(\vec{x},\vec{v})\) along the orbit is given by \begin{align} \frac{\mathrm{d} f(\vec{x},\vec{v})}{\mathrm{d} t} & = \frac{\partial f(\vec{x},\vec{v})}{\partial t} + \dot{\vec{x}}\,\frac{\partial f(\vec{x},\vec{v})}{\partial \vec{x}} + \dot{\vec{v}}\,\frac{\partial f(\vec{x},\vec{v})}{\partial \vec{v}}= 0\,,\label{eq-liouville} \end{align} because the right-hand side is simply the expression in the collisionless Boltzmann equation (5.29). Thus, we find Liouville’s theorem: the phase-space density is conserved along orbital trajectories. This is a very useful property of Hamiltonian systems.
Let’s look at an example of Liouville’s theorem. In a simple one-dimensional potential (one-dimensional such that we can directly visualize the two-dimensional phase space), we start with a set of \(n\) initial points \((x_i,v_i)\) along the edge of a small circular area. We then integrate all of these orbits for some time and plot the evolution of the phase-space area spanned by the \((x_i[t],v_i[t])\) at different times \(t\) in Figure 5.3.
[14]:
from galpy.orbit import Orbit
from galpy.potential import LogarithmicHaloPotential
from galpy.util import plot as galpy_plot
from matplotlib import pyplot
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
# Simple one-dimensional potential
lp= LogarithmicHaloPotential(normalize=1.).toVertical(1.)
# Trace out initially circular phase-space area
npoints= 101
theta= numpy.linspace(0,2.*numpy.pi,npoints,endpoint=True)
phase_space_radius= 0.15
xs= phase_space_radius*numpy.cos(theta)+1.
vs= phase_space_radius*numpy.sin(theta)+1
# Integrate all points along the edge of the area
os= [Orbit([x,v]) for x,v in zip(xs,vs)]
nts= 11
ts= numpy.linspace(0.,10.5,nts)
[o.integrate(ts,lp) for o in os]
# Create matplotlib polygons to represent area at each time
patches= []
for ii,t in enumerate(ts):
polygon= Polygon(numpy.array([[o.x(t),o.vx(t)] for o in os]),
closed=True)
patches.append(polygon)
p= PatchCollection(patches,alpha=1.)
fig, ax= pyplot.subplots(figsize=(10,5))
ax.set_aspect('equal')
ax.add_collection(p)
xlim(-3.,3)
ylim(-1.5,1.5)
xlabel(r'$x$')
ylabel(r'$v$')
pyplot.text(2.1,1.2,r'$t=0$',size=18.)
pyplot.show();
Figure 5.3: Liouville’s theorem.
The initially circular phase-space area quickly gets stretched out in one direction and squashed in the other and then gets deformed as time progresses. But it is clear that the total area is conserved: the area becomes progressively thinner, but also gets longer in such a way that the total area does not change.
.
.
.