Splashback Radius in a Spherical Collapse Model

Abstract

It was shown several years ago that dark matter halo outskirts are characterized by very steep density profiles in a very small radial range. This feature has been interpreted as a pile-up of different particle orbits at a similar location, namely, splashback material at half an orbit after collapse. Adhikari et al. (2014) obtained the location of the splashback radius through a very simple model by calculating a dark matter shell trajectory in the secondary infall model while it crosses a growing NFW profile-shaped dark matter halo. Because they imposed a halo profile instead of calculating it from the trajectories of the shells of dark matter, they were not able to find the dark matter profile around the splashback radius. In the present paper, we use an improved spherical infall model taking into account shell crossing as well as several physical effects such as ordered and random angular momentum, dynamical friction, adiabatic contraction, etc. This allows us to determine the density profile from the inner to the outer region and to study the behavior of the outer density profile. We compare the density profiles and their logarithmic slope of with the simulation results of Diemer and Kravtsov (2014), finding a good agreement between the prediction of the model and the simulations.

1. Introduction

The problem of determining the structure of dark matter haloes is an old one, and has been studied from analytical, numerical, and observational points of view. The first efforts were based on analytical models, in particular, on the spherical collapse model. The first trials to study the formation of virialized structure go back to the seminal paper by [1], after which several authors investigated the consequences of secondary infall and accretion onto proto-structures, studying in particular the structure of the density profiles. In [2], the authors studied the self-similar collapse of scale-free perturbations, determining important processes in the halo profiles. Several other papers have improved on these results, e.g., [3,4,5,6,7,8,9]. These and other studies improved the secondary infall model (SIM) by taking into account the effects of ordered and random angular momentum, adiabatic contraction of dark matter (DM) produced by baryons [10,11,12], the effects of dynamical friction, and many more effects. Further developments were obtained by means of DM-only N-body simulations with the determination of the Navarro–Frenk–White profile [13,14], the Einasto profile [15,16], and hydrodynamic simulations [17,18]. Much of the past work focused on the inner structure of haloes and was driven by efforts to understand and solve the cusp–core problem, that is, the discrepancy between the steep slopes predicted by DM-only simulations and the flat-cored profiles observed in low surface brightness (LSB) and dwarf galaxies. In the last several years, authors (e.g., [19]) have begun to study the outer part of haloes, finding that the outer profiles are inconsistent with typical fits such as the NFW and Einasto profiles. Findings include that the outer density profiles are characterized over a narrow range of radii by very steep logarithmic slopes of $dlog\rho /dlogr\le -4$. According to [19], the observed local steepening is due to a caustic related to the splashback of material accreted by the halo. The presence of caustics is due to the pile-up of different particle orbits at a similar location. The splashback radius corresponds to the outermost caustic associated with the first apoapse after collapse. As shown in [20], the splashback location is provided by the relation $dlog\rho /dlogr\le -3$. Caustics are not a rare phenomenon; they are present in the [2] similarity solutions, in 3D similarity solutions of the collapse of triaxial peaks [21], and in real galaxies, where they appear as radial shells [22]. However, detecting density enhancement related to caustics in N-body simulations of dark haloes is not easy. This is mainly due to the presence of small-scale structures that can smear out caustics [23,24]. As material accumulates, a steepening in the outer profile is observed at the splashback radius. A few years ago, [20] estimated the location of the splashback radius using a very simplified secondary collapse model to obtain the secondary infall trajectory of DM shells by means of a growing DM halo with an NFW profile. In their calculation, a shape of the halo profile was imposed instead of being computed from the trajectories of the DM shells. As a consequence, in their calculation they were not able to obtain the full shape of the DM profile around the splashback radius. In the following, we introduce a much more complex spherical collapse model which allows us to simultaneously calculate the trajectories and the DM halo profile.

The paper is organized as follows. In Section 2, we discuss the model used to determine the density profile and which allows us to determine of the features of the outer profile. In Section 4 we discuss the results following from the model and compare them with the results of [19]’s simulations. Section 5 is devoted to further discussion.

2. Model

In this section, we discuss the model that allows us to determine the density profile. It was first introduced in [8], followed by several applications: to density profiles universality studies [25,26], to galaxies [27,28] and clusters [28,29] density profiles, and to galaxies’ inner surface-density distributions [30].

The semi-analytical model (SAM) used in this work encompasses several upgrades on the SIM (e.g., [1,3,4,5,6,7,9,31,32]). In contrast to the anterior avatars of the SIM, it comprises non-radial collapse effects from random angular momentum (RAM [5,6,33,34,35,36,37,38,39,40])1, ordered, tidal angular momentum [41,42], and the impact of dynamical friction (e.g., [8,43,44,45]) and of baryonically-induced DM adiabatic contraction [10,11,12].

