Dark matter haloes and subhaloes

The development of methods and algorithms to solve the $N$-body problem for classical, collisionless, non-relativistic particles has made it possible to follow the growth and evolution of cosmic dark matter structures over most of the Universe's history. In the best studied case $-$ the cold dark matter or CDM model $-$ the dark matter is assumed to consist of elementary particles that had negligible thermal velocities at early times. Progress over the past three decades has led to a nearly complete description of the assembly, structure and spatial distribution of dark matter haloes, and their substructure in this model, over almost the entire mass range of astronomical objects. On scales of galaxies and above, predictions from this standard CDM model have been shown to provide a remarkably good match to a wide variety of astronomical data over a large range of epochs, from the temperature structure of the cosmic background radiation to the large-scale distribution of galaxies. The frontier in this field has shifted to the relatively unexplored subgalactic scales, the domain of the central regions of massive haloes, and that of low-mass haloes and subhaloes, where potentially fundamental questions remain. Answering them may require: (i) the effect of known but uncertain baryonic processes (involving gas and stars), and/or (ii) alternative models with new dark matter physics. Here we present a review of the field, focusing on our current understanding of dark matter structure from $N$-body simulations and on the challenges ahead.

dimensionless power spectrum, and the spectral index n s = 0.965 [1]. 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 [22]) 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 5 , 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 6 (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 [24]. Another example is collisional damping due to interactions between dark matter and relativistic particles in the early universe (either photons or neutrinos, e.g. [25,26], or, in non-standard models, dark radiation in hidden dark sector models, e.g. [20,21]). 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 7 . 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. 1. 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 8 . Measurements of galaxy clustering on scales larger than individual galaxies, together with constraints from the flux spectrum of the Ly-α forest (e.g. [28]) constrain the power spectrum to be like CDM to the left of the hashed area in Fig. 1. On smaller scales the power spectrum could have a damping cutoff due to either collisionless (as in WDM models; red 5 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 B T 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. 6 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 [23]. 7 Note that acoustic oscillations are also present in WIMP-CDM models (e.g. [27]), but they occur at much smaller scales than in relevant hidden dark sector models where they can be of galactic scale. 8 We use M = 4π/3ρ(π/k) 3 , where ρ is the mean dark matter density today. 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 [29]). Figure 1. Dimensionless linear dark matter power spectrum in different dark matter models. In the current paradigm, cold dark matter (CDM), the power spectrum keeps on rising to well below subgalactic scales. Alternative models such as warm dark matter (WDM) or interacting dark matter (DAOs) have a cutoff at or sligthly below galactic scales, which determines the abundance and structure of small-mass dark matter haloes and subhaloes and the galaxies within. In the black hashed area, the dark matter is constrained by the observed large-scale distribution of galaxies [e.g. 30,31] and the Ly-α forest constraints on WDM [28] to behave as CDM. Figure adapted from [32].
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. [33]). 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 [34,35]. (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 [36] and Ellipsoidal collapse [37] models which link a primordial perturbation to the final equilibrium configuration: the dark matter halo. (iii) The halo model (for a review see [38]), 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 ( [39] see also [40]), 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 [41][42][43]. (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 [44]. 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 [45].
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): where d/dt is the Lagrangian derivative. Cosmological N-body simulations 9 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 10 distribution function: 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 ε 11 , 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 12 ; 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 13 . Since each simulation particle represents a region of phase space containing a very large number, m i /m χ , of real dark matter particles, 9 For a review see e.g. Section 3 of [46]. 10 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. 11 In principle each particle can have an individual softening, see e.g. Section 4 of [47]. 12 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. 13 In Fourier space, Eq. 6 is simply a multiplicationΦ( k) =ĝ( k)ρ( k). 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. [48,49]) or direct integration of the N 2 interactions (e.g. [50]). 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 3 M method (e.g. [51,52]), 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 [53] 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 [54], which uses the treePM algorithm, whereby short-range forces are computed with the tree method and long-range forces with Fourier techniques 14 .
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 (1-6) 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 2.3 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. 1, 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. [55][56][57][58][59]; see Appendix A of [58] for a detailed derivation). An alternative to the N-body approach is the "gravothermal fluid" approximation [60], 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. [61]). 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. [62,63]).

