11institutetext:  Instituto de Astrofísica de Canarias, s/n, E-38205, La
Laguna, Tenerife, Spain 22institutetext: Departamento de Astrofísica,
Universidad de La Laguna, E-38206, La Laguna, Tenerife, Spain 33institutetext:
Departamento de Física, Universidad Técnica Federico Santa María, Casilla
110-V, Avda. España 1680, Valparaíso, Chile

# Secondary halo bias through cosmic time I: Scaling relations and the
connection with the cosmic web

Andrés Balaguera-Antolínez 1,2<EMAIL_ADDRESS>0000-0001-5028-3035 Antonio D.
Montero-Dorta3 and Ginevra Favole1,2

(Received /Accepted)

###### Abstract

Context. The spatial distribution of dark matter halos carries cosmological
and astrophysical information. Cosmological information can be considered to
be contained in the connection between halo main properties and the large-
scale halo bias, while the astrophysical information would be encoded in the
scaling relations between halo properties. The combination of these two
contributions leads to the effect of secondary halo bias.

Aims. To measure the signal of secondary halo bias as a function of a variety
of intrinsic and environmental halo properties, and characterize its
statistical significance as a function of cosmological redshift.

Methods. Using fixed and paired $N$-body simulations of dark-matter halos -the
UNIT simulation - with masses above $10^{11}M_{\odot}h^{-1}$ identified over a
wide range of cosmological redshifts ($0<z<5$), we explore the behavior of the
scaling relations among different halo properties. We include novel
environmental properties based on the halo distribution as well as the
underlying dark-matter field. We implement an object-by-object estimator of
large-scale effective bias and test its validity against standard approaches.
With a bias assigned to each tracer, we perform a statistical analysis aiming
at characterizing the distribution of bias and the signal of secondary halo
bias.

Results. We show how the halo scaling relations linking direct probes of the
halo potential well do not depend on the environment. On the contrary, links
between halo mass and the so called set of secondary halo properties are
sensitive to the cosmological environment. We show that the signal of
secondary bias derives statistically from secondary correlations beyond the
standard link to halo mass.

Conclusions. We show that the secondary bias arise through non-local and/or
environmental properties related either to the halo distribution or to the
properties of the underlying dark-matter field. In particular, properties such
as the tidal field (measure of the anisotropy of the density field) and the
local Mach number (a measure of the local kinetic temperature of the halo
distribution) generates the signals of secondary bias with the highest
significance. We propose applications of the assignment of individual bias for
the generation of mock catalogs containing the signal of secondary bias, as
well as a series of cosmological analyses aiming at mining large galaxy data
sets.

###### Key Words.:

Cosmology: theory, large-scale structure of the Universe, dark matter –
Galaxies:

##  1 Introduction

It is well known that the clustering of dark matter halos strongly depends on
their internal properties. In this context, halo mass is usually considered to
be responsible for one of the primary dependencies, as a direct manifestation
of the more fundamental dependence of halo bias on the peak height of density
fluctuations, $\nu$ (e.g. Press & Schechter, 1974; Sheth & Tormen, 1999a;
Sheth et al., 2001; Sheth & Tormen, 2002). More massive halos are more
strongly clustered than less massive halos, as indeed described by the
$\Lambda$-cold dark matter ($\Lambda$-CDM) structure formation formalism.
However, a number of additional secondary dependencies at fixed halo mass have
been discovered mostly using cosmological simulations (see e.g. Sheth &
Tormen, 2004; Gao et al., 2005; Wechsler et al., 2006; Gao & White, 2007;
Pandey & Bharadwaj, 2007; Dalal et al., 2008a; Angulo et al., 2008; Li et al.,
2008; Faltenbacher & White, 2010; Lazeyras et al., 2017; Salcedo et al., 2018;
Han et al., 2019; Mao et al., 2018; Sato-Polito et al., 2019; Johnson et al.,
2019; Ramakrishnan et al., 2019; Montero-Dorta et al., 2020; Tucci et al.,
2021; Montero-Dorta et al., 2021). As an example, lower mass halos that
assemble a significant portion of their mass early on are more strongly
clustered than those that form at later times, an effect named halo assembly
bias. This type of secondary dependencies have been measured for a number of
internal properties, including halo concentration, spin, shape, anisotropic
velocity dispersion, and a variety of definitions of halo age and proxies of
assembly history (see e.g. Gao & White, 2007; Dalal et al., 2008a, a;
Faltenbacher & White, 2010; Sato-Polito et al., 2019; Montero-Dorta et al.,
2021, just to name a few).

Throughout this work, we refer to the complete set of secondary dependencies
of halo bias at fixed halo mass as secondary halo bias. These also include
environmental properties and parameters characterizing the tidal field around
halos, which are known to display such signals111Note that it is common to
refer to all secondary dependencies as halo assembly bias, instead of just the
dependence on assembly history. We avoid this nomenclature since it is not
established that all secondary dependencies of halo bias originate from the
same set of physical mechanisms. (see e.g. Paranjape et al., 2018). Despite
the abundance of measurements, the physical origins of secondary halo bias are
still not completely characterized, although several results point precisely
towards tidal fields as the main drivers. Dalal et al. (2008b) introduced the
notion that low-mass halo assembly bias, as opposed to its high-mass
counterpart, could be due to a sub-population of halos whose accretion was
halted early on. Hahn et al. (2009) subsequently claimed that tidal effects
produced by a neighbouring massive halo could be the dominant driver for this
suppressed growth. These results are align well with the works of
Borzyszkowski et al. (2017) and Musso et al. (2018). Borzyszkowski et al.
(2017), based on zoom-in simulations of a small number of halos, concluded
that low-mass assembly bias is due with the existence of same-mass of
“stalled” and “accreting” halos, typically associated with filaments and
nodes, respectively. Musso et al. (2018) extended this analysis to take into
account the whole geometry of tides. It has been also shown that the secondary
bias signals (i.e., not only assembly bias) correlate with the anisotropy of
the tidal tensor (parametrized as $\alpha$) at fixed halo mass (see e.g.
Paranjape et al., 2018; Ramakrishnan et al., 2019). Importantly, the exact
relations between these mechanisms and the different halo properties are still
unclear.

The measurement of halo bias (primary or secondary) is usually performed on
the basis of the two-point correlation function, either in configuration or in
Fourier space (see e.g. Kravtsov & Klypin, 1999; Pollack et al., 2012;
Montero-Dorta et al., 2021). In this paper, however, we implement an object-
by-object estimator of halo bias (Paranjape et al., 2018), which allows for a
robust statistical characterization of halo bias and its dependencies on a
number halo properties. This bias estimation is employed to explore, for the
first time, the signal of secondary bias in the UNIT simulation (Chuang et
al., 2019), which followed the evolution of a cosmological volume of
$1($Gpc$/h)^{3}$ with outputs from redshift $z\sim 5$. We dissect secondary
bias by analyzing a wide range of both intrinsic halo and environmental
properties. The assessment of secondary bias performed in this work includes
novel halo properties such as the _relative Mach number_ and the _neighbour
statistics_ , as well as the evaluation of properties which are inherent to
the dark-matter field. Our analysis is not only relevant in the context of
secondary bias. We also perform several numerical tests to assess the strength
of the dependencies of the halo scaling relations on secondary properties and
their impact on the signal of halo primary and secondary bias.

This work is the first of a series of papers where we intend to shed light on
several key aspects of secondary bias, by taking advantage of the singular
cosmological characteristics of the paired-fixed amplitude UNIT simulation,
and the novel angle that the object-by-object bias estimate brings. The
outline of this paper is the following. In § 2 we introduce the UNIT
simulation an describe the different halo properties employed in the analysis.
In particular, we characterize the correlation among the different halo
properties and the impact of the environment over their scaling relations. In
§4 we show the measurements of secondary bias based on an object-by-object
estimator and show the scaling relations between large-scale bias and the halo
properties. Section 5 presents an analytical characterization of the bias
dependencies. In §6, we discuss forthcoming research related to the analysis
of halo bias at the tracer level. Finally, in §7, we discuss our results and
present the main conclusions of this work.

##  2 The UNIT simulation

The UNIT suite of simulations (UNITsim hereafter) 222http://www.unitsims.org/
(Chuang et al., 2019) is a set of $N$-body simulations of dark matter in
cosmological volumes evolved from $z=99$ until $z=0$. Initial conditions are
generated using the FastPM algorithm (Feng et al., 2016) imposing a linear
matter power spectrum in a $\Lambda$CDM Universe with cosmological parameter
taken from Planck Collaboration et al. (2016). The UNITsim provides three
different volumes. We use a cosmological volume of
$1(\mathrm{Gpc}\,h^{-1})^{3}$, run with $4096^{3}$ dark matter particles with
mass resolution of $1.2\times 10^{9}M_{\odot}h^{-1}$. The dark matter field is
represented by a mesh of $2048^{3}$ cells, which we have reduced to $512^{3}$,
yielding a spatial resolution of $\sim 1.95$Mpc$h^{-1}$.

Figure 1: Slices of $\sim 80$ Mpc$h^{-1}$ thickness in the UNITsim (at $z=0$)
showing the dark matter density field (upper left), the local halo Mach
$\mathcal{M}_{5}$ number (upper right, see §3.2.1), the neighbour statistics
$\mathcal{D}_{5}$ (lower left, see §3.2.2) and the tidal anisotropy
$\mathcal{T}_{A}$ (see §3.3.2. ) Figure 2: Probability distribution of dark
matter density in the UNITsim for two redshifts, in the four cosmic-web types
defined in §3.3.1. As we have used a threshold $\lambda_{th}=0$ in the cosmic-
web classification, we identify under-densities with voids and overdensities
with knots.

Importantly, the UNITsim is a paired-fixed amplitude simulation (Angulo &
Pontzen, 2016). On the one hand, fixed-amplitude means that, while the phases
$\phi_{k}$ of a Fourier mode $\delta(\@vec{k})=A_{k}e^{i\phi_{k}}$ are
randomly distributed in the range $[0,2\pi)$, the amplitudes $A_{k}$ follow a
distribution of the form $\delta^{D}(A_{k}-\sqrt{VP_{0}(k)/(2\pi)^{3}})$,
instead of a Rayleigh distribution with mode $P_{0}(k)$, as it is assumed in
standard simulations based on Gaussian random field (GRF hereafter). Pairing,
on the other hand, means that one realization is in reality a set of two
simulations, in which the phases of the initial density fields are shifted by
an amount $\pi/2$ (or the initial overdensity differing by a minus sign). This
means that two paired simulations display the same initial power spectrum by
construction, with full lack of correlation towards large scales.
Gravitational evolution induces mode coupling, which translates into a
variance in the power spectrum that is shown to be much smaller than that
obtained from standard simulations with the same initial conditions, when
analyzed from a variety of tracers (Villaescusa-Navarro et al., 2018; Chuang
et al., 2019). This feature makes this type of simulations suitable for the
modeling of cosmological probes and the assessment of systematic errors in
large-scale structure analyses (see e.g. Garrison et al., 2018; Zhai et al.,
2019; Maksimova et al., 2021). It has also been demonstrated that these types
of simulations do not introduce bias (as compared to standard simulations) in
the one- or two-point statistics (Angulo & Pontzen, 2016; Villaescusa-Navarro
et al., 2018) nor in galaxy/halo properties and scaling relations (means and
scatter). The scatter in probes such as density distributions, void
abundances, and halo bias is still compatible with that observed in standard
simulations (see e.g. Villaescusa-Navarro et al., 2018; Klypin et al., 2020).
This point is of particular interest for this paper, as we will explore the
behavior of halo scaling relations and derive many of our conclusions based on
the scatter observed from one simulation out of the pair. All the calculations
shown in this paper are performed using the same simulation, and we have
verified that the same results are obtained using the paired one. Fixed-paired
simulations have been also shown to be suitable for the generation of halo
mock catalogs based on learning algorithms (Balaguera-Antolínez et al., 2023).

##  3 Halo and environmental properties: correlation and scaling relations

In this section, we present the different properties employed in the analysis
(intrinsic to halos, dark matter, and the environment), and we discuss the
intricate relations between them as a function of redshift. The information
displayed and discussed in this section is organized in terms of several
figures, which combine multiple aspects of the analysis and were hence
designed to be self-contained. Since we will refer to these figures back and
forth throughout the section, and in order to facilitate the reading, we will
first provide a brief description of the figures. This description will be
expanded with detail in the subsections. The following figures are presented:

  * •

Figure 1: the spatial distribution of dark matter and halo properties (defined
in §3.1) as a function of several properties. The examples of distributions
shown correspond to slices of $\sim 80$Mpc$h^{-1}$ in thickness, colour coded
by dark matter density field, represented by $\log(\delta_{\rm dm}+1)$ (where
$\delta_{\rm dm}$ is the overdensity field), the Mach number (to be defined in
§3.2.1), the neighbour statistics (to be defined in §3.2.2) and the tidal
anisotropy (to be defined in §3.3.2).

  * •

Figure 2: the distribution of dark matter overdensity for two different
UNITsim snapshots. This plot is included for guidance, as in forthcoming
sections we will describe the behavior of halo bias and halo properties in
high and low density regions, defined through a cosmic-web classification
presented in §3.3.1.

  * •

Figure 3: the distribution of halo properties at redshift $z=0$ for different
cosmic-web types (to be defined in §3.3.1).

  * •

Fig. 4: the mean scaling relations between halo mass and the rest of the halo
properties as a function of redshift. We show this rather standard plots since
there is a clear connection between secondary bias, halo scaling relations and
the location of halos in the cosmic-web types, as we discuss throughout this
paper.

  * •

Figure 5: the volume filling fraction (VFF) for different cosmic-web
environments (§3.3.1).

  * •

Figures 6 and 7: the Spearman’s rank correlation coefficients for halo
properties computed at different redshifts and cosmic-web types (§3.4).

  * •

Figure 8: the scaling relations (as a function of halo mass) at a fixed
redshift, but this time measured in different cosmic-web environments.

###  3.1 Intrinsic halo properties

Halo catalogs and their properties have been built using the ROCKSTAR halo
finding algorithm (Behroozi et al., 2013). The robustness of the halo
properties has been explored in a number of comparison projects (see e.g.
Evrard et al., 2014; Mansfield & Avestruz, 2021; Knebe et al., 2011).
Throughout this paper, we shall refer to _scaling relations_ (or link) between
halo properties. This term must be understood as the probability that a halo
has a property $a$ conditional to a property $b$ in the range $b,b+{\rm d}b$,
i.e, $\mathcal{P}(a|b){\rm d}b$.

####  3.1.1 Halo Mass

We have chosen the virial mass as mass proxy, obtained in the halo finder
using spherical overdensities of $\Delta\sim M_{\rm vir}R_{\rm
vir}^{-3}\bar{\rho}_{\rm mat}(z)$ (Bryan & Norman, 1998; Behroozi et al.,
2013). We select parent halos with a minimum host-halo mass (virial mass) of
$2\times 10^{11}M_{\odot}h^{-1}$, which corresponds to $\sim 300$ dark matter
particles for a minimum halo mass. The halo mass distribution within the
UNITsim is shown at $z=0$ in Fig.3.

####  3.1.2 Maximum Circular Velocity

