Open AccessReview

Dark Matter Haloes and Subhaloes

Jesús Zavala

^1,*

and

Carlos S. Frenk

Center for Astrophysics and Cosmology, Science Institute, University of Iceland, Dunhagi 5, 107 Reykjavík, Iceland

Institute of Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK

Author to whom correspondence should be addressed.

Galaxies 2019, 7(4), 81; https://doi.org/10.3390/galaxies7040081

Submission received: 29 July 2019 / Revised: 19 September 2019 / Accepted: 20 September 2019 / Published: 25 September 2019

(This article belongs to the Special Issue The Role of Halo Substructure in Gamma-Ray Dark Matter Searches)

Download

Browse Figures

Versions Notes

Abstract

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.

Keywords:

dark matter; structure formation; cosmological N-body simulations

Contents
1 Introduction	2
2 Formation of Dark Matter Haloes	3
2.1 Initial Conditions: The Primordial Power Spectrum in the Linear Regime	3
2.2 The Non-Linear Regime: N-Body Simulation Methods	5
2.3 The Non-Linear Regime: Initial Conditions and The Emergence of the CosmicWeb	7
2.4 The Structural Properties of Dark Matter Haloes	11
3 Halo Mergers and the Emergence of Subhaloes	18
3.1 Halo Mass Assembly: Smooth Accretion vs Mergers	18
3.2 Evolution of Subhaloes: Initial Conditions	19
3.3 Dynamics of Subhaloes	21
3.4 The Abundance, Spatial Distribution and Internal Structure of Dark Matter Subhaloes	28
3.5 The Impact of the Nature of the Dark Matter	32
4 Outlook	33
4.1 The Impact of Baryonic Physics on Dark Matter Structure	34
4.2 Astrophysical Tests of the Nature of the Dark Matter	36
References	39

1. Introduction

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

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

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

Dark matter could become non-relativistic at sufficiently late times to suppress, by free-streaming, the formation of low-mass galactic-scale haloes. This case is, in fact, one of the best-studied alternatives to CDM, known as warm dark matter (WDM). In contrast to WIMPs, these particles have masses of

O

(1 keV). A sterile neutrino, included as a part of a model that accounts for neutrino masses and for the baryon asymmetry of the universe, is the favorite WDM candidate (for a recent review see [16]). Another possibility is that dark matter is made of extremely light bosons with a

O

(1 kpc) de Broglie wavelength, in which case quantum effects would be relevant on galactic scales (such possibility falls in the category of “fuzzy dark matter”; for a review see [17].)

Although the interactions between dark matter and standard model particles are severely constrained, the interactions among the dark matter particles themselves are not. It is possible that dark matter may have its own rich phenomenology hidden from the ordinary matter. This hidden dark matter sector might possess new forces and particles, some of which could be viable dark matter particles that are strongly self-interacting2. These collisional particles fall under the category of self-interacting dark matter (SIDM; for a review see [18]3). Some of the hidden particles might be light enough that they effectively act as dark radiation that prevents the gravitational collapse of dark matter on subgalactic scales4 (e.g., [20,21]). As mentioned earlier, the CDM hypothesis is only supported to some extent: astronomical data allow a variety of models in which dark matter behaves significantly differently from CDM.

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

2. Formation of Dark Matter Haloes

2.1. Initial Conditions: The Primordial Power Spectrum in the Linear Regime

A theory of structure formation aims to explain the evolution of the universe from a nearly homogeneous initial state, with tiny matter density perturbations,

δ ρ / ρ

, seeded by inflation, which grow to leave an imprint on the CMB (emitted at the time of recombination,

z \sim 1100

, when

δ ρ / ρ \sim 10^{- 3}

), through the emergence of the self-gravitating dark matter haloes where galaxies form (

δ ρ / ρ ≫ 1

), to the universe we observe today characterized by a web of filamentary large-scale structure (

δ ρ / ρ \sim 1

The starting point is the end of cosmic inflation when dark matter perturbations are predicted to have a nearly scale-invariant power spectrum,

Δ^{2} \propto k^{3 + n_{s}}

, where

Δ^{2} (k) = k^{3} P (k) / 2 π^{2}

is the 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_{χ} \sim 100

GeV,

k_{fs} \sim 2.5 \times 10^{6}

h/Mpc); warm (

m_{χ} \sim 1

keV,

k_{fs} \sim 3.8

h/Mpc); and hot (

m_{χ} \sim

30 eV,

k_{fs} \sim 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 \times 10^{5}

–

1.2 \times 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 Figure 1. On the largest scales, not affected by the Mézáros effect, the power spectrum is nearly scale-invariant,

Δ^{2} \propto 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}

_{⊙}

, roughly corresponding to a (co-moving) wavenumber ∼12 h/Mpc8. 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 Figure 1. On smaller scales the power spectrum could have a damping cutoff due to either collisionless (as in WDM models; red line) or to collisional (as in models with DAOs; blue line) processes. We have not included the cutoff characteristic of fuzzy dark matter models, but we note that it is also oscillatory like the DAOs models (but due to quantum rather than collisional effects; see e.g., Figure 2 of [29]).

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) \propto 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 (

δ ρ / ρ \sim 0.1

), perturbation theory breaks down since gravity couples perturbations on different scales and their evolution can no longer be calculated as independent modes.

2.2. The Non-Linear Regime: N-Body Simulation Methods

To follow the evolution of dark matter density perturbations beyond the linear regime, several 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 virialized 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,40]), which assumes that the number of neighboring dark matter particles within a fixed physical separation remains constant, and can be used to study the deeply non-linear regime; a recent reformulation 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 discretization 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 artificial 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 (\vec{x}, \vec{v}; t)

, which in principle is given by the collisionless Boltzmann equation coupled with the Poisson equation for the gravitational field,

Φ (\vec{x})

(a combination known as the Vlasov-Poisson equation):

\begin{matrix} \frac{d f}{d t} & = & \frac{\partial f}{\partial t} + \sum_{i} v_{i} \frac{\partial f}{\partial x_{i}} - \sum_{i} \frac{\partial Φ}{\partial x_{i}} \frac{\partial f}{\partial v_{i}} = 0 \end{matrix}

(1)

\begin{matrix} ρ_{χ} (\vec{x}; t) & = & \int f (\vec{x}, \vec{v}; t) d^{3} \vec{v} \end{matrix}

(2)

\begin{matrix} \nabla^{2} Φ (\vec{x}) & = & 4 π G ρ_{χ} (\vec{x}) \end{matrix}

(3)

where

d / d t

is the Lagrangian derivative. Cosmological N-body simulations9 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”,

{\vec{x_{i}}, \vec{v_{i}}}

, with

i = 1, . . ., N

. Since the collisionless Boltzmann equation implies that the phase-space distribution remains constant in time along any trajectory

{\vec{x} (t), \vec{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-grained10 distribution function:

\begin{matrix} \tilde{f} (\vec{x}, \vec{v}) & \sim & \sum_{i} m_{i} W (| \vec{x} - \vec{x_{i}} |; ε) δ^{3} (\vec{v} - \vec{v_{i}}); \frac{d \tilde{f}}{d t} = 0 \end{matrix}

(4)

\begin{matrix} \tilde{ρ} (\vec{x}) & = & \int \tilde{f} (\vec{x}, \vec{v}) d^{3} \vec{v} \sim \sum_{i} m_{i} W (| \vec{x} - \vec{x_{i}} |; ε) \end{matrix}

(5)

\begin{matrix} \tilde{Φ} (\vec{x}) & = & \int g (\vec{x} - {\vec{x}}^{'}) \tilde{ρ} ({\vec{x}}^{'}) d^{3} {\vec{x}}^{'} \end{matrix}

(6)

where

m_{i}

is the mass of the simulation particle,

δ^{3} (\vec{v} - \vec{v_{i}})

is the DiRAC delta function 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 distribution12; 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 function13. Since each simulation particle represents a region of phase space containing a very large number,

m_{i} / m_{χ}

, of real dark matter particles, the information in an N-body simulation is always incomplete, limited by the phase-space resolution and the softening length.

With the discretization method employed in an N-body simulation, calculating the evolution of the phase-space distribution is reduced to following self-consistently the dynamics of a system of N particles (usually in terms of the Hamiltonian of the system in the co-moving frame) according to the potential derived from the particle distribution. Modern codes used to solve this problem employ efficient methods for computing the gravitational potential and integrating the Hamiltonian system forward in time. Early cosmological simulation codes used the particle-mesh (PM) technique in Fourier space (e.g., [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 techniques14.

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 Section 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 Equation (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 neighbors 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]).

2.3. 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,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:

(i): create a realization of an unperturbed cube of side L by distributing N particles homogeneously in a lattice or in a glass-like configuration15 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.

An illustration of the end result of this procedure is shown in Figure 2. The initial conditions generator, MUSIC [65], was used to construct the particle distribution on the right, which is a statistical realization 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 length scale that can be simulated is determined by the (co-moving) side of the computational cube16, and the minimum length scale that can represented in the initial conditions is set by the Nyquist frequency of the particle distribution17. 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 Figure 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].

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 Figure 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 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 performed19, and the most expensive simulations to date have reached the one trillion particle milestone [80].

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, 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 Figure 4.

Dark matter N-body simulations are the cornerstone of the current understanding of how galaxies form and evolve and, as illustrated in Figure 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 Figure 1 and Figure 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 Section 1 and Section 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.

2.4. The Structural Properties of Dark Matter Haloes

As a consequence of gravitational instability, 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 virialized 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

\sim 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 \sim 0.2

in units of the mean interparticle separation [14]; (ii) a spherical overdensity mass,

M_{Δ}

, contained within a sphere centered on the halo (with the center 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,

Δ \sim 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 characterized 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:

\frac{d n}{d M} \propto M^{α}, where α \sim - 1.9 for low masses,

(7)

with an overall normalization 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 virialized haloes as modeled 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 cutoff 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]).

Figure 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 a well-known example of the free-streaming effect (see Figure 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 conditions20 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]:

ρ_{NFW} (x) = \frac{ρ_{s} (c)}{c x {(1 + c x)}^{2},}

(8)

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} = 3 H^{2} / 8 π G

is the critical density of the universe, and:

δ_{c} = \frac{Δ c^{3}}{3 K_{c} (c)}

(9)

where

K_{c} (c) = \ln (1 + c) - c / (1 + c)

. Although recent simulations have shown that a different profile, the so-called Einasto profile, which has three parameters, is a slightly better fit to simulations [121], the asymptotic behavior of the NFW profile for

ρ (r \to 0) \sim r^{- 1}

remains a remarkably good approximation to the inner structure of CDM haloes (see top left panel of Figure 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 Equation (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 Figure 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.

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 [141]:

\frac{d (ρ σ_{r}^{2})}{d r} + 2 \frac{β}{r} ρ σ_{r}^{2} = - ρ \frac{d Φ}{d r},

(10)

where

Φ (r)

is the gravitational potential related to the density by Poisson’s equation,

σ_{r} (r)

is the radial velocity dispersion profile, and

β (r) = 1 - (σ_{θ}^{2} + σ_{ϕ}^{2}) / σ_{r}^{2}

is the velocity anisotropy profile, which quantifies the degree of anisotropy of the particle orbits in the dark matter halo. Haloes tend to be isotropic only in their centers and are radially anisotropic in their outskirts [121,142], with a velocity anisotropy that is related to the logarithmic slope of the density profile [143]. The velocity structure of dark matter haloes in equilibrium is thus intimately linked to their spatial distribution (see top right panel of Figure 6), which is strikingly evident in the so-called pseudo-phase-space density,

Q \equiv ρ / σ^{3}

, where

σ^{2} = σ_{r}^{2} + σ_{θ}^{2} + σ_{ϕ}^{2}

is the square of the 3D velocity dispersion. This quantity has been found to be an almost perfect power law,

Q \propto r^{- 1.875}

, over several orders of magnitude [121,144], and is in remarkable agreement with the self-similar solution for infall onto a point mass in an Einstein-de Sitter universe [145]. The radial behavior 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]:

f (E) = \frac{1}{\sqrt{8} π^{2}} \int_{0}^{\sqrt{E}} d u \frac{d^{2} ρ}{d Ψ^{2}} (r (Ψ (u))),

(11)

where

u = \sqrt{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 Equation (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,

ρ \propto r^{- 2}

, results in a Maxwellian distribution in Equation (11)). In other words, haloes are non-Maxwellian simply by virtue of their NFW density profiles. In fact, CDM haloes in simulations have local25 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 ([140]; see bottom right panel of Figure 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,139,147,148]; see bottom left panel of Figure 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 [139,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, i.e., 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, 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 characterized 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 Figure 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 Equation (39) of Reference [157], which connects the concentration to the cutoff scale in the power spectrum28. This lower concentration is also reflected in the velocity dispersion profile (see bottom-right panel of Figure 7).

It is interesting to note that the pseudo-phase-space density,

Q \equiv ρ / σ^{3}

, in WDM haloes scales with radius in the same way as in CDM,

Q \propto 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}

t_{rel, s} \sim {(ρ (r_{s}) 〈 v_{rel, s} 〉 σ_{T} / m_{χ})}^{- 1}

, where

〈 v_{rel, s} 〉

is the characteristic local relative velocity, becomes comparable to the age of the inner halo. The original CDM density cusp turns into a core within the region where this condition is satisfied. The interaction cross-section thus introduces a new scale in structure formation—the mean free-path for particle collisions—which breaks the near universality of CDM haloes. The transformation of the cusp into a core due to elastic collisions at the halo center 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 Figure 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 center (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

^{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.2 r_{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 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 such as this is shown in the bottom panels of Figure 7). Within the thermalized region, the orbits of dark matter particles are isotropized 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.

3. 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.

3.1. 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,172]). 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 Figure 8). One could argue for a more physical definition, based for instance on 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 [173,174,175], and thus, the subhalo population today extends to radii a few times the current virial radius ([175]; those systems beyond the virial radius are commonly known as backsplash haloes [173]).

3.2. 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 Figure 8). Thus, a merger tree contains information about the accretion times and masses of the haloes that eventually become subhaloes. Both 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,

d N_{m} / d ξ / d z

[178], which gives the mean number of mergers,

d N_{m}

, per mass ratio,

d ξ

(relative to the main progenitor at the time of accretion), per redshift interval,

d z

. This quantity has been found to have a functional form that is nearly universal [178,179]:

\frac{d N_{m}}{d ξ d z} (M_{0}, ξ, z) = A {(\frac{M_{0}}{10^{12} M_{⊙}})}^{α} ξ^{β} \exp [{(\frac{ξ}{\tilde{ξ}})}^{γ}] {(1 + z)}^{η_{z}},

(12)

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 realizations 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—Equation (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 idealized 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} = \sqrt{G M_{200} / r_{200}}

. Another common choice is to characterize 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 accretion32.

Figure 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}

_{⊙}

. 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 (Equation (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 Figure 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_{r}^{2} + V_{θ}^{2})}^{1 / 2}

, is determined by the potential of the host halo. Therefore, the ridge line of the bivariate distribution, (

V_{r}, V_{θ}

), shown in the right panel of Figure 9, is close to circular (see also [187]).

3.3. 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 [138]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.

Figure 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 color-coded according to the subhalo infall time. This figure reveals the path of a subhalo as it orbits in the host halo: it first falls in radially and 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 center 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; [138]). The color in this case encodes the number of caustics that a given particle traverses37. 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 Figure 10 since particles that are part, or were part, of a subhalo have undergone more particle orbits in their earlier host. This plot thus shows the richness of the substructure present in haloes today. As we mentioned earlier, most of the matter in a halo today has been accreted through mergers and consists of material that was stripped from subhaloes.

Identifying substructure–-Several algorithms are in common use to identify subhaloes in N- body simulations and define their boundaries and properties. These subhalo finders are based on different techniques. Here we will list only the most popular ones; for a comprehensive comparison study see [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 center. 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]. A similar procedure is used in SURV [176,199], which uses the merger history tree of a halo across all its branches to identify subhalos.

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 localized in real space, but have clear signatures in phase space. Examples of phase-space finders are the Hierarchical Structure Finder [200] and ROCKSTAR [201].

Current subhalo finders can 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 [202]38:

r_{t} = R {[\frac{m_{sub} (< r_{t})}{(2 - dln M / dln r) M (< R)}]}^{1 / 3} .

(13)

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., [203]), 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 modeled as: [204,205,206]:

d m_{sub} = α_{t} m_{sub} (> r_{t}) \frac{d t}{t_{orb} (R)},

(14)

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. Equation (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) / d R

[207]. The left panel of Figure 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 [208]). As pointed out by [209], a combination of the ill-defined tidal radius and uncertainty in the parameter

α_{t}