The non-linear regime: initial conditions and the emergence of the cosmic web
The techniques of Section 2.2 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 [14] and [64] and have been refined over the years (e.g. [65,66]; for a review see [67] or Appendix C1.1.4 of [68]). 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: 14 For a review of the force computation methods see Section 3.5 of [46]. 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 15 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 [70], or the more accurate second-order Lagrangian perturbation theory (e.g. [71]), 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. Figure 2. Illustration of the initial conditions for an N-body simulation. Left: the dimensionless linear CDM power spectrum. The vertical dashed lines mark the modes corresponding to the maximum and minimum scales that can be represented in the initial conditions: the fundamental mode, 2π/L, and the Nyquist mode, π/d, where L and d are the cube length and interparticle separation, respectively. Right: a realisation of the dark matter density field generated from the power spectrum on the left at redshift z = 127 using N = 1024 3 particles in a cosmological cube of co-moving side, L = 40 Mpc/h. The code MUSIC [65] was used to generate the particle distribution and the Pynbody package [72] to create the image.
An illustration of the end result of this procedure is shown in Fig. 2. The initial conditions generator, MUSIC [65], 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 16 , and the minimum lengthscale that can represented in the initial conditions is set by the Nyquist frequency of the particle 15 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 [69]. 16 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. distribution 17 . 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. 2 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. [73]). 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 [74] that requires careful treatment to either remove small-scale artificial clumps [75] or avoid their formation altogether by using non-standard simulation techniques [45,76]. non-linear CDM linear CDM Figure 3. Emergence of the cosmic web. Left: evolution of the (projected) dark matter density field in a slab of length L = 100 Mpc/h and thickness 15 Mpc/h from the Millennium II simulation [77].
The redshift corresponding to each snapshot is shown on the top right. Right: The dimensionless dark matter power spectrum (solid lines) at the redshifts shown on the left. For comparison, also shown are: the linear power spectrum (thin grey lines) and the non-linear power spectrum for the lower resolution but larger scale (500 Mpc/h) Millennium I simulation (in dotted lines; [4] 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 [77], is illustrated in Fig. 3. 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 17 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, ε. 18  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 [14,78] 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 O(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 O(10 9 ) particles are routinely performed 19 , and the most expensive simulations to date have reached the one trillion particle milestone [80].  [4]. The small slice at the top shows the CfA2 "Great Wall" [81], with the Coma cluster at the centre. Just above is a section of the Sloan Sigital Sky Survey in which the "Sloan Great Wall" [82] is visible. The wedge on the left shows one-half of the 2-degree-field galaxy redshift survey [83]. At the bottom and on the right, mock galaxy surveys constructed using a semi-analytic model applied to the simulation [84] are shown, selected to have geometry and magnitude limits matching the corresponding real surveys. Adapted from [85].
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" [14,86]. Physically based models of galaxy formation that could be grafted onto N-body simulations were developed in the early 1990s [87]. 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, 19 For a review on the state of cosmological simulations circa 2012 see [79]. 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 3.1). Semi-analytic models have been extremely successful in linking the distribution of dark matter computed in an N-body simulation to the observed Universe [88][89][90][91] and have become very sophisticated in predicting visible galaxy properties over a large range of wavelengths (e.g. [92]). An example based on the Millennium simulation is shown in Fig. 4.
Dark matter N−body simulations are the cornerstone of the current understanding of how galaxies form and evolve and, as illustrated in Fig. 4, have been very successful in explaining the large-scale structure of the Universe [85]. 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 [93], 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 [14]. By now, it is firmly established that whatever the dark matter is, it must behave as CDM on large scales (see Figs. 1 and 4). 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 1 and 2.1. 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 [36]. 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 3. 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 [94]. A recent update of this work (applied also to WDM) puts the fraction even lower, at the percent level [95].
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 [14,96,97]. 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 [14]; (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 [98]. The overdensity, ∆(z), is a redshift-dependent function of cosmology [99,100], 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 [98]. 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 [e.g. 77, [101][102][103][104][105][106], and is now well determined over the full range of epochs and masses relevant to galaxy formation, from O(10 8 M ) dwarf-size haloes to O(10 15 M ) cluster-size haloes. The number density of haloes per unit mass scales as: with an overall normalisation that correlates with the large scale environment, with denser environments having a larger halo abundance [78,107]. 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 [108] and extended Press-Schechter (EPS) formalisms [109,110], 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 [37]. 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 [e.g. 45,75,111,112] and interacting dark matter [e.g. 20,113,114] models have now been fairly well characterized with N-body simulations (with appropriate corrections for spurious fragmentation due to particle discreteness near the cutoff [74,75]). The Press-Schechter approach can be readily extended to provide a reasonable approximation to the halo mass function in these models as well [e.g. [115][116][117]. Fig. 5 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 [118], is an example of a model with dark acoustic oscillations, while WDM is an well-known example of the free streaming effect (see Fig. 1). 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 20 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 [119,120]: 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: 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 [121], 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. 6). 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 [122][123][124][125]. 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. 8: 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 [120,[126][127][128][129][130][131][132][133][134][135][136]. 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. 1). 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 [137] 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 [120,131,135], 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 [132]. 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.  [143]. The top panels show the spherically averaged radial density (left; [143] 21 ) and velocity dispersion (right; [121] 22 ) profiles, which are nearly universal for haloes in dynamical equilibrium. The bottom panels show the halo shape (left: moment of inertia axis ratios, and triaxiality: T = (a 2 − b 2 )/(a 2 − c 2 ); [144] 23 ) and local dark matter velocity distribution near the solar circle: 2 kpc < r < 9 kpc (right; [145] 24 ).
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 [138]: 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 φ )/σ 2 r is the velocity anisotropy profile, which quantifies the degree of anisotropy of the particle orbits in the dark matter halo. Haloes tend to be 21  isotropic only in their centres and are radially anisotropic in their outskirts [121,139], with a velocity anisotropy that is related to the logarithmic slope of the density profile [140]. The velocity structure of dark matter haloes in equilibrium is thus intimately linked to their spatial distribution (see top right panel of Fig. 6), 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 [121,141], and is in remarkable agreement with the self-similar solution for infall onto a point mass in an Einstein-de Sitter universe [142]. 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 (E ), and is fully given by the Eddington formula [146]: where u = √ E − Ψ and E 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. 11. 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. 11). In other words, haloes are non-Maxwellian simply by virtue of their NFW density profiles. In fact, CDM haloes in simulations have local 25 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 ( [145]; see bottom right panel of Fig. 6).
Halo shapes.-Although to first order, CDM haloes are well described by the spherical NFW profile, they are, in fact, triaxial [78]. In general, CDM haloes are prolate in the inner parts and oblate in the outskirts ( [78,144,147,148]; see bottom left panel of Fig. 6). 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 [144,149]. 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 [150,151] 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 [152].
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, Milky-Way-size halo . Bottom left: spherically averaged density profiles. WDM haloes are well described by an NFW profile, but have lower concentrations than their CDM counterparts of the same mass; SIDM haloes develop flat density cores during a transient stage that inevitably ends with the collapse of the core once the gravothermal catastrophe is triggered. Bottom right: spherically averaged velocity dispersion profiles. WDM haloes still obey the universal scaling for the pseudo-phase-space density, ρ/σ 3 ∝ r −1.875 , at most radii, except in the very centre, which results from a similar velocity dispersion profile to that in CDM but shifted downwards and to the right as a result of the lower concentration. SIDM haloes develop isothermal density cores of size of the order of the scale radius.
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 [75,112,[154][155][156][157].
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 [158] which is therefore reflected in the concentration of the haloes. An example of this is shown in the bottom left panel of Fig. 7, 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. [157], which connects the concentration to the cutoff scale in the power spectrum 28 . This lower concentration is also reflected in the velocity dispersion profile (see bottom-right panel of Fig. 7).
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 [159]. 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, O(10 pc) for ∼ keV thermal WDM relics in dwarf-size haloes [160,161]; 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 [112].
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 , 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 formationthe 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. 7). 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. [162]). 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 [61,163] 29 .
For 0.1 σ T /m χ 10 cm 2 /g, cosmological N-body simulations have shown that SIDM haloes today should have cores of size of the order of the scale radius [56][57][58]162,167]. At the lower end of this range of cross-sections, SIDM cores are small, ∼ 0.2r s [58], 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 [168]. At the higher end of the 28 This functional form has been corroborated by [112], but the parameters in the two studies are different. The formula is nevertheless a good approximation to the general behaviour. 29 The gravothermal collapse [164] 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. [165]). 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 [166].
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. 7). 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 [169]. This makes haloes centrally rounder than their CDM counterparts [170] and causes them to have local velocity distributions that are close to Maxwellian [171]. Since the onset of the gravothermal collapse phase is expected to be ∼ (250 − 400)t rel,s [165], 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.
CDM Milky-Way size halo at z=0

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 [23,24], 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. [137,173]). 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 [94,137]. These calculations show that the amount of smooth accretion onto present-day CDM Milky Way-size haloes is in the range ∼ (10 − 20)% [94,137]. 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 [137].
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. 8). 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 [174][175][176], and thus the subhalo population today extends to radii a few times the current virial radius ( [176]; those systems beyond the virial radius are commonly known as backsplash haloes [174]).

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 [177]). 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. 8). 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 [178], 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 [178,179] where M 0 is the mass of the descendant and the fitting parameters (for the Millennium simulation) are given in Table 1 of [179] (see also [180]). 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. [181]. This analytical approach, calibrated to simulations [182], 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. [183]). Distribution of accretion times and orbital properties.-Eq. 12 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 [184], 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 [185]. 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 32 . Fig. 9 shows a sample of the orbital parameters of the subhalo population at the time of infall, calculated from high-resolution N-body simulations [186]. 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. 12), 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 [187]. The distributions of radial and tangential velocities of infalling satellites for Milky Way-size haloes are shown in the middle panels of Fig. 9. 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 2 r + 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. 9, is close to circular (see also [187]).

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 3.1, 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 [143] 34 . 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. 10 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 [188] is shown colour-coded according to the subhalo infall time. Way-size halo simulation (Via Lactea II [188]; figure adapted from [192] 35 ). Subhaloes are colour-coded according to their infall time (measured from z = 0). Subhaloes that are just being accreted are radially infalling, while those that were accreted earlier and have completed many orbits lose energy through dynamical friction and sink towards the halo centre. Right: the 2D radial phase-space structure of simulation particles in a different Milky Way-size halo simulation (Aquarius [143]; figure adapted from [190] 36 ). Each particle is color-coded according to the number of caustics it passes (roughly proportional to the number of orbits executed by a given particle). The top panel includes bound subhaloes, while the bottom one does not. In the latter, tidal streams from disrupted subhaloes are more clearly visible. 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; [143]). The colour in this case encodes the number of caustics that a given particle traverses 37 . 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. 10 37 Caustics represent folds in the fine-grained phase-space distribution function, which in CDM evolves according to the collisionless Boltzmann equation (Eq. 1). 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. 1, 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 2.4), 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. [189][190][191]).
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 Nbody 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 [193]. 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 [91]; Bound Density Maxima [194]; VOBOZ [195]; Amiga Halo Finder [196]. 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 [197,198].
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 [199] and ROCKSTAR [200].
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 [193]. 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 [201] 38 : 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. [202]), 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: [203][204][205]: 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. 14 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 t dφ host (R)/dR [206]. The left panel of Fig. 11 shows an example of tidal stripping in this slow (adiabatic) mode for a subhalo in a circular orbit around a static host potential (from [207]). As pointed out by [208], 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. 14 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. 14) 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 [203,208,209]. This process is well described by the "impulsive approximation" (see [210,211] 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 [212,213]). The approximation gives the specific energy change suffered by particles in the subhalo due to a tidal shock as 39 : 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. 11 (taken from [208]). Although the energy injection in Eq. 15 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 [209,214]. 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. [214], 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. 11).
For a general subhalo orbit, a combination of the adiabatic (Eq. 14) and impulsive (Eq. 15) 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 39 It is possible to evaluate Eq. 15 for a given fixed potential and a given subhalo orbit (see e.g Eq. 20 of [208] for an NFW halo). 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 [208]).

