Multi-stream radial structure of cold dark matter haloes from particle trajectories: deep inside splashback radius

Our analysis is based on two cosmological $N$-body simulations (LR and HR) performed in a periodic comoving box with a side length of $41\,h^{-1}\mathrm{Mpc}$, loaded with different numbers of particles as listed in Table 1. We assume a flat-geometry $\Lambda$CDM universe consistent with the recent result from the Planck satellite: $\Omega_{m}=0.3089,\Omega_{\Lambda}=0.6911,h=0.6774,n_{s}=0.9667,\sigma_{8}=0.8259$ (Planck Collaboration et al., 2016). The initial conditions are generated by adding displacements to particles arranged in a regular lattice based on second-order Lagrangian perturbation theory (Scoccimarro, 1998; Crocce et al., 2006), sourced by a Gaussian random field drawn from the linear matter power spectrum computed using the CLASS Boltzmann solver (Blas et al., 2011). The LR and HR simulations share the same random realization to allow for straightforward comparison. We then evolve the positions and velocities of the particles using a TreePM code GINKAKU (Nishimichi, Tanaka & Yoshikawa in prep.). This code is developed based on a public library, the Framework for Developing Particle Simulators (FDPS; Iwasawa et al., 2016; Namekata et al., 2018), which is aimed at efficient particle simulations in modern supercomputers. Working in a hybrid MPI-openmp parallel mode, the library allows an efficient domain decomposition as well as communication between processors. The short-range tree force is accelerated by the use of SIMD instructions implemented in the Phantom-GRAPE library (Nitadori et al., 2006; Tanikawa et al., 2012, 2013), while the details of the long-range PM force can be found in Yoshikawa & Fukushige (2005) (see also Ishiyama et al. 2009, 2012).

To track the trajectories of particles, we utilized 1001 snapshots of LR, which were uniformly sampled in redshift from $z=5$ to $z=0$. Since LR comprises approximately five times the number of snapshots used in S20, we anticipate that we can analyze particle orbits in greater detail, particularly those that undergo apocentre passages more than five times. In order to assess the convergence of the density profile considering the limited softening length, we rely on the $z=0$ snapshot of HR. Note that we have not retained particle snapshots at higher redshifts for the HR simulation due to constraints in disk storage capacity. Therefore, we investigate the evolution of the phase-space structure based on LR data, and then verify the accuracy of the density profile at $z=0$ through HR data.

We identified haloes at $z=0$ using a 6D phase-space temporal friends-of-friends halo finder Rockstar (Behroozi et al., 2013). As explained in more detail in Section 3.1, we tracked both the centre and the surrounding particles of the haloes identified at $z=0$ backward in time to establish the most massive main progenitor branch.

While Rockstar can identify haloes undergoing major mergers, the positions, velocities, and particle memberships of such haloes are often highly uncertain. Additionally, as we will discuss detail in Section 3.1, it is typically challenging to rigorously track the main progenitor branch of merger trees for these haloes, especially when they are in close proximity to other massive haloes. To circumvent these challenges, we concentrate on relaxed haloes that are less influenced by recent mergers, using two parameters: the spin parameter $\lambda$ and the offset parameter $X_{\mathrm{off}}$, both calculated by the Rockstar halo finder. Here, the parameter $X_{\mathrm{off}}$ represents the distance between the density peak location and the centre of mass of particles, normalized by the virial radius $R_{\mathrm{vir}}$. On the other hand, $\lambda$ denotes the amplitude of angular momentum divided by $\sqrt{2}R_{\mathrm{vir}}V_{\mathrm{vir}}M_{\mathrm{vir}}$, and is referred to as the Bullock spin parameter in Behroozi et al. (2013). Following Klypin et al. (2016), we consider haloes that meet either of the following conditions as undergoing major mergers and remove them from our halo catalogue:

Fig. 1 illustrates the distribution of $X_{\mathrm{off}}$ and $\lambda$, and the outcomes of applying the criteria in equation (1). The distribution closely resembles that shown in Fig. 4 of Klypin et al. (2016).

In addition to the haloes undergoing major mergers, we also exclude subhaloes. Subhaloes are typically surrounded by particles that belong to a distinct nearby host halo. This complicates the verification of whether the DM particles are bound by subhaloes or not, and it makes it challenging to accurately determine their density profiles. Thus, we exclude subhaloes from our halo catalog through the following procedure. First, we calculate the mass, denoted as $M_{\mathrm{vir,all}}$, which represents the total mass of particles within the virial radius $R_{\mathrm{vir}}$ as calculated by the Rockstar halo finder ($R_{\mathrm{vir,ROCK}}$). The Rockstar halo finder computes $R_{\mathrm{vir,ROCK}}$ and the virial mass $M_{\mathrm{vir,ROCK}}$ after removing unbound particles, ensuring that $M_{\mathrm{vir,all}}$ is always greater than $M_{\mathrm{vir,ROCK}}$. For host haloes, most of the particles within $R_{\mathrm{vir,ROCK}}$ are bound to the host halo, so $M_{\mathrm{vir,ROCK}}\approx M_{{\mathrm{vir,all}}}$. In contrast, for subhaloes, many particles around them are bound to the host haloes and not to the subhaloes themselves, resulting in $M_{{\mathrm{vir,all}}}\gg M_{\mathrm{vir,ROCK}}$. Here, there is no theoretical value that strictly distinguishes between host haloes and subhaloes. To address this, we identify haloes that satisfies the condition:

as subhaloes, and remove them from our halo catalogue. Fig. 2 displays the distribution of $M_{\mathrm{vir,all}}/M_{\mathrm{vir,ROCK}}$ and the outcomes of applying equation (2). We retain those plotted by the black dots below the horizontal dashed line after the two conditions (685 haloes). We can see that the two conditions are almost independent: the ratio of the haloes removed by the second condition does not depend on the first condition.

In addition to the criteria discussed above, we introduce two additional criteria to exclude two more haloes that are considered to be undergoing major mergers, as discussed in Section 3.1. Table 2 provides a summary of the resulting halo catalog. We also present the mass ranges from S to XL that we will use for stacking analysis after Section 4.

In summary, even after excluding subhaloes and the haloes that are undergoing major mergers, nearly $70\%$ of all haloes remain in our catalog, with equations (1) and (2) serving as the primary criteria for constructing our halo catalog.

For a particle with position $\boldsymbol{r}$ and peculiar velocity $\boldsymbol{v}$, the radial velocity $v_{r}$ is defined as:

We monitor the sign of this quantity over different snapshots to determine $p$. Therefore, to calculate $p$, we must determine the values of $\boldsymbol{x}_{\mathrm{h}}$ and $\boldsymbol{v}_{\mathrm{h}}$ for each halo in our catalog at each snapshot. While these values can be found by applying the halo finder to the snapshots at higher redshifts, they change discontinuously over time because different simulation particles are used to determine them in each snapshot. This leads to discontinuous variations in $v_{r}$ between snapshots, making it challenging to accurately count $p$.

To overcome this challenge, we calculate $\boldsymbol{x}_{\mathrm{h}}$ and $\boldsymbol{v}_{\mathrm{h}}$ as the average position and velocity of a fixed list of particles over time containing 1000 particles that can be considered the oldest progenitors. These progenitors are identified by tracking the massive branch of merger trees for each halo. By averaging the positions and velocities of these 1000 particles, we can calculate continuously varying $\boldsymbol{x}_{\mathrm{h}}$ and $\boldsymbol{v}_{\mathrm{h}}$ between snapshots while reducing the impact of individual particle noise.

The particles in the oldest progenitor are gravitationally well-bound within the halo at $z=0$ and move in tandem with it. Thus, these particles are suitable for representing the position and bulk motion of the halo centre. Here, the specific number 1000 is chosen to ensure numerical noise is minimised and can be adjusted as long as it provides a sufficient particle count. In Appendix A we vary this number and confirm that the results for $p\leq 40$ remain consistent quantitatively. Therefore, for the subsequent analysis, we focus solely on the particles with $p\leq 40$.

To construct merger trees, we initiate the process by determining the particles that compose each halo at $z=0$. This is achieved through a series of steps. First, we identify the centre of the halo using the shrinking sphere method. We then expand this sphere around the fixed centre until the overdensity within it decreases to the virial overdensity $\Delta_{\mathrm{vir}}$, as defined in Bryan & Norman (1998). The particles found within this last sphere are designated as members of the halo. In the shrinking-sphere method, we initially set the centre of the halo, as calculated by the Rockstar finder, as the starting point for the sphere. We then systematically reduce the radius of the sphere in a logarithmically equally manner through 120 bins, ranging from $R_{\mathrm{vir,ROCK}}$ down to the radius closest to the first centre (typically corresponds to $\mathcal{O}(10^{-3})\times R_{\mathrm{vir,ROCK}}$). The shrinking process concludes when the number of particles contained within the sphere falls below 100. We denote the centre-of-mass position of these 100 particles as $\boldsymbol{x}_{\mathrm{h,ss}}$. Note that the shrinking sphere method converges towards the primary peak of the halo. Consequently, this approach may not be suitable for haloes classified as subhaloes or undergoing major mergers, which exhibit another peak within $R_{\mathrm{vir}}$. This is another reason why we focus on relaxed and host haloes.

