Dark matter haloes: a multistream view

TL;DR

This paper introduces a new method for identifying dark matter haloes in cosmological simulations by analyzing the multistream field, which captures the collisionless nature of dark matter and its nonlinear growth, leading to improved halo detection.

Contribution

The paper proposes a novel framework utilizing the multistream field to detect dark matter haloes, offering a new perspective beyond traditional particle-based methods.

Findings

Effective identification of halo substructures using multistream maxima.

Convexity conditions isolate potential halo regions.

Enhanced understanding of dark matter clustering in simulations.

Abstract

Mysterious dark matter constitutes about 85% of all mass in the Universe. Clustering of dark matter plays the dominant role in the formation of all observed structures on scales from a fraction to a few hundreds of Mega-parsecs. Galaxies play a role of lights illuminating these structures so they can be observed. The observations in the last several decades have unveiled opulent geometry of these structures currently known as the cosmic web. Haloes are the highest concentrations of dark matter and host luminous galaxies. Currently the most accurate modeling of dark matter haloes is achieved in cosmological N-body simulations. Identifying the haloes from the distribution of particles in N-body simulations is one of the problems attracting both considerable interest and efforts. We propose a novel framework for detecting potential dark matter haloes using the field unique for dark matter -- multistream field. The multistream field emerges at the nonlinear stage of the growth of perturbations because the dark matter is collisionless. Counting the number of velocity streams in gravitational collapses supplements our knowledge of spatial clustering. We assume that the virialized haloes have convex boundaries. Closed and convex regions of the multistream field are hence isolated by imposing a positivity condition on all three eigenvalues of the Hessian estimated on the smoothed multistream field. In a single-scale analysis of high multistream field resolution and low softening length, the halo substructures with local multistream maxima are isolated as individual halo sites.

Tables

Table 1. Table 1: Statistics of the Hessian eigenvalues in the halo candidates

Statistics	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$
Minimum	$1.5$	$0.5$	$1.3\times {10}^{-2}$
Maximum	$1.7\times {10}^3$	$1.5\times {10}^3$	$1.1\times {10}^3$
Median	$10.5$	$5.5$	$1.9$

Table 2. Table 2: Number of λ 3 subscript 𝜆 3 \lambda_{3} -haloes identified at smoothing of n s ​ t ​ r ​ ( 𝒙 ) subscript 𝑛 𝑠 𝑡 𝑟 𝒙 n_{str}(\mn@boldsymbol{x}) at different scales.

$N_p$	$0.20h^{-1}\text{Mpc}$	$0.39h^{-1}\text{Mpc}$	$0.78h^{-1}\text{Mpc}$
${128}^3$	$5181$	$4492$	$2923$
${256}^3$	$27929$	$18221$	$7897$

Table 3. Table 3: Local maxima and minima of n s ​ t ​ r subscript 𝑛 𝑠 𝑡 𝑟 n_{str} in each of 4492 haloes. The highest n s ​ t ​ r subscript 𝑛 𝑠 𝑡 𝑟 n_{str} values in the interior of haloes span over a large range of values. Low values of n s ​ t ​ r subscript 𝑛 𝑠 𝑡 𝑟 n_{str} in haloes, which are generally near halo boundaries, have a median of 3.

Statistics	$n_{str}^{high}\left(H_i\right)$	$n_{str}^{low}\left(H_i\right)$
Minimum	$7$	$3$
Maximum	$2831$	$459$
Median	$17$	$3$

Table 4. Table 4: Number of haloes, N H subscript 𝑁 𝐻 N_{H} detected by the three halo finder algorithms in the two simulations of L = 100 ​ h − 1 ​ Mpc 𝐿 100 superscript ℎ 1 Mpc L=100h^{-1}\text{Mpc} with different mass resolutions. Values shown for λ 3 subscript 𝜆 3 \lambda_{3} -haloes are done in the multistream fields with refinement factor of 2, and smoothing scale equal to the diagnostic grid size.

$N_p$	$N_H^{\mathrm{AHF}}$	$N_H^{{\lambda }_3}$	$N_H^{\mathrm{FOF}}$
${128}^3$	3374	4492	5440
${256}^3$	24710	27929	35765

Equations

$${\text{H}}_{ij}(f)=\frac{{\partial }^{2}f}{\partial x_i\partial x_j}$$

Full text

††pagerange: Dark matter haloes: a multistream view–15

1 Introduction

