Relating the Structure of Dark Matter Halos to Their Assembly and Environment

The N-body simulation used here is the ELUCID simulation carried out by Wang et al. (2016) using L-GADGET code, a memory-optimized version of GADGET-2 (Springel 2005). The simulation uses 3072³ dark matter particles, each with a mass $3.08\times {10}^{8}\,{h}^{-1}{M}_{\odot }$, in a periodic cubic box of 500 comoving $\,{h}^{-1}\mathrm{Mpc}$ on a side. The cosmology parameters used are those based on WMAP5 (Dunkley et al. 2009), a ΛCDM universe with density parameters ${{\rm{\Omega }}}_{{\rm{K}},0}=0$, ${{\rm{\Omega }}}_{{\rm{M}},0}=0.258$, ${{\rm{\Omega }}}_{{\rm{B}},0}=0.044$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }},0}=0.742$, a Hubble constant ${H}_{0}\,=100\ h\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ with h = 0.72, and a Gaussian initial density field with power spectrum $P(k)\propto {k}^{n}$ with n = 0.96 and an amplitude specified by ${\sigma }_{8}=0.80$. A total of 100 snapshots, uniformly spaced in $\mathrm{log}(1+z)$ between z = 18.4 and z = 0, are taken and stored.

Halos and subhalos with more than 20 particles are identified by the friends-of-friends (FoF; see, e.g., Davis et al. 1985) and SUBFIND (Springel 2005) algorithms. Halos and subhalos among different snapshots are linked to build halo merger trees.⁸ Halo virial radius ${R}_{\mathrm{vir}}$ is related to halo mass, ${M}_{\mathrm{halo}}$, through

$$\begin{eqnarray}&&{M}_{\mathrm{halo}}=\displaystyle \frac{4\pi {R}_{\mathrm{vir}}^{3}}{3}{{\rm{\Delta }}}_{\mathrm{vir}}\overline{\rho },\end{eqnarray} \tag{ 1 }$$

where $\overline{\rho }$ is the mean density of the universe, and Δ_vir is an overdensity obtained from the spherical collapse model (Bryan & Norman 1998). The center of a halo is assumed to be the position of the most bound particle of the main subhalo. Halo mass ${M}_{\mathrm{halo}}$ is computed by summing over all particles enclosed within ${R}_{\mathrm{vir}}$. The virial velocity, ${V}_{\mathrm{vir}}$, is defined as ${V}_{\mathrm{vir}}=\sqrt{{{GM}}_{\mathrm{halo}}/{R}_{\mathrm{vir}}}$, where G is the gravitational constant. We use halos with masses $\geqslant {10}^{10}\,{h}^{-1}{M}_{\odot }$ directly from the simulation, and we use Monte Carlo–based merger trees to extend the mass resolution of trees down to ${10}^{9}\,{h}^{-1}{M}_{\odot }$ (see Chen et al. 2019).

From the halo merger tree catalog constructed above, we form four samples according to halo mass: S₁ contains 2000 halos with ${10}^{11}\leqslant {M}_{\mathrm{halo}}/\,{h}^{-1}{M}_{\odot }\leqslant {10}^{11.2};$ S₂ contains 1000 halos with ${10}^{12}\leqslant {M}_{\mathrm{halo}}/\,{h}^{-1}{M}_{\odot }\leqslant {10}^{12.2};$ S₃ contains 500 halos with ${10}^{13}\leqslant {M}_{\mathrm{halo}}/\,{h}^{-1}{M}_{\odot }\leqslant {10}^{13.2};$ and S₄ contains 500 halos with ${10}^{14}\leqslant {M}_{\mathrm{halo}}/\,{h}^{-1}{M}_{\odot }\leqslant {10}^{14.5}$. We also construct a complete sample, S_c, which contains 10,000 halos with ${M}_{\mathrm{halo}}/\,{h}^{-1}{M}_{\odot }\geqslant {10}^{11}$ to represent the total halo population. All halos in samples S₁, S₂, S₃, S₄, and S_c are randomly selected from simulated halos at z = 0. When halos need to be divided into subsamples according to some properties, the sample size of S_c may be insufficient. In this case, we construct a larger sample by the following steps. Starting from all simulated halos at z = 0 with ${M}_{\mathrm{halo}}/\,{h}^{-1}{M}_{\odot }={10}^{11}$ –${10}^{14.6}$, we divide them into mass bins of width 0.3 dex. If the number of halos in a bin exceeds 10,000, we randomly choose 10,000 from them; otherwise all halos in this mass bin are kept. This gives a sample of 94,524 halos, which is referred to as sample S_L. Due to the mass resolution of the ELUCID simulation, some properties of small halos cannot be derived reliably. Whenever these properties are needed, we use another halo sample, ${{\rm{S}}}_{{\rm{c}}}^{{\prime} }$, which contains all halos with ${M}_{\mathrm{halo}}\geqslant 5\,\times {10}^{11}\,{h}^{-1}{M}_{\odot }$ in sample S_c.

Note that halos with recent major mergers may have structural properties that are very different from virialized halos, and including them in our sample will significantly increase the variance of halo properties, thereby affecting the statistics derived from the sample. We use the criteria described in Appendix C to exclude those "unrelaxed" halos.

All of the samples used in this paper are summarized in Table 1. Note that we use only a fraction of all halos in a given mass range available in the simulation to save computational time. We have made tests using larger samples to confirm that the samples we use are sufficiently large to obtain robust results.

Table 1.  Samples Used in Our Analysis

Sample	Nhalo	${M}_{\mathrm{halo}}$/($\,{h}^{-1}{M}_{\odot }$)	Usage
S1	2000	10[11, 11.2]	Samples with constrained halo masses. Used in Section 3.1.
S2	1000	10[12, 12.2]
S3	500	10[13, 13.2]
S4	500	10[14, 14.5]
Sc	10000	≥1011	Mass-limited sample. Used in Section 3.1.
${{\rm{S}}}_{{\rm{c}}}^{{\prime} }$	2335	≥5 × 1011	Subsample of Sc, with all halo properties well defined. Used in Sections 3.2, 4.1, 4.2.
SL	94524	10[11, 14.6]	The "larger" sample for binning statistics. Used in Sections 4.3, 4.4.

Note. The four columns are sample identifier, number of halos, halo mass range, and usage, respectively. The exact sample definitions can be found in Section 2.1.

Download table as:  ASCIITypeset image

Following the literature (e.g., van den Bosch 2002; Wechsler et al. 2002; Zhao et al. 2003a, 2003b), we define the MAH of a halo as the main-branch mass ${M}_{z}$ as a function of redshift z in the halo merger tree rooted in that halo. Based on a theoretical consideration of halo formation (e.g., Zhao et al. 2009), we use the following quantity as the mass variable:

$$\begin{eqnarray}&&s(z)=\sigma ({M}_{0})/\sigma ({M}_{z}),\end{eqnarray} \tag{ 2 }$$

where $\sigma (M)$ is the rms of the z = 0 linear density field at the mass scale M. Similarly, we use

$$\begin{eqnarray}&&{\delta }_{{\rm{c}}}(z)={\delta }_{{\rm{c}},0}/D(z)\end{eqnarray} \tag{ 3 }$$

