Halo concentrations from extended Press-Schechter merger histories

TL;DR

This paper investigates the relationship between dark matter halo concentrations and their formation histories using merger trees from excursion set theory, identifying missing environmental effects and proposing modifications to improve scatter predictions.

Contribution

It demonstrates that incorporating environmental dependence into excursion set merger trees enhances the accuracy of halo concentration scatter predictions.

Findings

Model predicts median concentration-mass relation accurately.

Environmental dependence increases scatter in concentration predictions.

Modified merger trees show correct correlation between environment and formation epoch.

Abstract

We apply the model relating halo concentration to formation history proposed by Ludlow et al. to merger trees generated using an algorithm based on excursion set theory. We find that while the model correctly predicts the median relation between halo concentration and mass, it underpredicts the scatter in concentration at fixed mass. Since the same model applied to N-body merger trees predicts the correct scatter, we postulate that the missing scatter is due to the lack of any environmental dependence in merger trees derived from excursion set theory. We show that a simple modification to the merger tree construction algorithm, which makes merger rates dependent on environment, can increase the scatter by the required amount, and simultaneously provide a qualitatively correct correlation between environment and formation epoch in the excursion set merger trees.

Figures

Tables

Table 1. Table 1: Maximum posterior probability values and the region containing 68% of the posterior probability for all parameters used in our model (the nuisance parameters C 0 subscript 𝐶 0 C_{0} , C m subscript 𝐶 m C_{\mathrm{m}} and C a subscript 𝐶 a C_{\mathrm{a}} are not listed). For each parameter we also list the previously determined value where available.

	This work	Previous work
$a$	$+{0.791}_{-0.023}^{+0.074}$	$+{0.874}_{-0.005}^{+0.005}$	Benson (2017b)
$p$	$+{0.218}_{-0.130}^{+0.045}$	$-{0.031}_{-0.005}^{+0.005}$	Benson (2017b)
$A$	$+{0.302}_{-0.003}^{+0.013}$	$+{0.332}_{-0.0002}^{+0.0002}$	Benson (2017b)
$G_0$	$+{0.591}_{-0.010}^{+0.009}$	$+{0.635}_{-0.0002}^{+0.011}$	Benson (2017b)
${\gamma }_1$	$+{0.253}_{-0.019}^{+0.019}$	$+{0.176}_{-0.015}^{+0.002}$	Benson (2017b)
${\gamma }_2$	$+{0.124}_{-0.019}^{+0.017}$	$+{0.041}_{-0.009}^{+0.001}$	Benson (2017b)
$\alpha$	$+{0.077}_{-0.009}^{+0.017}$	—
$C$	$+{625.0}_{-15.0}^{+77.0}$	$+400.0$	Ludlow et al. (2016)
$f$	$+{0.061}_{-0.005}^{+0.020}$	$+0.020$	Ludlow et al. (2016)

Table 2. Table 2: Mean and scatter in log 10 ⁡ c 200 ​ c subscript 10 subscript 𝑐 200 c \log_{10}c_{\mathrm{200c}} for different concentration models, and as measured from COCO N-body halos, for halos in the mass range 9.402 ≤ log 10 ⁡ ( M 200 ​ c / M ⊙ ) < 9.902 9.402 subscript 10 subscript 𝑀 200 c subscript M direct-product 9.902 9.402\leq\log_{10}(M_{\mathrm{200c}}/\mathrm{M}_{\odot})<9.902 at z = 0 𝑧 0 z=0 .

Model	$⟨{\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}{𝒄}_{\mathrm{𝟐𝟎𝟎}𝐜}⟩$	${𝝈}_{{\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}{𝒄}_{\mathrm{𝟐𝟎𝟎}𝐜}}$
COCO (N-body)	1.084	0.152
Ludlow et al. (2016), $\alpha =+0.000$	1.104	0.103
Ludlow et al. (2016), $\alpha =+0.071$	1.077	0.114
Ludlow et al. (2014), $\alpha =+0.000$	1.077	0.133
Ludlow et al. (2014), $\alpha =+0.071$	1.051	0.140

Equations

$$r$$

$$\kappa (r)$$

$${\varphi }^2=\sum _{i=1}{\log }^2\left(\frac{N_i^{(\mathrm{n})}}{N_i^{(\mathrm{e})}}\right),$$

$${\log }_{10}{\sigma }_{{\log }_{10}{c}_{200c}}$$

$$a(c_{200c})$$

$$b$$

$${\log }_{10}{\sigma }_{{\log }_{10}{c}_{200c}}$$

$$a(c_{200c})$$

$$b$$

$$\frac{⟨{\rho }_{-2}⟩}{{\rho }_0}=C\frac{{\rho }_{\mathrm{crit}}(t_{\mathrm{c}})}{{\rho }_0},$$

$$n_{\mathrm{i}}=\frac{1}{M_{\mathrm{i}}}{\int }_{\sqrt{M_{\mathrm{i}-1}{M}_{\mathrm{i}}}}^{\sqrt{M_{\mathrm{i}}{M}_{\mathrm{i}+1}}}Mn(M,{\delta }_{\mathrm{e}})\mathrm{d}M.$$

$$w(t;{\delta }_{\mathrm{e}})={\delta }_{\mathrm{c}}(t)D^{\alpha {\delta }_{\mathrm{e}}-1}(t),$$

$$C=C_0{\left(\frac{M_1}{M_0}\right)}^{C_{\mathrm{m}}}{\left(\frac{a_1}{a_0}\right)}^{C_{\mathrm{a}}},$$