Computed as $V_{\rm max}\equiv{\rm Max}(V_{c})$, where
$V_{c}(r)=\sqrt{GM(r)/r}$ is the halo circular velocity at a radius $r$ and
$M(r)$ is the mass enclosed within that radius. A simple version of the virial
theorem predicts $V_{\rm max}\propto M_{\rm vir}^{1/3}$, confirmed in the
UNITsim at different redshifts (see also Shaw et al., 2006; Cataldi et al.,
2023). We have also read the velocity dispersion $\sigma_{v}$, computed from
velocity statistics of dark matter in halos. As for $V_{\rm max}$, velocity
dispersion scales approximately as $\sigma_{\rm v}\propto M_{\rm vir}^{1/3}$
(see e.g. Ragone-Figueroa et al., 2010). Figure 3 contains the distribution of
$V_{\rm max}$ at $z=0$, while its scaling relation with mass is shown in
Fig.4.

####  3.1.3 Scale radius and halo concentration

The scale radius $R_{s}$ is provided by the halo-finder after fitting the halo
density profiles to a Navarro-Frenk-White profile (Navarro et al., 1996).
Fig.4 shows the scaling of this property (in kpc$/h$) with respect to halo
mass. Besides of shaping the density profile (it marks the scale at which the
logarithmic slope of the profile equals $2$), this scale has been related to
halo anisotropy parameters (see e.g. Binney & Tremaine, 2008; Barnes et al.,
2005).

We compute the halo concentration as $C_{\rm vir}\equiv R_{\rm vir}/R_{s}$,
using the virial radius $R_{\rm vir}$ provided by the halo finder (see also
Klypin et al., 2016). Halo concentration scales with halo mass (see e.g.
Ludlow et al., 2016; Child et al., 2018) and has been shown to be a good proxy
for halo formation time, at fixed halo mass (Wechsler et al., 2002; Ludlow et
al., 2016; Wang et al., 2020). Figure 4 displays the mean scaling relation
with halo mass, showing a positive slope at at high redshift, i.e, more
massive halos are more concentrated. The slope in this link changes towards
low redshift, fully inverted for redshifts below $z\sim 1$ in agreement with
previous results in the literature (see e.g. Prada et al., 2012; Diemer &
Joyce, 2019).

Figure 3: Distribution of intrinsic and environmental halo properties in the
UNITsim at $z=0$, evaluated in different cosmic-web types as defined in
§3.3.1.

####  3.1.4 Dimensionless spin parameter

ROCKSTAR computes dimensionless spin parameters using both the Peebles (1971)
definition, as well as the (Bullock et al., 2001) prescription, i.e.,
$\lambda_{b}=J/(\sqrt{2}M_{\rm vir}V_{c}(R_{\rm vir})R_{\rm vir})$ where $J$
is the total halo angular momentum. We use the latter in this analysis. Figure
3 presents the spin distribution for different cosmic web types, while Fig.4
displays its behavior with respect to halo mass at different redshifts.

####  3.1.5 The virial parameter

Defined as the ratio between kinetic to gravitational potential energy,
$\mathcal{V}=T/|U|$, the virial parameter is used as a proxy for the degree of
relaxation of a halo (Knebe & Power, 2008). Note that the virial theorem (see
e.g. Collins, 1978) predicts $\mathcal{V}\sim 1/2$ for relaxed isolated halos.
Since the internal dynamics of dark matter halos involve pressure terms (from
velocity anisotropies) and surface terms (inflow-outflow of material),
deviations from the value $\mathcal{V}\sim 0.5$ do not strictly imply a lack
of relaxation (Henriksen & Widrow, 1999; Bett et al., 2007; Klypin et al.,
2016). Halos in the UNITsim have $\mathcal{V}>1/2$.

####  3.1.6 Halo triaxiality and ellipticity

Ellipsoidal shape parameters are defined in terms of the ratio between the
halo semi-axes $(a,b,c)$ as

$\mathcal{T}_{h}\equiv\frac{1-q^{2}}{1-s^{2}},\,\,\,\,\,\,\,\,\mathcal{E}_{h}\equiv\frac{1-s^{2}}{1+s^{2}+q^{2}},$
(1)

where $q=b/a$ and $s=c/a$ ($a>b>c$) 333The set $(q,s)$ is originally provided
by the halo finder.. The triaxiality parameter $\mathcal{T}_{h}$ defines
ellipsoidal geometries such as the oblate ($\mathcal{T}_{h}\to 0$) and prolate
($\mathcal{T}_{h}\to 1$) spheroids. On the other hand, halo ellipticity
$\mathcal{E}_{h}$ has as special limits the spherical ($\mathcal{E}_{h}\to 0$)
and non-spherical (or aspherical) configurations (see e.g. Cole & Lacey, 1996;
Binney & Tremaine, 2008).

####  3.1.7 Peak height

We assign to each halo a value of peak height
$\nu(M,z)=\delta_{sc}/\sigma(M,z)$, where $\delta_{sc}=(3/20)(12\pi)^{2/3}$ is
the critical collapse threshold linearly extrapolated to redshift $z=0$ and

$\sigma^{2}(M,z)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}{\rm
d}kk^{2}P(k,z)|W(k;R)|^{2}$ (2)

is the variance in the matter distribution on scales $R\propto(M_{\rm
vir}/\bar{\rho})^{1/3}$. $P(k,z)$ represents the linear matter power spectrum
at the corresponding redshift (obtained from the power spectrum defining the
initial conditions of the UNITsim, scaled by the corresponding growth factor
(Heath, 1977)). $W(k;R)$ is a top-hat filter with scale $R$ and
$\bar{\rho}(z)$ is the mean matter density at redshift $z$. In this work, we
use $\log\nu$ as a primary halo property, that can be exchange with the halo
mass according its definition. The mass scale set by $\nu(M_{\star},z)=1$
defines the mass of collapsing halos at redshift $z$. At $z=0$, $M_{\star}\sim
3.1\times 10^{12}M_{\odot}h^{-1}$. The peak height $\nu$ is the natural
variable in which the peak-formalism can predict halo abundances and bias Mo &
White (1996); Sheth & Tormen (1999b). These quantities display approximately a
universal shape (independent of mass and redshift) when expressed in terms of
it (Tinker et al., 2008, 2010; More et al., 2011). Also, trends in scaling
relations (such as the Mass-concentration relation) can be well described with
the same functional form in terms of this parameter (see e.g. Diemer & Joyce,
2019). Importantly, it is common to assume that the redshift evolution of
secondary bias scales with $\nu(M,z)$ (so it vanishes at fixed $\nu(M,z)$
Wechsler et al. (see e.g. 2006); Gao & White (see e.g. 2007)). Peak height
$\nu(M,z)$ is, therefore, a cornerstone for theoretical models of halo bias,
including its secondary dependencies.

Figure 4: Mean scaling relations between halo virial mass and a number of halo
properties described in §3.1, at different redshifts. The last row is
dedicated to non local and environmental halo properties. The error bars
denote the error in the mean in each mass bin. Figure 5: Volume filling
fraction (VFF) and the fraction of tracers $f_{\rm tr}$ in each cosmic-web
type defined in §3.3.1 as a function of redshift.

###  3.2 Non-local halo properties

The following set of properties are computed based on the halo distribution on
small scales. Throughout this paper we refer to these as _non-local_ , as they
are computed based on the phase-space properties of tracers enclosed on
spheres of radius $R=5$Mpc $h^{-1}$ computed around each halo. As such, this
set provide a good characterization of small-scale clustering.

####  3.2.1 Relative Mach number

The Mach number $\mathcal{M}_{R}$, introduced in the cosmological context by
Ostriker & Suto (1990), is a measure of the _temperature_ of the cosmic fluid.
In its original definition (see also Ostriker & Suto, 1990; Suto et al., 1992;
Ma et al., 2012; Agarwal & Feldman, 2013; Meriot et al., 2022), this quantity
measures the amount of bulk motion of patches of the Universe (characterized
by a scale $R$) with respect to the velocity dispersion of the random motions
inside those regions. This parameter contains cosmological information, as it
can used to probe the shape of the matter power spectrum (Strauss et al.,
1998). It has also been shown to display correlations with galaxy and halo
properties (Nagamine et al., 2001; Meriot et al., 2022).

In this work, we define a _relative Mach number_ by computing the relative
velocity $\@vec{v}_{ij}=\@vec{v}_{i}-\@vec{v}_{j}$ of a dark matter halo
(labelled $i$) with peculiar velocity vector $\@vec{v}_{i}$, and its
neighbouring halos with velocities $\@vec{v}_{j}$ inside spheres of radius
$R$, namely:

$\mathcal{M}^{(i)}_{R}=\frac{\langle|\@vec{v}_{ij}|\rangle_{j,R}}{\sigma_{i}},$
(3)

where $\langle\rangle_{j,R}$ denotes averages inside spheres of radius $R$ and
$\sigma_{i}$ is the dispersion in the variable $|\@vec{v}_{ij}|$ inside those
spheres. As in the original definition of the Mach number, high (low) values
of $\mathcal{M}^{(i)}_{R}$ can be associated with warm
$\mathcal{M}^{(i)}_{R}\geq 1$ (cool, $\mathcal{M}^{(i)}_{R}<1$) halo
environments. Under the assumption of no velocity bias (i.e, dark matter halos
and dark matter particles share the same velocity field), this measures the
_local temperature_ of the dark matter density field. In Fig.3, we present the
distribution of local halo Mach number at $z=0$ for different cosmic-web
types. In addition, Fig. 4 displays the mean scaling relation between halo
mass and Mach number at different redshifts. The trends at various redshift
imply that tracers live in _warm_ environments, and the scaling relations can
in principle be described by two different power-laws with a break at mass
scales varying with redshift, as well at the overall amplitude (the high-
redshift population of halos are warmer than low-redshift one). In particular,
at $z=0$, the relative Mach number for masses below $\sim 10^{13}M_{\odot}/h$
varies little. The information encoded in Fig.2 reveals that the halo
population, up to such mass scale, is mainly living in filaments, above which
the fraction of halos living in knots dominate the halo sample. Hence, the
relative Mach number varies faster with mass in high density regions. This is
also seen in Fig.8, where the scaling relation is split in different cosmic-
web environments (to be defined in §3.3.1).

####  3.2.2 Neighbour statistics

Along with the relative Mach number previously defined, we introduce the
_neighbour statistics_ $\mathcal{D}_{R}$, as a probe for the statistics of
pair separations around each tracer. This estimator is computed as

$\mathcal{D}^{(i)}_{R}=\frac{\langle\lambda_{ij}\rangle_{j,R}}{\sigma_{i,\lambda}},$
(4)

where $\langle\lambda_{ij}\rangle_{j,R}$ is the mean separation of halos at a
position $\@vec{r}_{j}$ and within spheres of radius $R$ from the main tracer,
while $\sigma_{\lambda_{R}}$ is the corresponding population variance. The
neighbour statistics aims at compressing the two first moments of the
distribution of pair separations, and, as such, contains information on the
small-scale clustering of halos. Figure 3 shows the distribution of
$\mathcal{D}_{R}$ across the halo population at $z=0$ in different cosmic-web
environments, while Fig. 4 displays the mean scaling relation between halo
mass and $\mathcal{D}$ at different redshifts. Figure 8 shows $z=0$ scaling
relation as a function of cosmic-web environments.

Similar quantities have been explored in the context of secondary bias, such
as a _scale dependent correlation_ (Kerscher, 2018) or _neighbour distance_
(Salcedo et al., 2020).

####  3.2.3 Halo overdensity

We compute a _local halo overdensity $\Delta_{R}$_ by collecting the number of
tracers $N_{R}$ used to measure $\mathcal{M}_{R}$ (or $\mathcal{D}_{R}$) and
defining $\Delta_{R}=\log(N_{R}/\bar{N})$, where $\bar{N}$ is the expected
number of randomly distributed tracers in spheres of volume $\propto R^{3}$.

###  3.3 Environmental properties

By environmental properties, we refer to quantities computed from the
underlying dark matter field. These properties can be classified as local,
such as the dark matter density in a cell $\delta_{\rm dm}$444For clarity in
the presentation of some results, we use the label $\Delta_{\rm
dm}=\log(1+\delta_{\rm dm})$., and non-local, like the tidal field (see e.g.
Balaguera-Antolínez et al., 2020).

####  3.3.1 Cosmic-web classification

The cosmic web is a well recognized pattern (see e.g. Gott et al., 2005)
arising from the dynamical evolution of density fluctuations. Several methods
have been designed to define different environments in this cosmic structure
(see e.g. Schaap & van de Weygaert, 2000; González & Padilla, 2010; Aragón-
Calvo et al., 2010; Libeskind et al., 2018, and references therein). Here we
use the approach by Hahn et al. (2007) and compute the Hessian of the
gravitational potential (or tidal field) $\mathcal{T}_{ij}=\partial_{ij}\Phi$
(where $\Phi$ is the comoving gravitational potential obtained from Poisson’s
equation $\nabla^{2}\Phi=\delta_{\rm dm}$) 555We compute the tidal field
directly from the overdensity field $\delta_{\rm dm}(\textbf{r})$, as the
inverse Fourier transform $\mathcal{T}_{ij}=F^{-1}[(k_{i}k_{j})\delta_{\rm
dm}(\textbf{k})/k^{2}]$ using the FFTW algorithm (Frigo & Johnson, 1998). The
eigenvalues have been computed using the gsl public library (Contributors,
2010). and its eigenvalues $\lambda_{i},(i=1,2,3)$ (see e.g. Hahn et al.,
2007; van de Weygaert et al., 2009; Forero-Romero et al., 2009; Aragon-Calvo,
2016; Yang et al., 2017; Paranjape et al., 2018) satisfying
$\sum_{i}\lambda_{i}=\delta_{\rm dm}$. Using these eigenvalues, each cell is
classified as part of a knot ($\lambda_{1}>\lambda_{\rm th}$,
$\lambda_{2}>\lambda_{\rm th}$, $\lambda_{3}>\lambda_{\rm th}$), filament
($\lambda_{1}<\lambda_{\rm th}$, $\lambda_{2}>\lambda_{\rm th}$,
$\lambda_{3}>\lambda_{\rm th}$), sheet ($\lambda_{1}<\lambda_{\rm th}$,
$\lambda_{2}<\lambda_{\rm th}$, $\lambda_{3}>\lambda_{\rm th}$), or void
($\lambda_{1}<\lambda_{\rm th}$, $\lambda_{2}<\lambda_{\rm th}$ and
$\lambda_{3}<\lambda_{\rm th}$), where $\lambda_{\rm th}$ is an arbitrary
threshold. In this work, we use $\lambda_{\rm th}=0$, which allows us to
associate knots (voids) with over-(under-) dense regions. The relevance of the
cosmic-web classification resides in the dynamical properties of each
environment and the density these are characterized with (Zel’dovich, 1970;
Doroshkevich & Shandarin, 1978; Bond et al., 1996), as well as the connection
with perturbation theory (see e.g. Kitaura et al., 2022). Halo bias has been
observed to change within cosmic-web environment (see e.g. Yang et al., 2017;
Balaguera-Antolínez et al., 2019; Bonnaire et al., 2022), a claim we support
with this paper. In future work we shall also explore the dependencies of halo
properties with respect to the geometrical properties of peaks in the density
field (Peacock & Heavens, 1985), which has been shown to correlate with halo
bias as well (see e.g. Sinigaglia et al., 2021).