as the time variable, where ${\delta }_{{\rm{c}},0}=1.686$ is the critical overdensity for spherical collapse, and D(z) is the linear growth factor at z. We use the transfer function given by Eisenstein & Hu (1998) and the linear growth factor D(z) from Carroll et al. (1992). These definitions for the mass and time variables are well motivated by the self-similar behavior of halo formation expected in the Press–Schechter formalism (e.g., Press & Schechter 1974; Mo et al. 2010).

Thus, in our definition, the MAH of a halo is a vector ${\boldsymbol{s}}={(s({z}_{1}),s({z}_{2}),...,s({z}_{{\rm{M}}}))}^{{\rm{T}}}$, with each of its elements being the main-branch mass at a snapshot in the merger tree.⁹ Such a high-dimensional vector is obviously too complex to be useful in characterizing the formation of a halo. To overcome this problem, a common practice is to characterize the full MAH by a set of formation times (e.g., Li et al. 2008). In our analysis, we will use both the formation times and the principal components (PCs; see Section 3) to reduce the dimension of the MAH.

The MAH introduced above only uses the main branch of a halo merger tree, and thus it may potentially lose important information about the formation of a halo. However, our tests including the side branches (e.g., by modeling the assembly history with the progenitor mass distribution, as used in Parkinson et al. 2008) showed that it does not provide important information about halo structural properties. We therefore only use the main-branch assembly history for our analysis.

Dark matter halo profiles are usually modeled by a universal two-parameter form—the Navarro–Frenk–White (NFW) profile (Navarro et al. 1997),

$$\begin{eqnarray}&&{\rho }_{\mathrm{NFW}}(r)=\displaystyle \frac{\delta {\rho }_{\mathrm{crit}}}{(r/{r}_{{\rm{s}}}){\left(1+r/{r}_{{\rm{s}}}\right)}^{2}},\end{eqnarray} \tag{ 4 }$$

where ${\rho }_{\mathrm{crit}}=3{H}^{2}/(8\pi G)$ is the critical density of the universe. This profile is specified by a dimensionless amplitude, δ, and a scale radius, r_s. The scale radius is usually expressed in terms of the virial radius, ${R}_{\mathrm{vir}}$, through a concentration parameter, $c\equiv {R}_{\mathrm{vir}}/{r}_{{\rm{s}}}$. Since both $\overline{\rho }$ and Δ_vir in Equation (1) are known for a given cosmology, the profile of a halo can be specified by the parameter pair $({M}_{\mathrm{halo}},c)$, where the halo mass ${M}_{\mathrm{halo}}$ is defined in Section 2.1.

We obtain the concentration parameter of a halo through the following steps (see Bhattacharya et al. 2013):

• 1.  
  We divide the volume centered on a halo into N_r = 20 radial bins equally spaced between 0 and ${R}_{\mathrm{vir}}$ (see Section 2.1 for definitions of halo center and radius), and we calculate the mass within each bin i, M_i, using the number of particles in the bin.
• 2.  
  We compute the mass expected from the NFW profile in this bin, ${M}_{i,\mathrm{NFW}}(c)$, assuming a concentration parameter c.
• 3.  
  We define an objective function, ${\chi }^{2}(c)$, to be minimized as

  $$\begin{eqnarray}&&{\chi }^{2}(c)=\displaystyle \sum _{i=1}^{{N}_{{\rm{r}}}}\displaystyle \frac{{\left[{M}_{i}-{M}_{i,\mathrm{NFW}}\left(c\right)\right]}^{2}}{{M}_{i}^{2}/{n}_{i}}.\end{eqnarray} \tag{ 5 }$$
• 4.  
  We minimize the objective function to find the best-fit concentration parameter ${c}_{\mathrm{fit}}={\mathrm{argmin}}_{c}{\chi }^{2}(c)$.

In what follows, we will drop the subscript "fit" and use c to denote the concentration parameter obtained this way. As tested by Bhattacharya et al. (2013), changing the radius range and binning scheme in the fitting only introduces a negligible difference in c. According to our tests, our results presented in the following sections are also insensitive to such changes.

We model the mass distribution in a dark matter halo with an ellipsoid and use its axis ratio to characterize its shape. Our modeling consists of the following steps (see MacCiò et al. 2007):

• 1.  
  We start from a FoF halo and calculate the spatial position $d{{\boldsymbol{r}}}_{i}$ relative to the center of mass for each particle linked to the halo. The inertia tensor ${ \mathcal M }$ of the halo is given by the dyadic of $d{{\boldsymbol{r}}}_{i}$ summing over all of the ${N}_{{\rm{p}}}$ particles in that halo:

  $$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \sum _{i=1}^{{N}_{{\rm{p}}}}{m}_{i}d{{\boldsymbol{r}}}_{i}d{{\boldsymbol{r}}}_{i}^{{\rm{T}}},\end{eqnarray} \tag{ 6 }$$

  where m_i is the particle mass.
• 2.  
  We calculate the eigenvalues ${\lambda }_{{ \mathcal M },i}$ and eigenvectors ${v}_{{ \mathcal M },i}(i=1,2,3)$ of ${ \mathcal M }$ and rank the eigenvalues in descending order: ${\lambda }_{{ \mathcal M },1}\geqslant {\lambda }_{{ \mathcal M },2}\geqslant {\lambda }_{{ \mathcal M },3}$. So defined, ${v}_{{ \mathcal M },i}$ gives the axis direction of the inertia ellipsoid, and ${a}_{{ \mathcal M },i}$=$\sqrt{{\lambda }_{{ \mathcal M },i}}$ gives the length of the corresponding axis.
• 3.  
  We define the axis ratio q_axis of the halo as

  $$\begin{eqnarray}&&{q}_{\mathrm{axis}}=\displaystyle \frac{{a}_{{ \mathcal M },2}+{a}_{{ \mathcal M },3}}{2{a}_{{ \mathcal M },1}}.\end{eqnarray} \tag{ 7 }$$

  So defined, ${q}_{\mathrm{axis}}=1$ for a spherical halo, and close to zero if the halo is very elongated along the major axis.

Following Bullock et al. (2001), we define the spin parameter of a halo as

$$\begin{eqnarray}&&{\lambda }_{{\rm{s}}}\equiv \displaystyle \frac{\parallel {\boldsymbol{j}}\parallel }{\sqrt{2}{M}_{\mathrm{halo}}{R}_{\mathrm{vir}}{V}_{\mathrm{vir}}},\end{eqnarray} \tag{ 8 }$$

where ${M}_{\mathrm{halo}}$, ${R}_{\mathrm{vir}}$, and ${V}_{\mathrm{vir}}$ are, respectively, the halo mass, virial radius, and virial velocity defined in Section 2.1. The total angular momentum ${\boldsymbol{j}}$ is defined as

$$\begin{eqnarray}&&{\boldsymbol{j}}=\displaystyle \sum _{i=1}^{{N}_{{\rm{p}}}}{m}_{i}d{{\boldsymbol{r}}}_{i}\times d{{\boldsymbol{v}}}_{i},\end{eqnarray} \tag{ 9 }$$

where m_i is the particle mass, and $d{{\boldsymbol{r}}}_{i}$ and $d{{\boldsymbol{v}}}_{i}$ are particle position and velocity vectors relative to the center of mass, respectively. The summation is over all of the N_p particles linked in the FoF halo.