$$\log \mathcal{L}=-\frac{1}{2}\Delta {\text{C}}^{-1}{\Delta }^{\mathrm{T}},$$

Full text

††pubyear: 2018††pagerange: Halo concentrations from extended Press-Schechter merger histories–A

1 Introduction

It is now very well established that the density profiles of cold dark matter (CDM) halos can be well described by simple, universal forms such as the NFW profile Navarro et al. (1996, 1997) or, somewhat more accurately, by the Einasto profile (Gao et al., 2008; Ludlow & Angulo, 2017). The density profile is thought to play a significant role in shaping the properties of galaxies which form in dark matter halos, as it contributes significantly to the rotation curve, and gravitational potential.

In semi-analytic models of galaxy formation, the density profile of the halo directly affects determinations of galaxy sizes (Cole et al., 2000; Jiang et al., 2018), and in many models affects the mass loading of outflows from galaxies (Cole et al., 2000; Benson, 2012). While semi-analytic models of galaxy formation are most often applied to dark matter halo merger trees extracted from N-body simulations, from which dark matter profiles can be measured directly, they can also be applied to merger trees generated using other techniques, such as those based on excursion set theory (Bond et al., 1991; Lacey & Cole, 1993; Cole et al., 2000; Parkinson et al., 2008). Such approaches have some advantages over the use of N-body merger trees, for example allowing much finer time resolution to be achieved (which can affect the results of semi-analytic models; Benson et al. 2012), rapid exploration of different dark matter physics (Benson et al., 2013), and avoidance of numerical noise issues which occur in N-body halos at low particle number (Benson, 2017b).

However, merger trees built from excursion set theory have, so far, only provided halo masses—they say nothing about other key physical properties of halos such as their density profiles, spins, or environment. Typically, these quantities are assigned to halos in such trees by appealing to an empirical correlation with halo mass and redshift, or by drawing at random from some measured distribution. In the case of halo concentrations for example, the usual approach is to assign a concentration based on a concentration-mass-redshift relation measured from an N-body simulation. This approach has two significant disadvantages. First, concentration is not uniquely determined by halo mass and redshift—there is significant scatter in concentration at fixed mass and redshift. As such, this scatter will be missing from calculations which depend on halo structure (e.g. galaxy sizes; Jiang et al. 2018). This scatter cannot be incorporated by simply adding a random perturbation around the median concentration relation as the offset from the median is expected to be correlated across time in any given halo. The second disadvantage of this approach is that it is known that halo structure correlates with the formation history of the halo (Navarro et al., 1997; Bullock et al., 2001; Ludlow et al., 2013).

Based on these correlations, Ludlow et al. (2014) (see also Correa et al. 2015) and Ludlow et al. (2016) developed models which relate the density profiles of dark matter halos to their assembly histories. The model of Ludlow et al. (2016) is based upon the time evolution of the total mass of progenitor halos collapsed (the ‘collapsed mass history’, or CMH), and is able to correctly predict concentrations for both CDM and warm dark matter (WDM) power spectra, while the Ludlow et al. (2014) model, which is based on the ‘mass accretion history’ (MAH; the mass of the main progenitor halo) fails to reproduce the concentration mass relation in WDM. As shown by Ludlow et al. (2014) and Ludlow et al. (2016), these models work well for halos which are relaxed and in equilibrium, for which concentration is mainly set by formation time and is independent of other factors (such as substructure, recent collapse times, etc.). As such, the models are relevant to a biased set of simulated merger trees, since selecting relaxed halos naturally biases trees to those that form early (on average). While these models were developed and tested primarily on assembly histories extracted from N-body simulations, they were also shown to be applicable to simple, spherical collapse models. As such, they can be similarly applied to assembly histories derived from excursion set-based merger trees.

It has also been shown (e.g. Zehavi et al., 2018) that the formation times of halos correlate with environment. Any model which relates concentration to halo assembly history (and, therefore, formation time) should also consider the relation with environment. In excursion set theory, the large scale environment of a halo has no effect on its assembly history (Bond et al., 1991), at least in the standard case where the power spectrum is filtered using a window function that is sharp in $k$-space.

In this work we examine the outcome of applying the Ludlow et al. (2016) model for halo concentration to merger trees built using the Parkinson et al. (2008) algorithm. We further explore how a simple model for the effects of environment on halo assembly can be introduced into this algorithm, and show that this improves agreement with the distribution of concentrations derived from N-body simulations.

The remainder of this paper is organized as follows. In §2 we describe how we measure the concentrations of N-body halos, and estimate the uncertainties on these measurements. We then describe our algorithm for introducing an environmental dependence into the Parkinson et al. (2008) tree-building algorithm, and how the parameters of this model are constrained. In §3 we show the results of this model and compare it to those from N-body simulations. Finally, in §4 we discuss the implications of our results.

2 Methods

Our goal in this work is two-fold: to introduce a model for halo concentrations into the Parkinson et al. (2008) tree-building algorithm which is dependent upon the formation history of each halo, and to compare the results of that model to those from cosmological N-body simulations. This will require determining the full distribution of concentration parameter as a function of halo mass. To that end, we first describe how we measure concentrations of a sample of N-body halos and, importantly, how we determine the uncertainties in these measurements arising from the finite number of particles with which each halo is represented (Benson, 2017a, b).

2.1 Concentrations of N-body halos

2.1.1 Fitting procedure

