Predicting the density profiles of the first halos

TL;DR

This paper develops and validates universal models to predict the density profiles of the first dark matter halos from initial density peaks, enabling their use as probes of early Universe physics and dark matter properties.

Contribution

The authors introduce models that accurately predict the density profiles of the first halos from peak properties, validated across diverse cosmologies and applicable to various scenarios.

Findings

Models predict the $ ho ightarrow r^{-3/2}$ asymptote and halo properties like $v_{max}$ and $r_{max}.

Models validated with six diverse $N$-body simulations.

Effective in narrow power spectrum scenarios, with potential for broader spectra understanding.

Abstract

The first dark matter halos form by direct collapse from peaks in the matter density field, and evidence from numerical simulations and other analyses suggests that the dense inner regions of these objects largely persist today. These halos would be the densest dark matter structures in the Universe, and their abundance can probe processes that leave imprints on the primordial density field, such as inflation or an early matter-dominated era. They can also probe dark matter through its free-streaming scale. The first halos are qualitatively different from halos that form by hierarchical clustering, as evidenced by their $\rho\propto r^{-3/2}$ inner density profiles. In this work, we present and tune models that predict the density profiles of these halos from properties of the density peaks from which they collapsed. These models predict the coefficient $A$ of the $\rho=Ar^{-3/2}$ small-radius asymptote of the density profile along with the maximum circular velocity $v_\mathrm{max}$ and associated radius $r_\mathrm{max}$. These models are universal; they can be applied to any cosmology, and we confirm this by validating them using six $N$-body simulations carried out in wildly disparate cosmological scenarios. We find that these models can even predict the full population of halos with reasonable accuracy in scenarios with narrowly supported power spectra, although for broader power spectra, an understanding of the impact of halo mergers is needed. With their connection to the primordial density field established, the first dark matter halos will serve as probes of the early Universe and the nature of dark matter.

Tables

Table 1. Table 1: The simulation list with basic parameters. Also listed are σ 0 subscript 𝜎 0 \sigma_{0} , the rms density variance, and R ∗ subscript 𝑅 R_{*} , the comoving characteristic correlation length, both calculated in linear theory. When matter dominates, σ ​ ( a ) ∝ a proportional-to 𝜎 𝑎 𝑎 \sigma(a)\propto a , which motivates the scaling in this column (with a = 1 𝑎 1 a=1 today).

Simulation	$z_{\mathrm{init}}$	$z_{\mathrm{final}}$	Box (kpc)	${\sigma }_0\left(a\right)/a$	$R_{\ast }$ (kpc)
$w=0.1$	${10}^4$	$50$	$7.4$	$164$	$0.25$
$w=0.3$	${10}^4$	$50$	$7.4$	$164$	$0.19$
$w=0.5$	${10}^4$	$50$	$7.4$	$164$	$0.12$
EMDE	$3\times {10}^4$	$100$	$1.5\times {10}^{-3}$	$330$	$4.1\times {10}^{-5}$
$100$ GeV	$500$	$9$	$0.15$	$30$	$1.9\times {10}^{-3}$
$3.5$ keV	$150$	$2$	$4.4\times {10}^3$	$7.6$	$66$

Table 2. Table 2: For each simulation, we list the halo count N halo subscript 𝑁 halo N_{\mathrm{halo}} above the mass cutoff at the final redshift, the number of peaks N peak subscript 𝑁 peak N_{\mathrm{peak}} in the linear density field, and the number N match subscript 𝑁 match N_{\mathrm{match}} of successful halo-peak matches. Of these matches, also listed are the number N A subscript 𝑁 A N_{\mathrm{A}} with density profiles resolved at small radii, the number N r subscript 𝑁 r N_{\mathrm{r}} with r max < r vir subscript 𝑟 max subscript 𝑟 vir r_{\mathrm{max}}<r_{\mathrm{vir}} , and the subsets N A < 3 ​ M ​ M superscript subscript 𝑁 A absent 3 M M N_{\mathrm{A}}^{<3\mathrm{MM}} and N r < 3 ​ M ​ M superscript subscript 𝑁 r absent 3 M M N_{\mathrm{r}}^{<3\mathrm{MM}} that also underwent fewer than three major mergers.

Simulation	$N_{\mathrm{halo}}$	$N_{\mathrm{peak}}$	$N_{\mathrm{match}}$	$N_A$		$N_r$
$w=0.1$	696	439	411	408	350	404	347
$w=0.3$	801	933	565	561	481	537	457
$w=0.5$	1011	3896	795	790	683	738	634
EMDE	610	730	424	416	366	388	339
100 GeV	216	6108	157	153	141	82	73
3.5 keV	581	4401	463	455	430	322	303

Table 3. Table 3: A summary of the predictive models along with their simulation-tuned coefficients. In each model, the halo quantity is proportional to the peak quantity with the given coefficient between them. The mean and rms deviation σ 𝜎 \sigma in the proportionality coefficient are computed in log space. We express the spread as a ratio: in particular, the middle 50% spread is the ratio between the 75th and 25th percentiles, while the 1 σ 𝜎 \sigma spread is the ratio corresponding to the log-space σ 𝜎 \sigma . We also list the number of halos contributing to the coefficient statistics.

Halo quantity	Peak quantity	Model	Median coefficient	Middle 50% spread	Mean coefficient	1$\sigma$ spread	Sample size	Section ref.
$A$	${\overline{\rho }}_0{\delta }^{9/4}{\left|{\nabla }^2\delta \right|}^{-3/4}$	Spherical collapse	$8.76$	$1.43$	$8.50$	$1.44$	$2451$	III.1
$A$	${\overline{\rho }}_0{\delta }^{9/4}{\left|{\nabla }^2\delta \right|}^{-3/4}{f}_{\mathrm{ec}}^{-3/2}(e,p)$	Ellipsoidal collapse	$12.1$	$1.31$	$12.2$	$1.36$	$2450$171717One peak is so ellipsoidal that Eq. (6) has no solution, so this peak is discarded from analysis.	III.2
$r_{\max }$	$q_{\max }/\Delta \left(q_{\max }\right)$	Turnaround,	$0.131$	$1.34$	$0.132$	$1.39$	$1759$	IV.1
$r_{\max }$	$q_v/\Delta \left(q_v\right)$	Turnaround, $q_{\max }>q_v$	$0.088$	$2.18$	$0.085$	$1.75$	$712$	IV.1
$r_{\max }$	$q_{\max }/\left[\Delta \left(q_{\max }\right)X\left(q_{\max }\right)\right]$	$s=0$ contraction	$0.414$	$1.40$181818The $s=0$ coefficient has spread $1.33$ and $1.74$ over halos for which the turnaround model predicts and $q_{\max }>q_v$, respectively, implying that the contraction model is superior or comparable to the turnaround model in both cases.	$0.406$	$1.52$	$2471$	IV.2
$r_{\max }$	$q_{\max }/\left[\Delta \left(q_{\max }\right)X\left(q_{\max }\right)\right]$	$s=1$ contraction	$0.846$	$1.39$191919Similarly, the $s=1$ coefficient has spread $1.33$ and $1.68$ in these two cases.	$0.848$	$1.49$	$2471$	IV.2
$r_{\max }$	$q_{\max }/\kappa \left(q_{\max }\right)$	Virialization,	$0.042$	$1.36$	$0.042$	$1.37$	$1578$	IV.3
$r_{\max }$	$q_v/\Delta \left(q_v\right)$	Virialization, $q_{\max }>q_{\mathrm{ta}}$	$0.092$	$1.95$	$0.087$	$1.70$	$893$	IV.3
$M\left(r_{\max }\right)$	$\left(4\pi /3\right)q_{\max }^3{\overline{\rho }}_0$	Turnaround,	$0.273$	$1.64$	$0.258$	$1.84$	$1514$	IV.1
$M\left(r_{\max }\right)$	$\left(4\pi /3\right)q_v^3{\overline{\rho }}_0$	Turnaround, $q_{\max }>q_v$	$0.134$	$2.97$	$0.123$	$3.13$	$639$	IV.1
$M\left(r_{\max }\right)$	$\left(4\pi /3\right)q_{\max }^3{\overline{\rho }}_{0}X\left(q_{\max }\right)$	$s=0$ contraction	$0.441$	$1.76$	$0.396$	$2.30$	$2153$	IV.2
$M\left(r_{\max }\right)$	$\left(4\pi /3\right)q_{\max }^3{\overline{\rho }}_{0}X\left(q_{\max }\right)$	$s=1$ contraction	$0.658$	$1.73$	$0.619$	$2.23$	$2153$	IV.2
$M\left(r_{\max }\right)$	$\left(4\pi /3\right)q_{\max }^3{\overline{\rho }}_0$	Virialization,	$0.150$	$1.59$	$0.145$	$1.71$	$1341$	IV.3
$M\left(r_{\max }\right)$	$\left(4\pi /3\right)q_v^3{\overline{\rho }}_0$	Virialization, $q_{\max }>q_{\mathrm{ta}}$	$0.143$	$2.79$	$0.129$	$2.98$	$812$	IV.3

Table 4. Table 4: This table shows how accurately our models predict aggregate halo signals. We consider simplified aggregate signals ∑ A 2 superscript 𝐴 2 \sum A^{2} for annihilation and ∑ M ​ ( r max ) ​ r max 𝑀 subscript 𝑟 max subscript 𝑟 max \sum M(r_{\mathrm{max}})r_{\mathrm{max}} for lensing, each summed over the population of predicted or simulated halos. In each case, we list the ratio of the signal predicted from the linear density field to the signal aggregated over halos in the simulation box. We use the ellipsoidal collapse and s = 0 𝑠 0 s=0 contraction models.

Simulation	$\sum A^2$	$\sum M\left(r_{\max }\right)r_{\max }$
$w=0.1$	0.91	0.56
$w=0.3$	0.89	0.69
$w=0.5$	0.91	2.2
EMDE	0.78	0.66
100 GeV	1.15	33
3.5 keV	1.16	14

Equations

$$\mathcal{P}(k)=\frac{\mathcal{A}}{\sqrt{2\pi }w}\exp \left[-\frac{1}{2}{\left(\frac{\ln (k/k_s)}{w}\right)}^2\right]$$

$${\sigma }_j$$

$$R_{\ast }$$

$$\tilde{\delta }(\mathbf{x})\equiv \delta (\mathbf{x},a)/a,$$

$$A=\alpha {\bar{\rho }}_0{\delta }^{9/4}|{\nabla }^2\delta {|}^{-3/4},$$

$$f_{\mathrm{ec}}(e,p)=1+0.47{\left[5(e^2-p|p|)f_{\mathrm{ec}}^2(e,p)\right]}^{0.615}.$$

$$A=\alpha {\bar{\rho }}_0{\delta }^{9/4}|{\nabla }^2\delta {|}^{-3/4}{f}_{\mathrm{ec}}^{-3/2}(e,p),$$

$$\Delta (q)=\frac{3}{q^3}{\int }_0^q\delta (q^{\prime })q^{\prime 2}\mathrm{d}{q}^{\prime },$$

$$r=\beta q/\Delta (q),$$

$$M(q)={\beta }_M\frac{4\pi }{3}{q}^3{\bar{\rho }}_0$$

$$\frac{\mathrm{d}\ln M}{\mathrm{d}\ln r}=\frac{3}{1+3ϵ(q)},$$

$$ϵ(q)\equiv -\frac{1}{3}\frac{\mathrm{d}\ln \Delta }{\mathrm{d}\ln q}=1-\frac{\delta (q)}{\Delta (q)}.$$

$$r_{\max }={\beta }^{\prime }{q}_{\mathrm{v}}/\Delta (q_{\mathrm{v}}),$$