Many definitions can be found in the literature to characterize the environment of a halo at different scales, traced by different objects, and including or excluding the halo itself (see Haas et al. 2012 for a review). Here we define two quantities to describe the environment in which a halo resides: the density contrast as specified by the bias factor, and the tidal tensor. These quantities are computed directly from the N-body simulation.

First we follow Mo & White (1996) to define the halo bias b as the ratio between the halo–matter cross-correlation function and the matter–matter autocorrelation function. For each simulated halo i, we calculate its bias b_i by

$$\begin{eqnarray}&&{b}_{i}=\displaystyle \frac{{\xi }_{\mathrm{hm},i}(R)}{{\xi }_{\mathrm{mm}}(R)},\end{eqnarray} \tag{ 10 }$$

where ${\xi }_{\mathrm{hm},i}(R)$ is the overdensity centered at halo i at a radius R, and ${\xi }_{\mathrm{mm}}(R)$ is the matter–matter autocorrelation function at the same radius in the simulation at the redshift in question. On linear scales, the bias factor is expected to depend only on halo mass, independent of R. We have checked the values of b on different scales and found that the bias factor is almost constant at $R\gt 5\,{h}^{-1}\mathrm{Mpc}$. In the following, we compute both ${\xi }_{\mathrm{hm},i}$ and ${\xi }_{\mathrm{mm}}$ for R between 5 and $15\,{h}^{-1}\mathrm{Mpc}$ at z = 0, and we obtain the corresponding local linear bias b_i for each halo using the above equation.

For the tidal field, we follow Ramakrishnan et al. (2019) and define the halo environment in the following steps: we first divide the simulation box into a sufficiently fine grid (of ${N}_{\mathrm{cell}}^{3}$ grid points) and compute the density field $\rho ({\boldsymbol{x}})$ on each grid point using the clouds-in-cell method (Hockney & Eastwood 1988). To describe the environment of a given halo, the density field is smoothed on some scale R_sm with a Gaussian kernel. We then compute the potential field ${\rm{\Phi }}({\boldsymbol{x}})$ by solving the Poisson equation

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}({\boldsymbol{x}})=4\pi G\overline{\rho }\delta ({\boldsymbol{x}}),\end{eqnarray} \tag{ 11 }$$

where $\delta ({\boldsymbol{x}})=\rho ({\boldsymbol{x}})/\overline{\rho }-1$ is the overdensity and $\overline{\rho }$ the mean density of the universe. Next we obtain the tidal field ${ \mathcal T }({\boldsymbol{x}})$ through

$$\begin{eqnarray}&&{ \mathcal T }({\boldsymbol{x}})={\rm{\nabla }}{\rm{\nabla }}{\rm{\Phi }}({\boldsymbol{x}}).\end{eqnarray} \tag{ 12 }$$

Finally, we solve the eigenvalue problem of the tidal tensor at each grid point to find the eigenvalues ${\lambda }_{1}({\boldsymbol{x}})$, ${\lambda }_{2}({\boldsymbol{x}})$, ${\lambda }_{3}({\boldsymbol{x}})$ (ranked in descending order).

With all these, we obtain four environmental quantities for each halo: the local bias factor and the three eigenvalues of the tidal tensor: λ_i (i = 1, 2, 3). We follow the definition in Ramakrishnan et al. (2019) to define the local tidal anisotropy at each halo's position as

$$\begin{eqnarray}&&{\alpha }_{{ \mathcal T }}=\sqrt{{q}^{2}}/(1+\delta ),\end{eqnarray} \tag{ 13 }$$

where ${q}^{2}=\tfrac{1}{2}\left[{\left({\lambda }_{1}-{\lambda }_{2}\right)}^{2}+{\left({\lambda }_{2}-{\lambda }_{3}\right)}^{2}+{\left({\lambda }_{1}-{\lambda }_{3}\right)}^{2}\right]$ is the halocentric tidal shear, and $\delta ={\lambda }_{1}+{\lambda }_{2}+{\lambda }_{3}$ is the local overdensity. For our applications, we choose ${R}_{\mathrm{sm}}=4{R}_{\mathrm{vir}}$, although we have checked that our conclusion does not change significantly by using other values of R_sm. We use ${N}_{\mathrm{cell}}=2560$, which is sufficiently fine for computing ${\alpha }_{{ \mathcal T }}$ for the smallest halo in our sample.

In Table 2, we summarize the halo assembly, structure, and environmental properties used in our analysis.

Table 2.  Summary of Halo Properties Studied in This Paper

Property Type	Notation	Meaning
Assembly History	${\mathrm{PC}}_{\mathrm{MAH},i}$	ith PC of MAH
Structure	c	Concentration parameter of the NFW profile
	λs	Dimensionless spin parameter
	qaxis	Axis ratio of inertia momentum
Environment	b	Halo bias factor
	${\alpha }_{{ \mathcal T }}$	Tidal anisotropy

Note. Detailed descriptions of halo MAH PCs and structural and environmental quantities can be found in Sections 2 and 3.1.

Download table as:  ASCIITypeset image

For a given sample of N halos, the MAHs are represented by ${\boldsymbol{S}}={({{\boldsymbol{s}}}_{1},{{\boldsymbol{s}}}_{2},\ldots ,{{\boldsymbol{s}}}_{{\rm{N}}})}^{{\rm{T}}}$, where ${{\boldsymbol{s}}}_{i}$ is the MAH of the ith halo, as defined in Section 2.2. In simulation, the MAH of a halo is not traced below the mass resolution of the simulation. Halos of different mass therefore may be traced down to different snapshots, resulting in different lengths of the vector ${{\boldsymbol{s}}}_{i}$. For a given sample, we choose a snapshot where 90% of the MAHs can be traced back to this time. MAHs extending beyond this snapshot are truncated, and MAHs that are terminated before this snapshot are padded with 0. In this way, all of the ${{\boldsymbol{s}}}_{i}$ will have the same length, M, suitable for PCA. After applying the PCA to ${\boldsymbol{S}}$, the ${\boldsymbol{s}}$ for each halo is transformed into a new coordinate system, producing a new vector that consists of a series of PCs, ${{\bf{pc}}}_{\mathrm{MAH}}=({\mathrm{PC}}_{\mathrm{MAH},1},{\mathrm{PC}}_{\mathrm{MAH},2},\,\ldots ,\,{\mathrm{PC}}_{\mathrm{MAH},{\rm{M}}})$. In terms of the capability of capturing sample variance and reducing reconstruction errors, PCA is theoretically the optimal linear method. Thus, if the assembly histories of halos have some dominating modes, they are expected to be captured by the PCA.

The upper panel of Figure 1 shows the PVE and CPVE curves (see Appendix A.1 for definitions) for the PCs of the five samples, S₁, S₂, S₃, S₄, and S_c, defined in Section 2.1 (see also Table 1). Since the proportional explained variance (PVE) measures the fraction of sample variance explained by a PC, it is clear that the first several PCs, among all cases, can capture most of the variance in the halo MAH (>80% by using the first three PCs). This demonstrates that a strong degeneracy exists in the MAH of individual halos, suggesting that the assembly history can be described by using only a few eigenmodes. As shown by the CPVE curve, using a single parameter is insufficient to describe the MAH. It can at most explain as much variance as PC₁, which is 55% for the most massive halos ($\approx {10}^{14}\,{h}^{-1}{M}_{\odot }$) and 70% for the smallest halos ($\approx {10}^{11}\,{h}^{-1}{M}_{\odot }$). In principle, we can add more PCs, which typically leads to better capture of the subtle structures in the MAHs. The lower panel of Figure 1 shows the distribution of the error e in each sample when MAHs are reconstructed with the first three PCs (see Appendix A.1 for the reconstruction algorithm). The reconstruction errors are almost all below 10%, demonstrating that the MAHs of halos can be represented well by a small number of PCs.