We make use of the Copernicus Complexio (“COCO”) simulations of Hellwing et al. (2016) to determine the distribution of N-body halo concentrations. As described by Hellwing et al. (2016), halos were identified in the COCO simulations using the friends-of-friends algorithm (Davis et al., 1985), with a linking length parameter of $b=0.2$. For each halo a mass $M_{\mathrm{200c}}$ is determined as the mass within a sphere centred on the halo particle with the minimum gravitational potential, and enclosing a mean density equal to 200 times the critical density. As shown by Ludlow et al. (2012), halos of a given mass that collapse very recently (i.e. double their mass in less than a crossing time) have concentrations that are biased relative to the normal concentration-formation time relation. As such, we follow Ludlow et al. (2016) and exclude any halos which more than doubled their mass in the last 1.25 crossing times, or approximately the last 3.7 Gyr.

To measure concentrations of halos via their N-body representations we construct histograms of the number of particles in the halo in 31 spherical shells with logarithmically-spaced radii between $3.46\times 10^{-2}r_{\mathrm{200c}}$ and $0.764r_{\mathrm{200c}}$, where $r_{\mathrm{200c}}$ is the radius enclosing a mean density equal to 200 times the critical density, and then fit these using an Einasto profile (Einasto, 1965), which has been found to be a good fit to the density profile of cosmological N-body halos (Gao et al., 2008). In performing the fit we include only those bins which satisfy the inequalities:

[TABLE]

where $\epsilon$ is the softening length (230pc for the COCO simulations), and $\kappa(r)=t_{\mathrm{relax}}(r)/t_{\mathrm{circ}}(r_{\mathrm{200c}})$ is a convergence criterion defined by Power et al. (2003; their eqn. 20). We choose a value of $\kappa_{\mathrm{conv}}=1$ to ensure the profile is minimally affected by numerical relaxation. (Power et al. 2003 advocate $\kappa=0.6$, so we are more conservative in defining convergence as larger $\kappa$ implies better convergence, Navarro et al. 2010). We find that this relaxation criterion exlcudes only 1.2% of halos. We furthermore retain only those halos for which at least 16 bins satisfy the above inequalities to ensure that the fit is well constrained.

To find the best fit profile for each halo we minimize a goodness-of-fit measure:

[TABLE]

where $N_{i}^{\mathrm{(n)}}$ is the number of particles in the $i^{\mathrm{th}}$ bin of the N-body profile, and $N_{i}^{\mathrm{(e)}}$ if the mean number of particles expected in that bin assuming an Einasto profile111We also considered an alternative approach in which we maximized the Poisson likelihood $\log\mathcal{L}=\sum_{i}-N^{\mathrm{(e)}}_{i}+N^{\mathrm{(n)}}_{i}\log N^{\mathrm{(e)}}_{i}-\log\Gamma(N^{\mathrm{(n)}}_{i}+1)$, which may be expected to be valid if the number of particles in each bin obeys Poisson statistics (Benson, 2017a). We found that this leads to small but significant differences in the resulting distribution of concentration parameters. Specifically, the scatter in $\log_{10}c_{\mathrm{200c}}$ at fixed mass is decreased by around 0.02 dex. However, we find that the Poisson model is not a good description of the distribution of particle number in each bin, presumably because of fluctuations introduced by subhalos and other internal structure. This does, however, point to the need for a more careful understanding of this issue for precision measurements of halo concentrations..

Our Einasto profiles are described by three parameters: $m_{0}$ (the total mass of the halo in units of the nominal $M_{\mathrm{200c}}$ mass reported by subfind), $r_{-2}$ (the radius at which the logarithmic slope of the density profile equals $-2$, measured in units of the nominal $r_{\mathrm{200c}}$ radius reported by subfind), and $\alpha$ (the shape parameter of the Einasto profile). We choose to fix $\alpha=0.18$ (consistent with the typical shape of halos reported by Gao et al. 2008) as it is generally not well constrained by the data, and explore broad ranges of the remaining two parameters: $m_{0}=(0.5,1.5)$, $r_{\mathrm{s}}=(0.01,0.50)$. The range for $m_{0}$ is chosen to be sufficiently broad to encapsulate all plausible profiles, and that on $r_{\mathrm{s}}$ is chosen to encompass the plausible range of concentrations (Gao et al., 2008). Once the best-fit parameters for the Einasto profile, $(\hat{m}_{0},\hat{r}_{-2})$, have been determined we use them to compute the value of $\hat{r}_{\mathrm{200c}}$ (i.e. the radius enclosing a mean density equal to 200 times the critical density in the best-fit profile in units of the nominal $r_{\mathrm{200c}}$ radius reported by subfind, and which may differ slightly from unity) for this profile, and then compute a concentration $c_{\mathrm{200c}}=\hat{r}_{\mathrm{200c}}/\hat{r}_{-2}$.

2.1.2 Uncertainties in N-body concentration estimates

As has been explored by Trenti et al. (2010) and Benson (2017a, b), measurements of halo properties from N-body simulations are subject to noise arising from the discrete sampling of the underlying distribution function by a finite number of particles. To explore the effects of the discreteness noise on our concentration measures we perform two Monte Carlo experiments.

In the first, we generate spherically symmetric Einasto density profiles with a range of different concentrations, and sample from each profile using different numbers of particles via a Poisson process (Benson, 2017a). These N-body representations are then fit using the procedure described above (without any application of the Power et al. (2003) convergence criterion, which does not apply in this simple experiment). This is repeated for many different random realizations of the particle distribution, and for many different total numbers of particle in the halo to quantify the uncertainty and bias in the recovered concentration.

