Cosmic Web Type Dependence of Halo Clustering

The following analysis is based the Millennium simulation (Springel et al., 2005). The first in a series of publicly available, high resolution, cosmological N-body simulations performed within the Millennium Simulation Project (Lemson & Virgo Consortium, 2006). For the original Millennium Simulation Run a $\Lambda$CDM cosmological model with the following parameters has been adopted: $\Omega_{dm}=0.205$; $\Omega_{b}=0.045$; Hubble parameter $h_{100}=0.73$; primordial power spectrum index $n=1$ and $\sigma_{8}=0.9$. The evolution of $2160^{3}$ dark matter simulation particles in a 500 Mpc/h cube is followed resulting in a mass resolution of $8.6\times 10^{8}\,h^{-1}M_{\odot}$.

To identify the dark matter halos, a friends-of-friends (FOF) (Davis et al., 1985) group finder with a linking length b=0.2 is employed and a minimum number limit of 20 particles is imposed. FOF groups are then decomposed into a distinct set of self-bound subhalos using SUBFIND (Springel et al., 2001), the most massive of which is then assigned as the main halo. The current analysis will be restricted to main halos.

The degree of clustering is quantified using the spatial two - point correlation function (e.g., Peebles, 1980). Commonly, one distinguishes between two variations of the two point correltaion function: (i) the two-point auto-correlation function which compares counts of pairs of points drawn from the same sample with the counts drawn from a corresponding random sample as a function of pair separation, while (ii) the two-point cross-correlation function compares counts of pairs drawn from a main and a reference sample with the counts from the main and a corresponding random sample as a function of pair separation. In this paper, we utilise halo samples, defined by mass and web type, as main samples ($Q$-sample) and an a random subset of $\sim 10$ million simulation particles as reference sample ($R$-sample). Thus, the $R$-sample represents the overall matter distribution. Due to the simple geometry of the simulation box, pair counts involoving the random sample ($\mathcal{R}$-sample) can be determined analytically.

In the following we utilise the auto-correlation for the $R$-sample which reflects the clustering behaviour of the overall matter distribution and the matter-halo cross-correlation function. The relation between the two correlation functions is comonly used to determine the bias factor of the halo sample. The auto-correlation function of the $R$-sample is computed by:

where $RR(r)$ indicates the number of pairs from the reference sample with a spatial separation of $r$. $\mathcal{RR}(r)$ is the corresponding number of pair counts drawn from a random sample which is determined analytically. The halo-matter cross-correlation is computed for pairs from the $Q$ and the $R$ sample:

where $QR(r)$ are the mixed pair counts drawn from the halo sample ($Q$) and the reference sample ($R$). As for the auto correlation function, the correctly normalised counts for the main and the random sample, $Q\mathcal{R}(r)$, are computed analytically.

It is worth mentioning, that the cross-correlation functions, discussed here, can be interpreted as the averaged density contrast profiles, scaled by the cosmic mean density, for a given halo sample:

where $\rho$ is the local and $\rho_{\rm mean}$ cosmic mean density. This is a consequence of the fact that our $R$-sample reflects the overall matter distribution. Consequently, the halo-matter cross-pair counts for a given distance bin correspond to the cumulative mass found in all shells centered on the halos of that sample. The division by the halo-random cross-pair counts for the same shells is effectively a normalisation by the mean cosmic density. Finally, subtraction of one gives the density contrast.

We determine the cross-correlation functions on distance scales from 0.1 $h^{-1}$Mpc to 50 $h^{-1}$Mpc. The amplitudes at distances $\lesssim 2h^{-1}$Mpc (depending on the mass of the halo) reflect the density within the halo. At distances of $\sim 10h^{-1}$Mpc outside of the halo volumes, the comparison of the halo - matter cross-correlation function, $\xi(r)_{hm}$ (Eq. 2), and the matter auto-correlation function, $\xi(r)_{mm}$ (Eq. 1), determines the clustering strength of the halo sample relative to the overall clustering of matter. If the amplitudes of the halo correlation function are larger than that of the mass correlation function, the halo sample is said to show positively biased clustering. The reverse behaviour is referred to as anti-biased halo clustering.

Classifying individual halos into web types is done using the Lagrangian classifier given by Fisher et al. (2016) applied to a given halo sample. This is accomplished by constructing the velocity shear tensor (see Hoffman et al., 2012) using techniques borrowed from smoothed particle hydrodynamics (SPH) given as:

where the Greek indices refer to the coordinate system $(x,y,z)$ and the Latin indices represent the halo labels, $n_{i}$ is the halo number density at the location of halo $i$²²2Note: More commonly, SPH approaches utilise the density, $\rho_{i}$, instead of the number-density, $n_{i}$, used here. This requires to include the masses of the particles (halos), $m_{j}$, within the summation. Since we interpret halos as test particles of the underlying velocity field we bypass the reliance on generally ill-determined halo masses by using $n_{i}$ instead of $\rho_{i}$., $\mathbf{v}$ are the halo velocities and $W(|\mathbf{r}_{i}-\mathbf{r}_{j}|,h_{i})$ is the smoothing kernel. We adopt a cubic spline kernel (Monaghan & Lattanzio, 1985), given by Springel (2010) as:

with $q=r/h$ and $W(r,h)=w(q)/h^{3}$. The adaptive smoothing length, $h$, is set individually to ensure exactly 32 neighbours are accounted for. The adaptive approach largely prevents local density and velocity field sampling errors.

Once the 9 elements of the velocity shear tensor have been constructed according to Eq. 4, the velocity shear tensor is diagonalised to obtain its eigenvalues which are then sorted. Defining an eigenvalue threshold value $\lambda^{th}$ allows one to classify the halo web type by counting how many of its eigenvalues are greater than the threshold value. Nodes, filaments, sheets and voids have 3, 2, 1 and 0 eigenvalues greater than the threshold value respectively.

An eigenvalue greater than the threshold value is interpreted as a collapse along its corresponding eigenvector (Zel’dovich, 1970; Hoffman et al., 2012). Expansion along a given eigenvector occurs when the corresponding eigenvalue is less than the threshold value. Nodes collapse and voids expand along all 3 axes. Sheets collapse along one of its axis and expand along two axes while filaments are observed to collapse along two axes and expand along the third. In this work we adopt a threshold value of zero.

Hoffman et al. (2012) considered the threshold value as a free parameter and determined visually as $\lambda^{th}=0.44$. Fisher et al. (2016), based on a Lagrangian approach for the determination of the velocity shear tensor, showed that adding the Hubble flow to the peculiar velocity field (Hoffman’s analysis is based on the peculiar velocity field) allows one to set $\lambda^{th}=0$. Zero seems to be the most natural discriminator between contraction ($\lambda^{th}>0$) and expansion ($\lambda^{th}<0$) along the eigenvectors of the velocity shear tensor.

One subtle aspect of the classification algorithm that needs clarification is the in/exclusion of the Hubble flow when sampling the velocity field for the subsequent determination if the velocity shear tensor. It has been shown in Fisher et al. (2016) that the exclusion of the Hubble flow contaminates the classification scheme by the artificial shrinking of virialised structures. Therefore, inclusion of the Hubble flow is particularly important if the classification scheme is applied to dynamical information from within virialised structures. This is the case if simulation particles or subhalo populations (i.e. not only main halos) are employed for sampling the velocity field.

On the other side, if one focuses on the classification of the main halos based on the large scale environment outside of virialised structures the above mentioned shrinking of the virialised regions due to the exclusion of the Hubble flow does not contaminate the classification scheme. On the contrary, the restriction of the sample to main halos causes sph-kernels to extend to rather large distances where the Hubble expansion begins to dominate. In this case the inclusion of the Hubble flow causes a problem as for distances $\gtrsim 5$ Mpc the ambient main halos show an outward drift due to the Hubble flow resulting in a bias of the classification scheme towards more ’expanding environments’, for instance a filament my be classified as sheet or void or a sheet may be determined as void. The expansion measured this way is not a result of the local gravitational field.

The classification scheme used for the present analysis employs the peculiar velocities of main halos without including the Hubble flow to reveal the impact of gravity on the large scale velocity field. The utilisation of peculiar velocities in combination with the sampling of the velocity field over relatively large distances allows the threshold value for the classification to be set to $\lambda^{th}=0$. This approach produces number fractions for the halo - web type correlations comparable to Hahn et al. (2007a). A zero threshold value permits a natural interpretation as positive and negative eigenvalues correspond to gravitational contraction and expansion respectively.

According to the classification scheme discussed above, the web type number fractions for the classification of the Millennium main halo sample (at $z=0$), including all main halos with top-hat masses $\geq 9\times 10^{10}\,h^{-1}{\rm M}_{\odot}$, are: 9.15% node, 62.61% filament, 27.55% sheet and 0.91% void halos. For the simulation output at $z=1$ using all main halos with top-hat masses $\geq 5\times 10^{9}\,h^{-1}{\rm M}_{\odot}$ (the scaling factor between the two mass limits (at $z=0$ and $1$) is obtained by comparing the characteristic masses of the spherical collapse model at the two redshifts, also see Figure in the Appendix) the corresponding number fraction are: 10.22% node, 64.76% filament, 24.51% sheet and 0.91% void halos. Filament halos constitute the largest subset of halos in the mass range considered here. The statistics indicate a skewness towards web types with a larger number of positive eigenvalues of the velocity shear tensor. It is suggestive that environments of higher infall dimensionality (ID), such as nodes (ID=3) and filaments (ID=2), are more conducive to the formation of halos. On the other hand, one does not expect to find many halos in void region as shear forces are causing an expansion along all 3 dimensions (ID=0). The relatively small number of node halos indicates that only a small volume fraction of the universe shows infall along all three dimensions.

Fig. 1 shows the two point cross and auto-correlation functions for the halo and the overall matter distributions, respectively. The matter distribution is sampled by a subset of  10 Million randomly selected simulation particles (or  0.1% of the total sample). The classification is based on the velocity field sampled by all main halos with masses $\geq 9\times 10^{10}\,h^{-1}{\rm M}_{\odot}$. We use a linear scale below $\xi=1$ to be able to display cases where under-clustering is occurring, i.e. where $\xi<0$. The level of clustering varies significantly with web type as the cross-correlation functions are staggered according to the ID of the web type.

For small $r$ the halo-matter cross correlation functions effectively give the averaged density profiles of the halo sample ($Q$-sample). Thus, the tiered behaviour observed on small scales $\lesssim 1\,h^{-1}\rm Mpc$ reveals a size and mass segregation of the halos corresponding to web type. Void halos are preferentially found among the smallest, least massive halos in the entire sample. The opposite holds for node halos which occupy the largest and most massive halo ranges. Filament and sheet halos occupy the intermediate mass and size ranges.

The staggered arrangement observed at distances $\gtrsim 1\,h^{-1}\rm Mpc$ is a seamless continuation of the behaviour on small scales. At about 10$h^{-1}$ Mpc, where commonly the bias for a given halo population is measured, the comparison of the halo-matter cross-correlation functions with the matter auto-correlation function reveal that halo populations associated with lower ID web types exhibit larger anti-bias. Void halos are shown to be the least clustered web type. Interpreting the halo-matter cross-correlation function for the void halos as averaged density profiles suggests that they reside in under-dense environments of 10s of Mpc in size.

The thin dotted line with the ’+’ markers in Fig. 1 shows the halo-matter cross-correlation function for the entire halo sample. One finds a noticeable similarity to the filament halo cross-correlation function which can be explained by the fact that filament type halos constitute the dominant fraction in the entire halo sample (62.61%). The entire halo sample shows a marginal anti-bias which is a result of the mass distribution within the entire sample, i.e. in numbers the anti-biased halos at the lower end of the mass range dominate over the strongly biased halos at the high mass end.

The tiered behaviour of the halo-matter cross-correlation functions on all scales indicates that halos classified as web types with the highest IDs (IDs = 3 and 2 for nodes and filaments respectively) are usually found in higher density environments. Similarly, low ID (0 and 1 for voids and sheets respectively) corresponds to low density regions. Thus, as expected, there is a direct correlation between ID and density of the environment. Since our web type classification scheme is based on halo positions and velocities alone, this agreement is an independent confirmation of the good correspondence between the cosmic web classification based on density and that based on velocity shear reported by Cautun et al. (2014).

So far the sampling of the velocity field, and ultimately the computation of the eigenvalues of the velocity shear tensor, is based on the entire halo sample, comprised of all main halos $\geq 9\times 10^{10}\,h^{-1}{\rm M}_{\odot}$ in the simulation box. The determination of the velocity gradient uses the 32 nearest neighbour main halos. The most massive halos in the simulation show masses above $10^{15}\,h^{-1}M_{\odot}$. In general, it is expected that high mass halos are surrounded by a plethora of low mass halos (10 - 10000 smaller in mass). Thus, the velocity shear tensor associated with high mass halos may frequently be based on a set of low mass halos close to, but outside of, the virial radius.

These regions are know to be contaminated by backsplash halos, i.e. halos which have been part of the main halo before being ejected (e.g., Moore et al., 2004; Warnick et al., 2008). To check the effect of the backsplash halos we have computed the shear tensor excluding neighbour halos within two times the virial radius and obtained very similar results, with the encouraging exception that the very high mass void halos disappear. In the following section, which presents the main results, we only use narrow mass bins for the sampling of the velocity field and calculation of the velocity shear tensor which suppresses the contribution of backsplash halos as they are typically orders of magnitude less massive than the main halo. Thus, we do not investigate the impact of backsplash halos any further.

Akin to the work in Gao et al. (2005b), we examine the mass dependence of halo clustering in four narrow mass bins. Within each mass bin we subdivide the halo sample according to web type and investigate the differences in the associated halo-matter cross-correlation functions. The sampling of the velocity field and consequently the determination of the cosmic web type is confined to halos in the same mass bin. For the highest mass bin the velocity shear tensor is computed based on halos more massive than $2\times 10^{13}\,h^{-1}\rm M_{\odot}$ resulting in an average distance to the 32nd neighbour of $35\,h^{-1}\rm Mpc$. The confinement to narrow mass bins takes the different sizes of the halos into account, i.e. the determination of the web type of low mass halos is based on low mass neighbour halos which on average have smaller separations, whereas for the determination of the web type of high mass halos only the halos in the high mass bin are used which show on average larger separations. Somewhat counterintuitive is the appearance of a small fraction (0.3%) of void halos in the highest mass bin. This average expansion of the ambient halo-velocity filed as reliably measured by the SPH-kernel may be caused by halos of similar size which just experienced a close encounter and presently are on a separation course or by the gravitational impact of massive structures nearby ($\sim 35\,h^{-1}\rm Mpc$).

As a cross check, we have used the entire halo sample to compute the web types for the halos in the various mass bins and find qualitative agreement despite the contamination by backsplash halos in this case. Since the former approach determines web type and clustering based on the same halo sample we have decided to present the results of this, methodologically more consistent approach.

Fig. 2 shows the halo-mass cross-correlation functions for all halos in the indicated mass bins, again on linear scales for $\xi(r)\leq 1$. As expected, clustering increases with mass. The correlation functions are tiered. The charatersitic mass derived from the sperical collapse model within the adopted cosmology is $M_{*}(z=0)approx7\times 10^{12}h^{-1}M_{\odot}$. In the lower mass bins, less than the characteristic mass, a anti-bias, measured at scales from $5$ to $20\,h^{-1}\rm Mpc$, can be observed. Conversely, the highest mass bin is defined greater than $M_{*}$ and one observes a large bias factor. Halos in the second highest bin have masses most similar to $M_{*}$ and they display very little bias.

In Fig. 3, we calculate the web type clustering per mass bin. All panels of Fig. 3 display the same behaviour as in Fig. 1 whereby the web type cross-correlation functions are tiered by web type. At small scales, the halo-matter cross-correlation function measures the average density profiles of the halo sample and it can be seen that the profiles converge independent of web type. The deviation seen in the highest mass bin is a result of the large mass range covered by this mass bin, which causes a behaviour reminiscent of Fig. 1 where the correlation function for the entire mass range is displayed.

Arguably the most striking feature in Fig. 3 is the similarity of the correlation functions for a given web type independent of mass. In all four mass bins the node halos are strongly biased at scales $\gtrsim 5\,h^{-1}\rm Mpc$, the filament halos follow the matter auto-correlation function closely and sheet and void halos are strongly anti-biased with the zero-crossing occurring at much smaller scales, when compared to the entire halo sample. The exact scale of the zero-crossing seems to be influenced by the average halo radius in the specific sample but the maximal negative amplitudes are very similar for a given web type independent of mass.

On linear scales, $\sim 10h^{-1}$ Mpc, the clustering of low mass halos ($\lesssim 10^{13}M_{\odot}$) has been shown to depend substantially on halo formation time (Gao et al., 2005b). The clustering of halos which currently have the same mass depends on their formation time. Halos that assembled at high redshift are substantially more clustered than those that assembled more recently. Subsequently, many other properties of halos, such as concentration and velocity anisotropy parameter, have been shown to have an impact on halo clustering besides the halo mass (see introduction). Fig. 3 clearly demonstrates that the web type of halos also has a strong impact on halos clustering. Very simular conclusions can be drawn from the correponding $z=1$ results displayed in Fig. 5 in the Appendix.

The above analysis has revealed that same mass halos show different clustering behaviour depending on their web type. In order to gain insight into the mass dependence of the halo-matter cross-correlation functions for the entire halo sample, as seen in Fig. 2, we now investigate which web types are responsible for the overall clustering per mass bin.

Fig. 4 illustrates the web type fractions per mass bin (for visual impression, the last bin mass value is set to $10^{14}$, otherwise plotted using the midpoint of the respective bin). The fraction of halo web types is largely mass dependant. Node halo fractions are monotonically increasing with mass so that in the highest mass bin they account for the second largest fraction of all halos, followed by sheets and voids. Filaments are dominant in all mass bins. In the low mass regime, filament halos are seen to occupy the largest fraction of halos followed by sheet halos. Sheet and void fractions decrease with mass, the latter being less than 1% for all mass bins.

Fig.4 suggests that filaments and sheets determine the overall clustering in the lower mass ranges ( $\sim 10^{11}$ to $\sim 10^{12}\,h^{-1}\rm M_{\odot}$) causing anti-bias. For intermediate mass ranges ( $\sim 10^{12}$ to $\sim 10^{13}\,h^{-1}\rm M_{\odot}$), filament halos dominate which results in marginal bias. Finally, the substantial fraction of node halos along with filament halos in the high mass ranges ($\gtrsim 5\times 10^{14}\,h^{-1}\rm M_{\odot}$) generates positively biased halo clustering. The thin lines in Fig.4 display the web type fractions for $z=1$, where the mass bins are rescaled by the ratio of the characteristic masses of the spherical collapse model, $M_{*}(z=0)\approx 7\times 10^{12}h^{-1}M_{\odot}$ and $M_{*}(z=1)\approx 3.5\times 10^{11}h^{-1}M_{\odot}$. The agreement between $z=0$ and $z=1$ shows that the contribution of web types to the overall halo clustering is relatively constant over time.