This SAM evolves perturbations from their linear expansion with the Hubble flow to turn-around before collapsing, with adiabatic central potential variations including “shell-crossing” [2,46].

Spherical SIMs2 in the filiation from [1] describe a bound mass shell expansion from a comoving initial radius $x_i$ to its maximum (turn-around or apapsis) radius $x_m=x_{ta}$,

$$x_m=g\left(x_i\right)=x_i/{\overline{\delta }}_i,$$

(1)

with the linearly grown mean overdensity inside the shell extrapolated at current epoch $t_0$ and with $\delta \left(y\right)$ obtained as in Appendix B of [8], results from

$${\overline{\delta }}_i=\frac{3}{x_i^3}{\int }_0^{x_i}\delta \left(y\right)y^{2}dy.$$

(2)

More generally [47], Equation (1) extends to a pure dust Universe, with the density parameter ${\Omega }_i$ as

$$x_m=g\left(x_i\right)=x_i\frac{1+{\overline{\delta }}_i}{{\overline{\delta }}_i-({\Omega }_i^{-1}-1)}.$$

(3)

This generalisation represents the core of the SIM, for which a Lagrangian shell’s time averaged radius remains proportional to its initial radius. Using Equation (3), the final radius x can be written as proportional to the turn around radius, $x_m$:

$$x=f\left(x_i\right)x_m$$

(4)

with the scaling fitted by [48]:

$$\begin{array}{cc}\hfill f=f\left(\alpha \right)=& 0.186+0.156\alpha +0.013{\alpha }^2+0.017{\alpha }^3\hfill \\ & -0.0045{\alpha }^4+0.0032{\alpha }^5.\hfill \end{array}$$

(5)

Mass conservation at turn around radius yields the density profile [3,35,47]

$${\rho }_{ta}\left(x_m\right)={\rho }_i\left(x_i\right){\left(\frac{x_i}{x_m}\right)}^2\frac{d{x}_i}{d{x}_m}.$$

(6)

Beyond turnaround, shell crossing effects are bypassed with the Virial theorem, yielding a collapse factor $f=0.5$, resulting in the final density profile through mass conservation:

$$\rho \left(x\right)x^{2}dx={\rho }_i{x}_i^{2}d{x}_i.$$

(7)

This produces the power-law density profile from [3] from the shell’s initial density approximation

$${\rho }_i\left(x_i\right)={\rho }_{b,i}[1+{\delta }_i\left(x_i\right)],$$

(8)

and the linear ${\delta }_i$ expansion of Equation (3), reading

$$\rho \left(x\right)\propto x^{-3(n+3)/(n+4)}.$$

(9)

However, the Virial theorem relies on energy conservation, and oscillations of collapsing shells through their inner shells break up the energy integral of motion and vary the value of f.

This modifies the SIM dynamics, which assume a “gentle” collapse, as follows: with the conjecture of adiabatic variation of the central potential [2,46], shells near the centre oscillate many times without significant changes in potential. In other words, inner shells’ orbital period can be neglected compared with outer shells’ collapse time [49]. In this case, the inner shells admit the radial action $\oint v\left(r\right)dr$, with the radial velocity $v\left(r\right)$ as adiabatic invariant. The collapse of outer shells slowly changes the potential, shrinking the inner shells via the radial action invariant.

The mass inside a shell with initial radius $x_i$ at its apocenter (apapsis radius) $x_m$ can be decomposed in the sum of its inner shells masses, i.e., with apocenters inside $x_m$, the permanent component $m_p$; with the contribution from its temporarily crossing outer shells masses, the additional mass is $m_{add}$. Mass conservation yields the first component,

$$m_p\left(x_m\right)=m\left(x_i\right)=\frac{4}{3}\pi {\rho }_{b,i}{x}_i^3(1+{\overline{\delta }}_i),$$

(10)

with initial time constant density ${\rho }_{b,i}$ of the homogeneous Universe. From Equation (6), the distribution of mass $m\left(x\right)=m\left(x_m\right)$ and the system radius R, its outer shell apapsis, together with the probability $P_{x_m}\left(x\right)$ to find the shell with apapsis x inside the radius $x_m$, follows the additional component

$$m_{add}\left(x_m\right)={\int }_{x_m}^R{P}_{r_m}\left(x\right)\frac{dm\left(x\right)}{dx}dx,$$

(11)

with the resulting total mass reading

$$m_T\left(x_m\right)=m_p\left(x_m\right)+m_{add}\left(x_m\right).$$

