Dark matter haloes and subhaloes
Jesús Zavala Carlos S. Frenk
astro-ph.CO hep-ph







§ INTRODUCTION
 

The current theory of the formation and evolution of cosmic
structure in the Universe is based on the dark matter hypothesis in
which ∼84% of the mass-energy density of the Universe
<cit.> is in the form of a new type of particle, or
particles, with negligible electromagnetic interactions. The evidence
for the existence of dark matter is varied and compelling. It comes
from cosmic structures on all scales and across all epochs: from the
smallest, dark-matter-dominated dwarf galaxies
(e.g. <cit.>), through the largest clusters of galaxies
(e.g. <cit.>), to the large-scale structure of the Universe
(e.g. <cit.>) and back to the very seeds of cosmic
structure reflected in the temperature of the cosmic background
radiation (CMB; e.g. <cit.>). This body of evidence,
accumulated over the past three decades, can be accounted for within a
coherent theory of structure formation in which the gravity of the
dark matter amplifies primordial density perturbations imprinted
during an early period of cosmic inflation <cit.>. Empirical evidence for the existence of dark matter
comes purely from its gravitational effect: despite significant
efforts, experimental searches for dark matter particles in
accelerators (e.g. for a review in LHC searches see <cit.>), and
dedicated detectors on Earth (e.g. <cit.>) and in
space (e.g. <cit.>) so far remain unsuccessful.
Until the particles are discovered, dark matter will remain a
hypothesis, albeit one with strong empirical support.

In addition to the dark matter hypothesis, the standard theory of
structure formation makes a specific assumption about the nature of
dark matter, which is only partially supported by observations. This
is that the dark matter consists of classical, non-relativistic,
collisionless particles which had negligible thermal velocities at
early times.  This “cold dark matter” (CDM) is assumed to behave as
a fluid throughout most of the Universe's history, except at very
early times when this assumption breaks down in different ways
depending on the specific mechanism of dark matter production. The
most common hypothesis is that the dark matter particles are thermal
relics from the Big Bang (e.g. <cit.>). In this case,
dark matter was symmetric[Equal amounts of dark matter and
  anti-dark matter.]  and in thermal equilibrium with the
photon-baryon plasma through interactions with standard model
particles. As the Universe cools down, dark matter decouples from the
standard model particles, its creation and annihilation stops and its
co-moving density freezes out. If the strength of the
interactions is assumed to be on the scale of the weak force, then the
thermal relic abundance of these weakly interacting massive particles
(WIMPs) is quite close to the observed abundance of dark matter. This
remarkable coincidence, haplessly known as the WIMP miracle, has
enshrined WIMPs as the most popular dark matter candidates, especially
since new physics at the weak scale (and with them the emergence of
WIMP-like particles) was anticipated by Supersymmetric theories in
order to solve the hierarchy problem
(e.g. <cit.>). Moreover, WIMPs are the quintessential CDM
candidate because once they decouple, they are nearly collisionless
and, since they are massive (∼10 GeV – 1 TeV), they behave as a
classical (non-quantum) fluid that becomes non-relativistic very early
on.

The combination of the WIMP miracle with the success of the CDM
model in explaining the observed large-scale structure of the Universe
in the mid 1980s <cit.> established the current paradigm of
structure formation in which gravity is the only dark matter
interaction. This model has been widely adopted by the community
working on galaxy formation and evolution and, as a result, most of
our understanding on how cosmic structure emerges comes from studies
that assume the CDM model. This is a relevant remark in the context of
this review because the properties of dark matter haloes and their
substructure depend on the nature of dark matter (see
Section <ref>). In reality, the range of allowed
dark matter models, motivated to varying degrees by particle physics
considerations, is vast.  In this landscape of models, only a fraction
fall in the CDM category alongside WIMPs, e.g.  the QCD axion
(motivated by a proposed solution to the strong CP problem in particle
physics <cit.>).

Dark matter could become non-relativistic at sufficiently late times
so as to suppress, by free streaming, the formation of low-mass
galactic-scale haloes. This case is, in fact, one of the best studied
alternatives to CDM, known as warm dark matter (WDM). In contrast to
WIMPs, these particles have masses of 𝒪(1 keV). A sterile
neutrino, included as a part of a model that accounts for neutrino
masses and for the baryon asymmetry of the Universe, is the favourite
WDM candidate (for a recent review see <cit.>).  Another
possibility is that dark matter is made of extremely light bosons with
a 𝒪(1 kpc) de Broglie wavelength, in which case quantum
effects would be relevant on galactic scales (such possibility falls in
the category of “fuzzy dark matter”; for a review see
<cit.>.)  

Although the interactions between dark matter and standard model
particles are severely constrained, the interactions among the dark
matter particles themselves are not. It is possible that dark matter
may have its own rich phenomenology hidden from the ordinary
matter. This hidden dark matter sector might possess new forces
and particles, some of which could be viable dark matter particles
that are strongly self-interacting[By strong, we mean that
  the cross-section for self-interaction is of the order of the
  nuclear cross-section for visible matter (set by the strong
  force).]. These collisional particles fall under the category of
self-interacting dark matter (SIDM; for a review see
<cit.>[Some SIDM models are motivated by the baryon
  asymmetry; in these models, dark matter, unlike traditional WIMPs,
  shares this asymmetry (for a review of asymmetric dark matter see
  <cit.>).]). Some of the hidden particles might be light
enough that they effectively act as dark radiation that prevents
the gravitational collapse of dark matter on subgalactic
scales[In contrast to WDM, the damping of small structures is
  not due to free-streaming, but to a collisional, Silk-like,
  damping.]  (e.g. <cit.>). As mentioned
earlier, the CDM hypothesis is only supported to some extent:
astronomical data allow a variety of models in which dark matter
behaves significantly differently from CDM.

The goal of this paper is to provide a review of the formation,
evolution and dynamics of dark matter haloes and subhaloes, as
revealed primarily by N-body simulations.  Although no
account of the properties of haloes based purely on gravitational
dynamics can be complete since baryonic processes play a
significant role in galaxy formation, and new dark matter physics could
also do so, we focus on the standard CDM paradigm of
structure formation in part because the subfield of cosmological
N-body simulations has historically been developed in this context,
and also because the emergence and properties of dark matter
structures are most simply understood in the context of
CDM. Alternative dark matter models with additional physical
ingredients to gravity, albeit appealing, are more complicated. In various parts of this
review, we will explore how different assumptions for the nature of
dark matter can lead to different predictions from CDM.



§ FORMATION OF DARK MATTER HALOES
 



 §.§ Initial conditions: the primordial power spectrum in the linear regime


A theory of structure formation aims to explain the
evolution of the Universe from a nearly homogeneous initial state,
with tiny matter density perturbations, δρ/ρ, seeded by
inflation, which grow to leave an imprint on the CMB (emitted at the
time of recombination, z∼1100, when
δρ/ρ∼10^-3), through the emergence of the
self-gravitating dark matter haloes where galaxies form
(δρ/ρ≫1), to the Universe we observe today
characterised by a web of filamentary large-scale structure
(δρ/ρ∼1).

The starting point is the end of cosmic inflation when dark matter
perturbations are predicted to have a nearly scale-invariant power
spectrum, Δ^2∝ k^3+n_s, where
Δ^2(k)=k^3P(k)/2π^2 is the dimensionless power spectrum, and
the spectral index n_s=0.965 <cit.>. The growth of
dark matter perturbations in the expanding Universe is driven by
self-gravity. As long as the perturbations are small,
δρ/ρ≪1 (the linear regime), this growth can be
described by linear perturbation theory in which each perturbation
evolves independently of all others.

Two important processes occur in the linear regime, which modify the
primordial power spectrum. The first (known as the Mészáros effect
<cit.>) operates during the period when the energy
density in the Universe is dominated by radiation: the growth of dark
matter perturbations on scales smaller than the horizon stagnates,
while super-horizon scales continue to grow. This situation pertains
until matter overcomes radiation as the dominant component of the
energy density, after which all perturbations grow at the same rate.
The transition introduces a characteristic scale in the power
spectrum, the size of the horizon at the time of matter-radiation
equality. On scales smaller than this, the power spectrum
flattens. The second important scale, a cutoff in the power spectrum,
is of non-gravitational origin and reflects the particle nature of
dark matter. The physical mechanism that imposes this cutoff is model
dependent. For thermal relics (like many WIMP models and certain types
of WDM), the mechanism is free-streaming, a form of collisionless
(Landau) damping, whose scale is given by the horizon size at the
epoch when the dark matter particles become non-relativistic; the more
massive the particle, the earlier this epoch, and thus the smaller the
(co-moving) free-streaming scale is[The (co-moving)
  free-streaming scale is given by:
  l_ fs=2ct_ nr/a_ nr[1+ ln(a_
      eq/a_ nr)], where t_ nr is the age of the
  Universe at the time when the dark matter particles become
  non-relativistic (at a temperature 3k_BT_ nr∼ m_χ c^2);
  a_ nr=1/(1+z_ nr) is the scale factor at t_ nr
  (a∝ t^1/2 in the radiation-dominated era); and
  a_ eq is the scale factor at the time of matter-radiation
  equality.], k_ fs=2π/l_ fs.  This is the best-known
cutoff mechanism, which has been traditionally used to classify dark
matter into three categories (where m_χ denotes the mass of the
particle): cold[For cold particles, we have assumed CDM
  WIMPs, which requires taking into account the kinetic decoupling
  temperature and epoch; specifically we took Eq. 43 of
  <cit.>. ] (m_χ∼100 GeV,
k_ fs∼2.5×10^6 h/Mpc); warm (m_χ∼1 keV,
k_ fs∼3.8 h/Mpc); and hot (m_χ∼ 30 eV,
k_ fs∼0.3 h/Mpc).

A different type of damping is collisional damping,
which 
prevents the gravitational collapse of small structures, resulting in
an effective cutoff in the power spectrum. An example is kinetic
coupling of WIMPs, which effectively keeps dark matter coupled to the
photon-baryon plasma until the Universe cools enough that the
interactions become inefficient, damping perturbations beyond a scale
in the range (2.6×10^5 - 1.2× 10^8) h/Mpc
<cit.>. Another example is collisional damping due to
interactions between dark matter and relativistic particles in the
early universe (either photons or neutrinos,
e.g. <cit.>, or, in non-standard models, dark
radiation in hidden dark sector models,
e.g. <cit.>). The relativistic particles
create an effective radiation pressure that counteracts the
gravitational collapse, driving oscillations in the density
perturbations, akin to the well-known baryon acoustic oscillations
(BAOs), but on much smaller scales; by analogy they are called dark
acoustic oscillations, DAOs [Note that acoustic oscillations are also present in WIMP-CDM models (e.g. <cit.>), but they occur at much smaller scales than in relevant hidden dark sector models where they can be of galactic scale.]. 

Once the Universe cools down, the dark
matter decouples from the relativistic particles, imprinting a
characteristic scale (the size of the sound horizon at the time of
decoupling) in the power spectrum, followed by a Silk-like damping
cutoff.

The main features of the clustered dark matter distribution during the
linear regime are illustrated in Fig. <ref>. On the largest
scales, not affected by the Mézáros effect, the power spectrum is
nearly scale-invariant, Δ^2∝ k^3+n_s; on smaller scales
it bends to increasingly shallower slopes. For CDM (black line), the
power spectrum remains featureless well below galactic scales. For
reference, a dark matter halo today hosting a typical dwarf galaxy
would have a mass ∼10^10 M_⊙, roughly corresponding to a
(co-moving) wavenumber ∼12 h/Mpc[We use
  M=4π/3ρ(π/k)^3, where ρ is the mean dark
  matter density today.]. Measurements of galaxy clustering on scales
larger than individual galaxies, together with constraints from the
flux spectrum of the Ly-α forest (e.g. <cit.>)
constrain the power spectrum to be like CDM to the left of the hashed
area in Fig. <ref>. On smaller scales the power spectrum
could have a damping cutoff due to either collisionless (as in WDM
models; red line) or to collisional (as in models with DAOs; blue
line) processes. We have not included the cutoff characteristic of
fuzzy dark matter models, but we note that it is also
oscillatory like the DAOs models (but due to quantum rather than
collisional effects; see e.g. Fig. 2 of <cit.>).



As long as the dark matter perturbations remain linear
(δρ/ρ≪1), they grow at a rate that does not depend on
their co-moving scale, Δ^2(k;t)∝ D^2(t), where D(t) is
the growth factor, which depends only on the mean density of matter
and dark energy (see e.g. <cit.>). Once the density
contrast is no longer small (δρ/ρ∼0.1), perturbation
theory breaks down since gravity couples perturbations on different
scales and their evolution can no longer be calculated as independent
modes. 



 §.§ The non-linear regime: N-body simulation methods


To follow the evolution of dark matter density perturbations
beyond the linear regime, a number of approaches are possible
depending on the problem of interest. (i) High order perturbation
theory which can be used to study the quasi-linear regime
(δρ/ρ≲1), particularly in a modern reformulation
such as the Effective field theory of large-scale structure
<cit.>.  (ii) Analytical models with
simplified assumptions for the growth, turnaround (i.e., decoupling
from the expansion of the Universe), collapse and virialization (i.e.,
the formation of a gravitationally self-bound structure) of individual
perturbations. The best known examples are the Spherical
  collapse <cit.> and Ellipsoidal collapse
<cit.> models which link a primordial perturbation to the
final equilibrium configuration: the dark matter halo.  (iii) The halo model (for a review see <cit.>), which combines the
analytical models in (ii) with the assumption of a Gaussian density
field and can be used to compute the abundance of virialised haloes as
a function of their mass (the halo mass function); together with a
model for the dark matter distribution within haloes, it can be used
to model the non-linear dark matter power spectrum on all
scales. (iv) Models based on the Stable clustering hypothesis
(<cit.> see also <cit.>), which assumes that the
number of neighbouring dark matter particles within a fixed physical
separation remains constant, and can be used to study the deeply
non-linear regime; a recent re-formulation in phase space has been
shown to be a promising alternative to the halo model
<cit.>.  (v) Numerical N-body simulations, which solve ab initio the gravitational
evolution in phase space of a distribution of N particles sampled
from an initial power spectrum. This is the most general and powerful
approach to study the clustering evolution on all scales and is the
focus of this review.  (vi) Techniques that avoid the particle
discretisation inherent in N-body simulations by following the
phase-space distribution function directly <cit.>. These are
particularly useful to study evolution from truncated power spectra
such as for hot or warm dark matter for which standard N-body
techniques suffer from artifical fragmentation <cit.>. 