$$X(q)\equiv 1+\frac{3}{q^3}{\int }_q^{q_{\mathrm{v}}}{q}^{\prime 2}\mathrm{d}{q}^{\prime }f\left[r(q),r(q^{\prime })\right].$$

$$r(q)=\beta \frac{q}{\Delta (q)}\frac{1}{X(q)}$$

$$\frac{\mathrm{d}\ln X}{\mathrm{d}\ln q}=-\frac{3-\left[3(3-s)ϵ(q)-s\right](X-1)}{1+(4-s)(X-1)}$$

$$M(q)={\beta }_M(4\pi /3)q^3{\bar{\rho }}_{0}X(q)$$

$$\frac{\mathrm{d}\ln M}{\mathrm{d}\ln r}=\frac{3+\mathrm{d}\ln X/\mathrm{d}\ln q}{1+3ϵ(q)-\mathrm{d}\ln X/\mathrm{d}\ln q}$$

$$r(q)=\beta \frac{q}{\kappa (q)},\kappa (q)\equiv \frac{2}{q^5}{\int }_0^q{q}^{\prime 4}\mathrm{d}{q}^{\prime }\Delta (q^{\prime })$$

$$\delta (\mathbf{k})\equiv \int {\mathrm{d}}^3\mathbf{x}{\mathrm{e}}^{-i\mathbf{k}\cdot \mathbf{x}}\delta (\mathbf{x}).$$

$${\partial }_i{\partial }_j\delta =-\int \frac{{\mathrm{d}}^3\mathbf{k}}{(2\pi {)}^3}{\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}}{k}_i{k}_j\delta (\mathbf{k}).$$

$${\partial }_i{\partial }_j\varphi =4\pi G{\bar{\rho }}_0\int \frac{{\mathrm{d}}^3\mathbf{k}}{(2\pi {)}^3}{\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}}\frac{k_i{k}_j}{{\mathbf{k}}^2}\delta (\mathbf{k}).$$

$$e_{\varphi }\equiv \frac{{\lambda }_1-{\lambda }_3}{2({\lambda }_1+{\lambda }_2+{\lambda }_3)}\mathrm{and}{p}_{\varphi }\equiv \frac{{\lambda }_1+{\lambda }_3-2{\lambda }_2}{2({\lambda }_1+{\lambda }_2+{\lambda }_3)}.$$

$$\left\{\begin{array}[]{c}\delta(r)\\ \Delta(r)\\ \zeta(r)\\ \end{array}\right\}=\int\!\frac{\mathrm{d}^{3}\bm{k}}{(2\pi)^{3}}\mathrm{e}^{i\bm{k}\cdot\bm{x}}\delta(\bm{k})\left\{\begin{array}[]{c}\operatorname{sinc}(kr)\\ W(kr)\\ \cos(kr)\\ \end{array}\right\},$$

$$\frac{\mathrm{d}\delta }{\mathrm{d}\ln r}=\zeta (r)-\delta (r)\text{and}\frac{\mathrm{d}\Delta }{\mathrm{d}\ln r}=3\left[\delta (r)-\Delta (r)\right]$$

$$\frac{{\mathrm{d}}^{2}n}{\mathrm{d}\nu \mathrm{d}x}=\frac{{\mathrm{e}}^{-{\nu }^2/2}}{(2\pi {)}^2{R}_{\ast }^3}f(x)\frac{\exp \left[-\frac{1}{2}(x-\gamma \nu {)}^2/(1-{\gamma }^2)\right]}{[2\pi (1-{\gamma }^2){]}^{1/2}},$$

$$f(x)\equiv$$

$$+$$

$$\frac{\mathrm{d}n}{\mathrm{d}x}=\frac{f(x)}{8{\pi }^2{R}_{\ast }^3}{\mathrm{e}}^{-x^2/2}\left[1+\mathrm{erf}\left(\frac{x\gamma }{\sqrt{2(1-{\gamma }^2)}}\right)\right],$$

$$f(e,p|\nu )=\frac{1125}{\sqrt{10\pi }}e(e^2-p^2){\nu }^5\exp \left[-\frac{5}{2}{\nu }^2(3e^2+p^2)\right].$$

$$\begin{array}{rl}F(e|\nu )=& {\mathrm{e}}^{-\frac{15}{2}{e}^2{\nu }^2}\left(1-15e^2{\nu }^2\right)\mathrm{erf}\left(\sqrt{5/2}e\nu \right)\\ & -3\sqrt{10/\pi }e\nu {\mathrm{e}}^{-10e^2{\nu }^2}+\mathrm{erf}\left(\sqrt{10}e\nu \right)\end{array}$$

Full text

I Introduction

Decades of work have been devoted to understanding the halos that form by the gravitational collapse of collisionless dark matter. In the cold dark matter (CDM) model, $N$-body simulations have demonstrated that the radial density profiles $\rho(r)$ of these halos are well described by a remarkably universal form, known as the Navarro-Frenk-White (NFW) profile Navarro et al. (1996, 1997), which is a double power law that transitions from $\rho\propto r^{-1}$ at small radii to $\rho\propto r^{-3}$ at large radii. This profile has received only minor corrections since its introduction (e.g., Ref. Navarro et al. (2010)). In the CDM picture, all halos form by hierarchical clustering of smaller halos, and the NFW profile appears to be the generic consequence.

CDM represents an idealized scenario, however. In reality, dark matter particles are expected to have a nonzero temperature, and the corresponding random particle motions wash out density fluctuations smaller than a characteristic free-streaming scale. The first halos form by direct collapse of overdense regions at this scale, and $N$-body simulations with sufficient resolution show that these halos possess a markedly different density profile that asymptotes to $\rho\propto r^{-3/2}$ at small radii Ishiyama et al. (2010); Anderhalden and Diemand (2013a); *anderhalden2013erratum; Ishiyama (2014); Polisensky and Ricotti (2015); Ogiya and Hahn (2018); Angulo et al. (2017). This profile is stable; it does not relax to the NFW profile, at least in the absence of halo mergers.

The first halos subsequently merge to produce successively larger halos, and Refs. Ogiya et al. (2016); Angulo et al. (2017) find that these mergers gradually drive the halos’ density profiles toward the NFW form. As these works argue, the shallowing of the inner density profile is likely a consequence of violent relaxation Lynden-Bell (1967) during merger events. However, there are multiple lines of evidence suggesting that such relaxation does not erase the memory of prior states. For nearly equal-mass mergers, the density profile of the merger remnant depends sensitively on the profiles of its progenitors Polisensky and Ricotti (2015); Angulo et al. (2017). Successive mergers either raise or leave unaltered a halo’s characteristic density Ogiya et al. (2016); Angulo et al. (2017); Drakos et al. ; since the first halos are the densest, this trend preserves them. Finally, for highly unequal-mass mergers, the smaller halo is generally expected to survive as a subhalo of the larger Berezinsky et al. (2008). The smaller halo’s central density profile may even be mostly unaffected by the merger Hayashi et al. (2003); Penarrubia et al. (2010); Ogiya et al. (2019).

In this way, it is broadly plausible that the dense central regions of the first halos survive the hierarchical clustering process. These halos, forming during a denser epoch, should be the densest dark matter objects in the Universe, and such density leads to observational prospects. Signals from dark matter annihilation are dramatically enhanced by the high density within these halos (e.g., Refs. Diemand et al. (2007); Strigari et al. (2007); Pieri et al. (2008); Jeltema and Profumo (2008); Bartels and Ando (2015); Ando et al. ). Gravitational signatures, whether through microlensing (e.g., Refs.  Boden et al. (1998); Dominik and Sahu (2000); Erickcek and Law (2011); Van Tilburg et al. (2018)), timing delays (e.g., Ref. Siegel et al. (2007)), or stellar dynamics (e.g., Refs. Erkal and Belokurov (2015); Erkal et al. (2016); Buschmann et al. (2018)), are also enhanced. Meanwhile, since these halos form by direct collapse from overdense patches, they carry sensitive information about the primordial density field on scales that are inaccessible to other probes. Our goal is to develop a model that extracts this information; we aim to connect the properties of the first halos to the statistics of the primordial matter density field.

The matter density field is intimately connected to some of the most fundamental questions of cosmology. The power spectrum $\mathcal{P}(k)$, which quantifies the power in density fluctuations at scale wave number $k$, has been precisely measured at scales above approximately $1$ Mpc using the cosmic microwave background Hlozek et al. (2012) and the Lyman-$\alpha$ forest Bird et al. (2011). However, with a few exceptions Josan et al. (2009); Chluba et al. (2012), it is largely unconstrained at submegaparsec scales. Since density fluctuations are thought to have been seeded during inflation, $\mathcal{P}(k)$ serves as a valuable probe of inflationary models Lidsey et al. (1997). The power spectrum is also sensitive to the thermal history of the Universe after inflation. An early matter-dominated era (EMDE), driven by an unstable heavy relic, would amplify small-scale density fluctuations Erickcek and Sigurdson (2011); Barenboim and Rasero (2014); Fan et al. (2014); Erickcek (2015), as would a period of domination by a fast-rolling scalar field Redmond et al. (2018). The power spectrum $\mathcal{P}(k)$ thereby supplies one of the few windows into the Universe prior to big bang nucleosynthesis. Separately, its sensitivity to the free-streaming scale makes it a probe of the nature of dark matter.

Numerous previous works have explored the prospects of using the first dark matter halos to probe the small-scale primordial power spectrum Josan and Green (2010); Bringmann et al. (2012); Li et al. (2012); Yang et al. (2013a, b); Yang (2014); Clark et al. (2016); *clark2017erratumII; Yang and Qin (2017); Nakama et al. (2018). However, most treatments assumed that these halos, if they form sufficiently early, possess a particularly compact $\rho\propto r^{-9/4}$ density profile Ricotti and Gould (2009) derived from self-similar theory Bertschinger (1985). References Gosenca et al. (2017) and Delos et al. (2018a) (Paper I) used $N$-body simulations to show that this profile does not arise from a realistic formation scenario. In Ref. Delos et al. (2018b), hereafter Paper II, we showed that, despite possessing shallower $\rho\propto r^{-3/2}$ density profiles, the first halos can still supply competitive constraints on the power spectrum through nondetection of their observable signals. In that work, we used scaling arguments to model the population of the first halos given a particular family of power spectra. Each peak in the density field was mapped to a collapsed halo at later time.

Our present work represents a natural extension of that model to arbitrary power spectra. It is based on the notion that the density profile of a halo forming by direct collapse is uniquely related to the properties of its precursor density peak. Numerous prior works have explored the problem of explaining a halo’s density profile in terms of the structure of the density peak whence it collapsed. In work that pioneered the so-called spherical infall model, Refs. Gunn and Gott III (1972); Gott (1975) approximated the density profile of a collapsed halo by employing the simplifying assumption that particle orbits are unaltered after accretion. Subsequent works extended this model by including the contraction of orbiting material due to new accretion Gunn (1977); Fillmore and Goldreich (1984); Bertschinger (1985); Hoffman and Shaham (1985); Ryden and Gunn (1987); Zaroubi and Hoffman (1993); Łokas and Hoffman (2000); Dalal et al. , relaxing the assumption of spherical symmetry Ryden (1993); Lithwick and Dalal (2011), and modeling nonradial motions White and Zaritsky (1992); Nusser (2001); Ascasibar et al. (2004); Lu et al. (2006); Ascasibar et al. (2007); Zukin and Bertschinger (2010a, b). Because of the difficult, nonlinear nature of the matter infall problem, every treatment employs simplifying assumptions. Exact solutions only exist for the self-similar case Fillmore and Goldreich (1984); Bertschinger (1985); Ryden (1993); Lithwick and Dalal (2011), in which the primordial mass excess is a power law in radius. More general treatments employ Ansätze related to angular momenta and orbital contraction. As an alternative to tracking each orbit, other works have employed an Ansatz related to virialization Salvador-Solé et al. (2012); Juan et al. (2014).