(12)

The ratio of the time spent below $x_m$ to the outer shell’s x oscillation period allows $P_{x_m}\left(x\right)$ to be computed from the outer shell’s pericenter $x_p$ and $v_x\left(\eta \right)$, its radial velocity at radius $\eta$, as follows:

$$P_{x_m}\left(x\right)=\frac{{\int }_{x_p}^{x_m}\frac{d\eta }{v_x\left(\eta \right)}}{{\int }_{x_p}^x\frac{d\eta }{v_x\left(\eta \right)}}.$$

(13)

The outer shell’s radial velocity derives from integrating its equation of motion, including the tidal torque-generated ordered specific angular momentum3 $h(r,\nu )$, the random angular momentum $j(r,\nu )$ (see [8,33] and follow ups), the gravitational potential acceleration $G\left(r\right)$, $\Lambda$ (the cosmological constant), and the dynamical friction coefficient $\mu$:

$$\frac{d{v}_r}{dt}=\frac{h^2(r,\nu )+j^2(r,\nu )}{r^3}-G\left(r\right)-\mu \frac{dr}{dt}+\frac{\Lambda }{3}r.$$

(14)

Computations of $\mu$ and the angular momenta are explained in ([8], Appendices C and D). In the restriction where $\mu =0$, Equation (14) integrates into the square of velocity:

$$v{\left(r\right)}^2=2\left[ϵ-G{\int }_0^r\frac{m_T\left(y\right)}{y^2}dy+{\int }_0^r\frac{h^2}{y^3}dy+\frac{\Lambda }{6}{r}^2\right],$$

(15)

where the shell’s specific binding energy results from the turnaround value at $v\left(r\right)=dr/dt=0$ at $r=x_m$ in Equation (15).

The $\mu \ne 0$ case requires the numerical integration of Equation (14) for v after its quadrature into

$$\frac{d{v}^2}{dt}+2\mu v^2=2\left[\frac{h^2+j^2}{r^3}-G\frac{m_T}{r^2}+\frac{\Lambda }{3}r\right]v.$$

(16)

With the above computations complete, following [2,46] we obtain $f\left(x_i\right)$ the shell’s collapse factor, which starts at radius $x_i$ and reaches apapsis $x_m$

$$f\left(x_i\right)=\frac{m_p\left(r_m\right)}{m_p\left(r_m\right)+m_{add}\left(r_m\right)},$$

(17)

from Equations (6) and (4), the corresponding density profile at Virialisation is

$$\rho \left(x\right)=\frac{{\rho }_{ta}\left(x_m\right)}{f^3}{\left[1+\frac{dlnf}{dlng}\right]}^{-1}.$$

(18)

The calculations above allow us to evaluate the $f\left(x_i\right)$ variations from energy integral break up, confirmed from N-body simulations, finding its relation to the initial density perturbation profile and its increase with initial radius. The case $x_i\to 0$ and $f\to 0$ recovers a radial collapse [50]. The computation from Equation (17) of f via integration of Equation (11) can proceed numerically after variable change to express them in terms of initial radius [38,50]. Such variable change turns Equation (11) into

$$m_{add}\left(r_m\right)=4\pi {\rho }_{b,i}{\int }_{x_i}^{x_b}{P}_{x_i}\left(x_i^{\prime }\right)[1+{\delta }_i\left(x_i^{\prime }\right)]x_i^{\prime 2}\mathrm{d}{x}_i^{\prime },$$

(19)

with $P_{x_i}\left(x_i^{\prime }\right)=I\left(x_i\right)/I\left(x_i^{\prime }\right)$,

$$I\left(r\right)={\int }_{x_p^{\prime }}^r\frac{1}{v_{g\left(x_i^{\prime }\right)}\left(g\left(\eta \right)\right)}\frac{\mathrm{d}g\left(\eta \right)}{\mathrm{d}\eta }\mathrm{d}\eta ,$$

(20)

taking $r_m=g\left(x_i\right)$, $x_p^{\prime }=g^{-1}\left(x_p\right)$ and the initial $x_i^{\prime }$ shell’s pericenter $x_p$. The upper bound $x_b$ of Equation (19)’s integration corresponds to the presently collapsed sphere’s initial radius. A similar variable change in Equation (16) leads to the determination of the radial velocity v for a shell characterised with apapsis $x=g\left(x_i\right)$ at radius $r=g\left(r_i\right)$ through

$$\frac{d{v}_x^2\left(r\right)}{dt}+2\mu v_x^2=2\left[\frac{h_x^2+j_x^2}{r^3}-\Psi \left(r\right)+\frac{\Lambda }{3}r\right]v_x\left(r\right),$$