The results are shown in Fig. 3. Most of the haloes that passed the criteria of equations (1) and (2) exhibit close agreement between $\boldsymbol{x}_{\mathrm{h,ss}}$ and $\boldsymbol{x}_{\mathrm{h,ROCK}}$ to within $5\%$ of $R_{\mathrm{vir,ROCK}}$. However, there is one halo, despite satisfying the two previous criteria, with a mass of $\sim 5\times 10^{11}M_{\odot}$ and a relatively large positional difference ($\sim 20\%$ of $R_{\mathrm{vir}}$). In Appendix C, we carefully investigate this halo and find that it possesses a secondary peak, with both the peaks meeting the previously defined criteria. In general, haloes exhibiting significant deviations in $|\boldsymbol{x}_{\mathrm{h,ss}}-\boldsymbol{x}_{\mathrm{h,ROCK}}|$ can be considered as subhaloes or haloes undergoing major mergers. Therefore, we introduce an additional criterion based on Fig. 3 to exclude haloes that meet the following condition:

Next, we select the member particles of haloes in a snapshot one step before $z=0$ from those retained at $z=0$, employing the same method. This involves determining the centre of the halo using the shrinking-sphere method and expanding the sphere until the overdensity within it decreases to $\Delta_{\mathrm{vir}}$. These steps are repeated, tracking back to $z=5$, until the number of particles constituting the halo is below 1000. When the number goes below 1000, we slightly expand the sphere to ensure that it encompasses exactly 1000 particles, and we define the starting time of the counting of $p$, $t_{s}$, by the cosmic time of this snapshot. Note that when we move to the next snapshot (i.e., the one preceding it in cosmic time) in this procedure, we exclusively take into account the particles that are retained in the previous snapshot for the shrinking-sphere process. This guarantees that the list of member particles can only diminish as we move towards higher redshifts. Consequently, we define the particles within the oldest progenitor as consistently located near the halo centre from $t_{s}$ until the present.

As a result of this procedure, we can define the virial mass $M_{\mathrm{vir}}(z)$ and the virial radius $R_{\mathrm{vir}}(z)$ at each redshift as the mass contained within the sphere with an overdensity of $\Delta_{\mathrm{vir}}$ and its corresponding radius. These values are utilised in the calculation of the accretion rate in Section 5.1. For clarity, we distinguish these values from those computed by Rockstar at $z=0$ by denoting the latter simply as $R_{\mathrm{vir}}$ and $M_{\mathrm{vir}}$ (i.e., without the argument $z$) in what follows.

Fig. 4 illustrates the distribution of $t_{s}$, representing the cosmic time when the particle count in the progenitor of each halo falls below 1000. The majority of haloes exhibit $t_{s}$ values below $4\mathrm{Gyr}$, allowing us to trace the particles for over $\sim 10\mathrm{Gyr}$. However, note that in the case of 97 out of the 684 haloes in our catalogue, the tracking process extends all the way back to $z=5$, which is our first snapshot, without dropping below 1000 particles. While we can still define the centres of these 97 haloes using all the particles from their progenitor at $z=5$, the varying number of particles used to determine their centres may introduce systematic effects.

To address this issue, we choose to select 1000 particles from the progenitors at $z=5$, with a particular focus on their phase space. Initially, we calculate the position and velocity dispersions, denoted as $\sigma_{x}$ and $\sigma_{v}$, respectively, for the particles in the progenitor at $z=5$. Subsequently, we define the phase space distance, $d_{\mathrm{ps}}$, for each particle with respect to the average position $\boldsymbol{x}_{\mathrm{ave}}$ and velocity $\boldsymbol{v}_{\mathrm{ave}}$ for each particle using the following equation:

where $\boldsymbol{x}$ and $\boldsymbol{v}$ represent the positions and velocities of the particles, respectively. Finally, we select 1000 particles with the smallest values of $d_{\mathrm{ps}}$. This selection process ensures that we are focusing on particles that are well bound within the halo progenitor.

Fig. 5 illustrates the difference between the positions of halo centres at $z=0$, calculated by averaging over the positions of 1000 progenitor particles ($\boldsymbol{x}_{\mathrm{h,pro}}$), and those determined by the Rockstar halo finder ($\boldsymbol{x}_{\mathrm{h,ROCK}}$). Apart from one halo with a mass of $\sim 10^{12}M_{\odot}$, which exhibits a substantial difference ($\sim R_{\mathrm{vir}}$), we observe that all 683 $\boldsymbol{x}_{\mathrm{h,pro}}$ align with $\boldsymbol{x}_{\mathrm{h,ROCK}}$ within an accuracy of less than $10\%$ of $R_{\mathrm{vir}}$, with 564 $\boldsymbol{x}_{\mathrm{h,pro}}$ falling below $1\%$. This demonstrates that we can faithfully track the motion of the progenitor centres from $t_{s}$ to the present using the fixed list of 1000 particles. These particles are the oldest population that occupies the central region of a halo throughout its evolution. We utilise these particles to define the position and bulk motion of the centres, denoted as $\boldsymbol{x}_{\mathrm{h,pro}}$ and $\boldsymbol{v}_{\mathrm{h,pro}}$, respectively, at intermediate redshifts.