The web-like distribution of matter initially revealed by redshift surveys (with less than 300 galaxies by Gregory & Thompson 1978 and around 1000 galaxies by de Lapparent et al. 1986) and numerical modeling (using N-body simulations of around 30000 particles by Shandarin 1983a and Klypin & Shandarin 1983) pioneered morphological investigations of the cosmic web structures (see Bond et al. 1996, also reviews by Shandarin & Zeldovich 1989 and van de Weygaert & Bond 2008). Detailed mapping of the Universe has crossed three million objects today, by catalogues such as the Sloan Digital Sky Survey (SDSS; Albareti et al. 2016). The upcoming Large Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. 2009) is expected to probe the nature of dark matter using several billion galaxies. On the other hand, cosmological simulations have improved immensely in several aspects – numerical techniques, parallelization schemes, inclusion of various physical processes, volume and resolution (some of these developments are summarized in Bertschinger (1998) and Bagla & Padmanabhan (1997). Modern state-of-the-art simulations like the Illustris Project Vogelsberger et al. (2014), the EAGLE project Schaye et al. (2015) and Q-Continuum Heitmann et al. (2015) use more than a billion dark matter particles. Finally, the ever improving data analysis techniques have resulted in new and sophisticated density estimators, geometrical and topological indicators. A plethora of algorithms for identifying and characterizing dark matter structures have emerged in last two decades (a summary on cosmological data analysis is highlighted in van de Weygaert & Schaap 2009). Considering all these improvements, it is worth noting that the proto-structures detected in the modern simulations are qualitatively similar to the quasi-linear description of clustering by Zeldovich Approximation (ZA; Zeldovich 1970). Location and properties of these structures, i.e., the voids, walls, filaments and haloes maybe inconsistent across different structure finding algorithms, but that is primarily due to varied definitions.

Most of structure finders are halo finders only and majority of them are stemmed from three underlying algorithms. One of them is the SO (Spherical Over-density) halo finder that defines halos as spherical regions whose mass density exceeds the mean density by a specified factor (Press & Schechter, 1974). Another is the FOF (Friends-of-Friends) halo finder describing haloes as the groups of particles separated less than a specified linking length often chosen as 0.2 times the mean particle separation (Davis et al., 1985). The FOF can be also used for identifying filaments and walls/pancakes by increasing the linking length (Zeldovich et al. 1982, Shandarin 1983b, Shandarin et al. 2010). Finally the DENMAX (DENsity MAXimum) halo finder assumes that the halos are the peaks of the density fields and thus selects the particles concentrated in the vicinity of the density maxima (Bertschinger & Gelb, 1991). One of the common features of these techniques is that all three are based on density, in one form or another. And all of them depend on free parameters that are chosen chiefly on the ‘merits principle’ (Forero-Romero et al., 2009) rather than on physics. Over the years all three kinds of the halo finders have been experiencing various modifications and improvements. A few examples from a long list of these modifications may include:

(i) Improvised techniques of generation of the density field from the particle positions, and finding spherically bound over-densities (Lacey & Cole 1994, Jenkins et al. 2001, Evrard et al. 2002, Weinberg et al. 1997, Neyrinck et al. 2005, Knollmann & Knebe 2009, Sutter & Ricker 2010, Planelles & Quilis 2010 etc.)

(ii) Adaptive methods controlling the linking length in methods using FOF (Davis et al. 1985, van Kampen 1995, Gottlober et al. 1999, Springel et al. 2001b, Habib et al. 2009, Rasera et al. 2010 etc.)

(iii) Adaptive methods for searching the positions of density maxima (For example, Klypin et al. 1999.)

(iv) Generalization of FOF and DENMAX techniques to six-dimensional phase space, and many others (such as 6DFOF by Diemand et al. 2006 and ROCKSTAR by Behroozi et al. 2013 use velocity-position space with parameters analogous to linking-length.)

(v) Computing hierarchical tree of clusters in the phase-space such as the Hierarchical Structure Finder MacIejewski et al. (2009), and the 6-D Hierarchical Over-density Tree Ascasibar (2010). (vi) Hybrid algorithms: frameworks such as the Hierarchical Bound-Tracing algorithm Han et al. (2012) and SURV Giocoli et al. (2010) identify haloes at multiple time steps from the simulation to find prospective sub-haloes. In addition, there are HOP methods by Eisenstein & Hut (1998), Tweed et al. (2009) and Skory et al. (2010).

A detailed comparisons of several halo/sub-halo finders is provided in Knebe et al. (2011), Knebe et al. (2013), Onions et al. (2012). In a nice summary discussing these developments as well as describing a few new suggestions they concluded that there was no general consensus for a precise definition of a halo or a sub-halo. Consequently, there were different estimates of number of haloes, halo mass functions, halo centre locations, boundaries and other parameters.

There are significant concerns with SO, DENMAX and FOF algorithms - both in terms of underlying mechanisms of halo formation and the parameters used in halo identification. SO is motivated by the analytical toy model of the collapse of a top-hat spherical density perturbation. Parameters of the virial radii $r_{vir}$ and virial mass $M_{vir}$ are determined by the regions with density $\rho_{vir}\geq\Delta_{vir}\times\rho_{b}$, where $\rho_{b}$ is the background density of the simulation box. $\Delta_{vir}$ is generally taken around $180$ or $200$, derived for an isolated spherical collapse model, and it varies for different cosmologies and redshift. The peaks in CDM models not only aspherical, but their collapse is subject to tidal forces, mergers and presence of sub-structures - none of these complexities are weighed in the spherical collapse model.

For FOF, the free parameter of linking length is generally taken as $b=0.2$ times the mean separation of particles at $z=0$. This inter-particle separation corresponds to $\Delta_{vir}\approx 180$ if the halo has an isothermal density profile, $\rho\propto r^{-2}$. Using percolation theory, More et al. (2011) argued that this corresponds to a rather wide range of over-densities depending on halo mass and density profiles. They found out that $b=0.2$ corresponds to local over-density $\delta$ within the enclosed halo to be in the range of 250 to 600. Moreover, the resulting FOF-haloes need not have a compact geometry: often the haloes are irregularly shaped, which is unlikely if the halos are virialized. Hence modern algorithms re-define the halo boundaries by excluding particles using post-processing techniques. In recent simulations with clear delineation of walls and filaments Angulo et al. (2013), $b=0.05$ was used for finding FOF-haloes since the traditional value of $b=0.2$ corresponded to structures that percolate into the web structure.

Absence of dynamical traits in the FOF and SO algorithms are arguably more crucial. In phase-space, the halo collapse models show collisionless DM particles in oscillatory motions about a core, at successive foldings of the phase-space sheet. The velocity field within each oscillatory spiral is multi-valued in physical space. Incrementing multistream shells, separated by caustic surfaces sequentially trace the structures of the cosmic web - walls, filaments and the haloes. Majority of the mass accretion into the haloes along the filaments: from lower multistreams into higher. Thus the DM haloes are not independent of filaments around them, and the hierarchical layers of multistreams represents this portrait precisely. This picture of structure formation was initially theorized using ZA Zeldovich (1970) and in context of caustics Arnold et al. (1982) as well as in the Adhesion Approximation (Gurbatov, S.N., Saichev, A.I., Shandarin et al. 1989, Kofman et al. 1992). Shandarin & Zeldovich (1989) reviewed gravitational evolution of density perturbations in the context.

It has been demonstrated that the multistream field in Eulerian space can be computed directly from the Lagrangian sub-manifold (Shandarin et al. 2012 and Abel et al. 2012). About 90% of the field is single-streaming voids, and the rest of the volume comprises of multistream walls, filaments and haloes. Ramachandra & Shandarin (2015) found the multistream value of $n_{str}\approx 90$ corresponds to virial density $\Delta_{vir}=200$. On the other hand, DM particles are identified by (Falck et al., 2012) as belonging to haloes if they undergo flip-flop along 3 orthogonal axes. These analyses have opened up a new avenue in studies of halo formation, both qualitative and quantitative. Re-investigations of halo spins, physical radii of the halo, sub-structure in the light of streaming phenomena have shown that the halo structures and formations are more complicated than previously believed. Vogelsberger & White (2011) investigated the distribution of streams in small haloes at various redshifts. They concluded that tracking caustics and streams is better than density, since density fields are noisy in the dense inner regions of haloes. In another study, More et al. (2015) argued that the ‘splashback radius’ - distance from the halo core to the first caustic enumerated from outside - is a better physical indicator of DM halo boundary than the virial radius (also see the discussion on turn-around radius of bound objects by Lee & Yepes 2016). Angulo et al. (2013) also agree with the view that the locally overdense regions correspond better with the volumes within the first caustic than the virialized DM clumps. Recent toy model of anisotropic halo collapse by Neyrinck (2016) considers intersecting multistream filaments forming spinning nodes. Ramachandra & Shandarin (2017) showed that the virial surfaces of FOF haloes have varying number of streams, including single-streams. Study by Shandarin & Medvedev (2017) delineated the rich sub-structure of haloes using another derivative parameter from the Lagrangian sub-manifold called the ‘flip-flop’ defined on the Lagrangian space.

In this paper, we identify potential haloes by utilizing multistreaming as the governing dynamical phenomena. A review of the DM particle clustering in a one-dimensional dimensional universe is made in Section 2, and the concept of multistream field is extended to higher dimensions. The multistream field is computed on the cosmological simulations described in Section 3. The halo identification framework in this field is described in Section 4. This algorithm isolates convex regions of the multistream field using Hessian eigenvalues, each enclosing a local multistream maximum. Without employing any non-local thresholds that several halo finders generally use, these convex multistream regions are identified as potential halo sites. We also illustrate the significance of multistream refinement and softening scales in finding subhaloes. However, this paper does not focus on adaptive multi-scale analyses for substructure studies. A few properties of the multistream haloes are discussed in Section 5, and comparison of these haloes with AHF and FOF algorithms is done in Section 6. We also discuss the spatial distribution of the dark matter haloes along the percolating web structure.

2 Phase-space representation of gravitational clustering

We begin with a simple illustration showing the formation of a few haloes in a one-dimensional simulation. Dark matter clustering in a (1+1)-dimensional phase-space $(\mn@boldsymbol{p},\mn@boldsymbol{x})$ (where $p$ is the momentum and $x$ is the co-moving Eulerian coordinate) at four successive time steps is shown in the top panels of Figure 1. The lower panels show the corresponding multistream field (Shandarin et al. 2012 and Abel et al. 2012) $n_{str}(\mn@boldsymbol{x},z)$ (red) and density field $\rho(\mn@boldsymbol{x},z)$ (gray). At $z_{1}$ (left-most panel), velocity is single-valued in Eulerian co-ordinates shown, except at a small three-stream region near $\mn@boldsymbol{x}=5\pi/4$. This is the first instance of multistreaming in the region, which was previously had $n_{str}=1$ throughout. The interface of $n_{str}=1$ and $n_{str}=3$ regions is also the location of the first caustic. On the other hand, the density calculated at a high resolution shows variations, even in the mono-streaming regions. The variations are especially more pronounced around the caustic (near $\mn@boldsymbol{x}=5\pi/4$).

The gravitational clustering is more evolved in the two center panels ($z_{2}$ and $z_{3}$) with three prominent phase-space spirals. The regions between the spirals have sparsely distributed dark matter particles, and have $n_{str}=1$. Each spiral corresponds to a location of gravitational collapse with $n_{str}>1$ region, and higher density. A few of these regions within three-streaming regions are elevated to $n_{str}=5$. The corresponding density field is not only noisier, but also reaches very high values at the caustics. This is also a primary distinguishing feature between mass density fields and multistream fields: At the locations of caustic, the density (regardless of how it is calculated) is not smooth Vogelsberger & White (2011). Computational limitations on simulation resolutions and refinement of density calculations soften the fields, exceptionally at the zero volume measure regions of caustic surfaces. On the other-hand, multistream values are increased by finite values at caustic surface locations - This property is preserved at higher simulation resolutions and any refinements of multistream field calculations - although $n_{str}$ may be resolved enough to have intermediate even-values. Multistream fields are also intrinsically discrete valued, which is not true with density fields. Discreteness of multistream fields is discussed in more detail in Ramachandra & Shandarin (2017).

The right-most panel in Figure 1 shows the final structure at $z=0$. Two large spirals have spatially merged. These collapse environments are naturally very complex, with an increased number of successive caustic formation and merging. The corresponding velocity streams also show a more complicated structure. Clearly, the multistream field has a saddle point that is not as low as $n_{str}=1$. This poses a bigger problem in the context of most of halo detection algorithms, and we discuss this in Appendix A.

2.1 Collapse in higher dimensions

Extending the above results of one-dimensional collapse into higher dimension is vital, primarily in the context of halo formation. The individual spiral collapses in the one-dimension happen at a small region (left-most panel in Figure 1), and the region grows by volume, whilst increasing the spiral twists within. This is in contrast with the theoretical top-hat spherical model of halo formation when the shell crossing would not happen until the final moment of the collapse of the entire halo into a point-like singularity. Thus all shell crossings happen at a single point and at a single instant of time. The collapse of an isolated, spherically symmetric density peak is a very exceptional case, because every spherical shell feels only the force due to interior mass until it collapses into the caustic region. The collapse of the real peak proceeds in the field generated by the mass distribution - in both the mass within the forming halo, and the mass outside the halo.

The collapse of a uniform ellipsoid also results in a simultaneous collapse of the entire ellipsoid however this time not into a point but into a two-dimensional ellipse (Lin et al. 1965, Icke 1973, Eisenstein & Loeb 1995). Another customarily used spherical model of halo formation by Fillmore, J. A., Goldreich (1984) and Bertschinger (1985) does not consider the initial collapse at all. Instead it assumes self-similar initial conditions and the halo at advanced stage with formally infinite number of spherical caustic shells.

The ‘core’ in a collision-less dark matter collapse (in Figure 1) is a region where a multistream region is first formed due to caustic formation. This is conceptually similar to a shell crossing. However, there are successive caustic formations that follow the first shell crossing, and they are not limited to the halo cores. Each caustic increases the multistream value within by a finite number. The cores of the multistream haloes obviously have the local maxima of velocity streams in Eulerian coordinates. On the contrary, mass densities have infinite values at the caustics surfaces, including the core. Discontinuities in densities at these regions of sharp multistream transitions are clearly seen if the mass and spatial resolutions were sufficiently high( see two-dimensional simulations by Melott & Shandarin 1989 as well as in three-dimensional simulations by Hahn et al. 2013, Angulo et al. 2016, Hahn & Angulo 2016 etc.).

In three-dimensional simulations, the Lagrangian sub-manifold twists in complicated ways in a six-dimensional phase space. This is due to complexities involving caustic formations in higher dimensions, which is true even in the ZA, see Arnold et al. 1982 and Hidding et al. 2014 for detailed analyses of caustic formation. The resulting intricate geometrical structures can be characterized by the $n_{str}$ field. Nearly $90\%$ of the volume in N-body simulations are single-streamed voids at $z=0$ (Shandarin et al. 2012, also see Falck & Neyrinck 2015 for a percolation analysis of single-streaming voids). From the visualizations in Ramachandra & Shandarin (2015) and percolation analysis of Ramachandra & Shandarin (2017), we also know that the $n_{str}=3$ regions mostly form connected wall-like structures, unlike the isolated patches as seen in one-dimensional simulations of Figure 1. The structures become predominantly filamentary at higher thresholds of $n_{str}\gtrsim 17$. Subsequently, the regions around the multistream maxima have isolated closed surfaces (for example, in Figure 2), which may be identified as halo locations.

Caustic formations and mass accretion are also seen to occur more along the higher streams, which makes the haloes non-spherical, with the alignment generally determined by a complicated interplay of the intensities of the streams from neighboring filamentary structures. Number of streams corresponding to the dark matter halo also has a local environment dependence. The three small haloes in Figure 2, where the number of streams are higher than the neighbouring filaments, are aligned along three intersecting filaments. Halo environment studied in Ramachandra & Shandarin (2015) show similar hierarchical variation in $n_{str}$ values. The halo environments are thus very complicated, and empirical thresholds (on multistream or density fields) may not account for all the haloes uniformly. Hence we use a local geometrical method to identify potential haloes in multistream fields.

The first non-linear structures in the web have $n_{str}=3$. By visual inspection, these regions generally form a fabric-like open structures that resemble walls. The surface contours of higher $n_{str}$ are embedded within the walls. Figure 2 shows a filamentary structure of the web at $n_{str}\gtrsim 17$. The figure also shows regions around local maxima of the multistream field, which are generally located at the intersections of filaments.

3 Simulations and tools

The emphasis of this paper is to demonstrate the use of multistream field in identifying potential dark matter haloes, and not a full statistical analysis of halo properties. For this purpose, we have run simulations at two different mass resolutions (number of particles $N_{p}=128^{3}$ and $256^{3}$, and respective mass of particles, $m_{p}=3.65\times 10^{10}h^{-1}M_{\sun}$ and $4.57\times 10^{9}h^{-1}M_{\sun}$), with the same periodic side length $L=100h^{-1}Mpc$. The gravitational softening length $\epsilon=20h^{-1}kpc$ and $10h^{-1}kpc$ for low and high resolution simulation respectively. The initial conditions are generated by MUSIC (Hahn & Abel, 2011) with the transfer function from Eisenstein & Hu (1998) at a redshift of $z_{ini}=80$. The $\Lambda$CDM cosmological simulation is run using GADGET-2 (Springel 2005 and Springel et al. 2001a) is similar to the ones used in Ramachandra & Shandarin (2017). The cosmological parameters used in the simulation are $\Omega_{m}=0.276$, $\Omega_{\Lambda}=0.724$, the Hubble parameter, $h=0.703$, the power spectrum normalization, $\sigma_{8}=0.811$ and the spectral index $n_{s}=0.961$.

Multistream field $n_{str}(\mathbf{x})$ is calculated on the GADGET-2 snapshots at $z=0$ using the tessellation scheme by Shandarin et al. (2012). The multistream field can be computed at the native resolution of the Lagrangian grid of the simulation, i.e., at refinement factor of $l_{l}/l_{d}=1$ (where $l_{l}$ is the inter-particle separation in Lagrangian grid and $l_{d}$ is the side length resolution of diagnostic grid). Arbitrarily high refinement factors can be utilized in computing multistream fields as well, for example $l_{l}/l_{d}=8$ for the halo multistream environment shown in Figure 2. For analysis of full simulation boxes, we restrict $l_{l}/l_{d}$ to $1$ and $2$.

Two halo finders are also used to identify potential haloes with 20 or more particles at $z=0$: a classic Friends-of-Friends method (FOF-Davis et al. 1985) using a popular linking length, $b=0.2$ (e.g. Frenk et al. 1988 and Lacey & Cole 1994) and the Adaptive Mesh Investigations of Galaxy Assembly (AMIGA halo finder or AHF-Knollmann & Knebe 2009, Gill et al. 2004). Halo catalogue from these halo finders are used to compare with our implementation of halo detection in the multistream field. The haloes candidates from AHF and FOF algorithms are hereafter referred to as AHF-haloes and FOF-haloes respectively.

4 Haloes in the multistream field

We intend to identify haloes in the $n_{str}(\mn@boldsymbol{x})$ field instead of using just the position coordinate data. While the eigenvalue analysis itself is done at a chosen time, the multistream field inherently has data from six-dimensional Lagrangian space $(\mn@boldsymbol{q},\mn@boldsymbol{x})$ that contains the full dynamical information,similar to the phase-space sheet albeit in a different form. Dynamical history that is embedded in the multistream field is crucial in understanding the strength of gravitational binding of the particles. A physically motivated distinction between void and gravitationally collapsed regions – voids are the regions with a single stream – is a unique feature of multistream analysis (Shandarin et al. 2012 and Ramachandra & Shandarin 2017). Thus the haloes detected from local maxima of the $n_{str}$ field can be ensured to be away from the mono-streaming voids. Methods based on linking-length or density fields may not be able to ensure that all the particles in haloes are away from voids (as shown for FOF haloes in Ramachandra & Shandarin 2017) .

Numerical analyses of scalar fields generally depend on resolution as opposed to particle coordinates analysis tools like FOF. The multistream field conveniently has an advantage of being less noisy than mass density (Shandarin et al. 2012, also see the Appendix in Ramachandra & Shandarin 2017 ).

4.1 Hessian of multistream fields

Hessian matrix $\textbf{{H}}(f)$ of a scalar field $f$ involves local second-order variations in three orthogonal directions. Each element of the Hessian matrix $\textbf{{H}}_{ij}(f)$ (where $i$ and $j$ can be any of $x$, $y$ or $z$ directions) is given by Equation 1.

[TABLE]

By choosing a $f=-n_{str}(x)$ (smoothened if necessary), local multistream variations can be diagnosed. The Hessian matrices at each point on the configuration space are always symmetric matrices, resulting in real eigenvalues. The Hessian eigenvalues in multistream fields differ from that in density, gravitational potential or velocity shear tensor. We refer the readers to Ramachandra & Shandarin (2017) for an extensive analysis on multistream Hessians and their geometrical significance. Some of the salient features of Hessian eigenvalues of multistream field are as follows:

(i) Every element of Hessian matrices $\textbf{{H}}(-n_{str})$, and consequently the eigenvalues $\lambda_{i}$’s are zero in single-stream voids. Even if the multistream field is a smoothed, the eigenvalues peak at zero. This property is unique to multistream fields. Eigenvalues of Hessians of density Aragon-Calvo et al. (2007), velocity shear tensor Libeskind et al. (2013) do not peak at zero, and the eigenvalues of deformation tensor are negative in voids as a result of continuity equation (shown in Zel’dovich formalism as well).

(ii) The eigenvalues of these Hessian matrices are always real, and depending on if their values are positive or negative, one may infer local geometrical features in the field. For our choice of $-n_{str}(x)$ as the domain of Hessian, at least in principle, the conditions for geometric criteria are: $\lambda_{1}>0>\lambda_{2}\geq\lambda_{3}$ for locally flat regions, $\lambda_{1}\geq\lambda_{2}>0>\lambda_{3}$ for locally tubular structures and $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}>0$ for clumped blobs.

