# The Hierarchical Cosmic Web and Assembly Bias

J. M. Coloma-Nadal  F.-S. Kitaura  J. E. García-Farieta  F. Sinigaglia  G.
Favole  and D. Forero Sánchez

###### Abstract

Accurate modeling of galaxy distributions is paramount for cosmological
analysis using galaxy redshift surveys. However, this endeavor is often
hindered by the computational complexity of resolving the dark matter halos
that host these galaxies. To address this challenge, we propose the
development of effective assembly bias models down to small scales, i.e.,
going beyond the local density dependence capturing non-local cosmic
evolution. We introduce a hierarchical cosmic web classification that
indirectly captures up to third-order long- and short-range non-local bias
terms. This classification system also enables us to maintain positive
definite parametric bias expansions. Specifically, we subdivide the
traditional cosmic web classification, which is based on the eigenvalues of
the tidal field tensor derived, with an additional classification based on the
Hessian matrix of the negative density contrast. We obtain the large-scale
dark matter field on a mesh with $\sim 3.9\,h^{-1}$ Mpc cell side resolution
through Augmented Lagrangian Perturbation Theory. To assess the effectiveness
of our model, we conduct tests using a reference halo catalogue extracted from
the UNIT project simulation, which was run within a cubical volume of 1
$h^{-1}$ Gpc side. The resulting mock halo catalogs, generated trough our
approach, exhibit a high level of accuracy in terms of the one-, two- and
three-point statistics. They reproduce the reference power spectrum within
better than 2 percent accuracy up to wavenumbers $k\sim 0.8\,h$ Mpc-1 and
provide accurate bispectra within the scales that are crucial for cosmological
analysis. This effective bias approach provides a forward model appropriate
for field-level cosmological inference and holds significant potential for
facilitating cosmological analysis of galaxy redshift surveys, particularly in
the context of projects such as DESI, EUCLID, and LSST.

##  1 Introduction

The primary objective of cutting-edge galaxy surveys like DESI [1] or EUCLID
[2], is to meticulously chart the vast cosmic web, revealing the intricate
three-dimensional distribution of matter. This endeavor serves the fundamental
purpose of probing the elusive nature of dark energy and its influence on the
Universe’s expansion. Achieving this demands a comprehensive understanding of
cosmological parameters, the evolution of cosmic structures, the historical
trajectory of cosmic expansion, and even theories that challenge conventional
notions of gravity. To accomplish these feats, these surveys rely heavily on
theoretical frameworks and sophisticated galaxy mock catalogues, which are
essential for conducting robust analyses and estimating covariance matrices.

Generating mock catalogues necessitates the faithful recreation of cosmic
volumes observed in actual galaxy surveys, along with the replication of
clustering observables. Traditionally, simulating structure formation has
relied on computationally intensive $N$-body codes. However, with the immense
volumes involved in contemporary surveys, the computational burden becomes
prohibitive. As a pragmatic alternative, approximate gravity solvers come into
play. Among these, some leverage Lagrangian Perturbation Theory (LPT) [3] to
compute the dark matter density field. Noteworthy examples of these approaches
include in chronological order: PEAKPATCH [4, 5, 6, 7], PINOCCHIO [8, 9],
PThalos [10, 11, 12], PATCHY [13, 14, 15, 16], EZmocks [17, 18], HALOGEN [19],
ICE-COLA [20], and BAM [21, 22, 23, 24], each offering distinct methodologies
for generating galaxy mock catalogues (see comparison projects, [25, 26, 27,
28]).

These methods can be classified into two categories: the ones trying to
reproduce the halo distribution based on analytical and numerical
approximations (PEAKPATCH, PINOCCHIO, later versions of PThalos, HALOGEN, ICE-
COLA); and those trying to reproduce the distribution of large-scale structure
tracers from detailed reference calculations with effective field-level bias
models (e.g., the earliest PThalos version, PATCHY, EZmocks, and BAM).

The latter category offers substantial savings in computational resources, as
demonstrated in a comparative study by [27]. Moreover, it enables a detailed
dissection of contributions from individual bias components, as elucidated in
works by [14, 23]. Importantly, this methodology can be seamlessly extended to
field-level inference studies (see, e.g.,[29, 30]).

In this study, our approach to modelling structure formation hinges on
augmented Lagrangian Perturbation Theory (ALPT) [31], an advancement over
linear LPT. ALPT refines the displacement field by decomposing it into long-
and short-range components, drawing from second-order LPT (2LPT) and the
spherical collapse approximation, respectively. This nuanced treatment allows
for a more precise characterisation of the evolving cosmic structures within a
one-step gravity solver. We could, however, rely on more accurate and
expensive approaches (e.g., eALPT[32], FASTPM [33], and COLA [34]). The
findings of the present study will yield to improved results when using more
sophisticated gravity solvers, as demonstrated in [16].

The advantage of PATCHY over EZmocks and BAM is its explicit, clear analytic
bias description, i.e., the description of the intricate relationship between
the dark matter density field and its tracers such as halos, galaxies, galaxy
clusters, etc. However, PATCHY faced challenges in accurately modelling the
higher-order statistics of low-mass tracers, particularly halos with masses
above approximately $10^{8}M_{\odot}\,h^{-1}$ in a cubical volume of 400
$h^{-1}$ Mpc side, despite incorporating second-order non-local bias effects
[35].

The study by [23], using high mass halos (halos with $\gtrsim
10^{12}M_{\odot}\,h^{-1}$ in a 1,5 $h^{-1}$ Gpc side cubical volume), showed
that a proper modelling of the higher order statistics requires going to at
least third-order non-local bias, and that this can be done through the cosmic
web classification. These findings were exploited in the context of effective
baryonic physics bias mapping and Lyman-$\alpha$ forest modelling (see [36,
37, 38]).

This paper extends previous PATCHY versions with a novel hierarchical cosmic
web classification, indirectly accounting for both long- and short-range non-
local bias at least up to third order. This bias description turns out to be
crucial in acquiring a high fidelity forward model, as tested on a $N$-body
based halo catalogue with a mass cut of $\gtrsim 10^{11}M_{\odot}\,h^{-1}$ in
a 1 $h^{-1}$ Gpc side cubical volume with variance suppression [39].

Our method offers the advantage of streamlining the parametric bias model,
consolidating it into a maximum of six parameters that encompass both local
deterministic and stochastic bias components. Moreover, the indirect
incorporation of non-local bias terms through the hierarchical cosmic web
classification enhances the model’s comprehensiveness and efficiency. Our
approach avoids the explicit inclusion of non-local bias terms, a practice
that can lead to the introduction of negative densities when truncated to a
specific order. Instead, we segment a large-scale structure catalogue into
well-defined subsets, akin to dividing a halo catalogue into mass bins. This
segmentation allows us to treat each subset as a distinct catalogue,
simplifying the application of a straightforward local bias model to
accurately reproduce its characteristics. In this way, we manage to have a
solid description of the large scale bias on large scales as stated by
renormalised perturbation theory [40, 41], and extend it towards small scales
with short-range non-local bias terms. The bias model presented in this study
could potentially be combined with effective field theory based approaches
[42, 43, 44, 45, 46, 47, 48, 49] within field-level frameworks (see, e.g.,
[50]).

While our approach concentrates solely on number counts within a real-space
mesh, it is worth noting that a companion paper introduces a sub-grid model
designed to allocate positions and velocities within the cells (see the
CosmoMIA method [51]). The bias method presented in this paper together with
CosmoMIA are being applied to a variety of large-scale structure tracers, such
as Luminous Red Galaxies (LRGs), Emission Line Galaxies (ELGs), quasars (QSOs)
and Bright Galaxies (BGs), and will be presented in forthcoming publications.

The paper is structured as follows. In section 3, we briefly describe the
structure formation theory we use to evolve the density field. Subsequently,
we present the bias of dark matter tracers trough perturbation theory (PT) and
define our analytic model. At the end of this section we hierarchically
classify the cosmic web. In section 4, we present and statistically analyse
the test on a $N$-body reference catalog. Finally, we discuss the results in
section 5.

##  2 Theory

In this section, we provide an overview of the analytical expression of the
displacement field necessary to calculate the density field. Then we recap the
theoretical background in modelling the bias of large scale structure tracers.

###  2.1 Structure formation

On the large scales, we rely on Lagrangian perturbation theory (LPT) to
simulate structure formation [3]. This is an approximate solution to evolve
particles from the primordial density fluctuations until the desired epoch
(redshift). In this theory, the action of gravity across cosmic evolution is
encoded in the displacement field $\mathbf{\Psi}(\mathbf{q},z)$, which relates
initial Lagrangian positions $\mathbf{q}$ to final Eulerian comoving ones
$\mathbf{x}(z)$ through the following mapping relation

$\mathbf{x}(z)=\mathbf{q}+\mathbf{\Psi}(\mathbf{q},z).$ (2.1)

However, it is well known in the literature that LPT approaches suffer from
artificial shell-crossing, as the particles do not interact with each other
[52, 53, 54, 55]. This problem can be mitigated with Augmented LPT (ALPT)
[31], a combination of LPT with the spherical collapse (SC) model [56, 57,
58], which models virialisation towards small scales. In particular, we rely
on second-order LPT (2LPT) [52, 53], an improvement of the linear
approximation with non-local tidal field corrections. The ALPT displacement
field is then written as sum of two components

$\mathbf{\Psi}(\mathbf{q},z)=\underbrace{\overbrace{\mathcal{K}(\mathbf{q},r_{s})}^{\text{long-
range kernel}}\circ\mathbf{\Psi}_{\mathrm{2LPT}}(\mathbf{q},z)}_{\text{long-
range
interaction}}+\underbrace{\overbrace{\left(1-\mathcal{K}(\mathbf{q},r_{s})\right)}^{\text{short-
range kernel}}\circ\mathbf{\Psi}_{\mathrm{SC}}(\mathbf{q},z)}_{\text{short-
range interaction}}.$ (2.2)

The long-range component is described by the linear growth factor $D$ and the
first- and second-order gravitational potential terms $\phi^{(1)}$,
$\phi^{(2)}$:

$\mathbf{\Psi}_{\mathrm{2LPT}}=-D\nabla\phi^{(1)}+D_{2}\nabla\phi^{(2)}.$
(2.3)

The short-range component reads

$\mathbf{\Psi}_{\rm
SC}=\nabla\nabla^{-2}\left[3\left(\left(1-\frac{2}{3}D\delta^{(1)}\right)^{1/2}-1\right)\right].$
(2.4)