, makes the modeling of this adiabatic case quite complicated, with the end result that models based on Equation (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 (Equation (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 [204,209,210]. This process is well described by the “impulsive approximation” (see [211,212] 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 [213,214]). The approximation gives the specific energy change suffered by particles in the subhalo due to a tidal shock as39:

{(Δ E)}_{i, tid} = {(Δ v)}_{i, tid}^{2} \approx {|\int_{orbit} {\vec{a}}_{i, tid} (t) d t|}^{2},

(15)

where

{\vec{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 Figure 11 (taken from [209]).

Although the energy injection in Equation (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 center expand, reducing the inner density, although this process is not strong enough to create a central core [210,215]. 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 Reference [215], using idealized simulations in a static external halo potential, found a profile of the form

ρ \propto r^{γ} \exp (- r / r_{b})

, with a central slope,

γ \sim 1

, and a cutoff radius due to tidal shocks,

r_{b}

(see bottom right panel of Figure 11).

For a general subhalo orbit, a combination of the adiabatic (Equation 14) and impulsive (Equation (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 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 [209]).

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; [216]). A similar impulsive approach to the one above can be used to estimate the strength of this form of tidal heating, but the calculation is more complicated because, among other things, the distribution of subhaloes in the host and the encounter rate need to be modeled. A recent study finds that tidal shocking from encounters is subdominant (by a factor of several) compared to shocking during pericentric passages [209].

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 Figure 12).

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 analyzed as a series of uncorrelated sequential encounters between the object of mass,

M_{s}

, and velocity,

{\vec{v}}_{s}

, and particles randomly taken from the ambient medium with velocity distribution,

f_{m} (\vec{u})

. These interactions occur on a timescale much shorter than the variations in the velocity,

{\vec{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 [219], which, for the drag force,

{\vec{F}}_{df}

, takes the form:

\begin{matrix} {\vec{F}}_{df} = M_{s} \frac{d {\vec{v}}_{s}}{d t} & = & - 16 π^{2} M_{s}^{2} m \ln Λ [\int_{0}^{| {\vec{v}}_{s} |} f_{m} (| \vec{u} |) u^{2} d | \vec{u} |] \frac{{\vec{v}}_{s}}{| {\vec{v}}_{s} |} \\ = & - 4 π {(\frac{G M_{s}}{v_{s}})}^{2} \ln Λ ρ (< | {\vec{v}}_{s} |) \frac{{\vec{v}}_{s}}{| {\vec{v}}_{s} |}, \end{matrix}

(16)

where

ρ (< | {\vec{v}}_{s} |)

is the density of ambient particles with speed less than

| {\vec{v}}_{s} |

and

\ln Λ \equiv \ln [1 + {(b_{\max} / b_{90})}^{2}]

is the Coulomb logarithm, with

b_{90} = G (M_{s} + m) / v_{\infty}^{2}

v_{\infty}

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 center 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]):

\frac{t_{df}}{t_{dyn}} \approx \frac{1.17}{\ln (M_{h} / M_{s})} (\frac{M_{h}}{M_{s}}),

(17)

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, Equation (17) requires a correction that scales with the circularity,

η

, as

t_{df} \propto η^{γ_{η}}

, where

γ_{η} \sim 0.53

[220], implying that more eccentric orbits decay more rapidly (see right panel of Figure 12).

Equation (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. An example of this is provided in [218] with a series of idealized N-body simulations of a subhalo infalling into a host, both described by a Hernquist density profile45. This study 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 Figure 12):

\frac{t_{df}}{t_{dyn}} = A \frac{{(M_{h} / M_{s})}^{b}}{\ln (1 + M_{h} / M_{s})} \exp [c η (E)] {[\frac{r_{c} (E)}{r_{200}}]}^{d},

(18)

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 one46. Equation (18) was found to be valid over a wide range of orbital parameters; the most relevant restriction is

0.025 \leq M_{s} / M_{h} \leq 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 time47.

3.4. 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 [222]. 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.

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 [138,225,226,227,228,229]:

d n_{sub} / d m_{sub} \propto m_{sub}^{α}

, where

α \sim - 1.9

(see Equation (7)). This similarity to the halo mass function is partly because most subhaloes identified at a given time were accreted relatively recently and thus tidal effects have not had time to act; see Figure 9. The normalization of the subhalo mass function depends on the mass of the host halo, with more massive haloes having, on average, larger subhalo populations [137,227,230]. This reflects the earlier assembly of low-mass haloes, which allows tidal effects more time to act and disrupt subhaloes. For similar reasons other properties of the host halo can have second-order effects on the amplitude of the subhalo mass function, e.g., at fixed mass, more concentrated haloes (which assemble earlier) have fewer subhaloes.

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

N (> μ \equiv m_{sub} / M_{h}) = {(\frac{μ}{{\tilde{μ}}_{1}})}^{1 + α} \exp [- {(\frac{μ}{μ_{cut}})}^{b}],

(19)

where the exponential cutoff accounts for the increasing rarity of subhaloes of mass close to that of the host halo mass49. The parameter,

{\tilde{μ}}_{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 [232]. 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 normalization parameter,

{\tilde{μ}}_{1}

. This dependence is amplified, and perhaps better expressed, if the subhalo velocity function is used instead, i.e., if the abundance is given in terms of the maximum circular velocity instead of the mass (see left panel of Figure 13). In this case, an accurate approximation is given by [223]:

N (> x \equiv V_{sub} / V_{h}) \propto V_{h}^{1 / 2} x^{- 3}, x < 0.7,

(20)

which implies, for instance, that a cluster-sized halo (

V_{h} \sim 1000

km/s) has ∼2.2 times more substructure of a given velocity ratio than a Milky Way-size halo (

V_{h} \sim 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}

_{⊙}

) 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 \times 10^{4}

_{⊙}

) to date of a Milky Way-size halo [138]. 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 [138,225,231,234,235,236,237] 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 Figure 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 [138,175,231]. 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 [222] for an analytical model of the subhalo distribution). It is interesting that the radial distribution of subhaloes with maximum circular velocity >

V_{sub}

is steeper than that of subhaloes with mass >

m_{sub}

(see e.g., [238]), 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 normalized to the virial radius as defined in [224] and shown in the right-hand panel of Figure 13), is approximated quite accurately by the following functional form [224]:

ϕ (x \equiv r / r_{vir}) = \frac{d \tilde{N} / d x^{3} |_{sub}}{d \tilde{N} / d x^{3} |_{NFW}} = 4 \frac{x^{4}}{{(1 + x)}^{2}}

(21)

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 Figure 14), while the outer regions show a steep truncation at a radius approximately equal to the tidal radius given in Equation (13) (see Figure 15 of [138]). 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 [239]50:

ρ_{E} (x_{E} \equiv r / r_{- 2}) = ρ_{- 2} e^{- 2 (x_{E}^{α_{E}} - 1) / α_{E}},

(22)

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

α_{E} \sim 0.22

[241,242]. Although for isolated haloes the parameters

α_{E}

and

r_{- 2}

can be related to the virial mass of the halo,

M_{200}

[135,240], in a similar way as the halo (NFW) concentration is connected to the virial mass, the situation is less clear for subhaloes [241], 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 [138,226]:

δ_{V} = \frac{\bar{ρ} (< r_{\max})}{ρ_{crit}} = 2 {(\frac{V_{\max}}{H r_{\max}})}^{2},

(23)

where H is the Hubble constant. Equation (23) can be related to the standard scale density of the NFW profile (

δ_{c}

in Equation (9)), and thus to the NFW concentration, in a straightforward way [226]:

δ_{V} = {(\frac{c}{2.163})}^{3} \frac{K_{c} (2.163)}{K_{c} (c)} Δ,

(24)

where

K_{c}

was defined just after Equation (9). We note that for the NFW profile,

r_{\max} / r_{s} = 2.163

Since the concentration of haloes (and subhaloes) is tightly correlated with their mass (see Section 2.4), Equation (23) implies a tight correlation between

V_{\max}

and

r_{\max}

, which indeed has been found and characterized in simulations (see right panel of Figure 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 mass52:

\begin{matrix} V_{\max} & = & 10 km / s {(\frac{m_{sub}}{3.37 \times 10^{7} M_{⊙}})}^{0.29} \\ \frac{δ_{V} (z = 0)}{2} = {(\frac{V_{\max}}{H_{0} r_{\max}})}^{2} & = & 2.9 \times 10^{4} {(\frac{m_{sub}}{10^{8} M_{⊙}})}^{- 0.18} . \end{matrix}

(25)

The fitting function for the

r_{\max} - V_{\max}

relation implied by Equation (25) is shown as a dotted line in the upper right panel of Figure 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 Figure 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} \sim 2.6

larger for subhaloes than for haloes of the same

V_{\max}

(lower right panel of Figure 14), which roughly corresponds to a

30 %

increase in the NFW concentration. This relative increase in concentration is larger for subhaloes nearer the center 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 [210]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 center (e.g., [244] 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 equilibrium 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 [245]. 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 [246]. 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 center outwards [246], 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 behavior of CDM field haloes, but with a slightly shallower slope,

Q_{sub} \sim r^{- 1.6}

[246], 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 behavior 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 Equations (19) and (20), is broken [247] 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 [222,247].

In many SIDM models the subhalo mass function remains largely unchanged as long as the interaction cross-section,

σ_{T} / m_{χ} < 10

^{2} /

g [57,58,248]. 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 Figure 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 [247] 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 [248]. 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 [249].

4. 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 years55. This strong statement carries a couple of major caveats, which define today’s frontier in N-body simulations of cosmological structure formation.

First, 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 SIDM56, 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 SIDM57. 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 summarize these in Section 4.1 below.

Finally, we make two remarks concerning the limited resolution of current N-simulations: (i) there have been recent claims that subhaloes can be artificially disrupted in cosmological simulations due to discreteness effects and inadequate force resolution [208,209]; 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}

_{⊙}

, 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.

4.1. 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., [255]), 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: gas dynamics, 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,256]58. Subsequently, the gas can radiatively cool and condense towards the center 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}

_{⊙}

) below which cooling becomes highly inefficient. At higher masses (∼

10^{13}

_{⊙}

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,258]. A hot, quasi-hydrostatic corona forms from which gas can subsequently cool at the center. 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 modeled assuming an adiabatic response leading to the contraction of the halo [259,260]. Even though the assembly of baryonic matter by mergers is not, in general, adiabatic, the simple adiabatic model remains a reasonable approximation [261]. 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., [262,263]).

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 reionizes this gas and heats it up, suppressing gas cooling into low-mass haloes and subsequent star formation [264,265]. This heating mechanism moves the minimum scale for galaxy formation from the atomic cooling limit to larger halo masses of order

10^{9}

_{⊙}

today 59 [250,267,268,269,270,271,272,273,274]. 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 [250,275].

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 [256,276,277] and inner structure of low-luminosity galaxies. Acting in conjunction with reionization, such winds strongly suppress galaxy formation in small haloes, reducing the abundance of luminous low-mass galaxies [269,271,272].

Energy injection from supernova can potentially alter the inner structure of dark matter haloes: if gas becomes gravitationally dominant in the center 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 center to expand, turning the original NFW cusp into a core. This mechanism, first proposed in the 1990s [278], became fashionable again several years later [279,280,281,282,283,284,285,286] when tentative observational evidence for cores, particularly in dwarf galaxies, began to emerge [287]. This evidence, however, is controversial [288,289]. While the proof of concept in [278] 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 [282]. This core-formation mechanism depends on the details of the baryon physics implemented in the simulation [290] and not all cosmological simulations produce cores in dwarf galaxies [291]. 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 [292,293,294,295,296].

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 [297,298] 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 center [299,300,301].

There has been great progress in the past decade in incorporating baryonic physics into full cosmological simulations; today galaxy formation and evolution can be modeled in unprecedented detail [302,303,304,305,306]. 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. Despite 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 [304] for a discussion of the limitations of gas dynamics simulations).

4.2. 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 among 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 Center [307], the latter through a 3.5 kev decay line in the x-ray spectra of galaxies and clusters [308,309] 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}

_{⊙}

, are dark; (ii) the presence (in the absence of the baryon effects discussed in the preceding section) of a steep cusp (

ρ \propto 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 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 [312,313]. 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 [314]. 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., [315]). 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 [316,317,318].

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}

_{⊙}

was originally calculated using semi-analytic techniques [264,265,267], as were the implications for the abundance of galactic satellites in the CDM model [269,271,272]. This solution has been repeatedly confirmed by modern gas dynamic simulations (e.g., [250,266,319,320]). 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 [250]. 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 cores61. 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., [317,323]), once the statistics are calculated rigorously, taking into account the “look elsewhere effect” [324]62. The origin of these planes is almost certainly the anisotropic nature of the accretion of satellites along filaments of the cosmic web [326,327] although the exact mechanism is still unclear as is the expected frequency of these structures.

While it is now generally agreed among 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 length scale 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 [328] 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 conflict with the data if the Milky Way halo mass is smaller than

1.1 \times 10^{12}

_{⊙}

[329]. Using a similar approach, [330] 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 \times 10^{12}

_{⊙}

, a value that may already be in conflict with the most recent determinations of the Milky Way halo mass [331]. 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,332].

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

^{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_{χ} \sim 1

^{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,249,251,252] 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: most haloes in CDM (and in many alternative dark matter models) are dark, i.e., unaffected or almost unaffected by baryons. It is the existence of a vast population of such small-mass haloes (

m ≲ 5 \times 10^{9}

_{⊙}

) 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 \times 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}

_{⊙}

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 center 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}

_{⊙}

and

10^{7}

–

10^{8}

_{⊙}

respectively (although the smaller gap could also be due to the impact of a giant molecular cloud) [339,340].

However, 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 multiple-lensed quasars; the second are small distortions of Einstein rings and large arcs.

In a multiple-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 quadruple-lensed quasars but dark substructures alone are insufficient to explain the observed anomalies [345], implying that other effects such as inadequate lens modeling may be at work. With better modeling 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 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 center 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}

_{⊙}

, 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}

_{⊙}

[349] 63 and could detect haloes as small as ∼

10^{7}

_{⊙}

[353,354].

Detecting the small signal generated by individual projected haloes or subhaloes requires accurate modeling of the source and the lens (the “macro” model; (e.g., [348,355])) and sophisticated statistical techniques to analyze 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}

_{⊙}

, are completely dark: they have never been modified in any way by baryons. 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.

Funding

This research received funding from the Icelandic Research Fund (Rannís), Grant of Excellence number 173929, the European Research Council (ERC) Advanced Investigator grant DMIDAS (GA 786910) and the Science and Technology Facilities Council (STFC) [grant number ST/F001166/1, ST/I00162X/1, ST/P000541/1].

Acknowledgments

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

Conflicts of Interest

The authors declare no conflict of interest.

References