In Appendix D, we show that the exceptional halo we mentioned above possesses a significant secondery density peak that hinders us from accurately tracing the primary peak. This halo can be considered to be currently undergoing a major merger. Therefore, we introduce a final criterion below and remove haloes that satisfy it:

In summary, we constructed merger trees for each halo and selected 1000 particles in their progenitor; then we defined $\boldsymbol{x}_{\mathrm{h}}$ and $\boldsymbol{v}_{\mathrm{h}}$ as the average position and velocity of these 1000 particles at each snapshot. As demonstrated in Fig. 5, $\boldsymbol{x}_{\mathrm{h}}$ matches those calculated by Rockstar at $z=0$.

We can now calculate the radial velocity, $v_{r}$, of each particle at each of the 1001 snapshots using $\boldsymbol{x}_{\mathrm{h}}$, $\boldsymbol{v}_{\mathrm{h}}$ defined above, and count the number of apocentre passages from the change in sign of $v_{r}$. However, before proceeding with the actual counting of $p$, we need to account for particles residing in subhaloes. Our treatment of this issue is based on the method used in S20.

Generally, particles consisting a subhalo follow orbits relative to the its centre, and the sign of their radial velocities relative to the host halo can sometimes transition from negative to positive (and vice versa) due to this orbital motion within the subhalo. While we define $p$ as the number of apocentre passages relative to the host halo, this effect inadvertently increases the count of $p$. To mitigate this effect, we take into consideration the direction of each particle, denoted as $\hat{\boldsymbol{r}}\left(t\right)$, from the centre of the host halo as follows.

Suppose that a particle undergoes an apocentre passage at time $t_{1}$ and experiences a change in the sign of its radial velocity from positive to negative. Later, at time $t_{2}$, the positive-to-negative sign transition occurs for the first time after $t_{1}$. If $\hat{\boldsymbol{r}}\left(t_{1}\right)\cdot\hat{\boldsymbol{r}}\left(t\right)>0$ always holds for $t_{1}<t\leq t_{2}$, it implies that the particle remains on the same side of the halo, not completing a full orbit before the next apocentre passage. In such cases, the sign change at $t=t_{2}$ is likely due to orbital motion within the subhalo. We can show that a particle within a static spherical potential undergoes an apocentre passage after moving at least $180^{\circ}$ around the centre of that potential since the previous apocentre passage (see Sec.3.1 of Binney & Tremaine (2008)). Therefore, we count the number of apocentre passages only when $\hat{\boldsymbol{r}}\left(t_{1}\right)\cdot\hat{\boldsymbol{r}}\left(t\right)<0$ holds at least once for $t_{1}<t\leq t_{2}$. Note that this condition ensures that the particle has orbited at least $90^{\circ}$ from the previous apocentre passage. Nevertheless, we confirmed that this works in practice to remove fake apocentre passages due to the internal motion within a subhalo.

As an example, Fig. 6 shows the counting of $p$ for an orbiting particle using the method described above. In this figure, the colour of the points changes as the particle passes through its orbital apocentre, demonstraining the effectiveness of the method. Around $t=9\,\mathrm{Gyr}$, we can see that the sign of $v_{r}$ changes from positive to negative. However, the value of $p$ is not incremented because the particle has not traveled much from the previous apocentre passage (i.e., $\hat{\boldsymbol{r}}\left(t_{1}\right)\cdot\hat{\boldsymbol{r}}\left(t\right)$ remains positive since $t_{1}$). Indeed, the time since the previous apocentre passage is short compared to the typical intervals between other apocentre passages, suggesting that the sign change is due to internal motion with the subhalo. We can also see that the interval between the data points in the figure is sufficiently dense to count the apocentre passages. If the number of snapshots is insufficient, the algorithm cannot detect the apocentre passages, potentially leading to miscounts of $p$. To ensure accuracy, we conduct a verification in Appendix A, confirming that the $1001$ snapshots we used offer sufficient time resolution to avoid miscounting $p$ for $p\leq 40$. Furthermore, in Appendix B, we examine a convergence study and check if our simulation setup has enough mass resolution to fairly trace the multi-stream flow for a large number of apocenter passages, again confirming that the present mass resolution provides a reliable estimate of $p$ for $p\leq 40$.

