ON PHYSICAL SCALES OF DARK MATTER HALOS

Within the currently favored cosmological framework where the energy content is dominated by dark energy and cold dark matter and space is flat (ΛCDM), peaks of the initial matter density field become gravitationally unstable, collapse, then virialize and finally form what we call dark matter halos. Embedded within those dark matter halos is the baryonic matter that forms the galaxies, which we can observe. For a recent detailed and pedagogical review on the theory of cosmological structure formation, see Knobel (2012) and references therein. Therefore, for a comprehensive understanding of the universe, it is of fundamental importance to understand how structures form in cosmology and how dark matter halos and their enclosed galaxies evolve with cosmic time.

In the simplest conceivable scenario, one can study the collapse of a spherically symmetric density perturbation in an otherwise homogeneous universe (Gunn & Gott 1972). This is called the spherical collapse model and we review its essentials in Section 2. When combined with the spectrum of (known) perturbations, this results in an analytic and elegant model for structure formation in the universe. This approach was pioneered in the seminal work of Press & Schechter (1974) and led to the fundamental concepts of how we describe the size and mass of dark matter halos and their abundance (i.e., mass function), which are used to describe the properties of our universe.

The increase of computer power over the last few decades has enabled the possibility of simulating the formation of structure in the universe with ever better and more detailed models. This has led to a good understanding of the large-scale structure of dark matter (e.g., Springel et al. 2005; Boylan-Kolchin et al. 2009; Teyssier et al. 2009; Kim et al. 2011; Klypin et al. 2011), as well as the detailed small-scale and phase-space structures of individual dark matter halos (e.g., Diemand et al. 2007, 2008; Diemand & Kuhlen 2008; Springel et al. 2008; Stadel et al. 2009; Vogelsberger et al. 2009; Zemp et al. 2009). Currently, the major challenge for computational structure formation simulations is to model the complex physics of star formation and feedback in a robust fashion.

In order to describe the properties and distribution of dark matter halos in such simulations, it is common practice to use concepts that are rooted in the simplistic analytic models from the early days of the era of cosmological structure formation. Concretely, inspired by the spherical collapse model, the sizes and masses of dark matter halos are defined with respect to a density threshold that in general evolves with cosmic time. Such definitions lead to pseudo-evolution: an artificial growth of the size and mass of a dark matter halo without much happening within fixed physical scales (Diemand et al. 2007; Cuesta et al. 2008; Diemer et al. 2013b, see also Appendix A for an illustration of the pseudo-evolution mechanism). In this paper, we expand on earlier studies and shed further light on the issue from various previously unexplored angles. We also discuss potential solutions and more physical descriptions.

The structure of this paper is as follows. In Section 2, we review the ideas and models that have inspired dark matter halo size/mass definitions based on spherical overdensities. In Section 3, we study the evolution of various properties of dark matter halos in a cosmological structure formation simulation over the last 10 Gyr. This is then followed by a detailed discussion in Section 4, where we discuss the implications on models, descriptions, and relations within the realm of structure and galaxy formation. We then conclude and summarize in Section 5.

The simple picture of the formation of halos in the universe within the framework of the spherical collapse model served as an inspiration for definitions of sizes and masses of halos commonly used. Classically (see, for example, Binney & Tremaine 2008, Section 9.2.1, or Knobel 2012; Section 2.3.1), this is derived for an Einstein-de Sitter universe (i.e., Ω_M = 1, Ω_Λ = 0, Ω_K = 0, ρ_mean = ρ_crit at all times) where a spherical region is overdense with respect to the rest of the universe. The matter is assumed to be cold (i.e., no shell crossing) and to have no angular momentum (i.e., only radial motions). Therefore, the dynamics of the overdensity is solely determined by the total interior mass. This overdensity will initially expand, but this expansion will slow down due to gravity and will come to a halt at the turnaround radius at a time t_ta. At this turnaround point, the density contrast of the sphere with respect to the mean matter density of the universe will be ρ_sphere(t_ta)/ρ_mean(t_ta) = 9π²/16 ≈ 5.552. Theoretically, the oversimplified spherical collapse model would reach a singularity at t_vir = 2t_ta, but in reality the collapsing dark matter will virialize and settle into an equilibrium configuration which we call a halo. Energy conservation and the virial theorem then imply that the size of the sphere is r_vir = r(t_ta)/2. Since for the mean matter density in an Einstein-de Sitter universe ρ_mean∝t⁻², the mean matter density drops by a factor of four and the density within the sphere increases by a factor of eight. Thus, we get for the density contrast at virialization

$$\begin{equation} \Delta _\mathrm{vir} \equiv \frac{\rho _\mathrm{sphere}(t_\mathrm{vir})}{\rho _\mathrm{mean}(t_\mathrm{vir})} = 18 \pi ^2 \approx 177.7, \end{equation} \tag{ 1 }$$

i.e., the average density within the virialized region of size r_vir is given by

$$\begin{equation} \rho _\mathrm{vir} = 18 \pi ^2 \rho _\mathrm{mean} \end{equation} \tag{ 2 }$$

and the mass is given by

$$\begin{equation} M_\mathrm{vir} = (4 \pi /3) r_\mathrm{vir}^{3} \rho _\mathrm{vir} \,. \end{equation} \tag{ 3 }$$

The classic derivation for the Einstein-de Sitter universe can be generalized for an open universe (i.e., Ω_M < 1, Ω_Λ = 0, Ω_K > 0) or a ΛCDM universe (i.e., Ω_M + Ω_Λ = 1, Ω_K = 0) (Lacey & Cole 1993; Eke et al. 1996, 1998; Bryan & Norman 1998). For the more relevant case of a ΛCDM universe, there are two common formulae that are used. Eke et al. (1998) fit their calculations with

$$\begin{equation} \rho _\mathrm{vir} = 178 \, \Omega _\mathrm{M}^{0.45} \rho _\mathrm{crit}, \end{equation} \tag{ 4 }$$

whereas Bryan & Norman (1998) find the following fit

$$\begin{equation} \rho _\mathrm{vir} = (18 \pi ^2 - 82 q - 39 q^2) \rho _\mathrm{crit}, \end{equation} \tag{ 5 }$$

where

$$\begin{equation} q \equiv \frac{\Omega _{\Lambda,0}}{\Omega _\mathrm{M,0} a^{-3} + \Omega _{\Lambda,0}} \end{equation} \tag{ 6 }$$

and a is the scale factor.

These definitions of size and mass of a halo inspired by the spherical collapse model can be put in a more general form

