17.3. Virialization¶
The formation of dark-matter halos in a quasi-equilibrium state—“virialized”, because they satisfy the virial theorem—from dynamically-cold initial conditions is a complex process that needs to be studied numerically. Of particular interest is the question of how relaxed, quasi-equilibrium systems form through the effects of self-gravity and what the properties of the relaxed system are. In Chapter 5.1, we saw that the gravitational dynamics of galactic systems is collisionless—collisions between individual bodies play no role in its evolution—and that the distribution function therefore satisfies the collisionless Boltzmann equation. In Chapter 5.3 we demonstrated that this implies Liouville’s theorem: the phase-space density is conserved along orbital trajectories. “Relaxation” in the thermodynamical sense therefore cannot happen: The six-dimensional phase-space distribution function does not truly change (its total time derivative is zero) and there can therefore be no loss of information or increase of entropy associated with its evolution.
The conservation of the phase-space distribution function holds only on a microscopic level, however, not on the macroscopic, observable level. The evolution of a gravitational system causes the phase-space distribution to wrap up in complicated ways, and when averaged—or “coarse-grained”—over a macroscopic volume of phase-space, the distribution function will typically appear to have decreased. On the macroscopic level, systems therefore will appear to approach a relaxed state, even though they are no more relaxed than they were initially at the microscopic level. This type of relaxation that is due to the coarse-grained distribution function evolving while the microscopic one remains constant is called collisionless relaxation.
In the following subsections, we discuss the processes through which dark-matter halos and other self-gravitating systems virialize and relax. First, in Section 17.3.1, we use the virial theorem to understand the properties of the end state of the virialization process and we derive the standard result that the average density within collapsed dark-matter halos is \(\approx 200\) times the mean background matter density of the Universe. In the next sections, we discuss how dark-matter halos attain equilibrium through the growth and damping of instabilities and various collisionless relaxation mechanisms. As part of this discussion, we will touch upon a number of standard results about the stability of self-gravitating fluids and collisionless systems that are useful, for example, in studies of star formation and the stability of globular clusters.
While we focus here on the initial path towards a relaxed, quasi-equilibrium dark-matter halo, it is important to keep in mind that dark-matter halos continue to evolve after their initial collapse through the accretion of additional dark matter through smooth accretion and through minor and major mergers (see Chapter 19). The relaxation processes that we discuss also lead to the quick re-establishment of a relaxed state after such accretion events.
17.3.1. Overall properties of virialized dark-matter halos¶
The spherical-collapse picture from Section 17.2 above has serious limitations. In the spherical-collapse approximation, the radial-only collapse without shell crossing would continue until all shells reach \(r=0\) at \(t=t_c\). In reality this does not happen, because deviations from spherical symmetry will lead to shell-crossing resulting in non-linear dynamical interactions. Exactly what happens is a matter for \(N\)-body simulations, but it is plausible to assume that the end state is a collapsed mass distribution that is close to virial equilibrium. We can then compute the actual physical overdensity of the collapsed, virialized object (which we refer to as a “halo”) with respect to the background density. We work in a matter-only Universe for simplicity and because this is approximately the case for our Universe during the main epoch of halo formation.
We again consider the top-hat overdensity from the spherical-collapse calculation above. In the spherical-collapse picture, all shells expand until they reach their maximum radius \(r_\mathrm{max}\) at the same time \(t=t_c/2\) (or \(\eta = \pi\) in the parameterization of Equations 17.81); this time is known as the turn-around time. We assume that this dynamics still holds and that virialization happens at \(t > t_c/2\). At turn-around, all of the energy of the collapsing object is purely in the form of potential energy, because all shells turn around at the same time and their velocities are thus zero. As discussed in Chapter 14.1.2, the total potential energy of a density distribution \(\rho(\vec{x})\) is \begin{equation} W = \frac{1}{2}\,\int\mathrm{d}\vec{x}\,\rho(\vec{x})\,\Phi(\vec{x})\,. \end{equation} For a homogeneous density sphere of radius \(r_\mathrm{max}\), the potential energy is given by (see Chapter 14.1.3): \begin{align} W_{\mathrm{max}} & = -\frac{3}{5}\,\frac{GM^2}{r_\mathrm{max}}\,. \end{align}
After the collapsed, virialized halo has formed, the virial theorem holds which in the scalar form derived in Chapter 14.1.2 can be written as that the total energy of the system \(E = W_v/2\), where \(W_v\) is the potential energy after virialization. The total energy of the system is conserved during the collapse and therefore we have that \begin{equation} W_v = 2\,W_\mathrm{max}\,. \end{equation} Approximating the mass distribution of the collapsed object as a homogeneous density sphere with a radius \(r_\mathrm{vir}\)—the virial radius—this equation implies that \begin{equation} r_\mathrm{vir} = \frac{r_\mathrm{max}}{2}\,. \end{equation} A virialized halo should therefore have a radius that is half of the turn-around radius.
Next we compute the physical overdensity \(\Delta_v\) of the collapsed halo with respect to the cosmological critical density. From Equations (17.81) for spherical collapse, we see that in the spherical collapse picture, \(r = r_\mathrm{max}/2\) at \(\eta = 3\pi/2\) (\(\eta > \pi\), because we are interested in the collapsing part of the solution). This corresponds to \(t = (3/4+1/[2\pi])\,t_c\approx0.9\,t_c\). In reality it will take somewhat longer to attain virial equilibrium and we will instead assume that virialization is achieved at \(t = t_c\). We then compute the virial overdensity as \begin{equation} \Delta_v = \frac{\rho(t_v)}{\rho_c(t_v)} = \frac{M_\mathrm{vir}/[4\pi r_\mathrm{vir}^3/3]}{\rho_c(t_v)}\,, \end{equation} where \(\rho_c(t_v)\) is the critical density from Equation (17.3) at the time of virialization. Using the fact that the total mass is conserved during the collapse and the mass \(M_\mathrm{vir}\) is the eventual virial mass contained within the virial radius (we assume no mass escapes the system, even after spherical-collapse fails to be a good approximation of the virialization process), we can compute this by expressing everything in terms of the initial quantities \(r_i\) and \(t_i\) from the previous section \begin{align} \Delta_v & = \frac{3}{4\pi\,r_\mathrm{vir}^3}\,\frac{4\pi\,r_i^3}{3}\,\bar{\rho}_m(t_i)\,(1+\delta_i)\,\frac{1}{\rho_c(t_v)} = \left(\frac{r_i}{r_\mathrm{vir}}\right)^3\,(1+\delta_i)\,\frac{\bar{\rho}_m(t_i)}{\rho_c(t_v)}\,. \end{align} Now we use that \(r_\mathrm{vir} = r_\mathrm{max}/2\) and Equation (17.88) to obtain \begin{equation} \Delta_v = 8\,\left(1-\left[1-\frac{5}{3}\,\delta_i\right]\,\Omega_{i,m}^{-1}\right)^3\,(1+\delta_i)\,\frac{\bar{\rho}_m(t_i)}{\rho_c(t_v)}\,. \end{equation}
To go further, we need to know the behavior of the background density with time. For a general matter-only Universe with \(\Omega_m \neq 1\), this is complicated, but in the Einstein–de-Sitter Universe model that we are using we have that \(\rho_c(t_v) = \bar{\rho}_m(t_v)\) and from Equation (17.9) we find \({\bar{\rho}_m(t_i)/\rho_c(t_v)} = \left({t_v/t_i}\right)^2\). Plugging this into the equation for \(\Delta_v\), remembering that \(t_v = t_c\), and using Equation (17.89), this then simplifies to \begin{equation}\label{eq-dmhaloform-deltavir-eds} \Delta_v = 18\,\pi^2\,. \end{equation}
Thus, collapsed, virialized halos in a matter-only, \(\Omega_m = 1\) Universe have an overdensity with respect to the background of \(\Delta_v = 18\pi^2 \approx 178\). This value is therefore often used to define the boundary of a dark matter halo: the virial radius \(r_\mathrm{vir}\) is such that the mean overdensity within the virial radius is \(\approx180\) times the critical density. Comparing this \(\approx180\) number to the linear overdensity \(\delta_c \approx 1.686\) at the collapse time derived above, we get a sense of how fast overdensities grow once they enter the non-linear regime: the actual physical overdensity of the halo that forms is more than a hundred times larger than what the overdensity would have been under linear evolution alone.
The general result for an open matter-only Universe with \(\Omega_m \leq 1\) is \begin{equation}\label{eq-dmhaloform-deltavir-open-exact} \Delta_v = 8\pi^2\,\left[{1\over x}-{\left(1-x\right)\,\mathrm{sinh}^{-1}\left(\sqrt{{x\over 1-x}}\right) \over x^{3/2}}\right]^{-2}\,, \end{equation} where \(x = 1-\Omega_m(t_c)\) (Lacey & Cole 1993). This is well approximated by the following form (Bryan & Norman 1998) \begin{equation}\label{eq-dmhaloform-deltavir-open-approx} \Delta_v = 18\,\pi^2-60x-32x^2\,. \end{equation} For a flat Universe consisting of matter and dark energy (\(\Omega_m + \Omega_\Lambda = 1\))—which is the Universe that we live in—one finds that (Bryan & Norman 1998) \begin{equation}\label{eq-dmhaloform-deltavir-flat-de} \Delta_v = 18\,\pi^2-82x-39x^2\,, \end{equation} with again \(x = 1-\Omega_m(t_c)\). Note that if we consider the virial overdensity with respect to the mean matter density instead of the critical density, there is an additional \(\Omega_m^{-1}(t_c)\) factor in the expressions for \(\Delta_v\).
For our own Universe, \(\Omega_m \approx 0.31\) and \(\Omega_\Lambda \approx 0.69\) at present (Equation 17.5), so halos that form today have \begin{equation} \Delta_v \approx 100\,,\quad (\mathrm{today})\,, \end{equation} with respect to the critical density and an overdensity of \(\approx 330\) with respect to the mean matter density. Many galaxies that we see today formed when the Universe was much younger (say \(z \gtrsim 2\)) and then \(\Omega_m \approx 1\), such that the standard value applied at that time \begin{equation} \Delta_v \approx 178\,,\quad (z\gtrsim2)\,. \end{equation}
It is important to remember that these overdensities are with respect to the critical density at the time of collapse, not with respect to the density today. In particular, halos that form earlier when the Universe was denser have higher densities than those that form later, because even though their \(\Delta_v\) is similar, the mean density was higher.
Finally, we note that the virial mass and the virial radius that emerge from the simple spherical-collapse+virialization picture presented here do not appear to be as fundamental as they appear, because the ongoing accretion onto halos (see Chapter 19.2) causes the mass and radial boundary to change over time. An alternative halo boundary that has received much attention in recent years is the splashback radius (Diemer & Kravtsov 2014; Adhikari et al. 2014; More et al. 2015), the first apocenter radius of halo particles orbiting within the halo. In the self-similar collapse solutions with radial orbits from Fillmore & Goldreich (1984) and Bertschinger (1985), this radius cleanly distinguishes infalling particles from particles orbiting in the halo and it corresponds to a steep density drop, which is also observed in simulations (Diemer & Kravtsov 2014). The location of the splashback radius depends on the accretion rate and is \(\approx r_\mathrm{vir}\) for rapidly-accreting halos and \(\approx 1.5\,r_\mathrm{vir}\) for slowly-accreting halos, when we define the virial radius using \(\Delta_v = 200\) with respect to the mean matter density (More et al. 2015).
17.3.2. Phase mixing¶
Phase mixing is perhaps the most important collisionless relaxation process that occurs in galaxies, because it is involved in creating the initial virialized dark-matter halo, but is also instrumental every time a a dark matter halo or a galactic component is perturbed out of its quasi-equilibrium state (e.g., by a minor merger) in driving the system back to quasi-equilibrium.
Phase mixing refers to the process of driving the “phases” of a set of bodies to a uniform distribution, where the “phases” refers to the time-dependent component of an orbit. For example, in the epicycle approximation of orbits in galactic disks that we discussed in Chapter 9.3.1, an orbit consists of three sinusoidal oscillations and each oscillation has an associated amplitude that is conserved and an associated phase angle that increases linearly in time. More generally, all regular orbits in a galaxy can be described in action-angle coordinates (see Chapter 3.4.3), where the actions are three integrals of motion and the angles are three phase angles that increase linearly in time according to three frequencies. Phase mixing occurs when the phases of a set of orbits that are initially clustered around a given phase become uniform.
Phase mixing occurs because the frequencies of nearby points in phase-space generally differ and after long periods of time, the points therefore find themselves at very different phases, even if they started out close together. On a microscopic level, the phase-space distribution remains as featured as it was in the initial state—no information is lost during the evolution of regular orbits—but on a macroscopic level the distribution will appear to become more uniform. As the phases become uniform, the system will appear to become time-independent, because time dependence for regular orbits can only come from the evolution of the phase distribution, which ceases once it is uniform. Phase mixing therefore drives galactic systems to a quasi-equilibrium state.
As an example, we consider the evolution of a disk of initial conditions in a one dimensional system clustered around an initial \((x,v)\). First we compute the evolution of the phase-space volume of this disk and display it in Figure 17.5.
[5]:
from galpy.orbit import Orbit
from galpy.potential import LogarithmicHaloPotential
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 (in x = 2v)
npoints= 3001
theta= numpy.linspace(0,2.*numpy.pi,npoints,endpoint=True)
phase_space_radius= 0.25
xs= phase_space_radius*2.*numpy.cos(theta)+1.
vs= phase_space_radius*1.*numpy.sin(theta)+1
# Integrate all points along the edge of the area, making use
# of galpy's capability to handle multiple orbits at once
os= Orbit([[x,v] for x,v in zip(xs,vs)])
nts= 1001
ts= numpy.linspace(0.,220.5,nts)
os.integrate(ts,lp)
ts= [[0,9],[70],[220]]
# Plot
fig, axes= subplots(1,3,figsize=(10,3))
axes= axes.flatten()
for ii,(ax,tii) in enumerate(zip(axes,ts)):
for t in tii:
# Create matplotlib polygons to represent area at each time
polygon= Polygon(numpy.array([os.x(t),os.vx(t)]).T,
closed=True)
p= PatchCollection([polygon],alpha=1.)
sca(ax)
ax.add_collection(p)
xlim(-3.5,3.5)
ylim(-1.75,1.75)
xlabel(r'$x$')
ylabel(r'$v$')
ax.tick_params(labelbottom=False,labelleft=False)
galpy_plot.text(rf'$t={",".join(["{:d}".format(t) for t in tii])}$',top_left=True,size=18.)
tight_layout(w_pad=0.2);
Figure 17.5: The evolution of the fine-grained phase-space distribution.
Because the frequencies of different orbits present in the initial disk differ, the disk quickly starts to spread. Eventually the disk turns into a thin ribbon that wraps around the center many times. Because the microscopic phase-space volume is conserved, the total area of the ribbon is the same as the area of the initial disk, but now spread over a much wider area of phase-space. If we coarse-grain the distribution function, the area of the region covered by the ribbon would be larger and the distribution function would therefore be lower. To illustrate that phase-mixing happens in this system at the coarse-grained level, we instead sample points uniformly within the initial disk and integrate their orbits. In the top row of Figure 17.6, we plot their phase-space positions at the same times as in Figure 17.5, showing individual points with large symbols to perform a visual coarse-graining.
[7]:
# Simple one-dimensional potential
lp= LogarithmicHaloPotential(normalize=1.).toVertical(1.)
# sample points within initially circular phase-space area (in x = 2v)
npoints= 3001
from numpy.random import default_rng
rng= default_rng(1)
theta= rng.random(size=npoints)*2.*numpy.pi
phase_space_radius= 0.25*numpy.sqrt(rng.random(size=npoints))
xs= phase_space_radius*2.*numpy.cos(theta)+1.
vs= phase_space_radius*1.*numpy.sin(theta)+1
# Integrate all points, making use
# of galpy's capability to handle multiple orbits at once
os= Orbit([[x,v] for x,v in zip(xs,vs)])
nts= 3001
ts= numpy.linspace(0.,520.5,nts)
os.integrate(ts,lp)
ts= [[0,9],[70],[220],[320],[420],[520]]
# Contour the allowed energies
Es= vs**2./2.+lp(xs)
Emin, Emax= numpy.amin(Es), numpy.amax(Es)
Exs, Evs= numpy.meshgrid(numpy.linspace(-3.5,3.5,51),
numpy.linspace(-1.75,1.75,52),indexing='ij')
Egrid= Evs**2./2.+lp(Exs)
# Plot
fig, axes= pyplot.subplots(2,3,figsize=(10,5.8),dpi=200)
axes= axes.flatten()
for ii,(ax,tii) in enumerate(zip(axes,ts)):
sca(ax)
for jj,t in enumerate(tii):
line = ax.plot(os.x(t),os.vx(t),'o',
color=None if jj == 0 else line[0].get_color(),
alpha=0.3,ms=8,mec='none',rasterized=True)
# Contour Emin and Emax, allowed energies
ax.contour(Exs,Evs,Egrid,levels=[Emin,Emax],colors='0.75',lw=1.5)
xlim(-3.5,3.5)
ylim(-1.75,1.75)
xlabel(r'$x$')
ylabel(r'$v$')
ax.tick_params(labelbottom=False,labelleft=False)
galpy_plot.text(rf'$t={",".join(["{:d}".format(t) for t in tii])}$',top_left=True,size=18.)
tight_layout(w_pad=0.2);
Figure 17.6: The evolution of the coarse-grained phase-space distribution.
It is clear that at the final time of the top row, the distribution of the points appears to uniformly occupy the part of phase-space allowed by the initial energies (the range of the initial energies is shown by the gray contours). The phase-space distribution is therefore phase-mixed. If we integrate further, we see in the bottom row of Figure 17.6 that the coarse-grained phase-space distribution essentially stops evolving. Phase mixing occurs in all non-uniform mass density distributions. The time scale of phase-mixing depends on the initial frequency distribution as \(\tau \approx 2\pi/\Delta \Omega\), where \(\Delta \Omega\) is the typical spread in the initial frequencies. If the initial distribution has \(\mathcal{O}(1)\) frequency differences, phase-mixing occurs on a dynamical time, because \(\tau \approx 2\pi/\Delta \Omega \approx 2\pi/\Omega = t_\mathrm{dyn}\). If the frequency distribution is narrow—or “colder”—the phase-mixing time scale is longer and can be many Gyr for very cold systems (for example, globular cluster streams phase-mixing into the halo of the Milky Way have \(\Omega / \Delta \Omega \approx 100\), so \(\tau \approx 100\,t_\mathrm{dyn}\)).
17.3.3. The (in)stability of self-gravitating systems¶
Phase mixing is a relatively benign process that happens predictably, at a steady rate, and within both evolving and static mass distributions. Because phase mixing only changes the phases of the orbits, but not their integrals of motion like the energy, it leads to only small changes to the global structure of a system. In static systems, the three integrals of the motion—the actions in action-angle coordinates—are conserved for regular orbits and the only way to mix the integrals of the motion themselves is therefore to have a time-dependent gravitational potential. Time-dependence occurs naturally during the formation and evolution of dark-matter halos, both because the initial state of the dark-matter halo is highly non-phase-mixed and, thus, the phase-mixing evolution from the previous section leads to a strongly time-dependent density and gravitational potential, and because perturbations to the halo cause time-dependent fluctuations in the halo’s gravitational field.
A perturbation to a system for which self-gravity is unimportant—e.g., a low-mass disk galaxy embedded within a massive dark-matter halo—would lead to a perturbed state that relaxes through the phase-mixing process described in the previous section. But if self-gravity is important—as it is for the dark-matter halo as a whole—the perturbation gives rise to collective effects that can amplify or dampen the perturbation. The collective response of a self-gravitating system to different types of perturbations is an important topic in the study of fluids (e.g., gas disks) and collisionless systems (e.g., stellar clusters or dark-matter halos), because it allows us to determine which fluid or collisionless systems are stable and, thus, likely to be observed in nature. Yet despite its importance, useful mathematical results are few and far between, because a numerical approach is necessary to deal with the complexity of realistic systems and of the \(N\)-body dynamics. We will discuss this topic with more mathematical rigor in Chapter 20.1 to determine whether and when disk galaxies are dynamically stable. Here we focus on a qualitative understanding of the evolution of perturbed gravitational systems.
In response to a perturbation, a self-gravitating collisionless or fluid system can either start oscillating with fixed amplitude or the perturbation can grow or decay. While stable oscillations are common in fluids (e.g., stars, which have many oscillation modes that give rise to the field of asteroseismology), oscillations are rare for collisionless systems and any oscillatory modes are generally damped. Typically, there is a spatial scale at which perturbations become unstable, in a Fourier decomposition of the perturbation along the lines of Section 17.1.4.1, which for astrophysical systems is the Jeans length from Equation (17.55). We can re-use our mathematical framework describing the evolution of baryonic and dark-matter perturbations from Section 17.1 to qualitatively understand the possible behaviors.
Considering the evolution of perturbations in the cosmological baryonic fluid, we derived Equation (17.53) that gives the evolution of perturbations on Fourier scale \(\vec{k}\) and we found that the Jeans length defined in Equation (17.55) corresponding to the Jeans wavenumber \(k_J\) (Equation 17.54) separates two regimes: on large spatial scales \(k < k_J\), the fluid pressure is unimportant and the perturbation evolves just as it does in the pressureless, collisionless regime discussed in Section 17.1.2, while on small spatial scales \(k > k_J\), the effect of pressure becomes important. On large-scales, the evolution equation is given by Equation (17.25) and, as we saw in Section 17.1.4.2, generally this equation has two solutions: a growing mode and a decaying mode. The existence of the growing mode implies that perturbations grow on large spatial scales \(k < k_J\) and perturbations are thus unstable. This is referred to as gravitational instability or Jeans instability. That this is the case even when considering systems that are not themselves expanding like the Universe as a whole (e.g., a formed dark-matter halo), is clear from considering Equation (17.25) with the expansion of the Universe turned off, which corresponds to setting \(H = \dot{a}/a = 0\) and \(a=1\). Then we are left with \begin{equation}\label{eq-ddotdelta-dm-noexpansion} \frac{\partial^2 \delta_{\vec{k}}}{\partial t^2} = 4\pi\,G\,\overline{\rho}\,\delta_{\vec{k}}\,, \end{equation} which has as solutions the growing mode \(\propto e^{\omega t}\) and the decaying mode \(\propto e^{-\omega t}\) where \(\omega = \sqrt{4\pi G\,\overline{\rho}}\). The expansion of the Universe leads to a damping of this mode, known as the Hubble drag, but the drag only fully dampens the growing mode if the expansion of the Universe has the fast exponential form of a dark-energy dominated Universe.
On small scales \(k > k_J\), the existence of the \(k^2/k_J^2\) term on the right-hand side of Equation (17.53) leads to different behavior. Again considering a fluid in a non-expanding Universe, e.g., a gas cloud in a galaxy, Equation (17.53) becomes \begin{equation}\label{eq-ddotdelta-with-pressure-ink-baryons-noexpansion} \frac{\partial^2 \delta_{\vec{k}}}{\partial t^2}= 4\pi\,G\,\overline{\rho}\,\left(1-\frac{k^2}{k_J^2}\right)\,\delta_{\vec{k}}\,. \end{equation} Because \(k > k_J\), the right-hand side is negative and the solutions are therefore oscillating modes \(\propto e^{\pm i\omega t}\) where \(\omega = \sqrt{4\pi\,G\,\overline{\rho}\,\left(k^2/k_J^2-1\right)}\). The effect of the Hubble drag in an expanding Universe is to dampen this oscillation. Our derivation has not taken the self-gravity of the fluid, that is, its response to the perturbation, into account and is therefore an incomplete description of the behavior of fluids under perturbations. But self-gravity does not qualitatively change the behavior. Thus, perturbations to fluid systems are stable or damped on small spatial scales \(k > k_J\). The behavior of fluid perturbations that we have described in a static, non-expanding background are highly relevant for the evolution of gaseous astrophysical systems, playing for example a large role in star formation through the gravitational collapse of gas in galaxies. The Jeans wavenumber and length in this case are given by Equations (17.54) and (17.55) with \(a=1\).
For collisionless systems such as dark-matter halos, the fluid description so far does not hold, because there are no collisions to lead to pressure gradients. However, as we discussed in Section 17.1.2, if the velocity dispersion is not small, the velocity dispersion acts as an anisotropic pressure. In our derivation of the equivalent of the Euler equation (17.22) for collisionless, cold dark matter, we ignored the velocity dispersion term in Equation (17.21). But keeping this term and assuming for simplicity an isotropic, constant velocity dispersion (which we refer to as isothermal), we would eventually arrive at a version of Equation (17.24) that is the same as Equation (17.42) when we set \(\vec{v} \rightarrow \langle \vec{v}\rangle\) and \(c_s \rightarrow \sigma\). Thus, an isothermal collisionless system is in many ways equivalent to a fluid with a sound speed given by the velocity dispersion. In Fourier space, the evolution of perturbations is then given by the fluid form of Equation (17.53), with the Jeans wavenumber defined now as \begin{equation}\label{eq-dmhaloform-jeansk-collisionless} k^2_J = \frac{4\pi\,G\,\overline{\rho}\,a^2}{\sigma^2}\,, \end{equation} and the Jeans length as \begin{equation}\label{eq-dmhaloform-jeanslength-collisionless} \lambda_J = \frac{2\pi a}{k_J} = \sigma\,\sqrt{\frac{\pi}{G\,\overline{\rho}}}\,. \end{equation} When considering a collisionless system in a static, non-expanding background, again set \(a=1\) in these expressions.
Because of the correspondence between fluids and collisionless systems with non-zero velocity dispersions, the qualitative description from the previous paragraphs applies to collisionless systems as well: Instabilities occur on spatial scales similar to or larger than the Jeans length, while perturbations on scales smaller than the Jeans length oscillate or dampen. In fact, the effect of self-gravity is such that oscillations are typically damped and it is hard to find systems with undamped oscillations (but not impossible; Mathur 1990). The damping of oscillations in this context is known as Landau damping and it proceeds through the interaction of the perturbation and bodies in the collisionless system that are in resonance with the perturbation. This means that the difference between their orbital frequencies and the frequencies present in the perturbation are commensurate (see the discussion of resonant orbits in Section 13.1) . Physically, what happens is that these bodies “see” the perturbation always at the same phases of their orbits and, thus, strongly interact with it. While in principle it is possible for the bodies to both transfer energy to and from the perturbation, for realistic galactic distribution functions, the net energy transfer is always from the perturbation to the bodies and the perturbation is therefore damped. Effectively, the energy contained within a simple, ordered perturbation is dissipated as random energy of the bodies. Like phase mixing, Landau damping is important in establishing the original quasi-equilibrium of the dark-matter halo and in restoring equilibrium after subsequent perturbations.
Using the scalar virial theorem from Equation (14.18), we can show that the Jeans length for a self-gravitating collisionless system in equilibrium is always similar to the size \(R\) of the system: \(\lambda_J \sim R\). Assuming an absence of rotation, the scalar kinetic energy \(K \approx M\,\sigma^2/2\), where \(M = 4\pi \overline{\rho}R^3/3\) is the mass of the system, and approximating the system as a homogeneous density sphere we have that the scalar potential energy is \(W = -3 GM^2/[5R]\) (Equation 14.31). Using the scalar virial theorem \(2K+W = 0\), we then find that \(\sigma^2 \approx 4\pi G \overline{\rho}R^2/5\) and, thus, \(\lambda_J \approx 2\pi R/\sqrt{5} \approx 2.8 R\). Thus, the Jeans length of an equilibrium collisionless system is similar to the size of the system. This makes perfect sense, because the gravitational instability that creates the system must evolve to create a kinematically-hot system in such a way that it becomes insensitive to perturbations on scales smaller than the size of the system (or it would not be in equilibrium very long). The initial state of a dark-matter halo as it is collapsing is kinematically cold (see, e.g., the spherical collapse picture of non-crossing shells) with \(\sigma \approx 0\) and thus unstable on all scales. The growth of instabilities and their subsequent decay through phase mixing and Landau damping leads to a stable, virialized halo.
The importance of the Jeans length in the study of the stability of fluids and collisionless systems was first derived by Jeans (1902). To derive the Jeans stability criterion, Jeans (1902) considered an infinite, homogeneous system. Unfortunately, this system has appeared to be ill-defined, because in his derivation, Jeans (1902) had to ignore the gravitational potential generated by the smooth background density and posit that the Poisson equation only relates the perturbed potential and density. This is known as the Jeans swindle. We will not waste time discussing this mathematical curiosity, both because there are proper mathematical limits where the Jeans swindle is not necessary (Kiessling 2003) and, more importantly, because taking into account the background expansion removes the necessity of the Jeans-swindle assumption (e.g., Falco et al. 2013). This is clear in our own derivation of the Jeans length above in the cosmological context, where by working in comoving coordinates, it is natural that the Poisson equation describes the relation between the perturbed potential and density (Equation 17.17).
Our discussion of the stability of self-gravitating systems has been qualitative and even the mathematical derivation by Jeans (1902) only provides a qualitative understanding of real systems, which are far from infinite and inhomogeneous. Unfortunately, while they are mathematically better behaved than infinite, homogeneous systems (that is, no need for the Jeans swindle), more realistic systems such as spherical distribution functions or razor-thin disks are substantially more difficult to analyze. We will come back to this in Chapter 20, where we derive the stability criterion for stellar disks in detail, but to finish our discussion here, we state some useful results from the theory of the stability of spherical systems here without proof.
We discussed equilibrium distribution function models for spherical systems in Chapter 5.6 and discussed various models that can realistically represent galaxies, dark-matter halos, and stellar clusters. However, one question we did not ask is whether these models are stable against perturbations. In our discussion of spherical distribution functions, we distinguished between ergodic (isotropic) and anisotropic models and it turns out that isotropic models are generally fully stable against perturbations. This is because of a combination of two theorems, which discuss the stability of collisionless systems against radial and non-radial perturbations, in a decomposition of the angular dependence of perturbations in terms of spherical harmonics \(Y_l^m(\theta,\phi)\) (see Equation 12.63 in Chapter 12.3.1), where radial perturbations have \(l=0\) and non-radial ones \(l > 0\). For non-radial modes we have
Antonov’s second law (Antonov 1962b): Spherical collisionless systems with ergodic equilibrium distribution functions \(f(E)\) are stable against all non-radial perturbations if \(\mathrm{d} f(E) / \mathrm{d} E < 0\).
Similarly, for radial modes we have
Doremus-Feix-Baumann theorem (Doremus et al. 1971): Spherical collisionless systems with ergodic equilibrium distribution functions \(f(E)\) are stable against all radial perturbations if \(\mathrm{d} f(E) / \mathrm{d} E < 0\).
Thus, spherical systems with ergodic distribution functions are stable against all perturbations, as long as the distribution function decreases with increasing energy. This is typically the case and is in particular the case for all of the ergodic distribution functions that we discussed in Chapter 5.6.
For anisotropic systems, stability is not as easy to attain. For radial perturbations, a simple extension of the Doremus-Feix-Baumann theorem gives that anisotropic spherical systems are stable against all radial perturbations if \(\partial f / \partial E < 0\). However, there is no general result ensuring the stability of anisotropic spherical systems against non-radial perturbations. In fact, analytical (Antonov 1973) and numerical (Henon 1973) studies of spherical systems consisting preferentially of radial orbits—with constant orbital anisotropy \(\beta \gg 0\)—were found to be violently unstable through an instability known as the radial-orbit instability. For models with spatially-variable anisotropy, Polyachenko & Shukhman (1981) demonstrated this instability exists when the ratio of the kinetic energy in radial motions to that in tangential motions exceeds a threshold value of \(\approx 1.75/2\) (where an isotropic system has a value of \(1/2\)). In particular, systems with orbital anisotropy of the Osipkov-Merritt type from Chapter 5.6.4.2 are unstable if the anisotropy radius \(r_a\) is small. Because the formation of dark-matter halos starts out with the close-to-radial collapse of close-to-spherical shells, it is clear that the radial-orbit instability is highly important in the virialization process of dark-matter halos.
17.3.4. Violent relaxation¶
In the spherical-collapse model, dark-matter halos form when a spherical region of the Universe becomes so overdense with respect to the background density that it gravitationally decouples from the background expansion—known as “turn-around”—and spherical shells collapse towards the center without crossing each other. In reality, different shells of matter cross eventually—known as shell crossing—and the spherical symmetry is broken by the inhomogeneity of the local density distribution extending out to many Mpc (e.g., for the Milky Way, a significant force is the pull from the Virgo cluster, \(\approx 20\,\mathrm{Mpc}\) away). The roughly spherical initial collapse has high radial orbital anisotropy and the departure from spherical symmetry of the environment constitutes a non-radial perturbation. As such, the radial-orbit instability discussed in the previous section is an important ingredient in the evolution of dark-matter halos towards virialization. Numerical simulations of the evolution of spherical systems susceptible to the radial-orbit instability find that the end result of the instability is the formation of a strongly triaxial structure (e.g., Barnes et al. 1986). This is the likely explanation of why dark-matter halos in dissipationless simulations are strongly triaxial (see Chapter 12.1).
In the initial stage of collapse before shells cross and before the radial-orbit instability significantly affects the system, the gravitational potential felt by each gravitational body is constant and its specific energy is conserved. When shells start crossing and the radial-orbit instability leads to the breaking of spherical symmetry, the gravitational potential starts to evolve in time and the energy of individual bodies is no longer a conserved quantity, even though the total energy of the system is still conserved. This causes some bodies to lose energy and some to gain energy and the net result is a broadening of the energy distribution that drives the system towards a virialized state. A broadening of the energy distribution leads to a broadening of the distribution of orbital frequencies and phase-mixing therefore also becomes more rapid, leading to a relaxed phase-mixed state. Lynden-Bell (1967) termed this process violent relaxation, because, as he argued, it occurs on a dynamical time and is therefore fast compared to other collisionless relaxation processes (and certainly compared to collisional relaxation for galactic systems, which has \(\tau \approx N\,t_\mathrm{dyn}\)).
Lynden-Bell (1967)’s argument goes as follows. In a time-dependent gravitational potential, the time scale \(\tau_{\mathrm{vio}}\) for changing the energy is given by \begin{equation}\label{eq-dmhaloform-taurelax-1} \tau_{\mathrm{vio}} = \left\langle{E^2 \over \left({\mathrm{d} E / \mathrm{d} t }\right)^2 }\right\rangle^{1/2} = \left\langle{E^2/m^2 \over \left({\partial \Phi / \partial t}\right)^2 }\right\rangle^{1/2}\,, \end{equation} where \(E\) is the energy and we have used Equation (3.12) relating the time derivative of the energy to the partial derivative of the gravitational potential with respect to time. On its way towards virial equilibrium, the kinetic and potential energies of all bodies will oscillate with respect to their equilibrium value. As a whole, a self-gravitating system of bodies in virial equilibrium satisfies the virial theorem in the form of Equation (5.18), which states that twice the kinetic energy of the system equals the total potential energy of the system. Because the potential energy of the system as a whole is half of the sum of the individual potential energies, fluctuations in the potential energy are \(\approx 1/4\) of fluctuations in the kinetic energy; the total energy of bodies therefore oscillates by \(\approx 3/4\,m\Phi\). Plugging this estimate for \(E\) into Equation (17.111), we get \begin{equation}\label{eq-dmhaloform-taurelax-2} \tau_{\mathrm{vio}} = {3\over 4}\,\left\langle{\Phi^2 \over \left({\partial \Phi / \partial t}\right)^2 }\right\rangle^{1/2} = {3\over 4}\,\left\langle{1 \over \left({\partial \ln \Phi / \partial t}\right)^2 }\right\rangle^{1/2}\,. \end{equation} The relaxation time scale is therefore approximately that over which the potential fluctuates. As the collapsing dark-matter halo goes through shell crossing and the radial-orbit instability does its damage, the potential fluctuates on a time scale that is roughly the dynamical time from Equation (2.28), so \begin{equation} \tau_{\mathrm{vio}} \sim t_\mathrm{dyn}\,. \end{equation}
Potential fluctuations on the orbital time scale are able to inflict maximum damage on a gravitational system, because the perturbation’s time dependence is well matched with that of the orbits of many of the bodies in the system. Potential fluctuations on much longer time scales than the dynamical time do little to change the dynamical structure of a dynamical system, as the system is able to adjust its orbits without qualitatively altering them if the fluctuation is slow with respect to the dynamical time (this is known as adiabatic invariance). For slow perturbations, as soon as the perturbation goes away, the system goes back to its original state. For perturbations to the potential on the dynamical time scale, however, orbits are significantly and permanently altered by the perturbation.
To see the effect of a time-dependent potential on the energy distribution of a set of bodies, we can consider the following toy example. We set up a system of bodies with energies near the median of an NFW halo that is similar to the Milky Way’s and let it evolve without self-gravity in a fluctuating version of the underlying NFW potential, where the potential fluctuates between 50% and 150% of its equilibrium value on a time scale \(t_{\mathrm{osc}}\). To set \(t_{\mathrm{osc}}\), we estimate the dynamical time using Equation (2.28) for all bodies and set \(t_{\mathrm{osc}}\) to some fraction \(f\) of the dynamical time of the tenth percentile of the dynamical times in the sample. Because this toy example lacks self-gravity and instead applies the perturbation as an external force, the total energy of the \(N\) bodies is not conserved and the evolution of the system is not self-consistent. Nevertheless, we can qualitatively understand how the energy distribution evolves. We start by simulating the system for \(3\,t_{\mathrm{osc}}\) for \(f=1\) and show the result in Figure 17.7.
[14]:
from galpy.potential import (NFWPotential,
DehnenSmoothWrapperPotential)
from galpy.df import isotropicNFWdf
def violent_relaxation_demo(tosc_fac=1.,nsamp=1000,bins=31,
xmin=-1.39,xmax=-0.31,yScale='log',
ymin=0.,ymax=299.,
gcf=False,label=None):
np= NFWPotential(mvir=1.,conc=8.37)
dfh= isotropicNFWdf(pot=np,rmax=np.rvir())
dfh.turn_physical_off() # work in internal units
orig_samp= dfh.sample(n=int(1e6))
samp= orig_samp[numpy.fabs(orig_samp.E(pot=np)
-numpy.nanmedian(orig_samp.E(pot=np))) < 0.005]
# Use nsamp max
samp= samp[:numpy.amin([len(samp),nsamp])]
tosc= np.tdyn(numpy.percentile(samp.r(),10))/tosc_fac
halfpot= np/2.
oscpot= np\
+DehnenSmoothWrapperPotential(pot=halfpot,tform=0*tosc/2,tsteady=tosc/2)\
+DehnenSmoothWrapperPotential(pot=-1*halfpot,tform=1*tosc/2,tsteady=2*tosc/2)\
+DehnenSmoothWrapperPotential(pot=halfpot,tform=2*tosc/2,tsteady=3*tosc/2)\
+DehnenSmoothWrapperPotential(pot=-1*halfpot,tform=3*tosc/2,tsteady=4*tosc/2)\
+DehnenSmoothWrapperPotential(pot=halfpot,tform=4*tosc/2,tsteady=5*tosc/2)\
+DehnenSmoothWrapperPotential(pot=-1*halfpot,tform=5*tosc/2,tsteady=6*tosc/2)
ts= numpy.linspace(0.,6*tosc/2,3001)
samp.integrate(ts,oscpot,method='dop853_c')
if not gcf: figure(figsize=(13,5))
subplot(1,3,1)
# Time 0
line = hist(samp.E(),bins=bins,density=True,histtype='step',lw=2.)
yscale(yScale)
xlabel(r'$E$')
# Time 0.5
hist(samp.E(1*tosc/2),bins=101,density=True,histtype='step',lw=2.,color=line[2][0].get_edgecolor())
yscale(yScale)
xlabel(r'$E$')
xlim(xmin,xmax)
ylim(ymin,ymax)
galpy_plot.text(r'$t = 0.5\,t_{\mathrm{osc}}$',top_right=True,size=18.)
# Also label the t=0 line
if not gcf:
plot([-0.625,-0.45],[20.,10],color='k',ls='-')
galpy_plot.text(-0.5,7.,r'$t = 0$',size=18.)
subplot(1,3,2)
hist(samp.E(3*tosc/2),bins=101,density=True,histtype='step',lw=2.,label=label)
yscale(yScale)
xlabel(r'$E$')
xlim(xmin,xmax)
ylim(ymin,ymax)
galpy_plot.text(r'$t = 1.5\,t_{\mathrm{osc}}$',top_right=True,size=18.)
gca().tick_params(labelleft=False)
subplot(1,3,3)
hist(samp.E(6*tosc/2),bins=101,density=True,histtype='step',lw=2.)
yscale(yScale)
xlabel(r'$E$')
xlim(xmin,xmax)
ylim(ymin,ymax)
galpy_plot.text(r'$t = 3.0\,t_{\mathrm{osc}}$',top_right=True,size=18.)
gca().tick_params(labelleft=False)
return None
# Fast
violent_relaxation_demo(tosc_fac=1.,nsamp=5000,bins=31,yScale='log',label=r'$t_{\mathrm{osc}} = t_\mathrm{dyn}$')
# Faster
violent_relaxation_demo(tosc_fac=10.,nsamp=5000,bins=31,yScale='log',gcf=True,label=r'$t_{\mathrm{osc}} = t_\mathrm{dyn}/10$')
# Slow (adiabatic)
violent_relaxation_demo(tosc_fac=1./5.,nsamp=5000,bins=31,yScale='log',gcf=True,label=r'$t_{\mathrm{osc}} = 5\,t_\mathrm{dyn}$')
subplot(1,3,2)
handles, labels = gca().get_legend_handles_labels()
new_handles = [Line2D([], [], c=h.get_edgecolor()) for h in handles]
legend(handles=new_handles,labels=labels,frameon=False,fontsize=18.,loc='upper left')
tight_layout();
Figure 17.7: Violent relaxation.
We see that the initial narrow energy distribution quickly broadens and after a few oscillations reaches a quasi-steady state where the distribution retains its width. If the potential fluctuates on a time scale that is ten times shorter, the evolution is similar. Perturbing a system on time scales much smaller than the dynamical time does not increase the amount of relaxation induced in the system as can be seen by comparing the \(t_{\mathrm{osc}} = t_\mathrm{dyn}\) and \(t_{\mathrm{osc}} = 10\,t_\mathrm{dyn}\) histograms in Figure 17.7. If we instead perturb the system on a time scale that is much longer than the dynamical time, we see that while the energy of all bodies changes due to the perturbation, the system changes coherently and the energy distribution remains about as narrow as it was initially, as shown by the \(t_{\mathrm{osc}} = 5\,t_\mathrm{dyn}\) line in Figure 17.7.
Violent relaxation is a self-limiting process. Violent relaxation requires a time-dependent potential to operate, but drives self-gravitating systems towards a quasi-equilibrium state in which the potential no longer changes in time and violent relaxation therefore ceases. Simulations of the formation of dark-matter halos show that violent relaxation occurs, but is never complete, in the sense that the final energies of particles are correlated with their initial energies (e.g., van Albada 1982). The system therefore retains some memory of its initial condition.
The radial-orbit instability and the vaguer and broader concept of violent relaxation together with the other collisionless relaxation process that we described in the last few sections are part of the physical picture of what leads to the formation of virialized dark-matter halos with an end state that follows the NFW profile introduced in Chapter 2.4.6. However, the non-linear dynamics involved in the collapse and virialization of dark-matter halos formed from a cosmological density field are so complex that they need to be studied using \(N\)-body simulations. The discussion here has intentionally focused on providing some qualitative physical insight into what happens during this process, rather than discussing more detailed explanations for the form of the final density profile (but see, e.g., Lithwick & Dalal 2011; Pontzen & Governato 2013).
.
.
.