(iii) Convex neighbourhoods around local maxima of the multistream field are isolated by the positive definite Hessian matrices.

(iv) The resulting Hessian eigenvalues characterize the geometry in a four-dimensional hyper-space of $(-n_{str},x,y,z)$. The boundary of a region with $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}>0$ is a closed convex contour in this hyper-space, and thus it’s projection on to the three-dimensional Lagrangian space is also closed and convex.

Of the three geometries that are characterised by the eigenvalue conditions, we investigate the convexity of multistream contours in the context of halo finding in the section below.

4.2 Halo finder algorithm

Our goal is to isolate the locations of convex geometries in the multistream flow field. Prospective regions of the halo locations in the web structure are selected by positive definite condition of the Hessian $\textbf{{H}}(-n_{str})$: $\lambda_{1}>0$, $\lambda_{2}>0$ and $\lambda_{3}>0$, or simply the smallest of the eigenvalues, $\lambda_{3}>0$. We also filter out the regions if the multistream values inside do not suggest gravitational collapse into haloes. The sequence of our halo detection framework is listed below:

The multistream flow field is calculated on a diagnostic grid. The number of tetrahedra that encompass each vertex in the grid gives the $n_{str}$ field. Top left panel of Figure 3 shows the multistream web structure in a slice of the simulation with $n_{str}>1$ in gray and $n_{str}\geq 7$ in blue. 2. 2.

The discrete multistream flow field is smoothed in order to reduce numerical noise. We have used Gaussian kernel for smoothing in our analysis. Effect of smoothing scales in the halo identification is discussed in Section 4.3. 3. 3.

Second order variations of the smoothed $-n_{str}(\mn@boldsymbol{x})$ is computed at each point in the field. This gives symmetric Hessian matrices for this field whose eigenvalues are real. Ordered eigenvalues of the Hessian, $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$ are calculated. The $\lambda_{3}$ field is shown in the top right panel of Figure 3. 4. 4.

Using segmentation techniques, each region with $\lambda_{3}$ strictly greater than [math] within $n_{str}\geq 3$ regions of multistream field are isolated and labelled. This condition for each halo candidate guarantees that it is in the region where at least one gravitational collapse happened within the halo boundary. Mass particles belonging to these regions are shown shown as dark spots in in the top right panel of Figure 3. 5. 5.

The multistream field has a range of values within the isolated sites. We impose constraints on the isolated regions to rule out the labels with low multistream values. The local maxima of $n_{str}$ inside each halo must be at least 7. By this restriction, it is ensured that the halo sites with three Lagrangian sub-manifold foldings are selected. Bottom left panel of Figure 3 shows patches that are discarded in red. The resulting $\lambda_{3}$-haloes are shown in the bottom right. 6. 6.

In our comparisons with other halo finders in Section 6, we also used an additional constraint on the minimum number of mass particles in the haloes to be 20 - which is generally used as a criteria in several halo finders.

For the illustration halo detection framework in this section, we have calculated the number-of-streams at refinement factor of $2$ and smoothing scale of $0.39h^{-1}Mpc$ (equal to the grid length of the multistream field) for the simulation box of $128^{3}$ particles and size $L=100h^{-1}Mpc$ . Hessian matrices and eigenvalues are calculated on the same diagnosis grid. Results of the halo detection scheme for simulation box of higher mass resolution, and different smoothing factors are discussed in Sections 6 and 4.3. Hereafter we refer to the potential dark matter haloes detected from the Hessian analysis of the multistream field as $\lambda_{3}$-haloes for brevity.

Applying the above scheme on the simulation with side length of $100h^{-1}$ Mpc and $128^{3}$ particles (with cosmological parameters mentioned in Section 3), we detected approximately 50000 regions satisfying $\lambda_{3}>0$ within the non-void in the multistream field of refinement factor $l_{l}/l_{d}=2$ and smoothing scale of grid length, i.e, $0.39h^{-1}Mpc$. We filtered out the segments with local maxima of $n_{str}<7$, and around 14000 regions remained as prospective haloes. Majority of these regions have less than 20 particles, which are excluded in the halo catalogues. On the whole, our algorithm detected about 4500 haloes with more than 20 particles in the entire simulation box. We have not applied virialisation to define the halo boundaries. A more detailed study of halo edges, and comparison with that of FOF-haloes and AHF-haloes is done in Section 5. Here we concentrate on the three vital factors in our framework: local geometrical indicators $\lambda_{i}$’s, the softening scale of the field and multistream thresholds.

The maximum values of $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ in each of the haloes have peaks away from 0 as shown in Figure 4. The median values of $max(\lambda_{1})$ and $max(\lambda_{2})$ are in the range of 1-10 (Table 1), in-spite of the threshold for $\lambda_{3}$ being barely positive, by definition. Hence the interior of the potential halo segments is quite convex, with a local maxima inside. In some haloes, the local maxima of eigenvalue are in the order of thousands, as tabulated in Table 1.

With this algorithm, we obtain prospective dark matter haloes - regions with a local maximum of the multistream field in the interior of their closed convex surfaces. The haloes are detected without using density fields or linking lengths between particles. The parameters in the algorithm are entirely based on features of the multistream field and local geometry using Hessian matrices.

4.3 Effect of smoothing

In order to reduce noise, the field is smoothed for our analysis using a Gaussian filter. The effect of smoothing scale on the distribution of the eigenvalue $\lambda_{3}$ in the simulation of $128^{3}$ particles is shown in Figure 5. Effect of softening on the multistream fields does not alter the distribution of multistream distribution significantly (Seen in Figure 10 of Ramachandra & Shandarin 2017). However, the second order variation (and consequently the Hessian eigenvalues) is significantly changed due to the softening of the edges of structures. PDF of $\lambda_{3}$ at multistream smoothing scale of the half the side length of diagnostic grid , $0.5\times l_{d}=0.20h^{-1}\text{ Mpc}$ is noisier than in the scales of $l_{d}$ and $2\times l_{d}$. However, at every scale, the PDF peaks at [math]. The volume fraction of regions with $\lambda_{3}>0$ (i.e. with positive curvature) is $2.4\%$, $2.3\%$ and $2.5\%$ for scales $0.20h^{-1}\text{ Mpc}$, $0.39h^{-1}\text{ Mpc}$, $0.78h^{-1}\text{ Mpc}$ respectively. For the detection of haloes in Section 4, we only look at these regions.

In addition to reducing the numerical noise, smoothing of the multistream field also results in softening of the sub-structures Ramachandra & Shandarin (2017). Since our framework of detecting haloes isolates the multistream regions with local maxima, the closed curvatures are resolved separately. The halo or a sub-halo regions, that enclose the local maxima of the $n_{str}$ field, vary with the smoothing scale of the multistream field. By increasing smoothing of the multistream field, the number of haloes are reduced as shown in Table 2. In the simulation with $256^{3}$ particles, $27929$ $\lambda_{3}$-haloes are detected at smoothing scale equal to the diagnostic grid length, $l_{d}=0.20h^{-1}\text{Mpc}$. The number of haloes decreases to $18221$ and $7897$ at softening scales of two- and four times the grid lengths respectively.

Moreover, since the spatial resolution is higher at the low softening, more small haloes are detected, as shown in lower mass regime of halo mass functions in Figure 6. The tail of halo mass functions reveal that large haloes are more massive for higher softening scales. For instance, the largest haloes for the same simulation with multistream softening length of $0.20h^{-1}\text{Mpc}$, $0.39h^{-1}\text{Mpc}$ and $0.78h^{-1}\text{Mpc}$ have $30650$, $38333$ and $56257$ particles respectively.

The sub-halo finder methods (see Onions et al. 2012 and references therein) identify substructures within a large host halo. The sub-haloes are resolved individually as $\lambda_{3}$-haloes at different scales from our algorithm if the local maxima of the smoothed multistream field is enclosed within the boundary.

4.4 Effect of multistream thresholds

Environmental dependence of the haloes results in various multistream values for the halo core. Theoretical toy models of halo formation, such as the tetrahedral collapse model Neyrinck (2016) describes a three-dimensional halo with four filaments accreting mass into it, has 15 stream crossings. Ramachandra & Shandarin (2015) have previously showed that a high threshold of $n_{str}\geq 90$ is equivalent of virial density of $\rho_{vir}=200$, and filters most of the large haloes above $10^{13}M_{\sun}$.

The algorithm used for detecting multistream haloes initially detects all the closed regions in the multistream ($n_{str}>1$) regions of the cosmological simulation. In order to exclude some of the obvious non-halo sites, we impose a lower threshold of $n_{str}\geq 7$ on the multistream maximum (these regions were also seen as parts of walls or filaments in Ramachandra & Shandarin 2015 ), so that all the sites with three or more foldings in the Lagrangian sub-manifold are chosen. Combining this with the conditions on local eigenvalues, number of particles in haloes etc, we got a pretty good correspondence with other halo finders as demonstrated in section 6 .

Although this condition is by no means strict, it is necessary to check the validity of the assumption. Figure 7 shows the halo mass functions for the haloes detected with changing thresholds on the multistream values of the halo cores. The figure demonstrates that increasing the cut-off from $n_{str}\geq 3$ to $n_{str}\geq 25$ systematically excludes small mass haloes while the number of haloes with $M\gtrsim 2\times 10^{12}M_{\bigodot}$ remains the same.

5 Halo properties

Multistream environment of haloes can be very diverse. Ramachandra & Shandarin (2015) demonstrated that a majority of the FOF-haloes are in contact with the single-streaming voids. Illustration in Ramachandra & Shandarin (2017) also shows that a large number of FOF-haloes have more than 10 per cent void on the spherical surfaces with virial radii. The $\lambda_{3}$ haloes are significantly different: none of the $\lambda_{3}$-haloes are in contact with the regions where gravitational collapse has not occurred. This is guaranteed by the lower bound of $n_{str}=3$ on all potential halo candidates. Condition on the multistream field within the potential halo sites also ensures that there are collapses along more than one direction, which corresponds to $n_{str}=7$. Hence by definition, for any multistream halo $H_{i}$, highest and the lowest multistream value are $n_{str}^{high}(H_{i})\geq 7$ and $n_{str}^{low}(H_{i})\geq 3$ respectively.

The potential haloes $H_{i}$s selected by eigenvalue condition $\lambda_{3}>0$ have a local maxima of $n_{str}^{high}(H_{i})$ inside their boundaries. For a large number of these $\lambda_{3}$-halo candidates, the maximum $n_{str}^{high}$ is higher than the bound of $n_{str}\geq 7$, as shown in Table 3 and Figure 8. For simulation with $128^{3}$ particles, the median of this peak $n_{str}^{high}(H_{i})$ value is $17$. Unsurprisingly, the global maximum of the multistream field ($n_{str}=2831$) is also a local maximum for one of the haloes. On the other hand, the median of lowest multistream value $n_{str}^{low}(H_{i})$ in the haloes is 3 (Table 3), which is also the first stage of non-linearity.

An important feature of our halo detection method is that the detected $\lambda_{3}$-haloes do not have a global threshold on $n_{str}$ or mass density values. The local conditions may be more suited in identifying haloes in multistream fields, since the multistream environments around haloes are very diverse. For instance, regions with $n_{str}\geq 25$ are tubular around one of the the large haloes in Figure 9. Even the region with more than $75$ streams does not enclose a convex multistream region. Whereas, for $n_{str}\geq 200$ the region is convex and the particles detected by our method reside mostly within. We detect closed regions in the multistream field as long as they are not in void, and have at least three foldings in the Lagrangian sub-manifold.

However, the $\lambda_{3}$-halo boundary is different from any constant multistream contour. That is, from the function $n_{str}(\mn@boldsymbol{x})$, convex regions in the four-dimensional function space $(-n_{str},x,y,z)$ are projected onto three-dimensional co-ordinate space using eigenvalues. This is different from projecting ‘iso-multistream’ slice onto three-dimensional co-ordinate space. Appendix A illustrates the difference in the two approaches for a one-dimensional function.

The multistream field usually has concentric shells in the regions of successive gravitational collapses (as explained in Section 2 and Appendix in Ramachandra & Shandarin 2017). In a specific scenario of Figure 9, regions of lower number-of-streams ($n_{str}=25$ and lower) is are tubular and have regions of higher $n_{str}$ inside ($n_{str}=200$ and higher) that is closed. However, this transition from concavity to convexity of the multistream field does not occur at a constant value of $n_{str}$ throughout the field. Instead, it is a local geometrical change that occurs at $\lambda_{3}=0$. For the $\lambda_{3}$-haloes in our simulation ($N_{p}=128^{3}$), minimum of multistream values $n_{str}^{low}(H_{i})$ within each halo has a range of values shown in Table 3 and Figure 8 – this varies between $3\leq n_{str}^{low}(H_{i})\leq 459$. Hence a global condition on $n_{str}$ does not guarantee that the region is convex.

The particles in a massive $\lambda_{3}$-halo shown in Figure 9 form a spheroidal structure. The number of particles in similar massive haloes are in the order of $10^{3}-10^{4}$ particles. For instance, the most massive halo in the simulation (with $N_{p}=128^{3}$) has $5593$ particles. We have chosen a minimum threshold of 20 particles, which is an artificial parameter (may be cooked up or ad hoc) used by most halo finder methods. Majority of the $\lambda_{3}$-haloes have low number of particles; median of number of particles per halo is $49$. Each particle in this simulation is approximately $3.65\times 10^{10}h^{-1}M_{\sun}$. Hence the halo mass range varies in the order of $10^{11}M_{\sun}$ to $10^{14}M_{\sun}$. Combined mass of all the $\lambda_{3}$-halo candidates is about $31$ per cent of the total mass in the simulation. In contrast, the haloes occupy just $0.3$ per cent of the total volume. Thus the $\lambda_{3}$-haloes are extremely dense structures. Further analysis of halo mass function of $\lambda_{3}$-haloes and comparison with AHF- and FOF-haloes is done in the Section 6.

6 Correspondence with other halo finders

Comparison of haloes obtained from AHF and FOF method, along with our geometric analysis of the multistream field reveals several interesting features. The number of haloes ($N_{H}$) with at least 20 particles that were detected by all the algorithms is shown in Table 4. For both the simulations, FOF detects the highest number of haloes and AHF detects the least. By applying the Hessian algorithm on multistream fields smoothed at the scale of diagnostic grid size, $l_{d}$, we detected around $4500$ and $28000$ haloes in simulations with $128^{3}$ and $256^{3}$ particles respectively. The number of $\lambda_{3}$-halo is close to the mean of AHF- and FOF- haloes in each simulation – i.e., $N_{H}^{\lambda_{3}}$ is around $2$ per cent of mean of $N_{H}^{\rm AHF}$ and $N_{H}^{\rm FOF}$ for the $N_{p}=128^{3}$ simulation and $8$ per cent for the $N_{p}=256^{3}$ simulation. Multistream field both the simulations we calculated at a refinement factor of $l_{l}/l_{d}=2$.

The halo mass functions from all three finders are shown in Figure 10. For smaller haloes of order of $10^{13}M_{\sun}$, our method predicts a slightly higher number of haloes than FOF and AHF. For the most massive haloes of mass around $10^{14}M_{\sun}$, number of $\lambda_{3}$-haloes is fewer than the other 2 methods, albeit around the range of error of AHF-haloes.

By observing some of the massive haloes, like the one in Figure 11, we find that the $\lambda_{3}$-halo particles are within AHF- or FOF-halo region. This is generally observed in other massive haloes too: the large $\lambda_{3}$-haloes have fewer particles than corresponding AHF- or FOF-haloes. For haloes greater than $10^{14}M_{\sun}$, $\lambda_{3}$-haloes have boundaries slightly within the AHF virial radius. Without unbinding, the FOF-haloes can be very large compared to other methods, as seen in Figure 11. This results in a deviation in the $\lambda_{3}$-halo mass function (Figure 10) from the other two methods over halo mass of $10^{14}M_{\sun}$. Further discussion of size of the detected $\lambda_{3}$-haloes in the context of smoothing of the multistream is done in Section 4.3.

The particles identified by the AHF as belonging to haloes form spherical structures due to a series of processes (including virialization) applied to unbind the particles. Inherently, the iso-density contours at virial levels are not spherical or spheroidal. The virialized AHF-haloes on the web are shown in the top panel of Figure 12. However, the spherical AHF-haloes are fewer in number compared to the other methods.

The popular choice of linking length of $b=0.2$, although corresponding to virial density, does not ensure that the haloes have positive curvature. Most algorithms based on the FOF method re-define the halo boundaries by unbinding the particles outside a truncation radius. This truncation radius maybe the distance from the centre of mass of the halo to the farthest particle, rms distance, or an inflection point in the density field (For details on these methods, see Knebe et al. 2011 and references therein). Some halo finders define the virial radius, $r_{vir}$ at the distance from halo center where the density is $\Delta_{vir}$ times the background density. In the middle panel of Figure 12, the FOF-haloes are shown without any of the above post-processing schemes. Without any unbinding, the FOF-haloes are generally larger in size than $\lambda_{3}$-haloes in the centre panel of Figure 12. For a specific case of a massive halo, Figure 11, FOF identifies more particles as bounded, than AHF or our algorithm.

In contrast to the AHF and FOF algorithms, our halo method directly detects a closed, convex surface (approximately the largest one, since $\lambda_{3}>0$) for each of the haloes. There is no unbinding procedure on the particles identified within the halo site. The boundaries of $\lambda_{3}$-haloes themselves are not spherical or of any regular structure, but they are closed convex surfaces, as seen in the middle panel of Figure 12.

Haloes from the three finders in Figure 12 also show other differences in the halo boundaries. While all the AHF-haloes are spherical by definition, the FOF-haloes are irregular. The boundaries of the $\lambda_{3}$ haloes are not spherical either, but are more compact than FOF-haloes and in addition they are convex surfaces by design. At some junctions of the filaments, FOF identifies a large region as belonging to single halo, whereas AHF and our algorithm detect multiple isolated haloes. Each isolated $\lambda_{3}$-halo region enclose one maximum of multistream field, thus guaranteeing that multiple haloes are always resolved. On the other hand, a linking length cut-off or a constant threshold on scalar fields may enclose regions with multiple local maxima ( For one-dimensional fields, an illustration of this is shown in Appendix A).

For a simulation box with $N_{p}$ particles, each of mass $m_{p}$, the halo-mass fractions, $\displaystyle f_{h}=(\sum\limits_{i=1}^{N_{H}}m_{H}(i))/(m_{p}N_{p})$ (where $m_{H}$ is the mass of each halo and $N_{H}$ is total number of haloes) are shown in a Venn diagram in Figure 13. For the simulation with $N_{p}=128^{3}$ particles, AHF-, $\lambda_{3}$-, and FOF- haloes occupy $22$, $31$ and $35$ per cent of the total mass respectively. Nearly $19$ per cent of the total mass are concurrently detected as belonging to haloes by all the three algorithms. FOF (with highest halo mass fraction) detects virtually all the haloes that AHF (with least halo mass fraction). About $3$ per cent of the particles classified as belonging to haloes by both AHF and FOF are not classified as multistream halo particles. Our method also detected nearly $6$ per cent of mass particles as haloes, which neither FOF nor AHF classify as haloes. For simulation with $N_{p}=256^{3}$ particles, the corresponding halo mass fractions $f_{h}^{\rm AHF}=30$ per cent, $f_{h}^{\lambda_{3}}=32$ per cent, and $f_{h}^{\rm FOF}=42$ per cent respectively. Thus the mass fraction $f_{h}^{\lambda_{3}}$ remains fairly consistent over increasing mass resolution, as opposed to AHF and FOF. However, large fractions of these mass particles, nearly $23$ per cent of the $N_{p}=256^{3}$ (increased from $19$ per cent for low mass resolution simulation), are simultaneously detected as belonging to haloes by different methods, as shown in the right panel of Figure 13. For the simulation with $N_{p}=256^{3}$, we also see increase in agreement between any two pairs for halo finders. That is, the mass fraction of haloes simultaneously detected by $\lambda_{3}$ and FOF jumps from $25$ per cent to $27$ per cent. This correspondence increases from $19$ to $23$ per cent for $\lambda_{3}$-AHF, and $22$ to $30$ per cent for FOF-AHF pairs. For the same pair (in the $N_{p}=256^{3}$ simulation), $12.3$ per cent of particles are detected by FOF but not AHF, whereas almost all the particles ($>99.9$) for the AHF particles were also detected by FOF. For the pair $\lambda_{3}$-AHF, $9.5$ per cent of particles are detected by $\lambda_{3}$ but not AHF, and $7.3$ per cent of particles were detected by AHF but not by $\lambda_{3}$. Finally, for the $\lambda_{3}$-FOF pair, $5$ per cent of particles are detected by $\lambda_{3}$ but not FOF, and $15$ per cent of particles were detected by FOF but not by $\lambda_{3}$.

On the other hand, looking at the mass particles that were only detected as haloes by one method, but not by other two, we see that only the multistream haloes improve (i.e., the disagreement reduces from $5.9$ to $5$ per cent) with mass resolution. FOF detects $6.3$ and $7.8$ per cent of haloes in simulations of $128^{3}$ and $256^{3}$ particles respectively, that were not classified as haloes by the other two methods. AHF-halo particles, being sub-set of FOF-haloes for the most part, show less than $0.1$ per cent disagreement with other finders.

The discrepancies may have to be addressed on a case-by-case basis. One of the major difference between the haloes detected by isolating local multistream maxima regions and AHF/FOF is shown for a large halo in Figure 11. Without any unbinding procedure, FOF may detect very large irregular sized haloes, often consisting of multiple sub-haloes as shown in the bottom panel. On the other hand, the corresponding AHF-halo (top panel) is smaller spherical subset of FOF-halo. Furthermore, the $\lambda_{3}$-halo in the middle panel is smaller than both. Our multistream field detection technique selects convex regions with strictly one $n_{str}$ maxima within them. The sub-haloes detected by FOF (or AHF), may be detected as separate $\lambda_{3}$-haloes. Nevertheless, some of the mass particles between the two neighboring haloes (like ones along saddle regions of multistream fields) will not be included as belonging to the halo. This effect is seen in halo mass functions (Figure 10 for large haloes of mass more than $10^{14}M_{\bigodot}$ – number density of large $\lambda_{3}$-haloes is smaller than FOF. Similarly it causes a few discrepancies in mass fractions of potential haloes as well.

Other cause for differences in mass fraction is also rooted in the definition of haloes. Single-streaming regions are excluded from our halo search completely. Whereas, FOF and AHF employ no such mechanism to check for number of gravitational collapses. Ramachandra & Shandarin (2015) showed that a significant fraction (nearly $35$ per cent) of FOF-haloes are in contact with single-streaming voids. Particles within these regions would not be considered as potential $\lambda_{3}$-halo particles. This also contributes to the discrepancy in halo mass fraction by different halo finders.

6.1 Haloes in the percolating filament

The excursion set of multistreams above an $n_{str}$ threshold hosts a varying number of haloes. We compare the multistream halo candidates from our geometric method with the AHF and the FOF method in Figure 14

for the simulation with $N_{p}=256^{3}$. The regions in the co-ordinate space are classified into excursion set and non-excursion set regions based on whether the multistream is over or under the $n_{str}$ threshold. In the excursion set we also distinguish the largest structure from the rest of the structures because the largest region of the excursion set plays the crucial role in detecting the transition to percolation. Percolation takes place at thresholds $n_{str}\leq 17$ (Ramachandra & Shandarin 2017) to the right from the vertical dashed dashed line. Based on the coordinates of the halo particles, we check if a halo is in contact with the largest region of the excursion set or with rest of the excursion set.

The fraction of haloes in the non-excursion set are shown at various $n_{str}$ thresholds in the top panel of Figure 14. At thresholds greater than 17 streams (i.e. in non-percolating regime), a large fractions of the AHF-, FOF- and $\lambda_{3}$-haloes are in the non-excursion set, as shown in the top panel of Figure 14. The fraction of $\lambda_{3}$-haloes is slightly higher than FOF or AHF in this regime. At relatively high threshold of, say, $n_{str}=35$, about 65% of the AHF-haloes, about 75% of the FOF-haloes and about 80% of the $\lambda_{3}$-haloes are in the non-excursion set.

For multistream values lower than the percolation threshold of $n_{str}\leq 17$ (i.e. in percolating regime) the fractions of AHF-, $\lambda_{3}$- and FOF-haloes in the largest (i.e. percolating) region quickly grow with decreasing values of the threshold and surpass both the corresponding fractions in the non-percolating regions of the excursion set and that in the non-excursion set at $n_{str}\approx 10$. The majority of the haloes belong to the single percolating structure (shown for the simulation with $N_{p}=128^{3}$ in Figure 12. Similar spatial distribution of SUBFIND haloes Springel et al. 2001b in the multistream regions is shown by Aragon-Calvo et al. 2016) and at $n_{str}=3$, all the haloes are attached to the web.

At $n_{str}=3$, the filling fraction $f_{1}/f_{ES}$ is almost close to unity Ramachandra & Shandarin (2017). Most halo candidates from all three algorithms are at least in contact with the percolating cosmic web. Due to the threshold on $n_{str}$ in our halo detection method, the $\lambda_{3}$-haloes are not only in contact, but completely within this structure.

7 Discussion

The multistream field can have only integral values, moreover these values must be odd numbers practically at every point, because the number of streams can be even only on a set of measure zero. It means that in numerical simulations even values may occur on extremely rare occasions. We have analysed functional variation of the scalar field $-n_{str}(\mn@boldsymbol{x})$ using Hessian eigenvalues. The Hessian analysis is generally done for inherently continuous fields since it requires evaluation of the second derivatives. Geometries of web structures unveiled by Hessian signatures of smoothed density fields (such as Sousbie et al. 2008, Aragon-Calvo et al. 2007, Aragon-Calvo et al. 2010, Cautun et al. 2014, Bond et al. 2010a and many more), tidal shear or velocity shear tensor (Hahn et al. 2007, Forero-Romero et al. 2009, Hoffman et al. 2012, Hoffman et al. 2012, Libeskind et al. 2013, Cautun et al. 2014 etc.) and observational data from surveys (Sousbie et al. 2008, Bond et al. 2010a, Bond et al. 2010b, Pahwa et al. 2016 etc.).

Although the multistream field has discrete values by definition, it may be smoothed for numerical analysis at some scale, typically the scale of grid length of the field. The resulting Hessian eigenvalues approximately characterize the geometry in a four-dimensional hyper-space of ($-n_{str},x,y,z$). Our only assumption about the shape of the boundary of a virialized halo is that it is a convex surface. Therefore the boundary of a halo can be defined as a region with $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}>0$ since it is a closed convex contour in the ($-n_{str},x,y,z$) hyper-space, and thus it’s projection onto the three-dimensional $(x,y,z)$ space is also closed and convex.

Differentiating a smoothed $n_{str}(\mn@boldsymbol{x})$-field may still pose a problem in the regions where $n_{str}(\mn@boldsymbol{x})$ is close to a constant and therefore the eigenvalues represent noise about the zero level. Fortunately in the most of such regions the unsmoothed $n_{str}(\mn@boldsymbol{x})=1$ therefore they can be easily eliminated.

Our algorithm for detecting potential dark matter haloes is unique due to two important factors: the geometrical attribute and the choice of field. Local geometrical analysis on the multistream field conveniently delineates the non-void structures without any free parameters. The dark halo candidates have compact surfaces that enclose local maxima of the multistream field. We do not employ non-local thresholds that several halo finders use (see Knebe et al. 2011, Knebe et al. 2013 and Onions et al. 2012 for comparisons of various halo finders). Global thresholds (like a constant $n_{str}$ cut-off) might be unsuitable for detecting halo candidates since the halo multistream environments are generally not categorical. Secondly, the $n_{str}$ field enables us to mask out the regions belonging to mono-stream regions without a heuristic criteria. Our method guarantees that none of the $\lambda_{3}$-halo particles belong to void region.