Before presenting our main results, we provide a summary of some technical aspects related to calculating density profiles, stacking procedure, and fitting methodology. In Section 4.1, we showcase density profiles in physical coordinates. To construct these density profiles, we employ logarithmically spaced radial bins ranging from $2$ times the softening length to $2.5\,R_{\mathrm{vir}}$. Then, in Section 4.2, we utilise stacked density profiles, which involves averaging over individual profiles within specific subgroups. We define subgroups according to, for example, the number of particles with a certain value of $p$ and the mass. In our stacking procedure, we first normalise the radial distances of particles within individual haloes by their respective values of $R_{\mathrm{vir}}$, and create individual density profiles. Subsequently, we construct stacked density profiles by averaging over these individual profiles in the subgroup. For the stacked profiles for each mass bin, we define 80 logarithmically spaced radial bins spanning from $0.0025\,R_{\mathrm{vir}}$ to $2.5\,R_{\mathrm{vir}}$. However, we exclude bins that fall below $1.2$ times the maximum value of the softening length, scaled by $R_{\mathrm{vir}}$, for halos within the subgroup to avoid numerical artifact due to the limited force resolution. By comparing the LR and HR simulations, we confirm that the softening effect does not significantly impact the stacked profiles beyond this radius. These specific radii will be explicitly referenced in Section 4.2.2.

In Section 4, we will also present fitted curves for the density profiles. Our fitting procedure rely on the standard minimisation method for $\chi^{2}$, which is defined as follows:

Here, $\rho_{\mathrm{i,fit}}$ represents the fitted density in the $i$-th radial bin, while $\rho_{\mathrm{i,data}}$ denotes the individual or stacked density profile. $\sigma_{\mathrm{i}}$ is the uncertainty of $\rho_{\mathrm{i,data}}$, determined by the Poisson scatter for individual profiles and the root mean square deviation of the profiles stacked within subgroups. For our analysis, we consider all radial bins above the aforementioned lowest radial bins for each case.

Applying the method and algorithm in Section 3 to 683 halos in our simulations, summarised in Table 2, we are now able to characterise the radial phase-space structure by looking at the particle distributions classified with different numbers of apocentre passages. To do so, we first pick up three representative halos whose masses are $1.49\times 10^{14}$, $1.36\times 10^{13}$, and $1.45\times 10^{12}$ $h^{-1}M_{\odot}$, and plot their radial phase-space and density profiles in Figs. 7-9. In these figures, we show density profiles (top row) and phase-space distributions (bottom row) for particles with $p=1$ to $40$ (top left panel) and for $p=1,2,3,5,10,20,30,40$ (other panels).

Despite a wide range of mass scales, overall trends in the density profiles and phase-space structures are mostly similar for the three cases. Looking in particular at the outer part of the radial phase-space distribution, we see a substantial number of clumps with various spatial and velocity extents at $p=1$ and $2$. Although unrelaxed haloes are removed in our samples, this implies that the outer parts of haloes are generally not truly relaxed, having a large fluctuation in phase-space density due to the remnant of subhaloes that have recently accreted or undergone mergers. Note that this scatter has been reported in Diemer (2022) as the halo-to-halo scatter in characterizing the infalling profiles.

On the other hand, focusing on the inner structures (i.e., large values of $p$), the phase-space distribution becomes rather smooth, and we see fewer substructures. Furthermore, the density profile for a large value of $p$ generically exhibits double-power law features consisting of a shallow inner slope and a steep outer slope. At the transition scale where the slope suddenly changes, we see that particles in phase space are accumulated near zero velocity, associated with an apocentre passage for particles having the same $p$. In other words, this transition scale roughly corresponds to the radial caustic of each multi-stream flow. S20 reported that about $\sim 30$% of haloes have structures quantitatively similar to those predicted by the self-similar solution of Fillmore & Goldreich (1984). Although their self-similar solutions generally predict a spiky density structure at the caustics, such a spike is smeared in reality in the presence of non-zero tangential velocities and non-sphericity, leading to what is seen in Figs. 7 to 9, i.e., a smooth transition in slope.

Here, we investigate in detail the structure of the density profile for each $p$, and try to describe it with the following functional form:

which behaves like $\rho\propto r^{\alpha(p)}$ at $r\ll S(p)$ and $\rho\propto r^{\beta(p)}$ at $r\gg S(p)$, under the assumption of $\alpha(p)<0$ and $\beta(p)<0$. Here, the four parameters, $A(p)$, $S(p)$, $\alpha(p)$ and $\beta(p)$, are determined by fitting equation (8) to the measured density profile for each $p$ following the method introduced in Section 3.3.

The results are shown by black dashed lines in the upper panels of Figs. 7–9. Here, the fitting to equation (8) was performed for profiles with $p\geq 3$, since the profiles for $p=1$ and $2$ clearly show non-double-power law features due partly to the contributions of substructure, as seen in the radial phase-space distribution. Figs. 7–9 show that the fitting form of equation (8) describes the measured profile fairly well, especially around the transition scales. In particular, a better fitting result is obtained as increasing $p$. For a smaller value of $p$, a deviation from the double-power law function is found at the inner part of massive halos (see Figs. 7 and 8), and they tend to have a flat core. However, statistical variations seem large due to the lack of particles, and we cannot judge whether this deviation is significant or not on the basis of individual haloes. We will thus consider stacked halo profiles, and quantify statistically the goodness of fit to the double-power law function.