In Fig.5, we show the fraction of cells (for the fiducial resolution) in the
UNITsim volume (i.e., the volume filling fraction, VFF) classified in the
different cosmic-web environments, together with the fraction of tracers,
$f_{tr}$, which, at each redshift, are located in cells classified
accordingly. The effects of gravitational evolution and the formation of the
cosmic-web is condensed in this plot- As the Universe evolves, less cells are
classified as knots, while the fraction occupied by filaments increases
(filaments are, as expected, the dominant component across redshift). The
trends in the number of tracers follow that of the VVF, except for sheets:
while the fraction of cells classified as such decreases towards low redshift,
the fraction of tracers in such environments rises. This as a consequence of
the density regimes probed by sheets (see Fig.2) which reveals the interplay
between the growth of voids ($\delta_{\rm dm}<0$) and the formation of halos
in regions with $\delta_{\rm dm}>0$.

####  3.3.2 Tidal anisotropy

Based on the set of eigenvalues $\\{\lambda_{i}\\}$, the degree of anisotropy
of the tidal field can be quantified by computing geometrical quantities such
as the ellipticity $e=(\lambda_{1}-\lambda_{3})/2\delta_{\rm dm}$ or the
prolateness $p=(\lambda_{1}+\lambda_{3}-2\lambda_{2})/2\delta_{\rm dm}$. We
can also use the invariants of the tidal field, which have been shown to
correlate with halo statistics at the level of number counts (see e.g.
Balaguera-Antolínez et al., 2019; Kitaura et al., 2022). In this work, we
compute the so-called tidal anisotropy parameter

$\mathcal{T}_{A}\equiv\frac{\sqrt{(\lambda_{1}-\lambda_{2})^{2}+(\lambda_{1}-\lambda_{3})^{2}+(\lambda_{3}-\lambda_{2})^{2}}}{2(2+\lambda_{1}+\lambda_{2}+\lambda_{3})}.$
(5)

Notice that this definition differs from others often found in the literature
e.g. Hahn et al. (2009); Paranjape et al. (2018); Yan et al. (2013);
Ramakrishnan et al. (2019); Favole et al. (2022), where a tidal anisotropy is
computed based on the halo distribution and thus measures the level of
anisotropy of such field (smoothed on a certain scale). Also, notice that in
the denominator of the equation above we use a factor $2+\delta_{\rm dm}$
(instead of $1+\delta_{\rm dm}$), in order to avoid divergences. According to
the cosmic-web classification, high values of tidal anisotropy can be found in
filaments and sheets. Accordingly, high-mass halos, mainly found in knots, are
expected to be assigned a low value of anisotropy, as shown in Fig.4.

Figure 6: Spearman’s rank correlation coefficient $\rho_{s}$ between the halo
effective bias and a number of halo properties for two different redshifts.
Figure 7: Spearman’s rank correlation coefficient $\rho_{s}$ between the halo
effective bias and a number of halo properties at $z=0$. Upper diagonal shows
the correlation computed in knots, lower diagonal in voids.

###  3.4 Correlations among halo properties

We assess the level of correlation among the different halo properties using
the Spearman’s rank correlation coefficient. As an example, in Fig.6 we
present the correlation between all properties at high and low redshift ($z=5$
and 0, respectively). The highest correlations are found, as expected, among
properties probing the potential wells, as well as the correlation among
environmental properties. The transition from high to low redshift depicted in
this figure shows how the correlation between intrinsic and environmental
properties increases with decreasing redshift, a change that is made evident
for properties such as halo spin and concentration, as can be also seen from
Fig.4.

Halo properties display, in general, a different degree of variability across
the halo population. In appendix A, we perform a principal component analysis
across redshift to demonstrate that even though intrinsic properties linked to
the depth of the halo potential well (mass, velocity dispersion) contain most
of the variability, non-local and environmental properties can contain similar
levels of variability (understood here as information). In that line, halo
properties such as spin, concentration and geometry can be regarded as
_secondary_ due to their lower degree of variability. Statistically, such low
variability can be associated to a large scatter in their link with e.g. halo
mass, as can be inferred from Fig.6.

Figure 8: Mean scaling relations between halo virial mass and a number of halo
properties described in §3.1, evaluated in different cosmic-web environments
at $z=0$. Figure 9: Slice of thickness $\sim 80$ Mpc $h^{-1}$ of the UNITsim
(at $z=0$) showing the distribution of effective halo bias distribution,
obtained form the interpolation of the bias assigned to each halo according to
Eq.(9) on $512^{3}$ mesh.

Figure 10: Halo effective bias $b_{hm}$ computed with Eq.(9) at $z=0$. Left:
bias as a function of halo peak height; the points denote the mean bias in
different $\nu$ bins, and the error bars denote the standard error of the mean
in each bin. Dotted line shows the prediction of (Tinker et al., 2010). Right
panel shows the bias $b_{hm}$ as a function of halo virial mass. In both cases
the contours indicate region of equal number of tracers $N_{H}$.

###  3.5 Scaling relations conditional to cosmic-web environments

Based on the cosmic-web classification described in §3.3.1, we assess now the
level of correlation and the statistical behavior of the halo scaling
relations in different cosmological environments. Recall, first, that Fig.3
displays the distribution of halo properties (intrinsic and non-local)
evaluated in different cosmic-web environments, with the noticeable feature
that these distributions are much more sensitive to under-density regions
(voids). To complement these results, Figure 7 shows the Spearman’s rank
correlation coefficient $\rho_{s}$ at $z=0$ in extreme environments, i.e,
voids and knots, probing under- and over-dense regions, respectively. This
figure depicts again the two main regions in the parameter space (intrinsic
and environmental properties) displaying the largest correlations. The
variations in the correlation coefficient in these two cosmic-web types is a
signature that the links among halo properties can be sensitive to the
cosmological environment they live in. In Fig. 8 we preset the mean scaling
relations between halo mass and other properties at $z=0$, evaluated for
different cosmic-web types. Some conclusions of the main conclusions from
these results are:

  * •

Properties probing the depth of the halo potential well (e.g halo mass,
maximum circular velocity, or velocity dispersion) are linked through a rather
stable scaling relation against changes in the environment and across
redshift.

  * •

Statistically significant differences among scaling relations are found in
voids, specially for the spin-mass link, which not only displays a lower
amplitude (low mass halos in voids have lower spin than halos of the same mass
in higher density environments) but also a different trend (higher mass halos
in voids tend to have lower spin compared with halos of the same mass in knots
). This can be understood as a signature of the coupling between halo-spin and
tidal field, as suggested by the tidal-torque theory (see e.g. Doroshkevich,
1970; White, 1984; Porciani et al., 2002)(see also e.g. Codis et al., 2012;
Trowland et al., 2013; Aragon-Calvo & Yang, 2014, for an analysis of spin in
the cosmic-web). Spin has been shown not only to correlate with the cosmic
environments, but also with properties of its components such as the thickness
of filaments (Ganeshaiah Veena et al., 2021).

  * •

The virial $\mathcal{V}$ in voids is lower than in high density region. If we
assume that relaxation is represented as $T\to 2|U|$, this can be interpreted
as a signal that halos in voids are more relaxed as these tracers are less
prone to growing by merging.

  * •

Halos with masses up to $\sim 10^{12}M_{\odot}h^{-1}$ living in voids tend to
display higher ellipticity and triaxiality (i.e, halos tend to be more
prolate) than halos of the same mass in other environments. Above that mass
scale, ellipticity of halos becomes independent of the cosmic-web type. The
trend itself also changes: halos above that mass scale tend to be more oblate
than those in higher density regions.

  * •

The mean scaling relations between environmental properties and halo mass also
display (statistically) large differences among cosmic-web types, specially in
voids, with expected results for each of the shown properties. Namely, halos
in voids tend to live in a cooler environment (Mach number), with pairs at
larger separations (neighbour statistics) and in regions with a mean halo
density below the mean (local overdensity).

  * •

Environmental properties (local dark matter density and tidal anisotropy) show
large differences, which is a trivial consequence of their dependence on the
cosmic-web classification

These conclusions are in agreement with previous analysis (see e.g. Hellwing
et al., 2021), and highlight the importance of the environment in shaping the
link between different halo properties. Dependencies of halo scaling relations
might not only appear due to the environment, but also through the underlying
cosmological model (see e.g. Ho et al., 2018).

##  4 Halo bias across cosmic time

###  4.1 The large-scale halo bias

It is very well established that the bias of dark matter tracers needs to be
modelled beyond the linear scale-independent scheme (see e.g. Desjacques et
al., 2018, and references therein), in which the halo bias $b$ simply relates
dark matter $\delta_{\rm dm}$ with halo $\delta_{h}$ overdensities (determined
on a mesh characterized by some scale $R$) through an expression of the form
$\delta_{h}=b\delta_{\rm dm}$. Scale dependencies induced by the process of
halo formation, halo merging history, non-linear evolution of the dark matter
density field (see e.g. Matsubara, 1999; Sigad et al., 2000; Somerville et
al., 2001; Smith et al., 2007; Zentner, 2007; Tinker et al., 2010; Valageas,
2011; Pollack et al., 2012; Sheth et al., 2013; Baldauf et al., 2013; Ahn et
al., 2015; Han et al., 2019; Nasirudin et al., 2019, and references therein),
and the discrete presentation of dark matter halo and matter density fields
generalises the concept of halo bias to a non-local and stochastic quantity ,
characterized by a conditional probability distribution
$\mathcal{P}(\delta_{\rm h}|\delta_{\rm dm})_{R}$ (see e.g. Fry & Gaztanaga,
1993; Tegmark & Peebles, 1998; Dekel & Lahav, 1999; Blanton, 2000; Simon,
2005).

The effective large-scale halo bias, even if measured on scales where linear
perturbation theory applies, is not simply a link of the form
$\delta_{h}=b\delta_{\rm dm}$. Instead, it is the combination of terms
involved in the perturbative approach to describe the tracer density field,
(Heavens et al., 1998; McDonald & Roy, 2009; Desjacques et al., 2018). As
such, measurements of effective bias need to be complemented with other probes
(e.g. three-point statistics) in order to make connections with perturbation
theory (see e.g. Pollack et al., 2014; Werner & Porciani, 2020).

###  4.2 Estimators of effective bias

The large-scale effective halo bias can be measured via a number of methods,
such as number-count-in-cells , or two- and three-point statistics in Fourier
(see e.g. Kravtsov & Klypin, 1999; Smith et al., 2007; Tinker et al., 2010;
Pollack et al., 2012) and configuration space . These methods characterize the
clustering strength of a halo population either with respect to the underlying
dark matter (absolute bias) or to a reference population (relative bias) (see
e.g. Balaguera-Antolínez, 2014). In Fourier space, estimators of the large
scale bias can be based on auto-power spectrum, namely, ratio between halo
$P_{hh}(k)$ and dark matter power spectrum $P_{dm}(k)$:

$b_{hh}(k_{j})=\sqrt{\frac{P_{hh}(k_{j})}{P_{\rm dm}(k_{j})}},$ (6)

or based on the halo-dark matter cross-power spectra $P_{hm}(k)$:

$b_{hm}(k_{j})=\frac{P_{hm}(k_{j})}{P_{\rm dm}(k_{j})},$ (7)

the latter being free of shot-noise corrections and hence in principle more
suitable for sparse samples. In both cases, the estimate of effective large-
scale bias can be obtained by averaging the values of $b_{hm}(k_{j})$ (or
$b_{hh}(k_{j})$) over a range of wavenumbers $k_{j}<k_{max}$ in which the
ratio between the halo and the dark matter power spectra is constant, i.e.,

$b=\frac{\sum_{j}^{k_{j}<k_{max}}N_{k}^{j}b(k_{j})}{\sum_{j}^{k_{j}<k_{max}}N_{k}^{j}},$
(8)

where $N^{j}_{k}$ is the number of Fourier modes in the $j-$th spherical
shell. It is important to notice that for a fixed maximum wavenumber, these
estimators are not expected to provide the same large-scale bias, as non-
linear evolution affects differently the halo and dark matter field (see e.g.
Smith et al., 2007; Pollack et al., 2014).

Paranjape et al. (2018) used the bias estimator of Eq.(7) to write the
effective large-scale bias as a sample mean, obtained from a tracer catalog in
which each element contributes with an _object-by-object_ bias of the form

$b^{(i)}_{hm}=\frac{\sum_{j,k_{j}<k_{max}}N^{j}_{k}\langle{\rm
e}^{-i\textbf{k}\cdot\textbf{r}_{i}}\delta_{\mathrm{dm}}^{*}(\textbf{k})\rangle_{k_{j}}}{\sum_{j,k_{j}<k_{max}}N^{j}_{k}P_{\rm
dm}(k_{j})},$ (9)

where $\delta_{\mathrm{dm}}(\textbf{k})$ is the Fourier transform of the dark
matter density field (see also Ramakrishnan et al., 2019; Han et al., 2019;
Paranjape & Alam, 2020). The sum is carried over the range of wavenumbers in
which the ratio between the halo and the dark matter power spectra is
constant. We have taken such limit as $k\leq 0.08\,h$Mpc-1. This limit is
rather conservative for high redshift, where non-linear evolution of the halo
bias takes place at much smaller scales. At low-redshift ($z\sim 0$), the
ratio between halo and dark matter power spectrum is already scale-dependent
at $k\sim 0.1\,h$ Mpc-1.

The large-scale effective bias of a subsample of $N_{h}(\Delta\theta)$ tracers
in a bin of property $\Delta\theta$ can be obtained as a sample mean:

$\langle b_{hm}|\Delta\theta\rangle=\frac{1}{N_{\rm
h}(\Delta\theta)}\sum_{i\in\Delta\theta}^{N_{h}(\Delta\theta)}b^{(i)}_{hm}.$
(10)

where $N_{h}(\Delta\theta)$ denotes the number of halos in the $\Delta\theta$
bin. The effective bias estimator of Eq. (9) has a number of advantages, some
of which will be mined in this paper. The first of them is that, by assigning
a bias to each object, we can statistically analyze the effective bias in
terms of different (primary or secondary) halo properties without the explicit
need to divide the sample in bins and compute two-point statistics
(correlation function or power spectra). It also opens the possibility to
perform tests on models of redshift space distortion and applications in the
reconstruction of halo properties for mock catalogs (Balaguera-Antolínez et
al, in prep). Importantly, Eq. (9) has been already implemented in the context
of both halo and galaxy assembly bias (see e.g. Paranjape et al., 2018;
Ramakrishnan et al., 2019; Contreras et al., 2021).

Figure 9 shows a slice of thickness $\sim 80$ Mpc $h^{-1}$ through the UNITsim
at $z=0$, with the color code indicating the average (within the slice shown)
of the bias-weighted halo density field. It is interesting to notice how a
threshold in bias can define relatively regularly distributed structures with
typical sizes of $\lesssim 100$ Mpc. In a forthcoming paper we will explore
the possibility of using thresholds of individual bias to characterize large-
scale structures on scales larger than those of galaxy clusters.