In our second experiment, we utilize N-body halos from the Millennium Simulation (Springel et al., 2005), fitting them using the same method as described above. However, we sample particles uniformly at random from each halo (with replacement) to generate a large number of random realizations of the halo, and repeat the fitting procedure on each realization. This sampling is initially done at a rate of $\mu=1$ (where we define $\mu$ as the mean number of times that any given particle will appear in a random realization), and then repeated with lower sampling rates to generate realizations of the halos with fewer particles. This allows us to assess the uncertainty and bias in the measured concentrations for realistic, cosmological halos.

The results of these experiments are shown in Figure 1. The scatter, $\sigma_{\log_{10}c_{\mathrm{200c}}}$, is shown as a function of the number of particles in a halo, $N_{\mathrm{h}}$, on the $x$-axis, and as a function of concentration (shown by colour—with $c_{\mathrm{200c}}=4$ shown in blue, increasing to $c_{\mathrm{200c}}=20$ shown in red). Solid lines show the results of fitting Monte Carlo realizations of idealized, spherical, Einasto profiles. Dashed lines indicate results from fitting cosmological N-body halos, with each dashed line corresponding to the mean scatter measured over all cosmological halos of that concentration. Cosmological halos which deviate significantly from the expected scaling are excluded—in these cases the halos are not well-described by Einasto profiles and so we do not expect our fits to give meaningful results.

We fit the dependence of $\sigma_{\log_{10}c_{\mathrm{200c}}}$ on particle number and concentration using simple polynomial fits. For idealized halos we find that the scatter in concentrations can be described by the model

[TABLE]

while for cosmological halos the scatter is described by

[TABLE]

In both cases the dependence of the relation on particle number is consistent with the expected $N_{\mathrm{h}}^{-1/2}$ scaling due to Poisson sampling. The normalization is significantly higher for the cosmological halos than for the idealized halos.

Achieving an uncertainty in concentration parameter less than 0.1 dex requires $N_{\mathrm{h}}\gtrsim 2\times 10^{3}$ for idealized halos, and $N_{\mathrm{h}}\gtrsim 6\times 10^{3}$ for cosmological halos. When fitting our model to match N-body data, we will use the model of equation (5) to forward model the scatter in concentration parameter.

We also considered the bias in concentration—defined as the difference between the measured and true concentrations. (In the case of cosmological halos, the mean concentration found when fitting to the halo particles sampled at a rate of $\mu=1$ is taken to be our estimate of the true concentration.) Idealized halos show a trend of positive bias in fits, while cosmological halos show no significant trend in bias—the bias estimates for individual halos scatter around zero. As such, we assume no bias when modeling concentrations.

2.1.3 Concentrations of COCO halos

We fit profiles from the COCO simulations (Hellwing et al., 2016) using the procedure described in §2.1.1. Results are shown in Figure 2. Points show a random subsample of all halos fitted. The solid yellow line shows the median concentration as a function of mass, while the solid green lines show the $16^{\mathrm{th}}$ and $84^{\mathrm{th}}$ percentiles of the distribution. The expected scatter in concentration arising from finite particle number effects (as computed using the model of equation 5) is negligibly small, less than $0.035$ dex, in this figure as all halos shown have $N_{\mathrm{p}}>25,000$ particles.

2.2 Merger trees

Our goal is to apply the concentration model of Ludlow et al. (2016) to merger trees generated using an excursion set-based approach. Specifically, we will use the algorithm of Parkinson et al. (2008)—which is itself a modification of the algorithm of Cole et al. (2000)—to construct merger trees. In Appendix A we describe a small optimization to the tree building algorithm of Cole et al. (2000).

2.2.1 Concentrations

To compute concentrations of our merger trees we apply either the Ludlow et al. (2014) or Ludlow et al. (2016) algorithms. In the Ludlow et al. (2016) algorithm the mean density within the scale radius, $\langle\rho_{-2}\rangle$, for a halo is given by

[TABLE]

where $t_{\mathrm{c}}$ is the time at which the halo had first assembled a mass $M_{-2}=4\pi\langle\rho_{-2}\rangle r_{-2}^{3}/3$ into progenitors of mass greater than $fM$, where $M$ is the mass of the halo in question. To find $r_{-2}$ for any given halo in a merger tree we use an iterative approach. We make an initial guess for $r_{-2}$ (typically based on the median concentration-mass-redshift relation) and use this, together with the known mass and virial radius of the halo, to determine $M_{-2}$. We then search backward in time through the progenitors of the halo until we find the time, $t_{\mathrm{c}}$, at which the halo first had a total mass equal to $M_{-2}$ in progenitor halos of mass $fM$ or greater. An estimate of $\langle\rho_{-2}\rangle$ can then be made from eqn. (6), which in turn can be used (along with the mass and virial radius of the halo) to compute an updated estimate of $r_{-2}$. This new value of $r_{-2}$ is used in the next iteration of this procedure, which is continued until the value of $r_{-2}$ has converged to a desired level. In this way, scale radii, $r_{-2}$, are found for all halos in the merger tree222For halos close to the resolution limit of our merger trees progenitors of mass $fM$ may not be resolved. In such cases we can simply assign a concentration from a measured concentration-mass-redshift relation. In this work we do not care about such halos as we ensure that halos for which we wish to compute the concentration are sufficiently well resolved that this problem does not occur..

For the algorithm of Ludlow et al. (2014) we use a similar approach, except that the condition that the mass in all progenitors above mass $fM$ must equal $M_{-2}$ at $t_{\mathrm{c}}$, we instead require that the mass in the primary (most massive) progenitor equals $M_{-2}$.