In the case of classical, non-relativistic, collisionless particles, i.e., CDM, 
N-body simulations follow the evolution of the dark matter
phase-space distribution function, f(x⃗, v⃗; t), which in principle is given by the
collisionless Boltzmann equation coupled with the Poisson equation
for the gravitational field, Φ(x⃗) (a combination known as
the Vlasov-Poisson equation):

    df/dt   =   ∂ f/∂ t+∑_iv_i∂ f/∂ x_i-∑_i∂Φ/∂ x_i∂ f/∂ v_i=0
    ρ_χ(x⃗;t)   =   ∫ f(x⃗,v⃗; t)d^3v⃗
    ∇^2Φ(x⃗)   =   4π Gρ_χ(x⃗)

where d/dt is the Lagrangian derivative. Cosmological N-body
simulations[For a review see e.g. Section 3 of
  <cit.>.] solve this equation in an expanding universe
using a co-moving reference frame (with the expansion included
explicitly through the solution of the Friedmann equations for the
scale factor), discretizing the distribution function as an ensemble
of N phase-space elements or “particles”,
{x⃗_⃗i⃗,v⃗_⃗i⃗}, with i=1,...,N. Since the collisionless
Boltzmann equation implies that the phase-space distribution remains
constant in time along any trajectory {x⃗(t),v⃗(t)}, the
distribution obtained by following the N particles from initial
conditions sampled from the phase-space distribution at t=0,
constitute a representative Monte-Carlo sampling of the distribution
function at any subsequent time, t. The N particles are thus a
statistical representation of the coarse-grained[By this we
  mean an average of the fine-grained distribution function in the
  collisionless Boltzmann equation over the scales resolved in the
  simulation, typically several times the interparticle separation.]
distribution function:

    f̃(x⃗,v⃗)   ∼   ∑_im_iW(|x⃗-x⃗_⃗i⃗|;ε)δ^3(v⃗-v⃗_⃗i⃗);       df̃dt=0
    ρ̃(x⃗)   =   ∫f̃(x⃗,v⃗)d^3v⃗∼∑_im_iW(|x⃗-x⃗_⃗i⃗|;ε)
    
    Φ̃(x⃗)   =   ∫ g(x⃗-x⃗')ρ̃(x⃗')d^3x⃗'

where m_i is the mass of the simulation particle,
δ^3(v⃗-v⃗_⃗i⃗) is the Dirac delta funcion in 3D, W is
a kernel density with a softening length
ε[In principle each particle can have an
  individual softening, see e.g. Section 4 of <cit.>.],
introduced to obtain a smooth density field from the set of N
discrete particles; i.e., the kernel effectively models each
simulation particle as an extended mass distribution[The
  introduction of a softening scale in the density (or potential)
  suppresses gravitational two-body large-angle scatterings which are
  artificial for an approximately continuous dark matter density
  distribution.]; finally, the last equation for the potential is the
general solution to Poisson's equation as a convolution of the density
field with a suitable Green's function[In Fourier
    space, Eq. <ref> is simply a multiplication
    Φ̂̃̂(k⃗)=ĝ(k⃗)ρ̂̃̂(k⃗). ].


Since each simulation particle represents a region of phase space containing
a very large number, m_i/m_χ, of real dark matter particles, the
information in an N-body simulation is always incomplete, limited by
the phase space resolution and the softening length.

With the discretization method employed in an N-body simulation,
calculating the evolution of the phase-space distribution is reduced
to following self-consistently the dynamics of a system of N
particles (usually in terms of the Hamiltonian of the system in the
co-moving frame) according to the potential derived from the particle
distribution.  Modern codes used to solve this problem employ
efficient methods for computing the gravitational potential and
integrating the Hamiltonian system forward in time. Early cosmological
simulation codes used the particle-mesh (PM) technique in Fourier
space (e.g. <cit.>) or direct integration of the
N^2 interactions (e.g. <cit.>). The former is limited in
resolution by the size of the mesh while the latter is limited by
speed. These two shortcomings can be overcome by combining both
techniques in the P^3M method
(e.g. <cit.>), in which the long-range
forces acting on a particle are calculated on a PM grid and the
short-range forces by direct N^2 summation. An alternative approach
is the hierarchical tree method <cit.> in which an octree
is used to divide the volume recursively into cubic cells and
increasingly coarse cells are used to compute the forces on a particle
at increasingly large distances. The most widely used cosmological
simulation code is GADGET-2 <cit.>, which uses the treePM
algorithm, whereby short-range forces are computed with the tree
method and long-range forces with Fourier techniques[For a
  review of the force computation methods see Section 3.5 of
  <cit.>.].

If dark matter cannot be treated as CDM, then the fundamental
equations may need to be modified. For models that only deviate from
CDM because of a cutoff in the initial power spectrum (such as hot or
warm dark matter and certain DAO models), the N-body equations
(<ref>-<ref>) and methods used for CDM are still
valid as long as the dark matter behaves as a collisionless, classical
system, and the simulation starts well after the dark matter particles
have become non-relativistic ; all that is needed is a modification of
the initial conditions (see <ref> below).  On the other hand, if
dark matter is non-relativistic but no longer collisionless, like in
SIDM, then the collisionless Boltzmann equation needs to be replaced
by the full collisional Boltzmann equation, which has an extra term
(the collisional operator) in the right-hand-side of eq. <ref>, to
account for the effect of dark matter collisions according to a
self-scattering cross-section. It is possible to incorporate this new
term within the Monte Carlo approach of traditional N-body
simulations by adding “collisions” between each simulation particle
and its immediate neighbours in a probabilistic way that reflects the
effective scattering rate given by the cross-section
(e.g. <cit.>; see
Appendix A of <cit.> for a detailed derivation). An
alternative to the N-body approach is the “gravothermal fluid”
approximation <cit.>, which considers an SIDM dark
matter halo as a self-gravitating, spherically symmetric, ideal gas
with an effective thermal conductivity (related to the self-scattering
cross-section, see e.g. <cit.>). Although this approach is
restricted, it provides physical insight into the evolution of SIDM
haloes, and a degree of validation of SIDM N-body
simulations. Finally, if quantum effects are important for the dark
matter fluid, then the Vlasov-Poisson equation needs to be replaced by
the Schrödinger-Poisson equation, whose solution requires
numerical methods quite distinct from the N-body approach
(e.g. <cit.>).



 §.§ The non-linear regime: initial conditions and the
  emergence of the cosmic web


The techniques of Section <ref> can be used to
integrate forward in time a particle distribution starting from an
initial state, the initial conditions, usually taken to be in the
linear regime described by perturbation theory. The basic techniques
for generating general initial conditions were laid out in
<cit.> and <cit.> and have been refined over
the years (e.g. <cit.>; for a review see
<cit.> or Appendix C1.1.4 of <cit.>). They provide a
particle realization with the statistical properties of the linear
dark matter density field described by the power spectrum. In general
the procedure can be divided into two steps:


  i) create a realization of an unperturbed cube of side L by
  distributing N particles homogeneously in a lattice or in a
  glass-like configuration[The particles are initially placed
    at random in the simulation cube and then left to evolve under a
    repulsive force by reversing the sign of the gravitational force
    until they reach an equilibrium configuration that has no
    discernible grid pattern <cit.>.] to avoid imprinting a
  grid-like pattern in the simulation.

  ii) perturbations of wavelength λ down to the Nyquist
  frequency of the particle distribution are represented by plane
  waves of spatial frequency in Fourier space, k=2π/λ, whose
  amplitudes and phases are drawn at random from a Gaussian
  distribution with variance proportional to the desired linear power
  spectrum. The density field and its gravitational potential in real
  space are then obtained by an inverse Fourier transform. Using the
  Zel'dovich approximation <cit.>, or the more accurate
  second-order Lagrangian perturbation theory
  (e.g. <cit.>), these fields are used to compute the
  displacements needed to transform the uniform N-particle
  distribution in part i) into a distribution that has the desired
  power spectrum.




An illustration of the end result of this procedure is shown in
Fig. <ref>. The initial conditions generator, MUSIC
<cit.>, was used to construct the particle distribution on
the right, which is a statistical realisation of the CDM linear power
spectrum shown on the left. The main limitations for a cosmological
simulation are already set in the initial conditions: the maximum
lengthscale that can be simulated is determined by the (co-moving)
side of the computational cube[ A sufficiently large volume
  is needed to sample large-scale modes that remain approximately
  linear during the simulation where power is transferred from large
  to small scales; without appropriate large-scale sampling, the
  clustering is no longer accurate once perturbations on the scale of
  the cube become non-linear.], and the minimum lengthscale that can
represented in the initial conditions is set by the Nyquist frequency
of the particle distribution[In practice, power below the
  Nyquist frequency is generated non-linearly so the resolution of the
  simulation is not limited by the Nyquist frequency but rather by the
  gravitational softening scale, ε.]. The choice of cube
length and particle number depends on the science goal of the
simulation and on the computing resources available. We will come back
to this point below.