For more details, we refer the reader to [31]. The smooth dark matter density
field $\delta^{\rm ALPT}$ is finally obtained by applying Cloud-In-Cell (CIC)
interpolation to the final positions $\mathbf{x}(z)$ including tetrahedra-
tesselation [59]. The ALPT method presents a better cross-correlation with
$N$-body simulations than applying separately either LPT solutions until third
order or the spherical collapse model alone [31, 60]. Furthermore, compared to
$N$-body solutions, it is considerably more efficient, 111See table 1 in [27]:
comparison between ICE-COLA and PATCHY to reproduce the Minerva simulations
[61]: 33 hours CPU time with 340 Gb memory allocation vs 0.2 hours and 15 Gb,
respectively. while reproducing the power on the very large scales at
different redshifts.

Our studies have shown that the usage of improved gravity solvers makes only
sense when the resolution of the mesh is higher than about $2.5\,h^{-1}$ Mpc.
Since we want to explore at first, the maximal efficiency of the method, we
explore in this paper only resolutions of $\sim 3.9\,h^{-1}$ Mpc.
Conveniently, we opt for a resolution that is a multiple of two lower
($2^{4}$) than the original dark matter field, which serves as our reference
simulation from the UNIT project [39]. We leave an exploration of our
framework with higher resolutions based on an improved version of ALPT, called
eALPT [32] for future work.

###  2.2 The bias of dark matter tracers

Galaxies, galaxy clusters and halos are biased tracers of the dark matter
distribution (for review, see [62]). The relations between different tracers
and the underlying dark matter field, named bias relations, are well known to
be non-linear, non-local (hence, also scale-dependent) and stochastic (see
[13] and references therein). Bias relations also depend on the history
evolution of the Universe since the tracer formed, the so-called assembly bias
(see, e.g., [63, 64, 65]). Assembly bias is also studied through the
dependence of the tracer on the cosmic web (see, e.g., [66, 67]).

On the large scales, where density fluctuations are small, perturbation theory
(PT) can be used to expand the relation between the tracer $\delta_{\rm tr}$
and the dark matter $\delta$ density fields [68],

$\delta_{\rm
tr}(\mathbf{x})=f(\delta(\mathbf{x}))=\sum_{i=0}^{\infty}\frac{b_{i}}{i!}\delta^{i}(\mathbf{x}),$
(2.5)

where $b_{i}$ are the unknown bias parameters. This expansion, however,
presents only a dependence on the local density. In order to allow for history
dependence and long-range non-locality, one can also include a dependency on
the peculiar velocity through $\partial_{i}v_{j}$ and on the gravitational
potential through the tidal field tensor
$\mathcal{T}_{ij}\equiv\partial_{i}\partial_{j}\phi(\mathbf{x})$ [40, 23]. In
this work, we extend local bias models including long- and short-range non-
local bias components through a cosmic web classification. The short-range
non-local bias is based on the second-order derivatives of the density field,
$\Gamma_{ij}\equiv\partial_{i}\partial_{j}\delta(\mathbf{x})$. The
dependencies we investigate can be summarised as follows

$\delta_{\rm tr}\curvearrowleft P\left(\delta_{\rm
tr}|\nabla^{2}\phi(\mathbf{x}),\partial_{i}\partial_{j}\phi(\mathbf{x}),\partial_{i}\partial_{j}\delta(\mathbf{x}),\epsilon\right)\equiv
P\left(\delta_{\rm tr}|\delta,\mathcal{T},\Gamma,\epsilon\right),$ (2.6)

where $\epsilon$ describes stochastic uncertainty arising from the fact that
the tracer distribution is a discrete sample of the underlying dark matter
distribution. By convenience, following [40], the tidal field tensor can be
rewritten by removing its trace

$s_{ij}(\mathbf{x})\equiv\mathcal{T}_{ij}-\frac{1}{3}\delta_{ij}^{K}\delta(\mathbf{x})\,.$
(2.7)

Furthermore, the cosmological principle (homogeneity and isotropy) constrains
the bias parameters to be constant scalar quantities. Then, one should
construct scalars from $s_{ij}$. The first order scalars $s_{ii}$ are zero by
definition and the second and third order terms can be constructed by the
tensor contractions $s^{2}\equiv s_{ij}s_{ji}$ and $s^{3}\equiv
s_{ij}s_{jk}s_{ki}$, respectively. Analogously, we can define a traceless
tensor for the short-range non-local bias term,

$\mu_{ij}(\mathbf{x})\equiv\Gamma_{ij}-\frac{1}{3}\delta_{ij}^{K}\nabla^{2}\delta(\mathbf{x}),$
(2.8)

and its scalars $\mu^{2}\equiv\mu_{ij}\mu_{ji}$ and
$\mu^{3}\equiv\mu_{ij}\mu_{jk}\mu_{ki}$.

The resulting Taylor expansion of the bias model, in Eulerian coordinates
$\mathbf{x}$, has now local ($B_{\rm L}$), non-local long-range ($B^{\rm
LR}_{\rm NL}$), non-local short-range ($B^{\rm SR}_{\rm NL}$), velocity
($B^{v}_{\rm NL}$) and stochastic ($B^{\epsilon}$) bias contributions.
Following [23] and [40], the expansion can be written as

$\footnotesize\begin{split}\delta_{\rm tr}(\mathbf{x})=&B(\delta_{\rm
tr}|\delta)=\overbrace{B_{\rm
L}(\delta)}^{\text{local}}+\overbrace{\overbrace{B_{\rm NL}^{\rm
LR}(\delta)}^{\text{long-range}}+\overbrace{B_{\rm NL}^{\rm
SR}(\delta)}^{\text{short-range}}+\overbrace{B_{\rm
NL}^{v}(\delta)}^{\theta\neq\delta}}^{\text{non-
local}}+\overbrace{\overbrace{B^{\epsilon}(\delta)}^{\text{stochastic
noise}}}^{\text{local\;\&\;non-local}}\\\
=&\underbrace{\underbrace{\overbrace{c_{\delta}\delta(\mathbf{x})}^{\mathrm{local}}}_{\mathrm{first\;order}}+\underbrace{\overbrace{\frac{1}{2}c_{\delta^{2}}\left(\delta^{2}(\mathbf{x})-\langle\delta^{2}\rangle\right)}^{\mathrm{local}}+\overbrace{\frac{1}{2}c_{s^{2}}\left(s^{2}(\mathbf{x})-\frac{2}{3}\langle\delta^{2}\rangle\right)}^{\mathrm{non-
local\;long-
range}}}_{\mathrm{second\;order}}+\underbrace{\overbrace{\frac{1}{3!}c_{\delta^{3}}\delta^{3}(\mathbf{x})}^{\mathrm{local}}+\overbrace{\frac{1}{2}c_{\delta
s^{2}}\delta(\mathbf{x})s^{2}(\mathbf{x})+\frac{1}{3!}c_{s^{3}}s^{3}(\mathbf{x})}^{\mathrm{non-
local\;long-range}}}_{\mathrm{third\;order}}}_{\mathrm{curl-
free\;\&\;\theta=\delta\;terms}:\,B_{\rm L}|^{(1\to 3)}+B^{\rm LR}_{\rm
NL}|^{(2\to 3)}}\\\
&+\overbrace{\underbrace{\underbrace{\frac{1}{2}c_{\mu^{2}}\left(\mu^{2}(\mathbf{x})-\langle\mu^{2}(\mathbf{x})\rangle\right)}_{\mathrm{second\;order}}+\underbrace{\frac{1}{2}c_{\delta\mu^{2}}\delta(\vec{x})\mu^{2}(\mathbf{x})+\frac{1}{3!}c_{\mu^{3}}\mu^{3}(\mathbf{x})}_{\mathrm{third\;order}}}_{\mathrm{curl-
free\;\&\;\theta=\delta\;terms}:\,B^{\rm SR}_{\rm NL}|^{(2\to
3)}}}^{\mathrm{non-local\;short-
range}}+\underbrace{\underbrace{\overbrace{\underbrace{\mathcal{O}(F_{v}(\mathbf{x})|^{(3)})}_{\mathrm{third\;order+}}}^{\mathrm{non-
local}}}_{\mathrm{\theta\neq\delta\,terms}}}_{{\cal O}(B^{v}_{\rm
NL}|^{(3)})}+\underbrace{\underbrace{\overbrace{\underbrace{\mathcal{O}(F_{\delta,\mu,\mathcal{T}}(\mathbf{x})|^{(4)})}_{\mathrm{fourth\;order+}}}^{\mathrm{local\;\&\;non-
local}}}_{\mathrm{curl-free\,\theta=\delta\,terms}}}_{{\cal O}(\delta_{\rm
tr}|^{(4)})}+\underbrace{\underbrace{\underbrace{\overbrace{F_{\epsilon}(\epsilon(\mathbf{x}))}^{\mathrm{local\;\&\;non-
local}}}_{\mathrm{first\;order+}}}_{\mathrm{noise\;terms}}}_{{\cal
O}(B^{\epsilon}(\delta)|^{(1)})}\,,\end{split}$ (2.9)

where the $c$s are the bias factors and the superscripts $(1\to 3)$ denote
terms from first to third order.

In the large-scale limit of $k\to 0$, the observable, renormalised, large
scale bias $b_{\delta}$ as the sum of bias-like terms according to [40] is
given by

$b_{\delta}=c_{\delta}+\left(\frac{34}{21}c_{\delta^{2}}+\frac{68}{63}c_{s^{2}}+\frac{1}{2}c_{\delta^{3}}+\frac{1}{3}c_{\delta
s^{2}}-\frac{16}{63}c_{st}\right)\sigma^{2}+\frac{1}{2}c_{\delta\epsilon^{2}}\sigma^{2}_{\epsilon}\,,$
(2.10)

where $\sigma^{2}$ is the variance of the dark matter field at the considered
resolution, $\sigma^{2}_{\epsilon}$ is the corresponding variance of the
noise, and $c_{st}$ relates to a mixed non-local bias term including velocity
contributions for $\theta\neq\delta$.

##  3 Method

Equation (2.9) encapsulates the principal bias terms of the tracer density
field $\delta_{\rm tr}$ within a Taylor expansion. However, employing
truncated versions of an infinite expansion at a practical level is non-
trivial and carries inherent risks.

As depicted in (2.10), accurate reproduction of large-scale bias necessitates
accounting for at least third-order bias contributions. Yet, this does not
offer guidance on transitioning from large to small scales to achieve a
comprehensive field-level bias description.