Planck Collaboration. Planck 2018 results. VI. Cosmological parameters. arXiv 2018, arXiv:1807.06209. [Google Scholar]
Walker, M.G.; Mateo, M.; Olszewski, E.W.; Peñarrubia, J.; Evans, N.W.; Gilmore, G. A Universal Mass Profile for Dwarf Spheroidal Galaxies? Astrophys. J. 2009, 704, 1274–1287. [Google Scholar] [CrossRef]
Łokas, E.L.; Mamon, G.A. Dark matter distribution in the Coma cluster from galaxy kinematics: Breaking the mass-anisotropy degeneracy. Mon. Not. R. Astron. Soc. 2003, 343, 401–412. [Google Scholar]
Springel, V.; White, S.D.M.; Jenkins, A.; Frenk, C.S.; Yoshida, N.; Gao, L.; Navarro, J.; Thacker, R.; Croton, D.; Helly, J.; et al. Simulations of the formation, evolution and clustering of galaxies and quasars. Nature 2005, 435, 629–636. [Google Scholar] [CrossRef] [PubMed]
Guth, A.H. Inflationary universe: A possible solution to the horizon and flatness problems. Phys. Rev. D 1981, 23, 347–356. [Google Scholar] [CrossRef] [Green Version]
Linde, A.D. A new inflationary universe scenario: A possible solution of the horizon, flatness, homogeneity, isotropy and primordial monopole problems. Phys. Lett. B 1982, 108, 389–393. [Google Scholar] [CrossRef]
Kahlhoefer, F. Review of LHC dark matter searches. Int. J. Mod. Phys. A 2017, 32, 1730006. [Google Scholar] [CrossRef]
Xenon Collaboration. Dark Matter Search Results from a One Ton-Year Exposure of XENON1T. Phys. Rev. Lett. 2018, 121, 111302. [Google Scholar] [CrossRef] [Green Version]
ADMX Collaboration. Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment. Phys. Rev. Lett. 2018, 120, 151301. [Google Scholar] [CrossRef] [Green Version]
Fermi-LAT Collaboration; DES Collaboration. Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-Lat. Astrophys. J. 2017, 834, 110. [Google Scholar] [CrossRef]
Horiuchi, S.; Humphrey, P.J.; Oñorbe, J.; Abazajian, K.N.; Kaplinghat, M.; Garrison-Kimmel, S. Sterile neutrino dark matter bounds from galaxies of the Local Group. Phys. Rev. D 2014, 89, 025017. [Google Scholar] [CrossRef]
Gondolo, P.; Gelmini, G. Cosmic abundances of stable particles: Improved analysis. Nucl. Phys. B 1991, 360, 145–179. [Google Scholar] [CrossRef]
Jungman, G.; Kamionkowski, M.; Griest, K. Supersymmetric dark matter. Phys. Rep. 1996, 267, 195–373. [Google Scholar] [CrossRef] [Green Version]
Davis, M.; Efstathiou, G.; Frenk, C.S.; White, S.D.M. The evolution of large-scale structure in a universe dominated by cold dark matter. Astrophys. J. 1985, 292, 371–394. [Google Scholar] [CrossRef]
Preskill, J.; Wise, M.B.; Wilczek, F. Cosmology of the invisible axion. Phys. Lett. B 1983, 120, 127–132. [Google Scholar] [CrossRef] [Green Version]
Boyarsky, A.; Drewes, M.; Lasserre, T.; Mertens, S.; Ruchayskiy, O. Sterile Neutrino Dark Matter. arXiv 2018, arXiv:1807.07938. [Google Scholar] [CrossRef]
Hui, L.; Ostriker, J.P.; Tremaine, S.; Witten, E. Ultralight scalars as cosmological dark matter. Phys. Rev. D 2017, 95, 043541. [Google Scholar] [CrossRef] [Green Version]
Tulin, S.; Yu, H.B. Dark matter self-interactions and small scale structure. Phys. Rep. 2018, 730, 1–57. [Google Scholar] [CrossRef] [Green Version]
Zurek, K.M. Asymmetric Dark Matter: Theories, signatures, and constraints. Phys. Rep. 2014, 537, 91–121. [Google Scholar] [CrossRef] [Green Version]
Buckley, M.R.; Zavala, J.; Cyr-Racine, F.Y.; Sigurdson, K.; Vogelsberger, M. Scattering, damping, and acoustic oscillations: Simulating the structure of dark matter halos with relativistic force carriers. Phys. Rev. D 2014, 90, 043524. [Google Scholar] [CrossRef]
Cyr-Racine, F.Y.; Sigurdson, K.; Zavala, J.; Bringmann, T.; Vogelsberger, M.; Pfrommer, C. ETHOS: An effective theory of structure formation: From dark particle physics to the matter distribution of the Universe. Phys. Rev. D 2016, 93, 123527. [Google Scholar] [CrossRef]
Meszaros, P. The behaviour of point masses in an expanding cosmological substratum. Astron. Astrophys. 1974, 37, 225–228. [Google Scholar]
Green, A.M.; Hofmann, S.; Schwarz, D.J. The first WIMPy halos. J. Cosm. Part. Phys. 2005, 8, 3. [Google Scholar] [CrossRef]
Bringmann, T. Particle models and the small-scale structure of dark matter. New J. Phys. 2009, 11, 105027. [Google Scholar] [CrossRef]
Bœhm, C.; Riazuelo, A.; Hansen, S.H.; Schaeffer, R. Interacting dark matter disguised as warm dark matter. Phys. Rev. D 2002, 66, 083505. [Google Scholar] [CrossRef]
Bœhm, C.; Schewtschenko, J.A.; Wilkinson, R.J.; Baugh, C.M.; Pascoli, S. Using the Milky Way satellites to study interactions between cold dark matter and radiation. Mon. Not. R. Astron. Soc. 2014, 445, L31–L35. [Google Scholar] [CrossRef] [Green Version]
Loeb, A.; Zaldarriaga, M. Small-scale power spectrum of cold dark matter. Phys. Rev. D 2005, 71, 103520. [Google Scholar] [CrossRef]
Viel, M.; Becker, G.D.; Bolton, J.S.; Haehnelt, M.G. Warm dark matter as a solution to the small scale crisis: New constraints from high redshift Lyman-α forest data. Phys. Rev. D 2013, 88, 043502. [Google Scholar] [CrossRef]
Schive, H.Y.; Chiueh, T.; Broadhurst, T.; Huang, K.W. Contrasting Galaxy Formation from Quantum Wave Dark Matter, ψDM, with ΛCDM, using Planck and Hubble Data. Astrophys. J. 2016, 818, 89. [Google Scholar] [CrossRef]
Cole, S.; Percival, W.J.; Peacock, J.A.; Norberg, P.; Baugh, C.M.; Frenk, C.S.; Baldry, I.; Bland-Hawthorn, J.; Bridges, T.; Cannon, R.; et al. The 2dF Galaxy Redshift Survey: Power-spectrum analysis of the final data set and cosmological implications. Mon. Not. R. Astron. Soc. 2005, 362, 505–534. [Google Scholar] [CrossRef]
Percival, W.J.; Nichol, R.C.; Eisenstein, D.J.; Frieman, J.A.; Fukugita, M.; Loveday, J.; Pope, A.C.; Schneider, D.P.; Szalay, A.S.; Tegmark, M.; et al. The Shape of the Sloan Digital Sky Survey Data Release 5 Galaxy Power Spectrum. Astrophys. J. 2007, 657, 645–663. [Google Scholar] [CrossRef] [Green Version]
Zavala, J.; Lovell, M.R.; Vogelsberger, M.; Burger, J.D. The diverse dark matter density at sub-kiloparsec scales in Milky Way satellites:implications for the nature of dark matter. arXiv 2019, arXiv:1904.09998. [Google Scholar]
Carroll, S.M.; Press, W.H.; Turner, E.L. The cosmological constant. Annu. Rev. Astron. Astrophys. 1992, 30, 499–542. [Google Scholar] [CrossRef]
Baumann, D.; Nicolis, A.; Senatore, L.; Zaldarriaga, M. Cosmological non-linearities as an effective fluid. J. Cosmol. Astropart. Phys. 2012, 7, 51. [Google Scholar] [CrossRef]
Carrasco, J.J.M.; Hertzberg, M.P.; Senatore, L. The effective field theory of cosmological large scale structures. J. High Energy Phys. 2012, 9, 82. [Google Scholar] [CrossRef]
Gunn, J.E.; Gott, J.R., III. On the Infall of Matter Into Clusters of Galaxies and Some Effects on Their Evolution. Astrophys. J. 1972, 176, 1. [Google Scholar] [CrossRef]
Sheth, R.K.; Mo, H.J.; Tormen, G. Ellipsoidal collapse and an improved model for the number and spatial distribution of dark matter haloes. Mon. Not. R. Astron. Soc. 2001, 323, 1–12. [Google Scholar] [CrossRef] [Green Version]
Cooray, A.; Sheth, R. Halo models of large scale structure. Phys. Rep. 2002, 372, 1–129. [Google Scholar] [CrossRef] [Green Version]
Davis, M.; Peebles, P.J.E. On the integration of the BBGKY equations for the development of strongly nonlinear clustering in an expanding universe. Astrophys. J. Suppl. 1977, 34, 425–450. [Google Scholar] [CrossRef]
Smith, R.E.; Peacock, J.A.; Jenkins, A.; White, S.D.M.; Frenk, C.S.; Pearce, F.R.; Thomas, P.A.; Efstathiou, G.; Couchman, H.M.P. Stable clustering, the halo model and non-linear cosmological power spectra. Mon. Not. R. Astron. Soc. 2003, 341, 1311–1332. [Google Scholar] [CrossRef] [Green Version]
Afshordi, N.; Mohayaee, R.; Bertschinger, E. Hierarchy in the phase space and dark matter astronomy. Phys. Rev. D 2010, 81, 101301. [Google Scholar] [CrossRef]
Zavala, J.; Afshordi, N. Clustering in the phase space of dark matter haloes—II. Stable clustering and dark matter annihilation. Mon. Not. R. Astron. Soc. 2014, 441, 1329–1339. [Google Scholar] [CrossRef]
Zavala, J.; Afshordi, N. Universal clustering of dark matter in phase space. Mon. Not. R. Astron. Soc. 2016, 457, 986–992. [Google Scholar] [CrossRef] [Green Version]
Hahn, O.; Abel, T.; Kaehler, R. A new approach to simulating collisionless dark matter fluids. Mon. Not. R. Astron. Soc. 2013, 434, 1171–1191. [Google Scholar] [CrossRef] [Green Version]
Angulo, R.E.; Hahn, O.; Abel, T. The warm dark matter halo mass function below the cut-off scale. Mon. Not. R. Astron. Soc. 2013, 434, 3337–3347. [Google Scholar] [CrossRef]
Dehnen, W.; Read, J.I. N-body simulations of gravitational dynamics. Eur. Phys. J. Plus 2011, 126, 55. [Google Scholar] [CrossRef]
Dehnen, W. Towards optimal softening in three-dimensional N-body codes—I. Minimizing the force error. Mon. Not. R. Astron. Soc. 2001, 324, 273–291. [Google Scholar] [CrossRef]
Klypin, A.A.; Shandarin, S.F. Three-dimensional numerical model of the formation of large-scale structure in the Universe. Mon. Not. R. Astron. Soc. 1983, 204, 891–907. [Google Scholar] [CrossRef]
Melott, A.L. Massive neutrinos in large-scale gravitational clustering. Astrophys. J. 1983, 264, 59–78. [Google Scholar] [CrossRef]
Frenk, C.S.; White, S.D.M.; Davis, M. Nonlinear evolution of large-scale structure in the universe. Astrophys. J. 1983, 271, 417–430. [Google Scholar] [CrossRef]
Hockney, R.W.; Eastwood, J.W. Computer Simulation Using Particles; CRC Press: Boca Raton, FL, USA, 1988. [Google Scholar]
Efstathiou, G.; Eastwood, J.W. On the clustering of particles in an expanding universe. Mon. Not. R. Astron. Soc. 1981, 194, 503–525. [Google Scholar] [CrossRef]
Barnes, J.; Hut, P. A hierarchical O(N log N) force-calculation algorithm. Nature 1986, 324, 446–449. [Google Scholar] [CrossRef]
Springel, V. The cosmological simulation code GADGET-2. Mon. Not. R. Astron. Soc. 2005, 364, 1105–1134. [Google Scholar] [CrossRef]
Kochanek, C.S.; White, M. A Quantitative Study of Interacting Dark Matter in Halos. Astrophys. J. 2000, 543, 514–520. [Google Scholar] [CrossRef] [Green Version]
Davé, R.; Spergel, D.N.; Steinhardt, P.J.; Wandelt, B.D. Halo Properties in Cosmological Simulations of Self-interacting Cold Dark Matter. Astrophys. J. 2001, 547, 574–589. [Google Scholar] [CrossRef]
Vogelsberger, M.; Zavala, J.; Loeb, A. Subhaloes in self-interacting galactic dark matter haloes. Mon. Not. R. Astron. Soc. 2012, 423, 3740–3752. [Google Scholar] [CrossRef] [Green Version]
Rocha, M.; Peter, A.H.G.; Bullock, J.S.; Kaplinghat, M.; Garrison-Kimmel, S.; Oñorbe, J.; Moustakas, L.A. Cosmological simulations with self-interacting dark matter—I. Constant-density cores and substructure. Mon. Not. R. Astron. Soc. 2013, 430, 81–104. [Google Scholar] [CrossRef]
Robertson, A.; Massey, R.; Eke, V. Cosmic particle colliders: Simulations of self-interacting dark matter with anisotropic scattering. Mon. Not. R. Astron. Soc. 2017, 467, 4719–4730. [Google Scholar] [CrossRef]
Lynden-Bell, D.; Eggleton, P.P. On the consequences of the gravothermal catastrophe. Mon. Not. R. Astron. Soc. 1980, 191, 483–498. [Google Scholar] [CrossRef] [Green Version]
Koda, J.; Shapiro, P.R. Gravothermal collapse of isolated self-interacting dark matter haloes: N-body simulation versus the fluid model. Mon. Not. R. Astron. Soc. 2011, 415, 1125–1137. [Google Scholar] [CrossRef] [Green Version]
Schive, H.Y.; Chiueh, T.; Broadhurst, T. Cosmic structure as the quantum interference of a coherent dark wave. Nat. Phys. 2014, 10, 496–499. [Google Scholar] [CrossRef]
Mocz, P.; Vogelsberger, M.; Robles, V.H.; Zavala, J.; Boylan-Kolchin, M.; Fialkov, A.; Hernquist, L. Galaxy formation with BECDM—I. Turbulence and relaxation of idealized haloes. Mon. Not. R. Astron. Soc. 2017, 471, 4559–4570. [Google Scholar] [CrossRef] [PubMed]
Efstathiou, G.; Davis, M.; White, S.D.M.; Frenk, C.S. Numerical techniques for large cosmological N-body simulations. Astrophys. J. Suppl. 1985, 57, 241–260. [Google Scholar] [CrossRef]
Hahn, O.; Abel, T. Multi-scale initial conditions for cosmological simulations. Mon. Not. R. Astron. Soc. 2011, 415, 2101–2121. [Google Scholar] [CrossRef] [Green Version]
Jenkins, A. A new way of setting the phases for cosmological multiscale Gaussian initial conditions. Mon. Not. R. Astron. Soc. 2013, 434, 2094–2120. [Google Scholar] [CrossRef] [Green Version]
Sirko, E. Initial Conditions to Cosmological N-Body Simulations, or, How to Run an Ensemble of Simulations. Astrophys. J. 2005, 634, 728–743. [Google Scholar] [CrossRef]
Mo, H.; van den Bosch, F.C.; White, S. Galaxy Formation and Evolution; Cambridge University Press: Cambridge, UK, 2010. [Google Scholar]
White, S.D.M. Formation and Evolution of Galaxies. In Cosmology and Large Scale Structure; Schaeffer, R., Silk, J., Spiro, M., Zinn-Justin, J., Eds.; Cambridge University Press: Cambridge, UK, 1996; p. 349. [Google Scholar]
Zel’dovich, Y.B. Gravitational instability: An approximate theory for large density perturbations. Astron. Astrophys. 1970, 5, 84–89. [Google Scholar]
Jenkins, A. Second-order Lagrangian perturbation theory initial conditions for resimulations. Mon. Not. R. Astron. Soc. 2010, 403, 1859–1872. [Google Scholar] [CrossRef] [Green Version]
Pontzen, A.; Roškar, R.; Stinson, G.S.; Woods, R.; Reed, D.M.; Coles, J.; Quinn, T.R. Pynbody: Astrophysics Simulation Analysis for Python. Astrophysics Source Code Library. ascl:1305.002. 2013. Available online: https://pynbody.github.io/pynbody/#acknowledging-pynbody-in-scientific-publications (accessed on 24 September 2019).
Smith, R.E.; Markovic, K. Testing the warm dark matter paradigm with large-scale structures. Phys. Rev. D 2011, 84, 063507. [Google Scholar] [CrossRef]
Wang, J.; White, S.D.M. Discreteness effects in simulations of hot/warm dark matter. Mon. Not. R. Astron. Soc. 2007, 380, 93–103. [Google Scholar] [CrossRef] [Green Version]
Lovell, M.R.; Frenk, C.S.; Eke, V.R.; Jenkins, A.; Gao, L.; Theuns, T. The properties of warm dark matter haloes. Mon. Not. R. Astron. Soc. 2014, 439, 300–317. [Google Scholar] [CrossRef] [Green Version]
Hobbs, A.; Read, J.I.; Agertz, O.; Iannuzzi, F.; Power, C. NOVel Adaptive softening for collisionless N-body simulations: Eliminating spurious haloes. Mon. Not. R. Astron. Soc. 2016, 458, 468–479. [Google Scholar] [CrossRef]
Boylan-Kolchin, M.; Springel, V.; White, S.D.M.; Jenkins, A.; Lemson, G. Resolving cosmic structure formation with the Millennium-II Simulation. Mon. Not. R. Astron. Soc. 2009, 398, 1150–1164. [Google Scholar] [CrossRef] [Green Version]
Frenk, C.S.; White, S.D.M.; Davis, M.; Efstathiou, G. The formation of dark halos in a universe dominated by cold dark matter. Astrophys. J. 1988, 327, 507–525. [Google Scholar] [CrossRef]
Kuhlen, M.; Vogelsberger, M.; Angulo, R. Numerical simulations of the dark universe: State of the art and the next decade. Phys. Dark Univ. 2012, 1, 50–93. [Google Scholar] [CrossRef] [Green Version]
Potter, D.; Stadel, J.; Teyssier, R. PKDGRAV3: beyond trillion particle cosmological simulations for the next era of galaxy surveys. Comput. Astrophys. Cosmol. 2017, 4, 2. [Google Scholar] [CrossRef]
Geller, M.J.; Huchra, J.P. Mapping the universe. Science 1989, 246, 897–903. [Google Scholar] [CrossRef]
Gott, J.R., III; Jurić, M.; Schlegel, D.; Hoyle, F.; Vogeley, M.; Tegmark, M.; Bahcall, N.; Brinkmann, J. A Map of the Universe. Astrophys. J. 2005, 624, 463–484. [Google Scholar] [CrossRef]
Colless, M.; Dalton, G.; Maddox, S.; Sutherland, W.; Norberg, P.; Cole, S.; Bland-Hawthorn, J.; Bridges, T.; Cannon, R.; Collins, C.; et al. The 2dF Galaxy Redshift Survey: Spectra and redshifts. Mon. Not. R. Astron. Soc. 2001, 328, 1039–1063. [Google Scholar] [CrossRef]
Croton, D.J.; Springel, V.; White, S.D.M.; De Lucia, G.; Frenk, C.S.; Gao, L.; Jenkins, A.; Kauffmann, G.; Navarro, J.F.; Yoshida, N. The many lives of active galactic nuclei: Cooling flows, black holes and the luminosities and colours of galaxies. Mon. Not. R. Astron. Soc. 2006, 365, 11–28. [Google Scholar] [CrossRef]
Springel, V.; Frenk, C.S.; White, S.D.M. The large-scale structure of the Universe. Nature 2006, 440, 1137–1144. [Google Scholar] [CrossRef] [PubMed]
Bardeen, J.M.; Bond, J.R.; Kaiser, N.; Szalay, A.S. The statistics of peaks of Gaussian random fields. Astrophys. J. 1986, 304, 15–61. [Google Scholar] [CrossRef]
White, S.D.M.; Frenk, C.S. Galaxy formation through hierarchical clustering. Astrophys. J. 1991, 379, 52–79. [Google Scholar] [CrossRef]
Kauffmann, G.; Nusser, A.; Steinmetz, M. Galaxy formation and large-scale bias. Mon. Not. R. Astron. Soc. 1997, 286, 795–811. [Google Scholar] [CrossRef] [Green Version]
Kauffmann, G.; Colberg, J.M.; Diaferio, A.; White, S.D.M. Clustering of galaxies in a hierarchical universe—I. Methods and results at z = 0. Mon. Not. R. Astron. Soc. 1999, 303, 188–206. [Google Scholar] [CrossRef]
Benson, A.J.; Cole, S.; Frenk, C.S.; Baugh, C.M.; Lacey, C.G. The nature of galaxy bias and clustering. Mon. Not. R. Astron. Soc. 2000, 311, 793–808. [Google Scholar] [CrossRef] [Green Version]
Springel, V.; White, S.D.M.; Tormen, G.; Kauffmann, G. Populating a cluster of galaxies—I. Results at [formmu2]z = 0. Mon. Not. R. Astron. Soc. 2001, 328, 726–750. [Google Scholar] [CrossRef]
Lacey, C.G.; Baugh, C.M.; Frenk, C.S.; Benson, A.J.; Bower, R.G.; Cole, S.; Gonzalez-Perez, V.; Helly, J.C.; Lagos, C.D.P.; Mitchell, P.D. A unified multiwavelength model of galaxy formation. Mon. Not. R. Astron. Soc. 2016, 462, 3854–3911. [Google Scholar] [CrossRef]
White, S.D.M.; Frenk, C.S.; Davis, M. Clustering in a neutrino-dominated universe. Astrophys. J. Lett. 1983, 274, L1–L5. [Google Scholar] [CrossRef]
Angulo, R.E.; White, S.D.M. The birth and growth of neutralino haloes. Mon. Not. R. Astron. Soc. 2010, 401, 1796–1803. [Google Scholar] [CrossRef] [Green Version]
Stücker, J.; Busch, P.; White, S.D.M. The median density of the Universe. Mon. Not. R. Astron. Soc. 2018, 477, 3230–3246. [Google Scholar] [CrossRef] [Green Version]
White, M. The mass of a halo. Astron. Astrophys. 2001, 367, 27–32. [Google Scholar] [CrossRef] [Green Version]
Cuesta, A.J.; Prada, F.; Klypin, A.; Moles, M. The virialized mass of dark matter haloes. Mon. Not. R. Astron. Soc. 2008, 389, 385–397. [Google Scholar] [CrossRef] [Green Version]
Cole, S.; Lacey, C. The structure of dark matter haloes in hierarchical clustering models. Mon. Not. R. Astron. Soc. 1996, 281, 716. [Google Scholar] [CrossRef]
Eke, V.R.; Cole, S.; Frenk, C.S. Cluster evolution as a diagnostic for Omega. Mon. Not. R. Astron. Soc. 1996, 282, 263–280. [Google Scholar] [CrossRef]
Bryan, G.L.; Norman, M.L. Statistical Properties of X-Ray Clusters: Analytic and Numerical Comparisons. Astrophys. J. 1998, 495, 80–99. [Google Scholar] [CrossRef]
Jenkins, A.; Frenk, C.S.; White, S.D.M.; Colberg, J.M.; Cole, S.; Evrard, A.E.; Couchman, H.M.P.; Yoshida, N. The mass function of dark matter haloes. Mon. Not. R. Astron. Soc. 2001, 321, 372–384. [Google Scholar] [CrossRef] [Green Version]
Warren, M.S.; Abazajian, K.; Holz, D.E.; Teodoro, L. Precision Determination of the Mass Function of Dark Matter Halos. Astrophys. J. 2006, 646, 881–885. [Google Scholar] [CrossRef] [Green Version]
Lukić, Z.; Heitmann, K.; Habib, S.; Bashinsky, S.; Ricker, P.M. The Halo Mass Function: High-Redshift Evolution and Universality. Astrophys. J. 2007, 671, 1160–1181. [Google Scholar] [CrossRef]
Tinker, J.; Kravtsov, A.V.; Klypin, A.; Abazajian, K.; Warren, M.; Yepes, G.; Gottlöber, S.; Holz, D.E. Toward a Halo Mass Function for Precision Cosmology: The Limits of Universality. Astrophys. J. 2008, 688, 709–728. [Google Scholar] [CrossRef]
Trujillo-Gomez, S.; Klypin, A.; Primack, J.; Romanowsky, A.J. Galaxies in ΛCDM with Halo Abundance Matching: Luminosity-Velocity Relation, Baryonic Mass-Velocity Relation, Velocity Function, and Clustering. Astrophys. J. 2011, 742, 16. [Google Scholar] [CrossRef]
Hellwing, W.A.; Frenk, C.S.; Cautun, M.; Bose, S.; Helly, J.; Jenkins, A.; Sawala, T.; Cytowski, M. The Copernicus Complexio: a high-resolution view of the small-scale Universe. Mon. Not. R. Astron. Soc. 2016, 457, 3492–3509. [Google Scholar] [CrossRef] [Green Version]
Crain, R.A.; Theuns, T.; Dalla Vecchia, C.; Eke, V.R.; Frenk, C.S.; Jenkins, A.; Kay, S.T.; Peacock, J.A.; Pearce, F.R.; Schaye, J.; et al. Galaxies-intergalactic medium interaction calculation—I. Galaxy formation as a function of large-scale environment. Mon. Not. R. Astron. Soc. 2009, 399, 1773–1794. [Google Scholar] [CrossRef]
Press, W.H.; Schechter, P. Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation. Astrophys. J. 1974, 187, 425–438. [Google Scholar] [CrossRef]
Bond, J.R.; Cole, S.; Efstathiou, G.; Kaiser, N. Excursion set mass functions for hierarchical Gaussian fluctuations. Astrophys. J. 1991, 379, 440–460. [Google Scholar] [CrossRef]
Bower, R.G. The evolution of groups of galaxies in the Press-Schechter formalism. Mon. Not. R. Astron. Soc. 1991, 248, 332–352. [Google Scholar] [CrossRef]
Schneider, A.; Smith, R.E.; Reed, D. Halo mass function and the free streaming scale. Mon. Not. R. Astron. Soc. 2013, 433, 1573–1587. [Google Scholar] [CrossRef] [Green Version]
Bose, S.; Hellwing, W.A.; Frenk, C.S.; Jenkins, A.; Lovell, M.R.; Helly, J.C.; Li, B. The Copernicus Complexio: Statistical properties of warm dark matter haloes. Mon. Not. R. Astron. Soc. 2016, 455, 318–333. [Google Scholar] [CrossRef]
Schewtschenko, J.A.; Wilkinson, R.J.; Baugh, C.M.; Bœhm, C.; Pascoli, S. Dark matter-radiation interactions: The impact on dark matter haloes. Mon. Not. R. Astron. Soc. 2015, 449, 3587–3596. [Google Scholar] [CrossRef]
Vogelsberger, M.; Zavala, J.; Cyr-Racine, F.Y.; Pfrommer, C.; Bringmann, T.; Sigurdson, K. ETHOS—An effective theory of structure formation: Dark matter physics as a possible explanation of the small-scale CDM problems. Mon. Not. R. Astron. Soc. 2016, 460, 1399–1416. [Google Scholar] [CrossRef]
Benson, A.J.; Farahi, A.; Cole, S.; Moustakas, L.A.; Jenkins, A.; Lovell, M.; Kennedy, R.; Helly, J.; Frenk, C. Dark matter halo merger histories beyond cold dark matter—I. Methods and application to warm dark matter. Mon. Not. R. Astron. Soc. 2013, 428, 1774–1789. [Google Scholar] [CrossRef]
Leo, M.; Baugh, C.M.; Li, B.; Pascoli, S. A new smooth-k space filter approach to calculate halo abundances. J. Cosmol. Astropart. Phys. 2018, 4, 010. [Google Scholar] [CrossRef]
Sameie, O.; Benson, A.J.; Sales, L.V.; Yu, H.B.; Moustakas, L.A.; Creasey, P. The effect of dark matter-dark radiation interactions on halo abundance—A Press-Schechter approach. arXiv 2018, arXiv:1810.11040. [Google Scholar] [CrossRef]
Cyr-Racine, F.Y.; Sigurdson, K. Cosmology of atomic dark matter. Phys. Rev. D 2013, 87, 103515. [Google Scholar] [CrossRef]
Navarro, J.F.; Frenk, C.S.; White, S.D.M. The Structure of Cold Dark Matter Halos. Astrophys. J. 1996, 462, 563. [Google Scholar] [CrossRef] [Green Version]
Navarro, J.F.; Frenk, C.S.; White, S.D.M. A Universal Density Profile from Hierarchical Clustering. Astrophys. J. 1997, 490, 493–508. [Google Scholar] [CrossRef]
Navarro, J.F.; Ludlow, A.; Springel, V.; Wang, J.; Vogelsberger, M.; White, S.D.M.; Jenkins, A.; Frenk, C.S.; Helmi, A. The diversity and similarity of simulated cold dark matter haloes. Mon. Not. R. Astron. Soc. 2010, 402, 21–34. [Google Scholar] [CrossRef]
Anderhalden, D.; Diemand, J. Density profiles of CDM microhalos and their implications for annihilation boost factors. J. Cosmol. Astropart. Phys. 2013, 4, 009. [Google Scholar] [CrossRef]
Ishiyama, T. Hierarchical Formation of Dark Matter Halos and the Free Streaming Scale. Astrophys. J. 2014, 788, 27. [Google Scholar] [CrossRef]
Angulo, R.E.; Hahn, O.; Ludlow, A.D.; Bonoli, S. Earth-mass haloes and the emergence of NFW density profiles. Mon. Not. R. Astron. Soc. 2017, 471, 4687–4701. [Google Scholar] [CrossRef] [Green Version]
Delos, M.S.; Erickcek, A.L.; Bailey, A.P.; Alvarez, M.A. Density profiles of ultracompact minihalos: Implications for constraining the primordial power spectrum. Phys. Rev. D 2018, 98, 063527. [Google Scholar] [CrossRef] [Green Version]
Bullock, J.S.; Kolatt, T.S.; Sigad, Y.; Somerville, R.S.; Kravtsov, A.V.; Klypin, A.A.; Primack, J.R.; Dekel, A. Profiles of dark haloes: Evolution, scatter and environment. Mon. Not. R. Astron. Soc. 2001, 321, 559–575. [Google Scholar] [CrossRef]
Eke, V.R.; Navarro, J.F.; Steinmetz, M. The Power Spectrum Dependence of Dark Matter Halo Concentrations. Astrophys. J. 2001, 554, 114–125. [Google Scholar] [CrossRef]
Wechsler, R.H.; Bullock, J.S.; Primack, J.R.; Kravtsov, A.V.; Dekel, A. Concentrations of Dark Halos from Their Assembly Histories. Astrophys. J. 2002, 568, 52–70. [Google Scholar] [CrossRef] [Green Version]
Neto, A.F.; Gao, L.; Bett, P.; Cole, S.; Navarro, J.F.; Frenk, C.S.; White, S.D.M.; Springel, V.; Jenkins, A. The statistics of Λ CDM halo concentrations. Mon. Not. R. Astron. Soc. 2007, 381, 1450–1462. [Google Scholar] [CrossRef]
Zhao, D.H.; Jing, Y.P.; Mo, H.J.; Börner, G. Accurate Universal Models for the Mass Accretion Histories and Concentrations of Dark Matter Halos. Astrophys. J. 2009, 707, 354–369. [Google Scholar] [CrossRef]
Prada, F.; Klypin, A.A.; Cuesta, A.J.; Betancort-Rijo, J.E.; Primack, J. Halo concentrations in the standard Λ cold dark matter cosmology. Mon. Not. R. Astron. Soc. 2012, 423, 3018–3030. [Google Scholar] [CrossRef]
Ludlow, A.D.; Navarro, J.F.; Angulo, R.E.; Boylan-Kolchin, M.; Springel, V.; Frenk, C.; White, S.D.M. The mass-concentration-redshift relation of cold dark matter haloes. Mon. Not. R. Astron. Soc. 2014, 441, 378–388. [Google Scholar] [CrossRef]
Sánchez-Conde, M.A.; Prada, F. The flattening of the concentration-mass relation towards low halo masses and its implications for the annihilation signal boost. Mon. Not. R. Astron. Soc. 2014, 442, 2271–2277. [Google Scholar] [CrossRef]
Diemer, B.; Kravtsov, A.V. A Universal Model for Halo Concentrations. Astrophys. J. 2015, 799, 108. [Google Scholar] [CrossRef]
Klypin, A.; Yepes, G.; Gottlöber, S.; Prada, F.; Heß, S. MultiDark simulations: The story of dark matter halo concentrations and density profiles. Mon. Not. R. Astron. Soc. 2016, 457, 4340–4359. [Google Scholar] [CrossRef]
Pilipenko, S.V.; Sánchez-Conde, M.A.; Prada, F.; Yepes, G. Pushing down the low-mass halo concentration frontier with the Lomonosov cosmological simulations. Mon. Not. R. Astron. Soc. 2017, 472, 4918–4927. [Google Scholar] [CrossRef]
Wang, J.; Navarro, J.F.; Frenk, C.S.; White, S.D.M.; Springel, V.; Jenkins, A.; Helmi, A.; Ludlow, A.; Vogelsberger, M. Assembly history and structure of galactic cold dark matter haloes. Mon. Not. R. Astron. Soc. 2011, 413, 1373–1382. [Google Scholar] [CrossRef] [Green Version]
Springel, V.; Wang, J.; Vogelsberger, M.; Ludlow, A.; Jenkins, A.; Helmi, A.; Navarro, J.F.; Frenk, C.S.; White, S.D.M. The Aquarius Project: The subhaloes of galactic haloes. Mon. Not. R. Astron. Soc. 2008, 391, 1685–1711. [Google Scholar] [CrossRef]
Vera-Ciro, C.A.; Sales, L.V.; Helmi, A.; Frenk, C.S.; Navarro, J.F.; Springel, V.; Vogelsberger, M.; White, S.D.M. The shape of dark matter haloes in the Aquarius simulations: Evolution and memory. Mon. Not. R. Astron. Soc. 2011, 416, 1377–1391. [Google Scholar] [CrossRef]
Vogelsberger, M.; Helmi, A.; Springel, V.; White, S.D.M.; Wang, J.; Frenk, C.S.; Jenkins, A.; Ludlow, A.; Navarro, J.F. Phase-space structure in the local dark matter distribution and its signature in direct detection experiments. Mon. Not. R. Astron. Soc. 2009, 395, 797–811. [Google Scholar] [CrossRef]
Binney, J.; Tremaine, S. Galactic Dynamics: Second Edition; Princeton University Press: Princeton, NJ, USA, 2008. [Google Scholar]
Ludlow, A.D.; Navarro, J.F.; White, S.D.M.; Boylan-Kolchin, M.; Springel, V.; Jenkins, A.; Frenk, C.S. The density and pseudo-phase-space density profiles of cold dark matter haloes. Mon. Not. R. Astron. Soc. 2011, 415, 3895–3902. [Google Scholar] [CrossRef]
Hansen, S.H.; Moore, B. A universal density slope Velocity anisotropy relation for relaxed structures. New Astron. 2006, 11, 333–338. [Google Scholar] [CrossRef]
Taylor, J.E.; Navarro, J.F. The Phase-Space Density Profiles of Cold Dark Matter Halos. Astrophys. J. 2001, 563, 483–488. [Google Scholar] [CrossRef]
Bertschinger, E. Self-similar secondary infall and accretion in an Einstein-de Sitter universe. Astrophys. J. Suppl. 1985, 58, 39–65. [Google Scholar] [CrossRef]
Eddington, A.S. The distribution of stars in globular clusters. Mon. Not. R. Astron. Soc. 1916, 76, 572–585. [Google Scholar] [CrossRef]
Jing, Y.P.; Suto, Y. Triaxial Modeling of Halo Density Profiles with High-Resolution N-Body Simulations. Astrophys. J. 2002, 574, 538–553. [Google Scholar] [CrossRef]
Hayashi, E.; Navarro, J.F.; Springel, V. The shape of the gravitational potential in cold dark matter haloes. Mon. Not. R. Astron. Soc. 2007, 377, 50–62. [Google Scholar] [CrossRef] [Green Version]
Ganeshaiah Veena, P.; Cautun, M.; van de Weygaert, R.; Tempel, E.; Jones, B.J.T.; Rieder, S.; Frenk, C.S. The Cosmic Ballet: Spin and shape alignments of haloes in the cosmic web. Mon. Not. R. Astron. Soc. 2018, 481, 414–438. [Google Scholar] [CrossRef]
Bonamigo, M.; Despali, G.; Limousin, M.; Angulo, R.; Giocoli, C.; Soucail, G. Universality of dark matter haloes shape over six decades in mass: Insights from the Millennium XXL and SBARBINE simulations. Mon. Not. R. Astron. Soc. 2015, 449, 3171–3182. [Google Scholar] [CrossRef]
Vega-Ferrero, J.; Yepes, G.; Gottlöber, S. On the shape of dark matter haloes from MultiDark Planck simulations. Mon. Not. R. Astron. Soc. 2017, 467, 3226–3238. [Google Scholar] [CrossRef] [Green Version]
Despali, G.; Giocoli, C.; Tormen, G. Some like it triaxial: The universality of dark matter halo shapes and their evolution along the cosmic time. Mon. Not. R. Astron. Soc. 2014, 443, 3208–3217. [Google Scholar] [CrossRef]
Lovell, M.R.; Eke, V.; Frenk, C.S.; Gao, L.; Jenkins, A.; Theuns, T.; Wang, J.; White, S.D.M.; Boyarsky, A.; Ruchayskiy, O. The haloes of bright satellite galaxies in a warm dark matter universe. Mon. Not. R. Astron. Soc. 2012, 420, 2318–2324. [Google Scholar] [CrossRef] [Green Version]
Colín, P.; Avila-Reese, V.; Valenzuela, O. Substructure and Halo Density Profiles in a Warm Dark Matter Cosmology. Astrophys. J. 2000, 542, 622–630. [Google Scholar] [CrossRef] [Green Version]
Avila-Reese, V.; Colín, P.; Valenzuela, O.; D’Onghia, E.; Firmani, C. Formation and Structure of Halos in a Warm Dark Matter Cosmology. Astrophys. J. 2001, 559, 516–530. [Google Scholar] [CrossRef] [Green Version]
Colín, P.; Valenzuela, O.; Avila-Reese, V. On the Structure of Dark Matter Halos at the Damping Scale of the Power Spectrum with and without Relict Velocities. Astrophys. J. 2008, 673, 203–214. [Google Scholar] [CrossRef] [Green Version]
Schneider, A.; Smith, R.E.; Macciò, A.V.; Moore, B. Non-linear evolution of cosmological structures in warm dark matter models. Mon. Not. R. Astron. Soc. 2012, 424, 684–698. [Google Scholar] [CrossRef] [Green Version]
Ludlow, A.D.; Bose, S.; Angulo, R.E.; Wang, L.; Hellwing, W.A.; Navarro, J.F.; Cole, S.; Frenk, C.S. The mass-concentration-redshift relation of cold and warm dark matter haloes. Mon. Not. R. Astron. Soc. 2016, 460, 1214–1232. [Google Scholar] [CrossRef]
Dalcanton, J.J.; Hogan, C.J. Halo Cores and Phase-Space Densities: Observational Constraints on Dark Matter Physics and Structure Formation. Astrophys. J. 2001, 561, 35–45. [Google Scholar] [CrossRef]
Macciò, A.V.; Paduroiu, S.; Anderhalden, D.; Schneider, A.; Moore, B. Cores in warm dark matter haloes: A Catch 22 problem. Mon. Not. R. Astron. Soc. 2012, 424, 1105–1112. [Google Scholar] [CrossRef]
Shao, S.; Gao, L.; Theuns, T.; Frenk, C.S. The phase-space density of fermionic dark matter haloes. Mon. Not. R. Astron. Soc. 2013, 430, 2346–2357. [Google Scholar] [CrossRef]
Colín, P.; Avila-Reese, V.; Valenzuela, O.; Firmani, C. Structure and Subhalo Population of Halos in a Self-interacting Dark Matter Cosmology. Astrophys. J. 2002, 581, 777–793. [Google Scholar] [CrossRef]
Balberg, S.; Shapiro, S.L.; Inagaki, S. Self-Interacting Dark Matter Halos and the Gravothermal Catastrophe. Astrophys. J. 2002, 568, 475–487. [Google Scholar] [CrossRef]
Lynden-Bell, D.; Wood, R. The gravo-thermal catastrophe in isothermal spheres and the onset of red-giant structure for stellar systems. Mon. Not. R. Astron. Soc. 1968, 138, 495. [Google Scholar] [CrossRef]
Pollack, J.; Spergel, D.N.; Steinhardt, P.J. Supermassive Black Holes from Ultra-strongly Self-interacting Dark Matter. Astrophys. J. 2015, 804, 131. [Google Scholar] [CrossRef]
Vogelsberger, M.; Zavala, J.; Schutz, K.; Slatyer, T.R. Evaporating the Milky Way halo and its satellites with inelastic self-interacting dark matter. Mon. Not. R. Astron. Soc. 2019, 484, 5437–5452. [Google Scholar] [CrossRef]
Yoshida, N.; Springel, V.; White, S.D.M.; Tormen, G. Weakly Self-interacting Dark Matter and the Structure of Dark Halos. Astrophys. J. Lett. 2000, 544, L87–L90. [Google Scholar] [CrossRef]
Zavala, J.; Vogelsberger, M.; Walker, M.G. Constraining self-interacting dark matter with the Milky Way’s dwarf spheroidals. Mon. Not. R. Astron. Soc. 2013, 431, L20–L24. [Google Scholar] [CrossRef]
Brinckmann, T.; Zavala, J.; Rapetti, D.; Hansen, S.H.; Vogelsberger, M. The structure and assembly history of cluster-sized haloes in self-interacting dark matter. Mon. Not. R. Astron. Soc. 2018, 474, 746–759. [Google Scholar] [CrossRef]
Peter, A.H.G.; Rocha, M.; Bullock, J.S.; Kaplinghat, M. Cosmological simulations with self-interacting dark matter—II. Halo shapes versus observations. Mon. Not. R. Astron. Soc. 2013, 430, 105–120. [Google Scholar] [CrossRef]
Vogelsberger, M.; Zavala, J. Direct detection of self-interacting dark matter. Mon. Not. R. Astron. Soc. 2013, 430, 1722–1735. [Google Scholar] [CrossRef] [Green Version]
Genel, S.; Bouché, N.; Naab, T.; Sternberg, A.; Genzel, R. The Growth of Dark Matter Halos: Evidence for Significant Smooth Accretion. Astrophys. J. 2010, 719, 229–239. [Google Scholar] [CrossRef]
Gill, S.P.D.; Knebe, A.; Gibson, B.K. The evolution of substructure—III. The outskirts of clusters. Mon. Not. R. Astron. Soc. 2005, 356, 1327–1332. [Google Scholar] [CrossRef]
Sales, L.V.; Navarro, J.F.; Abadi, M.G.; Steinmetz, M. Cosmic ménage à trois: The origin of satellite galaxies on extreme orbits. Mon. Not. R. Astron. Soc. 2007, 379, 1475–1483. [Google Scholar] [CrossRef]
Ludlow, A.D.; Navarro, J.F.; Springel, V.; Jenkins, A.; Frenk, C.S.; Helmi, A. The Unorthodox Orbits of Substructure Halos. Astrophys. J. 2009, 692, 931–941. [Google Scholar] [CrossRef]
Giocoli, C.; Tormen, G.; Sheth, R.K.; van den Bosch, F.C. The substructure hierarchy in dark matter haloes. Mon. Not. R. Astron. Soc. 2010, 404, 502–517. [Google Scholar] [CrossRef]
Srisawat, C.; Knebe, A.; Pearce, F.R.; Schneider, A.; Thomas, P.A.; Behroozi, P.; Dolag, K.; Elahi, P.J.; Han, J.; Helly, J.; et al. Sussing Merger Trees: The Merger Trees Comparison Project. Mon. Not. R. Astron. Soc. 2013, 436, 150–162. [Google Scholar] [CrossRef] [Green Version]
Fakhouri, O.; Ma, C.P. The nearly universal merger rate of dark matter haloes in ΛCDM cosmology. Mon. Not. R. Astron. Soc. 2008, 386, 577–592. [Google Scholar] [CrossRef]
Fakhouri, O.; Ma, C.P.; Boylan-Kolchin, M. The merger rates and mass assembly histories of dark matter haloes in the two Millennium simulations. Mon. Not. R. Astron. Soc. 2010, 406, 2267–2278. [Google Scholar] [CrossRef] [Green Version]
Poole, G.B.; Mutch, S.J.; Croton, D.J.; Wyithe, S. Convergence properties of halo merger trees; halo and substructure merger rates across cosmic history. Mon. Not. R. Astron. Soc. 2017, 472, 3659–3682. [Google Scholar] [CrossRef]
Lacey, C.; Cole, S. Merger rates in hierarchical models of galaxy formation. Mon. Not. R. Astron. Soc. 1993, 262, 627–649. [Google Scholar] [CrossRef] [Green Version]
Parkinson, H.; Cole, S.; Helly, J. Generating dark matter halo merger trees. Mon. Not. R. Astron. Soc. 2008, 383, 557–564. [Google Scholar] [CrossRef]
Cole, S.; Lacey, C.G.; Baugh, C.M.; Frenk, C.S. Hierarchical galaxy formation. Mon. Not. R. Astron. Soc. 2000, 319, 168–204. [Google Scholar] [CrossRef]
Benson, A.J. Orbital parameters of infalling dark matter substructures. Mon. Not. R. Astron. Soc. 2005, 358, 551–562. [Google Scholar] [CrossRef] [Green Version]
Tormen, G. The rise and fall of satellites in galaxy clusters. Mon. Not. R. Astron. Soc. 1997, 290, 411–421. [Google Scholar] [CrossRef] [Green Version]
Jiang, L.; Cole, S.; Sawala, T.; Frenk, C.S. Orbital parameters of infalling satellite haloes in the hierarchical ΛCDM model. Mon. Not. R. Astron. Soc. 2015, 448, 1674–1686. [Google Scholar] [CrossRef]
Wetzel, A.R. On the orbits of infalling satellite haloes. Mon. Not. R. Astron. Soc. 2011, 412, 49–58. [Google Scholar] [CrossRef]
Diemand, J.; Kuhlen, M.; Madau, P.; Zemp, M.; Moore, B.; Potter, D.; Stadel, J. Clumps and streams in the local dark matter distribution. Nature 2008, 454, 735–738. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Rocha, M.; Peter, A.H.G.; Bullock, J. Infall times for Milky Way satellites from their present-day kinematics. Mon. Not. R. Astron. Soc. 2012, 425, 231–244. [Google Scholar] [CrossRef] [Green Version]
Vogelsberger, M.; White, S.D.M. Streams and caustics: The fine-grained structure of Λ cold dark matter haloes. Mon. Not. R. Astron. Soc. 2011, 413, 1419–1438. [Google Scholar] [CrossRef]
Natarajan, A.; Sikivie, P. Inner caustics of cold dark matter halos. Phys. Rev. D 2006, 73, 023510. [Google Scholar] [CrossRef] [Green Version]
Vogelsberger, M.; White, S.D.M.; Mohayaee, R.; Springel, V. Caustics in growing cold dark matter haloes. Mon. Not. R. Astron. Soc. 2009, 400, 2174–2184. [Google Scholar] [CrossRef] [Green Version]
Onions, J.; Knebe, A.; Pearce, F.R.; Muldrew, S.I.; Lux, H.; Knollmann, S.R.; Ascasibar, Y.; Behroozi, P.; Elahi, P.; Han, J.; et al. Subhaloes going Notts: The subhalo-finder comparison project. Mon. Not. R. Astron. Soc. 2012, 423, 1200–1214. [Google Scholar] [CrossRef]
Klypin, A.; Gottlöber, S.; Kravtsov, A.V.; Khokhlov, A.M. Galaxies in N-Body Simulations: Overcoming the Overmerging Problem. Astrophys. J. 1999, 516, 530–551. [Google Scholar] [CrossRef]
Neyrinck, M.C.; Gnedin, N.Y.; Hamilton, A.J.S. VOBOZ: An almost-parameter-free halo-finding algorithm. Mon. Not. R. Astron. Soc. 2005, 356, 1222–1232. [Google Scholar] [CrossRef]
Knollmann, S.R.; Knebe, A. AHF: Amiga’s Halo Finder. Astrophys. J. Suppl. 2009, 182, 608–624. [Google Scholar] [CrossRef]
Han, J.; Jing, Y.P.; Wang, H.; Wang, W. Resolving subhaloes’ lives with the Hierarchical Bound-Tracing algorithm. Mon. Not. R. Astron. Soc. 2012, 427, 2437–2449. [Google Scholar] [CrossRef]
Han, J.; Cole, S.; Frenk, C.S.; Benitez-Llambay, A.; Helly, J. HBT+: An improved code for finding subhaloes and building merger trees in cosmological simulations. Mon. Not. R. Astron. Soc. 2018, 474, 604–617. [Google Scholar] [CrossRef]
Tormen, G.; Moscardini, L.; Yoshida, N. Properties of cluster satellites in hydrodynamical simulations. Mon. Not. R. Astron. Soc. 2004, 350, 1397–1408. [Google Scholar] [CrossRef] [Green Version]
Maciejewski, M.; Colombi, S.; Springel, V.; Alard, C.; Bouchet, F.R. Phase-space structures—II. Hierarchical Structure Finder. Mon. Not. R. Astron. Soc. 2009, 396, 1329–1348. [Google Scholar] [CrossRef]
Behroozi, P.S.; Wechsler, R.H.; Wu, H.Y. The ROCKSTAR Phase-space Temporal Halo Finder and the Velocity Offsets of Cluster Cores. Astrophys. J. 2013, 762, 109. [Google Scholar] [CrossRef]
Tormen, G.; Diaferio, A.; Syer, D. Survival of substructure within dark matter haloes. Mon. Not. R. Astron. Soc. 1998, 299, 728–742. [Google Scholar] [CrossRef]
Tollet, É.; Cattaneo, A.; Mamon, G.A.; Moutard, T.; van den Bosch, F.C. On stellar mass loss from galaxies in groups and clusters. Mon. Not. R. Astron. Soc. 2017, 471, 4170–4193. [Google Scholar] [CrossRef] [Green Version]
Taylor, J.E.; Babul, A. The Dynamics of Sinking Satellites around Disk Galaxies: A Poor Man’s Alternative to High-Resolution Numerical Simulations. Astrophys. J. 2001, 559, 716–735. [Google Scholar] [CrossRef]
Zentner, A.R.; Bullock, J.S. Halo Substructure and the Power Spectrum. Astrophys. J. 2003, 598, 49–72. [Google Scholar] [CrossRef] [Green Version]
Zentner, A.R.; Berlind, A.A.; Bullock, J.S.; Kravtsov, A.V.; Wechsler, R.H. The Physics of Galaxy Clustering. I. A Model for Subhalo Populations. Astrophys. J. 2005, 624, 505–525. [Google Scholar] [CrossRef] [Green Version]
Johnston, K.V. A Prescription for Building the Milky Way’s Halo from Disrupted Satellites. Astrophys. J. 1998, 495, 297–308. [Google Scholar] [CrossRef]
Van den Bosch, F.C.; Ogiya, G. Dark matter substructure in numerical simulations: A tale of discreteness noise, runaway instabilities, and artificial disruption. Mon. Not. R. Astron. Soc. 2018, 475, 4066–4087. [Google Scholar] [CrossRef]
Van den Bosch, F.C.; Ogiya, G.; Hahn, O.; Burkert, A. Disruption of dark matter substructure: Fact or fiction? Mon. Not. R. Astron. Soc. 2018, 474, 3043–3066. [Google Scholar] [CrossRef]
Hayashi, E.; Navarro, J.F.; Taylor, J.E.; Stadel, J.; Quinn, T. The Structural Evolution of Substructure. Astrophys. J. 2003, 584, 541–558. [Google Scholar] [CrossRef] [Green Version]
Spitzer, L., Jr. Disruption of Galactic Clusters. Astrophys. J. 1958, 127, 17. [Google Scholar] [CrossRef]
Gnedin, O.Y.; Hernquist, L.; Ostriker, J.P. Tidal Shocking by Extended Mass Distributions. Astrophys. J. 1999, 514, 109–118. [Google Scholar] [CrossRef] [Green Version]
Aguilar, L.A.; White, S.D.M. Tidal interactions between spherical galaxies. Astrophys. J. 1985, 295, 374. [Google Scholar] [CrossRef]
Aguilar, L.A.; White, S.D.M. The Density Profiles of Tidally Stripped Galaxies. Astrophys. J. 1986, 307, 97. [Google Scholar] [CrossRef]
Kazantzidis, S.; Mayer, L.; Mastropietro, C.; Diemand, J.; Stadel, J.; Moore, B. Density Profiles of Cold Dark Matter Substructure: Implications for the Missing-Satellites Problem. Astrophys. J. 2004, 608, 663–679. [Google Scholar] [CrossRef]
Moore, B.; Katz, N.; Lake, G.; Dressler, A.; Oemler, A. Galaxy harassment and the evolution of clusters of galaxies. Nature 1996, 379, 613–616. [Google Scholar] [CrossRef]
Ogiya, G.; Burkert, A. Dynamical friction and scratches of orbiting satellite galaxies on host systems. Mon. Not. R. Astron. Soc. 2016, 457, 2164–2172. [Google Scholar] [CrossRef]
Boylan-Kolchin, M.; Ma, C.P.; Quataert, E. Dynamical friction and galaxy merging time-scales. Mon. Not. R. Astron. Soc. 2008, 383, 93–101. [Google Scholar] [CrossRef] [Green Version]
Chandrasekhar, S. Dynamical Friction. I. General Considerations: The Coefficient of Dynamical Friction. Astrophys. J. 1943, 97, 255. [Google Scholar] [CrossRef]
Van den Bosch, F.C.; Lewis, G.F.; Lake, G.; Stadel, J. Substructure in Dark Halos: Orbital Eccentricities and Dynamical Friction. Astrophys. J. 1999, 515, 50–68. [Google Scholar] [CrossRef] [Green Version]
Hernquist, L. An analytical model for spherical galaxies and bulges. Astrophys. J. 1990, 356, 359–364. [Google Scholar] [CrossRef]
Han, J.; Cole, S.; Frenk, C.S.; Jing, Y. A unified model for the spatial and mass distribution of subhaloes. Mon. Not. R. Astron. Soc. 2016, 457, 1208–1223. [Google Scholar] [CrossRef] [Green Version]
Klypin, A.A.; Trujillo-Gomez, S.; Primack, J. Dark Matter Halos in the Standard Cosmological Model: Results from the Bolshoi Simulation. Astrophys. J. 2011, 740, 102. [Google Scholar] [CrossRef]
Jiang, F.; van den Bosch, F.C. Statistics of dark matter substructure—III. Halo-to-halo variance. Mon. Not. R. Astron. Soc. 2017, 472, 657–674. [Google Scholar] [CrossRef]
Gao, L.; White, S.D.M.; Jenkins, A.; Stoehr, F.; Springel, V. The subhalo populations of ΛCDM dark haloes. Mon. Not. R. Astron. Soc. 2004, 355, 819–834. [Google Scholar] [CrossRef]
Diemand, J.; Kuhlen, M.; Madau, P. Formation and Evolution of Galaxy Dark Matter Halos and Their Substructure. Astrophys. J. 2007, 667, 859–877. [Google Scholar] [CrossRef] [Green Version]
Gao, L.; Navarro, J.F.; Frenk, C.S.; Jenkins, A.; Springel, V.; White, S.D.M. The Phoenix Project: The dark side of rich Galaxy clusters. Mon. Not. R. Astron. Soc. 2012, 425, 2169–2186. [Google Scholar] [CrossRef]
Garrison-Kimmel, S.; Boylan-Kolchin, M.; Bullock, J.S.; Lee, K. ELVIS: Exploring the Local Volume in Simulations. Mon. Not. R. Astron. Soc. 2014, 438, 2578–2596. [Google Scholar] [CrossRef]
Griffen, B.F.; Ji, A.P.; Dooley, G.A.; Gómez, F.A.; Vogelsberger, M.; O’Shea, B.W.; Frebel, A. The Caterpillar Project: A Large Suite of Milky Way Sized Halos. Astrophys. J. 2016, 818, 10. [Google Scholar] [CrossRef]
Gao, L.; Frenk, C.S.; Boylan-Kolchin, M.; Jenkins, A.; Springel, V.; White, S.D.M. The statistics of the subhalo abundance of dark matter haloes. Mon. Not. R. Astron. Soc. 2011, 410, 2309–2314. [Google Scholar] [CrossRef]
Angulo, R.E.; Lacey, C.G.; Baugh, C.M.; Frenk, C.S. The fate of substructures in cold dark matter haloes. Mon. Not. R. Astron. Soc. 2009, 399, 983–995. [Google Scholar] [CrossRef] [Green Version]
Boylan-Kolchin, M.; Springel, V.; White, S.D.M.; Jenkins, A. There’s no place like home? Statistics of Milky Way-mass dark matter haloes. Mon. Not. R. Astron. Soc. 2010, 406, 896–912. [Google Scholar] [CrossRef]
Rodríguez-Puebla, A.; Behroozi, P.; Primack, J.; Klypin, A.; Lee, C.; Hellinger, D. Halo and subhalo demographics with Planck cosmological parameters: Bolshoi-Planck and MultiDark-Planck simulations. Mon. Not. R. Astron. Soc. 2016, 462, 893–916. [Google Scholar] [CrossRef]
Ghigna, S.; Moore, B.; Governato, F.; Lake, G.; Quinn, T.; Stadel, J. Density Profiles and Substructure of Dark Matter Halos: Converging Results at Ultra-High Numerical Resolution. Astrophys. J. 2000, 544, 616–628. [Google Scholar] [CrossRef]
Diemand, J.; Moore, B.; Stadel, J. Velocity and spatial biases in cold dark matter subhalo distributions. Mon. Not. R. Astron. Soc. 2004, 352, 535–546. [Google Scholar] [CrossRef] [Green Version]
Nagai, D.; Kravtsov, A.V. The Radial Distribution of Galaxies in Λ Cold Dark Matter Clusters. Astrophys. J. 2005, 618, 557–568. [Google Scholar] [CrossRef]
Diemand, J.; Kuhlen, M.; Madau, P. Dark Matter Substructure and Gamma-Ray Annihilation in the Milky Way Halo. Astrophys. J. 2007, 657, 262–270. [Google Scholar] [CrossRef]
Gao, L.; De Lucia, G.; White, S.D.M.; Jenkins, A. Galaxies and subhaloes in ΛCDM galaxy clusters. Mon. Not. R. Astron. Soc. 2004, 352, L1–L5. [Google Scholar] [CrossRef]
Navarro, J.F.; Hayashi, E.; Power, C.; Jenkins, A.R.; Frenk, C.S.; White, S.D.M.; Springel, V.; Stadel, J.; Quinn, T.R. The inner structure of ΛCDM haloes—III. Universality and asymptotic slopes. Mon. Not. R. Astron. Soc. 2004, 349, 1039–1051. [Google Scholar] [CrossRef]
Gao, L.; Navarro, J.F.; Cole, S.; Frenk, C.S.; White, S.D.M.; Springel, V.; Jenkins, A.; Neto, A.F. The redshift dependence of the structure of massive Λ cold dark matter haloes. Mon. Not. R. Astron. Soc. 2008, 387, 536–544. [Google Scholar] [CrossRef]
Vera-Ciro, C.A.; Helmi, A.; Starkenburg, E.; Breddels, M.A. Not too big, not too small: The dark haloes of the dwarf spheroidals in the Milky Way. Mon. Not. R. Astron. Soc. 2013, 428, 1696–1703. [Google Scholar] [CrossRef]
Dutton, A.A.; Macciò, A.V. Cold dark matter haloes in the Planck era: Evolution of structural parameters for Einasto and NFW profiles. Mon. Not. R. Astron. Soc. 2014, 441, 3359–3374. [Google Scholar] [CrossRef]
Peñarrubia, J.; Navarro, J.F.; McConnachie, A.W. The Tidal Evolution of Local Group Dwarf Spheroidals. Astrophys. J. 2008, 673, 226–240. [Google Scholar] [CrossRef] [Green Version]
Moliné, Á.; Sánchez-Conde, M.A.; Palomares-Ruiz, S.; Prada, F. Characterization of subhalo structural properties and implications for dark matter annihilation signals. Mon. Not. R. Astron. Soc. 2017, 466, 4974–4990. [Google Scholar] [CrossRef]
Barber, C.; Starkenburg, E.; Navarro, J.F.; McConnachie, A.W. Galactic tides and the shape and orientation of dwarf galaxy satellites. Mon. Not. R. Astron. Soc. 2015, 447, 1112–1125. [Google Scholar] [CrossRef]
Vera-Ciro, C.A.; Sales, L.V.; Helmi, A.; Navarro, J.F. The shape of dark matter subhaloes in the Aquarius simulations. Mon. Not. R. Astron. Soc. 2014, 439, 2863–2872. [Google Scholar] [CrossRef]
Bose, S.; Hellwing, W.A.; Frenk, C.S.; Jenkins, A.; Lovell, M.R.; Helly, J.C.; Li, B.; Gonzalez-Perez, V.; Gao, L. Substructure and galaxy formation in the Copernicus Complexio warm dark matter simulations. Mon. Not. R. Astron. Soc. 2017, 464, 4520–4533. [Google Scholar] [CrossRef]
Dooley, G.A.; Peter, A.H.G.; Vogelsberger, M.; Zavala, J.; Frebel, A. Enhanced tidal stripping of satellites in the galactic halo from dark matter self-interactions. Mon. Not. R. Astron. Soc. 2016, 461, 710–727. [Google Scholar] [CrossRef] [Green Version]
Nishikawa, H.; Boddy, K.K.; Kaplinghat, M. Accelerated core collapse in tidally stripped self-interacting dark matter halos. arXiv 2019, arXiv:1901.00499. [Google Scholar]
Sawala, T.; Frenk, C.S.; Fattahi, A.; Navarro, J.F.; Theuns, T.; Bower, R.G.; Crain, R.A.; Furlong, M.; Jenkins, A.; Schaller, M.; et al. The chosen few: The low-mass haloes that host faint galaxies. Mon. Not. R. Astron. Soc. 2016, 456, 85–97. [Google Scholar] [CrossRef]
Sameie, O.; Yu, H.B.; Sales, L.V.; Vogelsberger, M.; Zavala, J. Self-Interacting Dark Matter Subhalos in the Milky Way’s Tides. arXiv 2019, arXiv:1904.07872. [Google Scholar]
Kahlhoefer, F.; Kaplinghat, M.; Slatyer, T.R.; Wu, C.L. Diversity in density profiles of self-interacting dark matter satellite halos. arXiv 2019, arXiv:1904.10539. [Google Scholar]
Arkani-Hamed, N.; Finkbeiner, D.P.; Slatyer, T.R.; Weiner, N. A theory of dark matter. Phys. Rev. D 2009, 79, 015014. [Google Scholar] [CrossRef]
Todoroki, K.; Medvedev, M.V. Dark matter haloes in the multicomponent model—I. Substructure. Mon. Not. R. Astron. Soc. 2019, 483, 3983–4003. [Google Scholar] [CrossRef]
Eisenstein, D.J.; Hu, W. Baryonic Features in the Matter Transfer Function. Astrophys. J. 1998, 496, 605–614. [Google Scholar] [CrossRef]
White, S.D.M.; Rees, M.J. Core condensation in heavy halos — A two-stage theory for galaxy formation and clustering. Mon. Not. R. Astron. Soc. 1978, 183, 341–358. [Google Scholar] [CrossRef]
Tseliakhovich, D.; Hirata, C. Relative velocity of dark matter and baryonic fluids and the formation of the first structures. Phys. Rev. D 2010, 82, 083520. [Google Scholar] [CrossRef]
Bower, R.G.; Benson, A.J.; Malbon, R.; Helly, J.C.; Frenk, C.S.; Baugh, C.M.; Cole, S.; Lacey, C.G. Breaking the hierarchy of galaxy formation. Mon. Not. R. Astron. Soc. 2006, 370, 645–655. [Google Scholar] [CrossRef] [Green Version]
Blumenthal, G.R.; Faber, S.M.; Flores, R.; Primack, J.R. Contraction of dark matter galactic halos due to baryonic infall. Astrophys. J. 1986, 301, 27–34. [Google Scholar] [CrossRef]
Mo, H.J.; Mao, S.; White, S.D.M. The formation of galactic discs. Mon. Not. R. Astron. Soc. 1998, 295, 319–336. [Google Scholar] [CrossRef] [Green Version]
Gnedin, O.Y.; Kravtsov, A.V.; Klypin, A.A.; Nagai, D. Response of Dark Matter Halos to Condensation of Baryons: Cosmological Simulations and Improved Adiabatic Contraction Model. Astrophys. J. 2004, 616, 16–26. [Google Scholar] [CrossRef]
Schaller, M.; Frenk, C.S.; Bower, R.G.; Theuns, T.; Jenkins, A.; Schaye, J.; Crain, R.A.; Furlong, M.; Dalla Vecchia, C.; McCarthy, I.G. Baryon effects on the internal structure of ΛCDM haloes in the EAGLE simulations. Mon. Not. R. Astron. Soc. 2015, 451, 1247–1267. [Google Scholar] [CrossRef]
Lovell, M.R.; Pillepich, A.; Genel, S.; Nelson, D.; Springel, V.; Pakmor, R.; Marinacci, F.; Weinberger, R.; Torrey, P.; Vogelsberger, M.; et al. The fraction of dark matter within galaxies from the IllustrisTNG simulations. Mon. Not. R. Astron. Soc. 2018, 481, 1950–1975. [Google Scholar] [CrossRef] [Green Version]
Efstathiou, G. Suppressing the formation of dwarf galaxies via photoionization. Mon. Not. R. Astron. Soc. 1992, 256, 43P–47P. [Google Scholar] [CrossRef]
Babul, A.; Rees, M.J. On dwarf elliptical galaxies and the faint blue counts. Mon. Not. R. Astron. Soc. 1992, 255, 346–350. [Google Scholar] [CrossRef]
Okamoto, T.; Gao, L.; Theuns, T. Mass loss of galaxies due to an ultraviolet background. Mon. Not. R. Astron. Soc. 2008, 390, 920–928. [Google Scholar] [CrossRef]
Thoul, A.A.; Weinberg, D.H. Hydrodynamic Simulations of Galaxy Formation. II. Photoionization and the Formation of Low-Mass Galaxies. Astrophys. J. 1996, 465, 608. [Google Scholar] [CrossRef]
Barkana, R.; Loeb, A. The Photoevaporation of Dwarf Galaxies during Reionization. Astrophys. J. 1999, 523, 54–65. [Google Scholar] [CrossRef] [Green Version]
Bullock, J.S.; Kravtsov, A.V.; Weinberg, D.H. Reionization and the Abundance of Galactic Satellites. Astrophys. J. 2000, 539, 517–521. [Google Scholar] [CrossRef] [Green Version]
Gnedin, N.Y. Effect of Reionization on Structure Formation in the Universe. Astrophys. J. 2000, 542, 535–541. [Google Scholar] [CrossRef]
Benson, A.J.; Lacey, C.G.; Baugh, C.M.; Cole, S.; Frenk, C.S. The effects of photoionization on galaxy formation—I. Model and results at z = 0. Mon. Not. R. Astron. Soc. 2002, 333, 156–176. [Google Scholar] [CrossRef]
Somerville, R.S. Can Photoionization Squelching Resolve the Substructure Crisis? Astrophys. J. Lett. 2002, 572, L23–L26. [Google Scholar] [CrossRef]
Hoeft, M.; Yepes, G.; Gottlöber, S.; Springel, V. Dwarf galaxies in voids: Suppressing star formation with photoheating. Mon. Not. R. Astron. Soc. 2006, 371, 401–414. [Google Scholar] [CrossRef]
Ocvirk, P.; Gillet, N.; Shapiro, P.R.; Aubert, D.; Iliev, I.T.; Teyssier, R.; Yepes, G.; Choi, J.H.; Sullivan, D.; Knebe, A.; et al. Cosmic Dawn (CoDa): The First Radiation-Hydrodynamics Simulation of Reionization and Galaxy Formation in the Local Universe. Mon. Not. R. Astron. Soc. 2016, 463, 1462–1485. [Google Scholar] [CrossRef]
Sawala, T.; Frenk, C.S.; Crain, R.A.; Jenkins, A.; Schaye, J.; Theuns, T.; Zavala, J. The abundance of (not just) dark matter haloes. Mon. Not. R. Astron. Soc. 2013, 431, 1366–1382. [Google Scholar] [CrossRef]
Larson, R.B. Effects of supernovae on the early evolution of galaxies. Mon. Not. R. Astron. Soc. 1974, 169, 229–246. [Google Scholar] [CrossRef]
Dekel, A.; Silk, J. The origin of dwarf galaxies, cold dark matter, and biased galaxy formation. Astrophys. J. 1986, 303, 39–55. [Google Scholar] [CrossRef]
Navarro, J.F.; Eke, V.R.; Frenk, C.S. The cores of dwarf galaxy haloes. Mon. Not. R. Astron. Soc. 1996, 283, L72–L78. [Google Scholar] [CrossRef]
Read, J.I.; Gilmore, G. Mass loss from dwarf spheroidal galaxies: The origins of shallow dark matter cores and exponential surface brightness profiles. Mon. Not. R. Astron. Soc. 2005, 356, 107–124. [Google Scholar] [CrossRef]
Gnedin, O.Y.; Zhao, H. Maximum feedback and dark matter profiles of dwarf galaxies. Mon. Not. R. Astron. Soc. 2002, 333, 299–306. [Google Scholar] [CrossRef] [Green Version]
Governato, F.; Brook, C.; Mayer, L.; Brooks, A.; Rhee, G.; Wadsley, J.; Jonsson, P.; Willman, B.; Stinson, G.; Quinn, T.; et al. Bulgeless dwarf galaxies and dark matter cores from supernova-driven outflows. Nature 2010, 463, 203–206. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Pontzen, A.; Governato, F. How supernova feedback turns dark matter cusps into cores. Mon. Not. R. Astron. Soc. 2012, 421, 3464–3471. [Google Scholar] [CrossRef] [Green Version]
Di Cintio, A.; Brook, C.B.; Macciò, A.V.; Stinson, G.S.; Knebe, A.; Dutton, A.A.; Wadsley, J. The dependence of dark matter profiles on the stellar-to-halo mass ratio: A prediction for cusps versus cores. Mon. Not. R. Astron. Soc. 2014, 437, 415–423. [Google Scholar] [CrossRef]
Chan, T.K.; Kereš, D.; Oñorbe, J.; Hopkins, P.F.; Muratov, A.L.; Faucher-Giguère, C.A.; Quataert, E. The impact of baryonic physics on the structure of dark matter haloes: The view from the FIRE cosmological simulations. Mon. Not. R. Astron. Soc. 2015, 454, 2981–3001. [Google Scholar] [CrossRef]
Tollet, E.; Macciò, A.V.; Dutton, A.A.; Stinson, G.S.; Wang, L.; Penzo, C.; Gutcke, T.A.; Buck, T.; Kang, X.; Brook, C.; et al. NIHAO - IV: Core creation and destruction in dark matter density profiles across cosmic time. Mon. Not. R. Astron. Soc. 2016, 456, 3542–3552. [Google Scholar] [CrossRef]
Read, J.I.; Agertz, O.; Collins, M.L.M. Dark matter cores all the way down. Mon. Not. R. Astron. Soc. 2016, 459, 2573–2590. [Google Scholar] [CrossRef] [Green Version]
Moore, B. Evidence against dissipation-less dark matter from observations of galaxy haloes. Nature 1994, 370, 629–631. [Google Scholar] [CrossRef]
Oman, K.A.; Navarro, J.F.; Fattahi, A.; Frenk, C.S.; Sawala, T.; White, S.D.M.; Bower, R.; Crain, R.A.; Furlong, M.; Schaller, M.; et al. The unexpected diversity of dwarf galaxy rotation curves. Mon. Not. R. Astron. Soc. 2015, 452, 3650–3665. [Google Scholar] [CrossRef]
Oman, K.A.; Marasco, A.; Navarro, J.F.; Frenk, C.S.; Schaye, J.; Benítez-Llambay, A. Non-circular motions and the diversity of dwarf galaxy rotation curves. Mon. Not. R. Astron. Soc. 2019, 482, 821–847. [Google Scholar] [CrossRef]
Benitez-Llambay, A.; Frenk, C.S.; Ludlow, A.D.; Navarro, J.F. Baryon-induced dark matter cores in the EAGLE simulations. arXiv 2018, arXiv:1810.04186. [Google Scholar] [CrossRef]
Bose, S.; Frenk, C.S.; Jenkins, A.; Fattahi, A.; Gómez, F.A.; Grand, R.J.J.; Marinacci, F.; Navarro, J.F.; Oman, K.A.; Pakmor, R.; et al. No cores in dark matter-dominated dwarf galaxies with bursty star formation histories. Mon. Not. R. Astron. Soc. 2019, 486, 4790–4804. [Google Scholar] [CrossRef]
Peirani, S.; Kay, S.; Silk, J. Active galactic nuclei and massive galaxy cores. Astron. Astrophys. 2008, 479, 123–129. [Google Scholar] [CrossRef]
Duffy, A.R.; Schaye, J.; Kay, S.T.; Dalla Vecchia, C.; Battye, R.A.; Booth, C.M. Impact of baryon physics on dark matter structures: A detailed simulation study of halo density profiles. Mon. Not. R. Astron. Soc. 2010, 405, 2161–2178. [Google Scholar] [CrossRef]
Teyssier, R.; Moore, B.; Martizzi, D.; Dubois, Y.; Mayer, L. Mass distribution in galaxy clusters: The role of Active Galactic Nuclei feedback. Mon. Not. R. Astron. Soc. 2011, 414, 195–208. [Google Scholar] [CrossRef]
Martizzi, D.; Teyssier, R.; Moore, B. Cusp-core transformations induced by AGN feedback in the progenitors of cluster galaxies. Mon. Not. R. Astron. Soc. 2013, 432, 1947–1954. [Google Scholar] [CrossRef]
Peirani, S.; Dubois, Y.; Volonteri, M.; Devriendt, J.; Bundy, K.; Silk, J.; Pichon, C.; Kaviraj, S.; Gavazzi, R.; Habouzit, M. Density profile of dark matter haloes and galaxies in the HORIZON-AGN simulation: The impact of AGN feedback. Mon. Not. R. Astron. Soc. 2017, 472, 2153–2169. [Google Scholar] [CrossRef]
D’Onghia, E.; Springel, V.; Hernquist, L.; Keres, D. Substructure Depletion in the Milky Way Halo by the Disk. Astrophys. J. 2010, 709, 1138–1147. [Google Scholar] [CrossRef]
Kazantzidis, S.; Łokas, E.L.; Callegari, S.; Mayer, L.; Moustakas, L.A. On the Efficiency of the Tidal Stirring Mechanism for the Origin of Dwarf Spheroidals: Dependence on the Orbital and Structural Parameters of the Progenitor Disky Dwarfs. Astrophys. J. 2011, 726, 98. [Google Scholar] [CrossRef]
Zolotov, A.; Brooks, A.M.; Willman, B.; Governato, F.; Pontzen, A.; Christensen, C.; Dekel, A.; Quinn, T.; Shen, S.; Wadsley, J. Baryons Matter: Why Luminous Satellite Galaxies have Reduced Central Masses. Astrophys. J. 2012, 761, 71. [Google Scholar] [CrossRef]
Sawala, T.; Pihajoki, P.; Johansson, P.H.; Frenk, C.S.; Navarro, J.F.; Oman, K.A.; White, S.D.M. Shaken and stirred: The Milky Way’s dark substructures. Mon. Not. R. Astron. Soc. 2017, 467, 4383–4400. [Google Scholar] [CrossRef]
Garrison-Kimmel, S.; Hopkins, P.F.; Wetzel, A.; Bullock, J.S.; Boylan-Kolchin, M.; Keres, D.; Faucher-Giguere, C.A.; El-Badry, K.; Lamberts, A.; Quataert, E.; et al. The Local Group on FIRE: Dwarf galaxy populations across a suite of hydrodynamic simulations. arXiv 2018, arXiv:1806.04143. [Google Scholar] [CrossRef]
Dubois, Y.; Pichon, C.; Welker, C.; Le Borgne, D.; Devriendt, J.; Laigle, C.; Codis, S.; Pogosyan, D.; Arnouts, S.; Benabed, K.; et al. Dancing in the dark: Galactic properties trace spin swings along the cosmic web. Mon. Not. R. Astron. Soc. 2014, 444, 1453–1468. [Google Scholar] [CrossRef]
Vogelsberger, M.; Genel, S.; Springel, V.; Torrey, P.; Sijacki, D.; Xu, D.; Snyder, G.; Nelson, D.; Hernquist, L. Introducing the Illustris Project: Simulating the coevolution of dark and visible matter in the Universe. Mon. Not. R. Astron. Soc. 2014, 444, 1518–1547. [Google Scholar] [CrossRef]
Schaye, J.; Crain, R.A.; Bower, R.G.; Furlong, M.; Schaller, M.; Theuns, T.; Dalla Vecchia, C.; Frenk, C.S.; McCarthy, I.G.; Helly, J.C.; et al. The EAGLE project: Simulating the evolution and assembly of galaxies and their environments. Mon. Not. R. Astron. Soc. 2015, 446, 521–554. [Google Scholar] [CrossRef]
Khandai, N.; Di Matteo, T.; Croft, R.; Wilkins, S.; Feng, Y.; Tucker, E.; DeGraf, C.; Liu, M.S. The MassiveBlack-II simulation: The evolution of haloes and galaxies to z 0. Mon. Not. R. Astron. Soc. 2015, 450, 1349–1374. [Google Scholar] [CrossRef]
Pillepich, A.; Springel, V.; Nelson, D.; Genel, S.; Naiman, J.; Pakmor, R.; Hernquist, L.; Torrey, P.; Vogelsberger, M.; Weinberger, R.; et al. Simulating galaxy formation with the IllustrisTNG model. Mon. Not. R. Astron. Soc. 2018, 473, 4077–4106. [Google Scholar] [CrossRef]
Goodenough, L.; Hooper, D. Possible Evidence For Dark Matter Annihilation In The Inner Milky Way From The Fermi Gamma Ray Space Telescope. arXiv 2009, arXiv:0910.2998. [Google Scholar]
Bulbul, E.; Markevitch, M.; Foster, A.; Smith, R.K.; Loewenstein, M.; Randall, S.W. Detection of an Unidentified Emission Line in the Stacked X-Ray Spectrum of Galaxy Clusters. Astrophys. J. 2014, 789, 13. [Google Scholar] [CrossRef]
Boyarsky, A.; Ruchayskiy, O.; Iakubovskyi, D.; Franse, J. Unidentified Line in X-Ray Spectra of the Andromeda Galaxy and Perseus Galaxy Cluster. Phys. Rev. Lett. 2014, 113, 251301. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Ruchayskiy, O.; Boyarsky, A.; Iakubovskyi, D.; Bulbul, E.; Eckert, D.; Franse, J.; Malyshev, D.; Markevitch, M.; Neronov, A. Searching for decaying dark matter in deep XMM-Newton observation of the Draco dwarf spheroidal. Mon. Not. R. Astron. Soc. 2016, 460, 1390–1398. [Google Scholar] [CrossRef]
Aharonian, F.A.; Akamatsu, H.; Akimoto, F.; Allen, S.W.; Angelini, L.; Arnaud, K.A.; Audard, M.; Awaki, H.; Axelsson, M.; Bamba, A.; et al. Hitomi Constraints on the 3.5 keV Line in the Perseus Galaxy Cluster. Astrophys. J. Lett. 2017, 837, L15. [Google Scholar] [CrossRef]
Klypin, A.; Kravtsov, A.V.; Valenzuela, O.; Prada, F. Where Are the Missing Galactic Satellites? Astrophys. J. 1999, 522, 82–92. [Google Scholar] [CrossRef] [Green Version]
Moore, B.; Ghigna, S.; Governato, F.; Lake, G.; Quinn, T.; Stadel, J.; Tozzi, P. Dark Matter Substructure within Galactic Halos. Astrophys. J. 1999, 524, L19–L22. [Google Scholar] [CrossRef]
Boylan-Kolchin, M.; Bullock, J.S.; Kaplinghat, M. Too big to fail? The puzzling darkness of massive Milky Way subhaloes. Mon. Not. R. Astron. Soc. 2011, 415, L40–L44. [Google Scholar] [CrossRef] [Green Version]
Walker, M.G.; Peñarrubia, J. A Method for Measuring (Slopes of) the Mass Profiles of Dwarf Spheroidal Galaxies. Astrophys. J. 2011, 742, 20. [Google Scholar] [CrossRef]
Kroupa, P.; Theis, C.; Boily, C.M. The great disk of Milky-Way satellites and cosmological sub-structures. Astron. Astrophys. 2005, 431, 517–521. [Google Scholar] [CrossRef] [Green Version]
Ibata, R.A.; Ibata, N.G.; Lewis, G.F.; Martin, N.F.; Conn, A.; Elahi, P.; Arias, V.; Fernando, N. A Thousand Shadows of Andromeda: Rotating Planes of Satellites in the Millennium-II Cosmological Simulation. Astrophys. J. Lett. 2014, 784, L6. [Google Scholar] [CrossRef]
Pawlowski, M.S.; Famaey, B.; Jerjen, H.; Merritt, D.; Kroupa, P.; Dabringhausen, J.; Lüghausen, F.; Forbes, D.A.; Hensler, G.; Hammer, F.; et al. Co-orbiting satellite galaxy structures are still in conflict with the distribution of primordial dwarf galaxies. Mon. Not. R. Astron. Soc. 2014, 442, 2362–2380. [Google Scholar] [CrossRef] [Green Version]
Wadepuhl, M.; Springel, V. Satellite galaxies in hydrodynamical simulations of Milky Way sized galaxies. Mon. Not. R. Astron. Soc. 2011, 410, 1975–1992. [Google Scholar] [CrossRef]
Simpson, C.M.; Grand, R.J.J.; Gómez, F.A.; Marinacci, F.; Pakmor, R.; Springel, V.; Campbell, D.J.R.; Frenk, C.S. Quenching and ram pressure stripping of simulated Milky Way satellite galaxies. Mon. Not. R. Astron. Soc. 2018, 478, 548–567. [Google Scholar] [CrossRef] [Green Version]
Mashchenko, S.; Couchman, H.M.P.; Wadsley, J. The removal of cusps from galaxy centres by stellar feedback in the early Universe. Nature 2006, 442, 539–542. [Google Scholar] [CrossRef] [PubMed]
Weinberg, M.D.; Katz, N. Bar-driven Dark Halo Evolution: A Resolution of the Cusp-Core Controversy. Astrophys. J. 2002, 580, 627–633. [Google Scholar] [CrossRef]
Müller, O.; Jerjen, H.; Pawlowski, M.S.; Binggeli, B. Testing the two planes of satellites in the Centaurus group. Astron. Astrophys. 2016, 595, A119. [Google Scholar] [CrossRef]
Cautun, M.; Bose, S.; Frenk, C.S.; Guo, Q.; Han, J.; Hellwing, W.A.; Sawala, T.; Wang, W. Planes of satellite galaxies: When exceptions are the rule. Mon. Not. R. Astron. Soc. 2015, 452, 3838–3852. [Google Scholar] [CrossRef]
Müller, O.; Pawlowski, M.S.; Jerjen, H.; Lelli, F. A whirling plane of satellite galaxies around Centaurus A challenges cold dark matter cosmology. Science 2018, 359, 534–537. [Google Scholar] [CrossRef] [Green Version]
Libeskind, N.I.; Frenk, C.S.; Cole, S.; Helly, J.C.; Jenkins, A.; Navarro, J.F.; Power, C. The distribution of satellite galaxies: The great pancake. Mon. Not. R. Astron. Soc. 2005, 363, 146–152. [Google Scholar] [CrossRef]
Shao, S.; Cautun, M.; Frenk, C.S. Evolution of galactic planes of satellites in the EAGLE simulation. arXiv 2019, arXiv:1904.02719. [Google Scholar] [CrossRef]
Wang, J.; Frenk, C.S.; Navarro, J.F.; Gao, L.; Sawala, T. The missing massive satellites of the Milky Way. Mon. Not. R. Astron. Soc. 2012, 424, 2715–2721. [Google Scholar] [CrossRef] [Green Version]
Kennedy, R.; Frenk, C.; Cole, S.; Benson, A. Constraining the warm dark matter particle mass with Milky Way satellites. Mon. Not. R. Astron. Soc. 2014, 442, 2487–2495. [Google Scholar] [CrossRef]
Lovell, M.R.; Bose, S.; Boyarsky, A.; Cole, S.; Frenk, C.S.; Gonzalez-Perez, V.; Kennedy, R.; Ruchayskiy, O.; Smith, A. Satellite galaxies in semi-analytic models of galaxy formation with sterile neutrino dark matter. Mon. Not. R. Astron. Soc. 2016, 461, 60–72. [Google Scholar] [CrossRef]
Callingham, T.M.; Cautun, M.; Deason, A.J.; Frenk, C.S.; Wang, W.; Gómez, F.A.; Grand, R.J.J.; Marinacci, F.; Pakmor, R. The mass of the Milky Way from satellite dynamics. Mon. Not. R. Astron. Soc. 2019, 484, 5453–5467. [Google Scholar] [CrossRef]
Lovell, M.R.; Gonzalez-Perez, V.; Bose, S.; Boyarsky, A.; Cole, S.; Frenk, C.S.; Ruchayskiy, O. Addressing the too big to fail problem with baryon physics and sterile neutrino dark matter. Mon. Not. R. Astron. Soc. 2017, 468, 2836–2849. [Google Scholar] [CrossRef] [Green Version]
Burger, J.D.; Zavala, J. The nature of core formation in dark matter haloes: Adiabatic or impulsive? Mon. Not. R. Astron. Soc. 2019, 485, 1008–1028. [Google Scholar] [CrossRef]
Rees, M.J. Lyman absorption lines in quasar spectra—Evidence for gravitationally-confined gas in dark minihaloes. Mon. Not. R. Astron. Soc. 1986, 218, 25P–30P. [Google Scholar] [CrossRef]
Benítez-Llambay, A.; Navarro, J.F.; Frenk, C.S.; Sawala, T.; Oman, K.; Fattahi, A.; Schaller, M.; Schaye, J.; Crain, R.A.; Theuns, T. The properties of ‘dark’ ΛCDM haloes in the Local Group. Mon. Not. R. Astron. Soc. 2017, 465, 3913–3926. [Google Scholar] [CrossRef]
Carlberg, R.G. Dark Matter Sub-halo Counts via Star Stream Crossings. Astrophys. J. 2012, 748, 20. [Google Scholar] [CrossRef]
Erkal, D.; Belokurov, V.; Bovy, J.; Sand ers, J.L. The number and size of subhalo-induced gaps in stellar streams. Mon. Not. R. Astron. Soc. 2016, 463, 102–119. [Google Scholar] [CrossRef]
Amorisco, N.C.; Gómez, F.A.; Vegetti, S.; White, S.D.M. Gaps in globular cluster streams: Giant molecular clouds can cause them too. Mon. Not. R. Astron. Soc. 2016, 463, L17–L21. [Google Scholar] [CrossRef]
Erkal, D.; Koposov, S.E.; Belokurov, V. A sharper view of Pal 5’s tails: Discovery of stream perturbations with a novel non-parametric technique. Mon. Not. R. Astron. Soc. 2017, 470, 60–84. [Google Scholar] [CrossRef]
Bovy, J.; Erkal, D.; Sanders, J.L. Linear perturbation theory for tidal streams and the small-scale CDM power spectrum. Mon. Not. R. Astron. Soc. 2017, 466, 628–668. [Google Scholar] [CrossRef]
Mao, S.; Schneider, P. Evidence for substructure in lens galaxies? Mon. Not. R. Astron. Soc. 1998, 295, 587. [Google Scholar] [CrossRef]
Schneider, P.; Weiss, A. The gravitational lens equation near cusps. Astron. Astrophys. 1992, 260, 1–13. [Google Scholar]
Metcalf, R.B.; Madau, P. Compound Gravitational Lensing as a Probe of Dark Matter Substructure within Galaxy Halos. Astrophys. J. 2001, 563, 9–20. [Google Scholar] [CrossRef] [Green Version]
Dalal, N.; Kochanek, C.S. Direct Detection of Cold Dark Matter Substructure. Astrophys. J. 2002, 572, 25–33. [Google Scholar] [CrossRef]
Xu, D.; Sluse, D.; Gao, L.; Wang, J.; Frenk, C.; Mao, S.; Schneider, P.; Springel, V. How well can cold dark matter substructures account for the observed radio flux-ratio anomalies. Mon. Not. R. Astron. Soc. 2015, 447, 3189–3206. [Google Scholar] [CrossRef]
Hsueh, J.W.; Enzi, W.; Vegetti, S.; Auger, M.; Fassnacht, C.D.; Despali, G.; Koopmans, L.V.E.; McKean, J.P. SHARP – VII. New constraints on warm dark matter free-streaming properties and substructure abundance from flux-ratio anomalous lensed quasars. arXiv 2019, arXiv:1905.04182. [Google Scholar]
Iršič, V.; Viel, M.; Haehnelt, M.G.; Bolton, J.S.; Cristiani, S.; Becker, G.D.; D’Odorico, V.; Cupani, G.; Kim, T.S.; Berg, T.A.M.; et al. New constraints on the free-streaming of warm dark matter from intermediate and small scale Lyman-α forest data. Phys. Rev. D 2017, 96, 023522. [Google Scholar] [CrossRef]
Vegetti, S.; Koopmans, L.V.E. Bayesian strong gravitational-lens modelling on adaptive grids: Objective detection of mass substructure in Galaxies. Mon. Not. R. Astron. Soc. 2009, 392, 945–963. [Google Scholar] [CrossRef]
Vegetti, S.; Lagattuta, D.J.; McKean, J.P.; Auger, M.W.; Fassnacht, C.D.; Koopmans, L.V.E. Gravitational detection of a low-mass dark satellite galaxy at cosmological distance. Nature 2012, 481, 341–343. [Google Scholar] [CrossRef] [PubMed]
Keeton, C.R. A Catalog of Mass Models for Gravitational Lensing. arXiv 2001, arXiv:astro-ph/0102341. [Google Scholar]
Vegetti, S.; Koopmans, L.V.E.; Bolton, A.; Treu, T.; Gavazzi, R. Detection of a dark substructure through gravitational imaging. Mon. Not. R. Astron. Soc. 2010, 408, 1969–1981. [Google Scholar] [CrossRef]
Vegetti, S.; Despali, G.; Lovell, M.R.; Enzi, W. Constraining sterile neutrino cosmologies with strong gravitational lensing observations at redshift z = 0.2. Mon. Not. R. Astron. Soc. 2018, 481, 3661–3669. [Google Scholar] [CrossRef]
Vegetti, S.; Koopmans, L.V.E.; Auger, M.W.; Treu, T.; Bolton, A.S. Inference of the cold dark matter substructure mass function at z = 0.2 using strong gravitational lenses. Mon. Not. R. Astron. Soc. 2014, 442, 2017–2035. [Google Scholar] [CrossRef]
Hezaveh, Y.D.; Dalal, N.; Marrone, D.P.; Mao, Y.Y.; Morningstar, W.; Wen, D.; Blandford, R.D.; Carlstrom, J.E.; Fassnacht, C.D.; Holder, G.P.; et al. Detection of Lensing Substructure Using ALMA Observations of the Dusty Galaxy SDP.81. Astrophys. J. 2016, 823, 37. [Google Scholar] [CrossRef]
Nightingale, J.; Dye, S.; Massey, R. AutoLens: Automated Modeling of a Strong Lens’s Light, Mass and Source. arXiv 2017. [Google Scholar] [CrossRef]
Brewer, B.J.; Huijser, D.; Lewis, G.F. Trans-dimensional Bayesian inference for gravitational lens substructures. Mon. Not. R. Astron. Soc. 2016, 455, 1819–1829. [Google Scholar] [CrossRef]
Diaz Rivero, A.; Cyr-Racine, F.Y.; Dvorkin, C. On the Power Spectrum of Dark Matter Substructure in Strong Gravitational Lenses. arXiv 2017, arXiv:1707.04590. [Google Scholar] [CrossRef]
Li, R.; Frenk, C.S.; Cole, S.; Wang, Q.; Gao, L. Projection effects in the strong lensing study of subhaloes. Mon. Not. R. Astron. Soc. 2017, 468, 1426–1432. [Google Scholar] [CrossRef] [Green Version]
Despali, G.; Vegetti, S.; White, S.D.M.; Giocoli, C.; van den Bosch, F.C. Modelling the line-of-sight contribution in substructure lensing. Mon. Not. R. Astron. Soc. 2018, 475, 5424–5442. [Google Scholar] [CrossRef]
Li, R.; Frenk, C.S.; Cole, S.; Gao, L.; Bose, S.; Hellwing, W.A. Constraints on the identity of the dark matter from strong gravitational lenses. Mon. Not. R. Astron. Soc. 2016, 460, 363–372. [Google Scholar] [CrossRef] [Green Version]