Whereas spherical infall models relate halo density profiles to the precursor mass distribution, a complementary paradigm empirically studies the distribution of halo density profiles as a function of cosmology and redshift. In this paradigm, cosmological $N$-body simulations are used to tune a parametric model that describes the density profile of a halo as a function of its mass, and models built in this paradigm are known as concentration-mass relations.111The concentration is a parameter in the NFW profile (and extended to other profiles) describing how centrally distributed the halo’s mass is. The distribution of halo masses can be subsequently obtained using Press-Schechter theory Press and Schechter (1974); Sheth et al. (2001); Sheth and Tormen (2002), so these models can predict the full population of halos and their density profiles. Concentration-mass relations have been studied extensively. The simplest models describe the halo distribution in a particular cosmological scenario and at a particular time of interest Avila-Reese et al. (1999); Jing (2000); Colin et al. (2004); Dolag et al. (2004); Avila-Reese et al. (2005); Neto et al. (2007); Duffy et al. (2008); *duffy2011erratum; Gao et al. (2008); Macciò et al. (2008); Klypin et al. (2011); Muñoz-Cuartas et al. (2011); Bhattacharya et al. (2013); Dutton and Macciò (2014); Heitmann et al. (2015); Klypin et al. (2016); Hellwing et al. (2016); Angel et al. (2016); Child et al. (2018). Other works have framed a halo’s density profile in terms of its age or assembly history Navarro et al. (1996, 1997); Bullock et al. (2001); Eke et al. (2001); Wechsler et al. (2002); Zhao et al. (2003, 2009); Giocoli et al. (2012); Ludlow et al. (2013, 2014); Bosch et al. (2014); Correa et al. (2015); Ludlow et al. (2016), and progress has been made in isolating the physical variables most directly relevant for predicting density profiles Prada et al. (2012); Okoli and Afshordi (2016); Diemer and Kravtsov (2015); Diemer and Joyce (2019). Nevertheless, due to their empirical nature, concentration-mass relations do not readily extend beyond the cosmological scenarios, times, and halo-mass ranges over which they are tuned.

Broadly, spherical infall models attempt to explain the structures of halos from first principles, while concentration-mass relations endeavor to predict these structures pragmatically. Our analysis constitutes a hybrid between these two procedures that is specialized to the first halos. Forming by direct collapse, these halos are well suited to the spherical infall description. At the smallest radii, we use ellipsoidal collapse arguments Bond and Myers (1996a); Sheth et al. (2001); Sheth and Tormen (2002) to predict the coefficient of the $\rho\propto r^{-3/2}$ inner profile. Beyond the inner asymptote, we employ the simplest spherical infall models to predict the larger profile, parametrized by the maximum circular velocity $v_{\mathrm{max}}$ and the radius $r_{\mathrm{max}}$ at which it is attained Ascasibar and Gottlöber (2008). By building from such first-principles descriptions, our predictive models are valid in any cosmological scenario; we demonstrate this by validating and tuning the models using six high-resolution cosmological $N$-body simulations carried out in wildly disparate cosmological scenarios. Our models nominally predict a halo’s density profile from the density peak whence it collapsed. However, modulo the influence of halo mergers (which we discuss), the statistics of peaks Bardeen et al. (1986) may be applied to thereby predict the full halo population at a given time.

This paper is structured as follows. In Sec. II, we detail our simulations and the procedure we use to connect halos with density peaks. Section III develops a model that predicts the small-radius asymptote $\rho\propto r^{-3/2}$ of a halo’s density profile, while Sec. IV compares models that predict the profile at larger radii. In Sec. V, we discuss prospects for predicting whole populations of halos, including discussion of halo mergers in Sec. V.2. Section VI concludes and discusses avenues for future work. Appendixes A and B further detail how we extract simulation data. Finally, in Appendix C, we present a procedure to directly sample halos from a power spectrum using our model.

II Simulations

We first build a halo catalogue on which to test our model. For this purpose, we simulated six different simulation boxes drawn from different initial power spectra. These power spectra are shown in Fig. 1. Three of these power spectra are constructed as “spikes” centered at the scale $k_{s}=6.8$ kpc*-1* with the form

[TABLE]

for different values of $w$. These spectra are primarily intended as artificial test beds for halo formation, but they do have qualitative motivations in inflationary phenomenology Salopek et al. (1989); Starobinskij (1992); Ivanov et al. (1994); Randall et al. (1996); Starobinsky (1998); Martin et al. (2000); Chung et al. (2000); Barnaby and Huang (2009); Barnaby (2010); Bugaev and Klimai (2011). A fourth power spectrum represents the impact of an EMDE with reheat temperature $T_{\mathrm{RH}}=100$ MeV and ratio $k_{\mathrm{cut}}/k_{\mathrm{RH}}=20$ between the free-streaming cutoff and the largest scale affected by the EMDE (see Ref. Erickcek (2015)). Finally, the last two power spectra include only the free-streaming cutoffs associated with cold dark matter with mass $m_{\chi}=100$ GeV and kinetic decoupling temperature $T_{\mathrm{kd}}=33$ MeV (corresponding to a typical weakly interacting massive particle Jungman et al. (1996)) and warm dark matter (WDM)222For WDM, we also take the dark matter particle to be fermionic with 2 degrees of freedom (so, $g_{X}=1.5$ in Ref. Bode et al. (2001)). with mass $m_{\chi}=3.5$ keV (close to lower bounds from the Lyman-$\alpha$ forest Viel et al. (2013)), respectively. The particular parameter choices for these power spectra are intended to represent very different cosmological scenarios so that we can test the broad applicability of our models.

We generate the three spiked power spectra from Eq. (1) at a starting redshift of $z=10^{4}$, normalizing $\mathcal{A}=2\times 10^{-7}$ at this time in order to effect abundant halo formation by $z\sim 100$. To generate the latter three power spectra, we begin with a primordial curvature power spectrum of amplitude $\mathcal{A}_{s}=2.142\times 10^{-9}$ and spectral index $n_{s}=0.9667$ Ade et al. (2016). We then use the Boltzmann solver Camb Sources Challinor and Lewis (2011); Lewis and Challinor (2007) to produce a matter power spectrum at $z=500$ in a scenario with no thermal dark matter motion. Next, we apply transfer functions from Refs. Erickcek (2015), Green et al. (2004), and Bode et al. (2001) to create, respectively, the EMDE, CDM, and WDM spectra. Finally, we use linear theory Hu and Sugiyama (1996) to evolve these power spectra to their initial redshifts. Note that we include radiation in our simulations, and the linear evolution accounts for this; see Paper II for details.333Modes that were subhorizon during an EMDE grow in an altered way during radiation domination Erickcek (2015), and we modify the linear growth function for this scenario accordingly. We also include a cosmological constant, although it is of minimal relevance.

The starting redshift $z_{\mathrm{init}}$, ending redshift $z_{\mathrm{final}}$, and periodic box size of each simulation are listed in Table 1. Also listed are the rms density variance $\sigma_{0}$ and comoving characteristic correlation length $R_{*}$ associated with each power spectrum in linear theory. These spectral parameters are defined as Bardeen et al. (1986)

[TABLE]

To illustrate the significance of these quantities, Fig. 1 also shows the correspondingly rescaled power spectra. This rescaling is useful because it factors out the scale differences between the spectra, leaving only their shapes. For this reason, we will use these spectral parameters as our units. The mass unit is $m_{\mathrm{unit}}\equiv\bar{\rho}_{0}R_{*}^{3}$, where $\bar{\rho}_{0}$ is the background matter density today. Meanwhile, the physical distance unit is $r_{\mathrm{unit}}\equiv[a/\sigma_{0}(a)]R_{*}$, where $\sigma_{0}(a)$ is evaluated using linear theory during matter domination (so $\sigma_{0}\propto a$) and $a=1$ today. For each simulation, the ending redshift is tuned so that $\sigma_{0}\sim 3$ (in linear theory) at simulation termination, and the comoving box size is tuned to be roughly 30 to 80 times $R_{*}$. The initial redshift is chosen so that the largest fractional density excesses ${\delta\equiv(\rho-\bar{\rho})/\bar{\rho}}$ are of order $0.2$ or smaller.

For each simulation, we draw a random density field from the matter power spectrum at $z_{\mathrm{init}}$ and generate initial conditions using the Zel’dovich approximation modified to account for radiation as described in Paper II.444For the EMDE, the methods in Paper II were adapted to the post-EMDE growth function. Finally, we carry out the simulations using a version of the cosmological simulation code Gadget-2 Springel (2005); Springel et al. (2001) that we modified to include the effects of radiation (see Paper II). All simulations employed $1024^{3}$ particles and a comoving force softening length set at $3\%$ of the initial interparticle spacing (with forces becoming non-Newtonian at 2.8 times this length).

Figure 2 shows the initial and final density fields for three of these simulations. This sample illustrates the characteristic differences between the density fields and halo populations that result from the different power spectra. The $w=0.1$ spectrum produces fluctuations only on a characteristic scale, and the resulting halo population is relatively uniform in size and separation. Meanwhile, the $w=0.5$ spectrum produces fluctuations on a wider range of scales, yielding very small halos as well as larger halo clusters and voids. Finally, the 100 GeV power spectrum includes fluctuations at scales up to and exceeding the box size.555In fact, the entire boxes of the 100 GeV and 3.5 keV simulations should be collapsing by the final redshift, a fact that cannot be reflected in these simulations. Thus, we do not expect these two simulations to accurately capture the large-scale dynamics; for example, halo mergers are probably under-represented. We focus on understanding the halos that form by direct collapse from peaks in the density field, so we do not expect this shortcoming to impact our main results. Accordingly, the bulk of the halos in the box reside within a few clusters, and some halos are growing exceedingly large.

At the final redshift of each simulation, we use the Rockstar halo finder Behroozi et al. (2013a) to identify every halo with mass larger than a cutoff $M_{\mathrm{cut}}=(4\pi/3)k_{\mathrm{max}}^{-3}\bar{\rho}_{0}$, where $k_{\mathrm{max}}$ is the wave number $k$ at which $\mathcal{P}(k)$ is maximized. This cutoff is intended to exclude most of the artificial fragments Angulo et al. (2013); Lovell et al. (2014) that are normally expected to appear in simulations with a small-scale cutoff in density fluctuations. Roughly speaking, fluctuations at the scale $k_{\mathrm{max}}$ are the first to collapse, so we do not expect to find many real halos at smaller scales. Table 2 lists the number $N_{\mathrm{halo}}$ of such halos found in each simulation.

Next, we wish to associate each halo with the density peak whence it collapsed. To do this, we use Rockstar to generate halo catalogues at intervals of roughly $5\%$ in the scale factor, and we use the Consistent Trees merger tree code Behroozi et al. (2013b) to improve the consistency of halo tracking through time. Starting at the final redshift, we trace each halo back to its formation time, defined as the time when its mass first exceeded $M_{\mathrm{cut}}$. In case of a merger, we follow the larger progenitor. We then identify all particles in the halo at its formation time, find their positions in the initial particle grid at $z_{\mathrm{init}}$, and compute the center of mass $\bm{x}_{\mathrm{CM}}$ and rms spread $r_{\mathrm{CM}}$ of these positions. Finally, we consider the initial grid $\delta(\bm{x})$ of fractional density contrasts and find the largest value of $\delta(\bm{x})$ within a sphere of radius $r_{\mathrm{CM}}$ about $\bm{x}_{CM}$. If this maximum is also a local maximum in $\delta(\bm{x})$ (a condition that could fail if the maximum lies at the sphere’s boundary), we take it to be the initial density peak for this halo. Otherwise, we discard the halo from further consideration, suspecting it to be an artificial fragment.666We explored the alternative procedure of following the density gradient to a local maximum, but the only consequence was the introduction of unphysically distant halo-peak associations. Moreover, when this procedure associates multiple halos in the final box with the same initial peak, we discard all but the most massive halo. Table 2 lists both the number $N_{\mathrm{peak}}$ of peaks in the initial density field and the number $N_{\mathrm{match}}$ of peak-halo identifications made through this procedure.