2.2.2 Environmental dependence

Since we wish to examine the effects of environment on our merger trees, we select for each tree an environmental overdensity (defined within some spherical region of radius $R_{\mathrm{e}}$) in the linear regime extrapolated to the present epoch, $\delta_{\mathrm{e}}$. This is drawn from the distribution of Mo & White (1996; their eqn. 9), that is from a normal distribution conditioned upon the fact that the region has not exceeded the threshold for collapse on any larger scale. This environmental overdensity is assumed to be the same for all progenitor halos in the given merger tree.

Each merger tree, with final halo of mass $M_{\mathrm{i}}$ (defined such that $M_{\mathrm{i}}>M_{\mathrm{i-1}}$), is considered to have an abundance

[TABLE]

That is, the abundance is chosen such that the number of such trees per unit volume contains the total mass expected from the halo mass function $n(M,\delta_{\mathrm{e}})$ within the mass interval associated with the tree.

To determine the environment-dependent halo mass function, $n(M,\delta_{\mathrm{e}})$, we use the peak-background split approach (Bond et al., 1991; Bond & Myers, 1996). This is formally derived within the excursion set theory, and predicts a mass function of the Press & Schechter (1974) form but with the mappings $\delta_{\mathrm{c}}\rightarrow\delta_{\mathrm{c}}-\delta_{\mathrm{e}}$ and $\sigma^{2}(M)\rightarrow\sigma^{2}(M)-\sigma^{2}(M_{\mathrm{e}})$, where $M_{\mathrm{e}}$ is the mass corresponding to the spherical region used to define the environment, and $\sigma^{2}(M)$ is the fractional variance in the mass in spheres containing mean mass $M$ in the linear theory density field extrapolated to the present day. In this work, we use the mass function of Sheth & Tormen (2002) with parameter values tuned as described in §2.2.3 rather than the Press-Schechter mass function. We retain this same mapping of $\delta_{\mathrm{c}}$ and $\sigma^{2}(M)$ suggested by the peak-background split.

In the case of the Press-Schechter mass function, when the peak-background split mass function is integrated over all environments (weighted by the appropriate distribution function—see Mo & White 1996) the result is the original Press-Schechter mass function (i.e. that unconditioned on environment). For halo masses which exceed the mass of the background, the unconditioned mass function applies. Therefore, in the case of the Press-Schechter mass function, the environmentally-averaged mass function is identical to the unconditioned mass function.

The same is not true for the Sheth-Tormen mass function. Instead, we find that averaging this mass function over environment leads to a result slightly larger than the corresponding unconditioned mass function. This overshoot peaks at around 5% close to the background mass, and approaches zero at low masses. This overshoot occurs because the distribution function we use for the background overdensity is derived for the constant barrier corresponding to the Press-Schechter mass function. In principle the correct distribution function could be found for the curved barrier corresponding to the Sheth-Tormen mass function, by numerical solution of the barrier crossing problem. Given the small magnitude of the overshoot we do not seek to find the correct background distribution here, but instead simply recalibrate the parameters of the Sheth-Tormen mass function such that when environmentally-averaged it agrees well with N-body measures of the mass function. Because of the overshoot, there is a discontinuity in the mass function at $M_{\mathrm{e}}$. To correct for this we multiply the mass function for $M>M_{\mathrm{e}}$ by a fixed factor to remove the discontinuity.

While the peak-background split can be used to accurately predict the environmental dependence of the halo mass function, it predicts no environmental dependence in halo formation histories. This is easy to see in excursion set theory, as the behaviour of the trajectory on scales larger than the collapse scale (at variances smaller than the collapse variance) is uncorrelated with the behaviour of the trajectory on smaller scales (larger variances)333This is true under the usual assumption of a sharp-$k$ filter—for other filters the trajectory is non-Markovian (Maggiore & Riotto, 2010)..

In the merger tree algorithm of Parkinson et al. (2008) the quantity $w(t)=\delta_{\mathrm{c}}(t)/D(t)$, where $\delta_{\mathrm{c}}$ is the critical linear theory overdensity for collapse at time $t$ and $D(t)$ is the linear theory growth factor, plays the role of a “time” variable, and sets the threshold which perturbations must reach to collapse and form halos. Therefore, to introduce an environmental dependence into this algorithm, we make an empirical modification such that

[TABLE]

where $\alpha$ is a parameter to be determined. This introduces a dependence on the mapping between merger tree “time”, $w(t)$, and true cosmic time, $t$, which depends on the environmental overdensity.

The choice of the radius used to define environment, $R_{\mathrm{e}}$, is somewhat arbitrary, but may influence the value of $\alpha$ in our model. We choose to use $R_{\mathrm{e}}=5h^{-1}$ Mpc (where $h=H_{0}/100\hbox{km~{}s}^{-1}\hbox{Mpc}^{-1}$) as this is larger than the collapse scale for the vast majority of halos, while still being a relatively “local” measure of environment. This point is discussed further in §3.2.

2.2.3 Calibration