(21)

with the shifted gravitational potential $\Psi$ from the initial mass profile $m\left(x_i\right)$ reading

$$\Psi \left[g\left(r_i\right)\right]=\frac{Gm\left(x_b\right)}{g\left(x_b\right)}+G{\int }_{r_i}^{x_b}\frac{m\left(x_i\right)}{g^2\left(x_i\right)}\frac{\mathrm{d}g\left(x_i\right)}{\mathrm{d}{x}_i}\mathrm{d}{x}_i.$$

(22)

In summary (see [50], Section 4), the equation of motions of a shell (21), given the angular momentum distribution, dynamical friction coefficient, and initial conditions, integrates to compute $P_{x_i}$, the probability from Equation (20), which then yields the transient part of the gravitating mass acting on the shell $m_{add}$ from Equation (19) and the collapse factor f from Equation (17) to obtain the final density profile through Equation (18) (as in [5], Section 2.1). Our model’s ordered angular momentum formation follows ([8], Appendix C1), its random angular momentum computation agrees with ([8], Appendix C2), its dynamical friction coefficient and the baryon dissipative collapse are described in ([8], Appendix D), the baryon’s adiabatic contraction is the object of ([8], Appendix E), and [Appendix B [8] describes initial condition generation. More explicitly, the model’s “ordered angular momentum” h derives from the tidal torque theory (TTT) [41,42,51,52,53], which uses processes from the tidal torques exerted on smaller-scale structures by larger-scale objects. On the other hand, the “random angular momentum” j is computed from the orbital axis ratio between the pericentric and apocentric radii $r_{\min }$ and $r_{\max }$; note that $e=\left(\frac{r_{\min }}{r_{\max }}\right)$ [54], modified according to the system’s dynamical state following simulations [5] into

$$e\left(r_{\max }\right)\simeq 0.8{\left(\frac{r_{\max }}{r_{\mathrm{ta}}}\right)}^{0.1},$$

(23)

a function of the ratio of $r_{\max }$ to the turnaround radius, where $r_{\max }<0.1r_{\mathrm{ta}}$. As for the dynamical friction effects, they are introduced into the equation of motion with a dynamical friction force, as computed in ([8], Appendix A; see Equation (A14)). Finally, for density profile steepening from adiabatic compression, the methods in [11] were followed.

3. Population of Haloes

The model is the principle that transforms from initial conditions to charateristics of the galaxy. Generating a population of galaxies requires use of the model as well as an initial range of parameters, a population of initial parameters, which provides the statistical distribution of galaxies given a sufficiently large number of the generated population. For the initial conditions and determination of the density profile of galaxies, it is necessary to calculate the initial overdensity, ${\overline{\delta }}_i\left(x_i\right)$, which can be calculated when the spectrum of perturbations is known. It is widely accepted that structure formation in the universe is generated through the growth and collapse of primeval density perturbations originating from quantum fluctuations in an inflationary phase of early Universe. The growth in time of small perturbations is due to gravitational instability. The statistics of density fluctuations originating in the inflationary era are Gaussian, and can be expressed entirely in terms of the power spectrum of the density fluctuations:

$$P\left(k\right)=\langle |{\delta }_{\mathbf{k}}{|}^2\rangle$$

(24)

where

$${\delta }_{\mathbf{k}}=\int d^{3}kexp(-i\mathbf{k}\mathbf{x})\delta \left(\mathbf{x}\right)$$

(25)

$$\delta \left(\mathbf{x}\right)=\frac{\rho \left(\mathbf{x}\right)-{\rho }_b}{{\rho }_b}$$

(26)

and ${\rho }_b$ is the mean background density. In biased structure formation theory it is assumed that cosmic structures of linear scale $R_f$ form around the peaks of the density field $\delta \left(\mathbf{x}\right)$, and are smoothed on the same scale. According to the hierarchical scenario of structure formation, haloes should collapse around maxima of the smoothed density field; see below. The statistics of peaks in a Gaussian random field were studied in the classic paper by ([55], hereafter BBKS). A well known result is the expression for the radial density profile of a fluctuation centered on a primordial peak of arbitrary height $\nu$:

$$\langle \delta \left(r\right)\rangle =\frac{\nu \xi \left(r\right)}{\xi {\left(0\right)}^{1/2}}-\frac{\vartheta (\nu \gamma ,\gamma )}{\gamma (1-{\gamma }^2)}\left[{\gamma }^2\xi \left(r\right)+\frac{R_{*}^2}{3}{\nabla }^2\xi \left(r\right)\right]·\xi {\left(0\right)}^{-1/2}$$

(27)

(Refs. [33,55], with [33] hereafter RG87), where $\nu =\delta \left(0\right)/\sigma$ (see the following for a definition of $\sigma$) is the height of a density peak and $\xi \left(r\right)$ is the two-point correlation function

$$\xi \left(r\right)=\frac{1}{2{\pi }^{2}r}{\int }_0^{\infty }P\left(k\right)ksin\left(kr\right)dk,$$

(28)

$\gamma$ and $R_{*}$ are two spectral parameters respectively provided by

$$\gamma =\frac{\int k^{4}P\left(k\right)dk}{{\left[\int k^{2}P\left(k\right)dk\int k^{6}P\left(k\right)dk\right]}^{1/2}}$$

(29)

$$R_{*}={\left[\frac{3\int k^{4}P\left(k\right)dk}{\int k^{6}P\left(k\right)dk}\right]}^{1/2}$$

(30)

and $\vartheta (\gamma \nu ,\gamma )$ is

$$\vartheta (\nu \gamma ,\gamma )=\frac{3(1-{\gamma }^2)+\left(1.216-0.9{\gamma }^4\right)exp\left[-\left(\frac{\gamma }{2}\right){\left(\frac{\nu \gamma }{2}\right)}^2\right]}{{\left[3\left(1-{\gamma }^2\right)+0.45+{\left(\frac{\nu \gamma }{2}\right)}^2\right]}^{1/2}+\frac{\nu \gamma }{2}}$$

(31)

Then, ${\overline{\delta }}_i$ is calculated from Equation (27) similarly to [5] (see their Section 2.2). In order to calculate $\delta \left(r\right)$ we need a power spectrum, $P\left(k\right)$. The CDM spectrum used in this paper is that of BBKS, with the transfer function

$$\begin{array}{cc}\hfill T\left(k\right)=& \frac{[ln\left(1+2.34q\right)]}{2.34q}{·[1+3.89q+\left(16.1q\right)}^2+{\left(5.46q\right)}^3\hfill \\ \hfill & +{\left(6.71\right)}^4{]}^{-1/4}\hfill \end{array}$$

(32)

where $q=\frac{k{\theta }^{1/2}}{{\Omega }_{\mathrm{X}}{h}^2{\mathrm{Mpc}}^{-1}}$. Here, $\theta ={\rho }_{\mathrm{er}}/\left(1.68{\rho }_{\mathrm{fl}}\right)$ represents the ratio of the energy density in relativistic particles to that in photons ($\theta =1$ corresponds to photons and three flavors of relativistic neutrinos). The spectrum is connected to the transfer function through the equation

$$P\left(k\right)=P_{CDM}{e}^{-1/2k^2{R}_f^2}$$

(33)

where $R_f$ is the smoothing (filtering) scale and $P_{CDM}$ is provided by

$$P_{CDM}=Ak{T}^2\left(k\right)$$

(34)

where A is the normalization constant. We normalized the spectrum by imposing the mass variance of the density field as

$${\sigma }^2\left(M\right)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }dk{k}^{2}P\left(k\right)W^2\left(kR\right)$$