The procedure illustrated in Fig. <ref> for CDM can
be readily applied to other dark matter models with different initial
power spectra. In fact, in models where dark matter only behaves
differently from CDM at very early times, e.g., in thermal-relic WDM
models, it is the different initial conditions (the lack of power on
small scales in WDM relative to CDM in the linear regime) that gives
rise to the main differences between these models (since the residual
thermal motions in WDM models of interest are negligible (see
e.g. <cit.>). In models with a truncated initial power
spectrum, the subsequent evolution is affected by particle
discreteness in the reconstruction of the density field which
introduces an irreducible (shot-noise) power. This results in spurious
clustering on scales close to the cutoff length <cit.> that
requires careful treatment to either remove small-scale artificial
clumps <cit.> or avoid their formation altogether by using
non-standard simulation techniques <cit.>.



Reproduced from Michael Boylan-Kolchin et al. Resolving cosmic structure formation with the Millennium-II Simulation. MNRAS (2009) 398 (3): 1150-1164, doi: 10.1111/j.1365-2966.2009.15191.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/398/3/1150/1261296link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com

Once the initial conditions are generated, an N-body simulation is
performed, most commonly in a computational cube with periodic
boundary conditions, to follow the evolution of the density and
velocity fields into the non-linear regime across all resolved scales.
An example, the Millennium II simulation <cit.>, is
illustrated in Fig. <ref>. The left set of panels shows the
projected dark matter density distribution at various snapshots
corresponding to the redshifts shown at the top right of each panel.  The emergence of the cosmic
web, the result of gravitational clustering, is apparent, with its
now familiar pattern of filaments over a range of scales surrounding
voids. The right panel shows the evolution of the power spectrum at
the same snapshots (solid lines). The hierarchical onset of non-linear
structure, from small to large scales is clearly apparent by
reference to the linear power spectrum (grey lines).

The first CDM cosmological N-body simulations in the 1980s
<cit.> already contained all the relevant physical
processes of gravitational clustering for collisionless dark matter,
but were computationally limited; they could follow the evolution of
only 𝒪(10^4) particles. In the decades since then, the
tremendous improvement in computational capabilities has been such
that cosmological (L≳100 Mpc/h) simulations with
𝒪(10^9) particles are routinely performed[For a
  review on the state of cosmological simulations circa 2012 see
  <cit.>.], and the most expensive simulations to date have
reached the one trillion particle milestone <cit.>.



To compare the simulations to astronomical data it is necessary to
make a correspondence between dark matter haloes and the galaxies that
would form within them. In the earliest simulations, galaxies were
identified with high peaks of the suitably filtered density field, an
assumption known as “biased galaxy formation”
<cit.>. Physically based models of galaxy
formation that could be grafted onto N-body simulations were
developed in the early 1990s <cit.>. These are known as
“semi-analytic models” because they encapsulate the relevant
physical models in a set of coupled differential equations that are
solved numerically. These equations assume spherical symmetry and
describe the cooling of gas, the formation of stars, chemical
evolution, the growth and merging of central supermassive black holes
and feedback effects arising from energy injected into the gas during
the course of stellar evolution and by active galactic nuclei
triggered by accretion of gas onto the central black hole. The model
is applied at every stage of the gravitational evolution of the
merging hierarchy of haloes, described by a merger
tree (see Section <ref>). Semi-analytic models have been extremely successful in linking
the distribution of dark matter computed in an N-body simulation to
the observed Universe <cit.> 
and have become very sophisticated in predicting visible galaxy properties
over a large range of wavelengths (e.g. <cit.>).

An example based on the Millennium simulation is shown in Fig. <ref>.

Dark matter N-body simulations are the cornerstone of the current
understanding of how galaxies form and evolve and, as illustrated in
Fig. <ref>, have been very successful in explaining the
large-scale structure of the Universe <cit.>. The latter
accomplishment is non-trivial and demands certain conditions about the
nature of dark matter. For instance, already in the 1980s, light
neutrinos were ruled out as the dominant component of dark matter by
their incompatibility with the observed large-scale structure
<cit.>, thus demonstrating the potential of N-body
simulations to test models for the nature of the dark matter.  By
contrast, the fact that CDM matched the observations available at the
time remarkably well contributed greatly to its establishment as the
standard model of cosmogony <cit.>. By now, it is firmly
established that whatever the dark matter is, it must behave as CDM on
large scales (see Figs. <ref> and <ref>). It is
important to recognize, however, that a wide range of dark matter
candidates behave just as CDM on large scales and thus are also
allowed by the large-scale structure data, as we discussed in
Sections <ref> and <ref>. In this sense, the
success of the CDM model in explaining the large-scale structure of
the Universe is shared by allowed WDM, SIDM and fuzzy dark
  matter models.



 §.§ The structural properties of dark matter haloes


As a consequence of gravitational clustering, the tiny density
perturbations present when the CMB was emitted grow in time,
eventually separating from the expansion of the Universe and becoming
self-gravitating bound structures known as dark matter haloes. This
process of forming virialised haloes can be understood from the
simple spherical collapse model <cit.>. Haloes become
increasingly more massive with time, smoothly by accreting mass
from their surroundings or merging with other, smaller
haloes. The latter thus become subhaloes, which is the topic of
Section <ref>. Although the large-scale environment and
spatial clustering of dark matter haloes are clearly relevant, here we
focus on the abundance (halo mass function) and internal structure of
dark matter haloes. These properties are the most useful when
attempting to differentiate dark matter models. Currently, however,
the halo mass function and halo structure on the key subgalactic
scales are only weakly constrained observationally. We should also
note that not all dark matter is contained within haloes. The fraction
of unclustered dark matter is naturally a strong function of
time reflecting the growth of collapsed objects by hierarchical
clustering. Even today at most ∼20% is expected to be
unclustered in the CDM model <cit.>. A recent update of
this work (applied also to WDM) puts the fraction even lower, at the
percent level <cit.>.

Definition of a halo.— Because of the dynamic nature of haloes
and their lack of spherical symmetry, precisely defining the boundary
of a halo, and thus its mass, is, to some extent, arbitrary
<cit.>. A variety of definitions exist
in the literature with the most common ones being: (i) the FOF
(friends-of-friends) mass, defined as the mass of the set of particles
that are linked together by a percolation scale, defined by a linking
length, b∼0.2 in units of the mean interparticle separation
<cit.>; (ii) a spherical overdensity mass, M_Δ,
contained within a sphere centred on the halo (with the centre placed
at the minimum of the gravitational potential of the halo), with a
radius given by the spherical collapse model, whereby the
collapsed region that defines a halo contains an average density
Δ(z) times the critical density for closure
<cit.>. The overdensity, Δ(z), is a redshift-dependent
function of cosmology <cit.>, but for the
Einstein-de Sitter cosmology, Δ∼178 at all times; (iii) the
viral mass, defined with Δ=200, which early simulations
identified as the radius that separates the region of the halo that is
in dynamical equilibrium from the surrounding region that is still
collapsing <cit.>.  Given the simplicity of the latter, its
relation to dynamical equilibrium, and its connection with the
Einstein-de Sitter spherical collapse overdensity, the radius,
r_200, and the enclosed mass, M_200, are widely used in the
field as the boundary and mass of dark matter haloes, respectively.

The halo mass function.— The mass function of dark matter
haloes, i.e., the number density of haloes of different mass, has been
characterised quite precisely in the last couple of decades by
N-body simulations
<cit.>,
and is now well determined over the full range of epochs and masses
relevant to galaxy formation, from 𝒪(10^8  M_⊙)
dwarf-size haloes to 𝒪(10^15  M_⊙) cluster-size
haloes.  The number density of haloes per unit mass scales as:

 
    dn/dM∝ M^α,    where  α∼-1.9    for low masses,

  with an overall normalisation that correlates with the large
 scale environment, with denser environments having a larger halo
 abundance <cit.>.



 The shape of the halo mass function is reasonably well understood
 from statistical arguments based on the properties of the initial
 Gaussian density field (described by the power spectrum) and the
 gravitational collapse of density perturbations into virialised
 haloes as modelled by the spherical collapse model. These arguments
 are the basis of the Press-Schechter <cit.> and extended
 Press-Schechter (EPS) formalisms <cit.>, which
 provide a good fit to the simulation results, particularly if the
 assumption of spherical symmetry for the collapse of overdensities is
 replaced by the assumption of ellipsoidal collapse <cit.>.
 At the small-mass end, the power-law form of the mass function is
 broken at a mass that depends on the nature of the dark matter. For
 example, a cutoff in the power spectrum, whether due to relativistic,
 collisional, or quantum effects, introduces a corresponding cut-off
 in the halo mass function. The mass function for WDM
 <cit.> and
 interacting dark matter
 <cit.> models have
 now been fairly well characterized with N-body simulations (with
 appropriate corrections for spurious fragmentation due to particle
 discreteness near the cutoff <cit.>). The
 Press-Schechter approach can be readily extended to provide a
 reasonable approximation to the halo mass function in these models as
 well <cit.>.

Fig. <ref> provides an example of the effect of
a cutoff in the primordial power spectrum on the halo mass function
relative to CDM. The “atomic dark matter model”, ADM_ sDAO
<cit.>, is an example of a model with dark acoustic
oscillations, while WDM is an well-known example of the free streaming
effect (see Fig. <ref>). These two models give rise to
qualitatively different types of suppression in the abundance of
low-mass haloes. The halo mass function thus contains a signature of
the type of primordial power spectrum cutoff.

The inner structure of dark matter haloes.— One of the
remarkable findings of the past few decades is that the spherically
averaged mass density profiles of dark matter haloes in dynamical
equilibrium have a nearly universal form which is independent 
of halo mass, initial conditions[This is only true
  if, on the scales of interest, the primordial power spectrum grows
  monotonically towards large k.]  and cosmological
parameters. These profiles are quite well described by a very simple
functional form with just two parameters, the so-called
Navarro-Frenk-White (NFW) profile <cit.>:
 
    ρ_ NFW(x)=ρ_s(c)/cx(1+cx)^2,

 where x=r/r_200, and c=r_200/r_s is the concentration of the
 halo; r_s is the scale length, which, for the NFW profile, coincides
 with the radius, r_-2, at which the logarithmic slope of the
 profile is equal to -2; finally
 ρ_s(c)=δ_cρ_ crit, where
 ρ_ crit=3H^2/8π G is the critical density of the
 Universe, and:

    δ_c=Δ c^3/3K_c(c)
  
where K_c(c)= ln(1+c)-c/(1+c). Although recent simulations have
shown that a different profile, the so-called Einasto profile, which
has three parameters, is a slightly better fit to simulations
<cit.>, the asymptotic behaviour of the NFW profile for
ρ(r→0)∼ r^-1 remains a remarkably good
approximation to the inner structure of CDM haloes (see top left panel
of Fig. <ref>). The physical origin of this divergent
cusp and the remarkably universal profile shape are not fully
understood.  It has been argued from N-body simulations of the early
stages of structure formation that the first CDM haloes to form, i.e.,
those near the free-streaming scale of CDM have a steeper cusp than
NFW, ∼ r^-1.5, which is subsequently flattened after a few
mergers to ∼ r^-1
<cit.>. More recent
simulations which follow the growth of the first mini-haloes all the
way to the present, seem to confirm this, suggesting that the
ubiquitous ∼ r^-1 slope develops at some point after the
formation of the halo and remains until z=0.


Halo concentration.—The remarkable simplicity of the NFW halo
density profile goes beyond Eq. <ref>: the profile is, in fact,
fully specified by a single parameter, halo mass, because the
concentration (or scale radius) correlates with mass, with lower-mass
haloes generally being more concentrated than higher-mass haloes
<cit.>. This
correlation is ultimately a consequence of the hierarchical nature of
structure formation by gravitational instability from a primordial
power spectrum that, as in CDM, grows monotonically towards small
scales (see Fig. <ref>). Lower-mass haloes form earlier,
when the mean density of the Universe is larger, and larger-mass
haloes form later when the mean density of the Universe is lower. The
inner regions of haloes collapse first <cit.> and their
density reflects the mean density of the Universe at that time. Hence,
smaller-mass haloes are more concentrated than larger-mass
haloes. Furthermore, since (at least for CDM), the power spectrum,
Δ^2(k), becomes flatter at larger k (due to the Mészáros
effect), haloes with a wide range of masses collapse in a short time
interval and this flattens the concentration-mass relation at low
masses. Models based on these simple arguments explain, at some level,
the concentration-mass relation measured in simulations
<cit.>, and also provide a natural
connection between the mass assembly of haloes in time and their
structure as a function of radius: each radial shell has the
characteristic density of the cosmic background density at the time
when it collapses <cit.>.  Random deviations around the
mean collapse time expected for haloes of a fixed mass and a
stochastic merger history introduce scatter in the concentration-mass
relation.



footnote-4
 footnoteReproduced from Volker Springel et al. The Aquarius Project: the subhaloes of galactic haloes. MNRAS (2008) 391 (4): 1685–1711, doi: 10.1111/j.1365-2966.2008.14066.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/391/4/1685/1747035link. 

 footnoteReproduced from Julio Navarro et al. The diversity and similarity of simulated cold dark matter haloes. MNRAS (2010) 402 (1): 21–34, doi: 10.1111/j.1365-2966.2009.15878.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/402/1/21/1028856link. 

 footnoteReproduced from Carlos Vera-Ciro et al. The shape of dark matter haloes in the Aquarius simulations: evolution and memory. MNRAS (2011) 416 (2): 1377–1391, doi: 10.1111/j.1365-2966.2011.19134.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/416/2/1377/1061105link. 
 
  footnoteReproduced from Mark Vogelsberger et al. Phase-space structure in the local dark matter distribution and its signature in direct detection experiments. MNRAS (2009) 395 (2): 797–811, doi: 10.1111/j.1365-2966.2009.14630.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/395/2/797/1747020link. 


The figures mentioned in footnotes 21-24 are not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com

Halo velocity distribution.— For a spherical,
self-gravitating, collisionless system in dynamical equilibrium, with
radial density profile, ρ(r), the Boltzmann equation reduces to
the well-known Jeans equation <cit.>:

    d(ρσ_r^2)dr+2β/rρσ_r^2=-ρdΦdr,

where Φ(r) is the gravitational potential related to the density
by Poisson's equation, σ_r(r) is the radial velocity dispersion
profile, and β(r)=1-(σ_θ^2+σ_ϕ^2)/σ_r^2
is the velocity anisotropy profile, which quantifies the degree of
anisotropy of the particle orbits in the dark matter halo. Haloes tend
to be isotropic only in their centres and are radially anisotropic in
their outskirts <cit.>, with a velocity
anisotropy that is related to the logarithmic slope of the density
profile <cit.>. The velocity structure of dark matter
haloes in equilibrium is thus intimately linked to their spatial
distribution (see top right panel of Fig. <ref>),
which is strikingly evident in the so-called pseudo-phase-space
  density, Q≡ρ/σ^3, where
σ^2=σ^2_r+σ_θ^2+σ^2_ϕ is the square of
the 3D velocity dispersion. This quantity has been found to be an
almost perfect power law, Q∝ r^-1.875, over several orders
of magnitude <cit.>, and is in remarkable
agreement with the self-similar solution for infall onto a point mass
in an Einstein-de Sitter universe <cit.>. The radial
behaviour of Q is a manifestation of the nearly universal structure
of dark matter haloes, which connects both their spatial and
kinematical distributions.

Besides having anisotropic particle orbits, CDM haloes have a 
non-Maxwellian velocity distribution. This may be seen in a simple way
by noting that for a purely isotropic spherical system (β=0), the
full velocity distribution function of the halo depends only on the
specific energy, f(ℰ), and is fully given
by the Eddington formula <cit.>:

    f(ℰ)=1/√(8)π^2∫_0^√(ℰ)dud^2ρ/dΨ^2(r(Ψ(u))),

where u = √(ℰ-Ψ) and ℰ and Ψ(r)
are the (negative) specific energy and gravitational potential,
respectively. If haloes were spherical and isotropic, their velocity
distribution would be given purely by the NFW density profile through
Eq. <ref>. Even at this level of approximation we can see
that haloes would not be described by a Maxwellian distribution
function (only the singular isothermal sphere, ρ∝ r^-2,
results in a Maxwellian distribution in Eq. <ref>). In
other words, haloes are non-Maxwellian simply by virtue of their NFW
density profiles. In fact, CDM haloes in simulations have local[Given the limited resolution of simulations, local
  in this sense refers to regions of at least 𝒪(10
  kpc^3) as in <cit.>.]  velocity distributions
that show significant departures from Maxwellian, related to the
dynamical assembly of the halo.  The features that appear in the local
velocity distribution are unique for a particular halo, and retain the
memory of its assembly history (<cit.>; see bottom
right panel of Fig. <ref>).

Halo shapes.— Although to first order, CDM haloes are well
described by the spherical NFW profile, they are, in fact, triaxial
<cit.>. In general, CDM haloes are prolate in the inner
parts and oblate in the outskirts
(<cit.>; see bottom left
panel of Fig. <ref>). This radial dependence of halo
shape seems to be related to the assembly history of the halo within
the cosmic web: the central regions, being assembled at earlier times
through accretion along narrow filaments, end up being more prolate,
while the outskirts, more recently assembled by less anisotropic
accretion end up more oblate <cit.>. Thus, the
halo shape profile at z=0 carries some memory of its assembly
history. Overall, more massive haloes are more aspherical than lower
mass haloes <cit.> because in the hierarchical
CDM model, the more massive haloes form more recently and thus their
shapes retain memory of the most recent accretion event
<cit.>.

Dependence on the nature of the dark matter.— There are
significant changes in the structure of dark matter haloes if the dark
matter particles do not behave as CDM. In currently allowed models, these
deviations are mostly confined to the central regions, that is within
the scale radius, r_s. By introducing a new scale in the process of
structure formation, be it in the initial conditions (e.g. a cutoff in
the linear power spectrum), or during the non-linear evolution phase
(e.g. a subgalactic scale mean free path due to self-interactions),
these models break the near universality of CDM haloes. 

In models with a galactic-scale cutoff in the primordial power
  spectrum, such as WDM and interacting DM, the main changes can be
understood from the later collapse of the first generation of haloes
in these models compared to CDM. In contrast to CDM, these galactic-size haloes are
not formed hierarchically from the assembly of smaller haloes but,
instead, by monolithic collapse.  Their characteristic density
therefore reflects the mean background density at the time of the
monolithic collapse. By contrast, 




footnote-2
 footnoteReproduced from Mark Vogelsberger et al. Subhaloes in self-interacting galactic dark matter haloes. MNRAS (2012) 423 (4): 3740–3752, doi: 10.1111/j.1365-2966.2012.21182.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/423/4/3740/1749150link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com

 footnoteReproduced from Mark Lovell et al. The haloes of bright satellite galaxies in a warm dark matter universe. MNRAS (2012) 420 (3): 2318–2324, doi: 10.1111/j.1365-2966.2011.20200.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/420/3/2318/979379link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com
 

a CDM halo of the same mass forms
from the assembly of smaller fragments that typically formed earlier
and are therefore denser. Simulations of WDM models have characterised
the spatial structure of WDM haloes quite accurately, showing that the
density profiles of allowed models are, in fact, well described by the
NFW profile but with a lower concentration at a given mass
<cit.>.
The concentration-mass relation for WDM haloes can then be modeled in
an analogous way to CDM, but taking into account the cutoff in the
power spectrum <cit.> which is therefore reflected in the
concentration of the haloes.  An example of this is shown in the
bottom left panel of Fig. <ref>, where the
density profile of a CDM halo is mapped into that of a 2.3 keV thermal
relic WDM halo by simply scaling down the concentration using Eq. 39
of Ref. <cit.>, which connects the concentration to the
cutoff scale in the power spectrum[This functional form has
  been corroborated by <cit.>, but the parameters in the two
  studies are different.  The formula is nevertheless a good
  approximation to the general behaviour.]. This lower concentration
is also reflected in the velocity dispersion profile (see bottom-right
panel of Fig. <ref>).

It is interesting to note that the pseudo-phase-space density,
Q≡ρ/σ^3, in WDM haloes scales with radius in the same
way as in CDM, Q∝ r^-1.9. In principle, Q can never exceed
its primordial value, Q=Q_ max, determined by the thermal
velocities of the unclustered dark matter particles
<cit.>. This is because for a collisionless system,
Liouville's theorem requires conservation of the fine-grained
phase-space density and the coarse-grained density, approximated by
Q, can never exceed this value. Thus, the central regions of WDM
haloes cannot exceed a maximum density, i.e., they form a central
density core. However, the value of Q_ max is so large in
allowed WDM models that the core size is tiny, 𝒪(10 pc)
for ∼ keV thermal WDM relics in dwarf-size haloes
<cit.>; WDM cores are thus irrelevant in
practice. Finally, WDM haloes are slightly less triaxial than CDM
haloes as a whole for a fixed mass, at masses near the cutoff scale
<cit.>.

In SIDM, collisions between the dark matter particles have an on 
impact in the inner structure of haloes once the timescale for
collisional relaxation at the characteristic radius of the halo r_s,
t_ rel, s∼(ρ(r_s) ⟨ v_ rel,
  s⟩σ_T/m_χ)^-1, where
⟨ v_ rel, s⟩ is the characteristic local relative
velocity, becomes comparable to the age of the inner halo. The
original CDM density cusp turns into a core within the region where
this condition is satisfied. The interaction cross-section thus
introduces a new scale in structure formation – the mean free-path
for particle collisions – which breaks the near universality of CDM
haloes. The transformation of the cusp into a core due to elastic
collisions at the halo centre is a transitory phase that leads to a
quasi-equilibrium configuration once the core has acquired its maximum
size, which is approximately the radius at which the velocity
dispersion profile peaks (see bottom panels in
Fig. <ref>). Prior to this, the transfer of
energy during elastic collisions proceeds from the outside in since
the velocity dispersion profile has a positive gradient in the inner
regions and so there is a net “heat flux” from the regions close to
the maximum of the velocity dispersion to the centre
(e.g. <cit.>). Once the core reaches its maximum size,
subsequent collisions can only result in a net heat flux from the
inside out since the velocity dispersion profile has a negative slope
in the outer regions. This condition triggers the gravothermal
collapse of the central parts of the SIDM halo, which results in the
contraction of the core to form a new cusp, ultimately collapsing into
a black hole <cit.>[The gravothermal
  collapse <cit.> is a familiar process in globular
  clusters, where the inner regions have negative specific heat that
  is smaller than the positive specific heat in the outer regions. In the
  case of globular clusters, the collapse can be prevented by the
  formation of binary stars at the centre. In the case of a SIDM
  halo, since the interactions are purely elastic, the process is expected
  to continue until a black hole forms. The black hole efficiently accretes the
  inner core of the SIDM halo (e.g. <cit.>). This
  discussion refers strictly to elastic self-scattering. If collisions
  are inelastic, then the energy released needs to be taken into
  account and, in fact, it could prevent the gravothermal collapse;
  see <cit.>.].

For 0.1≲σ_T/m_χ≲10cm^2/g, cosmological
N-body simulations have shown that SIDM haloes today should have
cores of size of the order of the scale radius
<cit.>. At
the lower end of this range of cross-sections, SIDM cores are small,
∼0.2r_s <cit.>, making SIDM haloes only slightly
different from their CDM counterparts. This is why below this
cross-section, SIDM is essentially indistinguishable from CDM as a
theory of structure formation <cit.>. At the higher end of
the cross-section range, the core sizes are slightly larger than the
scale radius and approach the full thermalization of the core, with a
maximum size bounded by the radius at which the velocity dispersion
peaks (a case like this is shown in the bottom panels of
Fig. <ref>). Within the thermalised region, the
orbits of dark matter particles are isotropised by collisions, erasing
most of the memory of the assembly of the central regions
<cit.>. This makes haloes centrally rounder than
their CDM counterparts <cit.> and causes them to have local
velocity distributions that are close to Maxwellian
<cit.>. Since the onset of the gravothermal collapse
phase is expected to be ∼(250-400)t_ rel,s
<cit.>, the core phase of SIDM haloes in this range of
cross-sections is relatively long-lived.



§ HALO MERGERS AND THE EMERGENCE OF SUBHALOES


In the previous section we reviewed the structural properties of dark
matter haloes in CDM and in well-known alternatives such as WDM and
SIDM. In this section we focus on dark matter subhaloes. Since these
exhibit similar structural properties to haloes, modified by a few
relevant physical processes, we draw extensively on the results of the
previous section. Our goal now is to describe these processes and how
they affect the abundance and structure of subhaloes.



footnote-2
  footnoteReproduced from Carlo Giocoli et al. The substructure hierarchy in dark matter haloes . MNRAS (2010) 404 (1): 502–517, doi: 10.1111/j.1365-2966.2010.16311.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/404/1/502/3101607link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com

 footnoteReproduced from Volker Springel et al. The Aquarius Project: the subhaloes of galactic haloes. MNRAS (2008) 391 (4): 1685–1711, doi: 10.1111/j.1365-2966.2008.14066.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/391/4/1685/1747035link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com



 §.§ Halo mass assembly: smooth accretion vs mergers


Haloes grow by accreting dark matter, either through mergers with
smaller haloes or by accretion of diffuse, smooth material. The
importance of each of these channels depends on the shape of the
primordial power spectrum and on the smallest mass halo that can be
formed, both of which, in turn, depend on the nature of dark matter
particles. For instance, in WIMP CDM models, the minimum halo mass is
in the range (10^-12-10^-6) M_⊙
<cit.>, which is many orders of magnitude
below the resolution of current cosmological simulations. Thus, in
reality, the amount of smooth accretion measured in a simulation
consists of a combination of true unclustered dark matter and
unresolved dark matter haloes
(e.g. <cit.>). However, estimates based, for
example, on the excursion set formalism can be used to extend the
results of simulations into the unresolved regime.  In the resolved
regime, the analytical expectations are in good agreement with
simulations <cit.>. These calculations show that the amount
of smooth accretion onto present-day CDM Milky Way-size haloes is in
the range ∼ (10-20)% <cit.>.  Thus, in this
hierarchical structure formation model, haloes today are mainly
composed of the remnants of disrupted smaller haloes.  Of these, major
mergers (i.e., those with progenitor with mass ratios greater than 1:10)
contribute, on average, less than 20% of the final mass
<cit.>.

When does a halo become a subhalo? We mentioned earlier that the
boundary of a halo is not sharply defined, but rather chosen
approximately to separate the region within which the dark matter is
in dynamical equilibrium from an outer region where the dark matter is
still mostly infalling. In a similar way, the moment at which a halo
becomes a subhalo, i.e., when it crosses this transition region for
the first time, is somewhat arbitrary. For simplicity, it is common to
use the virial radius, e.g. r_200, as the boundary of the halo and
thus to define a subhalo as a halo that has crossed the virial radius
of a larger halo at some point in the past (see left panel of
Fig. <ref>). One could argue for a more physical
definition, based for instance in the relevance of the tidal forces
exerted by the dominant halo host in the local environment but, for
simplicity, and because of its common usage, we will use the former,
simple definition of a subhalo. We should remark that it is not
uncommon for subhaloes to leave the boundary of the main halo at some
point after first crossing <cit.>, and
thus the subhalo population today extends to radii a few times the
current virial radius (<cit.>; those systems beyond the
virial radius are commonly known as backsplash haloes
<cit.>).



 §.§ Evolution of subhaloes: initial conditions


Halo merger trees and merger rates.– N-body simulations have
been instrumental in determining the mass assembly history of
haloes. A particular powerful tool are halo merger trees (for an
overview of different algorithms to construct these trees see <cit.>). A halo at
the redshift of interest is regarded as the trunk of the tree and the
merger tree structure consists of a catalogue of progenitors, which
constitute the secondary branches that eventually merge onto the main
branch of the tree (see left panel of Fig. <ref>). Thus,
a merger tree contains information about the accretion times and
masses of the haloes that eventually become subhaloes. Both of these
properties, together with the corresponding position and velocity
vectors, represent the initial conditions for the subsequent dynamical
evolution of the subhalo.

An interesting statistical property that can be extracted from a
merger tree is the merger rate per halo, dN_m/dξ/dz
<cit.>, which gives the mean number of mergers, dN_m, per mass
ratio, dξ (relative to the main progenitor at the time of
accretion), per redshift interval, dz. This quantity has been found to have a
functional form that is nearly universal
<cit.>:


    dN_mdξ dz(M_0,ξ,z)=A(M_0/10^12  M_⊙)^αξ^β exp[(ξ/ξ̃)^γ](1+z)^η_z,

where M_0 is the mass of the descendant and the fitting parameters
(for the Millennium simulation) are given in Table 1 of
<cit.> (see also <cit.>). In fact, η_z is
very small, which implies a very weak redshift dependence; α>0
and β>0, implying that the merger rate is higher in more massive
haloes and for small mass ratios (the expected outcome in a
hierarchical model). We should remark that halo merger trees can also
be constructed from Monte Carlo realisations based on merger rates
computed using the extended Press-Schechter formalism,
e.g. <cit.>. This analytical approach, calibrated to
simulations <cit.>, is also widely used to model the
assembly of dark matter haloes, particularly in the context of
semi-analytic models of galaxy formation (e.g. <cit.>).

Distribution of accretion times and orbital properties.–
Eq. <ref> can be used to compute the average number
of haloes that become subhaloes of a host at a given redshift in a
certain mass range. Halo merger trees can be used to compute other
statistics of the subhalo population that are directly linked to their
subsequent evolution, specifically: (i) the distribution of accretion
(or infall) times, and (ii) the distribution of orbital properties at
the time of accretion.  In a time-independent spherical potential only
two variables are needed to specify the orbit of a tracer particle (plus
the orientation of the orbit). Although subhaloes orbiting a host halo
are far from this idealised case, it is nevertheless useful to
describe the initial orbital parameters in this way since this
provides a point of comparison with the simple spherical
potential. These two parameters can be chosen to be the radial, V_r,
and tangential, V_θ, velocities of the subhalo at the time of
infall <cit.>, typically expressed in terms of the virial
velocity of the host halo, V_200=√(GM_200/r_200). Another
common choice is to characterise the orbital properties at infall in
terms of a circular orbit of the same energy, E, and the same magnitude
of the angular momentum, j <cit.>. The initial orbit is
then defined by the circularity, η=j(E)/j_ circ(E),
and the infall radius, r_ 200/r_ circ(E), at the time of
accretion[Another set of parameters that can be used to
  define the orbit are the apocentre and pericentre. Different
  parametrizations can be transformed into one another since they are
  all related to the potential, ϕ(r).].

Fig. <ref> shows a sample of the orbital parameters
of the subhalo population at the time of infall, calculated from
high-resolution N-body simulations <cit.>. The bottom left
panel shows the distribution of subhalo infall times for Milky
Way-size haloes, M_200(z=0)=10^12 M_⊙. This plot only
includes subhaloes accreted after the formation redshift
of the halo, z_ HF, defined as the time at which the main
progenitor had half the mass of the final halo. Since the merger rate
is higher for larger descendant masses (Eq. <ref>), and
since the halo is growing from z=z_ HF until today, we expect
the distribution of infall times to decrease with redshift down to a
minimum at z_ HF, independently of the mass ratio. Naturally,
recently accreted haloes will be found mostly near the virial radius
of the host, while haloes accreted long ago will be mostly found in
the central regions (we will return to this point below). For
subhaloes in bound orbits, we would expect orbital velocities at
infall to be close to the virial velocity of the host halo, V_200
<cit.>.  The distributions of radial and tangential
velocities of infalling satellites for Milky Way-size haloes are shown
in the middle panels of Fig. <ref>. Although broad,
the distributions of orbital velocities do indeed have median values
around V_200. In fact, these distributions are not independent
since the total velocity of the subhalo orbit,
(V_r^2+V_θ^2)^1/2, is determined by the potential of the host
halo. This is why the ridge line of the bivariate distribution,
(V_r, V_θ), shown in the right panel of
Fig. <ref>, is close to circular (see also
<cit.>).



Reproduced from Lilian Jian et al. Orbital parameters of infalling satellite haloes in the hierarchical ΛCDM model. MNRAS (2015) 448 (2): 1674–1686, doi: 10.1093/mnras/stv053. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/448/2/1674/1053529link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com



 §.§ Dynamics of subhaloes


The material content of a halo consists of: (i) a smooth component
made up mostly of the debris of disrupted subhaloes but also of
material that was accreted in diffuse form; (ii) gravitationally
self-bound substructure - the subhaloes. As mentioned in
Section <ref>, the contribution of truly smooth
accretion to the total mass of a halo at z=0 is subdominant, and
even though most dark matter is accreted by (minor) mergers, only a
small fraction, ∼10%, remains today in bound subhaloes, at least
in the CDM model <cit.>[This estimate is based
  on an extrapolation over many orders of magnitude of the subhalo
  mass function determined in simulations down to the free-streaming
  mass of WIMP CDM particles. We discuss this in more detail
  below.]. Thus, most of the mass in a halo consists of the remnants
of the environmental processes responsible for stripping mass from
subhaloes after infall.

Fig. <ref> provides an illustration of the
richness of information contained in the phase-space structure (shown
here in its 2D radial projection) of haloes today that is relevant to
these environmental processes. In the left panel, the location in
phase space (at z=0) of the subhaloes of the Via Lactea II
simulation of a Milky Way-size halo <cit.> is shown
colour-coded according to the subhalo infall time. This figure reveals
the path of a subhalo as it orbits in the host halo: it first falls in
radially and 




footnote-2 
 footnoteReproduced from Miguel Rocha et al. Infall times for Milky Way satellites from their present-day kinematics. MNRAS (2012) 425 (1): 231–244, doi: 10.1111/j.1365-2966.2012.21432.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/425/1/231/998181link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com
 footnoteReproduced from Mark Vogelsberger and Simon D. M. White. Streams and caustics: the fine-grained structure of Λ cold dark matter haloes. MNRAS (2011) 413 (2): 1419–1438, doi: 10.1111/j.1365-2966.2011.18224.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/413/2/1419/1070092link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com

then loses energy as it is subjected to dynamical
friction and tidal forces in the host halo. The former causes the
subhalo to sink towards the centre while the latter gradually strip
mass from it, creating tidal streams. This picture can be appreciated
with clarity in the right-hand panels of the figure where the (2D)
phase-space structure of the dark matter particles is shown for a
different Milky Way-size halo simulation (Aquarius;
<cit.>). The colour in this case encodes the number of
caustics that a given particle traverses[Caustics represent
  folds in the fine-grained phase-space distribution function, which
  in CDM evolves according to the collisionless Boltzmann equation
  (Eq. <ref>). Before the formation of non-linear structures, CDM
  particles are distributed nearly uniformly in space with small
  density and velocity perturbations and very small thermal
  velocities. CDM particles thus occupy a thin, approximately three
  dimensional, sheet in 6D phase-space volume. Since CDM particles are
  collisionless and evolve according to Eq. <ref>, the
  fine-grained phase-space density is conserved during gravitational
  evolution (this was discussed earlier in the context of the maximum
  phase-space density in WDM models in <ref>), which
  implies that the original thin sheet can be stretched and folded but
  it cannot be broken. Caustics appear where folds occurs, and have
  very large spatial densities, limited only by primordial thermal
  motions
  (e.g. <cit.>).]. Since
caustics occur near orbital turning points, the number of caustics is
roughly proportional to the number of orbits each particle traverses.
The caustic count is thus an excellent way to highlight substructure
in the 2D phase-space structure seen in the right-hand panel of
Fig. <ref> since particles that are part, or
were part, of a subhalo have undergone more particle orbits in their
earlier host. This plot thus shows the richness of the substructure
present in haloes today. As we mentioned earlier, most of the matter
in a halo today has been accreted through mergers and consists of
material that was stripped from subhaloes.

Identifying substructure.– Several algorithms are in common use
to identify subhaloes in N- body simulations and define their
boundaries and properties. These subhalo finders are based on
different techniques. Here we will list only the most popular ones;
for a comprehensive comparison study see <cit.>.  A common
approach consists of finding local density maxima in the parent halo
density field and then associating adjacent particles with this peak
using a binding energy criterion; a subhalo is thus defined as the
collection of particles that is gravitationally self-bound, with the
density peak at its centre. Examples are: SUBFIND <cit.>;
Bound Density Maxima <cit.>; VOBOZ <cit.>;
Amiga Halo Finder <cit.>. A different approach is the “time
domain subhalo finder” which follows the time evolution of haloes and
tracks them when they become subhalos by identifying their bound
particles, as in the Hierachical Bound-Tracing or HBT
<cit.>.

An improvement over this 3D spatial approach can
be made by including information on the particle velocities, which has
the advantage that subhaloes that are in close spatial proximity with
one another can be more easily disentangled. In this case a density
criterion is not enough, but the relative velocity between merging
subhaloes is a telltale sign of a merger. Modifications of these
algorithms can be used to identify the tidal streams that are the
remnants of the tidal stripping process (see below) and which are not
localised in real space, but have clear signatures in phase
space. Examples of phase-space finders are the Hierarchical Structure
Finder <cit.> and ROCKSTAR <cit.>.

Current subhalo finders are able to identify subhaloes down to
20-100 particles and different algorithms roughly agree with one
another on their location and main properties <cit.>. Below
we review the main physical processes that affect the evolution and
inner structure of subhaloes along their orbits, as well as relevant
lessons learned from N-body simulations.

Tidal stripping.– Once a halo reaches the outer boundary of the
host halo into which it will merge, tidal forces begin to act,
suppressing the accretion of matter into the merging halo and
stripping mass from its outer layers in a process known as tidal
stripping. Since the enclosed overdensity of a subhalo depends on its
position within the host, the virial radius is no longer a meaningful
concept. A more relevant scale is the tidal radius, r_t,
defined as the radius at which the differential tidal force of the
host halo is equal to the gravitational force due to the mass of the
subhalo, or equivalently, as the radius within which the enclosed mean
mass density of the satellite is comparable to the mean mass density
of the main halo interior to the distance, R, to the satellite. The
expectation is that the matter beyond the tidal radius will be removed from
the subhalo, reducing its mass as it orbits around the host. For a circular
orbit and assuming that the subhalo mass, m_ sub(<r_t), is much
smaller than the enclosed mass of the host, M(<R), and that
r_t≪ R, the tidal radius is given by
<cit.>[This equation ignores the effects of the
  centrifugal force on the satellite as it orbits around the
  host. Including this effect (assuming circular orbits) modifies
  Eq. <ref> slightly by changing the factor of 2 in
  the denominator to 3 <cit.>.]:

    r_t=R[m_ sub(<r_t)/(2- dlnM/ dlnr)M(<R)]^1/3.


For non-circular orbits the situation is more complex (in fact in the
most general cases, the tidal radius can be ill-defined; see
e.g. <cit.>), but the principle behind the concept of tidal
radius remains valid: the relevant physical quantity to determine the
boundary of the subhalo is the relative strength of the gravitational
attraction of the subhalo and the tidal forces of the host. In this
way, the tidal radius is commonly used to model tidal stripping in a
variety of collisionless systems, not only subhaloes, but also, for
example, globular clusters. In particular, for a slowly varying
tidal field (i.e., in the adiabatic approximation), the mass loss due
to stripping may be modelled as:
<cit.>:

    dm_ sub=α_t m_ sub(>r_t)dt/t_ orb(R),

where t_ orb(R) is the instantaneous orbital period at the
radius of the subhalo, and α_t is a tuning parameter, which
encapsulates departures from this simple approximation (e.g., the
details of subhalo structure); α_t is typically calibrated from
simulations but the values used vary significantly in the literature,
which is a major limitation. 
Eq. <ref> assumes that the relevant timescale for mass
loss is the orbital period of the subhalo, which is justified by
noting that the energy scale for the tidally stripped material is
given by the change in the potential of the host across the body of
the satellite ∼ r_tdϕ_ host(R)/dR
<cit.>. The left panel of Fig. <ref>
shows an example of tidal stripping in this slow (adiabatic) mode for
a subhalo in a circular orbit around a static host potential (from
<cit.>). As pointed out by <cit.>, a
combination of the ill-defined tidal radius and uncertainty in the
parameter α_t, makes the modelling of this adiabatic case quite
complicated, with the end result that models based in
Eq. <ref> do not, in general, match simulation results
accurately. In any case, it is rare for subhalo orbits to be nearly
circular; for realistic orbits, most of the tidal mass loss
happens near pericentre.

Tidal shock heating.– While tidal stripping
(Eq. <ref>) refers to the gradual loss of loosely bound
material from a subhalo due to a slowly-varying external potential, a
rapid (impulsive) variation in the potential causes a transfer of the
satellite's orbital energy to the internal energy of its
particles. These tidal shocks are most important at pericentre where
the impulsive condition is best satisfied. The redistribution of
internal energy produced by the shock alters the inner structure of
the subhalo and can unbind some of its particles
<cit.>.  This process is well
described by the “impulsive approximation” (see
<cit.> for the case of globular clusters) in
which tidal forces are assumed to act during a much shorter time than
the dynamical timescale of the satellite (see
<cit.>). The approximation gives the specific
energy change suffered by particles in the subhalo due to a tidal
shock as[It is possible to evaluate Eq. <ref>
  for a given fixed potential and a given subhalo orbit (see e.g
  Eq. 20 of <cit.> for an NFW halo).]:

    ( Δ E )_ i,tid=(Δ v)^2_ i,tid≈|∫_ orbita⃗_ i,tid(t) dt|^2,

where a⃗_ i,tid is the acceleration experienced by a
particle in the subhalo and the integral is performed along the orbit
of the subhalo. If ( Δ E )_ i,tid>E_i,b, where E_i,b
is the binding energy of the particle, we may assume that the particle
will become unbound. The mass fraction
of particles that satisfies the inequality is then assumed to be
removed instantaneously from the subhalo. The impulsive
approximation accurately captures the results of simulations for
radial orbits. An example is shown in the right-hand panel of
Fig. <ref> (taken from <cit.>).

Although the energy injection in Eq. <ref> from tidal
shocks might not be enough to unbind particles, it can still
affect the inner structure of a subhalo. As a result of the shock, the
orbits of dark matter particles in the centre expand, reducing the
inner density, although this process is not strong enough to create a
central core <cit.>. The resulting density
profile is, in fact, still well described by an NFW profile in the
inner regions, albeit with a higher concentration, while the outer
regions are considerably steeper than NFW due to stripping. For
instance Ref. <cit.>, using idealised simulations in a
static external halo potential, found a profile of the form
ρ∝ r^γ exp(-r/r_b), with a central slope,
γ∼1, and a cutoff radius due to tidal shocks, r_b (see
bottom right panel of Fig. <ref>).

For a general subhalo orbit, a combination of the adiabatic
(Eq. <ref>) and impulsive (Eq. <ref>)
approximations provides a good estimate of the amount of stripped
mass; the former is valid particularly near the apocentre of nearly
circular orbits while the latter is more appropriate near the
pericentre of eccentric orbits. The impulsive approximation reproduces
the results of simulations quite accurately for radial orbits, but the
adiabatic approximation is not very adequate for the reasons discussed
above (see also Section 4.3 of <cit.>).



footnote-3 
 footnoteReproduced from Frank C van den Bosch and Go Ogiya. Dark matter substructure in numerical simulations: a tale of discreteness noise, runaway instabilities, and artificial disruption. MNRAS (2018) 475 (3): 4066–4087, doi: 10.1093/mnras/sty084. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/475/3/4066/4797185link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com
 footnoteReproduced from Frank C van den Bosch et al. Disruption of dark matter substructure: fact or fiction? MNRAS (2018) 474 (3): 3043–3066, doi: 10.1093/mnras/stx2956. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/474/3/3043/4638541link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com
 footnoteAAS. Reproduced with permission. For the original article, please visit the following https://iopscience.iop.org/article/10.1086/420840link.

Subhalo-subhalo encounters.– Tidal heating can also be caused
by impulsive encounters with other subhaloes, which can add up to
produce a net effect on the subhalo inner structure and mass loss (a
process called galaxy harassment in the context of satellite
galaxies; <cit.>). A similar impulsive approach to the one
above can be used to estimate the strength of this form of tidal
heating, but the calculation is more complicated because, among other
things, the distribution of subhaloes in the host and the encounter
rate need to be modeled. A recent study finds that tidal shocking
from encounters is subdominant (by a factor of several) compared to
shocking during pericentric passages <cit.>.

Dynamical friction.– When an object of mass, M_s, moves
through an ambient medium of collisionless particles of mass
m≪ M_s, the object experiences a drag force known as dynamical
friction. This force may be thought of as the gravitational pull
exerted by a local enhancement in the ambient density formed behind
the moving object (a trailing wake) as the object gravitationally
focuses the surrounding particles (see left panels of
Fig. <ref>).



footnote-2
 footnoteReproduced from Go Ogiya and Andreas Burkert. Dynamical friction and scratches of orbiting satellite galaxies on host systems. MNRAS (2016) 457 (2): 2164–2172, doi: 10.1093/mnras/stw091. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/457/2/2164/970692link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com

 footnoteReproduced from Michael Boylan-Kolchin et al. Dynamical friction and galaxy merging time-scales . MNRAS (2008) 383 (1): 93–101, doi: 10.1111/j.1365-2966.2007.12530.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/383/1/93/1067887link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com


The net result of
dynamical friction is a transfer of orbital angular momentum and
energy from the moving object into the surrounding medium. The process
can be analysed as a series of uncorrelated sequential encounters
between the object of mass, M_s, and velocity, v⃗_s, and
particles randomly taken from the ambient medium with velocity
distribution, f_m(u⃗). These interactions occur on a timescale
much shorter than the variations in the velocity, v⃗_s, of the
object. If the ambient medium is assumed to have a homogeneous
density, ρ_m, then the changes to the velocity of the object
perpendicular to its motion average to zero by symmetry, while the
velocity changes parallel to the direction of motion are given by the
well-known Chandrasekhar dynamical friction formula
<cit.>, which, for the drag force, F⃗_ df,
takes the form: 

    F⃗_ df=M_sdv⃗_sdt   =   -16π^2M_s^2m lnΛ[∫_0^|v⃗_s| f_m(|u⃗|)u^2d|u⃗|]v⃗_s/|v⃗_s|
       =   -4π(GM_s/v_s)^2 lnΛ ρ(<|v⃗_s|)v⃗_s/|v⃗_s|,

where ρ(<|v⃗_s|) is the density of ambient particles
with speed less than |v⃗_s| and
lnΛ≡ ln[1+(b_ max/b_90)^2]
is the Coulomb logarithm, with b_90=G(M_s+m)/v_∞^2, 
v_∞ the initial relative velocity of an individual encounter,
and b_ max the maximum impact parameter
(b_ max≫ b_90). As a consequence of dynamical friction, the
orbit of the object decays in time sinking towards the centre of the
host halo. For circular orbits in a spherical singular isothermal host
halo (implying a Maxwellian velocity distribution, f_m) of mass,
M_h, the timescale for the orbit to decay to zero (i.e., the
dynamical friction time) is approximately given by <cit.>: 

    t_ df/t_ dyn≈1.17/
        ln(M_h/M_s)(M_h/M_s),

where t_ dyn=r_ vir/V_c is the dynamical timescale at the
virial radius of the halo, r_h, with V_c the circular velocity of
the host halo, which is independent of radius for a singular
isothermal sphere.  For more general orbits, Eq. <ref> requires
a correction that scales with the circularity, η, as
t_ df∝η^γ_η, where γ_η∼0.53
<cit.>, implying that more eccentric orbits decay more
rapidly (see right panel of Fig. <ref>).

Eq. <ref> is derived under the following assumptions:
(i) both the satellite and the particles that make up the ambient
medium can be treated as point masses; (ii) the self-gravity of the
ambient medium can be ignored, and (iii) the distribution of ambient
medium particles is infinite, homogeneous and isotropic. None of these
assumptions is strictly valid in realistic situations. Nevertheless,
Chandrasekhar's formula provides a reasonable description of dynamical
friction, particularly when modifications are included to account for
the density profile of the subhaloes and their orbits; in practice,
this can be done by regarding the Coulomb logarithm as a free
parameter that depends on these properties. For example, <cit.>
carried out a series of idealised N-body simulations of a subhalo
infalling into a host, both described by a Hernquist density
profile,[The Hernquist halo profile <cit.> has
  the same asymptotic behaviour at the centre as the NFW halo and has
  the advantage that the velocity distribution function in the
  isotropic case has an analytic form (see Eq. <ref>), which
  makes it particularly simple to set up initial conditions for
  simulating haloes in dynamical equilibrium.] and found the following fitting function 
  to the dynamical friction timescale
  (a few examples of the orbital evolution in this study are shown in the right
panel of Fig. <ref>):

    t_ df/t_ dyn=A(M_h/M_s)^b/
        ln(1+M_h/M_s) exp[cη(E)][r_c(E)/r_200]^d,

where E is the initial orbital energy of the satellite (we recall
that r_c(E) is the radius of a circular orbit of the same energy,
E), and the fitting parameters are of order one[The values
  for these parameters reported in <cit.> are: A=0.216,
  b=1.3, c=1.9, and d=1.0, but we point out that in this study
  both the halo and the subhalo were modelled as Hernquist
  profiles.]. Eq. <ref> was found to be valid over a wide
range of orbital parameters; the most relevant restriction is
0.025≤ M_s/M_h≤0.3. For smaller mass ratios, the dynamical
friction timescale becomes much larger than the age of the Universe,
while for larger mass ratios, the relevant timescale is just the
dynamical or free-fall time[Eq. <ref> was only
  explored for values of the circularity in the range
  0.3≤η≤1.0 and for 0.65≤ r_c/r_200≤1.0; the
  lower limits were imposed in order to avoid radial orbits that would
  take the subhalo so close to the centre of the halo in the first
  orbit that the tidal effects of the galaxy cannot be ignored.  So
  far we have not discussed baryonic effects, but it is worth
  mentioning them here since Eq. <ref> was not investigated
  outside this range and might not be valid there even in the absence
  of a central galaxy.].



 §.§ The abundance, spatial distribution and internal structure
  of dark matter subhaloes


The abundance, spatial distribution within the host halo and 
internal structure of subhaloes are determined by the combined effects
of the initial conditions at the time of accretion, which depend on
cosmology, and the dynamical processes described in the previous
section. These properties are best derived in full generality using
cosmological N-body simulations but analytic models can provide
valuable physical insights <cit.>.  In this section we present
some of the key structural properties of the subhalo population, as 
revealed by simulations. Naturally, these properties are closely
related to those of isolated haloes (discussed in
Section <ref>) with a few relevant modifications.



AAS. Reproduced with permission. For the original article, please visit the following https://iopscience.iop.org/article/10.1088/0004-637X/740/2/102link.

The subhalo mass function.– As in the case of isolated haloes,
the total CDM subhalo mass function (measured within the virial
radius) is remarkably close to universal and, in fact, has a similar
low-mass slope as the halo mass function
<cit.>:
dn_ sub/dm_ sub∝ m_ sub^α, where
α ∼ -1.9 (see Eq. <ref>). This similarity to the
halo mass function is partly due to the fact that most subhaloes
identified at a given time were accreted relatively recently and thus
tidal effects have not had time to act; see
Fig. <ref>. The normalisation of the subhalo mass
function depends on the mass of the host halo, with more massive
haloes having, on average, larger subhalo populations
<cit.>. This reflects the earlier assembly of
low-mass haloes, which allows tidal effects more time to act and
disrupt subhaloes. For similar reasons other properties of the host
halo can have second-order effects on the amplitude of the subhalo
mass function, e.g., at fixed mass, more concentrated haloes (which
assemble earlier) have fewer subhaloes.

When the subhalo mass function is scaled to the host halo
mass it becomes nearly universal across halo masses, with a functional
form that is well fitted by:

    N(>μ≡ m_ sub/M_h)=(μ/μ_1)^1+α exp[-(μ/μ_ cut)^b],

<cit.> where the exponential cutoff
accounts for the increasing rarity of subhaloes of mass close to that
of the host halo mass[The fitting parameters in
  Eq. <ref> in the case of the Millennium simulations may be
  found in <cit.>, where a redshift dependence is also
  provided.]. The parameter, μ_1, is the typical mass
fraction of the most massive subhalo (relative to the host halo mass),
which, for a Milky Way-size halo, is of order 0.01, but with a large
spread <cit.>. The universality of the subhalo mass function
is, however, not perfect; the remaining dependence on host halo mass
can be captured by allowing a relatively weak scaling of the
normalisation parameter, μ_1. This dependence is
amplified, and perhaps better expressed, if the subhalo velocity
  function is used instead, that is, if the abundance is given in
terms of the maximum circular velocity instead of the mass (see left
panel of Fig. <ref>). In this case, an accurate
approximation is given by <cit.>:

    N(>x≡ V_ sub/V_h)∝ V_h^1/2x^-3,      x<0.7,

which implies, for instance, that a cluster-sized halo
(V_h∼1000 km/s) has ∼2.2 times more substructure of a given
velocity ratio than a Milky Way-size halo (V_h∼200 km/s). Notice
that this difference is considerably weaker if the mass ratio is used
since, in this case, the abundance scales as a power law of exponent
∼-0.9 rather than ∼-3.

The fact that the power-law exponent of the subhalo mass function at
low masses, α, is greater than -2 is important; a steeper
slope would imply that the total mass in substructures diverges when
extrapolated to arbitrarily low masses. For a given particle dark
matter model we know, of course, that this extrapolation cannot be
continued beyond the truncation mass below which the properties of the
dark matter particle prevents the formation of smaller structures (due
to the suppression mechanisms mentioned in
Section <ref>). In the case of CDM WIMPs, the
extrapolation of the subhalo mass function down to the Earth's mass
(10^-6 M_⊙) implies that the fraction of mass contained in
unresolved subhaloes is ∼4.5%, in contrast to the
∼13% mass fraction found in the highest resolution simulation
(with a particle mass of 2× 10^4 M_⊙) to date of a Milky
Way-size halo <cit.>. As mentioned earlier, most of the
mass in haloes is not, in fact, in the form of self-bound subhaloes,
but in the remnants of the tidal stripping process accumulated over
the entire history of the halo.

The radial distribution of subhaloes.– The spatial distribution
of a subhalo population reflects the balance between accretion of new
subhaloes and tidal disruption of older ones. This distribution has
been studied extensively in N-body simulations
<cit.>
and the picture that emerges is that the radial distribution of
subhaloes is significantly less centrally concentrated than the dark
matter distribution (i.e., the smooth halo), and is relatively
independent of the host halo mass (see right panel of
Fig. <ref>). Most remarkably, when subhaloes are
selected according to mass (rather than maximum circular velocity) and
the distribution is normalized to the mean number density of subhaloes
of a given mass within the virial radius, there appears to be no trend
in the shape of the number density profile with subhalo mass
<cit.>.  A recent analysis using
the HBT finder, however, has shown that the most massive subhaloes are
actually more concentrated in the central regions than lower mass
subhaloes <cit.>, which seemingly reflects the resilience of
very massive subhaloes to tidal stripping despite suffering from
substantial dynamical friction. The near universality of the radial
distribution of subhaloes is then the result of a convolution of the
distribution of subhaloes before infall (sometimes called the unevolved radial distribution of subhaloes), which is nearly scale
free, with a tidal stripping process that is also nearly scale free,
except at the massive end (see <cit.> for an analytical model
of the subhalo distribution).








It is interesting that the radial distribution of subhaloes with
maximum circular velociy >V_ sub is steeper than that of
subhaloes with mass >m_ sub (see e.g. <cit.>), since
the latter is more heavily influenced by tidal stripping.






The ratio between the average, mass-selected subhalo radial
distribution and the average NFW mass density profile of their host
haloes (both normalised to the virial radius as defined in
<cit.> and shown in the right-hand panel of
Fig. <ref>), is approximated quite accurately by the
following functional form <cit.>:

    ϕ(x≡ r/r_ vir)=dN/dx^3|_ sub/dN/dx^3|_ NFW=4x^4/(1+x)^2


The inner structure of dark matter subhaloes.– Cosmological
N-body simulations have shown that the density profiles of subhaloes
retain the near universal properties of isolated field haloes but with
modifications that reflect the tidal effects discussed in
Section <ref>. These modifications are consistent with
expectations of analytical estimates and controlled simulations.  In
particular, for CDM, the subhalo radial density profile exhibits the
same central cusp as an isolated halo in equilibrium (left panel of
Fig. <ref>), while the outer regions show a steep
truncation at a radius approximately equal to the tidal radius given
in Eq. <ref> (see Fig. 15 of <cit.>). We
should remark that, as has been found for field haloes
<cit.>, a better fit to the density profile of subhaloes
is given by the 3-parameter Einasto profile
<cit.>[Although the introduction of a third
  parameter will obviously improve the quality of the fit, the Einasto
  profile is, in fact, a slightly better fit to simulations than the
  2-parameter NFW profile even after one of the parameters
  (α_E) is fixed to an appropriate value. For instance, fixing
  α_E∼0.16 gives a better fit than NFW to haloes across a
  range of halo masses <cit.>.]:

    ρ_E(x_E≡
    r/r_-2)=ρ_-2e^-2(x_E^α_E-1)/α_E,

where α_E is a shape parameter and ρ_-2 and r_-2 are
the density and radius at which the logarithmic slope of the density
profile is equal to -2.  The Einasto and NFW profiles are quite
similar, and both give good fits to the subhalo profiles in the range
0.01<x_E<100 if α_E∼0.22
<cit.>. Although for isolated haloes the
parameters α_E and r_-2 can be related to the virial mass
of the halo, M_200 <cit.>, in a similar way as the
halo (NFW) concentration is connected to the virial mass, the
situation is less clear for subhaloes <cit.>, and the
spread of the parameters across subhalo masses is large. Thus, for its
simplicity, the NFW profile remains a reasonable approximation to the
structure of both haloes and subhaloes.

Since for subhaloes the virial radius no longer has a proper meaning
as the “boundary” of the object, the concentration parameter,
defined as c=r_200/r_s commonly used to characterize NFW haloes,
is no longer appropriate. Instead, it is convenient to define the
concentration of a subhalo in a way that is independent of its size.
One such measure of concentration is the characteristic overdensity,
δ_V, defined as the mean density within the radius,
r_ max, where the circular velocity peaks, at a value of
V_ max, relative to the critical density
<cit.>:

    δ_V=ρ(<r_ max)/ρ_ crit=2(V_ max/Hr_ max)^2,

where H is the Hubble constant. Eq. <ref> can be related to
the standard scale density of the NFW profile (δ_c in
Eq. <ref>), and thus to the NFW concentration, in a
straightforward way <cit.>:

    δ_V=(c/2.163)^3K_c(2.163)/K_c(c)Δ,

where K_c was defined just after Eq. <ref>. We note that
for the NFW profile, r_ max/r_s=2.163.



 Reproduced from Volker Springel et al. The Aquarius Project: the subhaloes of galactic haloes. MNRAS (2008) 391 (4): 1685–1711, doi: 10.1111/j.1365-2966.2008.14066.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/391/4/1685/1747035link. This figure is not included under the CC-BY license of this publication. For permissions, please email: journals.permissions@oup.com

Since the concentration of haloes (and subhaloes) is tightly
correlated with their mass (see Section <ref>),
Eq. <ref> implies a tight correlation between V_ max
and r_ max, which indeed has been found and characterised in
simulations (see right panel of Fig. <ref>).  For
instance in the case of the Aquarius-A Milky Way-size halo, the
following fitting functions (to the mean relations in subhaloes)
provide a direct connection in terms of the subhalo mass[We
  note that there is a typo in the caption of Fig. 28 in
  <cit.>, which gives the fitting function for δ_V
  and m_ sub
  (5.8×10^8 M_⊙→5.8×10^4 M_⊙).]:

    V_ max   =   10  km/s(m_ sub/3.37×10^7 M_⊙)^0.29
    δ_V(z=0)/2=(V_ max/H_0r_ max)^2   =   2.9×10^4(m_ sub/10^8 M_⊙)^-0.18.

The fitting function for the r_ max - V_ max relation
implied by Eq. <ref> is shown as a dotted line in the upper
right panel of Fig. <ref>, while the corresponding
relation for isolated haloes in the Millennium I simulation is shown
as a solid line <cit.>[Obtained by taking the
  power-law concentration-mass relation in <cit.> (their
  Eq. 4), and using Eqs. <ref> and <ref>.]. The most
relevant result when comparing haloes and subhaloes in the
r_ max-V_ max plane is that both share the same relation,
but subhaloes have systematically higher concentrations at a given
V_ max <cit.>: in Fig. <ref> the
dotted line is a factor of 0.62 lower than the solid line,
i.e., subhaloes have on average r_ max values that are smaller
than haloes of the same V_ max by this factor. Equivalently,
the characteristic overdensity, δ_V, in subhaloes is roughly a
factor of (1/0.62)^2∼2.6 larger for subhaloes than for haloes of
the same V_ max (lower right panel of
Fig. <ref>), which roughly corresponds to a 30%
increase in the NFW concentration. This relative increase in
concentration is larger for subhaloes nearer the centre of the host,
as expected from the nature of the tidal forces experienced by the
subhaloes as described in Section <ref>: while tidal
stripping naturally reduces V_ max, it reduces r_ max
even further <cit.>[For a clear illustration of
  the evolutionary track of suhaloes in the r_ max-V_ max
  plane due to tidal stripping, see Fig. 8 of <cit.>.],
making the subhalo effectively more concentrated; the stronger the
mass loss, the stronger the effect, and hence the trend with
halocentric distance.

It is thus possible to model the inner density profile of the subhalo
population by assuming a model for the concentration-mass relation of
field haloes and making a simple correction to the subhalo
concentration depending on the location of the subhalo.  More
exhaustive studies of subhalo concentration exist that provide fitting
functions across a wide range of subhalo masses, host halo masses, and
distance to the halo centre (e.g. <cit.> for the case of
Milky Way-size haloes).

The shapes and internal kinematics of subhaloes.– The impact of
tidal forces in the structure of subhaloes is reflected also in their
shapes. Although tides tend to elongate objects, these distortions are
short-lived features accentuated during pericentric passages. Once the
tidal streams cease to be bound to the subhalo, simulations have shown
that the bound material remains in an equilibirum configuration that
is, in fact, more spherical than it was at the time of infall; the
stronger the tidal effects, the more spherical the subhalo becomes
<cit.>. Although these differences are significant for the
fraction of the subhalo population whose orbits are strongly
influenced by the tides of the host, the subhalo population as a whole
is only slightly affected and exhibits a small systematic shift
towards less aspherical shapes compared to field haloes
<cit.>. This is because the global subhalo population is
dominated in number by subhaloes near the virial radius of the host,
which have only recently fallen in.

When tidal effects are strong, the internal kinematics of subhaloes
are also substantially altered. In particular, the velocity anisotropy
of the dark matter particles becomes increasingly tangential
(β<0) from the subhalo centre outwards <cit.>, in
contrast to field haloes that are radially anisotropic at larger
radii. This is the result of the preferential stripping by tides (when
the subhalo is near pericentre) of subhalo particles with radial
orbits. On the other hand, the pseudo-phase space density,
Q_ sub, of subhaloes in equilibrium retains the universal
power-law behaviour of CDM field haloes, but with a slightly shallower
slope, Q_ sub∼ r^-1.6 <cit.>, compared to
∼ r^-1.9 for field haloes.



 §.§ The impact of the nature of the dark matter


Subhalo abundance.– By far the main difference in the subhalo
populations predicted in models with different kinds of dark matter is
the abundance of low-mass subhaloes. 



In particular, as we
discussed in Section <ref>, models in which the
primordial power spectrum of density perturbations has a cutoff at
relatively low k (such as WDM and interacting dark matter) have a corresponding cutoff in the mass
function of haloes and subhaloes. These models predict far fewer haloes
and subhaloes than CDM, and this offers the best prospect for
distinguishing between them and perhaps constraining the properties of
the particles themselves (such as the WDM particle mass).

A cutoff in the mass function breaks the universal behaviour of the
halo and subhalo mass functions at low masses in a way that also
depends on the nature of the dark matter particle. For example, the
self-similarity of the abundance of CDM subhaloes as a function of
relative mass, exhibited in Eqs. <ref> and <ref>, is
broken <cit.> because the cutoff scale expressed in terms of
the ratio, μ=m_ sub/M_h, occurs at larger values of μ for
smaller values of M_h. The radial distribution of subhaloes in WDM
models is quite similar to that in the CDM case, with only minor
differences explained by the enhanced tidal stripping of low-mass WDM
haloes resulting from their lower concentrations
<cit.>.

In many SIDM models the subhalo mass function remains largely
unchanged as long as the interaction cross-section,
σ_T/m_χ<10 cm^2/g
<cit.>. For higher values,
collisions between dark matter particles within subhaloes and in the
host are frequent enough to unbind material from the halo. This form
of subhalo evaporation is energetically efficient because the
energy transfer is determined by the relative velocity of the
colliding particles, which is of the order of the orbital velocity.
In this case, the mass loss in subhaloes is enhanced and the subhalo
abundance is depleted relative to the CDM case, particularly in the
central regions <cit.>.

Inner structure of subhaloes.– The inner structure of WDM and
SIDM subhaloes is rather similar to that of field haloes (see
Fig. <ref>), and the outer structure is altered
by tidal effects in a very similar way as in CDM. The main difference
is an enhancement in the concentration of subhaloes relative to their
field counterparts in WDM models <cit.> due to the increased
efficiency of tidal stripping of WDM subhaloes, which are less
concentrated than in CDM at the time of infall.  Tidal stripping also
plays a greater role in SIDM subhaloes <cit.>. Most
importantly, it can trigger a gravothermal catastrophe and this can
give rise to segregation according to particular orbits, with cuspy
profiles for subhaloes which have experienced substantial tidal mass
loss and central cores for those where tidal effects have been minimal
<cit.>.



§ OUTLOOK


It is fair to say that the evolution of the phase-space distribution
of classical, non-relativistic, collisionless dark matter (CDM) down
to galactic-scale haloes and subhaloes is now essentially a solved problem,
largely through the application of N-body simulations over the past
40 years[By galactic-scale haloes and subhaloes, we mean
  self-bound dark matter structures that can potentially host a
  galaxy, that is haloes of mass above ∼10^8 M_⊙, in which
  gas can cool by atomic processes
  (e.g. <cit.>).]. This strong statement carries
a couple of major caveats, which define today's frontier in N-body
simulations of cosmological structure formation.

Firstly, the statement above can still hold if, instead of CDM, most
of the dark matter consists of other types of particles, such as WDM
and SIDM[This is true only for elastic SIDM, and for cross
  sections that do not exceed the gravothermal collapse threshold,
  σ_T/m_χ∼10 cm^2/g, for dwarf-size haloes (see the
  last paragraphs of Section <ref>). Although the regime
  of gravothermal collapse has been known for a couple of decades
  <cit.>, a comprehensive analysis of this
  regime has yet to be carried out (see
  <cit.> for recent
  developments in this interesting regime).], for which N-body
simulations with appropriate modifications have been applied at a
similar level of detail as in CDM; in this review we have discussed
the most important changes in the dark matter phase-space structure
that occur in these alternative models. Nevertheless, there are still
dark matter models that remain unexplored, or only partially explored
with N-body simulations, e.g. hidden dark sector models with DAOs
<cit.> and inelastic
SIDM[There is a class of inelastic SIDM models in which the
  dark matter can have ground and excited states
  (e.g. <cit.>), and in which scattering between the
  excited and ground states can result in energy injection at the
  centre of dark matter haloes thus altering their structure. Only
  until very recently have these models began to be explored with
  simulations <cit.>.].  Secondly, and
crucially, the statement above does not take account of the interplay
between baryons and dark matter, which are dynamically coupled through
gravity. Several mechanisms that can radically modify the
dark-matter-only predictions of N-body simulations and which are, of
course, crucial for a complete theory of structure formation and its
connection to reality, have been studied extensively for several
decades. We briefly summarise these in Subsection <ref>
below.

Finally, we make two remarks concerning the limited
resolution of current N-simulations: (i) there have been recent
claims that subhaloes can be artificially disrupted in cosmological
simulations due to discreteness effects and inadequate force
resolution <cit.>; if correct, these effects
could alter some of the current results on the abundance and structure
of subhaloes, particularly at low masses; (ii) as we have seen, the
best current N-body simulations only resolve haloes of mass greater
that ∼ 10^5 M_⊙, many orders of magnitude larger than the
cutoff mass in the linear density power spectrum for CDM. Yet, if the
dark matter is made of Majorana particles, these, so far unresolved, 
haloes could be crucial for predicting the properties of their
annihilation radiation and thus for elucidating the nature of dark
matter. The first attempts at understanding the properties of haloes
down to the cutoff in the CDM primordial power spectrum have been made
<cit.> but new techniques will
be required to tackle this problem in full generality.



 §.§ The impact of baryonic physics on dark matter structure


In the linear regime, the (gravitational) impact of baryons (and
electrons and photons) in the dark matter distribution, of which
baryonic acoustic oscillations is perhaps the best known outcome
<cit.>, is fairly well understood. In the
non-linear regime, on the other hand, the complexity of baryonic
physics is much greater and the list of relevant processes is
extensive: gasdynamics, radiative processes, star formation and
evolution, supermassive black hole formation and evolution, etc. Here
we focus on some of the most important mechanisms that modify the
predictions for the abundance and structure of CDM haloes from
N-body simulations.

Condensation of baryons into haloes: adiabatic gas cooling and
  mergers.- In the classical theory of galaxy formation, gas
initially follows dark matter; as haloes collapse and virialize, the
associated gas heats up by shocks and adiabatic compression to the
virial temperature of the halo
<cit.>[Large relative velocities between
  gas and dark matter inherited from the photon-baryon coupling
  before recombination can impede the growth of gravitational
  perturbations and also stop gas from accreting into the first haloes
  <cit.>. This process, however, is only thought to
  be relevant for the formation of the first stars.].  Subsequently,
the gas can radiatively cool and condense towards the centre of the
halo if the cooling time is shorter than the free-fall time. The halo
mass threshold for effective cooling depends on the density,
temperature and metallicity of the gas; cooling is quite efficient in
low-mass haloes down to the atomic cooling limit (virial temperatures
∼10^4 K, corresponding to halo masses ∼10^8 M_⊙) below
which cooling becomes highly inefficient. At higher masses
(∼10^13 M_⊙ for gas with solar metallicity,
e.g. <cit.>), cooling is also suppressed because the cooling
time exceeds the free-fall time, limiting the condensation of baryons,
a process that can be exacerbated when the gas is heated by energy
input from Active Galactic Nuclei (AGN) <cit.>. A
hot, quasi-hydrostatic corona forms from which gas can subsequently
cool at the centre. Additional gas may be brought in by galaxy
mergers.  Regardless of the condensation mode, the assembly of the
central galaxy ultimately results in an enhancement of the central
gravitational potential, compared to the situation where the galaxy is
absent. The dark matter distribution reacts dynamically, becoming more
concentrated, a process first modelled assuming an adiabatic response
leading to the contraction of the halo
<cit.>. Even though the assembly of baryonic
matter by mergers is not, in general, adiabatic, the simple adiabatic
model remains a reasonable approximation <cit.>. In the
absence of heating processes, the general expectation is thus that
haloes should be cuspier than the NFW profile in the central regions,
as indeed is seen in cosmological hydrodynamics simulations
<cit.>.

Energy injection into haloes: UV background photoheating.- The
hydrogen emerging from recombination is, of course, neutral. However,
the UV radiation produced by stars in the first generations of
galaxies reionises this gas and heats it up, suppressing gas cooling
into low-mass haloes and subsequent star formation
<cit.>. This heating mechanism moves the
minimum scale for galaxy formation from the atomic cooling limit to
larger halo masses of order 10^9 M_⊙ today[This mass
  threshold is smaller at higher redshifts, see e.g. Fig. 3 of
  <cit.>.]
<cit.>.
This baryonic process is important also because, in conjunction with
the expulsion of gas from haloes by supernova feedback (see below) at
high redshift, it reduces the overall baryonic content, and thus, the
total mass content of low-mass haloes; this reduces the growth rate and
final masses of these haloes compared to their counterparts in
simulations without baryons <cit.>.

Energy injection into haloes: supernova and AGN feedback.- When
massive stars explode as supernovae in the final stages of their
evolution they release vast amounts of energy, a fraction of which
may couple effectively to the surrounding interstellar medium (ISM),
heating it and pushing it in a violent blowout. The combined impulsive
removal of baryonic outflows from several supernovae creates a
collective effect in the host galaxy known as supernova feedback,
which has a fundamental role in regulating star formation
<cit.>. Supernova feedback affects the evolution of the
galaxy population at all galactic masses, but is particularly
important in low-mass haloes which have shallow potential wells;
supernova-driven galactic winds affect both the abundance
<cit.> and inner structure of
low-luminosity galaxies.  Acting in conjunction with reionisation,
such winds strongly suppress galaxy formation in small haloes, reducing
the abundance of luminous low-mass galaxies
<cit.>.

Energy injection from supernova can potentially alter the inner
structure of dark matter haloes: if gas becomes gravitationally
dominant in the centre and most of it is removed suddenly, as could
happen in a starburst, energy can be transferred from the gas to the
dark matter and this can cause the centre to expand, turning the
original NFW cusp into a core. This mechanism, first proposed in the
1990s <cit.>, became fashionable again several years later
<cit.>
when tentative observational evidence for cores, particularly in dwarf
galaxies, began to emerge <cit.>. This evidence, however, is
controversial <cit.>. While the proof of concept in
<cit.> was based on a single explosive event, recent
simulations have shown that repeated outflows can create rapid
fluctuations in the gravitational potential which efficiently transfer
energy to the dark matter <cit.>. This core-formation
mechanism depends on the details of the baryon physics implemented in
the simulation <cit.> and not all cosmological
simulations produce cores in dwarf galaxies <cit.>.  On
scales larger than dwarfs, energy injection by AGN has been invoked as
a mechanism for core formation; however, the conditions required to
alter the deep potential wells of massive galaxies appear quite
extreme
<cit.>.

Energy injection into subhaloes: tidal effects from baryonic
  structures.- In Section <ref> we described the tidal
effects that the host halo induces on the dynamics and structure of
subhaloes. The presence of a central galaxy enhances these effects
both in subhaloes and in the satellite galaxies within them,
particularly when their orbits cross the region where the central
galaxy dominates the tidal field. Tidal shocking by a galactic disc
can result in the total disruption of subhaloes around the central
regions of the host <cit.> and other
structural changes. Current hydrodynamical simulations of
Milky-Way-like galaxies and their environment seem to agree that the
overall effect is a substantial reduction in the number of subhaloes
near the centre <cit.>.

There has been great progress in the past decade in incorporating
baryonic physics into full cosmological simulations; today galaxy
formation and evolution can be modelled in unprecedented detail
<cit.>. In
this way the effect on the dark matter phase-space distribution of the
complex interplay between the cooling and heating mechanisms of
baryons described above can be studied in their full cosmological
setting. In spite of this undeniable progress, many aspects of
baryonic physics remain poorly understood and, when they involve
processes on scales below the resolution of the simulation, they need to be
included as a subgrid model.  There are different approaches to
this problem which are often difficult to validate and this translates
into substantial uncertainties in some of the predictions of the
simulations (see <cit.> for a discussion of the limitations of
gasdynamic simulations).





 §.§ Astrophysical tests of the nature of the dark matter


Laboratory searches for dark matter have so far proved
unsuccessful. This, and the failure to find evidence for Supersymmetry, has generated gloom amongst proponents of the
lightest stable supersymmetric particle as the dark matter (even
though the mass of the Higgs boson suggests that the supersymmetry
scale is likely to be larger than a few tev, beyond the
reach of the lhc).  There have been, however, claims that both
CDM-WIMPs, and WDM particles in the form of sterile neutrinos of
mass 7 kev, have been discovered, the former through γ-ray
annihilation radiation from the Galactic Centre <cit.>,
the latter through a 3.5 kev decay line in the x-ray
spectra of galaxies and clusters
<cit.>[Contrary to some claims, xmm data for Draco, and Hitomi data for Perseus, are consistent
  with a 7 kev neutrino <cit.>.].
These claims are highly controversial but, since cosmogonic models
based on such particles have strong predictive power, they are
disprovable with appropriate astrophysical observations.

The standard CDM model has naturally come under the closest
scrutiny. Perhaps the two most important predictions of this model
(derived from N-body simulations) are: (i) the existence of a vast
population of haloes and subhaloes which, below a mass of order
10^9 M_⊙, are dark; (ii) the presence (in the absence of the
baryon effects discussed in the preceding section) of a steep cusp
(ρ∝ r^-1) in the density profile of dark matter haloes of
all masses. These two predictions are related to three of the much
publicized four problems of the CDM model on subgalactic
scales (often referred to as the “small-scale crisis” of
CDM): the (i) missing satellites; (ii)
too-big-to-fail and (iii) core-cusp problems. The fourth is the
so-called (iv) planes of satellites problem. Indeed some of
alternative dark matter particle models, such as SIDM, have been
proposed specifically to solve some or all of these perceived
astrophysical problems.

The missing satellites problem is the discrepancy between the
relatively small number of satellites observed around the Milky Way
and M31 and the many orders of magnitude larger number of halo
substructures predicted by CDM N-body simulations
<cit.>. The “too-big-to-fail” problem is the
existence in CDM N-body simulations of massive, dense galactic
subhaloes (maximum circular velocities, V_ max > 30 km/s) whose
kinematics appear inconsistent with those of the brightest Milky Way
satellites <cit.>.  The core-cusp problem is the
discrepancy between the cuspy universal NFW density profiles predicted
for pure CDM/WDM haloes and the inference of central cores in some
galaxies, particularly dwarfs <cit.>.  The “planes
of satellites” problem is the arrangement of the bright satellites of
the Milky Way, M31 (and a few others) on a thin plane in which the
satellites seem to be coherently rotating and which have been claimed
to be incompatible with CDM <cit.>.

The first three of the four perceived problems can be solved once the
effects of baryons discussed in Section <ref> are taken
into account. Perhaps paradoxically, the solution to what later became
known as the “missing satellites” problem was understood long before
it came to be regarded as a problem for CDM.  The strong suppression
of galaxy formation in haloes below a mass of ∼ 10^10M_⊙
was originally calculated using semi-analytic techniques
<cit.>, as were the implications
for the abundance of galactic satellites in the CDM model
<cit.>. This solution has been
repeatedly confirmed by modern gasdynamic simulations
<cit.>. Similarly,
the “too-big-to-fail” problem disappears when baryons are taken into
account, in this case through the more subtle effect of the reduced
growth of subhaloes arising from the early loss of baryons mentioned
above <cit.>. The “core-cusp problem”, if it exists at
all, can also be solved by the type of explosive baryonic effects
discussed in Section <ref>, which can transform NFW cusps
into cores[Other baryon effects that can transform cusps into
  cores have been proposed (e.g. <cit.>)
  but have been less studied.].  The existence of “planes of
satellites” in the Milky Way and M31 turns out not to be as unlikely
as has been claimed <cit.>, once the
statistics are calculated rigorously, taking into account the “look
elsewhere effect” <cit.>[See <cit.>
  for an opposed view.]. The origin of these planes is almost
certainly the anisotropic nature of the accretion of satellites along
filaments of the cosmic web <cit.> although the
exact mechanism is still unclear as is the expected frequency of these
structures.

While it is now generally agreed amongst practitioners of the field
that CDM is not afflicted by a “missing satellite” or a
“too-big-to-fail” problem, the data on the satellites of the Milky
Way can be used to constrain alternative dark matter models,
particularly those with a cutoff in the primordial power spectrum. In
WDM, the cutoff lengthscale varies roughly inversely with the mass of
particle. Thus, if this mass is too small, then too few small-mass
haloes would form and their abundance could be too low to account for
the observed number of satellites of the Milky Way. The expected
subhalo abundance increases roughly in proportion to the mass of the
parent halo <cit.> so, in reality, the observed abundance of
satellites constrains both the particle mass and the host halo mass
simultaneously.  For instance, using a semi-analytic model of galaxy
formation, the thermal WDM model was found to be in conflict with the
data if the Milky Way halo mass is smaller than
1.1 × 10^12 M_⊙ <cit.>. Using a similar
approach, <cit.> have ruled out a significant fraction of
the parameter space of sterile neutrinos and conclude that the models
that are in best agreement with the observed 3.5 keV line require 
the Milky Way halo to have a mass no smaller than
1.5× 10^12 M_⊙, a value that may already be in conflict
with the most recent determinations of the Milky Way halo mass
<cit.>. We should note that since the central densities
of WDM haloes are lower than those of their CDM counterparts, the
“too-big-to-fail problem” is easily avoided in WDM
<cit.>.

Although the strongest constraints on the SIDM cross section come from
the shapes and dynamics of massive haloes (particularly of galaxy
clusters, see e.g. Table 1 of <cit.>), the Milky Way
satellites are perhaps the best testbed for SIDM, since it is in these
systems that the model shows its greatest promise as an alternative to
CDM.  A few years ago it was suggested that the interesting range of
cross sections for the SIDM model to alleviate the “too-big-to-fail”
problem (without taking into account the baryonic processes just
mentioned) is 0.1≲σ_T/m_χ≲10cm^2/g
<cit.>.  Since then, several studies have taken a closer
look at the properties of the Milky Way satellites within the context
of SIDM and the picture that is emerging points to promising tests for
the near future which will either strengthen SIDM as an alternative to
CDM or narrow the range of allowed cross sections. For instance, the
diversity of dark matter densities on subkiloparsec scales in the
Milky Way satellites is difficult to accommodate for SIDM cross
sections σ/m_χ∼1cm^2/g <cit.>. The
inferred high dark matter densities in the ultra-faint satellites
(albeit uncertain due to possible systematic effects) are at first
sight difficult to explain within SIDM, which naturally predicts
cores, particularly in low-mass dark matter-dominated haloes. However,
a gravothermal collapse phase in SIDM haloes has recently been
proposed <cit.> as
a mechanism to create a diverse population of dwarf-size haloes, some
of which would be cuspy (those that collapse), and others that would
have cores. If cores are indeed shown to be present in (some) dwarf
galaxies, then dark matter self-interactions and the explosive baryon
effects in CDM mentioned above provide alternative explanations that
need to be contrasted. A promising way to achieve this, recently put
forward <cit.>, is to search for distinct signatures in the
detailed kinematics of the stellar population as they respond
differently to these two core formation mechanisms, one impulsive
(supernova feedback) and the other adiabatic (SIDM).

Since, as we have seen, the simplicity of the predictions of N-body
simulations can be easily obscured by the complexity of baryon
effects, testing dark matter models with astronomical observations
might, at first sight, seem a hopeless task. In fact, this is not the
case: the vast majority of haloes in CDM (and in many alternative dark
matter models) are dark, that is, unaffected or almost unaffected by
baryons. It is the existence of a vast population of such small-mass
haloes (m ≲ 5× 10^9 M_⊙) that is the hallmark of
the CDM model that distinguishes it from, for example, the WDM model.
Fortunately, nature has provided us with several tools to detect dark
objects in the Universe. One of these takes advantage of a side effect
of cosmic reionization which allows haloes in a small mass window
(10^8 ≲ 5× 10^9) M_⊙ to retain neutral hydrogen in
hydrostatic equilibrium with the dark matter potential and in thermal
equilibrium with the ionizing UV background, gas which is, however,
too diffuse to make stars <cit.>. These objects called
RELHICs (REionization-Limited HI Clouds, <cit.>)
may be detectable in forthcoming blind HI surveys and provide, in
principle, a critical test of CDM and related models in a regime that
has not been proved before.

An interesting idea that has been proposed to infer the existence of
small dark subhalos orbiting in the Milky Way halo is the disturbance
they cause when they cross a tidal stellar stream <cit.>.
When a subhalo crosses a stream it induces velocity changes along and
across the stream that can give rise to a visible gap, particularly in
cold streams such as those stripped from globular clusters. The
cross-section for gap creation is dominated by the smallest subhalos
so gaps can, in principle, constrain the identity of the dark matter.
The creation of gaps has been investigated with analytical treatments
or idealized numerical studies and it has been suggested that
perturbers of mass ∼ 10^7 M_⊙ could be detected in the gd-1 and Pal 5 globular cluster stellar streams <cit.>. A
complication of this method is that perturbations on the streams can
be induced not only by dark subhaloes but also by giant molecular
clouds and the bar at the centre of the Milky Way
<cit.>. Recent deep imaging around the Pal 5 stellar
stream does indeed reveal significant disturbances, in particular two
gaps which have been attributed to the impact of subhalos of mass in
the range 10^6-10^7 M_⊙ and 10^7-10^8 M_⊙ respectively
(although the smaller gap could also be due to the impact of a giant
molecular cloud) <cit.>.

But perhaps the most direct method to search for the ubiquitous
small-mass dark haloes is gravitational lensing.  There are two
specific instances where strong gravitational lensing could provide
the means to do this.  The first are the “flux-ratio anomalies” seen
in some multiply-lensed quasars; the second are small distortions of
Einstein rings and large arcs.

In a multiply-lensed image, the magnifications are determined by high
order derivatives of the lensing potential and are therefore
particularly sensitive to small changes in the potential such as those
produced by intervening small-mass structures. If the mass
distribution of the lens is smooth, the ratios of the fluxes of close
images (formed when the sources are close to a fold or a cusp of the
caustic) follow a certain asymptotic relation
<cit.>. These smooth-lens relations are violated
if there are intervening structures or substructures in the lens
giving rise to flux-ratio anomalies, which probe the total amount of
mass in structures along the line of sight to the lens
<cit.>.  Flux-ratio anomalies have been
observed in several quadruply-lensed quasars but dark substructures
alone are insufficient to explain the observed anomalies
<cit.>, implying that other effects such as inadequate lens
modelling may be at work. With better modelling of the lens (including
stellar discs and luminous satellites), it has been possible to set a
lower limit to the mass of a thermal WDM particle (see
<cit.> and Harvey et al., in preparation), similar to the
limits from satellite counts discussed above and also to those derived
from the observed inhomogeneity of the gas distribution at high
redshift probed by the Lyman-α forest <cit.>.



A more direct strategy for detecting dark structures and substructures
is to search for distortions in strongly lensed images. When the
source (a background galaxy), the lens (a massive halo) and the
observer are perfectly aligned, a circular feature near the centre of
the lens, an Einstein ring, is formed; if the alignment is not
perfect, then giant arcs are formed.  If the lens is a halo of mass
larger than ∼ 10^13 M_⊙, the radius of the Einstein ring
is generally larger than the image of the central galaxy and can thus
be studied in detail. If a halo or subhalo happens to be projected
onto the Einstein ring, it too will gravitationally lens the light
from the source producing a small distortion in the image of the
Einstein ring or giant arc <cit.>.  This strategy has
already yielded a halo of ∼ 10^8 M_⊙  <cit.>
[This halo mass was estimated assuming a truncated
  pseudo-Jaffe profile (see e.g. eq. 42 in <cit.>). The
  inferred mass is likely to be larger if an NFW profile is assumed
  instead. For instance, a similar  dark matter substructure detected
  with lensing was reported by <cit.> with a mass of
  ∼3.5×10^9 M_⊙ assuming a truncated pseudo-Jaffe
  profile, while assuming an NFW profile this substructure is estimated to
  have a mass of ∼10^10 M_⊙  <cit.>.] and 
could detect haloes as small as ∼ 10^7 M_⊙
<cit.>.

Detecting the small signal generated by individual projected haloes or
subhaloes requires accurate modelling of the source and the lens (the
“macro” model; <cit.>) and sophisticated
statistical techniques to analyse the image residuals. Dark haloes
imprint other observable features onto strong arcs. For example,
distortions to the lensing potential caused by the cumulative
contribution of many hundreds of projected structures produce unique
correlated residuals in the lensed image, the nature of which is
dependent on the abundance and mass distribution of the halo
population and, therefore, on the nature of the dark
matter <cit.>.  The mass function of dark haloes may
also be detectable through the n-point functions of the
projected density field or the substructure convergence power
spectrum <cit.>.

A very attractive feature of strong lensing as a means to detect
small-mass objects is that, for lens configurations of interest, the
dominant source of strong arc distortions are field haloes along the
line of sight, rather than subhaloes resident in the lens
<cit.>[This has been demonstrated
  explicitly for the case of Einstein ring distortions but it may hold
  true for other tests as well.].  This makes this test uniquely
powerful because, as we have seen, the haloes of interest, of mass less
than ∼ 10^8 M_⊙, are completely dark: they have never been
modified in any way by bayrons. Thus, the test depends mostly on the
abundance of pristine “field” dark matter haloes which we know very
well how to calculate rigorously and precisely with N-body simulations
for cosmological models of interest.

Approximately a few hundred high quality strong lens systems would
suffice to rule out either the 7 keV sterile neutrino model or
CDM itself <cit.>. Very high resolution imaging is the
primary requirement, either in the optical or UV, or using
interferometry at submillimeter and longer wavelengths <cit.>.  At least
several tens of systems with high quality data are already available
and future imaging facilities such as LSST and Euclid will increase
the number of suitable strong lenses by orders of magnitude.  By
bypassing the complications introduced by baryons, which have spoiled
all previous efforts to test the CDM model unambiguously and
distinguish it from alternative models, be they on small or large
scales, gravitational lensing offers a unique opportunity for a
breakthrough in this quest from astrophysics evidence alone.

We are very grateful to Alejandro Benitez-Llambay,
  Sownak Bose, Jiaxin Han,Mark Lovell and Simon White for their
  valuable comments.  Special thanks to Simon for his careful reading
  of the manuscript, insightful comments and suggestions. All of these
  have significantly improved our paper. We also thank Sebastian Bohr for
  his help in creating the right panel of Fig. 2.
  JZ acknowledges support by a
  Grant of Excellence from the Icelandic Research fund (grant number
  173929-052). CSF acknowledges support from the European Research
  Council (ERC) Advanced Investigator grant DMIDAS (GA 786910) and the
  Science and Technology Facilities Council (STFC) [grant number
  ST/F001166/1, ST/I00162X/1, ST/P000541/1]. Some of the simulations
  reported in this paper were carried out on the DiRAC Data Centric
  system at Durham University, operated by the ICC on behalf of the
  STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded
  by BIS National E-infrastructure capital grant ST/K00042X/1, STFC
  capital grant ST/H008519/1, and STFC DiRAC Operations grant
  ST/K003267/1 and Durham University. DiRAC is part of the National
  E-Infrastructure.
  
 
  
yes