To test our models, we must still make further cuts to the halo-peak population. For the small-radius asymptote of the density profile (Sec. III), we consider the subset of our peak-matched halos that have well-resolved inner density profiles777In particular, we require that at least 100 particles reside within the simulation softening length. Otherwise, the central density is small, usually implying that the density profile was not centered on the halo’s cusp, which is only true of a handful of halos not found by Rockstar but filled in by Consistent Trees. Due to the minuscule fraction of halos affected, we did not make further efforts to find the true center.. The number $N_{A}$ of such halos in each simulation is listed in Table 2. Moreover, Ref. Ogiya et al. (2016) found that successive major mergers alter the inner asymptote of a halo’s density profile, and we confirm this effect in Sec. V.2. We aim to model the initial density profile of a halo before it undergoes any disruptive dynamics, leaving a treatment of mergers for future work.888We prefer to exclude the impact of mergers in this work because mergers are a continuing process, so their impact is sensitive to the arbitrarily chosen simulation termination redshift. Thus, we also restrict our sample to halos that underwent fewer than three major mergers, which we define to be mergers between two halos with a mass ratio smaller than 3. The number $N_{A}^{<3\mathrm{MM}}$ of these halos is also listed in Table 2. The particular threshold of three major mergers is chosen as a compromise between minimizing the impact of mergers on our results and maximizing the sample of halos from which our results are drawn. As we will see in Sec. V.2, the first two major mergers only marginally alter halo density profiles.

Meanwhile, for the larger-radius density profile (Sec. IV), we consider the subset of our peak-matched halos that have $r_{\mathrm{max}}$, the radius at which the circular velocity is maximized, inside their virial radius, defined as the radius $r_{\mathrm{vir}}$ enclosing average density 200 times the background.999The requirement $r_{\mathrm{max}}<r_{\mathrm{vir}}$ culls a significant fraction of the halo population in the 100 GeV simulation, but this is not a serious concern. It turns out that the models we discuss predict very large $r_{\mathrm{max}}$ for these halos as well, and $r_{\mathrm{max}}>r_{\mathrm{vir}}$ only implies that the halo’s outermost profile has not yet steepened to the point that $\mathrm{d}\ln M/\mathrm{d}\ln r<1$. The number $N_{\mathrm{r}}$ of these halos is listed for each simulation in Table 2. Mergers only minimally alter $r_{\mathrm{max}}$ (see Sec. V.2), so when testing models for $r_{\mathrm{max}}$, it is not necessary to restrict to halos that underwent few mergers. However, mergers can significantly alter the mass $M(r_{\mathrm{max}})$ enclosed within $r_{\mathrm{max}}$ (or equivalently $v_{\mathrm{max}}$), so to test our mass predictions, we restrict the sample to the subset $N_{\mathrm{r}}^{<3\mathrm{MM}}$ of halos that underwent fewer than three major mergers. This number is also listed in Table 2. Altogether, it is evident from Table 2 that only a small fraction of the full peak and halo populations is used to test our models. Section V explores the impact of these restrictions by returning to the full peak and halo populations.

We now have a catalogue that matches peaks in an initial density field to collapsed halos at a much later time. All that remains is to collect the halo density profiles and the parameters of each peak, and we detail these processes in Appendixes A and B, respectively. With these data, we are now prepared to test any model relating the structure of a dark matter halo to its precursor density peak. In the following sections, we develop such a model.

III Density profile at small radii

We first study the coefficients $A$ of the $\rho=Ar^{-3/2}$ asymptotes of the density profiles of the first halos. It is of prime importance to accurately predict the density profiles in this regime because these radii source the bulk of the prospective signal from dark matter annihilation, and this remains true even if these halos relax toward ${\rho\propto r^{-1}}$ profiles due to mergers. Moreover, there is another reason to study small radii separately: the ${\rho\propto r^{-3/2}}$ inner density profile is established almost immediately after collapse. To illustrate this fact, we simulate101010This simulation employed about 9 million particles in a comoving vacuum-bounded sphere of radius 1.5 kpc. The starting redshift was $z=10^{6}$. the collapse of the isolated density peak shown in Fig. 3. This peak represents a typical $3\sigma$ peak drawn from the $w=0.3$ power spectrum using the statistics of peaks as described in Ref. Bardeen et al. (1986). The shallowing of this profile toward $r=0$ is associated with the absence of small-scale fluctuations, and this is the single feature common to density peaks drawn from all of the power spectra we are studying. Indeed, it has been suggested that this feature is responsible for the development of the ${\rho\propto r^{-3/2}}$ profile Ogiya and Hahn (2018).

Figure 4 shows the halo resulting from this initial peak in the moments after collapse. We also compute the expected collapse time from the initial peak using both the spherical collapse model, $\delta(a_{\mathrm{sc}})=1.686$, and the ellipsoidal collapse model described in Ref. Sheth et al. (2001). The redshifts $z_{\mathrm{sc}}$ and $z_{\mathrm{ec}}$ of spherical and ellipsoidal collapse, respectively, are indicated in Fig. 4. The critical observation is that the $\rho\propto r^{-3/2}$ asymptote of the density profile develops almost immediately after collapse, matching closely the late-time density profile as early as $a\simeq 1.15a_{\mathrm{ec}}$. Subsequently, the profile grows outward in radius alone, gradually steepening as accretion of new material slows. The interval $\Delta a/a\simeq 0.15$ over which the profile develops is not arbitrary. The characteristic dynamical time $\Delta t=\sqrt{3\pi/(16G\rho)}$ Binney and Tremaine (1987) of a virialized region with density $200$ times the background is $\Delta t=0.16/H$, where $H$ is the Hubble rate. This interval corresponds to $\Delta a/a=0.16$, implying that the inner density profile develops over a single dynamical time interval.

III.1 Spherical collapse

The rapid development of a halo’s inner $\rho=Ar^{-3/2}$ profile implies that the coefficient $A$ can only be influenced by the immediate neighborhood of the precursor density peak. To make this argument precise, we define111111After this section, we will omit the tilde, and any quantity related to the linear density field should be understood to be evaluated on the scaled linear density field $\tilde{\delta}(\bm{x})$ unless otherwise specified.

[TABLE]

where $\delta(\bm{x},a)$ is evaluated using linear theory during matter domination (when $\delta\propto a$). Here, $\bm{x}$ is a comoving coordinate, and we define $a=1$ today. We now claim that the inner $\rho\propto r^{-3/2}$ asymptote is only sensitive to the local properties $\tilde{\delta}$ and $\partial_{i}\partial_{j}\tilde{\delta}$ of the precursor density peak. If we neglect deviations from spherical symmetry, the peak reduces to two parameters: its amplitude $\tilde{\delta}$ and curvature $|\nabla^{2}\tilde{\delta}|$. During matter domination, spherical collapse theory implies that $a_{\mathrm{sc}}\propto\tilde{\delta}^{-1}$. Meanwhile, $|\nabla^{2}\tilde{\delta}|$ defines a comoving length scale $q_{\mathrm{pk}}\equiv(\tilde{\delta}/|\nabla^{2}\tilde{\delta}|)^{1/2}$ associated with the peak. At the time of collapse, this comoving length corresponds to the physical length scale $a_{\mathrm{sc}}q_{\mathrm{pk}}\propto(\tilde{\delta}|\nabla^{2}\tilde{\delta}|)^{-1/2}$. There is also a physical density scale $\bar{\rho}_{0}a_{\mathrm{sc}}^{-3}\propto\bar{\rho}_{0}\tilde{\delta}^{3}$, i.e., the cosmological background density at collapse. If these are the only scales, then up to a constant coefficient, there is a unique prediction for the coefficient $A$ of the $\rho=Ar^{-3/2}$ asymptote of the density profile,

[TABLE]

where $\alpha$ is a proportionality constant. We omit the tildes in Eq. (5), but all quantities related to the linear density field are understood to be evaluated on $\tilde{\delta}(\bm{x})$. This convention applies to the remainder of this work. The parameter $\alpha$ must be calibrated by simulations, but this calibration is only necessary once; it should be the same for any power spectrum.

We test this model on the $N_{A}^{<3\mathrm{MM}}$ peak-matched halos in each simulation that have well-resolved inner density profiles and underwent fewer than three major mergers (see Sec. II). For these halos, Fig. 5 plots the left-hand side against the right-hand side of Eq. (5) in order to test the model. Evidently, our model works well for how simple it is, predicting the asymptotes with reasonable success in all six simulations. There does, however, appear to be a correlated effect wherein halos with the densest predicted asymptotes tend to exceed that prediction, and vice versa at the less dense end. In fact, this effect is caused by the assumption of spherical collapse. Peaks of smaller amplitude tend to be less spherical, thereby collapsing later and forming less dense halos than their amplitude would suggest. We next account for this effect.

III.2 Ellipsoidal collapse

Equation (5) neglects deviations from spherical symmetry. However, it can be immediately generalized in the following way. Since $A\propto a_{\mathrm{sc}}^{-3/2}$, we can use the theory of ellipsoidal collapse Bond and Myers (1996a); Sheth et al. (2001); Sheth and Tormen (2002) to predict how the three-dimensional shape of a peak alters its collapse time. In particular, ellipsoidal collapse occurs later than spherical collapse by a factor $f_{\mathrm{ec}}(e,p)\equiv a_{\mathrm{ec}}/a_{\mathrm{sc}}$, which is a function of the ellipticity $e$ and prolateness $p$ of the gravitational potential $\phi$ in the vicinity of the peak (see Appendix B for definitions). To compute $f_{\mathrm{ec}}$, we use the approximation Sheth et al. (2001)

[TABLE]

Accounting for ellipsoidal collapse,

[TABLE]

where $\alpha$ is a proportionality constant.

Figure 5 compares the predicted asymptotes in the ellipsoidal collapse model to their measured values. The correlated scatter associated with the spherical model is no longer apparent, confirming that the spherical collapse assumption was the source. Figure 6 compares the overall scatter between these two models; the ellipsoidal collapse model is clearly superior. We remark that not all of the scatter depicted is physical; there is numerical noise in our simulated density profiles and finite density grids. Nevertheless, there are still clear sources of physical scatter. Absent an understanding of why the $\rho\propto r^{-3/2}$ profile develops, we constructed our model using purely scaling arguments, and these arguments can break down in two ways. First, the inner profile does not develop instantly, so it is sensitive to a small region about the density peak instead of only its immediate neighborhood. Also, as depicted in Fig. 4, the inner asymptote is not completely static but instead grows marginally after it develops, implying that it exhibits some sensitivity to the larger density field. Second, we only accounted for nonspherical peak shapes by altering the time of collapse. The sensitivity of the asymptotic coefficient $A$ to the shape of the peak may be more complicated.

Nevertheless, it is a significant success to obtain a mere approximately $30\%$ scatter across such a broad range of cosmologies from a model as simple as the one presented here. The proportionality coefficients and their statistics for our models are summarized in Table 3. Thus calibrated, the model developed in this section may be employed to predict the inner asymptotes of the first halos’ density profiles in any cosmological scenario.

IV Density profile at large radii

We next study the density profile beyond the inner $\rho\propto r^{-3/2}$ asymptote. These larger radii, containing the bulk of the halo mass, are relevant to gravitational lensing signatures as well as to the dynamical evolution of the halo through mergers and tidal stripping. In this section, we discuss and validate models that can predict the density profile at large radii. We focus on three physical models:

1. (1)

a “turnaround” model, tracing back to Ref. Gunn and Gott III (1972), in which each mass shell freezes at a fraction of its turnaround radius, or the radius of first apocenter; 2. (2)

a “contraction” model, expressed most concisely by Ref. Dalal et al. , that accounts for contraction of halo particle orbits due to subsequent infall; 3. (3)

a “virialization” model, put forward by Ref. Salvador-Solé et al. (2012), in which the final radius of a mass shell is determined by enforcing that the shell’s enclosed energy be distributed according to the virial theorem.

All of these models rely on the assumption that material is accreted gradually. Hence, they are not applicable to the density profiles at small radii, and none of them predicts the $\rho\propto r^{-3/2}$ asymptote. We will explore in Sec. IV.4 precisely where these models are accurate. In the meantime, we validate and tune the models by using them to predict the radius $r_{\mathrm{max}}$ at which the halo’s circular velocity is maximized along with the mass $M(r_{\mathrm{max}})$ within this radius. The maximum circular velocity itself follows as $v_{\mathrm{max}}=\sqrt{GM(r_{\mathrm{max}})/r_{\mathrm{max}}}$.

The radius $r_{\mathrm{max}}$ characterizes where the profile bends away from its inner asymptote.121212For example, if $\rho(r)=\rho_{s}(r/r_{s})^{-3/2}(1+r/r_{s})^{-3/2}$ Delos et al. (2018a); Moore et al. (1999) (where $\rho_{s}$ and $r_{s}$ are scale parameters), $r_{\mathrm{max}}=1.055r_{s}$. It is more common in the literature to study the scale radius $r_{s}$, often defined to be where $\mathrm{d}\ln\rho/\mathrm{d}\ln r=-2$, instead of $r_{\mathrm{max}}$, which is where $\mathrm{d}\ln M/\mathrm{d}\ln r=1$. However, we favor $r_{\mathrm{max}}$ over $r_{s}$ for several reasons. First, $r_{\mathrm{max}}$ can be read from a simulation more robustly than $r_{s}$, since the mass profile is less noisy than the density profile. Second, $r_{\mathrm{max}}$ turns out to be cleaner to predict from the linear density field. Finally, since $r_{\mathrm{max}}$ characterizes the total mass rather than the local density, it is the more relevant quantity for understanding halo dynamics Ascasibar and Gottlöber (2008) and for predicting gravitational lensing signatures. The combination of $A$ and $r_{\mathrm{max}}$ serves to mostly characterize the density profile of a halo, and the mass $M(r_{\mathrm{max}})$ (or velocity $v_{\mathrm{max}}$) supplies an additional constraint.

IV.1 Turnaround

To develop a model that can predict the full density profile of a collapsed halo, we must relate this profile to the spherically averaged fractional density excess profile $\delta(q)$, where $q$ is the comoving radius, about the corresponding peak in the linear density field. Following the convention established in Sec. III, we define $\delta(q)\equiv\delta(q,a)/a$ evaluated in linear theory during matter domination (with $a=1$ today). We will also use

[TABLE]

the fractional enclosed mass excess defined under the same convention. As in Ref. Gunn and Gott III (1972), we consider the simplified spherical infall model in which each spherical shell freezes at a fixed fraction of its turnaround radius. A mass shell initially at comoving radius $q$ turns around at physical radius $r_{\mathrm{ta}}=(3/5)q/\Delta(q)$. Hence, the final physical radius $r$ of this shell is

[TABLE]

where $\beta$ is a proportionality constant to be measured in simulations. But in this model, mass shells never cross, so the mass enclosed within this shell, now at $r$, is still

[TABLE]

(at zeroth order in $\Delta$) with $\beta_{M}=1$. This equation gives the mass profile $M(r)$ of the collapsed halo if $q$ is obtained from $r$ by inverting Eq. (9).

To find $r_{\mathrm{max}}$ (up to the proportionality constant $\beta$), we may maximize $M(r)/r$. Alternatively, we can write

[TABLE]

where we define131313We define $\epsilon(q)$ as a generalization of the index $\epsilon$ of $\delta M/M\propto M^{-\epsilon}$ in the self-similar theory Fillmore and Goldreich (1984). Consequently, Eq. (11) has exactly the same form as its analogue in the self-similar theory.

[TABLE]

In this case, $r_{\mathrm{max}}$ is obtained as the solution to $\mathrm{d}\ln M/\mathrm{d}\ln r=1$ or, equivalently, ${\epsilon(q_{\mathrm{max}})=2/3}$ with $r_{\mathrm{max}}$ computed from $q_{\mathrm{max}}$ using Eq. (9).

This model, at first glance, seems to be far divorced from a realistic description. However, it turns out to be a reasonable approximation of late accretion. At late times, halo density profiles are stable in time; see, e.g., Papers I and II. This observation is explained by noting that once the halo is established and accretion begins to slow, newly accreted mass only contributes significantly to the outskirts of the halo; there is too little new matter, spending too small a fraction of its orbital period at small radii, to substantially raise the interior density. Since the density profile is static, a newly accreted particle settles into a stable orbit with time-averaged radius proportional to its orbital apocenter, the turnaround radius.141414The orbital apocenter actually decays over time, an effect for which the contraction model accounts. The final radius in this model is to be understood as that orbital average.

We test this model on the $N_{\mathrm{r}}$ peak-matched halos for which $r_{\mathrm{max}}<r_{\mathrm{vir}}$ (see Sec. II). In the top panel of Fig. 7, we plot the measured $r_{\mathrm{max}}$ against the prediction from this model. The model appears to work well for the bulk of the halos. However, there is a significant population of halos, especially coming from the broader 100 GeV, 3.5 keV, EMDE, and $w=0.5$ power spectra, for which the predicted $r_{\mathrm{max}}$ is much larger than the measured value. In fact, in these simulations, there are many peaks for which the mass shell with $\epsilon=2/3$ has not yet accreted onto the halo by the final redshift. In these cases, it makes no sense to use the properties of this shell to predict $r_{\mathrm{max}}$.

A minimal correction, for these halos, is to instead relate $r_{\mathrm{max}}$ to the outer boundary of the halo. In the spherical collapse model, a mass shell virializes when it falls to half of its turnaround radius; at this point, its energy is distributed according to the virial theorem. The halo’s outer boundary may be taken to be the physical radius of the last shell to virialize in this way, and we define $q_{\mathrm{v}}$ to be the Lagrangian radius of this shell. If we are considering the halo population at scale factor $a$, then $q_{\mathrm{v}}$ is the smallest $q$ that satisfies $\Delta(q,a)=a\Delta(q)\leq\delta_{v}$, where $\delta_{v}=1.583$. If $r_{\mathrm{max}}$ is proportional to the halo boundary so defined, then

[TABLE]

where $\beta^{\prime}$ is another proportionality constant, and this model is understood to apply only to the halos of which the $\epsilon=2/3$ mass shells have not yet virialized, i.e., $q_{\mathrm{max}}>q_{\mathrm{v}}$. In Fig. 8, we plot the turnaround model for only those halos of which the $\epsilon=2/3$ shells have virialized, separately plotting the $q_{\mathrm{v}}$ model for the remainder of the halos. Evidently, the halos of which the $\epsilon=2/3$ shells had not yet virialized were indeed the population of which the $r_{\mathrm{max}}$ was severely overpredicted, and the $\epsilon=2/3$ model exhibits significantly less scatter with them excluded. Meanwhile, using $q_{\mathrm{v}}$ to predict $r_{\mathrm{max}}$ for these halos works reasonably well, although the scatter here is significantly larger. The top panel of Fig. 9 depicts these differences in scatter more transparently.

Finally, we test how well this model predicts $M(r_{\mathrm{max}})$ in similar fashion by using Eq. (10) and allowing $\beta_{M}$ to float. In this case, we employ the $N_{\mathrm{r}}^{<3\mathrm{MM}}$ halos with $r_{\mathrm{max}}<r_{\mathrm{vir}}$ that also underwent fewer than three major mergers (see Sec. II). As above, we set $q=q_{\mathrm{max}}$ when $q_{\mathrm{max}}<q_{\mathrm{v}}$ and $q=q_{\mathrm{v}}$ otherwise. The bottom panel of Fig. 9 shows the scatter in these predictions. For both $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$, Table 3 lists the tuned coefficients $\beta$ and $\beta_{M}$ and their statistics.

IV.2 Contraction

The turnaround model assumed that each new shell freezes at a fixed fraction of its turnaround radius, contributing mass to that radius alone. In reality, a shell contributes mass to a large range of radii. Figure 10 shows the density profiles151515This figure depicts a halo that collapsed from an isolated $3\sigma$ peak drawn from the $w=0.3$ power spectrum. This peak was simulated with about 70 million particles in a comoving vacuum-bounded sphere of radius 1 kpc. laid down by a range of initial mass shells. These shells consist of successive factors of 1.5 in initial radius, so that the lowest shell contains the mass initially in the comoving radius band 0.044 to 0.066 kpc; the second shell contains the mass initially in the band 0.066 to 0.10 kpc; and so on. Notably, each shell has a characteristic radius within the final halo below which it contributes roughly constant density. As Ref. Dalal et al. argues, the constant-density contribution follows from the notion that particles from large-radius mass shells cross the lower radii at such high velocity that their motions are effectively unaccelerated.

Consider a halo particle initially orbiting with apocenter radius $r$ that encloses halo mass $M(r)$. As time goes on, newly accreted shells contribute to the enclosed mass, increasing $M(r)$. In a spherically symmetric and self-similar system, the quantity $\oint v_{r}dr\propto[M(r)r]^{1/2}$ is an adiabatic invariant, implying that as the enclosed mass grows $r$ shrinks as $r\propto 1/M(r)$. While the quantity $Mr$ need not be conserved in a more general picture, Ref. Lithwick and Dalal (2011) found that it remains nearly invariant for each particle even with spherical symmetry relaxed.

For a mass shell with apocenter $r^{\prime}$, let $f(r,r^{\prime})$ be the fraction of the shell’s mass that is within $r$. It is readily seen that the mass enclosed within a shell with Lagrangian radius $q$ and final apocenter $r(q)$ increases by the factor

[TABLE]

due to the contribution of shells with $q^{\prime}>q$. As before, $q_{\mathrm{v}}$ is the Lagrangian radius of the latest shell to virialize; it is the smallest $q$ that satisfies $a\Delta(q)\leq\delta_{v}$, where $a$ is the scale factor at which we wish to characterize the halo population (and $a=1$ today). If orbital apocenters shrink according to $r\propto 1/M$, then

[TABLE]

describes the apocenter of the $q$ shell after contraction, with $\beta=3/5$.

As a simple model, let us assume that each shell contributes density $\rho(r)\propto r^{-s}$ below its apocenter radius and 0 above it, so $f(r,r^{\prime})=(r/r^{\prime})^{3-s}$. With this shell profile, Eq. (14) yields the ordinary differential equation

[TABLE]

with initial condition $X(q_{\mathrm{v}})=1$. This equation is equivalent to the simpler expression in Ref. Dalal et al. , but since it is expressed with respect to the variable $q$, it is more straightforward to integrate over a numerically tabulated peak profile $\epsilon(q)$. The enclosed mass profile after contraction is now

[TABLE]

(to zeroth order in $\Delta$) with $\beta_{M}=1$. Note that $X$ ranges from $1$ (at $q_{\mathrm{v}}$) to $\mathcal{O}(10)$. Using these equations, one may maximize $M(q)/r(q)$ to obtain $q_{\mathrm{max}}$, after which $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$ are obtained immediately. Alternatively,

[TABLE]

