A natural boundary of dark matter haloes revealed around the minimum bias and maximum infall locations

In our current understanding of structure formation, dark matter haloes are the building blocks of the large-scale structure in the Universe. The properties and evolution of dark matter haloes are fundamental to models that describe many aspects of the Universe ranging from the galaxy formation and evolution to the overall make-up and history of the Universe. In this halo model framework, the large-scale structure can be decomposed into the distribution of haloes on large scales, convolved by the internal structure of haloes on small scales (e.g., Scoccimarro et al., 2001; Cooray & Sheth, 2002). Galaxies form and evolve in the potential well provided by dark matter haloes, with many intrinsic galaxy properties such as the colour, morphology and mass largely determined by the structure and evolution of the haloes they reside in (e.g. Baugh, 2006; Benson, 2010; Somerville & Davé, 2015).

Despite the substantial work done on understanding the Universe using haloes as building blocks, practical characterisation of the fundamental properties of dark matter haloes is still subject to improvement. Most importantly, our understanding of what is the size of a halo is at best at a premature stage. Almost all of the studies about dark matter haloes so far are based on the classical virial definition of halo size derived from the spherical collapse model (Gunn & Gott, 1972), which is expected to be only an approximate description of these objects. In the spherical collapse model, haloes are modelled from the collapse of a spherical initial overdensity embedded in an otherwise uniform background Universe, and the final size of the halo is defined through virialization argument. However, in the real Universe haloes are live objects that do not necessarily have a clear separation from the background Universe and are constantly accreting from and interacting with the neighbouring non-uniform environment. Furthermore, the top-hat spherical collapse model assumes that collapsing overdensities are isolated, spherical, and with no initial velocity dispersion until relaxation. This leads to the model failing to predict the evolution of simulated haloes after their turn-around epochs (Suto et al., 2016). Such complications make the virial radius more of a practical working definition, with many variants of it defined through different overdensity criteria, such as 200 times the mean density or 200 times the critical density. Despite their popularity, such definitions do not always correctly separate the virialized part of haloes (e.g. Zemp, 2014) and can lead to apparent evolution of halo properties in absence of physical evolution (Diemer et al., 2013).

The recently proposed splashback radius marks a significant improvement over the classical model by incorporating halo accretion into the definition of a dynamical halo boundary (Diemer & Kravtsov, 2014; Adhikari et al., 2014; More et al., 2015; Shi, 2016) and hence has attracted substantial attentions (Snaith et al., 2017; Mansfield et al., 2017; Umetsu & Diemer, 2017; Baxter et al., 2017; Fong et al., 2018; Chang et al., 2018; Xhakaj et al., 2020; Sugiura et al., 2020; Aung et al., 2020; Murata et al., 2020; Deason et al., 2020b, a). Another classical boundary of physical importance is the turnaround radius. The turnaround radius also arises from the spherical collapse picture, but has been relatively less studied. We summarize these boundaries below:

• •

  The virial radius, $r_{\rm vir}$, is the expected radius of a virialized halo according to the spherical collapse model. Normally this is defined through the expected virialization density, which we take from the prediction of (Bryan & Norman, 1998) assuming a tophat initial density in a spherical collapse model.

• •

  The splashback radius, $r_{\rm sp}$, is practically determined to be where the density profile reaches its steepest slope. The steepening in the slope has been attributed to the build-up of particles during their first orbital apogees, where the particles have low radial velocities (Fillmore & Goldreich, 1984; Bertschinger, 1985). This radius has gathered significant interest since its discovery as it has been shown to probe the mass accretion rate of haloes (Diemer & Kravtsov, 2014; Adhikari et al., 2014; More et al., 2015; Shi, 2016). Note the splashback radii of individual halo particles can span a wide range, and hence there is a large freedom in defining the overall splashback radius of a halo from particle dynamics (Mansfield et al., 2017; Diemer, 2017; Diemer et al., 2017). In this work, we will refer to the splashback radius as the one estimated from the steepest slope location, unless explicitly specified otherwise.

• •

  The turnaround radius, $r_{\rm ta}$, in the spherical collapse model is located where a mass shell of a halo reaches zero radial velocity before collapsing back towards the halo at a given time (Mo et al., 2010; Pavlidou & Tomaras, 2014; Tanoglidis et al., 2015, 2016). For individual haloes the radial velocity profiles of particles are a combination of the peculiar and Hubble flow velocities, and the turnaround radius is located where the Hubble flow overcomes the peculiar velocity. This can potentially be used as a cosmological probe, as it reflects the competition between dark energy and gravity (Taruya & Soda, 2000; Falco et al., 2014; Faraoni et al., 2015; Lee & Li, 2017; Korkidis et al., 2019).

Typically $r_{\rm vir}\lesssim r_{\rm sp}<r_{\rm ta}$. Most studies of haloes use mass or radius definitions that are within or roughly around $r_{\rm vir}$. However, the influence of haloes can go well beyond the virial radius. For example massive clusters are found to impact on the star formation rates of galaxies that are well outside their host halo’s virial radius, after the galaxies have passed within the host halo’s virial radius (Wetzel et al., 2014; von der Linden et al., 2010; Wetzel et al., 2012; Kukstas et al., 2020; Adhikari et al., 2020). In Bahé et al. (2013); Wetzel et al. (2014) the authors show that a significant fraction of these galaxies are not infalling for the first time and remain bound to their host haloes (see also Ludlow et al., 2009). Even though one can always work with a given definition of halo size, it is not difficult to imagine that simple intrinsic physical relations may become complicated or obscured in absence of a correct physical description of haloes.

The different halo radii definitions also reflect our understanding of the halo structure in different dimensions. While the virial radius is defined through the average density expected from virialization, the splashback radius is defined through the slope of the density profile or more physically the population property of evolving orbits in an accreting halo (Diemer, 2017; Diemer et al., 2017). The turnaround radius, on the other hand, can be practically defined through the velocity profile around the halo.

In this work, we introduce two other dimensions to define a natural boundary of haloes in aim of providing a physical characterisation of haloes. We do this by first studying the bias profile, that is, the overdensity profile around a halo relative to the average overdensity profile around a random matter particle. We make use of a high resolution $N$-body simulation to extract bias profiles for haloes of different properties. The bias profile shows a typical trough around the halo boundary. This identifies the scale where the correlation between matter and haloes is the weakest, relative to the average clustering of matter, thus providing a natural dimension to define the boundary of a halo. Importantly, this bias trough is found to be close to the region of maximum mass inflow rate around the halo, signalling ongoing depletion around the bottom of the trough. These two characteristics combined lead to the interpretation of the bias trough as the location where matter is being depleted from the environment by the growing haloes over time. Correspondingly, we propose two characterisations of the trough scale: an inner depletion radius defined at the location of the maximum mass inflow rate, and a characteristic depletion radius at the minimum of the bias.

As a first paper in studying the depletion boundary, in this work we aim at establishing the concept and general characterisations of the depletion region, with some additional efforts on characterising the radius at the bias minimum. To this end we examine how the boundary features relate to the virial, splashback, and turnaround radii. By splitting haloes into bins of different halo properties, we also study how the characteristic depletion radius depends on multiple halo properties, and show that it has some simple properties such as a nearly constant enclosed density. In terms of particle orbits, we demonstrate that the inner depletion radius can be interpreted as the outermost splashback radius visually identifiable in the phase space structure around haloes. An immediate application of the inner depletion radius is to improve over the halo model description of the matter clustering around the quasi-linear scale (e.g. Tinker et al., 2005; Hayashi & White, 2008; van den Bosch et al., 2013), which we briefly discuss in this work. In a follow-up work, we will study the depletion process in the spatial and velocity distributions explicitly by tracking haloes over time.

This paper is organized as follows. In Section 2 we describe the halo sample used in this work. In Section 3 we introduce the characteristic depletion radius through the bias profile, examine its dependence on halo properties, and relate it with the splashback radius. In Section 4 we study the dynamical interpretation of the depletion trough and introduce the inner depletion radius through the velocity profile and the phase space structure around haloes, while relating the two radii with the turnaround radius. In Section 5 we compare these radii and their enclosed densities to other halo boundary characterisations. In Section 6 we briefly discuss the implications of the halo depletion radii in the halo model. In Section 7 we summarize and give our conclusions.

Note that $\log$ used in this work is $\log_{10}$; unless otherwise stated all $M_{X}$ and $r_{X}$ units are in ${\rm M_{\odot}}/h$ and ${\rm Mpc}/h$, respectively; the haloes in our work are located at $z=0$.

In this section we introduce the simulation and the halo sample used in this paper. Our halo and clustering data are extracted from one of the CosmicGrowth Simulations (Jing, 2019), a grid of high resolution $N$-body simulations run in different cosmologies using a P³M code (Jing & Suto, 2002). We use the $\Lambda$CDM simulation with cosmological parameters $\Omega_{\rm m}=0.268$, $\Omega_{\Lambda}=0.732$, and $\sigma_{8}=0.831$. The box size is $600~{}{\rm Mpc}/h$ with $3072^{3}$ dark matter particles and softening length $\eta=0.01~{}{\rm Mpc}/h$. Groups are identified with the Friends-of-Friends (FoF) algorithm with a linking length $0.2$ times the mean particle separation. The haloes are then processed with HBT+¹¹1https://github.com/Kambrian/HBTplus (Han et al., 2012; Han et al., 2018) to obtain subhaloes and their evolution histories. The resulting halo catalog has about $2\times 10^{6}$ distinct haloes with masses $10^{11.5}<M_{\rm vir}[{\rm M_{\odot}}/h]<3\times 10^{15}$, where the minimum mass corresponds to roughly $500$ particles within the virial radius to keep a reliable resolution on the structure of haloes.

The halo sample and the list of halo properties are the same as those used in Han et al. (2019), which we briefly describe below. The halo centres are located at the most bound particle of the central subhalo and all haloes are at $z=0$.

• •

  $M_{\rm vir}$: The virial mass of the halo, where $M_{\rm vir}$ is the mass within a spherical volume of radius $r_{\rm vir}$ that encloses the mean density $\Delta_{\rm c}$ times the critical density of the Universe, or $M_{\rm vir}=4\pi r_{\rm vir}^{3}\Delta_{\rm c}\rho_{\rm crit}/3$. The virial overdensity in units of the critical density, $\Delta_{\rm c}$, is predicted from the spherical collapse model (Bryan & Norman, 1998). Note the total mass distribution is used in the computation of the virial mass and radius, not just the bound mass.

• •

  $V_{\rm max}$: The maximum of the circular velocity function, $V_{\rm circ}=\sqrt{GM(<r)/r}$, of the central subhalo. Because $V_{\rm max}$ is tightly correlated with halo mass, for this work we will use $V_{\rm max}/V_{\rm vir}$ to factor out the mass dependency, where $V_{\rm vir}=\sqrt{GM_{\rm vir}/r_{\rm vir}}$, and $G$ is the gravitational constant. This has also been used as a proxy for halo concentration (e.g. Gao & White, 2007; Angulo et al., 2008; Sunayama et al., 2016), and is also a description of the shape of the density profile, where $V_{\rm max}/V_{\rm vir}=1$ correspond to isothermal haloes with flat rotation curves while $V_{\rm max}/V_{\rm vir}>1$ correspond to haloes with steeper outer profiles.

• •

  $j$: The spin of the central subhalo: $j=\frac{L\sqrt{|E|}}{GM^{5/2}}$ (Peebles, 1969), where $L$, $E$, and $M$ are the total angular momentum, energy, and mass of the central subhalo.

• •

  $e$: The shape parameter of the halo, $e=e_{1}$ in this work, where $e_{i}=\lambda_{i}/\sum_{j=1}^{3}\lambda_{j}$ and $e_{1}\geq e_{2}\geq e_{3}$. The subscripts $i=1,2,3$ correspond to the three eigenvalues, $\lambda_{i}$, of the inertial tensor, the square of the three principle axes lengths of the halo mass distribution. The inertial tensor is defined as $I_{w,ij}=\sum_{p}m_{p}x_{p,i}x_{p,j}/r_{p}^{2}$, where $m_{p}$ is the mass of particle $p$, and $\vec{x_{p}}$ and $r_{p}$ are the coordinate and distance to particle $p$ relative to the halo centre.

• •

  $a_{1/2}$: The scale factor of the Universe when the halo mass was half of its final mass $M_{0}$ at $z=0$. The masses used here are calculated using the bound particles to avoid complications due to ejected or fly-by haloes.

• •

  $\delta_{\rm e}$: The halo environment defined as the matter overdensity at a halo-centric distance $r_{\rm e}\approx 1-2{\rm Mpc}/h$ around each halo, $\delta_{\rm e}=\rho(r_{\rm e})/\rho_{\rm m}-1$ (see Han et al., 2019, for more details).

The most straightforward way to study the boundary of a halo is perhaps to examine the density profile around haloes. However, as the average density profile decreases around a halo even out to the largest scale, it is not obvious how a boundary could be defined from the density profile alone. So instead of studying the density profile itself, in this study we work with the normalized density profile of a halo, or more specifically the bias profile.

The bias profile of haloes is defined as the ratio between the halo-matter correlation function and the matter-matter correlation function. In the limit of a single halo, the halo-matter correlation function reduces to the overdensity profile of matter around the halo center. A bias profile can then be similarly defined taking the ratio between the overdensity profile and the matter-matter correlation function. Because the matter-matter correlation function is the average overdensity profile around a random matter particle, the resulting bias profile is a measure of the relative clustering around a halo compared to the clustering around a random matter particle. According to linear theory, the clustering around haloes follows the average matter clustering on large scales up to a constant bias factor. The bias profile thus factors out the background clustering profile of the Universe and leaves only the profile relevant for the haloes.

In this work we first extract individual bias profiles around each halo from the overdensity profile and the matter-matter correlation function. In principle we can proceed to analyze each bias profile individually. However, to suppress noises associated with individual profiles, we will instead focus on analyzing the stacked profiles of haloes binned in various halo properties,

where $\xi_{\rm hm}$ and $\xi_{\rm mm}$ are the halo-matter and matter-matter correlation functions, $\delta(r)=\rho(r)/\rho_{\rm m}-1$ is the overdensity profile of matter around each halo, $\rho_{\rm m}$ is the mean matter density of the Universe, and the averaging is over all the haloes in each halo property bin.

In Figure 1 we show the bias binned by the halo parameters explored in this paper, where each panel is binned by the halo parameter labelled on the top right of each panel. Though the forms of the bias can be complex, there are some common features. On large scales the bias flattens to a constant value, where we expect that the halo-matter cross-correlation, or the ensemble average of the overdensities around haloes, follows the shape of the matter-matter correlation to produce the so called linear bias. In contrast to the decreasing density profile, the constant bias profile on large scales thus provides a natural reference of the background density distribution around haloes. On small scales the bias is large and dominated by the density profile of the central halo. The intermediate scale marks the transition between the one-halo feature and the large-scale background, providing the natural scale to look for the boundary of haloes. Indeed, in nearly all cases there is a clear trough in the bias in the intermediate scale, representing the location where the correlation between haloes and matter is the weakest, with respect to the correlation of matter around a random particle in the density field. The location of this bias minimum thus defines a characteristic radius of the halo, which we will call the characteristic depletion radius, $r_{\rm cd}$, for reasons that will become clear when discussing dynamics in Section 4.

Beyond the bias minimum, the shape of the bias profile closely resembles the density profile of a void (Hamaus et al., 2014), which rises with radius until it reaches a maximum in many cases. One general interpretation for this is that the mass accretion of the central halo depletes matter around it, creating a characteristic trough. Beyond the trough the influences of neighbouring haloes and the background expansion become increasingly important, and the matter distribution can also be dominated by matter in neighbouring haloes, creating the wall region that resembles the wall around voids for the same reason. In this sense, the trough can be regarded as a “relative void" formed by the accretion of the central halo, alongside the competing mass accretion from neighbouring haloes. We will come back to this discussion in later sections.

In the language of a halo model, the profile within the depletion boundary is thus the one-halo term, while that outside the boundary is dominated by the two-halo term. The depletion trough is thus a direct manifestation of halo exclusion. This is also supported by the location and depth of the bias trough in different bins. Taking the mass dependence as an example, the highest mass haloes can only be surrounded by lower mass neighbours. Because it is easy to tightly pack small haloes around big ones, the boundaries of the largest haloes are thus smoothly connected to the neighbourhood, with no obvious troughs created by halo exclusion. By contrast, the smallest haloes typically only find neighbours relatively farther away on halo virial scales, reflecting the difficulty in packing larger haloes around them due to halo exclusion. As we will show later in section 6, the inner depletion radius that characterises the inner edge of ongoing depletion well matches the optimal halo exclusion radius in a halo model.

In order to accurately estimate the characteristic depletion radius defined in the bias profile, we first fit each bias profile with the following function,

with $r_{0}<r_{1}<r_{2}$. This function has four components describing the inner-most one-halo profile ($b\propto r^{-\alpha}$ for $r\ll r_{0}$) before the trough, the rise beyond the trough ($b\sim r^{\beta}$ for $r_{0}<r<r_{1}$), the decrease after the wall ($b\sim r^{-\gamma}$ for $r>r_{1}$) and the large-scale linear bias ($b\sim b_{0}$ for $r\gg r_{2}$). Note the inner-most one-halo profile is better described by a NFW-like (Navarro et al., 1997) double power-law profile, but as we focus on studying the intermediate to large-scale feature, we will not go to that complexity and only adopt an asymptotic single power-law for the one-halo term. We find this fitting function is universal for haloes binned in different or even multiple properties, with the exception of the lowest environmental overdensity bins.

With this parametrization of the bias profile, in the following sections we move on with more quantitative analysis of the depletion radius, relating it to various halo properties as well as to other halo radii.

To get more physical insights into the meaning of the depletion radius we study whether $r_{\rm cd}$ corresponds to any features in the dynamics of matter near the halo boundary. As the halo mass is the primary driver (except $\delta_{\rm e}$) for the characteristic depletion radius, we will focus on binning by halo mass in this section.

In Figure 11 we compare the various characteristic radii as functions of halo virial mass. As discussed above, the inner depletion radius is close to the characteristic depletion radius, and both are close to the turnaround for haloes below $M_{\rm vir}\sim 10^{\rm 13}{\rm M}_{\odot}/h$. The splashback radius, on the other hand, is typically close to but larger than the virial radius.

We do not include the lowest mass bin for the maximum inflow and turnaround locations here, due to the lack of an infall region in the corresponding $v_{\rm r}$ profile (see the discussion in Figure 7). We also emphasize that $r_{\rm cd}$ for the high mass bins may not be reliably estimated due to their weak or absent trough feature in the bias profile, and some more careful studies are needed to accurately measure the behavior in these mass ranges. We show the ratio of many measured radii to the characteristic depletion radius in the bottom panel of Figure 11. As we discussed before, the maximum mass inflow rate location is almost identical with the maximum infall velocity location, both of which are lower than $r_{\rm cd}$ by $10-20\%$ in the mass range studied.

The characteristic depletion radius has a constant ratio to the virial radius when binned by halo mass, with the fitted line of $r_{\rm cd}=2.5r_{\rm vir}$ in Figure 11 (or $r_{\rm cd}\approx 2.0r_{\rm 200m}$, for reference). However, the relation becomes more interesting when further binning by other halo parameters. This can be seen in Figure 12, where the ratios $r_{\rm cd}/r_{\rm vir}$ vary with halo properties other than mass when binned by combinations of internal halo properties. An outlier is again the halo spin, which barely affect the ratio. The trends in $r_{\rm cd}/r_{\rm vir}$ are significantly different from those in Figure 6, showing that the various halo radii can probe different properties of haloes. The dependence when $\delta_{\rm e}$ is involved becomes much different from other panels, but roughly consistent with the behavior in Figure 6, reflecting that the $\delta_{\rm e}$ mostly influences $r_{\rm cd}$ rather than $r_{\rm sp}$ or $r_{\rm vir}$.

Comparing this figure to Figure 6, we find that low $r_{\rm cd}/r_{\rm vir}$ ratios for highly concentrated or early forming haloes correspond to high ratios in $r_{\rm cd}/r_{\rm sp}$, or high ratios in $r_{\rm vir}/r_{\rm sp}$. These correspond to the older, low-mass, highly concentrated, low spin, and spherical haloes (see Figure 4). As discussed previously, the low-mass haloes are more likely to have completed their accretion phases, and part of the low-mass population can be made up of distinct post-pericentre splashback haloes, or “ejected subhaloes”. The ejected subhaloes are thought to be one of the causes of assembly bias for low-mass haloes (Dalal et al., 2008; Sunayama et al., 2016; Mansfield & Kravtsov, 2020), where the tidal gravitational forces from their massive neighbours can have significant impacts on their internal properties (Mansfield & Kravtsov, 2020; Tucci et al., 2020). Seeing how ejected subhaloes can have such incredible impacts on halo properties and bias, it would be interesting to trace the evolution of the sample of haloes in the context of Figure 4, which we leave to a future work.

In Figure 13 we plot the average density enclosed within each characteristic radii, for haloes of different virial masses. The average density within $r_{\rm cd}$ is consistent with a constant of $\rho(<r_{\rm cd})\approx 40.2\rho_{\rm m}$, independent of halo mass. In other words, this depletion radius can be found as a simple spherical overdensity radius containing 40.2 times the background density (or equivalently $\sim 10.8$ times the critical density) of the Universe in our simulation, when only the mass dependence is considered. Similarly, $r_{\rm id}$ can also be found enclosing an almost constant density of $\rho(<r_{\rm id})\approx 62.8\rho_{\rm m}$ over the mass range probed. By contrast, the splashback radius corresponds to varying overdensities and can not be simply defined through an characteristic overdensity criterion, as also found previously (Diemer, 2017). This is mostly due to the sharp increase in density at the low mass end. Excluding the lowest mass bin leads to an average enclosed density of $\rho(<r_{\rm sp})\approx 196\rho_{\rm m}$. For completeness, the average density enclosed in the turnaround radius is also shown. It increases with decreasing mass and becomes close to the depletion density at the low mass end. The higher enclosed density at the low mass end can be understood as their turnaround radii froze at higher redshifts where the density of the Universe was higher.

In Figure 14 we plot the spherical density, $\Delta=\rho(<r_{\rm cd})/\rho_{\rm m}$, for the halo depletion radius binned by two halo parameters, one being mass. Similar to the behavior of $r_{\rm sp}/r_{\rm vir}$, when we bin haloes by two halo parameters, we find that $\Delta$ is not universal but depends on halo properties other than mass. At a fixed mass, the density variation is largely consistent with the $r_{\rm cd}$ variation in the parameter space, where a smaller $r_{\rm cd}$ corresponds to a higher density contrast. The oldest, spherical, and the most concentrated haloes typically have the highest enclosed density within $r_{\rm cd}$.

Note that combining the radius ratio in Figure 12 and the density contrast in Figure 14, one can immediately estimate the ratio between the mass enclosed within $r_{\rm cd}$ and the virial mass. For example, taking the typical $r_{\rm cd}=2.5r_{\rm vir}$ and $\rho(<r_{\rm cd})=40.2\rho_{\rm m}$, the enclosed mass can be found to be $M(<r_{\rm cd})\approx 1.7M_{\rm vir}$.

We have discussed that the depletion trough can be understood as a manifestation of halo exclusion, such that the radii characterising the trough are expected to find direct applications in the halo modelling of the large-scale structure. A better definition and characterisation of the halo is expected to improve both in the simplicity and the accuracy of the resulting halo model. During the preparation of this work, a very recent outstanding work by (Garcia et al., 2020, hereafter G20) has just addressed this problem from an inverse approach compared to ours, by solving for an optimal halo radius definition while optimizing the halo model fitting to the halo matter correlation function. We thus take the leisure to just compare our results to their inferred optimal halo exclusion radius, to demonstrate the significance of our new radius for halo modelling.

In G20, the halo density profile and the optimal halo radius are parametrized through scaling relations with the optimal halo mass. The optimal mass and radius are expected to be different from the conventionally defined ones, so the parameters determining the scaling relations are allowed to vary freely and are solved for by matching the predicted halo-matter correlation to the measurement from simulations. The resulting best-fit scaling relation then provides the basis for defining the optimal mass and radius, by self-consistently interpreting the halo mass as the mass enclosed within the optimal halo radius.

G20 found that their optimal halo radius has a roughly constant ratio to the splashback radius, $r_{\rm G20}\approx 1.3r_{\rm sp87}$, where $r_{\rm sp87}$ is the radius containing 87% of the splashback apogees of all the halo particles according to Diemer et al. (2017). To compare against the G20 result, we will use $1.3r_{\rm sp87}$ as a proxy of $r_{\rm G20}$, where $r_{\rm sp87}$ can be computed from scaling relation with the virial mass. As shown in Figure 11, for the mass range covered by the G20 model, the G20 radius well matches the inner depletion radius, $r_{\rm id}$. Note that we have not adopted any tuning in defining our inner depletion radius. Rather our definition comes from complete physical intuition. The excellent agreement between our ad-hoc characterisation of the halo radius and the optimized exclusion radius required by the halo model is thus truly encouraging.

This agreement is also consistent with the phase space interpretation that $r_{\rm id}$ marks the outermost radius of the splashback particles, or the outer edge of particles orbiting in the halo. The good match indicates this high completeness splashback radius is approximately proportional to $r_{\rm sp87}$. In fact, an apparently better proportionality between $r_{\rm id}$ and the steepest slope radius, $r_{\rm sp}$ can be found as $r_{\rm id}\simeq 1.6r_{\rm sp}$ as seen in Figure 11.

It was shown in G20 that $r_{\rm G20}$ is located at the minima of $r^{2}\xi_{\rm hm}$ for haloes around $10^{13}\mathrm{M}_{\odot}/h$, which is the "by-eye" boundary of the halo. However, their $r_{\rm G20}$ starts to deviate from the "by-eye" boundary in higher mass haloes. As shown in the bottom panel of Figure 11, our $r_{\rm cd}$ defined at the minimum bias location is also located very near the minima of $r^{2}\xi_{\rm hm}$ for all mass bins. This is because the matter auto-correlation, $\xi_{\rm mm}$, is found to have a logarithmic slope close to $-2$ around these scales. As a result, the minimum bias location is also expected where $\xi_{\rm hm}$ has the same slope of $-2$, or where $r^{2}\xi_{\rm hm}$ reaches an extreme. As $r_{\rm cd}$ tends to be further away from $r_{\rm id}$ in high mass haloes, this explains the poorer agreement between the optimal exclusion radius and the "by-eye" boundary in more massive haloes. However, instead of using the $r^{2}\xi_{\rm hm}$ to define the “by-eye" halo boundary, we believe that our definition of $r_{\rm cd}$ using the bias profile is more physically meaningful as discussed in section 3.

In this paper we have proposed two new characterisations of the halo boundary using a large sample of haloes from a cosmological $N$-body simulation, using the spatial and velocity distribution of matter around haloes. In the spatial distribution, the “characteristic depletion radius", $r_{\rm cd}$, is most evident from the existence of a ubiquitous trough in the bias profile around this halo boundary, where the clustering around the halo is the weakest relative to the clustering around a random matter particle. In the velocity domain, an “inner depletion radius", $r_{\rm id}$, can be defined at the location of the maximum mass inflow rate. According to continuity, the maximum inflow radius marks the transition between an inner growing halo and the depleting environment, or the inner edge of the active depletion region.

These two radii are closely related, with $r_{\rm cd}$ slightly larger than $r_{\rm id}$ by $10-20\%$ over the mass range we studied. Matter is drained from the environment outside $r_{\rm id}$, leading to the formation of an accretion trough in the bias profile around $r_{\rm cd}$, which also corresponds to a relatively flat shoulder in the density profile beyond $r_{\rm cd}$. The bias profile and the velocity profile are also complementary for the practical identification of the depletion scale, with $r_{\rm cd}$ clear in the bias profile for low mass haloes but $r_{\rm id}$ easier to identify in the velocity profile for cluster haloes. As the halo density profile typically has a logarithmic slope close to $-2$ around these scales, this $r_{\rm id}$ can also be identified as the maximum infall velocity location.

We study how $r_{\rm cd}$ depends on multiple halo properties including halo mass, formation time, concentration, shape, spin and environment by binning the bias profiles according to one or two halo properties. $r_{\rm cd}$ in the binned profiles depends strongly on both mass and environment. The formation time and concentration dependence is also clear and non-redundant from the mass dependence, while the dependence on halo spin is the weakest. $r_{\rm cd}$ likely also depends on halo shape beside mass and formation time.

Comparing $r_{\rm cd}$ to the splashback radius defined at the steepest slope location, we find $r_{\rm cd}$ is $\sim 2-3$ times the splashback radius. The ratio between the two increases with decreasing halo size and depends on halo properties in a way similar to the dependence of the NFW concentration on halo properties, with low mass and early forming haloes also having a larger “outer concentration” according to the depletion to splashback ratio. By contrast, the ratio between $r_{\rm cd}$ and the virial radius is $\sim 1-3$ and depends on halo properties in a very different way. The depletion to virial radius ratio is independent of mass ($r_{\rm cd}=2.5r_{\rm vir}$ at fixed mass) but sensitive to other halo parameters, and the trends are different from those for the depletion to splashback ratio. These results reflect that the different radii carry different information about the haloes. As $r_{\rm cd}$ is on average much farther out than the virial and splashback radius, it is subject to more influence from the environmental density.

Comparing $r_{\rm cd}$ to the turnaround radius, which can be interpreted as the outer depletion edge encompassing the infall region around haloes, we find $r_{\rm cd}$ approaches the turnaround in low mass haloes ($M_{\rm vir}\sim 10^{12}\mathrm{M}_{\odot}/h$ and below) that have stopped mass accretion on turnaround scales and reached their maximum turnaround radius. These haloes have well-formed bias troughs but are not surrounded by an infall region due to their completed mass accretion, leading to the agreement between $r_{\rm cd}$ and the turnaround radius. By contrast, massive haloes are surrounded by a clear infall region within the turnaround radius, while the bias trough is much shallower or absent, reflecting their younger age in the mass accretion process. At the high mass end, $r_{\rm cd}$ is thus located within the turnaround radius.

We have further studied the depletion radius in the phase space diagram of haloes, and find it naturally marks the transition between an inner halo structure and the outer environment according to the distribution of radial velocity at different radius. In particular, the inner depletion radius, $r_{\rm id}$, can be interpreted as the radius enclosing a highly complete population of splashback orbits, or the radius of the splashback-bounded halo. Such high percentile splashback radius is not easy to robustly quantify from particle dynamics alone, while our definition via the maximum infall provides a natural alternative.

We find the average density enclosed in $r_{\rm cd}$ is approximately $40$ times the background density of the Universe, independent of the virial mass but dependent on other halo properties. Early forming, highly concentrated, spherical and high spin haloes tend to have a smaller $r_{\rm cd}$ at a fixed mass. Correspondingly, such haloes also tend to have a higher enclosed density, irrespective of mass. For haloes with $M_{\rm vir}>10^{12}\mathrm{M}_{\odot}/h$, the density enclosed in $r_{\rm id}$ is also approximately a constant of $\sim 63$ times the background density.

As the depletion radii mark the outer edges of matter associated with the halo, they are expected to be applicable in the halo model description of the large-scale structure. We demonstrate this by comparing our new radii to the optimal halo exclusion radius found in the recent work of Garcia et al. (2020), who solved for the optimal halo exclusion radius by fitting a flexible halo model to the halo-matter correlation function. We found that the inner depletion radius, $r_{\rm id}$, is in excellent agreement with their optimal exclusion radius, while $r_{\rm cd}$ is in agreement with their “by-eye" halo boundary. As our radii are defined according to physical interpretation of the spatial and velocity distribution around haloes without tuning, it is really exciting to see that our data-driven exploration of the halo radius converges with their model-driven definition of the halo boundary, signalling the convergence towards the beauty and power of a more physical characterisation of structures in the universe.

In this work, we focus on introducing these characterisations and studying their properties in the present day Universe, with some extra efforts in characterising $r_{\rm cd}$. In a follow-up work, we will study the evolution of the bias and velocity profiles to examine the depletion process explicitly, and to gain a dynamical understanding of these radii and their further connections, with more focuses on $r_{\rm id}$.

The authors are grateful to Yipeng Jing for providing access to the CosmicGrowth simulation, and to Hongyu Gao, Zhaozhou Li, Haojie Xu, Ji Yao, Jiaxin Wang, and Lindsay King for helpful discussions. This work is supported by NSFC (11973032, 11890691, 11621303), National Key Basic Research and Development Program of China (No.2018YFA0404504) and 111 project No. B20019. We gratefully acknowledge the support of the Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education. The computation of this work is partly done on the Gravity supercomputer at the Department of Astronomy, Shanghai Jiao Tong University.

The data underlying this article will be shared on reasonable request to the corresponding author.