In this subsection, using the selected $683$ halos, we analyse stacked density profiles for four halo mass bins S – XL summarised in Table 2. Fitting them to the double-power law function, statistical properties of the radial profile for different $p$ are elucidated, particularly focusing on their mass dependence.

In Section 4.2.1, fitting to the double-power law profile is examined for stacked profiles of each $p$, and we show that the fitting function with $\alpha=-1$ and $\beta=-8$ reproduces well the measured profiles for all mass bins. With these fixed slopes, we quantitatively investigate the dependence of the fitting parameters, i.e., characteristic density $A(p)$ and radius $S(p)$, on the number of apocentre passages $p$ and halo mass in Section 4.2.2.

So far, we have studied individual stream profiles for the mass-selected halo samples and characterised their properties based on the double-power law profile. Here, to assess the robustness of our findings, we analyse a subset of $460$ halos within a specific mass range $[4.10\times 10^{11},\,2.39\times 10^{12}]\,h^{-1}M_{\odot}$, which are divided into two sub-samples based on two different criteria. One is the concentration parameter $c_{\rm vir}$, defined by the ratio $R_{\rm vir}/R_{\rm s}$. Here, the radius $R_{\rm s}$ is the scale radius of the NFW profile (equation (12)), and we estimate it from Rockstar based on the maximum circular velocity (Klypin et al., 2011). Another quantity used for sample selection is the mass accretion rate defined by

where the quantity $t_{\rm dyn}$ represents the dynamical time estimated from halo masses (Diemer, 2017)¹¹1 We use the virial mass, $M_{\mathrm{vir}}$, to measure $\Gamma_{\mathrm{dyn}}$, whereas Diemer (2017) uses $M_{\mathrm{200m}}$.:

which gives $4.04$ Gyr at $z=0$ in our case.

In both cases, we divide the halos into two halves, one with high values of these indicators and the other with low values. Fig. 16 shows the distribution of the parameters, $c_{\rm vir}$ and $\Gamma_{\rm dyn}$, and $M_{\rm vir}$. Darker points are the halos used in the present analysis, and the boundary of the samples is shown in red lines in the upper and lower left panels. Clearly, the distribution in the measured parameters $c_{\rm vir}$ and $\Gamma_{\rm dyn}$ exhibits a tight correlation, and these parameters are anti-correlated. That is, halos with a higher concentration parameter tend to have a smaller $\Gamma_{\rm dyn}$. Therefore, we expect these two cuts to affect the results in a similar manner.

Fig. 17 shows the results of density profiles for the subsamples defined by $c_{\rm vir}$ (upper) and $\Gamma_{\rm dyn}$ (lower). In each panel, results for low and high values of $c_{\rm vir}$ or $\Gamma_{\rm dyn}$ are depicted as red and black colours, respectively. The stacked profiles for each $p$ are fitted to the double-power law form in equation (9), and the best-fit results are plotted in dotted and solid curves for the two subsamples, respectively. We again see a good agreement between the double-power law function and measured profiles over a wide range of $p$ and radius.

A close look at each stream profile in Fig. 17 reveals that halos with high concentration or low accretion rate tend to have a large amplitude $A(p)$ and a large characteristic scale $S(p)$, in particular for $p\gtrsim 14$. The opposite trends in the samples divided with the parameters $c_{\rm vir}$ and $\Gamma_{\rm dyn}$ simply come from their anti-correlation behavior seen in the lower right panel of Fig. 16. On the other hand, the result that high-concentration halos have large $S(p)$ seems somewhat counter-intuitive in the sense that a large value of $c_{\rm vir}$ implies, by definition, a small scale radius $R_{\rm s}$. Nevertheless, the mass inside the radii $r\leq R_{\rm s}$ actually increases as increasing $c_{\rm vir}$. Since the inner mass of the halo is mainly determined by the sum of stream profiles for large values of $p$, and the mass of each stream profile is proportional to $A(p)S(p)^{3}$, a high mass concentration leads to a large value of $A(p)S(p)^{3}$, consistent with what is seen in Fig. 17²²2Assuming the double-power law form of equation (9), the mass of each stream profile, $M(p)$, is analytically expressed as follows: $\displaystyle M(p)$ $\displaystyle=\int_{r=0}^{r=+\infty}\frac{4\pi r^{2}A(p)}{\bigl{\{}r/S(p)\bigr{\}}\,\Bigl{[}1+\bigl{\{}r/S(p)\bigr{\}}^{7}\Bigr{]}}dr$ $\displaystyle=\frac{4\pi}{49}\left(2\sin{\frac{\pi}{7}-\cos{\frac{\pi}{14}}+3\cos{\frac{3\pi}{14}}}\right)A(p)S(p)^{3},$ which gives $M(p)\approx 0.574A(p)S(p)^{3}$. .