with $\mathrm{d}\ln X/\mathrm{d}\ln q$ given in Eq. (16), so one may solve $\mathrm{d}\ln M/\mathrm{d}\ln r=1$ to obtain $q_{\mathrm{max}}$. Finally, since we do not expect to employ the correct form of $f(r,r^{\prime})$, we will allow $\beta$ and $\beta_{M}$ to vary, tuning them to simulations.

We test the contraction model using $s=0$ as well as $s=1$. While Fig. 10 suggests that $s=0$ at small radii, a nonzero choice of $s$ is motivated by noting that for a given mass shell $q$ most of the enclosed mass contributed by higher shells comes from shells just slightly above $q$. Meanwhile, Fig. 10 shows that the density profile contributed by a mass shell does not level off to a constant value until well below its apocenter radius. Therefore, the bulk of the mass contribution comes from shell density profiles $\rho\propto r^{-s}$ with $s>0$. The particular choice $s=1$ is partly arbitrary, but it is roughly the slope of the shell profiles in Fig. 10 slightly below their apocenters.161616The choice $s=1$ also corresponds roughly to the $\rho^{1/2}$ model in Ref. Dalal et al. , since $\rho$ is approximately proportional to $r^{-2}$ at ${r=r_{\mathrm{max}}}$.

Figure 11 plots $r_{\mathrm{max}}$ against its model prediction for both contraction models using the $N_{\mathrm{r}}$ halos with ${r_{\mathrm{max}}<r_{\mathrm{vir}}}$. Both models successfully predict, with some scatter, the values of $r_{\mathrm{max}}$ for the halos in all six simulations. In fact, both models work equally well, as Fig. 9 demonstrates. Moreover, comparing the results of these models to those of the turnaround models, which used the $\epsilon=2/3$ or $q_{\mathrm{v}}$ shells to predict $r_{\mathrm{max}}$, we see that the main difference is that the contraction models can handle all of the halos in a single model requiring just one tuned parameter. Adiabatic contraction can produce a bend in the density profile, and hence an $r_{\mathrm{max}}$, near the halo outskirts because the outskirts are uncontracted while the rest of the halo is contracted. Other than this, there is no significant advantage to the contraction models for predicting $r_{\mathrm{max}}$, as the top panel of Fig. 9 shows.

We also use the $s=0$ and $s=1$ contraction models to predict $M(r_{\mathrm{max}})$ using Eq. (17) with $\beta_{M}$ allowed to float. For this test, we employ the $N_{\mathrm{r}}^{<3\mathrm{MM}}$ halos with $r_{\mathrm{max}}<r_{\mathrm{vir}}$ that underwent fewer than three major mergers. The scatter in these predictions is depicted in the bottom panel of Fig. 9. The tuned values of $\beta$ and $\beta_{M}$ and their statistics are listed in Table 3.

IV.3 Virialization

Reference Salvador-Solé et al. (2012) developed a model for halo density profiles in which a mass shell freezes where its enclosed energy is virialized.202020Note that a model in which a shell freezes where its own energy is virialized is equivalent to the turnaround model above. In this model, the final radius of the $q$ shell is $r(q)=-(3/10)GM(q)^{2}/E(q)$, where $M(q)$ and $E(q)$ are, respectively, the mass and energy enclosed within the $q$ shell in the linear density field. This model may be expressed as

[TABLE]

with $\beta=3/10$ and $M(q)$ defined as in Eq. (10). For this model, $\mathrm{d}\ln M/\mathrm{d}\ln r=3/[6-2\Delta(q)/\kappa(q)]$. As before, $q_{\mathrm{max}}$ is obtained by maximizing $M(q)/r(q)$ or by solving $\mathrm{d}\ln M/\mathrm{d}\ln r=1$. We test this model in a similar fashion, allowing $\beta$ to float and plotting in Fig. 12 the simulated $r_{\mathrm{max}}$ against the model prediction (using the $N_{\mathrm{r}}$ halos with $r_{\mathrm{max}}>r_{\mathrm{vir}}$). Evidently, this model performs similarly to the turnaround model. Again, there is a multitude of halos of which the $q_{\mathrm{max}}$ shells have not yet accreted, leading to dramatically overpredicted values of $r_{\mathrm{max}}$. In this case, the $q_{\mathrm{max}}<q_{\mathrm{v}}$ condition turns out to be too restrictive; it eliminates too many halos. Instead, we require $q_{\mathrm{max}}<q_{\mathrm{ta}}$, where $q_{\mathrm{ta}}$ is Lagrangian radius of the last shell to turn around. In particular, $q_{\mathrm{ta}}$ is the smallest $q$ satisfying $a\Delta(q)\leq\delta_{\mathrm{ta}}$, where $\delta_{\mathrm{ta}}=1.062$ and $a$ is the scale factor at which we are studying the halo population. If $q_{\mathrm{max}}<q_{\mathrm{ta}}$, we use $q_{\mathrm{max}}$ to predict $r_{\mathrm{max}}$ using Eq. (19); otherwise, we use $q_{\mathrm{v}}$ to predict $r_{\mathrm{max}}$ using Eq. (13). We also use this model to predict $M(r_{\mathrm{max}})$ in a similar fashion (again using the $N_{\mathrm{r}}^{<3\mathrm{MM}}$ halos that also underwent fewer than three major mergers), allowing $\beta_{M}$ to float. The scatter and statistics of these predictions212121We remark that the optimal $\beta\simeq 0.042$ is much smaller than the exact value $3/10$ claimed in Ref. Salvador-Solé et al. (2012). are shown in Fig. 9 and Table 3, and we find that this model’s scatter is comparable to that of the turnaround model.

IV.4 Discussion

Table 3 summarizes the $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$ models and the statistics of their simulation-tuned parameters. Evidently, all of these models exhibit similar scatter. For instance, when $q_{\mathrm{max}}<q_{\mathrm{v}}$ or $q_{\mathrm{ta}}$, each model’s predicted $r_{\mathrm{max}}$ exhibits about 30% scatter in the middle half of its predictions. These similarities suggest there may be a statistical floor limiting the precision of all three models equally. We propose two sources of such a floor. The first is artificial discreteness noise in our simulations. As noted in Paper II, there is discreteness noise in each radial bin of the density profile that is significantly larger than Poisson noise and may be associated with the accretion of artificial fragments. Our procedure of averaging radial profiles over a finite time interval (see Appendix A) mitigates but does not eliminate this noise.222222There are other sources of discreteness noise in our analysis, but they are minor. The binning of our density profiles in factors of $1.1$ introduces artificial scatter, but this scatter does not exceed 5% in $r_{\mathrm{max}}$ and is further reduced by our use of interpolation (see Appendix A). There is also noise in the predicted values of $r_{\mathrm{max}}$ because they arise from a finite density grid, but since $q_{\mathrm{max}}$ is typically much larger than a grid cell [so $\delta(q_{\mathrm{max}})$ averages over a large number of cells], the resulting scatter is small.

The other probable source of a statistical floor is physical and owes to a simplifying assumption common to all three models: spherical symmetry. Observe in Fig. 3 that at the initial radius of the $q_{\mathrm{max}}$ shell the density profile of the peak that is depicted—a typical peak drawn from the $w=0.3$ power spectrum—is highly ellipsoidal. It is likely that a substantial fraction of the approximately $30\%$ scatter in $r_{\mathrm{max}}$ results from deviations from spherical symmetry in the initial density peaks. To correct this, models must be employed that move beyond the assumption of spherically symmetric initial conditions, perhaps employing ellipsoidal collapse arguments Bond and Myers (1996a); Sheth et al. (2001); Sheth and Tormen (2002) or drawing from nonspherical self-similar infall theory Ryden (1993); Lithwick and Dalal (2011). The principal difficulty in moving beyond spherical symmetry is that the three-dimensional shape of a peak is different at each radius, an effect for which no model has accounted (to our knowledge). However, the spherical models presented here exhibit just approximately $30\%$ scatter in $r_{\mathrm{max}}$ and approximately $60\%$ scatter in $M(r_{\mathrm{max}})$ over the full range of cosmologies we simulated, a success that should not be understated. As we will see in Sec. V, halo mergers present a larger source of error.

Beyond predicting $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$ (or equivalently $v_{\mathrm{max}}$), we may ask how well the spherical infall models can predict the full density profiles of the first halos. No model successfully predicts $\rho\propto r^{-3/2}$ at small radii232323An $s=3/2$ contraction model predicts $\rho\propto r^{-3/2}$ at small radii, but only in a contrived way; the adiabatic approximation is not valid when this part of the profile develops., and this is not surprising because in the moments after halo collapse the accretion of new material is not adiabatic with respect to the orbits of the already bound material. Nevertheless, we may explore the radial range over which the spherical infall models do accurately predict the density profiles. Figure 13 compares the mass profiles of a small halo sample to their predicted profiles using the median proportionality coefficients in Table 3. Evidently, the turnaround and virialization models do not accurately predict the density profile below $r_{\mathrm{max}}$, although they may succeed at larger radii. The $s=0$ and $s=1$ contraction models, and especially the latter, more accurately capture the shape of the profile at smaller radii. However, as shown in Fig. 10, the true density profile of a mass shell is more complicated than $\rho\propto r^{-s}$, so a more sophisticated model for a shell’s density profile should yield still better results.

In this section, we presented models that can predict the outer portions of the density profiles of the first halos in any cosmology. While the models themselves are not new, our calibrations enable their use as predictive tools. Moreover, our simulations demonstrate the universality of these models by showing that they succeed across wildly dissimilar cosmological scenarios.

V Predicting halo populations

So far, we have studied the relationship between a density peak and its resulting halo. However, our ultimate goal is to study the larger populations. Recall from Sec. II that not all peaks matched to halos in our simulations and not all halos matched to peaks. These discrepancies are reflected in Table 2 and may arise from physical processes, such as halo mergers, or numerical artifacts in our simulations. We have shown that if we know that a peak developed into a halo that persisted to some later redshift then we can predict the properties of that halo from the properties of the peak. Ultimately, however, we wish to predict an entire population of halos directly from a population of peaks. In light of the halo- and peak-count discrepancies, can our models proceed in this way?

V.1 Population comparisons

We first study the population of halos distributed in the asymptote $A$. In Fig. 14, we compare the entire halo population found in our simulations to the population predicted by accounting for every peak in the initial density fields that would have collapsed by simulation termination. We use the ellipsoidal collapse model. Generally, we see that for the narrower $w=0.1$, $w=0.3$, and EMDE power spectra the predicted population matches the simulated population reasonably well. In these scenarios, our model underpredicts halos at the low-density end, a surprising result that may be due to artificial fragmentation. In simulations with a small-scale power-spectrum cutoff (like ours; see Fig. 1), discreteness noise causes filamentary structures to fragment into halos even below the scales of the smallest density fluctuations. These halos have been shown to be unphysical simulation artifacts Angulo et al. (2013); Lovell et al. (2014), but they could contribute to the excess of less-dense halos in our simulations relative to model predictions in the cases of the narrower power spectra.

Otherwise, Fig. 14 shows a tendency for the model to overpredict halos at the middle-density range and underpredict at the highest-density range. This discrepancy can be attributed to mergers, for which our model does not account. Halo mergers reduce the number of halos while raising the central density of the merger remnants. For the broader $w=0.5$, 100 GeV, and 3.5 keV power spectra (bottom panel of Fig. 14), the merger-based discrepancy is amplified: the model underpredicts the densest halos and dramatically overpredicts the rest. Evidently, while our model can predict the density profiles of individual halos, an understanding of halo mergers is necessary to accurately predict halo populations in scenarios with more broadly supported power spectra. We will return to this point.