Tidal stripping (constant tidal field strength)
Tidal shocking (impulsive) model ___ t=0 density profile evolution Figure 11. Tidal effects in subhaloes. Left: evolution of a subhalo in a circular orbit in a static host halo potential. Since the tidal field strength is constant, the subhalo gradually loses mass (red particles are bound to the subhalo, black particles are unbound) as it orbits in the host halo creating characteristic tidal streams (figure adapted from [207] 40 ). Right: the effect of tidal shocks. For nearly radial orbits, the variations in the potential near pericentre are rapid (relative to the internal dynamical timescale of the subhalo) and this leads to an impulsive tidal shock, which causes a drastic removal of mass (upper right) and a change in the dark matter distribution (bottom right). In the upper panel the fraction of stripped mass as a function of circularity (see Section 3.2), given by the impulsive approximation, is compared with that in a controlled simulation (figure adapted from [208] 41 ). The model works quite well for radial orbits but it fails for circular orbits (as in the left panel), for which an adiabatic model is more appropriate (Eq. 14). In the lower panel, tidal shocking is seen to reduce the mass in the central regions but preserves the asymptotic NFW behaviour, while the outer regions become considerably steeper than NFW (figure adapted from [214] 42 ).
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; [215]). 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 40 Reproduced from Frank C van den Bosch and Go Ogiya. 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 [208]. 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. 12).