We note that present halo finders employ a variety of physical and numerical processes to identify dark matter halo candidates. Furthermore, there is no consensus in the definition of haloes itself (see discussion in Knebe et al. 2011). This is also the cause for the few differences between FOF-, AHF- and multistream haloes: FOF and AHF haloes only use Eulerian co-ordinates $\mathbf{x}(z)$ – either raw positions or in the form of mass density. On the other hand, we utilize a mapping on the Lagrangian sub-manifold $\mathbf{x}(\mathbf{q},z)$ to define the multistream field $n_{str}(\mathbf{x})$. The boundaries of haloes in FOF and AHF are defined by the free-parameter thresholds of linking length and density, and the halo-center is usually defined as the center of mass of the particles. Conceptually, the center of a $\lambda_{3}$-halo is the location of the local multistream maximum, and the boundary of the halo is the convex region surrounding it. At least for large haloes like the one in Figure 11, this convex region is well within the FOF boundary.

The applicability of non-local thresholds in detecting haloes deserves deeper investigation. Lower bounds on over-density or linking-length thresholds traditionally define halo regions in several halo finders. Values such as $\Delta_{vir}\approx 180$ or $b\approx 0.2$ correspond to virial theorem applied to isolated spherical collapse models. Recently More et al. (2011) demonstrated that depending on halo environment, cosmology and redshift the over-densities corresponding to $b=0.2$ have different values. The virial theorem itself is a good measure of equilibrium of a system. However, the global thresholds empirically derived from it may not be pertinent to diverse environment of dark matter haloes.

The algorithm prescribed in Section 4 lists out a set of physically motivated steps that filter out the noisy $\lambda_{3}>0$ regions that cannot be identified as haloes. The analysis in the simulation of $100h^{-1}$ Mpc side length and $128^{3}$ particles, with the multistream calculated on $256^{3}$ diagnosis box, approximately 40000 labelled segments satisfying $\lambda_{3}>0$ criterion in the non-voids were filtered out by a lower limit on multistreaming regions. One of the possible improvements in our algorithm would be to use information of number of flip-flops of each particle (For instance, using methods prescribed by Shandarin & Medvedev 2014, 2017). Such methods involving the Lagrangian sub-manifold may reveal rich sub-structure in the haloes.

The requirement that each halo should have closed convex surfaces with a multistream maximum inside may identify sub-haloes in large haloes but is too demanding because a halo with sub-haloes must have saddle points in the $n_{str}$ field. This may explain the shortage of massive haloes shown in Figure 10. Applying a filter for smoothing the $n_{str}$ field increases the number of massive haloes but reduces the number of low mass haloes. Although our present method does not currently perform an analysis simultaneously on multiple smoothing scales, such approaches done in density, lognormal density, tidal, velocity divergence or velocity shear fields (see MMF by Aragon-Calvo et al. 2007 and NEXUS+ by Cautun et al. 2013) have shown interesting multi-scale features of the cosmic web. Applying a more sophisticated procedure for linking $\lambda_{3}$-sub-haloes in a more massive halo will be done in the follow-up paper.

Dark matter haloes, being localised structures, are uniquely convenient for our local Hessian analysis. Conditions of $\lambda_{1}>0>\lambda_{2}\geq\lambda_{3}$ and $\lambda_{1}\geq\lambda_{2}>0>\lambda_{3}$ also give information about curvature. Hessian eigenvalue analysis at high resolution of multistream fields may be very interesting in understanding the tubular edges of filaments and surfaces of walls at smaller scales. However, in this study, Hessian analysis is only applied to haloes. Walls and filaments span large volumes in the dark matter simulations, and we employ topological tools to investigate them.

8 Summary

We studied certain geometrical of the multistream field in the context of halo formation. Findings from our analysis are summarized as follows:

Several aspects of halo formation in the Lagrangian sub-manifold are considerably different than that of reference models of spherical top-hat collapse and ellipsoidal collapse. Successive formations of caustics (and consequently multiple velocity streams) play a crucial role in the process of clustering. 2. 2.

We present a novel halo detection algorithm for identifying dark matter halo candidates in the multistream field. Conditions on the local geometric indicators of the field are used to ensure that each closed halo boundary hosts a local multistream maxima. The positive signs of all principal curvatures (please note that we use curvatures of $-n_{str}(\mn@boldsymbol{x})$ field) inside the boundary also guarantee that the boundary is convex. Bounds on $n_{str}$ guarantee that all the halo particles are in the non-void structure. We also ensure that the halo regions have foldings in the Lagrangian sub-manifold in more than one direction. 3. 3.

The multistream field within the halo boundaries may be very diverse. We do not detect halo candidates from a global lower bound on $n_{str}$. Instead, we look for closed convex regions in the multistream field. For the simulation with $128^{3}$ particles, minima of $n_{str}$ in each halo vary from 3 to nearly 450. Maxima of $n_{str}$ in the halo vary from 7 to about 2800. 4. 4.

Our multistream halo candidates had a reasonably good correspondence with haloes from AHF and FOF catalogues. One notable difference was found with massive haloes. Our algorithm predicted fewer particles than the FOF method. This is likely to be caused by the requirement that the multistream field in the regions of the $\lambda_{3}$-halo candidates is convex which may be a reasonable approximation for simple haloes (i.e. having no sub-haloes) but massive haloes are more likely to have sub-haloes and therefore the $n_{str}$ field in the corresponding regions must have saddle points ant therefore cannot be entirely convex. Our study of the smoothing effects has shown that the number of massive sub-haloes tend to increase with growing smoothing scale which seems to agree with the above explanation. We will address this problem in the following study. 5. 5.

Halo candidates were mostly embedded on the excursion set of the multistream field after percolation transition ($n_{str}=17$ in the simulation with $256^{3}$ particles). At lower thresholds (around $n_{str}=5\text{ to }11$), the largest percolating structure in the excursion set hosts most of the haloes.

In conclusion, the Lagrangian sub-manifold contains dynamical information of structure formation. We analysed the multistream field that contains the information of foldings in the sub-manifold. In addition, we demonstrated the use of geometrical features of the multistream field in identifying potential dark matter halo candidates in cosmological N-body simulations.

Acknowledgements

This work has been funded in part by the University of Kansas FY 2017 Competition General Research Fund, GRF Award 2301155. This research used resources of the Advanced Computing Facility, which is part of the Center for Research Computing at the University of Kansas. We thank Mark Neyrinck and Mikhail Medvedev for discussions and suggestions. We also thank the anonymous referee for insightful comments on improving this manuscript and the editor for help with the submission.

Appendix A Hessian signatures of the multistream field

Second-order local variations of a scalar field $f$ is described by a Hessian matrix, whose element in a three-dimensional domain is given by Equation 1. The geometry of the scalar field is classified by the Eigenvalues of the Hessian. The convex regions have at-most one maxima within the (3+1)-dimensional functional space. Projection of this closed region onto three-dimensional coordinate space also gives a closed surface in coordinate space. An illustration of the projection is shown in Figure 15 for a simpler function $f(x)$ in one-dimensional domain. The eigenvalue criteria for regions are simplified: for instance, $\frac{\partial^{2}f}{\partial x^{2}}<0$ for convex region. Projection of these regions onto coordinate space is shown in the shaded regions. This is different from regions within a contour, which is the projection of the curve along which the function has a constant value. Boundaries of these two regions may, but not necessarily, intersect.

In the case of cosmic fields, thresholds like $\Delta_{vir}$ are equivalent to the green dotted line in Figure 15. The over-dense regions (green shaded regions) are not constrained to be convex. Similarly structures selected based on $n_{str}$ thresholds do not universally result in convex structures either. Local geometry can be probed from the eigenvalue criteria instead, as shown by the red line on the curve and corresponding shaded area. The projected structures, albeit convex, may have very small values of $f(x)$ (like the red shaded area around $x=5$). In the framework of identifying potential haloes in multistream field, multistream thresholds are devised in so that some of these small peaks detected by the Hessian are not considered as potential halo sites.