$$\begin{equation} \rho _\mathrm{SO} \equiv \frac{M_\mathrm{SO}}{(4 \pi /3) r_\mathrm{SO}^{3}} = \Delta \, \rho _\mathrm{ref}, \end{equation} \tag{ 7 }$$

so that the halo is defined by the sphere where the average spherical overdensity (SO) reaches a specified threshold. In general, Δ can be constant or a function of time and the reference density ρ_ref is either chosen to be the mean matter density of the universe ρ_mean (often also denoted as the background density) or the critical density of the universe ρ_crit, which are both functions of time as well.

Often used definitions in practice are ρ_SO = 200 ρ_mean (200m), ρ_SO = 180 ρ_crit (180c) or one of the definitions already given in Equations (4) (E98) or (5) (BN98) above. In Figure 1, we show the evolution of those four commonly used definitions over a range of 10 Gyr in a ΛCDM universe with Ω_{M, 0} = 0.28 and Ω_{Λ, 0} = 0.72. At early times (high redshift), the contribution of a cosmological constant is vanishing, Ω_M → 1, and, as a consequence, ρ_mean → ρ_crit. Therefore, all four spherical overdensity definitions are similar at high redshift. Today at z = 0, there is a difference of a factor of around three between the 200m and the 180c definitions with the E98 and BN98 definitions in the middle.¹ The overall evolution of all four definitions over approximately the last 10 Gyr is quite evident from Figure 1. In the literature, any of these four definitions, and other similar definitions, are used as proxy to define a virialized region. For example, in the case 200m, the formal virial radius and mass are defined by r_vir = r_200m and M_vir = M_200m = M(r_200m), respectively. In Appendix B, we illustrate the evolution of the difference between the various formal virial mass definitions used in this study.

Zoom In Zoom Out Reset image size

Here, it is interesting to remark that not all of those definitions will be physically meaningful in the far future. In the future, the critical density will approach the asymptotic value $\rho _\mathrm{crit,\infty } = 3 H_0^2 \Omega _{\Lambda,0}/8 \pi G \, (= 97.91 \, M_\odot \, \mathrm{kpc}^{-3}$ in our cosmology), whereas the mean matter density ρ_mean drops to zero. Therefore, the 180c and BN98 definitions will reach asymptotic values and are well defined also in the far future. The 200m and E98 spherical overdensity thresholds approach zero and will enclose increasingly larger regions. The enclosed mass in these cases will converge but is typically larger than the ultimate bound mass within the turnaround radius (see, e.g., Busha et al. 2005). Of course, the E98 and BN98 fits to the generalized spherical collapse model are only valid for a certain range of scale factors (for details, consult those papers) and should not be extrapolated into the far future.

3.1. Simulation Data

3.2. Spherical Overdensity Mass and Size Evolution

In Figure 2, we show the mass evolution for the size/mass definition 200m, M_200m, over the last ca. 10 Gyr for different classes of halos as specified in Table 1. For comparison, we also plot the evolution of the mass within a fixed physical scale of r_{200m, 0}, the size each halo had at z = 0. For each halo track the masses are normalized by the mass of the halo at z = 0, M_{200m, 0}. It can be easily seen that the spherical overdensity size/mass definition overpredicts the growth when compared to the mass content within a fixed physical scale. For example, for halos that are in the mass range 10.0 ⩽ log₁₀(M_{200m, 0}/M_☉) < 10.5, one would assume that they have grown by a factor of two in the last 10 Gyr when looking at M_200m, but typically all of the mass was essentially already in place 10 Gyr ago. Similarly, for the mass range 12.0 ⩽ log₁₀(M_{200m, 0}/M_☉) < 12.5, one would traditionally come to the conclusion that those objects have typically grown by ca. a factor of three, whereas ≈75% of the mass was already in place within r_{200m, 0} 10 Gyr ago.

Zoom In Zoom Out Reset image size

This discrepancy is clearly due to the evolving spherical overdensity threshold ρ_SO, as shown in Figure 1. If we fix this threshold density to the value it has at z = 0 and then calculate the mass within the resulting sizes $r_\mathrm{200m_0}$, we can see in Figure 2 that this definition is essentially identical to the evolution of the mass within the fixed physical scale r_{200m, 0}. Only for halos in higher-mass bins do we see some deviation mainly at high redshift.

In Figure 3, a complementary version of Figure 2 for size, we present for a subset of the mass categories (for brevity) the evolution of the different size scales associated with the different definitions. The large growth in mass seen before is also reflected in the extensive increase in the associated size r_200m. We also see that the radius that contains the mass at z = 0, r(M_{200m, 0}) is still contracting for higher-mass halos at high redshift. Later on, however, and for low-mass objects that mass shell is essentially in place and does not vary much with time. Similarly, the scale $r_\mathrm{200m_0}$ resulting from keeping the spherical overdensity threshold constant shows only some growth early on for high-mass objects.

Zoom In Zoom Out Reset image size

This is a clear indication that for all halo classes there was much more matter already in place 10 Gyr ago than one would infer from the usual spherical overdensity definition. For small objects there does not seem to be much physical accretion at all and the amount of physical accretion increases with halo mass. However, for all halo-mass bins the amount of growth is overpredicted when compared to what happens within fixed physical scales. This pseudo-evolution due to the evolving density threshold leads to an illusion of growth for dark matter halos.

In order to present our findings in a compact way and to include the results from the other two size/mass definitions, we show the average relative growth rate of the various mass medians of Figure 2 for all three size/mass definitions as a function of halo class in Figure 4. For every mass bin we calculate for each median curve m in Figure 2

$$\begin{equation} \kappa \equiv \frac{(m_0-m_{1.592})/m_{1.592}}{\Delta t}, \end{equation} \tag{ 8 }$$

with Δt = 9.609 Gyr over the redshift range z = 1.592...0 in our cosmology. We show the same three types of medians as plotted in Figure 2 (same line style): radius fixed at z = 0 value (solid, e.g., M(r_{200m, 0})/M_{200m, 0}), evolving density threshold (dashed, e.g., M_200m/M_{200m, 0}), and density threshold fixed at z = 0 value (dash–dotted, e.g., $M(r_\mathrm{200m_0})/M_\mathrm{200m,0}$). Of course, the growth is, in general, not linear with cosmic time over this time interval (as can be seen in Figure 2), but this average linear growth rate serves as a good indicator to quantify the differences between the various definitions.

Zoom In Zoom Out Reset image size