Figure 11: Effective halo bias in the UNITSim at $z=0$, measured in bins of
halo viral mass. Symbols denote the mean halo-mass bias $\langle b_{hm}|\log
M_{\rm vir}\rangle$ measured with Eq.(9), the halo-to-matter ratio $b_{hh}$
and the cross power spectrum $b_{hm}$. The dashed lines corresponds to the
prediction of (Tinker et al., 2010), computed with an spherical overdensity at
virialization of $\Delta_{\rm vir}=330$. Figure 12: Mean halo bias as a
function of peak height at different redshifts in the UNITsim. The dashed line
shows the fitting formula of (Tinker et al., 2010). Error bars denote the
error of the mean in each mass bin.

Figure 10 shows the scaling relation between the effective halo bias computed
with Eq.(9) and both peak height, $\nu$ (left panel), and halo mass (right
panel). In both panels, we show the fitting formula of Tinker et al. (2010,
T10), which depends on the spherical overdensity used to define dark matter
halos, $\Delta_{\rm vir}$. We have employed $\Delta\sim M_{\rm vir}R_{\rm
vir}^{-3}\bar{\rho}_{\rm mat}(z)\sim 331$ (at $z=0$) (Bryan & Norman, 1998),
which correspond to the mean spherical overdensity of the halo population in
the UNITsim (Behroozi et al., 2013). The differences observed with respect to
the prediction from T10 are expected, as the latter are generated using a
different definition of halo mass (see e.g. White, 2001; More et al., 2011;
Balaguera-Antolínez et al., 2012).

We have compared the estimates of effective halo bias from the estimators of
Eq. (6) (7) and (9). This is presented in Fig. 11 666All power spectra have
been computed using a $512^{3}$ mesh with cloud-in-cell mass assignment
(Hockney & Eastwood, 1988) and its corresponding correction in Fourier space
and Poisson shot-noise subtraction. When averaging in Fourier bins, use
weighted averages with number of modes in each bin as the weights., showing
measurements obtained using Fourier modes to a maximum $k_{\rm max}={\rm
min}(0.08,k_{sn})\,h$Mpc-1, where $k_{sn}$ is the wavenumber at which the
signal of power spectrum equals the Poisson shot-noise. The results based on
Eq. (6) is in agreement with those based on Eq.(9), while the measurement of
$b_{hm}$ based on Eq. (7) differs by a factor $\sim 1\sigma$. 777Both the
estimates and error bars assigned to the bias obtained from the power-
spectrum-based estimators are computed as weighted averages using the number
of Fourier modes in each shell as weights. The overall agreement among these
estimates confirms the robustness of Eq. (9) in terms of measuring effective
halo bias. From now on, we shall refer to the effective halo bias with the
symbol $b_{h}$.

###  4.3 Characterization of the effective halo bias distribution

In Fig.12, we present the measurement of mean halo bias as a function of the
peak height for different redshifts. There is an approximate universality as a
function of redshift, spoiled only at high redshift (high $\nu$) with
deviations of the order of $1\sigma$ compared to T10. The universality of the
halo bias can be extended to the full bias distribution. To verify that,
notice that in Fig. 13 we present the estimates of the higher moments
(variance, skewness and kurtosis) of the $\mathcal{P}(b_{h}|\nu)$
distribution. These moments are (within the error bars) compatible with zero,
which added to the fact that the correlation between peak height and effective
bias is small ($\rho_{s}\leq 0.03$ at all redshifts explored, see Fig.6) leads
to the conclusion that effective halo bias follows a normal distribution, i.e,
$\mathcal{P}(b_{h}|\nu,z)\sim\mathcal{N}(\langle
b_{h}|\nu\rangle,\sigma_{\nu}(z))$, where the standard deviation
$\sigma_{\nu}(z)$ is independent of $\nu$ 888Note that we have determined that
the variance scales with redshift as
$\sigma_{\nu}(z)\approx(1+bz)^{c}\sigma_{\nu}(z=0)$ $b=1.4\pm 0.1$ and
$c=0.41\pm 0.06$, obtained using $\sim 20$ UNITsim snapshots. Figure 13 also
shows the higher moments computed in high-density regions (knots), with no
evident modifications with respect to the behavior in the full sample. We have
verified that this is also the case in other cosmic-web environments (i.e,
sheets, filaments). This behavior can complement the idea of an approximate
universality of halo bias when expressed in terms of the peak height, as the
main structure-formation formalisms predict.

Figure 13: Second and third order moments (variance, skewness and kurtosis) of
the halo scaling relation between effective bias and peak height
$P(b_{h}|\nu)$, computed for different redshifts. The right column shows the
behavior in high density regions (knots).

The scaling relations between halo bias and any other halo property $\theta$
(internal or external) can be expressed based on their link to peak height (or
halo mass) 999As the scaling relation between peak height and halo mass is
$\mathcal{P}(\nu|M_{\rm vir}){\rm d}M_{\rm vir}=\delta^{D}(\nu-\nu_{0}(M_{\rm
vir})){\rm d}M_{\rm vir}$ where $\nu_{0}(M_{\rm
vir})=\delta_{sc}/\sigma(M_{\rm vir})$ (see Eq. (2)), the statistical
properties of the halo bias as a function of halo mass are directly inherited
from those of peak height. using

$\mathcal{P}(b_{h}|\theta)=\int_{-\infty}^{\infty}\mathcal{P}(b_{h}|\nu)\mathcal{P}(\nu|\theta){\rm
d}\nu,$ (11)

such that the mean effective bias as a function of a property $\theta$ can be
written as

$\langle
b_{h}|\theta\rangle=\frac{\int_{-\infty}^{\infty}b_{h}\mathcal{P}(b_{h}|\theta){\rm
d}b_{h}}{\int_{-\infty}^{\infty}\mathcal{P}(b_{h}|\theta){\rm
d}b_{h}}=\frac{\int_{-\infty}^{\infty}\,\mathcal{P}(\nu)\mathcal{P}(\nu|\theta)\langle
b_{h}|\nu\rangle{\rm d}\nu}{\int{\rm
d}\nu\mathcal{P}(\nu)\mathcal{P}(\nu|\theta)},$ (12)

where $\mathcal{P}(\nu)=\int{\rm d}b_{h}\mathcal{P}(b_{h}|\nu)$ is the
distribution function of dark-matter halos as a function of their peak height.
Equations (11) and (12) (the latter is the continuous representation of
Eq.(9)) represent the link between large scale structure and astrophysics,
which is present at the level of dark matter halos and galaxy clusters: it
implies that in order to generate predictions on large-scale bias of a halo
population (containing cosmological information) as a function of a halo
property, we need to characterize the halo scaling relations, mainly shaped by
astrophysical processes (Cooray, 2006; Allen et al., 2011; Kravtsov & Borgani,
2012; Mantz et al., 2014; Balaguera-Antolínez, 2014; Evrard et al., 2014;
Singh et al., 2020). Notice that this is a symbolic separation, as the scaling
relations of halos contain cosmological information (Wechsler et al., 2002;
Ludlow et al., 2016; Diemer & Joyce, 2019)

Figure 14: Spearman’s rank correlation coefficient between the halo effective
bias and a number of halo properties as a function of redshift.

Figure 15: Spearman’s rank correlation coefficient between the halo effective
bias and a number of halo properties as a function of redshift.

Figure 16: Mean effective bias $\langle b_{h}|\theta,w_{i}\rangle$ measured in
bins of halo properties $\theta$ in different cosmic-web $\omega_{i}$ at
$z=0$.

Figure 17: Example of secondary halo bias: mean halo effective bias $\langle
b_{h}|M_{\rm vir}\rangle_{q}^{(\theta)}$ measured in bins of halo (log) mass
using a number of halo properties as secondary properties $\theta$. In each
mass bin, the sample has been divided in quartiles $q$ of the property
$\theta$. We show results from the lower (first) quartile (red symbols) and
the upper (fourth) quartile (blue error bars), along with the results from the
full sample (in each mass bin, yellow error bars), at redshifts $z=0$ (dashed
lines) and $z=2.2$ (solid lines).

Figure 18: Ratio between halo bias as a function of halo-mass (left column)
and peak-height (right column) obtained in two different quartiles of halo
concentration (top) and halo spin (bottom). Shaded area denotes the
uncertainty of the ratio computed from the error bars associated to each
effective bias.

A full characterization of the scaling relations between halo bias and halo
properties is beyond the scope of this paper (see e.g. Böhringer et al., 2012;
Rozo et al., 2014; Evrard et al., 2014, for a detailed analysis on halo and
cluster scaling relations). We nevertheless describe some relevant statistical
features, the first being the correlation among properties. In Fig.14, we show
the Spearman’s rank correlation coefficient 101010We have also computed the
Pearson correlation coefficient, leading to similar results. $\rho_{s}$
between the effective halo bias and a number of halo properties, as a function
of redshift. In general, the correlation coefficient is small, meaning that
for all properties, the link between effective halo bias and halo properties
cannot be described through a monotonic relation (as expected, given the large
scatter shown in e.g. Fig.10). Among the intrinsic properties, halo mass has
the largest correlation, except at low redshifts, where halo concentration
displays, interestingly, a larger correlation. The halo geometry,
characterized by the ellipticity $\mathcal{E}_{h}$, shows a significant anti-
correlation (with respect to intrinsic properties), similar in absolute value
to that displayed by mass.

The main message from this plot is that non-local and environmental properties
have higher-degree correlation with halo bias (compared to intrinsic
properties). While this might be expected for the non-local halo properties
(e.g. $\mathcal{D}_{5}$ and $\Delta_{5}$ are probes of small-scale clustering
and hence can be coupled via non-linear evolution with the large-scale modes
used to measure the effective bias), the behavior with respect to the
environmental properties (local density, tidal anisotropy) plots an scenario
in which halo bias is mainly driven by the environment where halos reside.
This will be key in order to explain the signal of secondary bias, as was
previously addressed by Paranjape et al. (2018); Ramakrishnan et al. (2019).

In Fig.15 we present the correlation coefficient $\rho_{s}$ for knots (high
density) and sheets (intermediate to low density) as a function of redshift.
The environmental properties display in general more correlation with bias in
high density regions, while the intrinsic halo properties mildly change from
one cosmic-environment to another. The behavior in voids is very noisy and
oscillates around $\rho_{s}=0$ for almost all properties thought the redshift
range.

Figure 19: Statistical significance of secondary halo bias as a function of
halo mass, computed using Eq. (13) for different redshifts and for a number of
secondary halo properties. The dark-shaded stripe denotes a strip of
$\Delta^{4-1}_{p-s}\pm 2$ and the horizontal line marks
$\Delta^{4-1}_{p-s}=0$. $\Delta^{4-1}_{p-s}>0$ ($<0$) implies that upper
(lower) quartiles display higher clustering strength than lower (upper)
quartiles.

Figure 20: Secondary halo bias in the cosmic-web. Mean halo bias in bins of
virial mass, with samples divided in quartiles of secondary properties and
evaluated in different cosmic–web environments, at $z=0$.

###  4.4 Effective halo bias and the cosmic web

In Fig.16, we show the measurements of mean effective bias as a function of
halo mass and cosmic-web types $\langle b_{h}|\theta,w_{i}\rangle$ at $z=0$.
Some conclusions based on these results are:

  * •

A significant statistically difference in the mean effective halo bias can be
observed as a function of comic-web types.

  * •

The mean relation between halo bias and intrinsic properties is positive for
all cosmic-web types, except for voids.

  * •

Non-local and environmental properties show negative bias not only in voids,
but also in filaments and knots. In general these set of properties display
the largest differences in their links with effective bias as a function of
the cosmic-web type.

  * •

The halo effective bias is mildly higher in environments probing intermediate
densities (i.e, filaments) than in extreme environments (knots and voids).

One key feature arising from these results is the fact that we can measure,
with statistical significance, a negative halo bias. A negative halo bias
should not be surprising, since it can arise from the comparison of halo
overdensities ($\delta_{h}>0$) in voids ($\delta_{\rm dm}<0$ as
$\lambda_{th}=0$). This does not implies that all objects living in voids
acquire a negative bias. In fact, nearly $40\%$ of the halos living in voids
have assigned a positive effective bias.

Figure 21: Significance of secondary bias (Eq.(13)), evaluated in different
cosmic-web types at $z=0$.

Figure 22: Comparison between the measurements of effective halo bias as a
function $V_{\rm max}$, $C_{\rm vir}$ and $\lambda_{B}$ as read from the
original catalog against the signal measured using after reassigning the same
set of properties using a scaling relation based solely on the virial mass.
Bottom panels show the ratio between the bias in bins of the “reassigned” to
that measured in bins of the “original” property.

###  4.5 Secondary halo bias

The signal of secondary halo bias is obtained by assessing the mean effective
halo bias in bins of a primary property, in which the sample has been divided
in quartiles of a secondary property (see e.g.Sato-Polito et al. 2019 for a
similar analysis, but using the correlation function computed from this type
of subsets). For simplicity, we concentrate here in 6 representative
properties: spin, concentration, ellipticity, Mach number, dark matter
density, and tidal anisotropy. Figure 17 shows effective halo bias as a
function of halo virial mass computed in the first (lower) and fourth (upper)
quartiles of all these properties at redshifts $z\sim 2,0$ for the entire
sample. Key features can be recognized from this type of plots, such as
inversions of trends for quartiles at different mass scales, and redshifts at
which such inversions take place. This plot illustrates a well-known result:
virtually all internal and external halo properties display a certain level of
secondary halo bias in some mass ranges.

As we aim at measuring signals of secondary bias for a great number of halo
properties in different cosmic-web environments and cosmological redshifts, we
need to condense the information shown in Fig. 17 in a simple but still
informative way. One possibility that has been employed in the literature is
to compute the the relative bias between effective halo bias in the first and
fourth quartiles of $\theta_{s}$. We present one example of this procedure in
Fig.18, which displays the secondary halo bias with respect to concentration
and spin, at three different redshifts. Here, some of the main features of the
signals are clearly visible. For concentration, there is a well-known
crossover of the signal at around $\log M_{\rm vir}[h^{-1}M_{\odot}]=13$, so
low-concentrated halos have higher (lower) bias than same-mass high-
concentrated halos above (below) this characteristic mass. Higher spin halos
are, nevertheless, always more tightly clustered than lower spin halos at
fixed halo mass. The inversion of the signal in the secondary bias induced by
halo spin that was previously measured (at $\log M_{\rm
vir}[h^{-1}M_{\odot}]\simeq 11.5$, Sato-Polito et al. 2019; Johnson et al.
2019) is not observed in UNITSim when the entire sample is analyzed (it is
found only for certain cosmic-web environments, see below). This inversion was
shown to be caused by the so-called _splashback halos_ in Tucci et al. (2021),
at a mass scale very close to the minimum mass adopted for this work (see e.g.
Figure $2$ in Tucci et al. (2021)). We notice however that the relative bias
tends to unity for low-mass halos at $z=0$.

Importantly, the right-hand panels of Fig.18 display the secondary halo bias
as a function of peak height. At least to a first approximation, these plots
show that the redshift dependence of these particular signals is approximately
captured by $\nu$, as previously shown by, e.g. Wechsler et al. (2006); Gao &
White (2007).

We can also assess the significance of the signal of secondary bias by
computing