Several predictions of our model—specifically the environment-averaged halo mass function, halo progenitor mass functions, and the distribution of halo concentrations as a function of halo mass—will be altered as a consequence of our introduction of environmental dependences into both the halo mass function, $n(M;\delta_{\mathrm{e}})$, by adopting the peak-background split model, and the halo collapse threshold, $w(t;\delta_{\mathrm{e}})$. Therefore, we recalibrate the parameters of these models to match measurements from N-body simulations. This calibration is carried out by running a Markov Chain Monte Carlo (MCMC) simulation. Our approach follows that of Benson (2017b) in detail, including utilizing the same MCMC algorithm and convergence criteria. Benson (2017b) constrained the parameters of the Sheth & Tormen (2002) mass function, and the Parkinson et al. (2008) merger tree algorithm to match the mass function and progenitor mass functions measured from the MultiDark Planck N-body (MDPL2) simulation (Klypin et al., 2016), after having removed splashback halos to avoid double-counting. We constrain these parameters in the same way here, except that we additionally average halo and progenitor mass functions over environment, and include the parameter $\alpha$ as an additional parameter in our MCMC simulation. Furthermore, when constructing progenitor mass functions we adopt the error distribution of Trenti et al. (2010) when convolving the intrinsic halo mass function with the expected error distribution. As the Trenti et al. (2010) errors are dominated by the effects of missing structure below the resolution limit, we expect that errors in halo masses will be correlated across time. This correlation should be measured directly from N-body simulations, using a methodology similar to that adopted by Trenti et al. (2010). Lacking such an analysis, we instead adopt a simple model to describe the covariance between the halo masses of parent and progenitor halos. Given a parent halo, “0”, and progenitor halo, “1” we assume that their masses are drawn from a distribution with covariance matrix S. The diagonal elements of this matrix are set to $\sigma^{2}_{0}$ and $\sigma^{2}_{1}$, where $\sigma_{0}$ and $\sigma_{1}$ are the root-variances in the parent and progenitor halo mass according to Trenti et al. (2010). For the off-diagonal elements we assume that the covariance is $C\sigma_{0}\sigma_{1}$, with

[TABLE]

where $a$ is expansion factor, $M$ is halo mass, and $C_{0}$, $C_{\mathrm{m}}$, and $C_{\mathrm{a}}$ are parameters to be determined. We include these three nuisance parameters in our MCMC simulation, adopting broad, uniform priors in the interval 0–1 for $C_{0}$ and 0–2 for both $C_{\mathrm{m}}$, and $C_{\mathrm{a}}$, and will marginalize over them in our final analysis.

The concentration model of Ludlow et al. (2016) has two free parameters, $f$ and $C$. We include the parameters $C$ and $f$ in our MCMC simulation, and constrain the model to match the distribution of concentrations found in the COCO simulations (see §2.1.1). Specifically, we construct histograms of the concentrations of halos in seven logarithmically-spaced mass bins spanning the range $9.4<\log_{10}(M/M_{\odot})<12.4$, and adopt a log-likelihood for each mass bin of

[TABLE]

where $\Delta$ is a vector of differences between the N-body and model concentration histograms, and C is a covariance matrix. The covariance matrix is taken to be the sum of that of the N-body halo histogram, and that of the model halo histogram, both assuming Poisson counting statistics. When constructing the histogram of model halo concentrations we exclude halos which have more than doubled their mass in the last 1.25 crossing times (3.7 Gyr; see §2.1.1; which we find excludes 2.6% of halos), and smooth by a Gaussian with width chosen to match the expected uncertainty in N-body halo concentration estimates (see eqn. 5)—this smoothing is taken into accounting when computing the covariance matrix.

For the parameters of the Sheth & Tormen (2002) halo mass function and Parkinson et al. (2008) merger tree algorithm we adopt the same priors as were used by Benson (2017b). For the parameter $\alpha$ we adopt a uniform prior in the range $0.0$ to $0.4$—we expect this parameter to be positive to induce early formation times in higher density environments (see §3), and we find that values larger than $0.4$ can lead to $w(t;\delta_{\mathrm{e}})$ becoming an increasing function of $t$, which would result in progenitor halos collapsing after their descendants and so is physically impossible. For the parameters $C$ and $f$ we adopt uniform priors of $(100,800)$ and $(0.01,0.10)$ respectively. These priors are broad and include the values found by Ludlow et al. (2016) to match results for both cosmological N-body halos, $(C,f)=(400,0.02)$, and simple spherical collapse models, $(C,f)=(650,0.02)$. Furthermore, Ludlow et al. (2016) find that $C$ and $f$ correlate with the slope of the $\langle\rho_{-2}\rangle$–$\rho_{\mathrm{crit}}(z_{-2})$ relation, when $f$ approaches the mass $M_{-2}/M_{200}$ (which is typically $0.1$–$0.2$). This correlation is currently neglected in our modeling, which further motivates our choice to restrict $f$ to values less than $0.1$. f restricted to small values, e.g. f=(0.01,0.1)

3 Results

The parameter constraints derived from our MCMC simulation are listed in Table 1, along with previously determined values. The maximum likelihood values of the Sheth & Tormen (2002) and Parkinson et al. (2008) model parameters are significantly shifted by the inclusion of environmental dependences, but the resulting halo and progenitor mass functions remain equally good matches to the MDPL2 N-body measurements. We find that the parameter $\alpha$ is strongly constrained to be non-zero, indicating that the N-body data are better fit by a model with environmental dependence in halo merger rates. For the parameters $(C,f)$ of the Ludlow et al. (2016) model we find best-fit values higher than those reported by Ludlow et al. (2016). This may reflect some difference in the structure of our merger trees compared to those extracted from N-body simulations, but we note that the posterior distribution over $(C,f)$ shows that these two parameters are strongly degenerate in the direction which includes the $(C,f)=(400,0.02)$ results favoured by Ludlow et al. (2016).

3.1 Concentrations