The results shown in Figs. 17 suggest that the double-power law feature in the radial stream profiles generically appears, irrespective of the halo sample selection. Although the characteristic density and scale, $A(p)$ and $S(p)$, depend generally on the selection criteria, their $p$ dependence can be translated from one to another once the relationship between different samples is established. In this respect, it might be useful to re-derive the fitting formulas of $A(p)$ and $S(p)$ expressed in terms of $c_{\rm vir}$ rather than $M_{\rm vir}$, since the concentration parameter is tightly correlated with the inner structure of haloes.

To do this, we further divide the haloes according to the value of $c_{\rm vir}$. We consider four subgroups: from $10$ to $30$, $30$ to $50$, $50$ to $70$ and $70$ to $90$ percentiles of $c_{\rm vir}$, each containing the same number of halos. We then look for a fitting function for $A(p)$ and $S(p)$, now with dependence on $c_{\rm vir}$. We find the following functions give reasonable fit to the data:

The performance of the new fitting formulas given above is examined in Fig. 18, where the total density profile and their logarithmic slope are again plotted for the four subgroups. Similarly to the case of mass-selected samples in Fig. 15, the sum of the profiles with the parameters given by the formulas in equations (16) and (17) reproduces the total profile measured from the HR run well beyond the resolution limit of the LR run, which we actually use to calibrate the fitting functions.

The results in Section 4 and 5.1 show that the inner structure of halos exhibits common features, and each stream profile is characterised by a simple double-power law function in a rather universal manner, suggesting that haloes evolve in a self-similar fashion. Here, we compare the radial phase-space structure and density profile for each stream with those predicted by self-similar solutions. To be strict, self-similar solutions are valid only in the Einstein-de Sitter universe, and hence a face-to-face comparison with simulations in the $\Lambda$CDM model makes little mathematical sense. Nevertheless, the secondary infall model of Bertschinger (1985), equivalent to the self-similar solution by Fillmore & Goldreich (1984) with a specific parameter choice (see below), has been shown to reproduce the power-law slope of the pseudo-phase-space density, $Q(r)\propto r^{-1.875}$, found in simulations in the $\Lambda$CDM model.

Along the lines of this, we focus on the qualitative trends in simulations discussed so far and compare them with predictions by self-similar solutions. Among various models considered in the literature, we adopt spherically symmetric solutions by Fillmore & Goldreich (1984) and Sikivie et al. (1997) (see also Nusser (2001)) as representative models of halo evolution under stationary matter accretion. Note that the latter model allows us to incorporate the non-zero angular momentum into the matter flow. The number of parameters characterising the halo structure is two. One is the slope of the initial linear overdensity, parameterised as $\delta_{\rm lin}(r)\propto r^{-3\epsilon}$ in the range $0\leq\epsilon\leq 1$, which is related to the accretion rate through $M_{\rm ta}(t)\propto\{a(t)\}^{1/\epsilon}$ with $M_{\rm ta}$ being the mass inside the turn-around radius. Another parameter is the dimensionless angular momentum $\mathcal{L}$, defined by $\mathcal{L}\equiv L/(GM_{*}r_{*})^{1/2}$, with $r_{*}$ and $M_{*}$ being respectively the turn-around radius and the total mass inside the turn-around radius for each spherical shell accreting onto the halo.

The bottom panels of Fig. 19 plot the radial density profile (upper) and phase-space distribution (lower) predicted by the self-similar solution for specific parameter choices, $(\epsilon,\,L)=(1,\,0)$ (left), $(1/6,\,0)$ (middle) and $(1/6,\,1/10)$ (right). Note that the second solution corresponds to the secondary infall model of Bertschinger (1985) as we mentioned above. These are obtained numerically based on the method described by Zukin & Bertschinger (2010). For reference, we also show in the upper panel the results of a representative halo taken from Fig. 8.

Because of the spherical symmetry, the radial phase space structure of the self-similar solution exhibits distinctive streams, and all of the plotted results show a qualitatively similar trend at the outer part. That is, for each stream, a spiky structure appears in density around the caustics, and this is shortly followed by a sharp drop. On the other hand, the inner structure of self-similar solutions generically shows a power-law behavior, and in the cases with zero angular momentum, the slope of the total density profile as well as each stream profile is predicted to be nearly $-2$. To be strict, the asymptotic slope of the solution with $(\epsilon,\,\mathcal{L})=(1,\,0)$ is slightly different from $-2$: the total and each stream approach $\rho\propto r^{-9/4}$ and $r^{-15/8}$, respectively³³3To derive the asymptotic slope of each stream profile, we first note that Fillmore & Goldreich (1984) and Nusser (2001) analytically derived the asymptotic inner slope of total density profiles $\rho_{\rm tot}$ for self-similar solutions, summarised as: $\displaystyle\rho_{\rm tot}(r)\propto\begin{cases}{r^{-9\epsilon/(1+3\epsilon)}\,(\frac{2}{3}\leq\epsilon\leq 1)}\\ {r^{-2}\,(0\leq\epsilon\leq\frac{2}{3})}\end{cases}$ for $L=0$, and $\displaystyle\rho_{\rm tot}(r)\propto r^{-9\epsilon/(1+3\epsilon)}\,(0\leq\epsilon\leq 1)$ for $L\neq 0$. Then, consider the interior mass of particles/shells having the number of apocentre passages $p$ inside the radius $r$, which we denote by $M_{p}(r)$. For a stationary halo composed of equal mass particles/shells, this is proportional to the time $\delta t(r)$ that a particle/shell with $p$ spends inside $r$ (Fillmore & Goldreich, 1984; Dalal et al., 2010; Lithwick & Dalal, 2011). Except for the radii near the apoapsis, the time $\delta t(r)$ is proportional to $r/v_{\rm r}$ with $v_{\rm r}$ being the radial velocity. Since $v_{\rm r}\propto\sqrt{\Phi}$, we have $M_{p}(r)\propto r/\sqrt{\Phi(r)}$. For the total profile that has a slope less than or equal to $-2$, corresponding to the solutions with $\epsilon\leq 2/3$, the potential becomes constant near the centre, and we obtain $M_{p}(r)\propto r$, which gives the stream profile of $r^{-2}$. On the other hand, if the inner slope of the total profile is steeper than $-2$, corresponding to the solutions with $\epsilon>2/3$, the potential diverges at the centre and we have $\Phi\propto r^{(2-3\epsilon)/(1+3\epsilon)}$. Taking this dependence into account, the interior mass with $p$ becomes $M_{p}(r)\propto r^{9\epsilon/2(1+3\epsilon)}$, which yields the stream profile of $\propto r^{-3(2+3\epsilon)/2(1+3\epsilon)}$. Hence, for $(\epsilon,\,{\mathcal{L}})=(1,\,0)$, we obtain the slope of the stream profile, $-15/8$.. Setting aside a subtle difference, the predicted steep slopes for the total and stream profiles are rather contrasted with those found in the simulations. In the case of $(\epsilon,\,\mathcal{L})=(1/6,\,1)$, the non-zero angular momentum yields a potential barrier near the centre in each mass shell, making the pericentre radius finite. As a result, a sharp cutoff appears at the inner radii for each stream profile, and this results in the total profile being shallower than the slope of $-2$. Thus, apart from small bumpy structures, the resultant total profile resembles the one obtained from simulations. However, the density profile of each stream still possesses a steep slope close to $-2$ above the pericentre radius, which contradicts the asymptotic inner slope seen in our simulation results.

These results indicate that the self-similar solutions considered here miss something essential or some complexities inherent in the cosmological $N$-body simulations. Therefore, a more comprehensive study is necessary taking into account properly the missing ingredients in the self-similar solutions. This may involve relaxing the symmetry assumptions (Ryden, 1993; Lithwick & Dalal, 2011) or exploring non-zero tidal torques (Zukin & Bertschinger, 2010). Also, the angular momentum distribution may be the key. Although we have considered the self-similar solution with angular momentum, this allows for the angular momentum to be provided in a very specific manner. Introducing a broad angular momentum distribution yields a new radial dependence of the particle trajectories that potentially leads to a shallower stream profile consistent with simulations. This could be achieved perhaps by further breaking the self-similarity, as studied by Lu et al. (2006). We leave these investigations to future work.

As a final discussion toward a better understanding of the origin of the universal double-power law nature, we focus on the halo sample M in Table 2, and select the particles with $p=2$, $4$, $6$, $8$, $10$, $20$, $30$ and $40$ at $z=0$. Then we trace their trajectories to higher redshifts and measure the density profiles for each value of $p$ stacked over different halos. Here, the values of $p$ always refer to those counted until $z=0$ instead of the redshift at which the profiles are plotted. Namely, we investigate the time evolution of the profile for the same set of particles over time. Fig. 20 overplots the results at $z=0.3$ (green) and $1.6$ (red), on top of those at $z=0$ already shown in Fig. 10 (blue). We observe that the amplitude of the curve increases as the redshift decreases. Interestingly, however, the evolution of the inner profiles becomes significantly weaker as the value of $p$ increases, and at $p=40$, the profiles almost converge even at the outermost part. This suggests that the double-power law nature was established at an early stage of halo formation and remains stable against matter accretion, which can only affect the outer part of the density profile represented by particles with $p$ smaller than or equal to $6$. Apart from the origin of the universal profile, this picture is consistent with previous studies that show that the accreting matter mainly accumulates in the outer region (e.g., Fukushige & Makino, 2001; Zhao et al., 2003; Wang et al., 2011), and partly explains why the characteristic scale $S(p)$ in equation (11) is a decreasing function of $p$; particles with larger $p$ have accreted earlier and their distribution tends to be relaxed in the inner part of haloes. In this respect, the dynamics at the early stage of halo formation would clarify the origin of the double-power law nature.