Gyr 7.3 Gyr
mass ratio and circularity dependence 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 43  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 [216], which, for the drag force, F df , takes the form: where ρ(< | v s |) is the density of ambient particles with speed less than | v s | and lnΛ ≡ ln , 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 [e.g., 68]: 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. 17 requires a correction that scales with the circularity, η, as t df ∝ η γ η , where γ η ∼ 0.53 [217], implying that more eccentric orbits decay more rapidly (see right panel of Fig. 12). Eq. 16 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, [219] carried out a series of idealised N-body simulations of a subhalo infalling into a host, both described by a Hernquist density profile, 45 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. 12): 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 46 . Eq. 18 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 47 .

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 [221]. 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 2.4) with a few relevant modifications.
radial distribution of subhaloes global subhalo abundance 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 [143,[224][225][226][227][228]: dn sub /dm sub ∝ m α sub , where α ∼ −1.9 (see Eq. 7). 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. 9. 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 [137,226,229]. 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 47 Eq. 18 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. 18 was not investigated outside this range and might not be valid there even in the absence of a central galaxy. 48 c AAS. Reproduced with permission. For the original article, please visit the following link. 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: [ [229][230][231][232] where the exponential cutoff accounts for the increasing rarity of subhaloes of mass close to that of the host halo mass 49 . 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 [231]. 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. 13). In this case, an accurate approximation is given by [222]: 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 2.1). 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 [143]. 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 [143,224,230,[233][234][235][236] 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. 13). 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 [143,176,230]. 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 [198], 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 [221] 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. [237]), 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 [223] and shown in the right-hand panel of Fig. 13), is approximated quite accurately by the following functional form [223]: 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 3.3. 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. 14), while the outer regions show a steep truncation at a radius approximately equal to the tidal radius given in Eq. 13 (see Fig. 15 of [143]). We should remark that, as has been found for field haloes [121], a better fit to the density profile of subhaloes is given by the 3-parameter Einasto profile [238] 50 : 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 [240,241]. Although for isolated haloes the parameters α E and r −2 can be related to the virial mass of the halo, M 200 [135,239], in a similar way as the halo (NFW) concentration is connected to the virial mass, the situation is less clear for subhaloes [240], 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 [143,225]: where H is the Hubble constant. Eq. 23 can be related to the standard scale density of the NFW profile (δ c in Eq. 9), and thus to the NFW concentration, in a straightforward way [225]: host halo subhaloes self-similarity of the (sub)halo density profile structural parameters: subhaloes vs haloes structural parameters: halocentric dependence …. subhaloes r < r 50 __ haloes …. subhaloes r < r 50 __ haloes Figure 14. The inner structure of subhaloes. Left: spherically averaged density profile of subhaloes (which is remarkably similar to that of isolated haloes). The plot shows the density profile of a Milky Way-size halo (solid black line) and eight of its largest subhaloes (colour lines). The vertical dotted line marks the radius beyond which the simulation results are converged. The self-similarity in the central region is better appreciated in the inset where the density and radius are scaled to their values at the scale radius, r s . The figure is for the Via Lactea II simulation and is adapted from [188]. Upper right: mean relation between the maximum circular velocity, V max , and the radius at which it is achieved, r max , for subhaloes within r 50 (the radius within which the mean enclosed density is 50 times the critical density) of one the Milky Way-size halo simulations in the Aquarius project, at different resolution levels (colour lines). The red dashed lines show the scatter (68% of the distribution) for the highest resolution level. The dotted line is a fit to the mean relation for subhaloes and lies systematically below the equivalent line for isolated haloes (solid line). Lower right: a measure of concentration for subhaloes (see Eq. 24) within different radial ranges, as given in the legend. The solid line corresponds to isolated haloes. Figures adapted from [143] 51 . Since the concentration of haloes (and subhaloes) is tightly correlated with their mass (see Section 2.4), Eq. 23 implies a tight correlation between V max and r max , which indeed has been found and characterised in simulations (see right panel of Fig. 14). 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 52 :  52 We note that there is a typo in the caption of Fig. 28 in [143], which gives the fitting function for δ V and m sub (5.8 × 10 8 M → 5.8 × 10 4 M ).
The fitting function for the r max − V max relation implied by Eq. 25 is shown as a dotted line in the upper right panel of Fig. 14, while the corresponding relation for isolated haloes in the Millennium I simulation is shown as a solid line [129] 53 . 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 [126]: in Fig. 14 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. 14), 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 3.3: while tidal stripping naturally reduces V max , it reduces r max even further [209] 54 , 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. [243] 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 [244]. 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 [245]. 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 [245], 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 [245], compared to ∼ r −1.9 for field haloes.
3.5. 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 2.4, 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. 19 and 20, is broken [246] 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 [221,246].
In many SIDM models the subhalo mass function remains largely unchanged as long as the interaction cross-section, σ T /m χ < 10 cm 2 /g [57,58,247]. 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 [57].
Inner structure of subhaloes.-The inner structure of WDM and SIDM subhaloes is rather similar to that of field haloes (see Fig. 7), 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 [246] 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 [247]. 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 [248].

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 55 . 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 56 , 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 [21,114] and inelastic SIDM 57 . 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 4.1 below. 55 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. [68,249]). 56 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 2.4). Although the regime of gravothermal collapse has been known for a couple of decades [162,163], a comprehensive analysis of this regime has yet to be carried out (see [32,248,250,251] for recent developments in this interesting regime). 57 There is a class of inelastic SIDM models in which the dark matter can have ground and excited states (e.g. [252]), 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 [166,253].
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 [207,208]; 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 [122][123][124] 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 [e.g. 254], 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 [87,255] 58 . 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. [87]), 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) [84,257]. 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 [258,259]. Even though the assembly of baryonic matter by mergers is not, in general, adiabatic, the simple adiabatic model remains a reasonable approximation [260]. 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 [e.g. 261,262].
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 [263,264]. 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 59 [249,[266][267][268][269][270][271][272][273]. 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 [249,274].
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 [87]. 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 [255,275,276] 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 [268,270,271].
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 [277], became fashionable again several years later [278][279][280][281][282][283][284][285] when tentative observational evidence for cores, particularly in dwarf galaxies, began to emerge [286]. This evidence, however, is controversial [287,288]. While the proof of concept in [277] 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 [281]. This core-formation mechanism depends on the details of the baryon physics implemented in the simulation [289] and not all cosmological simulations produce cores in dwarf galaxies [290]. 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 [291][292][293][294][295].
Energy injection into subhaloes: tidal effects from baryonic structures.-In Section 3.3 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 [296,297] 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 [298][299][300].
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 [301][302][303][304][305]. 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 [303] 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 [306], the latter through a 3.5 keV decay line in the X-ray spectra of galaxies and clusters [307,308] 60 . 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 [311,312]. 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 [313]. 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 [e.g. 314]. 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 [315][316][317].
The first three of the four perceived problems can be solved once the effects of baryons discussed in Section 4.1 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 10 M was originally calculated using semi-analytic techniques [263,264,266], as were the implications for the abundance of galactic satellites in the CDM model [268,270,271]. This solution has been repeatedly confirmed by modern gasdynamic simulations [e.g. 249,265,318,319]. 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 [249]. The "core-cusp problem", if it exists at all, can also be solved by the type of explosive baryonic effects discussed in Section 4.1, which can transform NFW cusps into cores 61 . The existence of "planes of satellites" in the Milky Way and M31 turns out not to be as unlikely as has been claimed [e.g. 316,322], once the statistics are calculated rigorously, taking into account the "look elsewhere effect" [323] 62 . The origin of these planes is almost certainly the anisotropic nature of the accretion of satellites along filaments of the cosmic web [325,326] 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 [327] 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 [328]. Using a similar approach, [329] 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 [330]. 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 [153,331].
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 [18]), 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 χ 10 cm 2 /g [168]. 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 χ ∼ 1 cm 2 /g [32]. 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 [32,248,250,332] 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 [333], 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 [334]. These objects called RELHICs (REionization-Limited HI Clouds, [335]) 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 [336]. 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 [337]. 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 [338]. 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) [339,340].
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 [341,342]. 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 [341,343,344]. Flux-ratio anomalies have been observed in several quadruply-lensed quasars but dark substructures alone are insufficient to explain the observed anomalies [345], 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 [346] 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 [347].
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 [348]. This strategy has already yielded a halo of ∼ 10 8 M [349] 63 and could detect haloes as small as ∼ 10 7 M [353,354].
Detecting the small signal generated by individual projected haloes or subhaloes requires accurate modelling of the source and the lens (the "macro" model; [e.g. 348,355]) 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 [356,357]. 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 [357].
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 [358,359] 64 . 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 [360]. Very high resolution imaging is the primary requirement, either in the optical or UV, or using interferometry at submillimeter and longer wavelengths [354]. 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.