Figure 3 shows the distribution of concentrations, $c_{\mathrm{200c}}$, for halos in the mass range $9.402\leq\log_{10}(M_{\mathrm{200c}}/\mathrm{M}_{\odot})<9.902$ at $z=0$. The model of Ludlow et al. (2014) (right panel) produces a distribution which matches N-body results (shown by green points) quite well when no dependence on environmental density is included in $w(t;\delta_{7})$. However, as shown by Ludlow et al. (2016), the Ludlow et al. (2014) has significant failings, such as failing to reproduce the concentrations of halos in WDM models. The model of Ludlow et al. (2016) (left panel) predicts a distribution which is too narrow when $w(t)$ is independent of environmental density, but matches the N-body results reasonably well when $w(t;\delta_{7})=\delta_{\mathrm{c}}(t)D^{+0.071\delta_{7}-1}(t)$. Table 2 summarizes the mean and scatter in concentration in a narrow range of halo mass as measured from the COCO N-body halos, and predicted by the Ludlow et al. (2014) and Ludlow et al. (2016) models applied to our merger trees. Even with an environmental dependence introduced into $w(t;\delta_{7})$ the Ludlow et al. (2016) model is unable to produce a scatter quite as large as that measured for N-body halos. This is largely driven by the tail of low-concentration halos seen in the N-body simulation which is not matched by the Ludlow et al. (2016) model applied to our excursion set-derived merger trees.

Figure 4 shows the concentration–mass relation as measured from N-body halo profiles, and as predicted by the model of Ludlow et al. (2016).

3.2 Formation epochs

N-body halos show a correlation between environmental overdensity and formation time (defined as the time at which a given halo had first assembled 50% of its final mass into a single progenitor halo) as shown, for example, by Zehavi et al. (2018). While our model of environmental dependence was constructed to explain the scatter in halo concentrations at fixed mass, we can also ask whether it produces the correct correlation between environment and formation time.

To assess this correlation, and to compare to the measurements of Zehavi et al. (2018), we must account for the fact that our environmental overdensity, $\delta_{\mathrm{e}}$, is the linear theory overdensity in the Lagrangian volume surrounding each halo, while that measured by Zehavi et al. (2018) is the non-linear overdensity in an Eulerian volume surrounding each halo. While a fitting function relating these quantities has been proposed by Roth & Porciani (2011) we find that it does not well describe the relation here as we are interested in the relation for special points in the density field, namely the locations of halos. We therefore measure the linear theory Lagrangian, and non-linear Eulerian overdensities around particles in the Millennium N-body simulation (as was used by Zehavi et al. 2018).

Zehavi et al. (2018) define environment based on a $5h^{-1}$ Mpc Gaussian smoothing of the N-body particle field (computed in $2h^{-1}$ Mpc cells), which they label $\delta_{5}$. Based on their Figure 2, the median value of $1+\delta_{5}$ is around 2 for most halo masses. As such, these regions have collapsed in radius by a factor of approximately $2^{1/3}$ relative to their initial, Lagrangian volume. We therefore choose to use a radius $R_{\mathrm{e}}=2^{1/3}5h^{-1}\,\hbox{Mpc}=6.3h^{-1}\,\hbox{Mpc}$ to define environment in our model444We did not alter the value of $\alpha$ used in constructing these merger trees, even though it was constrained using $R_{\mathrm{e}}=5h^{-1}$ Mpc rather than $R_{\mathrm{e}}=6.3h^{-1}$ Mpc. Since the dispersion in $\delta_{\mathrm{e}}$ will depend on $R_{\mathrm{e}}$ it is possible that $\alpha$ should be recalibrated for each $R_{\mathrm{e}}$. Alternatively, it is possible that one specific value of $R_{\mathrm{e}}$ is the optimal choice (given some suitable metric to judge optimality) for our model and should always be used to define environment (in which case comparisons to measures using a different definition of environment become more difficult). We leave a full investigation of these issues to a future work..

Using the particle distribution in the Millennium simulation we measure the Lagrangian overdensity in spheres of radius $R_{\mathrm{e}}$ in the initial conditions around each halo’s most-bound particle, and extrapolate this to $z=0$ assuming linear perturbation theory. We then determine $\delta_{5}$ in the $z=0$ density field for the same particles. The results are shown in Fig. 5. We model the distribution of non-linear, Eulerian overdensity at each linear, Lagrangian overdensity as a log-normal distribution, and determine the mean and dispersion of this distribution from the measured points. Then, for each halo in our model, we take its assigned linear, Lagrangian overdensity, and draw a non-linear, Eulerian overdensity from the appropriate log-normal distribution.

Figure 6 shows joint and marginalized distributions of density contrast and formation epoch in three different halo mass ranges, and can be compared directly to Fig. 2 of Zehavi et al. (2018) in which these same distributions were measured directly from the Millennium Simulation (Springel et al., 2005). To construct these distributions from our model we generated large samples of merger trees as described in §2.2, matching the cosmological parameters and power spectrum of the Millennium Simulation. The resulting distributions of overdensity and formation epoch are in qualitative agreement with those found by Zehavi et al. (2018). In particular, the marginalized distributions of overdensity agree well with Zehavi et al. (2018), and are in good quantitative agreement for the position of the mode of the distribution, and the width. Marginalized distributions for formation epoch also agree well. Finally, the joint distribution shows a correlation between formation epoch and environmental overdensity which agrees moderately well with that found by Zehavi et al. (2018)—without our overdensity-dependent modification of the collapse threshold no correlation would exist. A more quantitative comparison of the effects of environment on assembly history in excursion set and N-body merger trees is currently limited by the different definitions of environment which naturally apply in these two approaches.

4 Discussion and Summary

We have applied the model of Ludlow et al. (2016) to merger trees generated via the algorithm of Parkinson et al. (2008) to predict the concentration of halos in those merger trees from their formation histories. We find that, with the original Parkinson et al. (2008) algorithm this model results in insufficient scatter in concentration at fixed halo mass. Since the Ludlow et al. (2016) model predicts the correct amount of scatter when applied to merger trees extracted from N-body simulations this suggests that the failure lies within the merger tree construction algorithm itself. We hypothesize that the missing ingredient is the effects of environment on halo formation times, and introduce a simple, empirical model for this effect to the Parkinson et al. (2008) algorithm. With this modification our model correctly predicts the mean concentration as a function of halo mass, and comes closer to matching the measured scatter, while simultaneously matching the statistics of progenitor halo mass functions.

We also show that this simple model for the influence of environment on halo formation reproduces, at least qualitatively, the correlation between environment and formation epoch measured in N-body simulations. This brings some of the important effects of assembly bias into the extended Press-Schechter framework for merger tree construction.

This model has important consequences for semi-analytic models of galaxy formation built on such merger trees. The concentrations of halos are a key factor in determining the sizes of galaxies (Jiang et al., 2018), and determine key observable properties such as rotation curves. To illustrate the effects of this work we utilize the Galacticus model (Benson, 2012) to predict the distribution of disk stellar masses, scale lengths, and rotation speeds555In Galacticus galaxy sizes and rotation speeds are solved for assuming each disk is in rotationally-supported equilibrium in the combined gravitational potential of the dark matter halo and the baryonic content of the galaxy itself. The sizes and rotation speeds will therefore depend on the halo concentration through its effect on the dark matter density profile and potential. Stellar masses of galactic disks will also be affected by concentration since they depend on the sizes and rotation speeds of those disks (via the dependence of the algorithms for star formation and feedback on those quantities). See Benson (2012) for further details. for central, disk-dominated galaxies occupying halos in the mass range 1–$3\times 10^{12}\mathrm{M}_{\odot}$ at $z=0$. The resulting distributions are shown in Figure 7, which also shows the mean and standard deviation of each distribution. The yellow line shows the result obtained when the concentration of halos does not include any scatter (i.e. it is set to the mean concentration as a function of halo mass and redshift, specifically using the fitting formula given by Ludlow et al. (2016)). The blue line shows the result when scatter in concentration is included using the model of Ludlow et al. (2016) with parameter values as calibrated in this work. Incorporating scatter into halo concentrations has little effect on the mean of the distributions but does add to the scatter by a small but non-negligible increase. This is most readily apparent in the case of rotation speeds, in which the scatter is increased from 0.06 to 0.08 dex. Semi-analytic models which assign concentrations based on the mean concentration-mass relation will therefore miss a significant contribution to the scatter in these observables—as far as we are aware all other semi-analytic models either utilize a mean (or median) concentration-mass-redshift relation to assign concentrations to halos, or else do not make use of concentrations in deriving galaxy properties. Models which attempt to include this scatter by randomly drawing concentrations from the measured distribution of concentrations will, by construction, incorporate the scatter in observable quantities, but will miss any correlation with other galaxy properties that may be influenced by formation history.

Acknowledgements

ADL acknowledges financial support from the Australian Research Council (project number FT160100250). We thank Nick Crook, whose initial investigation of applying the Ludlow et al. (2016) model to PCH merger trees motivated this work. SMC acknowledges the support of the Science and Technology Facilities Council [ST/P000541/1].

Appendix A Optimization of the Merger Tree Construction Algorithm

We use the merger tree construction algorithm of Cole et al. (2000; see also Parkinson et al. 2008). In that algorithm, small steps in “time” (actually in $w(t)\equiv\delta_{\mathrm{c}}(t)/D(t)$, where $\delta_{\mathrm{c}}$ is the critical threshold for spherical collapse and $D(t)$ is the linear growth factor, which serves the role of a time variable) are taken. These steps, $\delta w$, are required to be sufficiently small that the probability, $P=R\delta w$ (where $R$ is the rate of branching events per unit interval of $w$), of a branching event in any timestep is small, typically $P<\epsilon$ with $\epsilon=0.1$. Additionally, the timesteps must be small enough to ensure that subresolution accretion onto the halo is not too large during the timestep, and that the approximations of the Cole et al. (2000) merger rate are valid. In the limit of high branching rate (which occurs when the mass of the branch in question is much larger than the mass resolution), the step will be limited by branching.

We make a small modification to this algorithm. We remove the limit on the timestep due to branching rate, retaining the other limits on the timestep. Then, in each timestep, we find the interval to the next branching event by drawing a value at random from a negative exponential distribution with rate parameter $R$. If this interval exceeds the maximum allowed timestep, no branching occurs, and the timestep proceeds. If the interval is less than the maximum allowed timestep, branching occurs at that point. In the regime of high branching rates this approach allows for larger (by a factor $1/\epsilon$ on average) timesteps to be taken. Note that we do not have to concern ourselves in the subsequent timestep with the fact that no branching occurred in the previous timestep because of the memorylessness nature of the negative exponential distribution. That is, the distribution of branching intervals conditioned on the fact that no branching occurred in the previous timestep is just the same negative exponential distribution.

The increase in speed of the tree building will depend on both the algorithm implementation, and the tree resolution. In the case of the specific implementation of the Cole et al. (2000) algorithm, with a mass resolution equal to $5\times 10^{-6}$ times the final halo mass, this optimization results in speed-ups by a factor of 1.2.