1	Equal amounts of dark matter and anti-dark matter.
2	By strong, we mean that the cross-section for self-interaction is of the order of the nuclear cross-section for visible matter (set by the strong force).
3	Some SIDM models are motivated by the baryon asymmetry; in these models, dark matter, unlike traditional WIMPs, shares this asymmetry (for a review of asymmetric dark matter see [19]).
4	In contrast to WDM, the damping of small structures is not due to free-streaming, but to a collisional, Silk-like, damping.
5	The (co-moving) free-streaming scale is given by: $l_{fs} = 2 c t_{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 $3 k_{B} T_{nr} \sim m_{χ} c^{2}$ ); $a_{nr} = 1 / (1 + z_{nr})$ is the scale factor at $t_{nr}$ ( $a \propto 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 Equation (43) of [23].
7	Please 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 \bar{ρ} {(π / k)}^{3}$ , where $\bar{ρ}$ is the mean dark matter density today.
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, Equation (6) is simply a multiplication $\hat{\tilde{Φ}} (\vec{k}) = \hat{g} (\vec{k}) \hat{\tilde{ρ}} (\vec{k})$ .
14	For a review of the force computation methods see Section 3.5 of [46].
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.
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	Reproduced from Michael Boylan-Kolchin et al. Resolving cosmic structure formation with the Millennium-II Simulation. MNRAS (2009) 398 (3): 1150-1164, doi: 10.1111/j.1365-2966.2009.15191.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following u. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
19	For a review of the state of cosmological simulations circa 2012 see [79].
20	This is only true if, on the scales of interest, the primordial power spectrum grows monotonically towards large k.
21	Reproduced from Volker Springel et al. The Aquarius Project: the subhaloes of galactic haloes. MNRAS (2008) 391 (4): 1685–1711, doi: 10.1111/j.1365-2966.2008.14066.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/391/4/1685/1747035.
22	Reproduced from Julio Navarro et al. The diversity and similarity of simulated cold dark matter haloes. MNRAS (2010) 402 (1): 21–34, doi: 10.1111/j.1365-2966.2009.15878.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/402/1/21/1028856.
23	Reproduced from Carlos Vera-Ciro et al. The shape of dark matter haloes in the Aquarius simulations: evolution and memory. MNRAS (2011) 416 (2): 1377–1391, doi: 10.1111/j.1365-2966.2011.19134.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/416/2/1377/1061105.
24	Reproduced from Mark Vogelsberger et al. Phase-space structure in the local dark matter distribution and its signature in direct detection experiments. MNRAS (2009) 395 (2): 797–811, doi: 10.1111/j.1365-2966.2009.14630.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/395/2/797/1747020. The figures mentioned in footnotes 21–24 are not included under the CC-BY license of this publication. For permissions, please email: [email protected].
25	Given the limited resolution of simulations, local in this sense refers to regions of at least $O$ (10 kpc $^{3}$ ) as in [140].
26	Reproduced from Mark Vogelsberger et al. Subhaloes in self-interacting galactic dark matter haloes. MNRAS (2012) 423 (4): 3740–3752, doi: 10.1111/j.1365-2966.2012.21182.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/423/4/3740/1749150. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
27	Reproduced from Mark Lovell et al. The haloes of bright satellite galaxies in a warm dark matter universe. MNRAS (2012) 420 (3): 2318–2324, doi: 10.1111/j.1365-2966.2011.20200.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/420/3/2318/979379. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
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 behavior.
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 center. 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].
30	Reproduced from Carlo Giocoli et al. The substructure hierarchy in dark matter haloes . MNRAS (2010) 404 (1): 502–517, doi: 10.1111/j.1365-2966.2010.16311.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/404/1/502/3101607. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
31	Reproduced from Volker Springel et al. The Aquarius Project: the subhaloes of galactic haloes. MNRAS (2008) 391 (4): 1685–1711, doi: 10.1111/j.1365-2966.2008.14066.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/391/4/1685/1747035. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
32	Another set of parameters that can be used to define the orbit are the apocentre and pericenter. Different parametrizations can be transformed into one another since they are all related to the potential, $ϕ (r)$ .
33	Reproduced from Lilian Jian et al. Orbital parameters of infalling satellite haloes in the hierarchical $Λ$ CDM model. MNRAS (2015) 448 (2): 1674–1686, doi: 10.1093/mnras/stv053. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/448/2/1674/1053529. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
34	This estimate is based on an extrapolation over many orders of magnitude of the subhalo mass function determined in simulations down to the free-streaming mass of WIMP-CDM particles. We discuss this in more detail below.
35	Reproduced from Miguel Rocha et al. Infall times for Milky Way satellites from their present-day kinematics. MNRAS (2012) 425 (1): 231–244, doi: 10.1111/j.1365-2966.2012.21432.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/425/1/231/998181. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
36	Reproduced from Mark Vogelsberger and Simon D. M. White. Streams and caustics: the fine-grained structure of $Λ$ cold dark matter haloes. MNRAS (2011) 413 (2): 1419–1438, doi: 10.1111/j.1365-2966.2011.18224.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/413/2/1419/1070092. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
37	Caustics represent folds in the fine-grained phase-space distribution function, which in CDM evolves according to the collisionless Boltzmann equation (Equation (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 Equation (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 Section 2.4), which implies that the original thin sheet can be stretched and folded but it cannot be broken. Caustics appear where folds occur, and have very large spatial densities, limited only by primordial thermal motions (e.g., [190,191,192]).
38	This equation ignores the effects of the centrifugal force on the satellite as it orbits around the host. Including this effect (assuming circular orbits) modifies Equation (13) slightly by changing the factor of 2 in the denominator to 3 [203].
39	It is possible to evaluate Equation (15) for a given fixed potential and a given subhalo orbit (see e.g., Equation (20) of [209] for an NFW halo).
40	Reproduced from Frank C van den Bosch and Go Ogiya. Dark matter substructure in numerical simulations: a tale of discreteness noise, runaway instabilities, and artificial disruption. MNRAS (2018) 475 (3): 4066–4087, doi: 10.1093/mnras/sty084. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/475/3/4066/4797185. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
41	Reproduced from Frank C van den Bosch et al. Disruption of dark matter substructure: fact or fiction? MNRAS (2018) 474 (3): 3043–3066, doi: 10.1093/mnras/stx2956. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/474/3/3043/4638541. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
42	©AAS. Reproduced with permission. For the original article, please visit the following https://iopscience.iop.org/article/10.1086/420840.
43	Reproduced from Go Ogiya and Andreas Burkert. Dynamical friction and scratches of orbiting satellite galaxies on host systems. MNRAS (2016) 457 (2): 2164–2172, doi: 10.1093/mnras/stw091. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/457/2/2164/970692. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
44	Reproduced from Michael Boylan-Kolchin et al. Dynamical friction and galaxy merging timescales . MNRAS (2008) 383 (1): 93–101, doi: 10.1111/j.1365-2966.2007.12530.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/383/1/93/1067887. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
45	The Hernquist halo profile [221] has the same asymptotic behavior at the center as the NFW halo and has the advantage that the velocity distribution function in the isotropic case has an analytic form (see Equation (11)), which makes it particularly simple to set up initial conditions for simulating haloes in dynamical equilibrium.
46	The values for these parameters reported in [218] are: $A = 0.216$ , $b = 1.3$ , $c = 1.9$ , and $d = 1.0$ , but we point out that in this study both the halo and the subhalo were modeled as Hernquist profiles.
47	Equation (18) was only explored for values of the circularity in the range $0.3 \leq η \leq 1.0$ and for $0.65 \leq r_{c} / r_{200} \leq 1.0$ ; the lower limits were imposed to avoid radial orbits that would take the subhalo so close to the center 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 Equation (18) was not investigated outside this range and might not be valid there even in the absence of a central galaxy.
48	©AAS. Reproduced with permission. For the original article, please visit the following https://iopscience.iop.org/article/10.1088/0004-637X/740/2/102.
49	The fitting parameters in Equation (19) in the case of the Millennium simulations may be found in [230], where a redshift dependence is also provided.
50	Although the introduction of a third parameter will obviously improve the quality of the fit, the Einasto profile is, in fact, a slightly better fit to simulations than the 2-parameter NFW profile even after one of the parameters ( $α_{E}$ ) is fixed to an appropriate value. For instance, fixing $α_{E} \sim 0.16$ gives a better fit than NFW to haloes across a range of halo masses [240].
51	Reproduced from Volker Springel et al. The Aquarius Project: the subhaloes of galactic haloes. MNRAS (2008) 391 (4): 1685–1711, doi: 10.1111/j.1365-2966.2008.14066.x. By permission of Oxford University Press on behalf of the Royal Astronomical Society. For the original article, please visit the following https://academic.oup.com/mnras/article/391/4/1685/1747035. This figure is not included under the CC-BY license of this publication. For permissions, please email: [email protected].
52	We note that there is a typo in the caption of Figure 28 in [138], which gives the fitting function for $δ_{V}$ and $m_{sub}$ ( $5.8 \times 10^{8}$ M $_{⊙} \to 5.8 \times 10^{4}$ M $_{⊙}$ ).
53	Obtained by taking the power-law concentration-mass relation in [129] (their Equation (4)), and using Equations (23) and (24).
54	For a clear illustration of the evolutionary track of subhaloes in the $r_{\max} - V_{\max}$ plane due to tidal stripping, see Figure 8 of [243].
55	By galactic-scale haloes and subhaloes, we mean self-bound dark matter structures that can potentially host a galaxy, i.e., is haloes of mass above ∼ $10^{8}$ M $_{⊙}$ , in which gas can cool by atomic processes (e.g., [68,250]).
56	This is true only for elastic SIDM, and for cross sections that do not exceed the gravothermal collapse threshold, $σ_{T} / m_{χ} \sim 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,249,251,252] 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., [253]), and in which scattering between the excited and ground states can result in energy injection at the center of dark matter haloes thus altering their structure. Only until very recently have these models began to be explored with simulations [166,254].
58	Large relative velocities between gas and dark matter inherited from the photon-baryon coupling before recombination can impede the growth of gravitational perturbations and stop gas from accreting into the first haloes [257]. This process, however, is only thought to be relevant for the formation of the first stars.
59	This mass threshold is smaller at higher redshifts, see e.g., Figure 3 of [266].
60	Contrary to some claims, xmm data for Draco, and Hitomi data for Perseus, are consistent with a 7 kev neutrino [310,311].
61	Other baryon effects that can transform cusps into cores have been proposed (e.g., [321,322]) but have been less studied.
62	See [325] for an opposed view.
63	This halo mass was estimated assuming a truncated pseudo-Jaffe profile (see e.g., Equation (42) in [350]). The inferred mass is likely to be larger if an NFW profile is assumed instead. For instance, a similar dark matter substructure detected with lensing was reported by [351] with a mass of ∼ $3.5 \times 10^{9}$ M $_{⊙}$ assuming a truncated pseudo-Jaffe profile, while assuming an NFW profile this substructure is estimated to have a mass of ∼ $10^{10}$ M $_{⊙}$ [352].
64	This has been demonstrated explicitly for the case of Einstein ring distortions but it may hold true for other tests as well.

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 slightly 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].