We next study the distribution in the radius $r_{\mathrm{max}}$. Many of the simulated halos had $r_{\mathrm{max}}>r_{\mathrm{vir}}$, and while we discarded those halos from previous analyses, doing so now would alter the populations. Instead, we account for this problem here by substituting $r_{\mathrm{vir}}$ for $r_{\mathrm{max}}$ in those cases. As Fig. 15 shows, the $s=0$ contraction model predicts the $r_{\mathrm{max}}$ distribution reasonably well for the narrow power spectra, with the bulk of the discrepancy arising from overprediction of the total halo count due to halo mergers. The predicted distributions are also more sharply peaked, an effect that may owe to the model’s neglect of spherical asymmetry (see Sec. IV.4). For the broader power spectra, the halo count discrepancy is magnified due to the much larger frequency of mergers. However, unlike in the case of the asymptote $A$, we will soon see that the radius $r_{\mathrm{max}}$ does not significantly increase due to mergers beyond what the models already predict. Hence, there is no significant underprediction of the largest-radius halos.

We can also ask how well our models predict aggregate observational signals. If we assume that all halos have $\rho\propto r^{-3/2}$ inner density profiles, then the total dark matter annihilation rate in halos is proportional to the sum $\sum A^{2}$ over all halos.242424There is a logarithmic sensitivity to $r_{\mathrm{max}}$, which we neglect. Also, note that there is a density cap imposed by annihilations Berezinsky et al. (1992), so the signal converges even for $\rho\propto r^{-3/2}$. In Table 4, we show how accurately the ellipsoidal collapse model predicts the aggregate annihilation signal. Remarkably, despite not accounting for halo mergers, we successfully predict the annihilation signal to within a factor of 1.3 for all simulations. The increased central density within merger remnants has compensated for the drop in halo count. One caveat is that in this calculation, we did not account for changes in the slope $\gamma$ of the $\rho\propto r^{-\gamma}$ inner profile resulting from mergers. This change must be accounted for separately and will significantly reduce the annihilation signal (see Paper II).

We also consider an aggregate microlensing signal. The apparent image of a background star is deflected by an angle proportional to $M_{\mathrm{2D}}(\xi)/\xi$, where $M_{\mathrm{2D}}(\xi)$ is the projected mass enclosed within impact parameter $\xi$. As an approximation, we claim that the deflection due to a halo is proportional to $M(r_{\mathrm{max}})/r_{\mathrm{max}}$ and a function related to the impact parameter. Integrating over impact parameters introduces a factor $r_{\mathrm{max}}^{2}$, the characteristic area of the halo. Hence, for a halo population, the aggregate lensing signal, considered as the expected deflection of a given image, is proportional to $\sum M(r_{\mathrm{max}})r_{\mathrm{max}}$ summed over all halos.

In Table 4, we show how accurately the $s=0$ contraction model predicts this signal. We find that our models predict the lensing signal significantly worse than the annihilation signal, although in four of the simulations, the prediction still agrees to about a factor of 2. For the simulations drawn from the 3.5 keV and 100 GeV power spectra, however, the model overpredicts the lensing signal by an order of magnitude. These power spectra yield the greatest prevalence of halo mergers, and unlike in the case of annihilations, mergers do not sufficiently boost the remnant’s lensing signal relative to its model prediction to compensate for the loss of halo count in mergers. Hence, the lensing signal in the simulations is much smaller than that predicted from the linear density field. We remark, however, that the story may change if instrument sensitivities are taken into account. Mergers predominantly destroy smaller halos, and their signals may have been beyond sensitivity limits regardless.

It is clear from these results that an accurate accounting of halo mergers is necessary to predict a halo population in any generality. However, for narrower power spectra such as our $w=0.1$, $w=0.3$, or EMDE cases, mergers are subdominant, and our models can predict the populations reasonably well. Moreover, when considering aggregate annihilation signals, halo mergers may not have a significant impact beyond altering the slope $\gamma$ of the $\rho\propto r^{-\gamma}$ inner density profiles. For these uses, we describe in Appendix C a method to sample the halo population directly from the matter power spectrum $\mathcal{P}(k)$, bypassing the step of drawing a density field $\delta(\bm{x})$. Halo populations predicted using this method are slightly different from those predicted using our density fields, but we propose that they are more accurate; the method we describe can easily sample a much larger number of peaks, and it is not subject to errors related to the finite size and finite grid spacing of a sampled density field.252525If sampling directly from the power spectrum is more accurate, one may wonder why we chose to compare our simulation results to those predicted from the less-accurate finite density fields (Figs. 14 and 15). The reason for this choice is that it leads to a more explicit test of our model. Since we compare a simulation’s results to model predictions using the same simulation’s initial density field, any discrepancy that is not a simulation artifact can be attributed to the model.

V.2 Halo mergers

A full treatment of mergers is beyond the scope of this work, but we are poised to make some observations. Our procedure for matching a halo to its predecessor density peak involved tracking each halo backward through time. During this process, we counted the number of major mergers this halo underwent, which we define to be a merger between two halos with mass ratio smaller than 3. Over half of our halos with well-resolved asymptotes experienced at least one major merger, and 12% experienced at least three.

In an effort to find a simple way to cull halos that end up merging, we explored cutting out density peaks that were too close to an earlier-collapsing density peak. We considered two characteristic comoving length scales below which to make these cuts: the scale $q_{\mathrm{pk}}$ associated with the density peak (see Sec. III) and the scale $q_{\mathrm{max}}$ where $\epsilon(q_{\mathrm{max}})=2/3$ (see Sec. IV.1). Unfortunately, neither of these cuts produced sensible results. The $q_{\mathrm{pk}}$ cut culled far too few halos, while the $q_{\mathrm{max}}$ cut culled far too many. Ultimately, we expect that a more sophisticated accounting of mergers will be necessary, possibly following along the lines of extended Press-Schechter theory Lacey and Cole (1993) or the Peak Patch algorithm Bond and Myers (1996a, b); Stein et al. (2019).

Another question is how halo density profiles change due to mergers. Reference Ogiya et al. (2016) found that successive mergers cause the slope $\gamma$ of the small-radius asymptote of the density profile, $\rho\propto r^{-\gamma}$, to become shallower than its initial value of $\gamma=3/2$. The same work also found that successive mergers increase the central density. Figure 16 shows that both of these results are borne out in our own simulations as well. There is a clear trend wherein more major mergers lead to successively denser262626While the coefficient $A$ is only well defined if $\gamma=3/2$, our procedure for measuring $A$ (see Appendix A) obtains a logarithmically averaged value of $\rho r^{3/2}$ over the inner density profile regardless. but shallower inner structures. On the other hand, we also see in Fig. 16 that major mergers do not significantly increase $r_{\mathrm{max}}$ beyond what our models already predict, although they do increase $M(r_{\mathrm{max}})$.

VI Conclusion

The first halos form by direct collapse of peaks in the primordial density field. If some of these halos survive the subsequent hierarchical clustering process, as evidence suggests Polisensky and Ricotti (2015); Angulo et al. (2017); Ogiya et al. (2016); Drakos et al. ; Berezinsky et al. (2008); Hayashi et al. (2003); Penarrubia et al. (2010); Ogiya et al. (2019), then they would be the densest dark matter objects in the Universe. In this work, we presented models that predict the density profiles of these halos directly from the properties of the density peaks that formed them, and we used high-resolution cosmological simulations in a variety of different scenarios to tune and validate these models. The models are described in Secs. III and IV and summarized in Table 3. They have simulation-tuned parameters, but these parameters are independent of cosmology and do not need to be retuned to accommodate different scenarios.

We treated small and large radii separately, with the understanding that these regimes form under radically different circumstances. At small radii, the first halos develop $\rho=Ar^{-3/2}$ density profiles Ishiyama et al. (2010); Anderhalden and Diemand (2013a); *anderhalden2013erratum; Ishiyama (2014); Polisensky and Ricotti (2015); Angulo et al. (2017); Ogiya and Hahn (2018); Delos et al. (2018a, b), and our model predicts the coefficient $A$ with remarkable success. Over the full range of cosmologies we explored, the scatter from our model is only in the vicinity of 30%-40%, and it is likely that a significant fraction of that scatter arises from numerical noise in our simulations. The density profiles at large radii are more varied, but we parametrize them using the radius $r_{\mathrm{max}}$ at which the circular velocity is maximized along with the mass $M(r_{\mathrm{max}})$ enclosed within that radius. For this regime, we employ models already present in the literature Gunn and Gott III (1972); Dalal et al. ; Salvador-Solé et al. (2012), but we supply calibrations that enable their use as predictive tools applicable to any cosmology. The scatter from our models is roughly 40%-50% in $r_{\mathrm{max}}$ and a factor of 2 in $M(r_{\mathrm{max}})$, and, similarly, not all of this scatter is physical.

In this way, the models we presented can predict the halo arising from a given density peak. Our goal, however, is to predict populations of halos. For a power spectrum of density fluctuations that is narrowly supported, such as the spectrum imprinted by an EMDE Erickcek and Sigurdson (2011); Barenboim and Rasero (2014); Fan et al. (2014); Erickcek (2015) or certain inflationary models Salopek et al. (1989); Starobinskij (1992); Ivanov et al. (1994); Randall et al. (1996); Starobinsky (1998); Martin et al. (2000); Chung et al. (2000); Barnaby and Huang (2009); Barnaby (2010); Bugaev and Klimai (2011), our models replicate the entire halo population reasonably well. Thus, our models can serve as a tool to predict the observational signals of such cosmologies. However, for more broadly supported power spectra, such as those arising from a scale-invariant initial spectrum, halo mergers dramatically alter the halo population, causing the populations that arise in our simulations to differ substantially from the populations our models predict. Interestingly, modulo changes in the slopes $\gamma$ of the $\rho\propto r^{-\gamma}$ inner density profiles, dark matter annihilation signals seem to be sufficiently close to additive in halo mergers that our models still predict the aggregate annihilation signal to within 30% in every cosmology we tested.

Nevertheless, it is not clear that this additivity in the annihilation signal should extend beyond the timescales spanned by our simulations. More broadly, our models predict the initial halo population, but a proper understanding of halo mergers is needed in order to robustly connect it to the population today. It is necessary to understand both how mergers are distributed across halos and time and how they impact halo density profiles. Methods exist that can predict the distribution of mergers, the most prominent of which is the extended Press-Schechter theory Lacey and Cole (1993). The Peak Patch algorithm Stein et al. (2019) represents a method that may be easier to adapt, among other advantages, if more computationally expensive to apply. Meanwhile, the larger challenge is to predict how halo mergers alter density profiles. References Ogiya et al. (2016); Angulo et al. (2017); Drakos et al. (for major mergers) and Hayashi et al. (2003); Penarrubia et al. (2010); Ogiya et al. (2019) (for minor mergers) represent steps toward this goal, but there is not, as yet, a sufficiently general model.

Moreover, there is room for improvement in our models themselves. In predicting the density profiles at large radii, we assume that each mass shell contributes density in a profile that is a single power law up to a maximum radius. In reality, shell profiles follow more complicated forms that are sensitive to the total density profile; see Fig. 10 and Refs Lithwick and Dalal (2011); Dalal et al. . Utilizing more accurate shell profiles would likely improve the model predictions for the density profile, especially in predicting its broader shape rather than only $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$. Also, our model for the density profile at large radii completely discounts any deviations from sphericity in the initial peak. Another avenue for improvement may be to incorporate ellipsoidal collapse. At small radii, the reason the $\rho\propto r^{-3/2}$ profile arises is not well understood, and its accuracy has only been confirmed down to the resolution limits of $N$-body simulations. A physical understanding must be developed of the mechanism by which this profile arises in order to confirm whether the profile truly extends to arbitrarily small radii.

The primordial density field and its power spectrum of density fluctuations comprise a valuable window into the early Universe and the nature of dark matter. Our work in this paper was carried out as part of an effort to use the observational signatures of dark matter halos to probe these fluctuations. Further research is still needed to understand the impact of mergers on halo populations before these signatures, or their nonobservation, can be employed to robustly constrain cosmology. Nevertheless, the models presented in this work, which predict the initial halo population, represent a step forward in our capacity to use this probe.