In general, the average relative growth rate κ increases with halo mass for all mass definitions. There is, however, quite a discrepancy between the evolving spherical overdensity definition (dashed) and the other two definitions where either the physical size or the spherical overdensity threshold was held constant. The difference is larger the more the spherical overdensity threshold (see Figure 1) varies with time, i.e., largest for the 200m and smallest for the 180c definition. Also, for small-mass halos the degree of pseudo-evolution is larger than for high-mass objects.³

3.3. Intrinsic Physical Scales

In this section, we would like to quantify the evolution of several intrinsic physical scales of dark matter halos.

3.3.1. Maximum Value of Circular Velocity Curve

The maximum value of the circular velocity curve is an intrinsic physical scale for an individual dark matter halo. The value of v_{c, max} is often used in the literature as an alternative to mass in order to characterize the "size" of a dark matter halo. In Figure 5, we show the evolution of v_{c, max} over the last 10 Gyr. The same mass bins as in Figure 3 are shown. We see that the overall evolution is rather moderate and better reflects the modest evolution of the size and mass of halos when fixed physical scales are considered (see Figures 2 and 3).

Zoom In Zoom Out Reset image size

As expected for a characteristic scale of a dark matter halo, the determination of the maximum value of the circular velocity curve is not sensitive to the size/mass definition used. Therefore, the evolution of v_{c, max} differs only marginally for the other two size/mass definitions from what is shown in Figure 5 for the 200m definition.

3.3.2. Scale Radius

The family of the αβγ models, given by

$$\begin{equation} \rho (r) = \frac{\rho _0}{(r/r_\mathrm{s})^\gamma [1+(r/r_\mathrm{s})^\alpha ]^\frac{\beta -\gamma }{\alpha }} \,, \end{equation} \tag{ 9 }$$

is often used to describe the spherically averaged density structure of various objects in the universe. Navarro et al. (1996) found that the parameter set (α, β, γ) = (1, 3, 1) provides excellent fits for the density profiles of dark matter halos in dissipationless structure formation simulations: the so-called Navarro–Frenk–White (NFW) profile. Recent work has shown, however, that the central density profile is shallower than ρ∝r⁻¹ (Springel et al. 2008; Stadel et al. 2009; Navarro et al. 2010) and that the Einasto (1965) profile

$$\begin{equation} \rho (r) = \rho _0 \exp \left(-\frac{2}{\alpha }[(r/r_\mathrm{s})^\alpha -1]\right) \end{equation} \tag{ 10 }$$

often provides a better fit.

A generic feature of both of these density profiles is that at a certain physical scale, the scale radius r_s, a transition between a shallower inner region to a more steeper outer regions occurs. The dividing logarithmic slope at r_s is $\mathrm{d} \ln (\rho) / \mathrm{d} \ln (r) |_{r_\mathrm{s}} = -2$ for the NFW as well as the Einasto profile, i.e., r_s = r₋₂, the location where the logarithmic slope is −2. Therefore, we would like to describe the evolution of this physical transition scale r_s with cosmic time.

For each dark matter halo, we fit an NFW profile where the fitting radial range was limited to 0.8 times the size the halo had at z = 0 (i.e., r_{200m, 0}, r_{E98, 0}, and r_{180c, 0}) at all times, i.e., the fitting range was constant in time. It is common practice to limit the fitting range in order to avoid complications with massive substructures in the outer parts of dark matter halos (see, e.g., Navarro et al. 2010). We have experimented with various other fitting ranges and our conclusions are not sensitive to the specific choice. The Einasto profile has one more free parameter than the NFW profile. We fixed the shape parameter to α = 0.16, which is similar to other recent work in the literature (e.g., Springel et al. 2008; Navarro et al. 2010). We have checked that using similar reasonable values for α do not dramatically affect the fit values for r_s. Therefore, our fitting procedure is quite robust and not sensitive to choices of parameters.⁴ The fitting was done with the non-linear least-squares minimization Python package lmfit,⁵ where we used the standard Levenberg–Marquardt algorithm in order to minimize the figure-of-merit or residual function

$$\begin{equation} R^{2} \equiv \sum _{i=1}^{N_\mathrm{bin}} \left(\ln \left(\rho _i^\mathrm{data}\right) - \ln \left(\rho _i^\mathrm{fit}\right)^{2}\right. \,. \end{equation} \tag{ 11 }$$

In Figure 6, we show the relative evolution of the scale radius for the NFW and Einasto profile fits for the 200m size/mass definition. We only show mass bins where the formal virial radius r_{200m, 0} contains at least ca. 1000 particles so that the internal structure of the halo is properly resolved. For the low-mass halos there is nearly no evolution present, whereas the 12.5 ⩽ log₁₀(M_{200m, 0}/M_☉) < 13.0 mass class shows moderate growth over the last 10 Gyr. The determination of the scale radius is not sensitive to the details of the size/mass definition, as expected for an intrinsic scale, and the results for the other two size/mass definitions are essentially the same as in Figure 6. Our results are in good agreement with previous results in the literature which also show that the scale radius is approximately constant in time for dark matter halos after their formation (e.g., Bullock et al. 2001; Zhao et al. 2003).

Zoom In Zoom Out Reset image size

An alternative way of describing the density profile of a dark matter halo is relating the scale radius r_s to the virial radius r_vir via the virial concentration parameter

$$\begin{equation} c_\mathrm{vir} \equiv r_\mathrm{vir}/r_\mathrm{s}, \end{equation} \tag{ 12 }$$

where any spherical overdensity proxy is generally used for the virial radius (see Section 2). In Figure 7, we show for brevity the evolution of the concentration parameter for roughly Milky-Way-sized objects (11.5 ⩽ log₁₀(M_{200m, 0}/M_☉) < 12.0). We show the median evolution for both the NFW as well as the Einasto profile fits. Since the scale radius essentially remains constant (see Figure 6), it is clear that nearly all of the evolution is driven by the artificial growth of the virial radius proxy r_200m for this halo class, and one would think that the concentration had tripled in the last 10 Gyr. In Diemer et al. (2013b), they find as well that the concentration evolution is mainly driven by pseudo-evolution especially for low-mass halos and we discuss this in more detail in Section 4.

Zoom In Zoom Out Reset image size

In Figure 7, we also present the evolution of an alternative concentration measure (Alam et al. 2002; Diemand et al. 2007)

$$\begin{equation} c_\mathrm{v} \equiv \frac{\bar{\rho }(r_{v_\mathrm{c,max}})}{\rho _\mathrm{crit,0}} = 2 \left(\frac{v_\mathrm{c,max}}{H_0 \, r_{v_\mathrm{c,max}}}\right)^{2} \,. \end{equation} \tag{ 13 }$$

This intrinsic concentration measure has the advantage that it is well defined both for isolated halos and subhalos. It does not make any assumptions about a specific shape of the density profile and does not depend on any arbitrarily defined outer edge of the halo. This concentration measure remains rather flat over the same time range though with quite a large scatter. Since the circular velocity curve is typically rather flat over some range, the location where the maximum circular velocity is reached, $r_{v_\mathrm{c,max}}$, can exhibit some scatter over time which is then reflected in the scatter of c_v. The large scatter is certainly a drawback for c_v, but the median behavior nonetheless indicates that a concentration measure based on intrinsic properties of the circular velocity curve of the dark matter halo, c_v, shows a qualitatively different behavior than the virial concentration c_vir.

3.3.3. Virialized Region

In this section, we aim to explore the possibility of estimating the virialized region via an operational definition. A natural choice is to look at the local virial ratio

$$\begin{equation} Q \equiv -\frac{2 K}{W} \end{equation} \tag{ 14 }$$

in each spherical shell of our profiles, where K is the kinetic and W is the potential energy of the matter in the shell (see Appendix C for more information on how we calculate Q in detail). From the scalar viral theorem, 2K + W = 0 follows that for virial equilibrium we get Q = 1 and the shell is unbound for Q > 2. Here, we define the virialized region where Q < Q_crit = 1.4. Of course, this operational definition depends on the choice of the Q_crit value. Values similar to those used here with Q_crit = 1.4 are used in the literature (e.g., Bett et al. 2007; Neto et al. 2007). We have checked that our qualitative conclusions do not change with similar choices for Q_crit.

In Figure 8, we show the evolution of this operationally defined virial radius r_vir with time. We compare with the E98 definition, which should be the appropriate collapse model for a ΛCDM cosmology and show the results for two different halo masses. It is evident that the virialized region in general physically grows with time, i.e., the median r_vir/r_{E98, 0} curve increases. However, the E98 definition is in general underestimating the virialized region.

Zoom In Zoom Out Reset image size

In Figure 9, the median ratio of the resulting virial mass M_vir with the mass from our three size/mass definitions are shown as a function of halo mass at z = 0. We find that for all masses and size/mass definitions, the operationally defined virial mass does not coincide with the spherical overdensity definition. The characteristic shape of this ratio with the initial increase with halo mass and then decrease for higher-mass halos is due to the typical radial outflow (low-mass halos) and infall (high-mass halos) patterns around those dark matter halos. For a more detailed discussion, see Cuesta et al. (2008; Figures 4 and 6) or the brief recapitulation below.

Zoom In Zoom Out Reset image size

Macciò et al. (2003) also used the differential virial ratio in order to define a boundary layer of a halo. They, however, define the potential energy differently than we do in this study (see Appendix C). They also find that virial masses determined via their method are larger than the E98 size/mass definition. Cuesta et al. (2008) used a different operational definition of the virialized region (see also Busha et al. 2005; Prada et al. 2006). They use the mean radial velocity as a function of radius for a class of objects. Typically, in central regions of dark matter halos that are in equilibrium, the radial velocity is close to zero. Further out, the radial velocity curve starts deviating from zero. For objects that are still accreting, the radial velocity is negative over some radial range further out due to the infalling matter. Even further away, cosmic expansion manifests itself and the radial velocity becomes positive and grows with distance. They found that small-mass objects no longer have an infall region at z = 0 and the radial velocity shows an outflow pattern beyond the central static region where the radial velocity is close to zero. At earlier times, there was a moderate inflow region present in between the inner equilibrium and the outer cosmic expansion region for those small-mass objects, which essentially has disappeared by z = 1. For more massive objects, such an infall region is still present today and it becomes more prominent with halo mass. In the far future, though, even for massive clusters, this infall pattern will disappear (Busha et al. 2005). We refer the reader to Figures 4 and 13 in Cuesta et al. (2008), which depict these patterns excellently. They then define the virialized region as the central radial range where the radial velocity remains within a certain threshold from zero. They call the size/mass associated with their definition the static radius r_static and static mass M_static, respectively. Their findings for the static size/mass are in good qualitative agreement with our results.⁶ For example, they find a similar behavior of their static mass as we do in Figure 9, as well as a general discrepancy between M_static and M_BN98, which they use as their standard spherical overdensity definition (see their Figures 6, 9, 14, and 15).

The aim here is certainly not to establish or promote a better way to define the virialized region of a dark matter halo—actually, it is rather the contrary, as we will discuss in Section 4. All operationally defined estimates depend on some choice of parameters as well as assumed simplifications and symmetries (e.g., spherical symmetry). Additionally, individual dark matter property profiles can be quite noisy if one is interested in local (differential) properties. This can lead to some scatter (as in our case) and one can, in general, only make meaningful statements in a statistical sense by averaging over a class of halos than for an individual halo. This is certainly a drawback of such methods when compared to the simple spherical overdensity size/mass definition that operates with a smoother, integrated quantity, namely, the cumulative mass profile. Also, such definitions can only be applied to simulation data where the full phase-space information is available. Ultimately, we must compare our results to observational data where such operational definitions of virialized regions are not feasible.

The main point we would make here is that, in general, a simple spherical overdensity criterion does not imply deeper fundamental physical properties like virialization. A physically more precise nomenclature would be to name such size/mass definitions by their characteristic definition (e.g., r_200m or M_200m) and refrain from the practice of calling them virial radius or virial mass, since this is a deeper inherent physical property that does not follow from simple overdensity-based size/mass definitions.

3.4. Halo Build-up

In this section, we aim to illustrate how a halo of a certain mass class typically builds up over time. In Figure 10, we plot the medians over all halo tracks within a halo class of mass fractions M(x)/M_{200m, 0} within various constant physical radii x ≡ r/r_{200m, 0} = 3, 2.75, 2.5, 2.25, 2, 1.75, 1.5, 1.25, 1, 0.75, 0.5, 0.25, 0.1 (from top to bottom). The x = 0.1 line is missing for the smallest two mass classes since this is below the resolution scale in our simulation box.

Zoom In Zoom Out Reset image size

As evident from Figure 10, the inner region of a halo was typically already in place 10 Gyr ago (see, e.g., the x = 0.1 and x = 0.25 median lines). The region out to r_{200m, 0} does not show much evolution as well for small mass halos. For more massive objects, the mass assembly in the outer region is clearly visible. At early times beyond r_{200m, 0}, the Hubble flow manifests itself since the mass within a fixed physical radius is decreasing. In general, closer mass shells stabilize earlier than more distant mass shells. However, during the last few gigayears the mass shells beyond r_{200m, 0} do not show much evolution for all halo classes except for the most massive objects we consider. Those objects are still in the assembly process today. This is a good illustration that the mass profiles around dark matter density peaks extend way beyond a formal virial radius and do not show much evolution once a specific mass shell stabilized after the decoupling from the cosmic expansion (see also the discussion in Section 4).

This picture of the halo build-up is nicely complemented by the studies of Prada et al. (2006), Cuesta et al. (2008), and Diemer & Kravtsov (2014), where they have studied the radial velocity patterns around dark matter halos (see also Section 3.3.3). Low-mass objects today (i.e., low sigma peaks in the density field) stopped accreting a while ago since their infall region disappeared and the matter distribution about those peaks has essentially stabilized even beyond the formal virial radius. More massive objects like, e.g., clusters (i.e., high sigma peaks) still have an infall region at z = 0 and are therefore physically growing today. Figures 4 and 13 of Cuesta et al. (2008) and Figure 8 in Diemer & Kravtsov (2014) depict these radial velocity patterns excellently.

In the previous Section 3, we have seen that the evolution of spherical overdensity sizes and masses can lead to an inaccurate and distorted picture of the assembly process of dark matter halos. Especially for low-mass objects, most of the physical matter distribution was already in place 10 Gyr ago (see also Appendix A). Aspects of this pseudo-evolution effect have been pointed out before in the literature. In Diemand et al. (2007), they showed that, for a Milky-Way-sized object, the Via Lactea I halo, after the last major merger at z ≈ 1.7, the density profile is essentially established and the mass within fixed physical scales is roughly constant, whereas the formal virial mass M_200m is still growing (see their Figures 3 and 4). They also point out that the small degree of physical accretion is typical for Milky-Way-sized halos by looking at a larger sample of halos in a cosmological volume. Cuesta et al. (2008) also show that the dark matter profile of a simulated Milky-Way-sized object typically does not evolve after z = 1 (see their Figures 2 and 16). However, their main focus is on the discrepancy between an operationally defined virial mass, which they call static mass, that is defined by the size of the region of a halo with zero radial velocity and the BN98 virial mass definition. They found that the BN98 virial mass definition is, in general, underestimating/overestimating the static mass at low/high redshift (see their Figures 15 and 16). The transition mass scale where M_static = M_BN98 decreases dramatically with redshift so that, especially for low-mass objects, the associated dynamic mass of a halo has already been underestimated for a long time. This is similar to our findings in Section 3.3.3 for an operational definition of the virial mass based on the local virial ratio. This effect can be understood by the typical shape of the radial velocity profile for the different halo classes (see also Section 3.3.3). Thus, in general, the static mass does not coincide with the BN98 virial mass proxy. Diemer et al. (2013b) looked at the pseudo-evolution effect in a statistical sense for a wide range of halo masses in the Bolshoi simulation (Klypin et al. 2011). They investigate the differences between the density profiles of dark matter halos between z = 1 and z = 0 and found that low-mass halos only grow by ca. 10% on average during that time span. By assuming a static halo density profile at z = 1/0, they also estimate the pseudo-evolution for each halo in two modes by forward/backward integration over this static profile with time. We are tracking the individual halos for several simulation snapshots in cosmic time and quantify the pseudo-evolution effect differently by looking at fixed physical and intrinsic scales of the dark matter halos. Their findings, however, are in good agreement with our results (see, e.g., their Figure 4). Weinmann et al. (2013) studied the evolution of the circular velocity function for objects in the Millennium-II simulation (Boylan-Kolchin et al. 2009), measured at a fixed physical radius of 20 kpc and found "remarkably little evolution since z = 4" (see their Figure 4). This also indicates that the central region of dark matter halos has already been in place since early times (see also our Figure 10).

The growth of formal virial quantities mainly reflects the underlying cosmic expansion and dilution of matter—though in a very convoluted way. The fundamental problem is more that a structural description of a dark matter halo via a single size/mass is not fully adequate. Within the framework of the spherical collapse model, the concept of virial mass is a perfectly suitable description of the halo. All the matter within the overdense perturbation ends up in a spherical ball of a fixed size and the rest of the universe is expanding away for eternity. Therefore, a single size/mass number describes the properties of a halo in this scenario appropriately. However, when we apply this concept to the much more complex configuration of dark matter halos that form in cosmological structure formation simulations, we get the artefact of pseudo-evolution.

The effect of pseudo-evolution is also present for another popular dark matter halo size/mass definition: the friends-of-friends (FoF) method (Davis et al. 1985). With this method, particles in a simulation are linked together if their mutual distance is smaller than some linking length, which is a free parameter of the method. Since it is common practice to keep this linking length constant in comoving coordinates, i.e., the linking length grows in physical coordinates, the size and mass definition of FoF halos grows accordingly. Therefore, the FoF method is also subject to a similar pseudo-evolution effect as in the case for the spherical overdensity size/mass definitions.

Before it was feasible to extract halo tracks from cosmological structure formation simulations, such merger trees were modeled analytically via the extended Press-Schechter formalism (Bond et al. 1991; Bower 1991; Lacey & Cole 1993). These analytic prescriptions are still used nowadays since merger trees can be generated much faster, with larger statistics and input parameters (e.g., for a different cosmology) can be varied more easily than for merger trees extracted from simulations. However, the underlying physical process is still the spherical (or often a more general ellipsoidal) collapse model. These analytic merger trees are overpredicting the degree of evolution since the mass evolution is very similar to merger trees extracted from simulations where a pseudo-evolving mass definition was used (see, e.g., Jiang & van den Bosch 2014).

In order to describe the density profiles of dark matter halos, we believe that it is physically more intuitive to anchor the parameterization on intrinsic properties. For example, in the case of an NFW profile (Equation (9)), one would specify the scale radius r_s and normalization ρ₀⁷ (see, e.g., Zemp et al. 2008), instead of linking a formal evolving outer edge r_vir to a scale radius r_s, which is rather constant via another evolving parameter, the virial concentration c_vir. Dark matter halos do not have a well defined outer edge and the density profiles are smooth and can be extended well beyond the formal virial radius (Prada et al. 2006; Tavio et al. 2008). Thus, anchoring the parameterization of density profiles on intrinsic physical properties seems to be more natural.

Is it then possible to define an edge or a virialized region of a halo with the full phase-space information available in simulations? The answer is certainly yes. In Section 3.3.3, we have used the local virial ratio to estimate a virial radius. Busha et al. (2005) and Cuesta et al. (2008) use the radial velocity profile to find the equilibrated central region of halos in order to define a size/mass of a halo. Inspecting the radial velocity patterns around dark matter peaks allows one to identify the static as well as the turnaround region of a halo. This is probably the most promising method in order to define an associated dynamic mass of a dark matter halo and gives detailed insight on how the physical mass accretion happens. In practice, this method can only be done in a satisfactory way by averaging over radial profiles of similar halos since individual halo profiles can be quite noisy. This is a minor drawback of this method. Falck et al. (2012) present the novel halo finder ORIGAMI, which tracks the folding of the dark matter sheet in phase-space, which allows for a dynamical definition of a halo. They also find that their halos are larger than corresponding FoF masses. All these phase-space based definitions, however, are much more complex to calculate than specifying the mass within a certain radius. Ultimately, one would like to compare theoretical results from simulations and models with observations, where such definitions are not practical since the necessary phase-space information is not available. Such definitions might be useful, however, to compare theoretical models with simulation results.

Is it essential to define a virialized region or an outer edge of a halo? We do not necessarily think so and it might be helpful to liberate ourselves from the Gedankenkorsett⁸ of the simplistic spherical collapse model. For example, for the dynamics or gravitational lensing it is not relevant whether or not the matter further out is virialized. The matter is there and acts gravitationally. This is also reflected in the fact that, for example, other (infalling) halos feel the influence of their future host way beyond the formal virial radius and can be physically affected (e.g., Hahn et al. 2009; Bahé et al. 2013; Behroozi et al. 2014). Diemand & Kuhlen (2008) have studied the dynamics in the outskirts of a Milky-Way-sized dark matter halo (Via Lactea I) and found that matter (particle) and subhalo orbits typically extend out to about 90% of their turnaround radius after their first pericenter passage, i.e., a typical orbit brings subhalos and dark matter particles back to beyond the formal virial radius of a halo on the first two passages (see also Moore et al. 2004; Gill et al. 2005; Ludlow et al. 2009; Bahé et al. 2013). Also, the process of virialization is much more complex. Diemand et al. (2007) showed that inner mass shells contract by more than the canonical value of a factor of two from the spherical collapse model, whereas outer mass shells contract by a factor smaller than two (see their Figure 1). The dynamically affected and collapsing region around a halo is much larger and therefore also the associated mass, which is based on the original definition of the spherical collapse model. Anderhalden & Diemand (2011) measured this total mass within the turnaround region of a dark matter halo and found that the mass of a halo is significantly higher than conventional virial masses and the mass function is better described by the original Press & Schechter (1974) form.

We see that size and mass is a somewhat ill-defined concept in order to describe the complex structure of dark matter halos. Physical processes (e.g., tidal mass loss or dynamical friction) do not stop or start at a specific outer edge of a halo but their effect strength is rather transitioning smoothly with the strength generally increasing with decreasing distance from the host. Of course, one could argue that a spherical overdensity (or FoF) definition is just a convenient operational definition of halo size/mass, but any such definition is in principle arbitrary and does not indicate a fundamental physical scale. It seems fine to us to use a spherical overdensity size/mass definition at a fixed time or redshift in order to classify or rank the objects as we have done in this study. However, such simple size/mass definitions do not imply deeper fundamental physical properties like being virialized. In order to study evolutionary effects, it is not an adequate method though and just obscures the real underlying physical processes.

An additional complication in this matter arises from baryonic physics. The processes of star formation and feedback can have the potential to affect the dark matter distribution on large scales (e.g., Dubinski 1994; Kazantzidis et al. 2004; Stanek et al. 2009; Abadi et al. 2010; Duffy et al. 2010; Kazantzidis et al. 2010; Pedrosa et al. 2010; Tissera et al. 2010; Cui et al. 2012; Zemp et al. 2012; Bryan et al. 2013). Today, the major difficulty in computational structure formation is to simulate these complex processes in a robust way so that one can make reliable predictions. Despite recent improvements (e.g., Gnedin & Kravtsov 2011; Guedes et al. 2011; Hopkins et al. 2011; Wise et al. 2012; Zemp et al. 2012; Agertz et al. 2013; Hopkins et al. 2013a, 2013b; Renaud et al. 2013; Stinson et al. 2013; Sales et al. 2014), this is still an issue (e.g., Scannapieco et al. 2012) and the improvement and shift from phenomenological to more physical models is the subject of current research in computational galaxy formation. The changes in the overall dark matter distribution (e.g., contraction or expansion, depending on the details of the modeling) also affect the details of the size or mass definition of a halo (e.g., Stanek et al. 2009; Cui et al. 2012). Also, the baryon fraction is in principle a function of radius (see, e.g., Figure 5 in Zemp et al. 2012) and a single number specified at a pseudo-evolving formal virial radius is not an adequate description of the complex distribution of the baryons.

Of course, this illusion of growth undermines to some degree any concept that involves a pseudo-evolving size or mass definition of dark matter halos. For example, the standard way to describe the abundance of dark matter halos in the universe is to use the mass function, i.e., the mean number of objects of a certain mass per unit volume and unit mass. The basic functional form of the mass function is inspired by the original Press & Schechter (1974) form, but there exist various parameterizations in order to account for different cosmologies, mass definitions, and evolutionary trends (e.g., Sheth & Tormen 1999; Jenkins et al. 2001; White 2002; Reed et al. 2003; Springel et al. 2005; Warren et al. 2006; Reed et al. 2007; Tinker et al. 2008; Pillepich et al. 2010; Bhattacharya et al. 2011; More et al. 2011; Murray et al. 2013). The pseudo-evolution of the formal virial mass also affects the classification of halos into host and subhalos, since, due to the artificial growth of the virial radius, a neighboring halo could be classified as a subhalo at a later time or for a different size/mass definition. It has become clear in recent studies that the mass function is not of universal shape when aiming at high precision, i.e., the functional form of the mass function will change as a function of time. Decoding and understanding the real physical evolution of the universe from a mass function that uses a pseudo-evolving mass definition is a difficult task. Another way to describe the abundance is the velocity function, where the maximum circular velocity serves as a "size" indicator for a halo. However, the issue of classification into host or subhalo remains in this case as well, since a formal virial radius is often used for this.

In the literature, formal virial quantities also play a major role in describing relations with various other quantities (e.g., Bullock et al. 2001; Wechsler et al. 2002; Zhao et al. 2003, 2009; Macciò et al. 2007, 2008; Evrard et al. 2008; Chae 2011; Moster et al. 2013) as well as in semi-analytic models of galaxy formation, where such relations are used to encapsulate complex physical processes (e.g., Mo et al. 1998; Cole et al. 2000; Baugh 2006; Croton et al. 2006; Benson & Bower 2010; Guo et al. 2011). Any correlation between two pseudo-evolving quantities is essentially trivial since the underlying common mechanism is the cosmic expansion. Diemer et al. (2013b) illustrate this for the virial concentration–virial mass (c_vir − M_vir) relation where they show that pseudo-evolution almost entirely accounts for the evolution of this relation at the low-mass end. It seems difficult to identify plausible processes where a physical quantity should correlate with quantities defined within a pseudo-evolving formal virial region, which just sweeps up matter that was essentially already in place before. It is potentially more insightful to anchor scaling relations to intrinsic physical quantities that do not show pseudo-evolution. This makes it much easier to recognize the relevant physical processes.

How should one define the "size" of a dark matter halo? There is probably no unique or best definition and this likely depends on the application. We believe the driving force behind this question should be how we compare the results from simulations and models with observations. It is important for such comparisons that one uses the same definitions so that one is not comparing apples and oranges. A simple solution to avoid pseudo-evolution would be to just keep the threshold density fixed in time as we have done and shown as well in Figures 2 and 3. Any such definition would, of course, still be arbitrary and not define a fundamental physical scale of a dark matter halo, but in this way one keeps the "meter stick" constant and one can better disentangle real physical growth from pseudo-evolution. This should work fine in recent times and we have done this in this study for the last 10 Gyr, i.e., up to ca. z ≈ 1.6. However, this method is doomed to fail at early cosmic times, depending on the exact threshold density used. For example, in the case of keeping the threshold density at the value of 200 times the mean matter density today (200m), the mean matter density of the universe has reached that value at the redshift of z = 200^1/3 − 1 = 4.848. Therefore, one cannot use the same definition in order to describe halos at high redshift and today. In practice, however, the time or redshift range of interest is often not that extreme (e.g., limited redshift range in an observational survey) and a definition with a fixed density threshold in time is potentially a good way to avoid the effect of pseudo-evolution. This would also be convenient in a practical sense, since one can use existing software tools with only modest modifications. Also, intrinsic physical scales like the maximum value of the circular velocity v_{c, max} (and $r_{v_\mathrm{c,max}}$), the concentration c_v, or the scale radius r_s, which is physically increasing with halo size/mass, could be used to describe the "size" of a dark matter halo and are potentially also easier to compare directly with observations. The usefulness of v_{c, max} and c_v as halo scales has been pointed out before and they are also better suited to describe subhalos (e.g., Klypin et al. 1999; Diemand et al. 2007; Ascasibar & Gottlöber 2008; Knebe et al. 2011, 2013).

In general, a more holistic approach and look at the overall radial profile of physical quantities and a parameterization based on intrinsic physical scales is potentially a promising way to better describe and understand the structure and evolution dark matter halos. For example, another important (derived) virial scale used in the literature is the virial velocity, $v_\mathrm{vir} \equiv \sqrt{G M_\mathrm{vir}/r_\mathrm{vir}}$. This is just the value of the circular velocity at the formal virial radius. The virial velocity is often used in order to set a characteristic energy scale (e.g., $T_\mathrm{vir} \propto v_\mathrm{vir}^2$) of the halo.

In Figure 11, we see that the formal virial velocity v_200m is also pseudo-evolving when compared to the value of $v_\mathrm{200m_0}$ (fixed background density at z = 0). The pseudo-evolution of v_200m is rather small when compared to the effect for the size or mass of a halo as we have seen in previous sections (see, e.g., Figures 2 and 3). Such a modest evolution of the virial velocity is also seen in Weinmann et al. (2013) for a large sample of dark matter halos from the Millennium-II simulation (see their Figure 4). This is due to the partial cancellation of the pseudo-evolution of size and mass for the circular velocity $v_\mathrm{c}(r) \equiv \sqrt{G M(r)/r}$. Viewed differently, it is a characteristic of the shape of the circular velocity curve for a density profile similar to the NFW shape that exhibits a relatively flat plateau after the peak is reached rather than a fundamental scale at the formal virial radius.

Zoom In Zoom Out Reset image size

Similarly, in Diemer et al. (2013a), those authors investigate the relation between mass and velocity dispersion (M–σ). They show that this relation for clusters not only holds at a formal virial scale but over a wide radial range. Thus, it is the result of tight equilibrium relations between radial profiles of physical quantities, rather than a fundamental property due to the virial theorem. This is also a nice illustration that the apparent fundamental relation at a formal virial scale was obscuring a more general relation between physical quantities that holds over a wider radial range. There are likely many more such cases.

Another derived virial scale is the dynamical time at the formal virial radius, $t_\mathrm{vir} \equiv \sqrt{r_\mathrm{vir}^{3}/G M_\mathrm{vir}} = r_\mathrm{vir}/v_\mathrm{vir}$, which often serves as a characteristic timescale for physical processes (e.g., cooling or feedback) in analytic models. In Figure 12, we compare the evolution of the timescale t_200m with the timescale $t_\mathrm{200m_0}$. We can clearly see that for no physical reason the characteristic time t_200m would become much shorter at early times than a more physically motivated timescale $t_\mathrm{200m_0}$. Thus, the general discrepancy between physical scales and formal virial scales could potentially account for some fundamental problems in (semi-)analytic models like, e.g., low-mass galaxy evolution (Weinmann et al. 2012).

Zoom In Zoom Out Reset image size

The current status quo of how we describe the structure and evolution of dark matter halos in our universe is somewhat unsatisfactory since often the underlying physical processes are obscured by an unfortunate, pseudo-evolving parameterization or description. The fundamental problem here is that simplistic analytic models of structure formation cannot capture the complex processes seen in cosmological structure formation simulations. Using then concepts based on formal virial quantities can give a distorted picture that does not properly reflect the actual physical evolution of dark matter halos.

In this study, we have critically investigated the pseudo-evolution effect of sizes and masses of dark matter halos due to a spherical overdensity definition. We summarize our findings as follows.

• 1.  
  The commonly used method of defining sizes and masses of dark matter halos based on an evolving overdensity threshold shows a large degree of pseudo-evolution when compared to the size or mass evolution within fixed physical scales or the evolution of intrinsic physical properties (see Figures 2, 3, 5, 6, and 10).
• 2.  
  The pseudo-evolution is more pronounced for low-mass objects, which essentially had their density profiles in place 10 Gyr ago even beyond a formal virial radius. However, also for group- and cluster-sized objects, which are still physically accreting at z = 0, the pseudo-evolution is substantial. The degree of pseudo-evolution also depends on how much the threshold density varies (see Figures 1 and 4).
• 3.  
  The FoF size and mass definitions are subject to a similar pseudo-evolution effect as the spherical overdensity definitions.
• 4.  
  Analytic merger trees are overestimating the mass growth since they agree well with merger trees extracted from cosmological structure formation simulations that use a pseudo-evolving mass definition. This indicates that it is difficult to capture the complex mass assembly of dark matter halos with simplistic models.
• 5.  
  The concept of size or edge of a dark matter halo is somewhat ill-defined and one can in practice only specify operationally defined scales. Such "size" definitions are in principle arbitrary and do in general not coincide with the virialized region.
• 6.  
  It is a physically more appropriate nomenclature to name these halo scales by their characteristics, e.g., r_200m/M_200m or r_180c/M_180c, and one should refrain from calling them virial radius or virial mass, since this is a deeper physical property that is not implied by such simple "size" definitions.
• 7.  
  In general, we question the widely accepted view of attributing pseudo-evolving formal virial quantities central importance for describing the structure and evolution of dark matter halos (and sometimes even of galaxies within the halos, as done in semi-analytic models). Concepts and relations based on pseudo-evolving formal virial quantities do not properly represent the actual evolution of dark matter halos and lead to an inaccurate picture of the physical evolution of our universe. In order to gain a better physical picture of our universe, it might be auspicious to parameterize such relations with respect to intrinsic physical scales of dark matter halos or with "size" measures that do not show pseudo-evolution.

It is a great pleasure to thank Jürg Diemand for many discussions during our common time in Zurich and Santa Cruz, where many issues presented in this study were debated over coffee. We also thank Jürg Diemand, Benedikt Diemer, Andrey Kravtsov, and Amber Kilgore for comments on a draft version of this paper. We also thank Yinghe Lu for interesting questions and discussions during her Bachelor thesis in the academic year 2012/2013 at Peking University, where she was working on aspects of this work. M.Z. is supported by a 985 grant from Peking University and the International Young Scientist grant 11250110052 by the National Science Foundation of China. The simulation in this work has been performed on the Flux cluster at the University of Michigan. This research has made use of NASA's Astrophysics Data System (ADS), the arXiv.org preprint server, and the Python plotting library Matplotlib (Hunter 2007).

In Figure 13, we illustrate the pseudo-evolution mechanism for an individual dark matter halo. For this object, the 200m size/mass definition would suggest that the formal virial radius r_200m has grown from 202.0 kpc to 674.6 kpc (ca. a factor 3.3 increase) and the virial mass M_200m from 4.580 × 10¹² M_☉ to 9.793 × 10¹² M_☉ (ca. a factor 2.1 increase) in the last 10 Gyr. However, the density profiles do not show much physical evolution in the same time interval which is also reflected in the only 6% increase in the peak circular velocity value. At redshift z = 1.592, the mass within the formal virial radius at z = 0 was M(r_{200m, 0}) = 9.130 × 10¹² M_☉. This means that 93% of the mass within the physical scale r_{200m, 0} was already present 10 Gyr ago for this dark matter halo.

Zoom In Zoom Out Reset image size

Here, we illustrate the evolution of the mass ratio of the three size/mass definitions used in this work. In Figure 14, we plot M_E98/M_180c and M_200m/M_180c over the last 10 Gyr for three mass categories. In general, the discrepancy between the various definitions becomes larger with time as can be intuitively understood from the evolution of the defining threshold density in Figure 1. The ratios do not only grow with time but become also larger for more massive objects. Therefore, there can be quite a difference in halo mass for the same object, depending on the definition used (see also White 2001; Hu & Kravtsov 2003). The exact definition of the size of a halo also influences the classification of objects as host or subhalos, which in turn affects the abundance (i.e., mass function) of the dark matter halos.

Zoom In Zoom Out Reset image size

In each spherical shell, we calculate the velocities of the dark matter particles with respect to the halo center. We convert them to physical velocities by adding the coordinate expansion (Hubble flow) since the simulation was performed in comoving coordinates. The kinetic term in the shell is then simply calculated by

$$\begin{eqnarray} &&K \equiv \frac{1}{2} \sum _{i} m_i \left(v_\mathrm{rad,i}^{2} + v_\mathrm{\varphi,i}^{2} + v_\mathrm{\theta,i}^{2} \right) \end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray} &&= \frac{1}{2} M \left(\bar{v}_\mathrm{rad}^{2} + \bar{v}_\mathrm{\varphi }^{2} + \bar{v}_\mathrm{\theta }^{2} + \sigma _\mathrm{rad}^{2} + \sigma _\mathrm{\varphi }^{2} + \sigma _\mathrm{\theta }^{2} \right), \end{eqnarray} \tag{ C2 }$$

where M = ∑_im_i is the total mass in the shell and the sum is over all particles in the shell. We have decomposed the velocity field in spherical coordinates. The second form allows us to use the mean velocities and dispersions for the different velocity components calculated by the profiling tool.

For an isolated, spherically symmetric structure, the potential at r is given by (Binney & Tremaine 2008, Equation (2.28))

$$\begin{equation} \Phi (r) = -\frac{G M(r)}{r} - G \int _{r}^{\infty } \frac{\mathrm{d} M(r^\prime)}{r^\prime } \,. \end{equation} \tag{ C3 }$$

In cosmology, the potential results from the density perturbations and we can remove the contribution of the constant background density (Jeans Swindle, see, e.g., Binney & Tremaine 2008; Falco et al. 2013). In practice, we only do the integration of the outer part of the potential out to a finite radius instead of infinity or until we reach the mean matter density. We have checked and found it not to be very sensitive on the exact value as long as it is of the order of a few megaparsecs, as in our case. We then interpolate on the resulting potential to get the value at the middle radius of each shell. The potential energy of the mass M in each shell with middle radius r_m is then simply given by

$$\begin{equation} W = \frac{1}{2} M \Phi (r_\mathrm{m}) \,. \end{equation} \tag{ C4 }$$

Here, a side remark is in place. Often one finds studies that use the local (differential), as in this study, or cumulative version of the virial ratio Q as a function of radius, where it shows an upturn towards the center of the halo with large values for Q that would mean the matter in the center is unbound, i.e., Q > 2 (e.g., Cole & Lacey 1996; Macciò et al. 2003; Shaw et al. 2006). Of course, this does not make physical sense and is an artifact of neglecting the contribution of the outer part of the potential (second term on the right-hand side in Equation (C3)) as we have verified. In Shaw et al. (2006), the authors remedy this situation by including a surface pressure term.