α

forest constraints on WDM [28] to behave as CDM. Figure adapted from [32].

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 realization 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.

2 π / L

, and the Nyquist mode,

π / d

, where L and d are the cube length and interparticle separation, respectively. Right: a realization 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.

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]). The dashed lines show the Poisson noise limit for the Millennium I (left) and Millennium-II (right) simulations. Figure adapted from [77]18.

Figure 3. Emergence of the cosmic web. Left: evolution of the (projected) dark matter density field in a slab of length

L = 100

Figure 4. The galaxy distribution in various redshift surveys and in mock catalogues constructed from the Millennium simulation [4]. The small slice at the top shows the CfA2 “Great Wall” [81], with the Coma cluster at the center. Just above is a section of the Sloan Digital 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].

Figure 5. Halo mass function for different dark matter models (adapted from [20]). Left: The large-scale dark matter distribution in a slab of a 64 Mpc/h cube in different cosmologies: CDM and WDM in the top left and bottom right, respectively; two interacting dark matter models in the other two panels. Right: The halo mass function at

z = 0

for the models on the left. The transparent light blue region marks the resolution limit of the simulations. The cutoff in the primordial linear power spectrum of the non-CDM models results in a lower abundance of low-mass haloes, visible in the panels on the left and quantified in the halo mass function on the right.