$\Delta^{i-j}_{p-s}\equiv\frac{\langle
b_{h}|\theta_{p}\rangle^{(s)}_{i}-\langle
b_{h}|\theta_{p}\rangle^{(s)}_{j}}{\sqrt{\sigma^{2}_{s,i}+\sigma^{2}_{s,j}}},$
(13)

where $\langle b_{h}|\theta_{p}\rangle^{(s)}_{i}$ is the mean sample bias
measured in bins of a primary property $\theta_{p}$ in the $i$-th quartile of
secondary property $\theta_{s}$. Similarly, $\sigma_{s,i}$ denotes the
standard error in the mean bias obtained in the $i$th quartile of the
secondary property. For a pair of primary and secondary properties, a value
$|\Delta^{i-j}_{p-s}|\gg 1$ can be used as an statistical signature of
secondary halo bias. It is important not to mistake this statistic with the
amplitude of the secondary bias signal itself, i.e., a particular property
might display a large signal at a given mass range, but at low statistical
significance.

In Fig. 19, we show the behavior of the variable defined by Eq.(13) as a
function of halo mass, for the selected set of secondary properties and as a
function of redshift. We include an uncertainty region of
$\Delta^{i-j}_{p-s}=\pm 2$ to guide the reader. The information displayed in
Figs. 19 is expanded in Figs. 20 and 21, where we evaluate the dependence of
secondary bias on the cosmic web. Also, the redshift evolution of secondary
bias based on peak height can be found in the Appendix. This set of plots
represents the measurement of secondary bias as a function of redshift and
cosmic-web type, thus containing a vast amount of relevance information. Our
main conclusions can be summarized as follows:

Figure 23: Same as Fig.22 focused on $V_{\rm max}$. In each panel, different
dependencies are used to reassign maximum circular velocity.

  * •

_Halo concentration_. The secondary bias based on concentration displays the
well-known trend explained above, dominated by a inversion of the signal at
high masses at $z=0$. Due to the intrinsic dependence of secondary bias on
$\nu$, this crossover moves towards lower halo masses at higher redshift.
Importantly, the significance of the concentration-based secondary bias
depends on cosmic-web environment: it is quite significant for knots and
filaments, but not so much for sheets and voids. Also, the ”crossover mass”
changes for different cosmic web types.

  * •

_Halo spin_. As mentioned before, spin bias conserves the link between
clustering strength and quartiles, i.e, upper quartile has a larger effective
bias at all mass scales in the full sample at $z=0$. It is only when the
signal is analyzed as a function of cosmic-web environments that an inversion
of the signal is observed for knots at low redshift (lower spin halos become
more highly clustered below $\log M_{\rm vir}[h^{-1}M_{\odot}]=12.2$ at
$z=0$). This might be due to the fact that splashback halos tend to be close
to knots in UNITSim.

  * •

_Halo ellipticity_. At fixed halo mass, spherical halos (lower quartile in
$\mathcal{E}_{h}$) are more clustered than non-spherical halos at $z=0$. This
trend appears to be inverted at $z\gtrsim 2$ at a mass scale of $\log M_{\rm
vir}\sim 13.2-13.4$. The $z=0$ result is in good agreement with Faltenbacher &
White (2010), where the redshift evolution was not measured. The secondary
bias signal for this geometry indicator is more significant than that of halo
triaxiality. We have checked that, at fixed halo mass, oblate halos are more
clustered than prolate objects, also in agreement with Faltenbacher & White
(2010). Again, the secondary bias signals produced by these parameters are
more significant in knots and filaments.

  * •

_Relative Mach number_. This property displays a significant positive
secondary bias signal, at all redshifts, particularly towards the lower halo
masses: at fixed halo mass, halos with higher Mach number have higher bias.
However, an inversion of the signal is detected, albeit with low significance,
for high-mass halos at lower redshift (i.e, low-$\mathcal{M}_{5}$ halos are
more clustered). In terms of cosmic-web types, again the presence of an
inversion depends on environment (e.g. it is not observed clearly in sheets).

  * •

_Neighbour statistics, local density, local halo overdensity, tidal
anisotropy_. All these halo external properties tend to produce very large and
significant positive secondary bias signal (larger values at fixed halo mass
are associated with higher bias). Again, the characteristic evolution dictated
by $\nu$ is partially observed (with the signals shifting towards lower
masses). There are particular redshifts and cosmic-web environments for which
small and low (statistical) - significance inversions are detected.

  * •

Among all properties considered, the one that produces the most statistically
significant low-mass secondary bias signal are, in this order, tidal
anisotropy, DM density, Mach number, and either ellipticity or concentration
(at a reference mass of $\log M_{\rm vir}[M_{\odot}h^{-1}]\sim 12.2$).

  * •

Secondary halo properties tend to produce secondary bias signals that decrease
in significance with halo mass (which could be partially due to the large
uncertainty at the high-mas end of the mass function). At low redshift, spin
remains the property that produces the largest signal in UNITSim.

  * •

The redshift evolution of secondary bias is approximately (but not completely)
captured by the implicit redshift dependence of $\nu$ for internal halo
properties (particularly concentration and ellipticity). That is not the case
for the external halo properties, which still maintain significant evolution
even when $\nu$ is chosen as a primary property (see also Fig. 32 from the
Appendix).

Figure 24: Same as Fig.22 using more information in the scaling relation to
assign the property $\tilde{\theta}_{s}$ (halo concentration in this case) as
shown in each panel.

Figure 25: Same as Fig.22 using more information in the scaling relation to
assign the property $\tilde{\theta}_{s}$ (halo spin in this case) as shown in
each panel.

##  5 Reconstructing the signal of secondary bias

In the previous section, we showed that secondary bias is a statistically
significant effect in the context of halo clustering, as has been abundantly
reported in the literature. In this section, we aim at a detailed statistical
description of the signal.

Let us start by exploring the link between halo properties and halo bias. In §
4, we expressed halo effective bias as a function of a halo property
$\theta_{a}$ (different from halo mass or peak height) as

$\mathcal{P}(b_{h}|\theta_{a})=\int_{0}^{\infty}\mathcal{P}(b_{h}|M_{\rm
vir})\mathcal{P}(M_{\rm vir}|\theta_{a}){\rm d}M_{\rm vir},$ (14)

which implies that the statistical properties of halo bias as a function of
$\theta_{a}$ are the result of the interplay between the secondary-primary
property link and the bias-primary property connection (in our case, halo mass
is taken as the primary property). As pointed out in §4, the latter can be
readily obtained from a more fundamental property such as the peak-height as
$\mathcal{P}(b_{h}|M_{\rm vir})=\mathcal{P}(b_{h}|\nu=\nu_{0}(M_{\rm vir}))$,
according to Eq.(2).

Each of the conditional probability distributions under the integral in Eq.
(14) is a marginalization over the halo properties; denoting with
$\eta=\\{\eta_{i}\\}$ the remaining set of halo properties (i.e, not including
mass or $\theta_{a}$), this is expressed as

$\mathcal{P}(M_{\rm vir}|\theta_{a})=\int\,\mathcal{P}(M_{\rm
vir}|\theta_{a},\eta_{1},\cdots,\eta_{i}){\rm d}\eta_{1}\cdots{\rm
d}\eta_{i},$ (15)

where $\mathcal{P}(M_{\rm vir}|\theta_{a},\eta_{1},\cdots,\eta_{i})$
represents a multivariate distribution of the halo mass conditional to the
properties $\\{\eta\\}$ (see e.g. Evrard et al., 2014). In general, the
moments of the scaling relation $\mathcal{P}(M_{\rm vir}|\theta_{a})$ will be
dependent of the correlation between viral mass and the set $\\{\eta\\}$. In
particular, the mean and the scatter around it are a consequence of the
propagation of the statistical properties of the links $\mathcal{P}(M_{\rm
vir}|\eta_{i})$.

Figure 26: Reconstruction of secondary bias: mean effective halo bias measured
in bins of halo mass using quartiles in halo concentration (upper panels) and
spin (bottom panels). The signal is measured from the original set of halo
properties and those re-assigned using Eqs.(16) with different combinations of
non-local or environmental properties, including the information of tidal
field $\mathcal{T}$ through the cosmic-web classification.

Figure 27: Same as Fig.26 for halo spin as secondary property.

We now evaluate the impact of partial marginalization over a reduced number of
secondary properties in the scaling relations involved in Eq.(14). To simplify
the analysis, we use the results of Figs. 6 and 14, from which we have
concluded that the correlation between effective halo bias and halo properties
is in general smaller than the correlation between halo properties. Based on
this, we fix the scaling relation $\mathcal{P}(b_{h}|M_{\rm vir})$ to the
information provided by the $N$-body simulation, such that the impact of
secondary properties (including environmental properties) in the
$\mathcal{P}(b_{h}|\theta_{a})$ link will be addressed through the scaling
relation $\mathcal{P}(M_{\rm vir}|\theta_{a})$.

We measure the scaling relation $\mathcal{P}(\theta_{a}|M_{\rm
vir},\\{\psi\\})$, where $\\{\psi\\}\in\\{\eta\\}$ is a subset of secondary
properties111111We define a mesh of $N_{b}^{p}$ bins, where $p$ denotes the
number of properties used to characterize the scaling relation (i.e, for a
single mass dependency of a property $\theta$, $p=2$) and $N_{b}=100$, the
same number of bins for each property. We count the tracers in each property
bin and normalize it in bins of $\theta$.. We use this measurement to reassign
a new version of the property $\theta_{a}$ (labeled $\tilde{\theta}_{a}$) to
each tracer of the simulated volume, i.e

$\tilde{\theta}^{i}_{a}\curvearrowleft\mathcal{P}(\theta_{a}|M_{\rm
vir}=M^{i}_{\rm vir},\psi=\psi^{(i)}).$ (16)

By construction, the new property $\tilde{\theta}_{a}$ obeys the same scaling
relation with respect to the halo mass and the set $\\{\psi\\}$, while the
correlations between $\tilde{\theta}_{a}$ and the rest of the halo properties
(i.e, the set $\\{\eta\\}-\\{\psi\\}$ ) are partially inherited through the
link between those properties and the halo mass. Again, it is key to notice
that _even if the scatter is the same, it does not contain the full
information of the scaling relations marginalized over in Eq.( 15)._

We compare the properties of $\mathcal{P}(b_{h}|\tilde{\theta}_{a})$ (in
particular, the mean) against those from the original scaling relation
$\mathcal{P}(b_{h}|\theta_{a})$ 121212This procedure is similar to the _galaxy
shuffling_ performed in the context of galaxy assembly bias (see e.g. Croton
et al., 2007; Contreras et al., 2019); it aims at removing secondary
correlations in order to assess their impact. Figure 22 shows the result of
reassigning the each of the properties $\theta_{s}=\\{V_{\rm max},C_{\rm
vir},\lambda_{B}\\}$ using only the virial mass in Eq. (16). The resulting
mean effective bias displays large deviations with respect to the true one,
especially for halo concentration and spin. Interestingly, $\sim 20\%$
deviations in the bias measured in bins of $V_{\rm max}$ with respect to the
original behavior leads to the conclusion that $V_{\rm max}$, despite of being
a probe of the potential well of dark matter halos, (and a relatively small
scatter against halo mass compared with other properties) has secondary non-
negligible correlations with other properties, specially at is low-end (low
masses). This can be due to the high degree of anisotropy of the environment
(see e.g. Fig. 4) at those mass-scales, which can be represented by tidal
forces disturbing the gravitational potential and thereby the _virial_
predictions for $V_{\rm max}$. On the other hand, the results for the halo
concentration indicate that virial mass is far from being its main primary
property, while halo spin is more receptive to a single-mass dependence, but
still with large variations with respect to the true halo bias.

We can extend the set $\\{\psi\\}$ to other properties beyond halo mass. In
that regard, the correlation analysis of §3.4 and the principal component
analysis described in appendix §A indicate that environmental properties
display large variability in the parameter space defined by the halo
properties, and as such, these are likely to be the main driver for halo bias
(primary and secondary). Inspired by this, we have generated a number new halo
catalogs with newly assigned properties using variations of the set
$\psi=\\{\mathcal{M}_{5},\Delta_{\rm
dm},\Delta_{5},T/U,\mathcal{E}_{h},\omega_{i}\\}$ (again, keeping halo
effective bias and halo virial mass as their original values). The resulting
mean halo bias measured as a function of $\\{V_{\rm max},C_{\rm
vir},\lambda_{B}\\}$ reassigned under different configurations (sets of
$\\{\psi\\}$) are shown in Figs.23, 24 and 25.

In general, these tests indicate that the signal of halo effective bias,
measured as a function of halo properties different of halo mass (e.g. spin,
concentration, or X-ray luminosity and temperature for galaxy clusters), can
be obtained without systematic bias when non-local or environmental properties
are included in the link between the observable and a primary property. Such
extra-properties depend on the halo property used as observable. For example,
the test shown in these figures reveals that the inclusion of the viral
$\mathcal{V}$ and halo ellipticity $\mathcal{E}_{h}$ can be key to retrieve
the effective bias as a function of $V_{\rm max}$, while for halo spin, the
local halo overdensity is more relevant to obtain unbiased results.

In Figs.26 and 27 we show the measurements of secondary bias (using halo spin
and concentration as secondary properties) based on reassignments of halo
properties previously described, again at $z=0$. Left panels (top and bottom)
in this plot show that _the single mass-dependency use to sample the secondary
property following Eq. ( 16) cannot account for the signal of secondary bias,
even if partial information is present in the form of scatter_. Explicitly
including more information such as local halo overdensity $\Delta_{5}$ or the
cosmic-web classification can account help to bring the set of reconstructed
properties to the true signal of secondary bias, as shown in the third column
of the same plots. Again, the extra-secondary properties that help to improve
this signal depends on the correlation between these and the secondary
property under scrutiny, i.e, the set can vary e.g. between spin and
concentration.

Figure 28: Effective halo bias measured in a synthetic halo catalog with bias
and mass reassigned according to the learn-and-map approach (§6.1). The
assignment of halo bias used the set of environmental properties
$\\{\Delta_{5},\Delta_{\rm
dm},\mathcal{M}_{5},\mathcal{D}_{5},\,\mathcal{T}_{A},\omega_{i}\\}$, while
halo masses are drawn from the scaling relation $\mathcal{P}(M_{\rm
vir}|b_{h},\Delta_{\rm dm},\Delta_{5},\omega_{i})$.

##  6 Forthcoming work

As commented in the introduction, this paper represents the first in a series
of documents dedicated to the understanding of secondary bias and cosmological
analysis based on the analysis of the link between halo properties and its
large scale structure. Let us here briefly comment on forthcoming work related
to the large-scale effective halo bias and its link to halo scaling relations.

###  6.1 Mock catalogs

One application of the possibility to have individual halo bias assigned to
each tracers finds place in the generation of halo mock catalogs based on
calibrated methods (see e.g. Bond & Myers, 1996; Avila et al., 2015; Chuang et
al., 2015; Kitaura et al., 2016). In particular, the bias assignment method
(BAM) (see e.g. Balaguera-Antolínez et al., 2019; Balaguera-Antolínez et al.,
2020) has been designed to generate precise halo catalogs (in terms of two and
three point statistics) while preserving the large-scale clustering signal as
a function of halo properties. While different approaches can be followed to
assign halo properties based on the features of the underlying dark matter
(Zhao et al., 2015; Balaguera-Antolínez et al., 2023), the assignment of halo
properties constraint to respect the a reference large-scale bias remains a
difficult task, as such procedure in principle can be done mapping intrinsic
scaling relations which do not carry clustering information. In that regard,
the inclusion of small scale proxies for halo clustering such as $\Delta_{5}$
or the relative Mach number $\mathcal{M}_{5}$ can help improve the accuracy of
the reconstruction of halo fields.

Figure 29: Slices of $10$ Mpc$h^{-1}$ thickness in an area of $\sim
200\,($Mpc$h^{-1})^{2}$ of the UNITsim volume, showing the bias, relative Mach
number and density filed (at $z=0$) interpolated on a $512^{3}$ mesh.

With the possibility to assign individual bias to tracers, we can select a
number of environmental properties which, according e.g. to a correlation
analysis of § 3.4 (of Fig.14), are most correlated with the object-by-object
halo bias. Let us give a brief description of how it would work. Labeling this
set as $\\{\eta\\}$, we can measure (learn) the scaling relation
$\mathcal{P}(b_{h}|\\{\eta\\})$ and re-assign to each tracer a new bias
$\tilde{b}\curvearrowleft\mathcal{P}(b_{h}|\\{\eta\\}=\\{\eta\\}_{new})$. For
practical applications, the set $\\{\eta\\}_{new}$ will represent non local
properties evaluated in independent realization of dark matter filed on top of
which halos are to be mapped (Balaguera-Antolínez et al., 2019). Once bias is
assigned, halo properties would then be assigned in a learn-and-map procedure.
For example, virial mass can be assigned as

$\tilde{M}^{(i)}_{\rm vir}\curvearrowleft\mathcal{P}(M_{\rm
vir}|b_{h}=\tilde{b}_{h}^{i},\\{\eta\\}=\\{\eta\\}_{new}^{(i)}),$

after which $V_{\rm max}$ can be assigned using

$\tilde{V}^{(i)}_{\rm max}\curvearrowleft\mathcal{P}(V_{\rm
max}|\tilde{M}^{(i)}_{\rm vir}=M_{\rm
vir},\tilde{b}_{h}^{i}=b_{h},\\{\eta\\}_{new}^{(i)}=\\{\eta\\}),$

and so forth for other halo properties, following the order established by
e.g. the importance according to a principal component analysis (see appendix
§A) of the halo properties. As an example, in Fig.28 we show the halo
effective bias as a function of mass after having reassigned both halos and
masses following the learn-and-map prescription, with
$\\{\eta\\}=\\{\eta_{new}\\}=\\{\Delta_{5},\Delta_{\rm
dm},\mathcal{M}_{5},\mathcal{D}_{5},\,\mathcal{T}_{A},\omega_{i}\\}$ Contrary
to previous approaches (see e.g. Zhao et al., 2015; Balaguera-Antolínez et
al., 2023), this methodology can guarantee that the large-scale clustering of
a mock catalog, as a function of main halo properties, is in good agreement
with respect to the expected results (taken from a reference $N$-body
simulation.

Halo catalogs with precise properties provide high flexibility at the time of
producing galaxy mock catalogues using e.g. a halo occupation distribution
(HOD) framework (e.g. Cooray, 2002; Cooray & Sheth, 2002; Berlind & Weinberg,
2002; Kravtsov et al., 2004) the sub-halo abundance matching (SHAM) (see e.g.
Vale & Ostriker, 2004; Kravtsov et al., 2004; Conroy et al., 2006; Favole et
al., 2016) This is optimal for multi-tracer analyses (see e.g. Hamaus et al.,
2012; Abramo & Leonard, 2013; Abramo et al., 2016; Wang & Zhao, 2020; Zhao et
al., 2021) as expected to be performed in many experiments Levi et al. (2013),
EUCLID Amendola et al. (2016), J-PAS (Benitez et al., 2014), and the Nancy
Grace Roman Space telescope Spergel et al. (2015). The approach presented in
this paper opens the possibility to assign halo properties along with halo
bias, printing the signal of secondary bias in halo mock catalogs. This will
be explored in a forthcoming paper (Balaguera-Antolínez et al., in
preparation) part of this series.

###  6.2 Effective bias in redshift-space

The object-by-object estimator of effective bias in Eq.(9) can be generalized
to provide estimates of redshift-space effective halo-bias for each tracer.
Such generalization can be used to test different subjects linked to redshift
distortions, such as _velocity bias_ , _measurement of the growth of
structures_ and marked statistics in redshift space. This will be explored in
a forthcoming publication, also part of this series.

###  6.3 Applications to galaxy and cluster surveys

The analysis of the large-scale bias can be extrapolated to redshift surveys.
One can assign relative (to the full galaxy population) effective bias or
absolute effective bias to galaxies (as done e.g. by Contreras et al., 2021),
provided that underlying dark matter density field has been reconstructed
using a number of techniques techniques (see e.g. Hada & Eisenstein, 2018;
Jasche & Lavaux, 2019; Kitaura et al., 2021).

Similarly, in the context of galaxy clusters, as the relative bias of clusters
contains cosmological information (see e.g. Böhringer, 1999; Balaguera-
Antolínez et al., 2011; Pillepich et al., 2012; Clerc & Finoguenov, 2023;
Balaguera-Antolínez, 2014), we can assign relative effective bias using a
single measurement of cluster power spectrum from the the full cluster
population.

###  6.4 Halo effective bias as proxy for super-clusters

Inspection of Fig. 9 shows that thresholds in halo bias can be use to identify
large-scale structures (super-clusters) with a rather distinct spatial
distribution. In particular, from Fig. 29 we foresee the applicability of a
criteria based on individual bias to isolate super-structures. In that figure
we see how, while the relative Mach number follows the filamentary structure
designed by the dark matter density field, the halo bias isolates selects
regions of $\sim 200($Mpc$h^{-1})^{2}$ separated by vast filaments with
$b_{h}\sim 0$. In a forthcoming letter we will present the properties of such
extra-large scale structures.

###  6.5 Extensions to galaxy assembly bias

Some of the properties that we investigate in this work (such as the Mach
number) have not been explored with hydro-dynamical simulations (e.g.
Illustris TNG131313http://www.tng-project.org) in the context of galaxy
assembly bias (i.e., the secondary dependencies of galaxy clustering and
occupancy at fixed halo mass). Since the individual object-to-object bias has
been successfully implemented in TNG (Paranjape et al., 2018), our work opens
new routes to explore these effects. Also, the connection between galaxy
assembly bias and the cosmic web has only been recently addressed with hydro-
dynamical simulations (see e.g. Montero-Dorta & Rodriguez 2023), so there are
a number of aspects regarding the physical origins of the effect to be
understood.

##  7 Discussion and conclusions

In this paper, we have measured the signal of secondary halo bias for the
first time in the UNIT simulation (UNITsim), a fixed-amplitude and paired set
of $N$-body simulations following the evolution of dark matter and dark-matter
halos across cosmic time, since $z\sim 6$ to $z=0$, in a comoving box of $1$
Gpc$h^{-1}$ side. The primary and secondary dependencies of bias has been
obtained using a new approach to halo bias, in which every tracer can be
assigned a value of effective large-scale bias, based on the estimator
introduced by Paranjape et al. (2018). We have verified the robustness of this
estimator against the results obtained using estimators of large-scale bias
based on the measurements of two point-statistics, in particular the auto- and
cross-power spectra. The object-by-object assignment of halo bias opens the
possibility to analyze the statistical behavior of halo effective bias as a
function of a number of halo properties classified as _intrinsic_ , _non-
local_ and _environmental_.

Intrinsic properties (mass, velocity dispersion, halo spin) are derived from
the halo finder producing the catalogs of the UNITsim. These can be divided in
a set of _primary_ and _secondary_ halo properties. This classification is
based, on the one hand, on physical grounds (e.g. primary properties are those
probing directly the depth of the halo potential well, i.e mass, velocity
dispersion, maximum circular velocity), and on the other hand, on statistical
grounds (a principal component analysis, shown in Appendix A shows that these
contain higher level of variability, and therefore information, across time).

_Non-local_ halo properties refer to small-scale statistical diagnostics
aiming at condensing information on halo clustering on small scales. Among
these, we have introduced the _relative Mach number_ , which is a measure of
the kinetic _temperature_ “observed” by each tracer in spheres of $5$
Mpc$h^{-1}$. Similarly, we have used the _neighbour statistics_ , which
condenses information of the distribution of pair separations with respect to
each tracers in the same spheres.

Environmental properties are defined in this paper as those related to the
underlying dark-matter density field. This can be local (e.g local density)
and non-local (i.e, tidal field). The computation of the tidal field of the
dark-matter distribution allows for the exploration of halo bias as a function
of key properties such as its invariants or functions generated with its
eigenvalues. In particular, in this paper we have explored the behavior of
halo properties and effective bias as a function of cosmic-web environments
and the tidal anisotropy.

In the first part of this paper we have reported the measurements of halo
scaling relations in different cosmological environments. Conclusions can be
summarized as follows:

  * •

The halo scaling relations do not change in the same way through cosmic time.
While links between the so-called _primary_ properties (halo mass, velocity
dispersion, maximum circular velocity) vary approximately along with
expectations from e.g. the virial theorem, the scaling relations between mass
and _secondary properties_ are more complex (e.g these are not simple power
laws) and, in general, display larger scatter around the mean relations. One
particular example is the mass-concentration relation, whose slope varies with
mass and redshift, ranging from positive values at high-redshift (i.e, more
massive halos are more concentrated, with a Spearman’s rank correlation of
$\rho_{s}\sim 0.1$, see Fig.6), passing through a nearly mass-independent
relation at $z\sim 2$ (with $\rho_{s}\sim 0.04$) and ending with a negative
correlation (i.e, more massive halos are less concentrated) towards lower
redshift (with $\rho_{s}\sim-0.3$). This in agreement with a number of
findings in the literature

  * •

The variation of halo spin with halo mass at different redshifts is weak for
low mass halos ($\log M_{\rm vir}\sim 11.5$). In particular, at $z=0$ (with
$\rho_{s}\sim-0.01$), only the very massive halos ($\log M_{\rm vir}>14$)
display a relevant (negative) correlation with halo mass. The dependency is
slightly stronger at higher redshift (e.g. at $z\sim 5$, $\rho_{s}\sim-0.02$).

  * •

The halo-scaling relations involving secondary properties vary across cosmic-
web environment. Although the definition adopted based on the Hessian of the
gravitational potential is arbitrary (as it depends on a threshold set to zero
in this paper), the clear distinction of the scaling relations in over- and
under-densities indicates how the evolution and growth of halos is heavily
influenced by its environment (see Figs. 7 and Fig.8).

The second part of this paper has been devoted to the assessment of the
statistical properties of large-scale halo bias as a function of the
properties described in the first section of this work. Conclusions can be
summarized as follows:

  * •

The object-by-object halo assignment proposed by Paranjape et al. (2018) yield
estimates of large-scale halo bias in agreement with estimates based on halo
power spectrum (see Fig.10). This approach has the great advantage of no
explicit need to compute two-point statistics for sub-samples divided in bins
of primary and/or secondary properties.

  * •

The mean effective bias expressed as a function of the peak-height retains a
degree of universality not only in this mean, but also in higher moments (see
Fig.12). The probability distribution of halo effective bias as a function of
the peak height can be expressed as a normal distribution with a mean value in
good agreement with predictions from conjunction of $N$-body simulations and
the ellipsoidal collapse (Tinker et al., 2010), and a standard deviation that
scales only with redshift. The universality of the mean-effective bias is no
longer conserved when the halo sample is split in terms of cosmic-web
environments.

  * •

In general, effective halo bias expressed as a function of halo properties is
quite sensitive to the cosmic-web environment (see Fig.16).

  * •

The signal of secondary bias is present when exploring the bias as a function
of halo mass in quartiles of secondary properties (see Figs.17, 18 and 19), as
widely reported in the literature.

  * •

Most of the secondary halo properties induce secondary bias signals that
preserves a hierarchy across cosmic time, i.e, the difference between the bias
in the higher quartile and that in the lower quartile does not change sigh
with redshift, though the trend might differ from property to property. For
instance, high spin halos have larger bias than low-spin halos for all masses
and redshifts explored (§4.5)

  * •

Among the secondary properties, the signal of secondary bias due to halo
concentration displays an well-known behavior, inverting the hierarchy at
different mass-scale in different snapshot, a fact seen in previous works (see
e.g. Wechsler et al., 2006; Sato-Polito et al., 2019). Such behavior can be
seen as a direct consequence of the evolution of the mass-concentration
scaling relation which in turn is inherited from that of the scaling radius
$R_{s}$ (see Fig.4).

  * •

Given the mass resolution used in this work, we cannot detect the inversion of
clustering trend in the secondary bias induced by the halo spin, as reported
by (e.g. Sato-Polito et al., 2019; Tucci et al., 2021). However, such
crossover is observed when the analysis is conditioned to halos located in
high density regions (knots).

  * •

Non-local and environmental properties induce the most significant signals of
secondary bias. In particular, the tidal anisotropy, the relative Mach number
and the local dark matter density display the highest statistical significance
(see Fig.19).

  * •

The redshift evolution of secondary bias is approximately (but not completely)
captured by the implicit redshift dependence of $\nu$ for internal halo
properties (particularly concentration and ellipticity). That is not the case
for the external halo properties, which still maintain significant evolution
even when $\nu$ is chosen as a primary property.

We have demonstrated that the signal of secondary bias based on halo
properties is highly dependent on the environment where these halos reside. As
we have used the cosmic-web classification based on the tidal field, such
environments are linked high (e.g. knots), intermediate (filaments, sheets) or
low (voids) density regions. Extending the cosmic-web classification to other
dynamical properties of the dark-matter field (such as the sheer of the
velocity field or the local curvature) can provide more hints towards the
physical origin of secondary bias. We have mentioned a number of forthcoming
research in this context, involving the applicability to the generation of
mock catalogs with precise halo properties, redshift-space analysis and
applications to analysis of secondary bias in hydro-simulations.

###### Acknowledgements.

We thank Chia-Hsun Chuang, Gustavo Yépez and F-S Kitaura for making available
the dark matter and halo catalogs of the UNITsim. ABA acknowledges the IAC
facilities and the Spanish Ministry of Economy and Competitiveness (MINECO)
under the Severo Ochoa program SEV-2015-0548 grants. The UNITsim has been run
under at the MareNostrum Supercomputer hosted by the Barcelona Supercomputing
Center, Spain, with computing time granted by PRICE under grant number
2016163937. ADMD thanks Fondecyt for financial support through the Fondecyt
Regular 2021 grant 1210612. GF is supported by a Juan de la Cierva
Incorporación grant n. IJC2020-044343-I. No AI application has been used to
generate neither, text, data, code nor figures for this paper.

## References

  * Abramo & Leonard (2013) Abramo, L. R. & Leonard, K. E. 2013, MNRAS, 432, 318
  * Abramo et al. (2016) Abramo, L. R., Secco, L. F., & Loureiro, A. 2016, MNRAS, 455, 3871
  * Agarwal & Feldman (2013) Agarwal, S. & Feldman, H. A. 2013, MNRAS, 432, 307
  * Ahn et al. (2015) Ahn, K., Iliev, I. T., Shapiro, P. R., & Srisawat, C. 2015, MNRAS, 450, 1486
  * Allen et al. (2011) Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409
  * Amendola et al. (2016) Amendola, L., Appleby, S., Avgoustidis, A., et al. 2016, ArXiv e-prints
  * Angulo et al. (2008) Angulo, R. E., Baugh, C. M., & Lacey, C. G. 2008, MNRAS, 387, 921
  * Angulo & Pontzen (2016) Angulo, R. E. & Pontzen, A. 2016, MNRAS, 462, L1
  * Aragon-Calvo (2016) Aragon-Calvo, M. A. 2016, MNRAS, 455, 438
  * Aragón-Calvo et al. (2010) Aragón-Calvo, M. A., van de Weygaert, R., & Jones, B. J. T. 2010, MNRAS, 408, 2163
  * Aragon-Calvo & Yang (2014) Aragon-Calvo, M. A. & Yang, L. F. 2014, MNRAS, 440, L46
  * Avila et al. (2015) Avila, S., Murray, S. G., Knebe, A., et al. 2015, MNRAS, 450, 1856
  * Balaguera-Antolínez (2014) Balaguera-Antolínez, A. 2014, A&A, 563, A141
  * Balaguera-Antolínez et al. (2023) Balaguera-Antolínez, A., Kitaura, F.-S., Alam, S., et al. 2023, A&A, 673, A130
  * Balaguera-Antolínez et al. (2020) Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibáñez, M., et al. 2020, MNRAS, 491, 2565
  * Balaguera-Antolínez et al. (2019) Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibáñez, M., Zhao, C., & Abel, T. 2019, MNRAS, 483, L58
  * Balaguera-Antolínez et al. (2012) Balaguera-Antolínez, A., Sánchez, A. G., Böhringer, H., & Collins, C. 2012, MNRAS, 425, 2244
  * Balaguera-Antolínez et al. (2011) Balaguera-Antolínez, A., Sánchez, A. G., Böhringer, H., et al. 2011, MNRAS, 413, 386
  * Baldauf et al. (2013) Baldauf, T., Seljak, U., Smith, R. E., Hamaus, N., & Desjacques, V. 2013, Phys. Rev. D, 88, 083507
  * Barnes et al. (2005) Barnes, E. I., Williams, L. L. R., Babul, A., & Dalcanton, J. J. 2005, ApJ, 634, 775
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109
  * Benitez et al. (2014) Benitez, N., Dupke, R., Moles, M., et al. 2014, arXiv e-prints, arXiv:1403.5237
  * Berlind & Weinberg (2002) Berlind, A. A. & Weinberg, D. H. 2002, ApJ, 575, 587
  * Bett et al. (2007) Bett, P., Eke, V., Frenk, C. S., et al. 2007, MNRAS, 376, 215
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Blanton (2000) Blanton, M. 2000, ApJ, 544, 63
  * Böhringer (1999) Böhringer, H. 1999, Astronomische Nachrichten, 320, 189
  * Böhringer et al. (2012) Böhringer, H., Dolag, K., & Chon, G. 2012, A&A, 539, A120
  * Bond et al. (1996) Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603
  * Bond & Myers (1996) Bond, J. R. & Myers, S. T. 1996, ApJS, 103, 1
  * Bonnaire et al. (2022) Bonnaire, T., Aghanim, N., Kuruvilla, J., & Decelle, A. 2022, A&A, 661, A146
  * Borzyszkowski et al. (2017) Borzyszkowski, M., Porciani, C., Romano-Díaz, E., & Garaldi, E. 2017, MNRAS, 469, 594–611
  * Bryan & Norman (1998) Bryan, G. L. & Norman, M. L. 1998, The Astrophysical Journal, 495, 80
  * Bullock et al. (2001) Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240
  * Cataldi et al. (2023) Cataldi, P., Pedrosa, S. E., Tissera, P. B., et al. 2023, MNRAS, 523, 1919
  * Child et al. (2018) Child, H. L., Habib, S., Heitmann, K., et al. 2018, ApJ, 859, 55
  * Chuang et al. (2015) Chuang, C.-H., Kitaura, F.-S., Prada, F., Zhao, C., & Yepes, G. 2015, MNRAS, 446, 2621
  * Chuang et al. (2019) Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, MNRAS, 487, 48
  * Clerc & Finoguenov (2023) Clerc, N. & Finoguenov, A. 2023, in Handbook of X-ray and Gamma-ray Astrophysics. Edited by Cosimo Bambi and Andrea Santangelo, 123
  * Codis et al. (2012) Codis, S., Pichon, C., Devriendt, J., et al. 2012, MNRAS, 427, 3320
  * Cole & Lacey (1996) Cole, S. & Lacey, C. 1996, MNRAS, 281, 716
  * Collins (1978) Collins, G. W., I. 1978, The virial theorem in stellar astrophysics
  * Conroy et al. (2006) Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2006, ApJ, 647, 201
  * Contreras et al. (2021) Contreras, S., Angulo, R. E., & Zennaro, M. 2021, MNRAS, 504, 5205
  * Contreras et al. (2019) Contreras, S., Zehavi, I., Padilla, N., et al. 2019, MNRAS, 484, 1133
  * Contributors (2010) Contributors, G. P. 2010, The GNU Operating System
  * Cooray (2002) Cooray, A. 2002, ApJ, 576, L105
  * Cooray (2006) Cooray, A. 2006, MNRAS, 365, 842
  * Cooray & Sheth (2002) Cooray, A. & Sheth, R. 2002, Phys. Rep., 372, 1
  * Croton et al. (2007) Croton, D. J., Gao, L., & White, S. D. M. 2007, MNRAS, 374, 1303
  * Dalal et al. (2008a) Dalal, N., White, M., Bond, J. R., & Shirokov, A. 2008a, ApJ, 687, 12
  * Dalal et al. (2008b) Dalal, N., White, M., Bond, J. R., & Shirokov, A. 2008b, ApJ, 687, 12
  * Dekel & Lahav (1999) Dekel, A. & Lahav, O. 1999, ApJ, 520, 24
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep, 733, 1
  * Diemer & Joyce (2019) Diemer, B. & Joyce, M. 2019, ApJ, 871, 168
  * Doroshkevich (1970) Doroshkevich, A. G. 1970, Astrofizika, 6, 581
  * Doroshkevich & Shandarin (1978) Doroshkevich, A. G. & Shandarin, S. F. 1978, Sov. Ast., 22, 653
  * Evrard et al. (2014) Evrard, A. E., Arnault, P., Huterer, D., & Farahi, A. 2014, MNRAS, 441, 3562
  * Faltenbacher & White (2010) Faltenbacher, A. & White, S. D. M. 2010, ApJ, 708, 469
  * Favole et al. (2016) Favole, G., Comparat, J., Prada, F., et al. 2016, MNRAS, 461, 3421
  * Favole et al. (2022) Favole, G., Montero-Dorta, A. D., Artale, M. C., et al. 2022, MNRAS, 509, 1614
  * Feng et al. (2016) Feng, Y., Chu, M.-Y., Seljak, U., & McDonald, P. 2016, MNRAS, 463, 2273
  * Forero-Romero et al. (2009) Forero-Romero, J. E., Hoffman, Y., Gottlöber, S., Klypin, A., & Yepes, G. 2009, MNRAS, 396, 1815
  * Frigo & Johnson (1998) Frigo, M. & Johnson, S. G. 1998, in Proc. 1998 IEEE Intl. Conf. Acoustics Speech and Signal Processing, Vol. 3 (IEEE), 1381–1384
  * Fry & Gaztanaga (1993) Fry, J. N. & Gaztanaga, E. 1993, ApJ, 413, 447
  * Ganeshaiah Veena et al. (2021) Ganeshaiah Veena, P., Cautun, M., van de Weygaert, R., Tempel, E., & Frenk, C. S. 2021, MNRAS, 503, 2280
  * Gao et al. (2005) Gao, L., Springel, V., & White, S. D. M. 2005, MNRAS, 363, L66
  * Gao & White (2007) Gao, L. & White, S. D. M. 2007, MNRAS, 377, L5
  * Garrison et al. (2018) Garrison, L. H., Eisenstein, D. J., Ferrer, D., et al. 2018, ApJS, 236, 43
  * González & Padilla (2010) González, R. E. & Padilla, N. D. 2010, MNRAS, 407, 1449
  * Gott et al. (2005) Gott, J. R., Jurić, M., Schlegel, D., et al. 2005, The Astrophysical Journal, 624, 463
  * Hada & Eisenstein (2018) Hada, R. & Eisenstein, D. J. 2018, Monthly Notices of the Royal Astronomical Society, 478, 1866
  * Hadzhiyska et al. (2020) Hadzhiyska, B., Bose, S., Eisenstein, D., Hernquist, L., & Spergel, D. N. 2020, MNRAS, 493, 5506
  * Hahn et al. (2007) Hahn, O., Porciani, C., Carollo, C. M., & Dekel, A. 2007, MNRAS, 375, 489
  * Hahn et al. (2009) Hahn, O., Porciani, C., Dekel, A., & Carollo, C. M. 2009, MNRAS, 398, 1742
  * Hamaus et al. (2012) Hamaus, N., Seljak, U., & Desjacques, V. 2012, Phys. Rev. D, 86, 103513
  * Han et al. (2019) Han, J., Li, Y., Jing, Y., et al. 2019, MNRAS, 482, 1900
  * Hearin et al. (2016) Hearin, A. P., Zentner, A. R., van den Bosch, F. C., Campbell, D., & Tollerud, E. 2016, MNRAS, 460, 2552
  * Heath (1977) Heath, D. J. 1977, MNRAS, 179, 351
  * Heavens et al. (1998) Heavens, A. F., Matarrese, S., & Verde, L. 1998, MNRAS, 301, 797
  * Hellwing et al. (2021) Hellwing, W. A., Cautun, M., van de Weygaert, R., & Jones, B. T. 2021, Phys. Rev. D, 103, 063517
  * Henriksen & Widrow (1999) Henriksen, R. N. & Widrow, L. M. 1999, MNRAS, 302, 321
  * Ho et al. (2018) Ho, A., Gronke, M., Falck, B., & Mota, D. F. 2018, A&A, 619, A122
  * Hockney & Eastwood (1988) Hockney, R. W. & Eastwood, J. W. 1988, Computer simulation using particles
  * Jasche & Lavaux (2019) Jasche, J. & Lavaux, G. 2019, A&A, 625, A64
  * Jeeson-Daniel et al. (2011) Jeeson-Daniel, A., Dalla Vecchia, C., Haas, M. R., & Schaye, J. 2011, MNRAS, 415, L69
  * Johnson et al. (2019) Johnson, J. W., Maller, A. H., Berlind, A. A., Sinha, M., & Holley-Bockelmann, J. K. 2019, MNRAS, 486, 1156
  * Jolliffe & Cadima (2016) Jolliffe, I. & Cadima, J. 2016, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374, 20150202
  * Kerscher (2018) Kerscher, M. 2018, A&A, 615, A109
  * Kitaura et al. (2021) Kitaura, F.-S., Ata, M., Rodríguez-Torres, S. A., et al. 2021, MNRAS, 502, 3456
  * Kitaura et al. (2022) Kitaura, F.-S., Balaguera-Antolínez, A., Sinigaglia, F., & Pellejero-Ibáñez, M. 2022, Monthly Notices of the Royal Astronomical Society, 512, 2245
  * Kitaura et al. (2016) Kitaura, F.-S., Rodríguez-Torres, S., Chuang, C.-H., et al. 2016, MNRAS, 456, 4156
  * Klypin et al. (2020) Klypin, A., Prada, F., & Byun, J. 2020, MNRAS, 496, 3862
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340
  * Knebe et al. (2011) Knebe, A., Knollmann, S. R., Muldrew, S. I., et al. 2011, MNRAS, 415, 2293
  * Knebe & Power (2008) Knebe, A. & Power, C. 2008, ApJ, 678, 621
  * Kravtsov et al. (2004) Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35
  * Kravtsov & Borgani (2012) Kravtsov, A. V. & Borgani, S. 2012, ARA&A, 50, 353
  * Kravtsov & Klypin (1999) Kravtsov, A. V. & Klypin, A. A. 1999, ApJ, 520, 437
  * Lazeyras et al. (2017) Lazeyras, T., Musso, M., & Schmidt, F. 2017, J. Cosmology Astropart. Phys., 2017, 059
  * Lee et al. (2017) Lee, C. T., Primack, J. R., Behroozi, P., et al. 2017, MNRAS, 4, 3834
  * Levi et al. (2013) Levi, M., Bebek, C., Beers, T., et al. 2013, ArXiv e-prints
  * Li et al. (2008) Li, Y., Mo, H. J., & Gao, L. 2008, MNRAS, 389, 1419
  * Libeskind et al. (2018) Libeskind, N. I., van de Weygaert, R., Cautun, M., et al. 2018, MNRAS, 473, 1195
  * Ludlow et al. (2016) Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214
  * Ma et al. (2012) Ma, Y.-Z., Ostriker, J. P., & Zhao, G.-B. 2012, J. Cosmology Astropart. Phys., 2012, 026
  * Maksimova et al. (2021) Maksimova, N. A., Garrison, L. H., Eisenstein, D. J., et al. 2021, MNRAS, 508, 4017
  * Mansfield & Avestruz (2021) Mansfield, P. & Avestruz, C. 2021, MNRAS, 500, 3309
  * Mantz et al. (2014) Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2014, MNRAS, 440, 2077
  * Mao et al. (2018) Mao, Y.-Y., Zentner, A. R., & Wechsler, R. H. 2018, MNRAS, 474, 5143
  * Matsubara (1999) Matsubara, T. 1999, ApJ, 525, 543
  * McDonald & Roy (2009) McDonald, P. & Roy, A. 2009, J. Cosmology Astropart. Phys., 2009, 020
  * Meriot et al. (2022) Meriot, R., Khochfar, S., Oñorbe, J., & Smith, B. 2022, MNRAS, 512, 27
  * Mo & White (1996) Mo, H. J. & White, S. D. M. 1996, MNRAS, 282, 347
  * Montero-Dorta et al. (2020) Montero-Dorta, A. D., Artale, M. C., Abramo, L. R., et al. 2020, MNRAS, 496, 1182
  * Montero-Dorta et al. (2021) Montero-Dorta, A. D., Chaves-Montero, J., Artale, M. C., & Favole, G. 2021, MNRAS, 508, 940
  * Montero-Dorta & Rodriguez (2023) Montero-Dorta, A. D. & Rodriguez, F. 2023, arXiv e-prints, arXiv:2309.12401
  * More et al. (2011) More, S., Kravtsov, A. V., Dalal, N., & Gottlöber, S. 2011, ApJS, 195, 4
  * Musso et al. (2018) Musso, M., Cadiou, C., Pichon, C., et al. 2018, MNRAS, 476, 4877
  * Nagamine et al. (2001) Nagamine, K., Ostriker, J. P., & Cen, R. 2001, ApJ, 553, 513
  * Nasirudin et al. (2019) Nasirudin, A., Iliev, I. T., & Ahn, K. 2019, arXiv e-prints, arXiv:1910.12452
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
  * Ostriker & Suto (1990) Ostriker, J. P. & Suto, Y. 1990, ApJ, 348, 378
  * Pandey & Bharadwaj (2007) Pandey, B. & Bharadwaj, S. 2007, MNRAS, 377, L15
  * Paranjape & Alam (2020) Paranjape, A. & Alam, S. 2020, MNRAS, 495, 3233
  * Paranjape et al. (2018) Paranjape, A., Hahn, O., & Sheth, R. K. 2018, MNRAS, 476, 3631
  * Peacock & Heavens (1985) Peacock, J. A. & Heavens, A. F. 1985, Monthly Notices of the Royal Astronomical Society, 217, 805
  * Peebles (1971) Peebles, P. J. E. 1971, A&A, 11, 377
  * Pillepich et al. (2012) Pillepich, A., Porciani, C., & Reiprich, T. H. 2012, MNRAS, 422, 44
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * Pollack et al. (2012) Pollack, J. E., Smith, R. E., & Porciani, C. 2012, MNRAS, 420, 3469
  * Pollack et al. (2014) Pollack, J. E., Smith, R. E., & Porciani, C. 2014, MNRAS, 440, 555
  * Porciani et al. (2002) Porciani, C., Dekel, A., & Hoffman, Y. 2002, MNRAS, 332, 325
  * Prada et al. (2012) Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018
  * Press & Schechter (1974) Press, W. H. & Schechter, P. 1974, ApJ, 187, 425
  * Ragone-Figueroa et al. (2010) Ragone-Figueroa, C., Plionis, M., Merchán, M., Gottlöber, S., & Yepes, G. 2010, MNRAS, 407, 581
  * Ramakrishnan et al. (2019) Ramakrishnan, S., Paranjape, A., Hahn, O., & Sheth, R. K. 2019, MNRAS, 489, 2977
  * Rozo et al. (2014) Rozo, E., Rykoff, E. S., Bartlett, J. G., & Evrard, A. 2014, MNRAS, 438, 49
  * Salcedo et al. (2018) Salcedo, A. N., Maller, A. H., Berlind, A. A., et al. 2018, MNRAS, 475, 4411
  * Salcedo et al. (2020) Salcedo, A. N., Zu, Y., Zhang, Y., et al. 2020, arXiv e-prints, arXiv:2010.04176
  * Sato-Polito et al. (2019) Sato-Polito, G., Montero-Dorta, A. D., Abramo, L. R., Prada, F., & Klypin, A. 2019, MNRAS, 487, 1570
  * Schaap & van de Weygaert (2000) Schaap, W. E. & van de Weygaert, R. 2000, A&A, 363, L29
  * Shaw et al. (2006) Shaw, L. D., Weller, J., Ostriker, J. P., & Bode, P. 2006, ApJ, 646, 815
  * Sheth et al. (2013) Sheth, R. K., Chan, K. C., & Scoccimarro, R. 2013, Phys. Rev. D, 87, 083002
  * Sheth et al. (2001) Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1
  * Sheth & Tormen (1999a) Sheth, R. K. & Tormen, G. 1999a, MNRAS, 308, 119
  * Sheth & Tormen (1999b) Sheth, R. K. & Tormen, G. 1999b, MNRAS, 308, 119
  * Sheth & Tormen (2002) Sheth, R. K. & Tormen, G. 2002, MNRAS, 329, 61
  * Sheth & Tormen (2004) Sheth, R. K. & Tormen, G. 2004, MNRAS, 350, 1385
  * Sigad et al. (2000) Sigad, Y., Branchini, E., & Dekel, A. 2000, ApJ, 540, 62
  * Simon (2005) Simon, P. 2005, A&A, 430, 827
  * Singh et al. (2020) Singh, P., Saro, A., Costanzi, M., & Dolag, K. 2020, MNRAS, 494, 3728
  * Sinigaglia et al. (2021) Sinigaglia, F., Kitaura, F.-S., Balaguera-Antolínez, A., et al. 2021, ApJ, 921, 66
  * Skibba & Macciò (2011) Skibba, R. A. & Macciò, A. V. 2011, MNRAS, 416, 2388
  * Smith et al. (2007) Smith, R. E., Scoccimarro, R., & Sheth, R. K. 2007, Phys. Rev. D, 75, 063512
  * Somerville et al. (2001) Somerville, R. S., Lemson, G., Sigad, Y., et al. 2001, MNRAS, 320, 289
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757
  * Strauss et al. (1998) Strauss, M. A., Ostriker, J. P., & Cen, R. 1998, ApJ, 494, 20
  * Suto et al. (1992) Suto, Y., Cen, R., & Ostriker, J. P. 1992, ApJ, 395, 1
  * Tegmark & Peebles (1998) Tegmark, M. & Peebles, P. J. E. 1998, ApJ, 500, L79
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709
  * Tinker et al. (2010) Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878
  * Trowland et al. (2013) Trowland, H. E., Lewis, G. F., & Bland-Hawthorn, J. 2013, ApJ, 762, 72
  * Tucci et al. (2021) Tucci, B., Montero-Dorta, A. D., Abramo, L. R., Sato-Polito, G., & Artale, M. C. 2021, MNRAS, 500, 2777
  * Valageas (2011) Valageas, P. 2011, A&A, 525, A98
  * Vale & Ostriker (2004) Vale, A. & Ostriker, J. P. 2004, MNRAS, 353, 189
  * van de Weygaert et al. (2009) van de Weygaert, R., Aragon-Calvo, M. A., Jones, B. J. T., & Platen, E. 2009, ArXiv e-prints
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Naess, S., Genel, S., et al. 2018, ApJ, 867, 137
  * Wang et al. (2020) Wang, K., Mao, Y.-Y., Zentner, A. R., et al. 2020, MNRAS, 498, 4450
  * Wang & Zhao (2020) Wang, Y. & Zhao, G.-B. 2020, Research in Astronomy and Astrophysics, 20, 158
  * Wechsler et al. (2002) Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52
  * Wechsler & Tinker (2018) Wechsler, R. H. & Tinker, J. L. 2018, ARA&A, 56, 435
  * Wechsler et al. (2006) Wechsler, R. H., Zentner, A. R., Bullock, J. S., Kravtsov, A. V., & Allgood, B. 2006, ApJ, 652, 71
  * Werner & Porciani (2020) Werner, K. F. & Porciani, C. 2020, MNRAS, 492, 1614
  * White (2001) White, M. 2001, A&A, 367, 27
  * White (1984) White, S. D. M. 1984, ApJ, 286, 38
  * Wong & Taylor (2012) Wong, A. W. C. & Taylor, J. E. 2012, ApJ, 757, 102
  * Yan et al. (2013) Yan, H., Fan, Z., & White, S. D. M. 2013, Monthly Notices of the Royal Astronomical Society, 430, 3432
  * Yang et al. (2017) Yang, X., Zhang, Y., Lu, T., et al. 2017, ApJ, 848, 60
  * Zehavi et al. (2019) Zehavi, I., Kerby, S. E., Contreras, S., et al. 2019, ApJ, 887, 17
  * Zel’dovich (1970) Zel’dovich, Y. B. 1970, A&A, 5, 84
  * Zentner (2007) Zentner, A. R. 2007, International Journal of Modern Physics D, 16, 763
  * Zentner et al. (2014) Zentner, A. R., Hearin, A. P., & van den Bosch, F. C. 2014, MNRAS, 443, 3044
  * Zhai et al. (2019) Zhai, Z., Benson, A., Wang, Y., Yepes, G., & Chuang, C.-H. 2019, MNRAS, 490, 3667
  * Zhao et al. (2021) Zhao, C., Chuang, C.-H., Bautista, J., et al. 2021, MNRAS, 503, 1149
  * Zhao et al. (2015) Zhao, C., Kitaura, F.-S., Chuang, C.-H., et al. 2015, MNRAS, 451, 4266
  * Zhou et al. (2023) Zhou, S., Zhang, P., & Chen, Z. 2023, MNRAS, 523, 5789

##  Appendix A Principal component analysis of halo properties

Halo mass is not the main halo property, as it can swap or share such stage
with other quantities such as $V_{\rm max}$ (see e.g. Skibba & Macciò 2011;
Wong & Taylor 2012), velocity dispersion or the peak velocity (as these are
also probes of the depth of the halos gravitational potential) (see e.g.
Jeeson-Daniel et al. 2011; Lee et al. 2017).

Figure 30: Cumulative variance $\mathcal{C}_{i}$ (Eq.(17) for each principal
component $PC_{i}$ at different redshifts and cosmic-web environments. The
horizontal lines mark the $30$, $60$ and $100\%$ of the full variance. Almost
$90\%$ of the information is retained by the $\sim 50\%$ of the highest
principal components. The first component encodes $\sim 30\%$ of the full
information.

The need for a _primary_ property of dark matter halos is key in the context
of galaxy clustering, where typical approaches to the population of halos such
as the halo occupation distribution (HOD) framework (e.g. Cooray 2002; Cooray
& Sheth 2002; Berlind & Weinberg 2002; Kravtsov et al. 2004), or sub-halo
abundance matching (see e.g. Vale & Ostriker 2004; Kravtsov et al. 2004;
Conroy et al. 2006; Zentner et al. 2014; Favole et al. 2016; Hearin et al.
2016; Wechsler & Tinker 2018), which are historically based on halo masses,
with extensions to other properties (see e.g. Zehavi et al. 2019; Hadzhiyska
et al. 2020). This question is also present in the context of the generation
of mock catalogs using calibrated methods (Balaguera-Antolínez et al. 2023).

Figure 31: Importance $\mathcal{I}_{i}(z)$ (Eq.(18) for each halo property as
a function of redshift and cosmic-web environment. We have divided the set of
halo properties in four panels in order to avoid clutter. The horizontal lines
show the $5$ $10$ and $15\%$ importance levels. The level of importance is
normalized to be unity when all the information from the shown properties is
added.

To determine the main halo property in the UNITsim, we perform a principal
component analysis (PCA) (see e.g. Jolliffe & Cadima 2016; Skibba & Macciò
2011; Wong & Taylor 2012; Zhou et al. 2023), transforming each halo property
$\eta_{\alpha}$ to its standardized form,
$\eta_{\alpha}\equiv(\theta_{\alpha}-\langle\theta\rangle)/\sigma_{\theta}$,
where $\langle\theta\rangle$ is the sample mean of the property and
$\sigma_{\theta}=\langle(\theta_{\alpha}-\langle\theta\rangle)^{2}\rangle^{1/2}$
its corresponding sample variance. We then compute the eigenvalues
$\tilde{\lambda_{i}}$ (sorted in descending order) and their associated
eigenvectors $\hat{\lambda}^{q}$ with elements $\hat{\lambda}^{(q)}_{i}$
($i,q=1,2\cdots,N_{p}$) of the covariance matrix whose elements are
$\mathcal{C}_{ij}\equiv\langle\eta_{i}\eta_{j}\rangle$. The set of
eigenvectors defines a new orthogonal basis in which $\hat{\lambda}^{(1)}$
represents the axis on which the data set displays the largest variation. A
new set of halo parameters (the principal components) can be then generated as
as
$\tilde{\eta}^{\alpha}_{q}=\sum_{i}(\hat{\lambda}^{(q)}_{i})\eta^{\alpha}_{i}$
($p,q=1.,N_{PCA}$) where $(\hat{\lambda}_{i})_{q}$ are the components of the
eigenvector associated to the $q-$th eigenvalue. The cumulative variance
retained in each principal component is

$\mathcal{C}_{i}(z)=\frac{1}{\mathcal{V}(z)}\sum_{j<i}\tilde{\lambda}_{j}(z).$
(17)

and shown in Fig.30 for each principal component, evaluated at number of
redshifts and cosmic-web environments. Trends as a function of the
cosmological redshift are consistent within the different cosmic-web
environment, indicating that more principal components are requested as long
as $z\to 0$ in order to retain the same amount of information on the halo
properties. For example, the change in the number of PC’s requested to have
$90\%$ of full variability goes from $\sim 6$ at $z\sim 5$ to $\sim 8$ at
$z=0$. Averaging over the values shown, the first PC contains $\sim 30\%$ of
the total amount of information, the second holding $\sim 15\%$. We have
verified that this conclusion are not modified when the cumulative variance is
computed in bins of halo mass.

We define the _importance_ of each halo property $\eta_{i}$ (as a function of
redshift) as:

$\mathcal{I}_{i}(z)\equiv\frac{1}{\mathcal{V}(z)}\sum_{p=1}^{N_{PCA}}\tilde{\lambda}_{p}(z)|\rho_{ip}(z)|,$
(18)

normalized such that $\sum_{i}\mathcal{I}_{i=1}^{N_{p}}=1$ for all redshifts
and cosmic web types. In Fig.31 we show the _importance_ of all halo
properties as a function of cosmological redshift (for this calculation we
have used all principal components) and in different cosmic-web environments.
For each redshift and cosmic-web type, Eq. (18) is normalized adding the
values from the properties shown in the plot. The main results can be
summarized as follows:

  * •

The set of intrinsic properties $\\{M_{\rm vir},V_{\rm max},\sigma_{v},\nu\\}$
display the highest level of importance within the set of halo properties,
very similar evolution across redshift ($10-15\%$ on the explored redshift
range). This behavior mildly changes in different cosmic-web environments.
This set of properties amasses $\sim 40\%$ of the full information encoded in
the set of halo properties explored.

  * •

The scale radius $R_{s}$ displays an increasing trend of importance from high
to low redshift, with $\mathcal{I}\sim 10\%$ at $z\sim 0$: The halo
concentration with $\sim 5\%$ of average importance on average, overcoming
that limit only in sheets at all redshifts.

  * •

The environmental properties
$\\{\mathcal{M}_{5},\mathcal{D}_{5},\Delta_{5},\mathcal{T}_{A}\\}$ display an
average $\mathcal{I}\sim 8\%$. Their higher value is reached in voids,
(specially $\mathcal{D}_{5}$). These properties amass $\sim 25\%$ of the total
level of importance.

  * •

The set conformed by properties such as spin, the virial and halo geometry,
$\\{\lambda_{B},T/|U|,\mathcal{T}_{h}\\}$ show a $3\%$ importance over the
redshift range explored and for all cosmic web types.

In conclusion, most of the information is retained by probes of the potential
well. In general, this claim is valid for all cosmic-web environments,
reinforcing the define the set $\\{M_{\rm vir},\nu,V_{\rm max},\sigma_{v}\\}$
as _primary properties_.

##  Appendix B Secondary bias as a function of peak-height

In Fig.32 we present the measurements of relative bias (secondary bias in the
fourth quartile divided by the same signal from the first quartile) using the
for the halo ellipticity, relative Mach number, local dark matter density and
tidal anisotropy as secondary properties.

Figure 32: Ratio between halo bias as a function of the peak-height obtained
in two different quartiles of halo ellipticity $\mathcal{E}_{h}$, relative
Mach number $\mathcal{M}_{5}$, dark matter density $\Delta_{\rm dm}$ and tidal
anisotropy $\mathcal{T}_{A}$, at different redshifts of the UNITsim.