Acknowledgements.

The authors thank Dragan Huterer for bringing Ref. Dalal et al. to their attention. The simulations for this work were carried out on the KillDevil and Dogwood computing clusters at the University of North Carolina at Chapel Hill. M. S. D. and A. L. E. were partially supported by NSF Grant No. PHY-1752752. Several key figures in this work employ the cube-helix color scheme developed by Ref. Green (2011).

Appendix A Collecting halo data

In Sec. II, we carried out six simulations and catalogued all halos present at the final redshift of each. The models we present in Secs. III and IV require data on the density profiles of these halos, and in this section, we detail our methods for collecting these data.

We first need to obtain the density profile and enclosed mass profile of each halo at the final redshift of each simulation. To reduce random noise, we use the procedure described in Paper II, wherein the profiles are averaged over a time interval. The profiles are binned at successive factors of $1.1$ in radius, but to mitigate noise associated with the binning scheme, we use a cubic spline to smoothly interpolate them. These profiles are only valid down to the radius $r_{\mathrm{soft}}$ corresponding to the separation below which simulation forces become non-Newtonian; for Gadget-2, $r_{\mathrm{soft}}$ is 2.8 times the force softening length. At large radii, we cut off the density profile at the radius $r_{\mathrm{vir}}$ inside which the mean enclosed density is 200 times the background density. Figure 17 shows a random sample of halo density profiles from each simulation.

The model developed in Sec. III predicts the coefficient $A$ of the $\rho=Ar^{-3/2}$ asymptote of the density profile at small $r$. Thus, we wish to extract $A$ from each density profile. Starting at $r_{\mathrm{soft}}$, we find the cumulative (logarithmic) average of $\rho r^{3/2}$ across radial bins, and we set $A$ to be the maximum of this cumulative average. The averaging procedure is intended to minimize noise resulting from employing a narrow range of radii while at the same time minimizing the influence of any bend in the density profile at larger radii. However, we also tested the alternative procedure of simply finding the average of $\rho r^{3/2}$ within a factor of 3 in radius above $r_{\mathrm{soft}}$ and found it to yield similar results.

We also wish to accommodate deviations from ${\rho\propto r^{-3/2}}$ at small radii. Unlike those of Paper II, our simulations do not have the resolution to clearly resolve the power-law index of the small-radius asymptote, but we may still hope to see statistical correlations in our large sample. For each halo, we measure the slope $\gamma$ of the small-radius asymptote $\rho\propto r^{-\gamma}$ by considering the radius range from $r_{\mathrm{soft}}$ to $3r_{\mathrm{soft}}$ and fitting a line in log space to the density profile within this range.

To test the models presented in Sec. IV, we also need to compute the radius $r_{\mathrm{max}}$ at which the circular velocity is maximized along with the mass $M(r_{\mathrm{max}})$. Since we already obtained the enclosed mass profile $M(r)$ for each halo, we simply find the maximum value of $M(r)/r$. The parameters $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$ are then defined to be the radius and enclosed mass at which $M(r)/r$ is maximized.

Appendix B Collecting peak parameters

The models we present in Secs. III and IV employ data on the peaks in the linear density field ${\delta(\bm{x})\equiv\delta(\bm{x},a)/a}$ (evaluated during matter domination). We obtain these data using Fourier methods, defining

[TABLE]

The derivatives of $\delta$ at a peak located at $\bm{x}$ immediately follow as

[TABLE]

Similarly, using Poisson’s equation, the derivatives of the (peculiar) gravitational potential are

[TABLE]

The ellipsoidal refinement described in Sec. III.2 requires the three-dimensional shape parameters $e_{\phi}$ and $p_{\phi}$ for the potential $\phi$ about the peak. Taking $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$ to be the eigenvalues of $\partial_{i}\partial_{j}\phi$, these parameters are defined

[TABLE]

Our model in Sec. IV requires the density profile $\delta(r)$ and mass profile $\Delta(r)$ about the peak. For a peak centered at $\bm{x}$, these profiles are computed as

[TABLE]

where $\operatorname{sinc}(x)\equiv\sin(x)/x$ and $W$ is the top-hat window function, $W(x)\equiv(3/x^{3})(\sin x-x\cos x)$. The third profile $\zeta(r)$ is useful because it is related to the derivative of $\delta(r)$. In particular, to numerically integrate Eq. (16), we interpolate $\delta(r)$ and $\Delta(r)$ with piecewise polynomials, using the relations

[TABLE]

to fix their derivatives.

Figure 18 shows the density profiles $\delta(r)$ of three peaks from each of the six simulations. The peaks displayed are those that later collapse into the halos depicted in Fig. 17.

Appendix C Predicting the halo population from the power spectrum

A primary goal of this work is to enable a prediction of the halo population given a power spectrum $\mathcal{P}(k)$ of density fluctuations. One option is to sample a density field from the power spectrum, use the methods described in Appendix B to characterize the peaks, and then apply the models developed in this work. However, it is possible to exploit the statistics of a Gaussian random field to compute the halo distribution more directly. In this section, we outline a practical procedure to perform this computation by sampling from the peak distribution. Similarly to earlier sections, we define $\mathcal{P}(k)\equiv\mathcal{P}(k,a)/a^{2}$, where $\mathcal{P}(k,a)$ is the dimensionless matter power spectrum evaluated using linear theory during matter domination. All quantities derived therefrom, such as $\sigma_{j}$ and $\delta(\bm{x})$, inherit similar scaling.

C.1 Number density of peaks

The first step is to find the total number density $n$ of peaks. As derived in Ref. Bardeen et al. (1986), the differential number density of peaks in a Gaussian random field, in terms of parameters $\nu\equiv\delta/\sigma_{0}>0$ and $x\equiv-\nabla^{2}\delta/\sigma_{2}>0$, is

[TABLE]

where $\gamma\equiv\sigma_{1}^{2}/(\sigma_{0}\sigma_{2})$, $\sigma_{j}$ and $R_{*}$ are defined in Eqs. (2) and (3), and

[TABLE]

The $\nu$ integral can be carried out analytically, leading to

[TABLE]

and this equation can be integrated numerically over ${x\geq 0}$ to obtain $n$.

C.2 Asymptote $A$

Next, we use Monte Carlo methods to compute the distribution of coefficients $A$ of the $\rho=Ar^{-3/2}$ small-radius asymptote. The model described in Sec. III predicts this asymptote from the amplitude $\delta$ and curvature $|\nabla^{2}\delta|$ of the density peak along with the shape parameters272727Since these parameters describe the potential, they are different from the parameters $e$ and $p$ in Ref. Bardeen et al. (1986). $e$ and $p$ associated with the potential about the peak. Equation (26) supplies the peak distribution in $\nu\equiv\delta/\sigma_{0}$ and $x\equiv-\nabla^{2}\delta/\sigma_{2}$. Meanwhile, the conditional distribution of $e$ and $p$ for a peak of height $\nu$, derived in Ref. Sheth et al. (2001), is

[TABLE]

To compute the distribution of $A$, we now employ a Monte Carlo procedure. We use rejection methods to sample $x$ from Eq. (28) and then sample $\nu$ from Eq. (26). Next, we employ the cumulative distributions of $e$ and $p$,

[TABLE]

and

[TABLE]

numerically inverting them to inverse transform sample $e$ and $p$. This procedure yields a sample of peaks with parameters $\nu$, $x$, $e$, and $p$. Finally, we use Eqs. (6) and (7) to convert this sample into a halo sample distributed in the asymptote $A$, and we multiply the distribution by the total number density $n$ to obtain the differential number density $\mathrm{d}n/\mathrm{d}\ln A$.

As a test, we sample 400,000 density peaks from each of our six power spectra, and for each peak, we compute its predicted asymptote $A$ (without the proportionality constant) using the ellipsoidal collapse model of Sec. III.2. We plot the resulting distributions in Fig. 19 superposed with the distributions we find in the randomly generated simulation boxes. We find that the two distributions match well for all of our power spectra except for the 100 GeV spectrum. For this power spectrum, directly sampling the power spectrum yields halos of higher predicted density that we find in the box. This discrepancy is explained by noting that the simulation box only samples fluctuation modes up to the size of the box. The 100 GeV power spectrum has sufficient power at larger scales that neglecting it significantly reduces the amplitudes of fluctuations and therefore the density of the resulting halos. In this way, sampling peaks directly from the power spectrum is more accurate than using the intermediate step of sampling a density field.

C.3 Outer profile: $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$

So far, we have obtained the halo population distributed in the small-radius asymptote $A$. The next step is to extend this computation to find the multivariate distribution in $A$ and the outer profile parameters $r_{\mathrm{max}}$ and $M(r_{\mathrm{max}})$ using the models discussed in Sec. IV. This calculation is more difficult because, instead of the finite number of Gaussian variables relevant to the neighborhood of each peak, we must now handle a Gaussian distribution in the infinitely many variables corresponding to the density profile about the peak at each radius $q$. Nevertheless, a numerical computation is tractable.

The idea is to Monte Carlo sample the full density profile $\delta(q)$ about the peak. To accomplish this, we must discretize the radial coordinate such that

[TABLE]

for some large $N$. The maximum radius may be initially set at some $q_{N}\gtrsim R_{*}$ and raised as needed. For shorthand, we will write $\delta_{i}=\delta(q_{i})$ and use $\bm{\delta}$ to represent the full vector of $\delta_{i}$s. We now seek the distribution of $\bm{\delta}$ conditioned on the data we already used to find the asymptote $A$. In fact, $\bm{\delta}$ has vanishing covariance with the three-dimensional shape of the peak, so we need only the conditional distribution $f(\bm{\delta}|\nu,x)$. This distribution is Gaussian with mean

[TABLE]

(recalling $\gamma\equiv\sigma_{1}^{2}/[\sigma_{0}\sigma_{2}]$) and covariance matrix

[TABLE]

with the necessary covariances given by

[TABLE]

It is helpful to diagonalize $C=PDP^{T}$ so that $P$ is an orthogonal matrix and $D=\mathrm{diag}(\lambda_{1},...,\lambda_{N})$ is diagonal. If we define a new vector ${\bm{\kappa}\equiv P^{T}(\bm{\delta}-\bar{\bm{\delta}})}$, then $\bm{\kappa}$ is distributed as

[TABLE]

We then sample each $\kappa_{i}$ from its respective univariate Gaussian distribution, and the density profile $\bm{\delta}$ immediately follows using $\bm{\delta}=\bar{\bm{\delta}}+P\bm{\kappa}$. From the density profile $\bm{\delta}\rightarrow\delta(q)$, we can find the mass profile $\Delta(q)$ using Eq. (8). With these profiles in hand, it is now a simple matter to apply any of the models detailed in Sec. IV.

To test this procedure, we sample density profiles for the same 400,000 peaks from each power spectrum for which we already sampled the peak parameters. Using the $s=0$ contraction model in Sec. IV.2, Fig. 20 compares the $r_{\mathrm{max}}$ distributions binned from our initial density fields to those computed using the Monte Carlo method. As before, the distributions from the 100 GeV power spectrum are discrepant, likely owing again to the finite size of the random density fields. Also, there is a tendency for the random density fields to produce less sharply peaked distributions than the Monte Carlo method, especially for the $w=0.1$, $w=0.3$, and EMDE power spectra. It is unclear where this discrepancy arises, but it could be connected to the finite grid resolution of the random density fields. For a second test, we also plot the combined $A$-$r_{\mathrm{max}}$ distribution for the 3.5 keV power spectrum in Fig. 21. The Monte Carlo method proposed in this section accurately reproduces the correlations between $A$ and $r_{\mathrm{max}}$ found in the density field.