Zoom In Zoom Out Reset image size

Figure 2 shows some examples of the reconstructed MAHs using the first three PCs. Here ${\rm{\Delta }}s(z)=s(z)-\overline{s}(z)$ is shown for several halos in each halo-mass-constrained sample (see Section 2.1 and Table 1), where $\overline{s}(z)$ is the mean MAH in the sample. The overall shape of the MAH is well captured by the reconstruction, although the MAHs of individual halos are quite diverse. Some fine structures in the MAH, caused by violent changes in the formation history due to merger events, are missed in the reconstruction. They can, in principle, be captured by including more PCs.

Zoom In Zoom Out Reset image size

The rapidly converged PVE, the sharply peaked distribution of the reconstruction, and the well-reconstructed MAHs of individual halos all indicate that PCs are effective in reducing the dimension of the halo MAH. In the following subsection, we will show the relation between PCs and some widely used halo formation times to gain more physical insights into different PCs.

The diversity of the assembly history shown in Figure 2 indicates that no single parameter can provide a complete description of the MAH. To reflect different aspects of the assembly history, different assembly indicators, for example, formation times, have been defined in the literature. These formation times are physically more intuitive compared with the more abstract PCs, although each of them only provides partial information about the MAH. Here we examine the relations between PCs and a number of halo formation times to gain some physical understanding of the PCs we obtain.

In Appendix B we list the formation times we use and their detailed definitions. Because of the large number of formation times and the potential nonlinear pattern in their relations with PCs, RF is an ideal tool for this task. In Appendix A.2 we describe the RF algorithm in detail. The two important outputs from the RF analysis are (1) the fraction of explained variance, R², and (2) the feature importance, ${ \mathcal I }(x)$, for any predictor variable x. These two quantities measure the performance of the regression and the contribution from each predictor variable in explaining the target variable, respectively. Figure 3 shows how different formation times contribute to the diversity of the halo MAH. Here we use sample ${{\rm{S}}}_{{\rm{c}}}^{{\prime} }$ in which all of the formation times are well defined, as described in Section 2.1 (see also Table 1), to regress the first three MAH PCs on all of the formation times. The performance, R², and the contribution ${ \mathcal I }(x)$ of each formation time x to each PC are plotted. For all PCs, the contribution from z_mb,0.04 is the most dominant, while z_mb,1/2 is also significant. However, the importance of both decreases in higher order PCs. Interestingly, the last major merger redshift, z_lmm, which contributes little to PC_MAH,1, is increasingly important as the PC order increases. For ${\mathrm{PC}}_{\mathrm{MAH},3}$, the contribution from z_lmm is comparable to those from the other two formation times, z_mb,0.04 and z_mb,1/2. Since mathematically higher order PCs are capable of capturing more subtle patterns in the feature space, this behavior of the importance curves means that assembly variables, such as z_mb,0.04 and z_mb,1/2, mainly describe the low-order, overall patterns of the halo assembly history, while major mergers are an important factor in producing the fine structure in MAHs.

Zoom In Zoom Out Reset image size

The contribution curves in Figure 3 automatically suppress variable competition where multiple degenerate feature variables compete in the prediction contribution to the target variable. In the RF, if two variables compete but one of them is slightly better, then the split algorithm prefers the better one and gives it a higher importance value ${ \mathcal I }(x)$. To show this more clearly, Figure 4 plots the correlation between z_mb,0.04, ${z}_{\mathrm{0.02,1}/2}$, and the first two MAH PCs. Strong and nearly linear correlations between ${\mathrm{PC}}_{\mathrm{MAH},1}$ and ${z}_{\mathrm{mb},0.04}$ and between ${\mathrm{PC}}_{\mathrm{MAH},1}$ and ${z}_{\mathrm{0.02,1}/2}$ are seen, indicating that the variances in both formation times contribute significantly to the variance in the MAH of the halos. Inspecting the contours, one can also see that the correlation between ${\mathrm{PC}}_{\mathrm{MAH},1}$ and ${z}_{\mathrm{mb},0.04}$ appears stronger. The larger contribution value to ${\mathrm{PC}}_{\mathrm{MAH},1}$ from ${z}_{\mathrm{mb},0.04}$ than from ${z}_{\mathrm{0.02,1}/2}$ shown in Figure 3 validates the strength of the correlation.

Zoom In Zoom Out Reset image size

Figure 4 also demonstrates that the description provided by a single halo formation time is incomplete. The large scatter seen in the contours of PC_MAH,1 versus formation times means that the direction of the largest scatter in the MAH space is not fully aligned with the scatter caused by formation time, and that other parameters must also contribute to the distribution of halos in the MAH space. Moreover, compared to ${\mathrm{PC}}_{\mathrm{MAH},1}$, ${\mathrm{PC}}_{\mathrm{MAH},2}$ has a much weaker correlation with the halo formation times. Thus, the information contained in the second PC of MAH, which contributes 10%–20% to the total variance as seen from the upper panel of Figure 1, is almost entirely missed when a single formation time is used to predict the MAH.

The degeneracy and completeness problems can, in principle, be overcome if we use PCs to describe the assembly history, as PCs are linearly independent of each other. Typically, the use of a small number of low-order PCs can solve most of the problems associated with the regression of halo MAHs. If more subtle information is needed, one can always add more PCs without introducing too much degeneracy into the problem.

It is well known that halo concentration is correlated with halo MAH (see e.g., Navarro et al. 1997; Jing 2000; Wechsler et al. 2002; MacCiò et al. 2008; Zhao et al. 2003a, 2003b, 2009). However, it is still unclear which single assembly parameter best predicts the concentration, and whether combinations of multiple parameters can improve the prediction precision. The same problem exists when we consider other halo structural properties, for example, the axis ratio q_axis and the spin parameter λ_s.

The difficulties involved here arise from the high dimension of the feature space, the degeneracy or correlation among predictors, a possible nonlinear effect from predictors to target, and the "bias–variance" trade-off in choosing model complexity. Again, Random Forest can be used to tackle these problems (see Appendix A.2). To this end, we build the regressor $y=\mathrm{RF}({\boldsymbol{x}})$, where y is one of the three structural properties: c, ${q}_{\mathrm{axis}}$ or ${\lambda }_{{\rm{s}}}$, and ${\boldsymbol{x}}$ are halo assembly indicators, either the first 10 MAH PCs, or all formation times. We also include halo mass in the predictor variables, because it is treated as one of the major parameters in many halo-related problems. All these regressors are built based on sample ${{\rm{S}}}_{{\rm{c}}}^{{\prime} }$ (see Section 2.1 and Table 1) in which all formation times are well defined for all halos.

The outputs of RF regressions, including the performance, R², and the contribution ${ \mathcal I }(x)$ from each predictor variable, are shown in Figure 5. As one can see from the red curves, when a large number of predictors are used, the upper limit in the prediction of c using assembly history is about 65%. This indicates that the concentration parameter c of a halo is largely determined by its assembly history. About 35% of the variance is still missing if one uses only the mean relation to predict c from assembly history. Furthermore, as seen from the left panel, the first MAH PC is by far the most important, accounting for about 67% of the total information provided by the MAH. As a comparison, in the right panel, the combination of the three formation times, z_mb,1/2, z_mp,1/2, and z_mb,0.04, contains about the same amount of information as PC_MAH,1. Thus, if a single parameter is to be adopted as the predictor of c, PC_MAH,1 is the preferred choice. We have also tried to combine halo formation times and PCs of the MAH as predictors, and we found that the overall performance changes little, indicating that the MAH PCs dominate the information content about halo concentration.

Zoom In Zoom Out Reset image size

For q_axis and λ_s, the upper-limit performances R², achieved by either PCs or formation times, are about 15% and 20%, respectively, much worse than c. No single variable seems to dominate the contribution, as indicated by the long and low tail in the ${ \mathcal I }(x)$ plot. These suggest that the axis ratio and spin can be affected by many factors related to halo mass assembly, but the effects are all small. The similarity in behavior between λ_s and q_axis suggests that these two quantities may share parts of their origins. Indeed, as we will show below (Section 4.2), λ_s and q_axis are strongly correlated, and both show a strong correlation with the anisotropy of the local tidal field.

Figure 6 shows the distribution of halos in the PC_MAH,1–structural parameter space. There is a strong trend that halos assembled late (large PC_MAH,1) tend to be less concentrated, and a weak trend that such halos tend to be more elongated and spin faster. These are all consistent with the output from the RF regressors and verify that the small contributions from assembly indicators to λ_s and q_axis are produced by the large scatter, rather than by variable competitions.

Zoom In Zoom Out Reset image size

We now add environmental predictors to the regression of the structural properties. To quantify the effects of variable competition, we adopt a commonly used approach called "growing," where a series of regressors is built up with an increasing number of predictors. (The approach is called "pruning" if the series runs reversely). Whenever there is a tight correlation between an added predictor and the predictors already used, competition will show up as changes in their importance values, ${ \mathcal I }(x)$, but the overall performance, R², will not be changed significantly by the addition.

As demonstrated in Section 4.1, among all halo assembly indicators, the first PC of the halo MAH is the dominating indicator for c. Even for q_axis and λ_s, it is still the most important although less dominating. We therefore start by building a Random Forest regressor $y={\mathrm{RF}}_{1}({\mathrm{PC}}_{\mathrm{MAH},1},{M}_{\mathrm{halo}})$, where y is one of the three structural quantities. Again, the inclusion of ${M}_{\mathrm{halo}}$ is motivated by the fact that halo mass is traditionally considered as one of the most important quantities distinguishing halos. We also calculate the performance ${R}_{1}^{2}$ as well as the contributions ${{ \mathcal I }}_{1}(x)$ of the regressor. We then add the two environmental quantities, the bias factor b and the tidal anisotropy ${\alpha }_{{ \mathcal T }}$, to build a second regressor, $y={\mathrm{RF}}_{2}({\mathrm{PC}}_{\mathrm{MAH},1},{M}_{\mathrm{halo}},b,{\alpha }_{{ \mathcal T }})$, and obtain the corresponding ${R}_{2}^{2}$ and ${{ \mathcal I }}_{2}(x)$.

The left panel of Figure 7 shows the contribution, ${ \mathcal I }(x)$, of these two regressors for each of the three structure properties described above, with the values of R² indicated, using sample ${{\rm{S}}}_{{\rm{c}}}^{{\prime} }$ (see Section 2.1 and Table 1). To estimate the uncertainty in the results, we generate 100 random subsamples, each consisting of half of the halos randomly selected from the original sample ${{\rm{S}}}_{{\rm{c}}}^{{\prime} }$ without replacement. The errors of R² and ${ \mathcal I }(x)$ are then estimated as the standard deviations among these subsamples. The reason we do not use the standard bootstrap technique is that the sampling with replacement will lead to an artificially large R² in the Random Forest regressor, because the repeated data points may appear in both the training set that shapes the decision trees and the out-of-bag (OOB) set (see Appendix A.2) that is used to test the performance, making the performance overestimated.

Zoom In Zoom Out Reset image size

In the case of c, the inclusion of environment only increases R² from 0.56 to 0.57, and the importance value ${ \mathcal I }({x}_{\mathrm{env}})$ is below 0.1. These two results indicate that environment has little impact on halo concentration c, that the concentration c is well determined by the first PC of the MAH, and that the weak dependence of c on environment is mainly through the dependence of halo assembly on environment (assembly bias). These are consistent with the finding of Lu et al. (2006) that the density profiles of individual halos can be modeled accurately from their MAHs. These are also consistent with the result obtained with the Gaussian process regression by Han et al. (2019), who found that the dependence of halo bias on halo concentration is mainly through the dependence of the bias on formation time.

For q_axis, the R² is doubled, from 0.05 to 0.11, when environment factors are included. From the importance curve, ${ \mathcal I }(x)$, one can see that these environment factors take away about half of the contribution from the halo mass and PC_MAH,1. These results together suggest that environment can affect halo shape significantly and is at least as important as halo mass and the PCs of the MAH.

For the spin parameter, the value of R² increases from 0.16 to 0.20 when environment factors are included. The increase, 0.04, is about 20%, suggesting that these environment factors do have a sizable effect on halo spin. The contribution, ${ \mathcal I }(b)+{ \mathcal I }({\alpha }_{{ \mathcal T }})=35 \%$, is larger than the 20% they contribute to R², suggesting that some of the contribution is actually taken from the halo assembly history and halo mass. Thus, when interpreting the dependence of halo spin on environment, one should remember that part of it may actually come from its degeneracy with halo assembly history.

Another difference between c and the other two structural parameters is in their values of R². Both q_axis and $\mathrm{log}{\lambda }_{{\rm{s}}}$ have fairly small R², much smaller than 50%. This indicates that the major contributors to these two structural parameters are not yet found. In general, factors that can affect halo structural properties can be classified into three categories: the initial conditions of halos, the intrinsic properties of halos (e.g., ${M}_{\mathrm{halo}}$, PC_MAH,1), and halo environment (e.g., ${\alpha }_{{ \mathcal T }}$ and b). The interaction of halos with their environment depends not only on the environment, but also on halo properties. For example, because a halo is coupled to the local tidal torque only through its quadrupole, we expect that the spin and shape of a halo are correlated. Motivated by this and the similar behavior of the shape and spin parameters revealed in Section 4.1, we add the spin parameter into the predictors of the shape, and vice versa. We denote the target structure parameter as y_α and the added structure parameter y_β. In addition, we also consider the "initial condition" of a halo by tracing all of the halo particles back to redshift z = 3 and using these particles to calculate the corresponding shape and spin parameters. This quantity for ${y}_{\alpha }$, denoted by ${y}_{\alpha ,z=3}$, is also added into the set of predictors. The right panel of Figure 7 shows the result when this set of variables are used to predict halo structures. Among all the predictor variables, y_β is the most important for y_α, indicating that q_axis and $\mathrm{log}{\lambda }_{{\rm{s}}}$ are strongly correlated. For λ_s, its initial value also matters, ranked as the second important predictor. The performance, R², for both λ_s and q_axis is now boosted significantly, to ∼0.3. However, even in this case, R² is still less than 50%, meaning that the causes of the main parts of the variances in both λ_s and q_axis are still to be identified. This also indicates that, unlike the concentration parameter c, which is determined largely by halo assembly, λ_s and q_axis may depend on the details of the initial conditions, assembly, and environment. Morinaga & Ishiyama (2020) provided a possible scenario where accretion from filaments may partly account for halo shape and orientation. However, the large scatter and weak trend in their results indicate that the driving factor of halo shape and orientation is still missing. All these suggest that many nuanced factors can contribute to the variances of λ_s and q_axis. Models for distributions of λ_s and q_axis have to take into account these nuances by assuming some random processes, such as a normal or log-normal process according to the central-limit theorem.