z = 0

Figure 6. The structure of CDM haloes. The different panels show several characteristics of the spatial (left) and dynamical (right) structure of a Milky Way-size CDM halo (

M_{200} \sim 1.8 \times 10^{12}

_{⊙}

;

r_{200} \sim 250

kpc) from the Aquarius project [138]. The top panels show the spherically averaged radial density (left; [138]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})

; [139]23) and local dark matter velocity distribution near the solar circle:

2 kpc < r < 9 kpc

(right; [140]24).

Figure 6. The structure of CDM haloes. The different panels show several characteristics of the spatial (left) and dynamical (right) structure of a Milky Way-size CDM halo (

M_{200} \sim 1.8 \times 10^{12}

_{⊙}

;

r_{200} \sim 250

T = (a^{2} - b^{2}) / (a^{2} - c^{2})

; [139]23) and local dark matter velocity distribution near the solar circle:

2 kpc < r < 9 kpc

(right; [140]24).

Figure 7. Structure of haloes in models with different types of dark matter: collisional (SIDM;

σ_{T} / m_{χ} ≳ 1

^{2}

/g) and with a galactic-scale free-streaming cutoff (WDM;

m_{χ} \sim 2.3

keV). Upper panels: projected dark matter distribution of a Milky Way-size halo in the SIDM model (left panel; [57]26) and in the WDM model (right panel; [153]27). 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} \propto r^{- 1.875}