(35)

convolved with the top hat window

$$W\left(kR\right)=\frac{3}{{\left(kR\right)}^3}(sinkR-kRcoskR)$$

(36)

of radius 8 $h^{-1}$$Mp{c}^{-1}$ being ${\sigma }_8=0.76$. Throughout this paper, we adopt a $\Lambda$CDM cosmology with WMAP3 parameters [56], ${\Omega }_m=1-{\Omega }_{\Lambda }=0.24$, ${\Omega }_{\Lambda }=0.76$, ${\Omega }_b=0.043$, and $h=0.73$, where h is the Hubble constant in units of 100 km $s^{-1}$$Mp{c}^{-1}$.

The mass enclosed in $R_f$ is calculated, as in RG87, as $M=4\pi /3{\rho }_b{R}_f^3$, such that for $R_f=0.12$ Mpc we can say that $M\simeq {10}^9{M}_{\odot }$.4 Structure such as Galaxies form from high peaks in the density field, high enough that they stand out above the “noise" and dominate the infall dynamics of the surrounding matter.

The amplitude of a given peak is expressed in terms of its $\sigma$ deviation, where $\sigma =\xi {\left(0\right)}^{1/2}$. Thus, the central density contrast of an $\nu \sigma$ peak is $\nu \xi {\left(0\right)}^{1/2}$ and the peak height is provided by $\nu =\delta \left(0\right)/\sigma$. As galaxies are rather common, they must have formed from peaks that are not very rare, say, 2–4 $\sigma$ peaks (RG87). In [8] (Figure 6 of Appendix B), the authors show the density profiles $\delta \left(x\right)$ plotted for $\nu =2,3,$ and 4. We generated a set of galaxies starting from the initial conditions and using the model. The different masses are related to the filtering radius $R_f$ and peak height $\nu$. As seen in Figure 6 of Appendix B in [8], the value of $\delta \left(x\right)$, and thus the final density profile, changes with changing $\nu$. The halo characteristics are modified by the tidal torque (ordered angular momentum), random angular momentum (Equation (C17) in [8]), dynamical friction (Appendix D in [8]), and adiabatic compression (Appendix E in [8]).

In summary, the statistics of our model’s halo populations are provided by the BBKS Gaussian random field fluctuations. The generated sample size is the result of computing power constraints versus stability of the median density profile, and was set after verifying that an increase in sample size would not significantly affect the results.

4. Results

Before describing the results, we provide several commonly used definitions here. The three-dimensional radius with respect to the center of the halo is indicated by r, while R is used to indicate radii related to the mass of the halo. The critical density is indicated by ${\rho }_{\mathrm{c}}$ and the mean matter density by ${\rho }_{\mathrm{m}}$. Masses at specific overdensities are denoted by $M_{\Delta \mathrm{m}}=M(<R_{\Delta \mathrm{m}})$. For example, the mass with overdensity at $\Delta =200$ is $M_{200\mathrm{m}}$, while that corresponding to the critical density $M_{200\mathrm{c}}$ reads $M_{\Delta \mathrm{c}}=M(<R_{\Delta \mathrm{c}})$. In addition, $M_{\mathrm{vir}}$ and $R_{\mathrm{vir}}$ are related to $\Delta \left(z\right)$, which at $z=0$ corresponds to ${\Delta }_{\mathrm{vir}}(z=0)\approx 358$. Instead of using masses, haloes are binned using the mean values of the peak height $\nu$, as at fixed $\nu$ the halo properties across redshifts should be similar. The definition of the peak height follows the usual

$$\nu \equiv \frac{{\delta }_{\mathrm{c}}}{\sigma (M,z)}=\frac{{\delta }_{\mathrm{c}}}{\sigma (M,z=0)\times D_{+}\left(z\right)},$$

(37)

where the critical density is provided according to the spherical top hat model by ${\delta }_{\mathrm{c}}=1.686$ and $D_{+}\left(z\right)$ represents the linear growth factor normalized to unity at $z=0$. For a sphere of radius R, the rms density fluctuation is provided by

$${\sigma }^2\left(R\right)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }{k}^{2}P\left(k\right){\left|\tilde{W}\left(kR\right)\right|}^{2}dk.$$

(38)

Here, W is the spherical top hat filter function and $\tilde{W}\left(kR\right)$ is its Fourier transform. The linear power spectrum $P\left(k\right)$ is provided by the formula of [57], with the normalization $\sigma \left(8h^{-1}\mathrm{Mpc}\right)={\sigma }_8$ = $0.82$; $\sigma \left(M\right)=\sigma \left(R\right[M\left]\right)$ indicates the variance at a given mass, and the calculation of $\nu$ uses $M=M_{\mathrm{vir}}$. The relation between $\nu$ and the virial mass is represented in ([19], Figure 1). Halo mass can be translated to $\nu$ by means of $R_{\mathrm{vir}}$, as that radius corresponds to the largest for which the density profile scatter remains relatively small at fixed mass. For this reason, we prefer $R_{\mathrm{vir}}$ to $R_{200\mathrm{m}}$ for computing the difference in mass accretion rate between two redshifts. The mean or median profiles in rescaled radial units are obtained by rescaling the individual halo profiles using the $R_{\Delta }$ of the halo. The mean and median are obtained from the rescaled profiles. The slope profiles are obtained using the fourth-order Savitzky–Golay smoothing algorithm over the fifteen nearest bins [58], and the functional fits are obtained by means of the Levenberg–Marquart algorithm. Median profiles are used because they are a good approximation of the typical profile and can be used to study trends in the density profiles.

In this paper, we use the model described in Section 2 to generate a population of haloes from which to build up median density profiles of haloes binned by peak height.

We compare the results of our model with those of the simulations of [19]. The goal is to show that the model provides a correct description of the outer density profile of haloes.

In Figure 1, we plot the median density profiles of two different halo samples, either from [19] (solid lines) or from our model (dotted lines), as well as their logarithmic slope profiles $\gamma \left(r\right)\equiv dlog\rho /dlogr$. The top panels show the median density profiles of low-mass (top left) and very massive (top right) haloes at $z=0$. The shaded band represents the interval around the median containing 68% of the profiles of individual halos in the $\nu$ bin from [19]. The bottom panel represents the logarithmic slope profiles corresponding to the top panels. The NFW and Einasto fits of the profiles are indicated as well. The low mass sample (left panels) is characterized by $0.5<\nu <0.7$, while the high mass sample contains haloes with $\nu >3.5$. The solid lines represent the results from [19]. The dotted lines in the top panels represent the median density profiles while the bottom panels show the logarithmic slope profile, both obtained with our model. The plot shows a good agreement between the profiles predicted by the model and the result of the simulated model [19], only differing by a very small offset and lower numerical noise in our model’s case. The very small difference in slope cannot be perceived in the density plots. The low-$\nu$ sample median profile’s slope changes slowly until r ≳ R_vir. A large scatter is visible at larger radii. The high-ν sample presents a steep profile at r ≳ 0.5R_vir. The slope varies from −2 to −4 over a radius range restricted to ≈4 rescaled radius, as shown in the bottom panel representing the slope profiles. The plot shows that the NFW and Einasto profiles provide a good description of the low-ν sample out to r ≈ R_vir, while the fast steepening of the slope of the high-ν sample is poorly described. The NFW and Einasto profiles present radically different behaviors for the outer density profiles of haloes. They are able to fit the low-ν out to r ≈ R_vir for the high-ν haloes, and r ≈ 0.5R_vir in the case of the high-ν haloes. In order to fit the logarithmic slope profile, it is necessary to use a different fitting formula, such as that presented by [19] (Section 3.3 and Appendix). In the case of both the low-ν and high-ν samples, the profiles flatten to a slope of ≈−1 at r ≳ 2R_vir.

The previous plots in Figure 1 show the profiles of a $\nu$ bin at $z=0$. In general, it could be expected that the profiles of a given $\nu$ when density and radii are rescaled in the correct way would be self-similar in shape. The problem is to find the radii and density to be used for rescaling.

Figure 2 shows the self-similarity of the redshift evolution of the profiles for both [19] (solid lines) and for our model (dotted lines). The left panel displays the redshift evolution of the median density profiles for a peak with $\nu >3.5$ as a function of the radius rescaled by $R_{200m}$, with the density rescaled by ${\rho }_m$. The central panel represents the same density profiles rescaled by $R_{\mathrm{vir}}$ and ${\rho }_{\mathrm{vir}}$, respectively. The right panel represents the same density profiles rescaled by $R_{\mathrm{c}}$ and ${\rho }_{\mathrm{c}}$, respectively. The black lines correspond to $z=0$ and the red lines to $z=6$.

We emphasize that the figure plot profiles are in physical units, rescaled in the three panel from left to right by $R_{200\mathrm{c}}$, $R_{\mathrm{vir}}$, and $R_{200\mathrm{m}}$, respectively, while the densities are rescaled by the corresponding quantities ${\rho }_{\mathrm{c}}$, ${\rho }_{\mathrm{vir}}$, and ${\rho }_{\mathrm{m}}$. The plot shows clearly that the halo structure is approximately self-similar after rescaling with $R_{\Delta }$. It is further evident that the self-similarity depends on the kind of rescaling chosen. The most self-similar inner structures of haloes are obtained by rescaling radii and densities with ${\rho }_{\mathrm{c}}$ and $R_{200\mathrm{c}}$. Conversely, the most-self-similar outer profiles are obtained by rescaling with $R_{200\mathrm{m}}$ and ${\rho }_{\mathrm{m}}$. In order to present a more readable figure, we only plotted the haloes at two redshifts, $z=0$, and $z=6$. We compared the result of our model (dotted line) in Figure 3 with [19], and found again that both sets of results are in agreement, confirming the discussed self-similarity. In order to obtain a better understanding of the self-similarity, it might be possible to use logarithmic slope profiles, as they clearly reveal the radii at which rapid changes in slopes happen.

Figure 3 plots the logarithmic slope of three $\nu$ bins at different redshifts rescaled by $R_{200\mathrm{m}}$ (top row) and $R_{200\mathrm{c}}$ (bottom row). The left panels of Figure 3 present results for the $\nu >3.5$ sample of Figure 2. The central and right panels refer to the $\nu$ bins shown in [19], Figure 5. The radii in the top panels are rescaled in units of $R_{200\mathrm{m}}$, while those in the bottom panel are expressed in units of $R_{200\mathrm{c}}$. The dotted lines represent the result of our model. The solid lines stand for the [19] simulations. The dashed lines provide fits using the Einasto profile. Black and red lines correspond to different redshifts indicated in the legends of the top panels. In all cases, the slopes show a sharp steepening followed by a flattening. In units of $R_{200\mathrm{m}}$, such sharp variations of the slope occur at the same radii with almost no sign of evolution of the transition. The radius of the steepest slope occurs around $\approx 1-1.2R_{200\mathrm{m}}$ for all $\nu$ and redshifts. Furthermore, for haloes rescaled in units of $R_{200\mathrm{m}}$ the outer flattening displays almost no evolution or variation with $\nu$. The situation is different in the case where $r<R_{200\mathrm{m}}$. In this case, a variation of the slopes of the profiles at a given $r/R_{200\mathrm{m}}$ with $\nu$ and z is observed. If radii and densities are rescaled by ${\rho }_c$ and $R_{200\mathrm{c}}$, the opposite is valid. While the shapes of low-$\nu$ and high-$\nu$ profiles are different, in any case they show a certain degree of uniformity at $r>R_{200\mathrm{m}}$ when rescaled by $R_{200\mathrm{m}}$ and at $r<R_{200\mathrm{c}}$ when rescaled by $R_{200\mathrm{c}}$. The previous discussion allows us to reach the conclusion that the inner profiles are self-similar in units of $r/R_{200\mathrm{c}}$ while the outer profiles are self-similar in units of $r/R_{200\mathrm{m}}$. It is interesting to note that one might expect a concentration of haloes to be more universal in terms of $\nu$ if the radius definition were related to the critical density. As in the previous figure, we only plotted two redshift dependencies in order to obtain a more readable plot. Again, the dotted lines represent the result of our model, which is in agreement with that of [19]. Another issue studied by [19] is the dependence on the mass accretion rate. As mentioned, our main goal was to show that the predictions of the simulated model [19] were in agreement with our model, that is, that our model is able to capture the main characteristics of the density profiles and their logarithmic slope to correctly describe the behavior of the outer region of the density profiles. For this reason, we do not discuss the dependence on the mass accretion rate. The comparisons shown here are sufficient to establish that our model can adequately describe the outer region of the density profile.

5. Conclusions

DM halo outskirts are characterized by very different density profiles than the inner density profile, which are usually fitted with models such as the NFW or the Einasto profile. The outer density profile is very steep over a small radial range. This kind of behavior has been interpreted as due to a pile-up of the orbits of particles and splashback of material located near half an orbit after collapse. Modeling spherical haloes, such radii provide a sharp separation of infalling matter from material just orbiting the halo, including its satellites. For exact spherical symmetry, a caustic, in this case an infinitely sharp density drop, characterises the splashback radius. In realistic halos, this caustic is smoothed out. The complexity of orbits in $\Lambda$CDM halos along with their non-sphericity and substructures make splashback radius identification in each halo non-trivial. However, it has been observed in weak lensing as well as around stacked clusters in satellite galaxy density profiles.

Theoretically, the location of the splashback radius can be obtained through the very simple model of [20] by calculating the secondary infall of DM shell trajectories within a growing NFW-profiled DM halo. As they imposed the halo profile a priori instead of calculating it from trajectories of DM shells, they were not able to find the DM profile around the splashback radius. A more complete model was proposed by [59], who extended the self-similar collapse model of [2] to a $\Lambda$CDM universe, allowing for simultaneous calculation of the trajectories and DM halo profile. Furthermore, [19] used simulations to study the density profiles of $\Lambda$CDM haloes, focusing on the outer regions of the halo. They found noteworthy deviations of the outer profile from the classical distributions, such as the NFW and Einasto profiles, in the form of sharper steepening of the density profile’s logarithmic slope than had been predicted. They found that the outermost density profiles at $r\ge R_{200\mathrm{m}}$ are self-similar if the radii are rescaled by $R_{200\mathrm{m}}$. At the same time, the inner density profiles are the most self-similar if the radii are rescaled by $R_{200\mathrm{c}}$.

In the present paper, we have proposed an improved spherical infall model that, unlike the model of [20], takes into account shell crossing and several other effects, including ordered and random angular momentum, dynamical friction, and adiabatic contraction. Using this model, we obtained the halo density profile, and studied the behavior of the profile in its outer regions. From Gaussian initial conditions, we generated populations of haloes that can be statistically compared with other models and observations. A comparison of the profile and the logarithmic slope profile with the results of the simulations in [19] shows good agreement. We were able to reobtain the results of [19] concerning the characteristics of the density profile and its logarithmic slope.

[custom]