The regressors for halo structure (presented in Sections 4.1 and 4.2) are all built using the mass-limited sample. The model training processes and performance measurements are thus dominated by low-mass halos, which are more abundant. However, halos with different masses may have different properties. For example, the halo assembly bias, as reflected by the correlation between halo formation time and the bias factor, is found to be significant only for low-mass halos (e.g., Gao et al. 2005; Gao & White 2007; Li et al. 2008). It is, therefore, interesting to see how the structural properties of massive halos depend on assembly and environment, and how environmental and assembly effects on these halos are related to each other.

Here we quantify such a halo mass dependence by building RF regressors for subsamples of a given halo mass, using the large sample S_L (see Section 2.1 and Table 1). For halos with a given mass, ${M}_{\mathrm{halo}}$, and for each of the three structural properties, y = c, q_axis, and $\mathrm{log}{\lambda }_{{\rm{s}}}$, we build three forest regressors $y=\mathrm{RF}({\boldsymbol{x}}| {M}_{\mathrm{halo}})$ with different sets of predictor variables ${\boldsymbol{x}}$:

• 1.  
  ${\boldsymbol{x}}=({\mathrm{PC}}_{\mathrm{MAH},1},{\mathrm{PC}}_{\mathrm{MAH},2},{\mathrm{PC}}_{\mathrm{MAH},3})$, the first three PCs of the halo MAH;
• 2.  
  ${\boldsymbol{x}}=({\mathrm{PC}}_{\mathrm{MAH},1},{\alpha }_{{ \mathcal T }},b)$, the first PC of the halo MAH and environmental parameters;
• 3.  
  ${\boldsymbol{x}}=({\mathrm{PC}}_{\mathrm{MAH},1},{\mathrm{PC}}_{\mathrm{MAH},2},{\mathrm{PC}}_{\mathrm{MAH},3},{\alpha }_{{ \mathcal T }},b)$, the first three PCs of the MAH plus the environmental parameters.

The reason for including ${\mathrm{PC}}_{\mathrm{MAH},2}$ and ${\mathrm{PC}}_{\mathrm{MAH},3}$ is that the MAHs of massive halos may be more complicated than low-mass ones, and high-order PCs may be needed to capture the more subtle components in their MAHs.

Figure 8 shows the contribution curves and performances of regressors $y=\mathrm{RF}({\boldsymbol{x}}| {M}_{\mathrm{halo}})$ for different halo structural properties, y, using different predictor variables, ${\boldsymbol{x}}$, and for halos of different masses. In the case where ${\boldsymbol{x}}$ is the first three MAH PCs (panels in the left column), the PC_MAH,1 is always the most important for the structures of low-mass halos. However, as the halo mass increases, the importance of PC_MAH,1 decreases and eventually is taken over by higher order PCs (the second or the third). Massive halos may have more diverse accretion histories (see, e.g., Obreschkow et al. 2020, who found that tree entropy increases with halo mass), so their structural properties may also be more complex. By using higher order PCs, the complex formation history can be captured so that a better prediction for halo structure properties can be achieved.

Zoom In Zoom Out Reset image size

As discussed in Section 4.2, for the total halo population, environmental effects are different for the three halo structural properties. Similar conclusions can be reached for halos of a given mass. The middle and right columns in Figure 8 show the results for regressors that combine the MAH and environment as predictors. For the halo concentration, the PCs of MAH always outperform environmental quantities, although there is a slight increase in ${ \mathcal I }(x)$ for environment quantities at the high-mass end. Compared to the regressor with only MAH PCs (upper left panel), the performance R² including the two environmental properties (upper right panel) only increases slightly, indicating again that the environmental effect on halo concentration is mainly through the dependence of halo MAH on environment.

The environmental effect on the shape parameter, q_axis, is totally different. As seen from the contribution curves, the environment is as important as MAH, and including the environment variables increases R² significantly. This implies that the environmental effect on the halo shape parameter is important, and that the effect is not degenerate with that of the MAH. The environmental contribution to the spin parameter, λ_s, is intermediate, larger than that to the concentration parameter but smaller than that to the shape parameter. The value of R² after including environment variables increases, but less significantly than in the case of the shape parameter. This suggests that the environmental effect does contribute to halo spins, but part of the contribution is taken from the assembly.

As a final demonstration of the application of the Random Forest regressor, we show how assembly parameters correlate with environment for halos of a given mass. Such a correlation is usually referred to as the halo assembly bias. The purpose here is to identify the best correlated pair of variables $(X,Y)$ at a given halo mass, where X is an assembly property and Y is an environmental property. The method to measure the correlation strength is straightforward. First, we bin halos in sample S_L (see Section 2.1 and Table 1) into subsamples according to the halo mass. Within each subsample, we build a Random Forest regressor for each pair of variables $(X,Y)$. The value of ${R}^{2}(X,Y| {M}_{\mathrm{halo}})$ then provides a measurement of the correlation strength between the two quantities. Figure 9 shows the results for different cases, where the environmental quantity Y is either the bias factor b or the tidal anisotropy parameter ${\alpha }_{{ \mathcal T }}$, and X is either PC_MAH,1 or $\mathrm{log}{\delta }_{{\rm{c}}}({z}_{\mathrm{mb},1/2})$. As one can see, the dependence of $\mathrm{log}{\delta }_{{\rm{c}}}({z}_{\mathrm{mb},1/2})$ on b is present only for halos with ${M}_{\mathrm{halo}}\lt {10}^{13}{h}^{-1}{M}_{\odot }$ and totally absent for more massive ones. The values of R² between b and the two assembly properties are both smaller than 2%, indicating that the correlation between b and assembly history is weak. These results are consistent with those obtained previously (e.g., Gao et al. 2005; Gao & White 2007; Li et al. 2008): the assembly bias is significant only for low-mass halos, and one has to average over a large number of halos to detect the weak trend.

Zoom In Zoom Out Reset image size

For comparison, we also build regressors between structure properties (c, λ_s, and q_axis) and b for halos of a given mass. The results are shown in the left panel of Figure 9. Clearly, the dependence of these properties on b is also weak. The results are consistent with those obtained previously by Mao et al. (2018), who found that b is better correlated with c and λ_s than with assembly properties for massive halos.

In a recent paper, Ramakrishnan et al. (2019) showed that the tidal anisotropy parameter, ${\alpha }_{{ \mathcal T }}$, is a good variable that correlates well with many structural properties. We present the correlation between ${\alpha }_{{ \mathcal T }}$ and other halo quantities in the right panel of Figure 9. It is clear that ${\alpha }_{{ \mathcal T }}$ shows a better correlation with halo intrinsic properties than the bias factor, as indicated by the larger values of R². In particular, the correlations between the assembly properties (${\delta }_{{\rm{c}}}({z}_{\mathrm{mb},1/2})$ and ${\mathrm{PC}}_{\mathrm{MAH},1}$) and ${\alpha }_{{ \mathcal T }}$ are significant, except at the very massive end. In Section 4.2, we demonstrate that part of the contribution from the environment to the structural properties is produced by the degeneracy between the environment and the halo assembly history. The strong relation between ${\alpha }_{{ \mathcal T }}$ and assembly history is a proof of this degeneracy. Note that q_axis shows a strong correlation with ${\alpha }_{{ \mathcal T }}$ for massive halos, indicating that the local tidal field plays an important role in determining the shape of the halo.

As mentioned above, the small values of R² between assembly and environment imply that assembly bias is a weak relation compared to the variance. Detecting such a bias thus needs the average over a large number of halos. To obtain the mean trend in the bias relation, we can build linear regression models between pairs of variables and use the slopes of the regression, k, to represent the mean trend between the two variables in question. Figure 10 shows the slopes of halo properties with b and ${\alpha }_{{ \mathcal T }}$ for halos with different masses. Given any two variables, a larger absolute value of k means a more significant linear correlation of the two variables for halos of a given mass. The relations between halo environment and different halo properties show different trends with halo mass. For example, the slope $k(\mathrm{log}{\delta }_{{\rm{c}}}({z}_{\mathrm{mb},1/2}),b| {M}_{\mathrm{halo}})$ is significant only at the low-mass end, while $k({\mathrm{PC}}_{\mathrm{MAH},1},b| {M}_{\mathrm{halo}})$ is significant at the low-mass end and becomes less significant at the high-mass end, suggesting that halo assembly bias is more important for halos of lower mass. In contrast, the slopes $k(c,b| {M}_{\mathrm{halo}})$, $k({q}_{\mathrm{axis}},b| {M}_{\mathrm{halo}})$, and $k(\mathrm{log}{\lambda }_{{\rm{s}}},b| {M}_{\mathrm{halo}})$ are all more significant at the high-mass end. As discussed in Wang et al. (2011), this halo mass dependence is a result of the competition between the environmental effect and the self-gravity of the halo. For example, a higher density environment not only provides more material for halos to accrete, but also is the location of a stronger tidal field that tends to suppress halo accretion. Depending on the halo mass, one of the two effects dominates. For example, for lower mass halos where self-gravity is weaker, the local tidal field may be more effective in making them more elongated and spin faster, and in preventing them from accreting new mass, so as to make their mass assembly earlier and concentration higher. This is indeed what can be seen from the lower panels of Figure 10. We note, however, that the physical units of the different curves in these panels are not the same, so these curves cannot be compared directly with one another. Note also that a significant value of k does not necessarily imply a significant R² in Figure 9, and vice versa, because a linear model may not represent faithfully the data of nonlinear relations. In addition, the value of R² depends not only on the value of k, but also on the absolute value of the variance in the regressed sample.

Zoom In Zoom Out Reset image size

PCA is an unsupervised, reduced linear Gaussian dimension-reduction method (Pearson 1901; Hotelling 1933). Consider a set of N vectors ${\boldsymbol{X}}={({{\boldsymbol{x}}}_{1},...,{{\boldsymbol{x}}}_{{\rm{N}}})}^{{\rm{T}}}$, each in an M–d space, V_M. The idea of the PCA is to find an M'–d subspace, ${V}_{{\rm{M}}^{\prime} }$ (M' ≤ M), in which the projection of ${\boldsymbol{X}}$,

$$\begin{eqnarray}&&{\boldsymbol{X}}^{\prime} ={\boldsymbol{XP}},\end{eqnarray} \tag{ A1 }$$

has maximal variance, where ${\boldsymbol{P}}$ is the projection operator. It can be shown that the problem to be solved is equivalent to solving the eigenvalue problem for the sample covariance matrix of ${\boldsymbol{X}}$, defined as

$$\begin{eqnarray}&&{\boldsymbol{S}}=\displaystyle \sum _{i=1}^{{\rm{N}}}({{\boldsymbol{x}}}_{i}-\overline{{\boldsymbol{x}}}){\left({{\boldsymbol{x}}}_{i}-\overline{{\boldsymbol{x}}}\right)}^{{\rm{T}}},\end{eqnarray} \tag{ A2 }$$

where $\overline{{\boldsymbol{x}}}={\sum }_{i=1}^{{\rm{N}}}{{\boldsymbol{x}}}_{i}/N$ is the sample mean. If we rank the eigenvalues in descending order, the first eigenvector of ${\boldsymbol{S}}$, ${{\boldsymbol{v}}}_{1}$, is the direction along which the sample has the maximum projected variance, ${\lambda }_{1}={{\boldsymbol{v}}}_{1}^{{\rm{T}}}{{\boldsymbol{Sv}}}_{1}$. This variance is exactly the first eigenvalue of ${\boldsymbol{S}}$. Similarly, the ith eigenvector and the ith eigenvalue are, respectively, the direction and value of the ith largest variance. Consider the space ${V}_{{\rm{M}}^{\prime} }$ spanned by the first M' eigenvectors. The linearity of the transformation can be used to prove that the projected variance in M' is ${\sigma }^{2}={\sum }_{i=1}^{{\rm{M}}^{\prime} }{\lambda }_{i}$. Thus, one can project each data point ${\boldsymbol{x}}$ into ${V}_{{\rm{M}}^{\prime} }$ by ${\boldsymbol{P}}=({{\boldsymbol{v}}}_{1},...,{{\boldsymbol{v}}}_{{\rm{M}}^{\prime} })$, ${\boldsymbol{x}}^{\prime} ={{\boldsymbol{P}}}^{{\rm{T}}}{\boldsymbol{x}}$, to find a lower-dimension representation for it. The ith component of ${\boldsymbol{x}}^{\prime}$ is called the ith PC of this data point in the sample.

In general, any dimension-reduction algorithm will lose information contained in the original data. In PCA, the proportional variance explained (PVE) by the ith PC, defined as

$$\begin{eqnarray}&&{\mathrm{PVE}}_{i}=\displaystyle \frac{{\lambda }_{i}}{\mathrm{Tr}({\boldsymbol{S}})},\end{eqnarray} \tag{ A3 }$$

is used to quantify the importance of the ith PC. The cumulative PVE, defined as

$$\begin{eqnarray}&&{\mathrm{CPVE}}_{{\rm{M}}^{\prime} }=\displaystyle \sum _{i=1}^{M^{\prime} }{\mathrm{PVE}}_{i},\end{eqnarray} \tag{ A4 }$$

can be used to quantify the performance of using the first M' PCs in the dimension reduction. Typically, if the data in question are generated from an intrinsic process of lower dimension, the CPVE should quickly converge to 1 as M' increases. We will see that halo assembly histories have this property (Section 3.1).

The inverse operation of the projection, $\tilde{{\boldsymbol{x}}}={\boldsymbol{Px}}^{\prime}$, allows one to reconstruct the original vector, but with information loss. We use the following quantity,

$$\begin{eqnarray}&&e=| | \tilde{{\boldsymbol{x}}}-{\boldsymbol{x}}| | /| | {\boldsymbol{x}}| | ,\end{eqnarray} \tag{ A5 }$$

to quantify the reconstruction error of the data point ${\boldsymbol{x}}$.

The RF regressor or classifier is a supervised, decision-tree-based, highly nonlinear, nonparametric model ensemble method in statistical learning (Breiman 2001). Here we first introduce the decision tree algorithm, and then we discuss how the trees are combined into a forest to make a regression or classification.

Given a set of observations $D={\{({{\boldsymbol{x}}}_{i},{y}_{i})\}}_{i=1}^{N}$, each consisting of a vector of predictor variables ${{\boldsymbol{x}}}_{i}\in {{\mathbb{R}}}^{{\rm{M}}}$ and a target random variable y_i (continuous in the regression problem, discrete in the classification problem), a decision tree T can be trained to fit the data by minimizing some error functional ${ \mathcal E }({\rm{T}}| D)$. In regression problems, a common choice for the error functional, which we adopt here, is the mean residual sum-of-square,

$$\begin{eqnarray}&&{ \mathcal E }({\rm{T}}| { \mathcal D })=\displaystyle \frac{1}{N}\displaystyle \sum _{n=1}^{N}{\left[f({{\boldsymbol{x}}}_{i})-{y}_{i}\right]}^{2},\end{eqnarray} \tag{ A6 }$$

where $f({{\boldsymbol{x}}}_{i})$ is the predicted value for the ith observation by the tree T. The tree can then be used to predict the target value for a future test observation, ${\boldsymbol{x}}$.

A decision tree is built by sequentially bipartitioning the feature space along some axes. At each partitioned region ${ \mathcal R }$ in the feature space, the tree fits the observations by a constant function,

$$\begin{eqnarray}&&f({\boldsymbol{x}})=\displaystyle \frac{1}{{N}_{{ \mathcal R }}}\displaystyle \sum _{n=1}^{{N}_{{ \mathcal R }}}{y}_{i},\end{eqnarray} \tag{ A7 }$$

where the summation is over all observations in the region ${ \mathcal R }$, and ${N}_{{ \mathcal R }}$ is the number of training observations in this region. In the first step, the variable x to be bipartitioned and the position of the partition plane are both chosen to minimize the error functional ${ \mathcal E }({\rm{T}}| D)$. After a partition, the feature space is split into two subregions, each of which can be bipartitioned further to minimize the error functional. This recursive process partitions the feature space into a tree-like structure and is continued until some stop criterion (e.g., the maximum tree height, or the maximum number of tree nodes, or the minimum number of data points in individual leaf nodes, defined as the nodes at the top of the tree) is achieved.

Once a tree is built, the amount of error reduced by partitioning variable x can be computed as

$$\begin{eqnarray}&&{{ \mathcal I }}_{{\rm{T}}| D}(x)=C\displaystyle \sum _{n}{N}_{n}\left({{ \mathcal E }}_{n}-\displaystyle \frac{{N}_{{n}_{{\rm{L}}}}}{{N}_{n}}{{ \mathcal E }}_{{n}_{{\rm{L}}}}-\displaystyle \frac{{N}_{{n}_{{\rm{R}}}}}{{N}_{n}}{{ \mathcal E }}_{{n}_{{\rm{R}}}}\right).\end{eqnarray} \tag{ A8 }$$

Here the summation is over all tree nodes partitioned by variable x; ${{ \mathcal E }}_{n}$ is the error functional computed at observations in the region represented by node n before partition; ${{ \mathcal E }}_{{n}_{{\rm{L}}}}$ and ${{ \mathcal E }}_{{n}_{{\rm{R}}}}$ are the errors computed at the left and right child nodes of node n, respectively, after partition; and N_n, ${N}_{{n}_{{\rm{L}}}}$, and ${N}_{{n}_{{\rm{R}}}}$ are the number of observations in node n and its two child nodes, respectively. The normalization factor C is chosen so that the summation of ${ \mathcal I }$ from all feature variables is one. So defined, the quantity ${{ \mathcal I }}_{{\rm{T}}| D}(x)$ is the amount of contribution from variable x in building the regressor and can therefore be viewed as the important value of the variable in explaining target y.

Such a tree model suffers from the overfitting problem: the more complicated the tree is, the less training error it will have. However, the tree will eventually be dominated by noise as its height increases. Many methods have been proposed to control such overfitting, for example, the cross-validation, bootstrap, and jackknife ensembles. Here we adopt the RF, an extension of the bootstrap ensemble, designed specifically to deal with overfitting in tree-like algorithms. The building of an RF involves two levels of randomness. First, one uses n_re bootstrap resamplings, and each is used to train a tree ${T}_{i}\ (i=1,\,\ldots ,\,{n}_{\mathrm{re}})$. Second, when training each of the n_re trees, only a random subset (size ${n}_{\mathrm{var}}\lt M$) of all M predictor variables is used at each partition step. The n_re trees are then combined, and the final prediction for a given feature ${\boldsymbol{x}}$, denoted as $\mathrm{RF}({\boldsymbol{x}})$, is then averaged among all trees (arithmetic average in regression, and majority voting in classification). Also, the importance of the predictor, ${{ \mathcal I }}_{{\rm{T}}| D}(x)$, is averaged among trees, which we denote as ${{ \mathcal I }}_{\mathrm{RF}| D}(x)$. The overall performance of the RF is represented by the explained variance fraction, R², defined as

$$\begin{eqnarray}&&{R}^{2}=1-\displaystyle \frac{{\displaystyle \sum }_{n=1}^{{N}_{{\rm{T}}}}{\left[{y}_{n}-\mathrm{RF}({{\boldsymbol{x}}}_{n})\right]}^{2}}{{\displaystyle \sum }_{n=1}^{{N}_{{\rm{T}}}}{\left[{y}_{n}-\overline{y}\right]}^{2}},\end{eqnarray} \tag{ A9 }$$

where in principle the summation should be computed over an independent test sample, T, of size N_T. The RF regressor has the advantage that the test performance can be directly estimated with the OOB sample in the bootstrap process, and therefore an extra test sample is not necessary. In addition, RF does not suffer from the issue of scaling or arbitrary transformation of predictors, which exists in many nonlinear approaches, such as K-nearest-neighbors and support vector machines.

The RF method has some free parameters to be specified, and we choose them based on the following considerations: (1) The number of trees in the forest, n_re, should be as large as possible to suppress overfitting. But a larger n_re is computationally more difficult. In our analysis, we choose ${n}_{\mathrm{re}}=100$, which is sufficiently large for most applications of RF. (2) The number of predictors randomly chosen in the partition, n_var, also controls the suppression of overfitting. We optimize the value of n_var by maximizing the OOB score through grid searching. (3) The tree termination criterion affects the complexity of each tree. We choose to control the number of data points in the leaf nodes, s_leaf. This choice makes the tree self-adaptive when more data points are available, and also reduces issues associated with transformations of target variables and the choice of the error functional. The value of s_leaf is also optimized by maximizing the OOB score, again through grid searching.