, at most radii, except in the very center, 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.

Figure 7. Structure of haloes in models with different types of dark matter: collisional (SIDM;

σ_{T} / m_{χ} ≳ 1

^{2}

/g) and with a galactic-scale free-streaming cutoff (WDM;

m_{χ} \sim 2.3

ρ / σ^{3} \propto r^{- 1.875}

Figure 8. Dark matter subhaloes. Left: schematic representation of a dark matter halo merger tree (taken from [176]30) at discrete redshifts. In a hierarchical model, haloes grow by the accretion of smaller neighboring haloes (A,B,C,D), which become subhaloes at the time when they first cross the virial radius of the host halo. The main branch of the tree represents the evolution of the main progenitor (shown in blue). Since this process occurs across the entire hierarchy of structures, there are subhaloes within subhaloes (sub-subhaloes; like a, b, c in system D) and so on. Right: a simulated Milky Way-size CDM halo from the Aquarius project (figure taken from [138]31; this is the same halo illustrated in Figure 6). The circles in the main image mark six subhaloes that are shown enlarged in the surrounding panels, as indicated by the labels. Sub-subhaloes are clearly visible (corresponding to the configuration illustrated in the last step,

z_{0}

, in the left panel). The bottom row shows several generations of sub-subhaloes contained within subhalo f.

z_{0}

, in the left panel). The bottom row shows several generations of sub-subhaloes contained within subhalo f.