As elucidated in [23], the contributions of velocity bias primarily impact
small scales and can generally be disregarded for the majority of baryon
acoustic oscillations (BAO) and redshift-space distortions (RSD) studies.
Consequently, we defer an investigation into extending our bias model to
incorporate terms where $\theta\neq\delta$ for future endeavors. Instead, we
focus on terms where information from velocity shear is inherently encoded in
the tidal field tensor.

Here, we adopt the approach proposed in the aforementioned study and
introduce, for the first time in the context of parametric bias models of
discrete tracer distributions, a cosmic web classification to indirectly model
non-local bias terms. We extend this methodology to similarly model short-
range non-local bias contributions for the first time, as elaborated below.

This approach has several advantages, as it permits us to ensure positive
definite tracer density fields, and directly or indirectly include all bias
terms ($B_{\rm L}$, $B_{\rm NL}^{\rm LR}$, $B_{\rm NL}^{\rm SR}$, $B_{\rm
NL}^{v}$ and $B^{\epsilon}$). Our method consists on applying a deterministic
(section 3.1) and stochastic (section 3.2) bias model to each of the different
regions of the cosmic web which is hierarchically classified in terms of the
gravitational potential $\phi$ and the density field $\delta$ (section 3.3).

The deterministic bias, consists firstly on a power law [69, 70, 13] and
secondly on a threshold bias [71, 13, 72], describing the non-linear behaviour
and allowing to suppress low density regions where not enough mass is
available to form halos. We extend this picture with a suppression of density
peaks which some galaxies tend to avoid as a result of halo exclusion (see,
e.g., [73]).

Finally, the connection between the expected number counts and the discrete
realisation of large scale structure tracers is established through a
likelihood, accounting for stochastic biasing [74, 75]. As advanced by [76], a
non-vanishing correlation function introduces deviations from Poissonity [77,
78]. This affects the statistics at scales below the mesh resolution and needs
to be modelled with non-Poisson likelihoods [79, 80, 13, 35]

###  3.1 Local deterministic bias

Following [13], we can generally write the expected number of a certain tracer
$\langle\rho_{\rm tr}\rangle$ in a given volume element d$V$ as function of
the underlying dark matter density $\delta$,

$\langle\rho_{\rm tr}\rangle_{\mathrm{d}V}=f_{\rm tr}(1+B(\delta_{\rm
tr}|\delta)),$ (3.1)

where $B(\delta_{\rm tr}|\delta)$ is the deterministic bias relation and
$f_{\rm tr}$ is the normalisation factor which controls the number density and
can be described as

$f_{\rm tr}=\frac{n_{\rm tr}}{\langle 1+B(\delta_{\rm
tr}|\delta)\rangle_{V}},$ (3.2)

i.e., by requiring the tracer density field to have the number density of the
reference sample $\langle\langle\rho_{\rm
tr}\rangle_{\mathrm{d}V}\rangle_{V}=n_{\rm tr}$.

In particular, we define the bias relation as the product of a power-law,
describing the non-linear behaviour, an exponential damping for high-density
regions ($\epsilon>0$) and a low-density threshold ($\epsilon^{\prime}<0$),

$1+B_{\rm L}(\delta)=\underbrace{(1+\delta)^{\alpha}}_{\text{power-law
bias}}\times\underbrace{\exp{\left(-\left(\frac{1+\delta}{\rho_{\rm\epsilon}}\right)^{\epsilon}\right)}}_{\text{high-
density threshold
bias}}\times\underbrace{\exp{\left(-\left(\frac{1+\delta}{\rho^{\prime}_{\rm\epsilon}}\right)^{\epsilon^{\prime}}\right)}}_{\text{low-
density threshold bias}}\,,$ (3.3)

where we have extended the model presented in [13] to also account for high
density suppression.

###  3.2 Stochastic bias

The distribution of halos on a dark matter density field defined on a mesh,
characterised by a lower resolution than that needed to resolve the halos,
leads to a stochastic process. This can be described by the probability
$P(N_{i}|\lambda_{i})$ of finding a given number of halos $N_{i}$ in a given
cell $i$ with volume $dV$. Therefore, $\lambda_{i}=\langle\rho_{\rm
tr}\rangle_{\rm dV}\times\mathrm{d}V$ is the average number of particles in a
cell. The Poisson distribution,
$P(N_{i}|\lambda_{i})=\frac{\lambda_{i}^{N_{i}}}{N_{i}!}\exp(-\lambda_{i})$,
can account for the discrete nature of galaxies [81, 82]. However, it only
approximately holds in intermediate density regions. In general, we have over-
dispersion. Non-Poissonian distributions, such as the Negative Binomial (NB),
where introduced to account for field-level bias modelling [13, 30, 50]:

$P(N_{i}|\lambda_{i},\beta)=\underbrace{\frac{\lambda_{i}^{N_{i}}}{N_{i}!}\exp{(-\lambda_{i})}}_{\text{Poisson
distribution}}\times\underbrace{\frac{\Gamma(\beta+N_{i})}{\Gamma(\beta)(\beta+\lambda_{i})^{N_{i}}}\frac{\exp(\lambda_{i})}{(1+\lambda_{i}/\beta)^{\beta}}}_{\text{deviation
from Poissonity}}\,.$ (3.4)

The Poisson distribution is recovered in the limit $\beta\rightarrow\infty$
for which the deviation tends to one [16]. One can also model a density
dependent $\beta$-parameter to account for arbitrary variance expressions (see
[35]). In this work, we model the stochastic bias by drawing samples of the
deterministic bias model (3.3) from the NB distribution.

###  3.3 Non-local deterministic bias: The hierarchical cosmic web
classification

Numerical simulations have shown that the gravitational growth of the initial
density perturbations leads to a characteristic pattern of structures, dubbed
cosmic web [83, 84], very similar to the one from observations of the large-
scale distribution of galaxies (see, e.g., [85]). Most of the Universe is
composed of low-density regions (voids) surrounded by filaments of matter
generating a web-pattern. Inspired by the formation of cosmic sheet-like
pancake structures [86], [87] proposed a theoretical criterion to classify the
different regions of the cosmic web in terms of the particle motion inside the
gravitational potential of a dark matter halo, $\ddot{x}=-\nabla\phi$. The
center of mass of a halo is an extremum, i.e. $\nabla\phi=0$, and hence,
there, the equation of motion can be linearised. At this points, the particle
motion is completely determined by the tidal field tensor
$T_{ij}\equiv\partial_{i}\partial_{j}\phi$ and, more concretely, by its
eigenvalues $\lambda_{1}>\lambda_{2}>\lambda_{3}$. Analogously to the
Zel’dovich pancake, we consider regions with three positive eigenvalues as
those where gravitational collapse cause matter inflows. On the other hand, we
consider unstable orbits (expanding regions) when the three eigenvalues are
negative. Attending to this argument, one can define four regions in the
cosmic web

  * •

Knots: $\lambda_{1}$, $\lambda_{3}$, $\lambda_{3}>$ $\lambda_{\rm th}$

  * •

Filaments: $\lambda_{1}$, $\lambda_{2}>$ $\lambda_{\rm th}$ $\&$
$\lambda_{3}<$$\lambda_{\rm th}$

  * •

Sheets: $\lambda_{1}>$ $\lambda_{\rm th}$ $\&$ $\lambda_{2}$, $\lambda_{3}<$
$\lambda_{\rm th}$

  * •

Voids: $\lambda_{1}$, $\lambda_{3}$, $\lambda_{3}<$$\lambda_{\rm th}$

Figure 1: ALPT Dark matter density field and the cosmic web. Slices of 20
$h^{-1}$ Mpc thick and 1000 $h^{-1}$ Mpc side at redshift $z$ = 1.032 are
shown for the logarithmic density field (left panel), the $\phi$-web (middle
panel) and the $\delta$-web (right panel). The darker the color, the highest
the density.

The threshold eigenvalue $\lambda_{\rm th}$, allowing for $\lambda_{\rm
th}\neq 0$, was introduced by [88] to reduce the dependence on the smoothing
scale, mesh resolution and mass assignment scheme. This classification can be
generalised to the full density field rather than limit to the centre of the
halos. As is essentially defined by the gravitational potential, we dubbed
this classification $\phi$-web (see figure 1). It has been used in simulations
both to study the properties of dark matter halos, galaxies and intergalactic
gas [87, 67, 89, 36] and to generate halos and galaxy mock catalogues [90, 21,
35]. Furthermore, it is used in observational analysis [91, 92, 93].

The gravitational potential is a key ingredient in both the bias model (2.9)
and the cosmic web classification. The distribution of halos can be associated
with various regions of the cosmic web by considering the invariants of the
tidal field tensor. These invariants, in turn, directly describe the non-local
bias terms within halo perturbation theory (for details see [23] and
references therein), and are:

  * –

$I_{1}=\lambda_{1}+\lambda_{2}+\lambda_{3}$,

  * –

$I_{2}=\lambda_{1}\lambda_{2}+\lambda_{1}\lambda_{3}+\lambda_{2}\lambda_{3}$,

  * –

$I_{3}=\lambda_{1}\lambda_{2}\lambda_{3}$.

Consequently, the invariants can be directly linked to the different regions
of the cosmic web (for simplicity, and without loss of generality, for
$\lambda_{\rm th}=0$):

  * •

Knots: $I_{3}>0\>\&\>I_{2}>0\>\&\>I_{1}>\lambda_{1}$

  * •

Filaments:
$\left(I_{3}<0\>\&\>I_{2}<0\right)\;||\;\left(I_{3}<0\>\&\>I_{2}>0\>\&\>\lambda_{3}<I_{1}<\lambda_{3}-\lambda_{2}\lambda_{3}/\lambda_{1}\right)$

  * •

Sheets:
$\left(I_{3}>0\>\&\>I_{2}<0\right)\;||\>\left(I_{3}<0\>\&\>I_{2}>0\>\&\>\lambda_{1}-\lambda_{2}\lambda_{3}/\lambda_{1}<I_{1}<\lambda_{1}\right)$

  * •

Voids: $I_{3}<0\>\&\>I_{2}>0\>\&\>I_{1}<\lambda_{1}$

Figure 2: Cross power spectra of the cosmic web density fields
$\delta_{\phi-\text{web }}$ (solid lines) and $\delta_{\delta-\text{web}}$
(dashed lines) with the halo reference catalogue for different threshold
values $\lambda_{\rm th}$.

It emerges that both the long-range non-local bias (2.9) and the cosmic web
classification are determined by the eigenvalues of the tidal field. In
particular, the invariants can be related to the different bias components in
halo perturbation theory (see [23])

  * –

$\delta=I_{1}$,

  * –

$s^{2}=\frac{2}{3}I_{1}^{2}-2I_{2}$,

  * –

$s^{3}=-I_{1}I_{2}+3I_{3}+\frac{2}{9}I_{1}^{3}$.

The cosmic web classification based on the gravitational potential (the
$\phi$-web) provides information on $I_{1}\sim\delta$, $I_{2}$ and $I_{3}$,
and thus, on $\delta$, $s^{2}$, and $s^{3}$. Thus, we propose applying the
deterministic (3.3) and stochastic (3.4) bias model to the different
$\phi$-web regions to indirectly model long-range non-local bias.
Additionally, since the short-range non-local bias (2.9) is determined by the
density field tensor, $\Gamma_{ij}=\partial_{i}\partial_{j}\delta$, we suggest
sub-classifying the $\phi$-web regions based on the eigenvalues of the density
field tensor. Figure 1 illustrates the density field (left) alongside the
$\phi$-web classification (middle) and the $\delta$-web classification
(right).

We have employed an information theory-based approach to represent the density
field using only four distinct values, corresponding to different cosmic web
regions. Following the eigenvalues, we assign the highest value to knots and
then to filaments sheets and voids in decreasing order. The specific values
themselves are not critical; what is important is ensuring that the difference
between subsequent regions remains consistent to prevent bias toward any
particular region, while preserving the hierarchical order. With this
definition, we can define a density contrast based only on those four values
for both long-range $\delta_{\phi-\rm web}$, and short-range
$\delta_{\delta-\rm web}$ components.

A visual examination of figure 1 reveals the correlation between both
classifications and the density field, with the $\delta$-web providing
information on small scales. Moreover, figure 3 depicts the hierarchical
cosmic web classification scheme for the halo reference catalogue used for
calibration in section 4. In the subsequent section, we demonstrate how
progressively incorporating $\phi$-web and $\delta$-web information
significantly enhances the 1-, 2-, and 3-point statistics, culminating in the
reproduction of the reference catalogue with remarkable accuracy.

We propose selecting the threshold value $\lambda_{\rm th}$ based on the
cross-correlation between the cosmic web density fields $\delta_{\phi-\rm
web}$, $\delta_{\delta-\rm web}$, and the halo reference catalogue (section
4). After exploring a wide range of threshold values, we have determined that
$\lambda_{\rm th}=0.05$ is the most suitable choice. In figure 2, we display
the correlation for $\lambda_{\rm th}=0.05$, along with lower ($\lambda_{\rm
th}=0.0$) and upper ($\lambda_{\rm th}=0.1$) limits.

The prominent feature of this plot is that, on small scales, the $\delta$-web
exhibits a higher correlation with the reference halos compared to the
$\phi$-web. This observation supports the notion that $\delta$-web
classification encompasses information on small scales, effectively modeling
short-range non-local behavior. Across different thresholds, the cross-
correlation with the $\phi$-web remains relatively consistent at both large
and small scales, with a slight increase in correlation at intermediate scales
for higher $\lambda_{\rm th}$ values (up to 0.1, beyond which it starts to
decrease). However, the correlation with the $\delta$-web displays significant
differences for various thresholds: at large and intermediate scales, the
correlation for $\lambda_{\rm th}=0.0$ is notably lower compared to other
cases. Across almost all scales, the $\delta$-web with $\lambda_{\rm th}=0.05$
demonstrates the highest correlation with the reference catalogue. In light of
these findings, we have opted to utilise $\lambda_{\rm th}=0.05$ in this
study.

As a consequence of our hierarchical cosmic web classification, various
quantities related to the local deterministic and stochastic bias components
have to be computed for each sub-region. In particular, the bias parameters
depend on the cosmic web region denoted by the subscript cw: $B(\delta_{\rm
tr}|\delta)|_{\rm cw}:\,\\{\alpha|_{\rm cw},{\rho_{\epsilon}}|_{\rm
cw},\epsilon|_{\rm cw},{\rho_{\epsilon}^{\prime}}|_{\rm
cw},\epsilon^{\prime}|_{\rm cw},\beta|_{\rm cw}\\}$ and, thus, also the
normalisation: $f_{\rm tr}|_{\rm cw}$ depending on the specific bias, volume:
$\langle 1+B(\delta_{\rm tr}|\delta)|_{\rm cw}\rangle_{V}$, and mean number
density: $n_{\rm tr}|_{\rm cw}$.

##  4 Calibration on N-body halo catalog

In this section, we present the application of the method to generate halo
mock catalogues.

###  4.1 Reference catalog: UNIT simulation

We employ a gravity-only simulation from the UNIT
project222http://www.unitsims.org, specifically a fixed-amplitude run denoted
as UNITSIM1 (U1) [39], which is consistent with the cosmological parameters of
Planck Collaboration 2016 [94], i.e., $\Omega_{m}=0.3089$,
$h\equiv\mathrm{H}_{0}/100=0.6774$, $n_{s}=0.9667$ and $\sigma_{8}=0.8147$
(see table 4 in [94]). The simulation tracks the dynamic evolution of
$4096^{3}$ dark matter (DM) particles within a periodic box of 1 Gpc $h^{-1}$
on a side, starting at scale factor $a\equiv 1/(1+z)=0.01$ ($z=99$) with
initial conditions generated with FastPM [95]. The DM particles undergo
evolution up to redshift $z=0$ with the particle-mesh and tree algorithms
implemented in the GADGET code [96]. This setup ensures a mass resolution of
$1.2\times 10^{9}~{}M_{\odot}h^{-1}$. Halos were identified with the phase-
space halo finder ROCKSTAR333https://bitbucket.org/gfcstanford/rockstar [97].
The CONSISTENT444https://bitbucket.org/pbehroozi/consistent-trees algorithm
was applied on top of the ROCKSTAR catalogues to enhance the consistency of
merger histories through the implementation of dynamic tracking for halo
progenitors (for details, please refer to [98]). As for the reference
catalogue used in this research, we focus in the halo distribution at
$z=1.032$ with spherical overdensity masses of $M_{200c}\gtrsim
10^{11}~{}M_{\odot}\,h^{-1}$. This criterion agrees with the expected mass cut
at which dark matter halos can host Emission Line Galaxies (ELGs) (as detailed
in studies such as [99]). For this reason, we have chosen a halo mass cut of
$10^{11}~{}M_{\odot}\,h^{-1}$ yielding a reference catalogue containing a
total of 8,278,311 distinct halos.

Full Halo Sample $\underbrace{\hskip 85.35826pt}$
---
$\delta_{\rm tr}$
$\underbrace{\hskip 113.81102pt}$
---
$\phi$-web
$B_{\rm NL}^{\rm LR}(\delta)$
$\underbrace{\hskip 113.81102pt}$
---
$\delta$-web
$B_{\rm NL}^{\rm SR}(\delta)$
Knots $\left\\{\begin{array}[]{c}V_{\rm fr}=2.6\\\ N_{\rm
fr}=15.8\end{array}\right.$Filaments $\left\\{\begin{array}[]{c}V_{\rm
fr}=25.1\\\ N_{\rm fr}=53.6\end{array}\right.$Sheets
$\left\\{\begin{array}[]{c}V_{\rm fr}=51.1\\\ N_{\rm
fr}=28.2\end{array}\right.$Voids $\left\\{\begin{array}[]{c}V_{\rm fr}=21.2\\\
N_{\rm fr}=2.4\end{array}\right.$Knots $\left\\{\begin{array}[]{c}V_{\rm
fr}=0.7\\\ N_{\rm fr}=6.3\end{array}\right.$ $n_{\rm
p}=2$$\begin{rcases}\vspace{17cm}\end{rcases}$ $B_{\rm L}(\delta)$ $+$
$B^{\epsilon}(\delta)$ Filaments $\left\\{\begin{array}[]{c}V_{\rm fr}=1.2\\\
N_{\rm fr}=7.5\end{array}\right.$ $n_{\rm p}=2$Sheets
$\left\\{\begin{array}[]{c}V_{\rm fr}=0.6\\\ N_{\rm fr}=1.9\end{array}\right.$
$n_{\rm p}=2$Voids $\left\\{\begin{array}[]{c}V_{\rm fr}=0.1\\\ N_{\rm
fr}=0.1\end{array}\right.$ $n_{\rm p}=2$Knots
$\left\\{\begin{array}[]{c}V_{\rm fr}=0.1\\\ N_{\rm fr}=1.2\end{array}\right.$
$n_{\rm p}=2$Filaments $\left\\{\begin{array}[]{c}V_{\rm fr}=4.1\\\ N_{\rm
fr}=19.6\end{array}\right.$ $n_{\rm p}=2$Sheets
$\left\\{\begin{array}[]{c}V_{\rm fr}=13.1\\\ N_{\rm
fr}=26.1\end{array}\right.$ $n_{\rm p}=4$Voids
$\left\\{\begin{array}[]{c}V_{\rm fr}=7.8\\\ N_{\rm fr}=6.7\end{array}\right.$
$n_{\rm p}=4$Knots $\left\\{\begin{array}[]{c}V_{\rm fr}<0.1\\\ N_{\rm
fr}<0.1\end{array}\right.$ $n_{\rm p}=0$Filaments
$\left\\{\begin{array}[]{c}V_{\rm fr}=0.1\\\ N_{\rm fr}=0.3\end{array}\right.$
$n_{\rm p}=2$Sheets $\left\\{\begin{array}[]{c}V_{\rm fr}=8.3\\\ N_{\rm
fr}=10.5\end{array}\right.$ $n_{\rm p}=4$Voids
$\left\\{\begin{array}[]{c}V_{\rm fr}=42.7\\\ N_{\rm
fr}=17.4\end{array}\right.$ $n_{\rm p}=4$Knots
$\left\\{\begin{array}[]{c}V_{\rm fr}<0.1\\\ N_{\rm fr}<0.1\end{array}\right.$
$n_{\rm p}=0$Filaments $\left\\{\begin{array}[]{c}V_{\rm fr}<0.1\\\ N_{\rm
fr}<0.1\end{array}\right.$ $n_{\rm p}=0$Sheets
$\left\\{\begin{array}[]{c}V_{\rm fr}<0.1\\\ N_{\rm fr}<0.1\end{array}\right.$
$n_{\rm p}=0$Voids $\left\\{\begin{array}[]{c}V_{\rm fr}=21.1\\\ N_{\rm
fr}=2.3\end{array}\right.$ $n_{\rm p}=2$ Figure 3: Hierarchical cosmic web
classification of the UNIT halo reference catalogue. The colour-code has been
chosen to match the one of the cosmic web classification shown in figure 1.
The volume fraction $V_{\rm fr}$ is the percentage number of cells in each
region w.r.t to the total box. The number fraction $N_{\rm fr}$ is the
percentage number of halos in each region w.r.t to the full sample. The number
of bias parameters $n_{\rm p}$ used in each region is also indicated. Since
there are only a few halos in voids, two parameters are enough to describe all
the $\delta$-web region. Regions with extremely low number densities ($n_{\rm
tr}|_{\rm cw}$) can be modelled with a constant bias solely determined by
$n_{\rm tr}|_{\rm cw}$. In the bottom of the plot we include the label of both
classifications and their relation with the non-local bias terms.

###  4.2 Mock catalog

The first step involves running ALPT with the WebON code [100] in a cubical
volume of 1 $h^{-1}$Gpc side to $z=1.032$. To this end, we employ a down-
sampled version of the original initial conditions (corresponding to the
UNITSIM1 simulation relying on $4096^{3}$ cells and particles) to a mesh of
2563 cells and particles. This results in a cell side resolution of
approximately 3.9 $h^{-1}$ Mpc. We consistently use the same cosmology as the
UNITSIM1. We define the dark matter density contrast after applying the CIC
mass assignment scheme combined with tetrahedra-tesselation [59]. We then use
the resulting dark matter field to calibrate the bias parameters.
Subsequently, we execute 100 realisations with random initial conditions to
ensure robust statistical analysis.

Furthermore, from the WebON code, we obtain the $\phi$-web and $\delta$-web
classifications for the ALPT dark matter density field using a threshold value
of $\lambda_{\rm th}=0.05$. This classification is also applied to the halo
reference catalogue to determine the number density in each cosmic web region.
The classification scheme is illustrated in figure 3, where the two numbers
$V_{\rm fr}$ and $N_{\rm fr}$ represent the volume and number fraction,
respectively, characterising each region of the cosmic web.

In the second step, we apply the bias model described in (3.3) to each of the
individual regions. Although initially there are six possible free parameters
$\\{\alpha,\beta,\rho_{\epsilon},\epsilon,\rho^{\prime}_{\epsilon},\epsilon^{\prime}\\}$,
we find that we can successfully fit most of the individual regions using only
two parameters, $\\{\alpha,\beta\\}$. For some regions, we require four
parameters by incorporating the exponential cut-off, but in neither case do we
need all six parameters.

We are conducting separate studies to investigate which parameters are
relevant for each galaxy type.

Figure 4: Visual comparison of the reference and a mock halo catalog. The
slices are $\sim$20 $h^{-1}$ Mpc thick and 1000 $h^{-1}$ Mpc side. The
logarithm of the halo density contrast $\delta_{\rm tr}$ is shown and darker
regions are related to lower densities. Both fields share the same initial
conditions down-sampled to resolutions of $\sim 3.9\,h^{-1}$ Mpc cell side.
Below that resolution, stochastic bias dominates and the halo distributions
differ.

###  4.3 Results

Figure 4 depicts a visual comparison between the reference and the mock halo
catalogue. The disparities observed stem from stochasticity induced by the
various realisations below the scale of the down-sampled field. Nevertheless,
both fields exhibit striking similarities, and the large-scale structure is
clearly discernible and consistent in both panels.

A quantitative analysis is presented in figure 5, where each row illustrates
the 1-, 2- and 3-point statistics or, equivalently, the probability density
function PDF (left), the power-spectrum $P(k)$ (middle), and the reduced
bispectrum $Q(\theta_{12})$ (right) for the various bias models considered in
this study. We progressively include more bias terms from the top to the
bottom rows. The summary statistics corresponding to the reference catalogue
is represented with the black solid lines or black dots and corresponding
error bars.

The PDF [101], the power spectrum [102] and the bispectrum [103, 104] are
powerful tools for extracting cosmological information from the large-scale
structure of the Universe by constraining cosmological parameters. The PDF is
understood as the number of cells containing a given number of objects. The
Fourier-space 2- and 3-point functions, i.e., the power spectrum and
bispectrum, respectively, are defined as

$\displaystyle\langle\delta_{\rm tr}(\boldsymbol{k})\delta_{\rm
tr}(\boldsymbol{k}^{\prime})\rangle$ $\displaystyle\equiv$
$\displaystyle(2\pi)^{3}P(k)\delta^{D}(\boldsymbol{k}-\boldsymbol{k}^{\prime})\,,$
(4.1) $\displaystyle\langle\delta_{\rm tr}(\boldsymbol{k}_{1})\delta_{\rm
tr}(\boldsymbol{k}_{2})\delta_{\rm tr}(\boldsymbol{k}_{3})\rangle$
$\displaystyle\equiv$
$\displaystyle(2\pi)^{3}B(k_{1},k_{2},k_{3})\delta^{D}(\boldsymbol{k}_{1}+\boldsymbol{k}_{2}+\boldsymbol{k}_{3})\,,$
(4.2)

where the ensemble average stands for $k$-shell averages, while the reduced
bispectrum is defined as

$Q(\theta_{12}|k_{1},k_{2})\equiv\frac{B(k_{1},k_{2},k_{3})}{P(k_{1})P(k_{2})+P(k_{2})P(k_{3})+P(k_{1})P(k_{3})}\,.$
(4.3)

We have focused on a particular configuration of the reduced bispectrum,
$k_{2}=2k_{1}=0.2\,h$ Mpc-1, which is typically studied, as it is relevant to
BAO and RSD analyses. We are using larger volume simulations based on the
AbacusSummit [105] simulation suite to thoroughly study a larger variety of
configurations for haloes and each galaxy type tracer relevant to the DESI and
EUCLID collaborations. Such studies will be presented in forthcoming papers.
The bispectrum projection is shown as a function of the angle $\theta_{12}$
between the wavevectors ${\bf\it k}_{1}$ and ${\bf\it k}_{2}$, given by
$k_{3}^{2}=k_{2}^{2}\sin^{2}\theta_{12}+(k_{2}\cos\theta_{12}+k_{1})^{2}$.555We
are grateful to Cheng Zhao for making his bispectrum computation code
available as presented in [18].

The resulting mock catalogues are based on the following bias models (figure
5):

  1. 1.

First row: Local power-law bias

  2. 2.

Second row: Local power-law bias combined with threshold bias

  3. 3.

Third row: Local power-law bias combined with threshold bias and long-range
cosmic web classification based on the gravitational potential ($\phi$-web)

  4. 4.

Fourth row: Local power-law bias combined with threshold bias and both long-
range and short-range cosmic web classification based on the gravitational
potential ($\phi$-web) and the density field ($\delta$-web).

As outlined in previous studies [16, 38], the determination of the bias
parameters is performed jointly by minimising the likelihood of the 1-point
PDF and the power spectrum. This joint optimisation is crucial as the skewness
of the PDF is constrained by the bispectrum, aiding in breaking degeneracies
[14]. In particular, we determine the optimal values for our parameters by
sampling their posterior distribution through the affine-invariant emcee
Markov Chain Monte Carlo (MCMC) sampler [106, 107].

Through our approach, it is evident how we progressively enhance the fit to
the 1-, 2-, and 3-point statistics. With the simplest model, we are
significantly divergent from reproducing the reference catalogue. There are
substantial deviations in the PDF and the bispectrum, and the power spectrum
fails particularly towards high wave numbers $k$ (small scales).

The second row of plots demonstrates how incorporating a threshold bias (3.3)
improves statistics, notably by reducing the significant deviation in the
bispectrum. However, the statistics are still not well reproduced. Therefore,
we decided to incorporate information on the different cosmic web environments
(knots, filaments, voids, and sheets) by individually applying the bias model
to each region rather than directly to the entire catalogue. This leads to
significant improvements (third row): the PDF is much closer to the reference,
the power spectrum is reproduced within a 5 percent deviation until $k=0.3\,h$
Mpc-1, and only the wings of the bispectrum are not entirely compatible with
the UNIT simulation catalogue.

Despite these improvements, we aim for percentage accuracy and ultimately
achieve it by incorporating small-scale information through hierarchical
cosmic web classification. Now, the four regions of the $\phi$-web are
subclassified based on the density field tensor, which, as demonstrated in
section 3.3, indirectly models the short-range non-local bias. The bias model
(3.3) is now applied to each of the 16 regions, and the statistical result of
combining them to create the entire catalogue is shown in the fourth row of
figure 5. The PDF is now reproduced with 5 percent accuracy, he power spectrum
is compatible within 2 percent deviations almost up to $k=0.8\,h$ Mpc-1, and
the bispectrum, including the wings, is reproduced within error bars.

We find that the low-density threshold bias, previously required in studies
utilising the PATCHY code, becomes unnecessary. This is because the short-
range $\delta$-web classification already provides a comprehensive enough
description of halo formation physics at the scales analysed in our work. As a
consequence, the hierarchical cosmic web classification allows us to
dramatically reduce the number of parameters $n_{\rm p}$ used in each region,
so that the total number of parameters employed in our study is of 32 instead
of 96 ($=16\text{ (regions) }\times 6\text{ (parameters)}$), as indicated in
figure 5.

It should be emphasised once more that the strength of this approach lies in
the ability to treat each subset of the halo catalogue, corresponding to a
particular sub-region of the cosmic web, as an independent catalogue. This
parallels the practice in clustering analyses where galaxies are divided into
categories such as blue and red galaxies. The accuracy of our approach is
undoubtedly constrained by the availability of sufficient statistics to
effectively constrain the bias parameters of each sub-region. Nonetheless,
this approach is less susceptible to over-fitting compared to machine learning
techniques, as the model (and its degrees of freedom) is constrained by its
analytical formulation for both the deterministic and stochastic bias
components. One crucial advantage over the BAM approach is that it does not
require a kernel to reproduce the power-spectrum. This capability enables us
to efficiently apply our field-level bias mapping technique in redshift shells
to a dark matter lightcone, such as the ones we have computed with the WebON
code.

We are exploring in separate studies transformations of the density field
which improve the description of the cosmic web towards smaller scales and
allows us to explore such scales in the 2- and 3-point statistics.

---
Figure 5: Statistical comparison of the UNIT simulation reference catalogue
(black) and the different mock catalogues (red). The probability density
function PDF (left column), the power spectrum $P(k)$ (middle) and the reduced
bispectrum $Q(\theta_{12})$ (right) are shown for a model describing a
deterministic and stochastic bias model (first row), to which we add a
threshold bias (second) and then, progressively, the $\phi$-web and
$\delta$-web information (third and fourth rows).The dashed red line is the
mean over 100 realisations and the shaded area represents the standard
deviation. The grey area in the subsets represent the 5 percent deviation (2
percent also included in the bottom row).

##  5 Summary and discussion

In this study, we have introduced an innovative approach for generating mock
catalogues which extends the deterministic and stochastic bias model by
incorporating a hierarchical classification of the cosmic web. This
augmentation allows us to effectively model non-local bias, ensuring the
positive definiteness of the expected number counts of tracers.

By integrating information from the cosmic web classification, specifically
the $\phi$-web and $\delta$-web, we significantly enhance the fidelity of our
bias model. This comprehensive framework enables us to capture the influence
of large-scale structure on tracer distributions and accurately reproduce
statistical properties such as the probability density function, power
spectrum and bispectrum.

One of the most remarkable findings of this study is the simplification of the
bias model as we incorporate the hierarchical cosmic web classification. We
observe that each region can typically be modelled with just two parameters,
and in some instances, only four parameters are required. This reduction in
the number of parameters not only streamlines the modelling process but also
enhances the interpretability and efficiency of our approach. It underscores
the effectiveness of leveraging cosmic web information to capture the
complexity of large-scale structure with minimal complexity in the bias model.

Our method offers several key advantages, including its computational
efficiency, high accuracy in reproducing reference statistics up to third
order, and versatility in its application to different tracers of large-scale
structure. In comparison to automatic non-parametric learning methods, e.g.,
[21], our approach offers the advantage of being less susceptible to
overfitting. This is because our bias model, including the stochastic
component, is analytical and determined by only a few parameters.

However, further investigation is needed in future work to verify the
robustness of the model assumptions, such as the specific deviation from
Poissonity set by the negative binomial distribution. This ongoing research
will help to ensure the reliability and validity of our approach across a
range of cosmological scenarios and observational data sets.

Also, we conducted a series of 200 large-volume light-cone dark matter
simulations using the WebON code within cubical volumes of ($10\,h^{-1}$ Gpc)3
[100]. We are currently preparing a series of papers presenting mock galaxy
catalogues generated using this method with various tracers of the large-scale
structure. We also plan to explore the 4-point statistics of our mock
catalogues. Previous studies have shown that accurately reproducing the 1-,
2-, and 3-point statistics lays the foundation for reproducing covariance
matrices [108]. Therefore, we believe that we are on the right track with our
current approach to become accurate in the higher-order statistics. The
hierarchical cosmic web classification presented in this study could
potentially be applied for environmental studies of galaxies to study assembly
bias. The assignment methods of halo or galaxy properties can benefit from the
findings of this paper extending previous approaches with the short-range
classification scheme (see [90]). The combination of this bias model with
effective field theory can be as straightforward as considering different
local bias parameters at various hierarchical cosmic web regions.

Looking ahead, we anticipate further enhancements by continuing to refine the
cosmic web classification with additional terms, such as those based on
velocity bias components obtained with fast gravity solvers (e.g., [32]).
Additionally, we plan to implement a nonlinear mapping of the density field to
imprint missing power from the approximate low-resolution gravity solver and
improve the description towards small scales.

This ongoing development promises to advance our understanding of the complex
interplay between cosmic web environments and tracer distributions in the
Universe. By incorporating these advancements, we aim to further refine our
mock catalogue generation process and utilise this framework for field-level
inference. This will potentially lead to even more accurate representations of
large-scale structure and facilitate deeper insights into cosmological
phenomena.

## Acknowledgements

JMCN, FSK, JGF, FS, and GF acknowledge the Spanish Ministry of Economy and
Competitiveness (MINECO) for financing the Big Data of the Cosmic Web project:
PID2020-120612GB-I00 under which this work has been conceived and carried out,
and the IAC for continuous support to the Cosmology with LSS probes project.
JMCN thanks the Research Summer Grant at the IAC for continuing the study of
the master thesis and the ICE for support during the writing phase of this
paper. GF thanks support from the IJC2020-044343-I grant. DFS thanks support
from the Swiss National Science Foundation (SNF) 200020_175751 and
200020_207379.

## References

  * [1] M. Levi, C. Bebek, T. Beers, R. Blum, R. Cahn, D. Eisenstein et al., _The desi experiment, a whitepaper for snowmass 2013_ , 2013.
  * [2] R. Laureijs, J. Amiaux, S. Arduini, J.L. Auguères, J. Brinchmann, R. Cole et al., _Euclid definition study report_ , 2011.
  * [3] F. Bernardeau, S. Colombi, E. Gaztañaga and R. Scoccimarro, _Large-scale structure of the universe and cosmological perturbation theory_ , _Physics Reports_ 367 (2002) 1.
  * [4] J.R. Bond and S.T. Myers, _The Peak-Patch Picture of Cosmic Catalogs. I. Algorithms_ , _ApJS_ 103 (1996) 1.
  * [5] J.R. Bond and S.T. Myers, _The Peak-Patch Picture of Cosmic Catalogs. II. Validation_ , _ApJS_ 103 (1996) 41.
  * [6] J.R. Bond and S.T. Myers, _The Peak-Patch Picture of Cosmic Catalogs. III. Application to Clusters_ , _ApJS_ 103 (1996) 63.
  * [7] G. Stein, M.A. Alvarez and J.R. Bond, _The mass-Peak Patch algorithm for fast generation of deep all-sky dark matter halo catalogues and its N-body validation_ , _MNRAS_ 483 (2019) 2236 [1810.07727].
  * [8] P. Monaco, T. Theuns, G. Taffoni, F. Governato, T. Quinn and J. Stadel, _Predicting the number, spatial distribution, and merging history of dark matter halos_ , _The Astrophysical Journal_ 564 (2002) 8.
  * [9] P. Monaco, E. Sefusatti, S. Borgani, M. Crocce, P. Fosalba, R.K. Sheth et al., _An accurate tool for the fast generation of dark matter halo catalogues_ , _MNRAS_ 433 (2013) 2389 [1305.1505].
  * [10] R. Scoccimarro and R.K. Sheth, _PTHALOS: a fast method for generating mock galaxy distributions_ , _MNRAS_ 329 (2002) 629 [astro-ph/0106120].
  * [11] M. Manera, R. Scoccimarro, W.J. Percival, L. Samushia, C.K. McBride, A.J. Ross et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: a large sample of mock galaxy catalogues_ , _MNRAS_ 428 (2013) 1036 [1203.6609].
  * [12] M. Manera, L. Samushia, R. Tojeiro, C. Howlett, A.J. Ross, W.J. Percival et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: mock galaxy catalogues for the low-redshift sample_ , _MNRAS_ 447 (2015) 437 [1401.4171].
  * [13] F.-S. Kitaura, G. Yepes and F. Prada, _Modelling baryon acoustic oscillations with perturbation theory and stochastic halo biasing_ , _Monthly Notices of the Royal Astronomical Society: Letters_ 439 (2013) L21 [https://academic.oup.com/mnrasl/article-pdf/439/1/L21/3421996/slt172.pdf].
  * [14] F.-S. Kitaura, H. Gil-Marín, C.G. Scóccola, C.-H. Chuang, V. Müller, G. Yepes et al., _Constraining the halo bispectrum in real and redshift space from perturbation theory and non-linear stochastic bias_ , _MNRAS_ 450 (2015) 1836 [1407.1236].
  * [15] F.-S. Kitaura, S. Rodríguez-Torres, C.-H. Chuang, C. Zhao, F. Prada, H. Gil-Marín et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: mock galaxy catalogues for the BOSS Final Data Release_ , _MNRAS_ 456 (2016) 4156 [1509.06400].
  * [16] M. Vakili, F.-S. Kitaura, Y. Feng, G. Yepes, C. Zhao, C.-H. Chuang et al., _Accurate halo-galaxy mocks from automatic bias estimation and particle mesh gravity solvers_ , _MNRAS_ 472 (2017) 4144 [1701.03765].
  * [17] C.-H. Chuang, F.-S. Kitaura, F. Prada, C. Zhao and G. Yepes, _Ezmocks: extending the zel’dovich approximation to generate mock galaxy catalogues with accurate clustering statistics_ , _Monthly Notices of the Royal Astronomical Society_ 446 (2014) 2621–2628.
  * [18] C. Zhao, C.-H. Chuang, J. Bautista, A. de Mattia, A. Raichoor, A.J. Ross et al., _The completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: 1000 multi-tracer mock catalogues with redshift evolution and systematics for galaxies and quasars of the final data release_ , _MNRAS_ 503 (2021) 1149 [2007.08997].
  * [19] S. Avila, S.G. Murray, A. Knebe, C. Power, A.S.G. Robotham and J. Garcia-Bellido, _HALOGEN: A tool for fast generation of mock halo catalogues_ , _Mon. Not. Roy. Astron. Soc._ 450 (2015) 1856 [1412.5228].
  * [20] A. Izard, M. Crocce and P. Fosalba, _ICE-COLA: towards fast and accurate synthetic galaxy catalogues optimizing a quasi-N-body method_ , _MNRAS_ 459 (2016) 2327 [1509.04685].
  * [21] A. Balaguera-Antolínez, F.-S. Kitaura, M. Pellejero-Ibáñez, C. Zhao and T. Abel, _BAM: bias assignment method to generate mock catalogues_ , _Monthly Notices of the Royal Astronomical Society: Letters_ 483 (2018) L58 [https://academic.oup.com/mnrasl/article-pdf/483/1/L58/54700852/mnrasl_483_1_l58.pdf].
  * [22] A. Balaguera-Antolínez, F.-S. Kitaura, M. Pellejero-Ibáñez, M. Lippich, C. Zhao, A.G. Sánchez et al., _One simulation to have them all: performance of the Bias Assignment Method against N-body simulations_ , _MNRAS_ 491 (2020) 2565 [1906.06109].
  * [23] F.S. Kitaura, A. Balaguera-Antolínez, F. Sinigaglia and M. Pellejero-Ibáñez, _The cosmic web connection to the dark matter halo distribution through gravity_ , _MNRAS_ 512 (2022) 2245 [2005.11598].
  * [24] A. Balaguera-Antolínez, F.-S. Kitaura, S. Alam, C.-H. Chuang, Y. Yu, G. Favole et al., _DESI mock challenge. Halo and galaxy catalogues with the bias assignment method_ , _A &A_ 673 (2023) A130 [2211.10640].
  * [25] C.-H. Chuang, C. Zhao, F. Prada, E. Munari, S. Avila, A. Izard et al., _nIFTy cosmology: Galaxy/halo mock catalogue comparison project on clustering statistics_ , _MNRAS_ 452 (2015) 686 [1412.7729].
  * [26] M. Lippich, A.G. Sánchez, M. Colavincenzo, E. Sefusatti, P. Monaco, L. Blot et al., _Comparing approximate methods for mock catalogues and covariance matrices - I. Correlation function_ , _MNRAS_ 482 (2019) 1786 [1806.09477].
  * [27] L. Blot, M. Crocce, E. Sefusatti, M. Lippich, A.G. Sánchez, M. Colavincenzo et al., _Comparing approximate methods for mock catalogues and covariance matrices II: power spectrum multipoles_ , _MNRAS_ 485 (2019) 2806 [1806.09497].
  * [28] M. Colavincenzo, E. Sefusatti, P. Monaco, L. Blot, M. Crocce, M. Lippich et al., _Comparing approximate methods for mock catalogues and covariance matrices - III: bispectrum_ , _MNRAS_ 482 (2019) 4883 [1806.09499].
  * [29] F.-S. Kitaura, M. Ata, S.A. Rodríguez-Torres, M. Hernández-Sánchez, A. Balaguera-Antolínez and G. Yepes, _COSMIC BIRTH: efficient Bayesian inference of the evolving cosmic web from galaxy surveys_ , _MNRAS_ 502 (2021) 3456 [1911.00284].
  * [30] M. Ata, F.-S. Kitaura and V. Müller, _Bayesian inference of cosmic density fields from non-linear, scale-dependent, and stochastic biased tracers_ , _MNRAS_ 446 (2015) 4250 [1408.2566].
  * [31] F.-S. Kitaura and S. Heß, _Cosmological structure formation with augmented lagrangian perturbation theory_ , _Monthly Notices of the Royal Astronomical Society: Letters_ 435 (2013) L78.
  * [32] F.S. Kitaura, F. Sinigaglia, A. Balaguera-Antolínez and G. Favole, _The cosmic web from perturbation theory_ , _A &A_ 683 (2024) A215 [2301.03648].
  * [33] Y. Feng, M.-Y. Chu, U. Seljak and P. McDonald, _FASTPM: a new scheme for fast simulations of dark matter and haloes_ , _Monthly Notices of the Royal Astronomical Society_ 463 (2016) 2273 [1603.00476].
  * [34] S. Tassev, M. Zaldarriaga and D.J. Eisenstein, _Solving large scale structure in ten easy steps with COLA_ , _Journal of Cosmoly and Astroparticle Physics_ 6 (2013) 036 [1301.0322].
  * [35] M. Pellejero-Ibañez, A. Balaguera-Antolínez, F.-S. Kitaura, R.E. Angulo, G. Yepes, C.-H. Chuang et al., _The bias of dark matter tracers: assessing the accuracy of mapping techniques_ , _Monthly Notices of the Royal Astronomical Society_ 493 (2020) 586 [https://academic.oup.com/mnras/article-pdf/493/1/586/32513186/staa270.pdf].
  * [36] F. Sinigaglia, F.-S. Kitaura, A. Balaguera-Antolínez, K. Nagamine, M. Ata, I. Shimizu et al., _The bias from hydrodynamic simulations: Mapping baryon physics onto dark matter fields_ , _The Astrophysical Journal_ 921 (2021) 66.
  * [37] F. Sinigaglia, F.-S. Kitaura, A. Balaguera-Antolínez, I. Shimizu, K. Nagamine, M. Sánchez-Benavente et al., _Mapping the Three-dimensional Ly $\alpha$ Forest Large-scale Structure in Real and Redshift Space_, _ApJ_ 927 (2022) 230 [2107.07917].
  * [38] F. Sinigaglia, F.S. Kitaura, K. Nagamine, Y. Oku and A. Balaguera-Antolínez, _Field-level Lyman- $\alpha$ forest modeling in redshift space via augmented nonlocal Fluctuating Gunn-Peterson Approximation_, _A &A_ 682 (2024) A21 [2305.10428].
  * [39] C.-H. Chuang, G. Yepes, F.-S. Kitaura, M. Pellejero-Ibanez, S. Rodríguez-Torres, Y. Feng et al., _UNIT project: Universe N-body simulations for the Investigation of Theoretical models from galaxy surveys_ , _MNRAS_ 487 (2019) 48 [1811.02111].
  * [40] P. McDonald and A. Roy, _Clustering of dark matter tracers: generalizing bias for the coming era of precision LSS_ , _Journal of Cosmology and Astroparticle Physics_ 2009 (2009) 020.
  * [41] M. Crocce and R. Scoccimarro, _Renormalized cosmological perturbation theory_ , _Phys. Rev. D_ 73 (2006) 063519 [astro-ph/0509418].
  * [42] D. Baumann, A. Nicolis, L. Senatore and M. Zaldarriaga, _Cosmological non-linearities as an effective fluid_ , _J. Cosmology Astropart. Phys_ 2012 (2012) 051 [1004.2488].
  * [43] J.J.M. Carrasco, M.P. Hertzberg and L. Senatore, _The effective field theory of cosmological large scale structures_ , _Journal of High Energy Physics_ 2012 (2012) 82 [1206.2926].
  * [44] E. Pajer and M. Zaldarriaga, _On the renormalization of the effective field theory of large scale structures_ , _J. Cosmology Astropart. Phys_ 2013 (2013) 037 [1301.7182].
  * [45] R.A. Porto, L. Senatore and M. Zaldarriaga, _The Lagrangian-space Effective Field Theory of large scale structures_ , _J. Cosmology Astropart. Phys_ 2014 (2014) 022 [1311.2168].
  * [46] R. Angulo, M. Fasiello, L. Senatore and Z. Vlah, _On the statistics of biased tracers in the Effective Field Theory of Large Scale Structures_ , _J. Cosmology Astropart. Phys_ 2015 (2015) 029 [1503.08826].
  * [47] L. Senatore and M. Zaldarriaga, _The IR-resummed Effective Field Theory of Large Scale Structures_ , _J. Cosmology Astropart. Phys_ 2015 (2015) 013 [1404.5954].
  * [48] L. Senatore, _Bias in the effective field theory of large scale structures_ , _J. Cosmology Astropart. Phys_ 2015 (2015) 007 [1406.7843].
  * [49] Z. Vlah, M. White and A. Aviles, _A Lagrangian effective field theory_ , _J. Cosmology Astropart. Phys_ 2015 (2015) 014 [1506.05264].
  * [50] M. Schmittfull, M. Simonović, V. Assassi and M. Zaldarriaga, _Modeling biased tracers at the field level_ , _Phys. Rev. D_ 100 (2019) 043514 [1811.10640].
  * [51] D. Forero-Sánchez, F.-S. Kitaura, F. Sinigaglia, J. María Coloma-Nodal and J.-P. Kneib, _CosmoMIA: Cosmic Web-based redshift space halo distribution_ , _arXiv e-prints_ (2024) arXiv:2402.17581 [2402.17581].
  * [52] T. Buchert, _Lagrangian Theory of Gravitational Instability of Friedman-Lemaitre Cosmologies - a Generic Third-Order Model for Nonlinear Clustering_ , _MNRAS_ 267 (1994) 811 [astro-ph/9309055].
  * [53] F.R. Bouchet, S. Colombi, E. Hivon and R. Juszkiewicz, _Perturbative Lagrangian approach to gravitational instability._ , _A &A_ 296 (1995) 575 [astro-ph/9406013].
  * [54] P. Catelan, _Lagrangian dynamics in non-flat universes and non-linear gravitational evolution_ , _MNRAS_ 276 (1995) 115 [astro-ph/9406016].
  * [55] M. Crocce, S. Pueblas and R. Scoccimarro, _Transients from initial conditions in cosmological simulations_ , _MNRAS_ 373 (2006) 369 [astro-ph/0606505].
  * [56] F. Bernardeau, _The Nonlinear Evolution of Rare Events_ , _ApJ_ 427 (1994) 51 [astro-ph/9311066].
  * [57] R. Mohayaee, H. Mathis, S. Colombi and J. Silk, _Reconstruction of primordial density fields_ , _Monthly Notices of the Royal Astronomical Society_ 365 (2006) 939 [https://academic.oup.com/mnras/article-pdf/365/3/939/2910100/365-3-939.pdf].
  * [58] M.C. Neyrinck, _Quantifying distortions of the lagrangian dark-matter mesh in cosmology_ , _Monthly Notices of the Royal Astronomical Society_ 428 (2013) 141.
  * [59] T. Abel, O. Hahn and R. Kaehler, _Tracing the dark matter sheet in phase space_ , _Monthly Notices of the Royal Astronomical Society_ 427 (2012) 61 [1111.3944].
  * [60] M.C. Neyrinck, _Truthing the stretch: non-perturbative cosmological realizations with multiscale spherical collapse_ , _MNRAS_ 455 (2016) L11 [1503.07534].
  * [61] J.N. Grieb, A.G. Sánchez, S. Salazar-Albornoz and C. Dalla Vecchia, _Gaussian covariance matrices for anisotropic galaxy clustering measurements_ , _MNRAS_ 457 (2016) 1577 [1509.04293].
  * [62] V. Desjacques, D. Jeong and F. Schmidt, _Large-scale galaxy bias_ , _Phys. Rep._ 733 (2018) 1 [1611.09787].
  * [63] L. Gao, V. Springel and S.D.M. White, _The age dependence of halo clustering_ , _MNRAS_ 363 (2005) L66 [astro-ph/0506510].
  * [64] D.J. Croton, L. Gao and S.D.M. White, _Halo assembly bias and its effects on galaxy clustering_ , _MNRAS_ 374 (2007) 1303 [astro-ph/0605636].
  * [65] N. Dalal, M. White, J.R. Bond and A. Shirokov, _Halo Assembly Bias in Hierarchical Structure Formation_ , _ApJ_ 687 (2008) 12 [0803.3453].
  * [66] R. Tojeiro, E. Eardley, J.A. Peacock, P. Norberg, M. Alpaslan, S.P. Driver et al., _Galaxy and Mass Assembly (GAMA): halo formation times and halo assembly bias on the cosmic web_ , _MNRAS_ 470 (2017) 3720 [1612.08595].
  * [67] X. Yang, Y. Zhang, T. Lu, H. Wang, F. Shi, D. Tweed et al., _Revealing the cosmic web-dependent halo bias_ , _The Astrophysical Journal_ 848 (2017) 60.
  * [68] J.N. Fry and E. Gaztanaga, _Biasing and Hierarchical Statistics in Large-Scale Structure_ , _ApJ_ 413 (1993) 447 [astro-ph/9302009].
  * [69] R. Cen and J.P. Ostriker, _Cold Dark Matter Cosmogony with Hydrodynamics and Galaxy Formation: Galaxy Properties at Redshift Zero_ , _ApJ_ 417 (1993) 415.
  * [70] S. de la Torre and J.A. Peacock, _Reconstructing the distribution of haloes and mock galaxies below the resolution limit in cosmological simulations_ , _MNRAS_ 435 (2013) 743 [1212.3615].
  * [71] N. Kaiser, _On the spatial correlations of Abell clusters._ , _ApJ_ 284 (1984) L9.
  * [72] M.C. Neyrinck, M.A. Aragón-Calvo, D. Jeong and X. Wang, _A halo bias function measured deeply into voids without stochasticity_ , _MNRAS_ 441 (2014) 646 [1309.6641].
  * [73] R. García and E. Rozo, _Halo exclusion criteria impacts halo statistics_ , _MNRAS_ 489 (2019) 4170 [1903.01709].
  * [74] A. Dekel and O. Lahav, _Stochastic Nonlinear Galaxy Biasing_ , _ApJ_ 520 (1999) 24 [astro-ph/9806193].
  * [75] R.K. Sheth and G. Lemson, _Biasing and the distribution of dark matter haloes_ , _MNRAS_ 304 (1999) 767 [astro-ph/9808138].
  * [76] P.J.E. Peebles, _The large-scale structure of the universe_ (1980).
  * [77] R.S. Somerville, G. Lemson, Y. Sigad, A. Dekel, G. Kauffmann and S.D.M. White, _Non-linear stochastic galaxy biasing in cosmological simulations_ , _MNRAS_ 320 (2001) 289 [astro-ph/9912073].
  * [78] R. Casas-Miranda, H.J. Mo, R.K. Sheth and G. Boerner, _On the distribution of haloes, galaxies and mass_ , _Monthly Notices of the Royal Astronomical Society_ 333 (2002) 730–738.
  * [79] W.C. Saslaw and A.J.S. Hamilton, _Thermodynamics and galaxy clustering - Nonlinear theory of high order correlations_ , _ApJ_ 276 (1984) 13.
  * [80] R.K. Sheth, _Press–Schechter, thermodynamics and gravitational clustering_ , _Monthly Notices of the Royal Astronomical Society_ 274 (1995) 213 [https://academic.oup.com/mnras/article-pdf/274/1/213/18539982/mnras274-0213.pdf].
  * [81] V. Wild, J.A. Peacock, O. Lahav, E. Conway, S. Maddox, I.K. Baldry et al., _The 2dF Galaxy Redshift Survey: stochastic relative biasing between galaxy populations_ , _Monthly Notices of the Royal Astronomical Society_ 356 (2005) 247 [https://academic.oup.com/mnras/article-pdf/356/1/247/3589649/356-1-247.pdf].
  * [82] S. de la Torre and J.A. Peacock, _Reconstructing the distribution of haloes and mock galaxies below the resolution limit in cosmological simulations_ , _Monthly Notices of the Royal Astronomical Society_ 435 (2013) 743.
  * [83] J.R. Bond, L. Kofman and D. Pogosyan, _How filaments of galaxies are woven into the cosmic web_ , _Nature_ 380 (1996) 603 [astro-ph/9512141].
  * [84] R. van de Weygaert and W. Schaap, _The cosmic web: Geometric analysis_ , in _Data Analysis in Cosmology_ , pp. 291–413, Springer Berlin Heidelberg (2008), DOI.
  * [85] I. Gott, J. Richard, M. Jurić, D. Schlegel, F. Hoyle, M. Vogeley, M. Tegmark et al., _A Map of the Universe_ , _ApJ_ 624 (2005) 463 [astro-ph/0310571].
  * [86] Y.B. Zel’dovich, _Gravitational instability: An approximate theory for large density perturbations._ , _A &A_ 5 (1970) 84.
  * [87] O. Hahn, C. Porciani, C.M. Carollo and A. Dekel, _Properties of dark matter haloes in clusters, filaments, sheets and voids_ , _Monthly Notices of the Royal Astronomical Society_ 375 (2007) 489.
  * [88] J.E. Forero–Romero, Y. Hoffman, S. Gottlöber, A. Klypin and G. Yepes, _A dynamical classification of the cosmic web_ , _Monthly Notices of the Royal Astronomical Society_ 396 (2009) 1815 [https://academic.oup.com/mnras/article-pdf/396/3/1815/5804803/mnras0396-1815.pdf].
  * [89] D. Martizzi, M. Vogelsberger, M.C. Artale, M. Haider, P. Torrey, F. Marinacci et al., _Baryons in the Cosmic Web of IllustrisTNG – I: gas in knots, filaments, sheets, and voids_ , _Monthly Notices of the Royal Astronomical Society_ 486 (2019) 3766 [https://academic.oup.com/mnras/article-pdf/486/3/3766/28548252/stz1106.pdf].
  * [90] C. Zhao, F.-S. Kitaura, C.-H. Chuang, F. Prada, G. Yepes and C. Tao, _Halo mass distribution reconstruction across the cosmic web_ , _Monthly Notices of the Royal Astronomical Society_ 451 (2015) 4266 [https://academic.oup.com/mnras/article-pdf/451/4/4266/3901283/stv1262.pdf].
  * [91] K.-G. Lee and M. White, _Revealing the z ~2.5 Cosmic Web with 3D Ly$\alpha$ Forest Tomography: a Deformation Tensor Approach_, _ApJ_ 831 (2016) 181 [1603.04441].
  * [92] A. Krolewski, K.-G. Lee, Z. Lukić and M. White, _Measuring Alignments between Galaxies and the Cosmic Web at z $\sim$ 2-3 Using IGM Tomography_, _ApJ_ 837 (2017) 31 [1612.00067].
  * [93] B. Horowitz, K.-G. Lee, M. White, A. Krolewski and M. Ata, _TARDIS. I. A Constrained Reconstruction Approach to Modeling the z $\sim$ 2.5 Cosmic Web Probed by Ly$\alpha$ Forest Tomography_, _ApJ_ 887 (2019) 61 [1903.09049].
  * [94] Planck Collaboration, P.A.R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont et al., _Planck 2015 results. XIII. Cosmological parameters_ , _A &A_ 594 (2016) A13 [1502.01589].
  * [95] Y. Feng, M.-Y. Chu, U. Seljak and P. McDonald, _FASTPM: a new scheme for fast simulations of dark matter and haloes_ , _MNRAS_ 463 (2016) 2273 [1603.00476].
  * [96] V. Springel, _The cosmological simulation code GADGET-2_ , _MNRAS_ 364 (2005) 1105 [astro-ph/0505010].
  * [97] P.S. Behroozi, R.H. Wechsler and H.-Y. Wu, _The ROCKSTAR Phase-space Temporal Halo Finder and the Velocity Offsets of Cluster Cores_ , _ApJ_ 762 (2013) 109 [1110.4372].
  * [98] P.S. Behroozi, R.H. Wechsler, H.-Y. Wu, M.T. Busha, A.A. Klypin and J.R. Primack, _Gravitationally Consistent Halo Catalogs and Merger Trees for Precision Cosmology_ , _ApJ_ 763 (2013) 18 [1110.4370].
  * [99] S. Alam, J.A. Peacock, K. Kraljic, A.J. Ross and J. Comparat, _Multitracer extension of the halo model: probing quenching and conformity in eBOSS_ , _MNRAS_ 497 (2020) 581 [1910.05095].
  * [100] F.-S. Kitaura, F. Sinigaglia, et al._, in preparation_ (2024) .
  * [101] C. Uhlemann, O. Friedrich, F. Villaescusa-Navarro, A. Banerjee and S. Codis, _Fisher for complements: extracting cosmology and neutrino mass from the counts-in-cells PDF_ , _MNRAS_ 495 (2020) 4006 [1911.11158].
  * [102] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J.A. Blazek et al., _The clustering of galaxies in the completed sdss-iii baryon oscillation spectroscopic survey: cosmological analysis of the dr12 galaxy sample_ , _Monthly Notices of the Royal Astronomical Society_ 470 (2017) 2617–2652.
  * [103] E. Sefusatti, M. Crocce, S. Pueblas and R. Scoccimarro, _Cosmology and the bispectrum_ , _Phys. Rev. D_ 74 (2006) 023522.
  * [104] N.S. Sugiyama, D. Yamauchi, T. Kobayashi, T. Fujita, S. Arai, S. Hirano et al., _New constraints on cosmological modified gravity theories from anisotropic three-point correlation functions of BOSS DR12 galaxies_ , _Monthly Notices of the Royal Astronomical Society_ 523 (2023) 3133 [https://academic.oup.com/mnras/article-pdf/523/2/3133/50540934/stad1505.pdf].
  * [105] N.A. Maksimova, L.H. Garrison, D.J. Eisenstein, B. Hadzhiyska, S. Bose and T.P. Satterthwaite, _ABACUSSUMMIT: a massive set of high-accuracy, high-resolution N-body simulations_ , _MNRAS_ 508 (2021) 4017 [2110.11398].
  * [106] J. Goodman and J. Weare, _Ensemble samplers with affine invariance_ , _Communications in Applied Mathematics and Computational Science_ 5 (2010) 65.
  * [107] D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, _emcee: The MCMC Hammer_ , _PASP_ 125 (2013) 306 [1202.3665].
  * [108] F. Baumgarten and C.-H. Chuang, _Robustness of the covariance matrix for galaxy clustering measurements_ , _MNRAS_ 480 (2018) 2535 [1802.04462].


* Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Glendenning and Kettner (2000) N. K. Glendenning and C. Kettner, Astron. Astrophys. 353, L9 (2000), arXiv:astro-ph/9807155 .
  * Alford and Han (2016) M. G. Alford and S. Han, Eur. Phys. J. A 52, 62 (2016).
  * Christian and Schaffner-Bielich (2022) J.-E. Christian and J. Schaffner-Bielich, Astrophys. J. 935, 122 (2022), arXiv:2109.04191 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Koehn  _et al._ (2024) H. Koehn _et al._ , (2024), arXiv:2402.04172 [astro-ph.HE] .
  * Kocsis (2022) B. Kocsis, in  _Handbook of Gravitational Wave Astronomy_ (2022) p. 15.
  * Tsang (2013) D. Tsang, Astrophys. J. 777, 103 (2013), arXiv:1307.3554 [astro-ph.HE] .