Figure 9. Initial conditions for the orbits of subhaloes infalling into haloes of present-day mass

M_{200} = 10^{12}

_{⊙}

(figures taken from [186]33; see that paper for similar plots for other host masses). Upper left: the distribution of formation redshifts (defined as the redshift at which the mass of the main progenitor of the halo was half its present value). These and the other histograms in this figure are normalized such that the integral over the distribution is unity. Lower left: distribution of infall (accretion) redshifts of subhaloes of different mass ratios,

ξ

(relative to the host halo at the time of accretion; see legend). Middle: distributions of radial (upper panel) and tangential (lower panel) subhalo orbital velocities at infall, relative to the virial velocity of the host, for the same halo mass and subhalo-to-halo mass ratios as in the lower-left panel. Right: bivariate distribution of orbital parameters for infalling haloes into hosts of mass

M_{200} (z = 0) = 10^{13}

_{⊙}

Figure 9. Initial conditions for the orbits of subhaloes infalling into haloes of present-day mass

M_{200} = 10^{12}

_{⊙}

ξ

M_{200} (z = 0) = 10^{13}

_{⊙}

Figure 10. Left: Distribution in the 2D radial phase-space plane of subhaloes identified in a Milky Way-size halo simulation (Via Lactea II [188]; figure adapted from [189]35). Subhaloes are color-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 center. Right: the 2D radial phase-space structure of simulation particles in a different Milky Way-size halo simulation (Aquarius [138]; 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.

z = 0

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 [208]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 [209]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 (Equation (14)). In the lower panel, tidal shocking is seen to reduce the mass in the central regions but preserves the asymptotic NFW behavior, while the outer regions become considerably steeper than NFW (figure adapted from [215]42).

Figure 12. Dynamical friction experienced by subhaloes. Left: simulation of a subhalo orbiting a Milky Way-size halo; the initial mass ratio and circularity of the orbit are 0.1 and 0.5, respectively. The images show the projected over- (or under-) density relative to the initial value at

t = 0

, at different times during the evolution. The thick solid line marks the subhalo orbit, which decays over time due to dynamical friction. This gravitational process induces a wake in the host halo trailing behind the satellite (most clearly visible in the top left panel). The dipole feature at the center of the host is caused by the tidal effect of the subhalo, which perturbs the position of the halo potential minimum. This effect diminishes with time as the satellite is stripped of mass (figure adapted from [217]43). Right: evolution of the radial distance of a simulated subhalo orbiting a Milky Way-size halo (figure taken from [218]44). The orbit decays by dynamical friction on a timescale that strongly depends on the initial mass ratio (different colors) and circularity of the orbit (dashed and solid lines). The timescales are well approximated by the fitting formula (Equation (18)), which is an improvement over the classical Chandrasekhar formula (Equation (17)).

t = 0

Figure 13. Subhalo abundance. Left: the subhalo velocity function at

z = 0

for haloes of different maximum circular velocity, from ∼150 km/s to ∼1000 km/s (bottom to top). In terms of the velocity ratio,

x = V_{sub} / V_{h}

, the velocity function is nearly universal, scaling as

x^{- 3}

(dashed line) with a scale-dependent normalization (see Equation (20); figure adapted from [223]48). Right: the number density of subhaloes as a function of halocentric distance in units of the virial radius for Milky Way-size haloes (triangles) and cluster-size haloes (circles). All subhaloes with

m_{sub} / M_{h} > 10^{- 5}

have been included. The dashed lines are the average NFW fits to the density profiles of the hosts. These functions have been normalized to unity at the virial radius. Figure adapted from [224].

Figure 13. Subhalo abundance. Left: the subhalo velocity function at

z = 0

for haloes of different maximum circular velocity, from ∼150 km/s to ∼1000 km/s (bottom to top). In terms of the velocity ratio,

x = V_{sub} / V_{h}

, the velocity function is nearly universal, scaling as

x^{- 3}

m_{sub} / M_{h} > 10^{- 5}

have been included. The dashed lines are the average NFW fits to the density profiles of the hosts. These functions have been normalized to unity at the virial radius. Figure adapted from [224].

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 (color 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 (color 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 Equation (24)) within different radial ranges, as given in the legend. The solid line corresponds to isolated haloes. Figures adapted from [138]51.

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}

68 %

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Zavala, J.; Frenk, C.S. Dark Matter Haloes and Subhaloes. Galaxies 2019, 7, 81. https://doi.org/10.3390/galaxies7040081

AMA Style

Zavala J, Frenk CS. Dark Matter Haloes and Subhaloes. Galaxies. 2019; 7(4):81. https://doi.org/10.3390/galaxies7040081

Chicago/Turabian Style

Zavala, Jesús, and Carlos S. Frenk. 2019. "Dark Matter Haloes and Subhaloes" Galaxies 7, no. 4: 81. https://doi.org/10.3390/galaxies7040081

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Dark Matter Haloes and Subhaloes

Abstract

1. Introduction

2. Formation of Dark Matter Haloes

2.1. Initial Conditions: The Primordial Power Spectrum in the Linear Regime

2.2. The Non-Linear Regime: N-Body Simulation Methods

2.3. The Non-Linear Regime: Initial Conditions and The Emergence of the Cosmic Web

2.4. The Structural Properties of Dark Matter Haloes

3. Halo Mergers and the Emergence of Subhaloes

3.1. Halo Mass Assembly: Smooth Accretion vs Mergers

3.2. Evolution of Subhaloes: Initial Conditions

3.3. Dynamics of Subhaloes

3.4. The Abundance, Spatial Distribution and Internal Structure of Dark Matter Subhaloes

3.5. The Impact of the Nature of the Dark Matter

4. Outlook

4.1. The Impact of Baryonic Physics on Dark Matter Structure

4.2. Astrophysical Tests of the Nature of the Dark Matter

Funding

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI