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

,


Introduction
It is well known that the clustering of dark matter halos strongly depends on their internal properties.In this context, the halo mass is usually considered to be responsible for one of the primary dependencies, as a direct manifestation of the more fundamental dependence of the halo bias on the peak height of density fluctuations, ν (e.g., Press & Schechter 1974;Kaiser 1984; Bardeen et al. 1986;Sheth & Tormen 1999a;Sheth et al. 2001;Sheth & Tormen 2002).More massive halos are more strongly clustered than less massive halos, as is indeed described by the Λ-cold dark matter (Λ-CDM) structure formation formalism.However, a number of additional secondary dependencies at a 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. ⋆ balaguera@iac.es2017; 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.These types 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 the halo age and proxies of the 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 the halo bias at a fixed halo mass as a secondary halo bias.These also include environmental properties and parameters characterizing the tidal field Article number, page 1 of 37 arXiv:2311.12991v2[astro-ph.CO] 10 Feb 2024 around halos, which are known to display such signals1 (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 toward 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 subpopulation of halos whose accretion was halted early on.Hahn et al. (2009) subsequently claimed that tidal effects produced by a neighboring 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).Based on zoomin simulations of a small number of halos, Borzyszkowski et al. (2017) concluded that low-mass assembly bias is due with the existence of same-mass "stalled" and "accreting" halos, typically associated with filaments and nodes, respectively.Musso et al. (2018) extended this analysis to take the whole geometry of tides into account.It has also been shown that the secondary bias signals (i.e., not only assembly bias) correlate with the anisotropy of the tidal tensor at a 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).For this paper, however, we implemented an object-by-object estimator of the halo bias (Paranjape et al. 2018), which allows for a robust statistical characterization of the halo bias and its dependencies on a number halo properties.This bias estimation was employed to explore, for the first time, the signal of the 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 ∼ 5. We dissected the secondary bias by analyzing a wide range of both intrinsic halo and environmental properties.The assessment of the secondary bias performed in this work includes novel halo properties such as the relative Mach number and the neighbor statistics, as well as the evaluation of properties that are inherent to the dark matter field, extending recent analysis on the same subject (Wang et al. 2021).Our analysis is not only relevant in the context of secondary bias.We also performed 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 the 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 and describe the different halo properties employed in the analysis.In particular, we char-acterize the correlation among the different halo properties and the impact of the environment on their scaling relations.In §4 we show the measurements of the secondary bias based on an object-by-object estimator and show the scaling relations between the large-scale bias and the halo properties.Section 5 presents a methodology for the reconstruction of a secondary bias and potential applications for the generation of halo mock catalogs.Finally, in §6 we discuss our results, present the main conclusions of this research, and comment on forthcoming research in line with the present work.

The UNIT simulation
The UNIT suite of simulations (UNITSim hereafter)2 (Chuang et al. 2019) is a set of N-body simulations of dark matter particles in cosmological volumes evolved from z = 99 until z = 0 using the code Gadget (Springel et al. 2001).Initial conditions are generated using the FastPM algorithm (Feng et al. 2016) imposing a linear matter power spectrum in a Λ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(Gpc h −1 ) 3 , run with 4096 3 dark matter particles with mass resolution of 1.2 × 10 9 M ⊙ h −1 .The dark matter field is represented by a mesh of 2048 3 cells, which we have reduced to 512 3 by averaging in configuration space 3 , yielding a spatial resolution of ∼ 1.95Mpch −1 .We shall refer to this as our fiducial resolution.
Importantly, the UNITSim is a paired-fixed amplitude simulation (Angulo & Pontzen 2016).On the one hand, "fixed-amplitude" means that, while the phases ϕ k of the complex overdensity field in Fourier space δ(k) = A k e iϕ k are randomly distributed in the range [0, 2π), their amplitudes A k follow a distribution of the form δ D (A k − V P 0 (k)/(2π) 3 ), instead of a Rayleigh distribution with mode P 0 (k), as it is the case of 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 π/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 toward 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 been also demonstrated that these types of simulations do not introduce a 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 the distribution of 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 almost indistinguishable 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).

Halo and environmental properties: Correlation and scaling relations
In this section we present the halo properties employed in the analysis, and we discuss the intricate relations between them as a function of redshift.The set of halo properties used in this work is divided in three sectors, namely i) intrinsic properties (i.e., those directly provided by the halofinder), ii) nonlocal properties (i.e., those computed from the phase-space properties of neighbor halos) and iii) environmental properties (i.e., computed from the underlying dark matter density field.

Intrinsic halo properties
In the UNITSim, halo catalogs and their properties have been obtained using the ROCKSTAR (Robust Overdensity Calculation using K-Space Topologically Adaptive Refine) halo finder algorithm (Behroozi et al. 2013).The robustness of this halo finder has been studied in a number of comparison projects (see e.g., Knebe et al. 2011;Evrard et al. 2014;Rodríguez-Puebla et al. 2016;Mansfield & Kravtsov 2020;Mansfield & Avestruz 2021).In particular, based on scalefree cosmologies, Leroy et al. (2021); Maleubre et al. (2024) assessed the level of convergence to the physical limit of the statistical properties of dark matter halos defined with the ROCKSTAR algorithm, showing its better performance with respect to other halo finders, a claim that is arguably extrapolated to non scale-free simulations as the UNITSim.
In general, we expect that our results can be extrapolated to those obtained with other halo-finders, although a thorough comparison in this direction is beyond the scope of this work.
Throughout this paper, we shall refer to "scaling relations" (or link) between halo properties.Such term must be understood as the probability that a halo has a property a conditional to a property b in the range b, b + db, i.e, P(a|b)db.

Halo mass
The halo finder provides a number of definitions of halo mass, based on a spherical overdensity ∆ (Behroozi et al. 2013).In this work we use the virial mass M vir = (4/3)π ρ(z)∆(z)R 3 vir , obtained by the halo finder using spherical overdensities of ∆(z) = (18π 2 + 82x − 39x 2 )/Ω mat (z) where x ≡ Ω mat (z) − 1 (Bryan & Norman 1998).Here R vir denotes the radius at which the enclosed mass equals the virial mass (see e.g., Klypin et al. 2016), a quantity similarly provided by the halo finder.We select parent halos with M vir > 2×10 11 M ⊙ h −1 (which corresponds to ∼ 300 dark matter particles).The selection of M vir as a mass-proxy is consistent with other halo properties such as spin, concentration and shape parameters, which are computed based on R vir and hence M vir .
The halo mass distribution at z = 0 is shown in the panel (a) of Fig. 3 as a function of the cosmic-web type (to be defined in §3.3.1).We have verified that the halo abundance as a function of virial mass is in good agreement with fitting functions found in literature (see e.g., Tinker et al. 2008).

Maximum circular velocity
This quantity is computed as is the halo circular velocity at a radius r and M(r) is the mass enclosed within that radius (and G N is the gravitational constant).A simple version of the virial theorem (see e.g., Collins 1978) predicts V max ∝ M 1/3 vir , a dependency that we have confirmed in the UNITSim at different redshifts (see also Shaw et al. 2006;Rodríguez-Puebla et al. 2016;Cataldi et al. 2023).Panel (b) of Fig. 3 contains the distribution of V max at z = 0, while its scaling relation with the virial mass is shown in panel (a) Fig. 4. The maximum circular velocity can be regarded, along with the halo mass, as direct probes of the potential well of dark matter haloes (see appendix §A).As such, this quantity is sometimes preferred over the mass in order to model the statistics of galaxies inside halos (see e.g., Zehavi et al. 2019;Lehmann et al. 2017;Mezini et al. 2022).
The halo finder also provides the velocity dispersion of dark matter halos, computed from velocity statistics of dark matter in halos.As for V max , this quantity scales approximately as ∝ M 1/3 vir (see e.g., Ragone-Figueroa et al. 2010).Even though we do not show results in terms of this quantity, it is important to indicate that the behavior of the halo velocity dispersion as a function of halo mass is key for the assignment of phase-space coordinates for galaxies in dark matter halos in the framework of the Halo model (see e.g Cooray 2006;Balaguera-Antolínez et al. 2023).

Scale radius and halo concentration
High-resolution N-body simulations have demonstrated that the density profile of dark matter halos can be described by the Navarro-Frenk-White profile (Navarro et al.

1996
).As such, ROCKSTAR uses this profile (see e.g., Knebe et al. 2011;Mansfield & Kravtsov 2020) to provide values of scale radius R s and the halo concentration C vir .For the scale radius, we have used the estimates based on fits to the NFW profile, although the halo finder also provides estimates based on the method by (Klypin et al. 2011) (which also relies on the NFW profile).Besides of shaping the density profile (it marks the scale at which the logarithmic slope of the NFW profile equals 2), the scale radius R s has been shown to correlate with halo anisotropy parameters (see e.g., Binney & Tremaine 2008;Barnes et al. 2005) and hence can encode information of the dynamical state of dark matter halos.Panel (c) of Fig. 3 shows the distribution of this property at z = 0 in the UNITSim.Panel (b) of Fig. 4 shows the scaling relation between the scale radius (in kpc/h) and the halo mass for different redshifts.
The halo concentration, computed as C vir ≡ R vir /R s , is a key property, not only for characterizing the halo density profile, but also because it has been shown to be a proxy for halo formation time at fixed halo mass (Wechsler et al. 2002;Ludlow et al. 2016;Wang et al. 2020).Panel (d) of Fig. 3 shows the distribution (at z = 0) of halo concentration in the UNITSim, while panel (c) of Fig. 4 displays the mean scaling relation in terms of halo mass (see e.g., Macciò et al. 2007;Ludlow et al. 2016;Child et al. 2018, for a thorough analyusis of this scaling relation).This massconcentration scaling relation shows interesting features, as it displays a positive slope at high redshift (i.e, more massive halos are more concentrated), decreasing toward low redshift, to become fully inverted for redshifts below z ∼ 1.This behavior can be best understood when the halo concentration is expressed as a function of the peak height ν (to be defined in §3.1.7)in conjunction with cosmological parameters describing the statistical properties of primordial density fluctuations; we have verified the C vir − ν relation displays an approximately universal shape as a function of redshift (see e.g., Prada et al. 2012;Diemer & Joyce 2019).

Dimensionless spin parameter
ROCKSTAR computes dimensionless spin parameters using both the Peebles (1971) and the (Bullock et al. 2001) prescription.The latter defined as λ B = J/( √ 2M vir V c (R vir )R vir ) where J is the total halo angular momentum.We use this version of halo spin, and we have verified that consistent results are obtained for the test carried within this paper when using the Peebles' definition.Panel (e) in Fig. 3 presents the distribution of halo spin (at z = 0) for different cosmic web types, from which we can observe deviations (in the form of kurtosis and skewness) from a log-normal distribution (Bullock et al. 2001).Panel (d) of Fig. 4 displays its behavior with respect to halo mass at different redshifts.The behavior observed for the spin at different redshifts is in principle compatible with previous analysis (see e.g., Macciò et al. 2007;Knebe & Power 2008;Ahn et al. 2015), displaying a weak correlation with the halo mass across redshifts, passing from a correlation coefficient of ∼ −0.02 at z ∼ 5 to ∼ 0.01 at z = 0. Nevertheless, we see that the behavior of the scaling relation, despite the small correlation, has a remarkable deviation from a constant behavior, contrary to what has been reported in the literature.This is a matter of future research.

The virial parameter
We refer as the "virial parameter" to the ratio between kinetic to gravitational potential energy, T/|U|4 .This quantity can be used as a proxy for the degree of relaxation of a halo (Knebe & Power 2008), as the virial theorem (see e.g., Collins 1978) predicts T/|U| ∼ 1/2 for relaxed isolated halos.Nevertheless, 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 T/|U| ∼ 1/2 do not strictly imply a lack of relaxation (Henriksen & Widrow 1999;Shaw et al. 2006;Bett et al. 2007;Klypin et al. 2016).Halos in the UNITSim have T/|U| > 1/2, as can be checked from panel (f) of Fig. 3. Panel (e) of Fig. 4 presents its dependency with respect to the virial mass for different redshifts.

Halo triaxiality and ellipticity
Ellipsoidal shape parameters are defined in terms of the ratio between the halo semi-axes (a, b, c) as the triaxiality parameter T h and the halo ellipticity: where q = b/a and s = c/a (a > b > c).The set (a, b, c) is computed by the halo finder as the sorted eigenvalues of the mass distribution tensor using particles inside the virial radius (see e.g., Zemp et al. 2011).The triaxiality parameter T h defines ellipsoidal geometries such as the oblate (T h → 0) and prolate (T h → 1) spheroids.On the other hand, halo ellipticity E h has as special limits the spherical (E h → 0) and nonspherical (or aspherical) configurations (see e.g., Cole & Lacey 1996;Binney & Tremaine 2008).An example of the distribution of this geometry-proxies at z = 0 is displayed in panels (g,h) of Fig. 3.The dependence with the virial mass at different redshifts is shown in panels (f,g) of Fig. 4, from where we see how, at a fixed mass, low redshift halos are inclined to be more rounded with oblate-shapes, in agreement with previous results (see e.g., Shaw et al. 2006).The behavior of these scaling relations are linked to the merging history of halos (see e.g., Allgood et al. 2006).The exact details on the connection between mergers and its impact on halo properties are to be addressed in forthcoming research.

Peak height
We assign to each halo a value of peak height ν(M, z) = δ sc /σ(M, z), where δ sc = (3/20)(12π) 2/3 is the critical collapse threshold linearly extrapolated to redshift z = 0 and is the variance in the matter distribution on scales R ∝ (M vir / ρ) 1/3 (see e.g., Peebles 1980).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 the Fourier transform of a top-hat filter with scale R and ρ(z) is the mean matter density at redshift z.In this work, we use log ν as a primary halo property, that can be exchanged with the halo mass according its definition.The mass scale set by ν(M ⋆ , z) = 1 defines the mass of collapsing halos at redshift z.At z = 0, M ⋆ ∼ 3.1 × 10 12 M ⊙ h −1 .The peak height ν is the natural variable in which the peak-formalism can predict halo abundances and bias (Kaiser 1984;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(Tinker et al. , 2010;;More et al. 2011).Also, trends in scaling relations (such as the mass-concentration relation, §3.1.3)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 ν(M, z) (so it vanishes at fixed ν(M, z) (see e.g., Wechsler et al. 2006;Gao & White 2007)).Peak height ν(M, z) is, therefore, a cornerstone for theoretical models of halo bias, including its secondary dependencies.

Nonlocal 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 "nonlocal", as they are computed based on the phase-space properties of tracers enclosed on spheres of radius R = 5Mpc h −1 around each halo.As such, this set can provide information of the small-scale clustering of dark matter haloes.

Relative Mach number
The Mach number 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 "relative Mach number" by computing the relative velocity v i j = v i −v j of a dark matter halo (labeled i) with peculiar velocity vector v i , and its neighboring halos with velocities v j inside spheres of radius R, namely: where ⟨⟩ j,R denotes averages inside spheres of radius R and σ i is the dispersion in the variable |v i j | inside those spheres.As in the original definition of the Mach number, high (low) values of M (i) R can be associated with warm 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 panel (i) of Fig. 3 we present the distribution of local halo Mach number at z = 0 for different cosmic-web types.In addition, panel (h) of 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 ∼ 10 13 M ⊙ /h varies little.The information encoded in Fig. 2 reveals that the halo population, up to the mentioned mass scale, is mainly living in filaments; above that mass, 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 panel (h) of Fig. 8, where the scaling relation is split in different cosmic-web environments (to be defined in §3.3.1).

Neighbor statistics
Along with the relative Mach number previously defined, we introduce the "neighbor statistics" D R , as a probe for the statistics of pair separations around each tracer.This estimator is computed as where ⟨λ i j ⟩ j,R is the mean separation of halos at a position r j and within spheres of radius R from the main tracer, while σ λ R is its population variance.The neighbor statistics aims at compressing the two first moments of the distribution of pair separations, and as such, is expected to contain information on the small-scale clustering of halos.Panel (k) of Fig. 3 shows the distribution of D R across the halo population at z = 0 in different cosmic-web environments, while panel (j) in Fig. 4 displays the mean scaling relation between halo mass and D at different redshifts.Similarly, panel (j) in Fig. 8 shows (at z = 0) the scaling relation as a function of cosmic-web environments (to be defined in §3.3.1),evidencing a quantitative difference between this statistics in low (voids) and high (density regions), specially at low halo masses.Similar quantities have been explored in the context of secondary bias, such as a scale dependent correlation (Kerscher 2018) or "neighbor distance" (Salcedo et al. 2022).The impact of massive neighbors in the signal of large-scale bias will be addressed in a forthcoming publication.

Halo overdensity
We compute a local halo overdensity ∆ R by collecting the number of tracers N R used to measure M R (or D R ) and defining ∆ R = log(N R / N), where N is the expected number of randomly distributed tracers in spheres of volume ∝ R 3 .The distribution of this variable within the halo population at z = 0 is shown in panel (j) of Fig. 3, while its dependency with the halo mas is shown in panel(i) of Fig. 4, reflecting the expectation that more massive tracers tend to be surrounded by more neighbors compared to low-mass halos.

Environmental properties
By "environmental properties", we refer to quantities computed from the underlying dark matter field from which the halo catalog has been built.These properties can be classified as local, such as the dark matter density in a cell   Spearman's rank correlation coefficient ρ s between the halo effective bias and a number of halo properties for two different redshifts.
link between halo mass and the local overdensity ∆ dm at different redshifts.

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 within the cosmic structure (see e.g., Schaap & van de Weygaert 2000;González & Padilla 2010;Aragón-Calvo et al. 2010;Sousbie et al. 2011;Libeskind et al. 2018, and references therein).Here we use the approach by Hahn et al. (2007) and compute the Hes- sian of the gravitational potential (or tidal field) T i j = ∂ i j Φ (where Φ is the comoving gravitational potential obtained from Poisson's equation ∇ 2 Φ = δ dm )6 and its eigenvalues λ 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) , where λ th is an arbitrary threshold.We use λ th = 0, which has the practical advantage of allowing 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 (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;Kitaura et al. 2022;Bonnaire et al. 2022;Wang et al. 2023) , 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 (with volume defined by the fiducial resolution, see §2) 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 toward 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 (δ dm < 0) and the formation of halos in regions with δ dm > 0.

Tidal anisotropy
Based on the set of eigenvalues {λ i } defined in §3.3.1, the degree of anisotropy of the dark matter density field can be quantified by computing geometrical quantities such as the ellipticity e = (λ 1 − λ 3 )/2δ dm or the prolateness p = (λ 1 + λ 3 − 2λ 2 )/2δ 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 Notice that this definition differs from others often found in the literature e.g., Hahn et 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 + δ dm (instead of 1 + δ dm ), in order to avoid divergences.According to the threshold chosen for the eigenvalues λ i , filaments and sheets are expected to display higher values of tidal anisotropy compared to knots and voids.Hence, very high-mass halos, mainly found in high density regions (knots), are expected to reside in regions with low anisotropy, as can be confirmed by panel (l) of Fig. 3 and panel (l) of Fig. 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 depth of the halo potential well (e.g virial mass, maximum circular velocity), as well as the correlation among environmental properties (e.g., Mach number, neighbor statistics, tidal anisotropy).
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 evident for properties such as halo concentration and spin, as can be also seen from panels (c) and (d) 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, maximum circular velocity) contain most of the variability (understood here as information), nonlocal and environmental properties can contain similar levels of variability.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 a low variability can be associated with a large scatter in their link with, for example, halo mass, as can be inferred from Fig. 6.

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 the panels in Fig. 3 display the distribution of halo properties (intrinsic and nonlocal) evaluated in different cosmic-web environments, with the noticeable feature that these distributions are much more sensitive to under-dense regions (voids).To complement these results, Fig. 7 shows the Spearman's rank correlation coefficient ρ 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 are a signature that the links among halo properties can be sensitive to the cosmological environment they live in.In Fig. 8 we present the mean scaling relations between halo mass and other properties at z = 0, evaluated in different cosmic-web types.Some conclusions based on this figure are (we have verified that these conclusions are representative of the behavior at different cosmological redshifts): -Properties probing the depth of the halo potential well (e.g., halo mass, maximum circular velocity) are linked through a rather stable scaling relation against changes in the cosmic web type, as shown in panel (a).We have checked that this is also the case for the velocity dispersion.
-Statistically significant differences among scaling relations of halo secondary properties are mainly found in voids, as can be seen for the scale-radius R s (panel b), the halo concentration C vir (panel c) and the spin λ B (panel d).In particular, halo spin scaling relation not only displays a lower amplitude in voids (i.e, 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) compared to other cosmic-web classifications.-The ratio T/|U| computed for tracers in voids is lower than those values seen in high density region, as shown in panel (e) of Fig. 8.If we assume that relaxation is represented as T → 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.-The dependency of ellipticity and triaxiality with halo mass are not (statistically speaking) sensitive to the cosmic-web type, except, again, in voids: halos with masses up to ∼ 10 12 M ⊙ h −1 in such low-density environments tend to display higher ellipticity and triaxiality (i.e, halos tend to be more prolate) than halos of the same mass in other environments, as shown in panels (f, g) of Fig. 8. -The mean scaling relations between nonlocal properties and halo mass also display (statistically) large differences among cosmic-web types.In particular, halos in voids found themselves in a rather cool environment compared with other cosmic-web types, as shown in (panel h).Also, voids are sparsely populated, as shown by the local halo overdensity (panel i) and the neighbor statistics (panel j) and in regions with a mean halo density below the mean (local overdensity).
Article number, page 11 of 37 -The comparison between the z = 0 curve in panel (k) of Fig. 3 and panel (k) of Fig. 8 shows how the population of halos, expressed through their halo mass, is sensitive to the cosmic-web environment.This is a key aspect exploited in bias mapping techniques aimed at generating mock catalogs (see e.g., Balaguera-Antolínez et al. 2023).
-Environmental properties such as local dark matter density (panel k) and tidal anisotropy (panel l) show large differences, which is a trivial consequence of their dependence on the cosmic-web classification.In particular, sheets and filaments are characterized by a large anisotropy.
These conclusions are in agreement with previous analysis (see e.g., Shaw et al. 2006;Hellwing et al. 2021), and highlight the importance of the environment in shaping the link between different halo properties.The behavior of the halo scaling relations might not only change due to the environment, but also through the underlying cosmological model (see e.g., Ho et al. 2018).Also, halo properties and scaling relations have been shown not only to correlate with the cosmic environments, but also with properties of its components such as the thickness of filaments (Liao & Gao 2019;Galárraga-Espinosa et al. 2020;Ganeshaiah Veena et al. 2021), or the mass of collapsing regions defined as knots (Zhao et al. 2015;Balaguera-Antolínez et al. 2023;Davies et al. 2021) As a general remark, we can argue that the differences seen in the scaling relations, specially those linked to spin and geometry, can be understood as a signature of the coupling between halo dynamical properties and the 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).Nevertheless, there could be other evolutionary effects (aside the growth in the strength of such coupling) that influence the shape of the scaling relations.We remind the reader that our results are based on a very particular cosmic-web classification, and in general, the properties of halos in particular environments such as voids can vary when using different techniques to define them (see e.g., Colberg et al. 2008;Cautun et al. 2018).
A deep understanding the halo scaling relations demands a complex analysis and is beyond the goal of this paper, as it deserves a project on its own.We shall pursue this in future work.In the coming sections we focus on the main subject of this paper, namely, the dependencies of large-scale halo bias on the halo properties and the signal of secondary bias.

The large-scale halo bias
It is very well established that the bias of dark matter tracers needs to be modeled beyond the linear scaleindependent scheme (see e.g., Desjacques et al. 2018, and references therein), in which the halo bias b simply relates dark matter δ dm with halo δ h overdensities (determined on a mesh characterized by some scale R) through an expression of the form δ h = bδ dm .Scale dependencies induced by the process of halo formation, halo merging history, nonlinear 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 generalizes the concept of halo bias to a nonlocal and stochastic quantity, characterized by a conditional probability distribution P(δ h |δ dm ) R (see e.g., Fry & Gaztanaga 1993;Tegmark & Peebles 1998;Dekel & Lahav 1999;Blanton 2000;Simon 2005).
The effective large-scale halo bias measured from twopoint statistics, even if assessed on scales where linear perturbation theory applies, is not simply the factor b defining a link of the form δ h = bδ dm .Instead, it is the combination (or renormalization) of linear-order coefficients involved in the perturbative expansion of the overdensity field, aimed to describe the tracer density field, (see e.g., Heavens et al. 1998;McDonald 2006;McDonald & Roy 2009;Desjacques et al. 2018;Werner & Porciani 2020).Measurements of higher order effective bias (e.g., quadratic bias) need to be accomplished with other probes such as cross-correlations (see e.g., Smith et al. 2009), higher order statistics (see e.g., Pollack et al. 2014), direct approaches such as the so-called "separate universe technique" (see e.g., Wagner et al. 2015;Lazeyras et al. 2016;Paranjape & Padmanabhan 2017) or effective filed theories of large scale structure (see e.g., Senatore 2015; Lazeyras et al. 2017).

Estimators of effective bias
The large-scale effective halo bias can be measured via a number of methods, such as number-count-in-cells (e.g., Peebles 1980;Bernardeau 1994;Fry et al. 2011), or twoand 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 (see e.g., Colín et al. 1999).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): or based on the halo-dark matter cross-power spectra P hm (k): 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 = Article number, page 12 of 37 Author: Balaguera-Antolínez, Montero-Dorta & Favole 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 nonlinear 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. ( 8) 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 where δ dm (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 ≤ 0.08 hMpc −1 .This limit is rather conservative for high redshift, where nonlinear evolution of the halo bias takes place at much smaller scales.At low-redshift (z ∼ 0), the ratio between halo and dark matter power spectrum is already scale-dependent at k ∼ 0.1 h Mpc −1 .The large-scale effective bias of a subsample of N h (∆θ) tracers in a bin of property ∆θ can be obtained as a sample mean: where N h (∆θ) denotes the number of halos in the ∆θ bin.
The effective bias estimator of Eq. ( 10) has a number of Article number, page 13 of 37 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 twopoint 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. ( 10) 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).Mean halo bias as a function of the 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 9 shows a slice of thickness ∼ 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 ≲ 60 Mpc.This can be considered as a consequence of setting a limit k max to compute the effective bias: in practice, such limit acts as a low-pass filter that reduces the spatial resolution of the dark matter and halo fields from ∼ 1.95 Mpc h −1 (see §2) to π/k max ∼ 40 Mpc h −1 .In appendix §C we show how the bias-weighted distribution of separation of tracers reflects this fact in the form of a ringing effect.
Figure 10 shows the scaling relation between the effective halo bias computed with Eq.( 10) and both peak height, ν (left panel), and halo mass (right panel).In both pan-Article number, page 14 of 37 Author: Balaguera-Antolínez, Montero-Dorta & Favole We have compared the estimates of effective halo bias from the estimators of Eq. ( 7) (8) and ( 10).This is presented in Fig. 11 7 , showing measurements obtained using Fourier modes to a maximum k max = min(0.08,k sn ) hMpc −1 , where k sn is the wavenumber at which the signal of power spectrum equals the Poisson shot-noise.The results based on Eq. ( 7) are in agreement with those based on Eq. ( 10), while the measurement of b hm based on Eq. ( 8) differs by a factor ∼ 1σ. 8 The overall agreement among these estimates confirms the robustness of Eq. ( 10) in terms of measuring effective halo bias.From now on, we shall refer to the effective halo bias with the symbol b h .
7 All 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 shotnoise subtraction.When averaging in Fourier bins, use weighted averages with number of modes in each bin as the weights. 8Both 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.Spearman's rank correlation coefficient between the halo effective bias and a number of halo properties as a function of redshift.

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 ν) with deviations on the order of 1σ compared to T10.The universality of the halo bias can be extended to the full bias distribution.To verify that, in Fig. 13 we present the estimates of the higher moments (variance, skewness and kurtosis) of the P(b h |ν) 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 (ρ s ≤ 0.03 at all redshifts explored, see Fig. 6) leads to the conclusion that effective halo bias follows a normal distribution, i.e, P(b h |ν, z) ∼ N(⟨b h |ν⟩, σ ν (z)), where the standard deviation σ ν (z) is independent of ν9 .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 formalism of structureformation predicts (Bardeen et al. 1986;Sheth & Tormen 1999a, 2002).
The scaling relations between halo bias and any other halo property θ (internal or external) can be expressed based on their link to peak height (or halo mass) 10 using such that the mean effective bias as a function of a property θ can be written as where P(ν) = db h P(b h |ν) is the distribution function of dark-matter halos as a function of their peak height.Equations ( 12) and ( 13) (the latter is the continuous representation of Eq.( 10)) represent the link between large scale 10 As the scaling relation between peak height and halo mass is P(ν|M vir )dM vir = δ D (ν − ν 0 (M vir ))dM vir where ν 0 (M vir ) = δ sc /σ(M vir ) (see Eq. ( 3)), the statistical properties of the halo bias as a function of halo mass are directly inherited from those of peak height.
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 Eq.( 12) is a symbolic separation, as the scaling relations of halos contain cosmological information (Wechsler et al. 2002;Ludlow et al. 2016;Diemer & Joyce 2019) It is also key to notice that previous attempts to assign an individual bias were based on Eq.( 13) (see e.g., Schuecker et al. 2001;Balaguera-Antolínez et al. 2011); such approach misses the information contained in the higher moments of the bias distribution.
A full characterization of the halo scaling relations 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 coefficient11 ρ 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 and halos spin display the largest correlations, a trend that varies at low redshifts where halo concentration displays, interestingly, a larger correlation.The halo geometry, characterized by the ellipticity E h , shows a significant anticorrelation (with respect to intrinsic properties), similar in absolute value to that displayed by mass.
The main message from this plot is that nonlocal and environmental properties have higher-degree correlation with halo bias (compared to intrinsic properties).While this might be expected for the nonlocal halo properties (e.g., D 5 and ∆ 5 are probes of small-scale clustering and hence can be coupled via nonlinear 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 In Fig. 15 we present the correlation coefficient ρ 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 ρ s = 0 for almost all properties thought the redshift range.Example of secondary halo bias: mean halo effective bias ⟨b h |M vir ⟩ (θ)  q measured in bins of the halo (log) mass using a number of halo properties as secondary properties θ, namely, concentration C vir , spin λ B , ellipticity E h , Mach number M 5 , local dark matter overdensity ∆ dm , and tidal anisotropy T A .In each mass-bin, the sample has been divided in quartiles q of the property θ.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).
Ratio between the halo bias as a function of the halo mass (left column) and peak height (right column) obtained in two different quartiles of halo concentration C vir (top) and halo spin λ B (bottom).The shaded area denotes the uncertainty of the ratio computed from the error bars associated with each effective bias.

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 ⟨b h |θ, w i ⟩ at z = 0. Some conclusions based on these results are: -A statistically significant difference in the mean effective halo bias can be observed as a function of comic-web types, especially for environmental properties.-The mean relation between halo bias and intrinsic properties is positive for all cosmic-web types, except for voids, as shown in panels (a) to (f) of Fig. 16.-Halo effective bias displays a rather weak dependence with intrinsic properties such as concentration (panel b), spin (panel c) and triaxiality (panel e).-Nonlocal and environmental properties show negative bias not only in voids, but also in filaments and knots (see panels g-i).In general these set of properties display the largest differences in their links with effective bias as a function of the cosmic-web type.
-When explored in different cosmic-web types, the universality of the bias-peak height relation is waned, as shown in panel (j) of Fig. 16.-The halo effective bias is mildly higher in environments probing intermediate densities (i.e, filaments) than in extreme environments (knots and voids), as shown in panel (k) of Fig. 16.
One key feature arising from these results is the fact that we can measure, with statistical significance, a negative halo bias.Notice that negative values of mean effective bias are not reachable from estimators such as Eq.( 7), as these is based on positive definite quantities.On the other hand, such values should not be surprising, since it can arise from halo overdensities (δ h > 0) in voids (δ dm < 0).This does not imply that all objects living in voids acquire a negative bias.In fact, we have verified that nearly 60% of the halos living in voids (as defined by our cosmic-web classification) have a negative effective bias at all redshifts explored in this paper.
Article number, page 19 of 37 Statistical significance of the secondary halo bias as a function of halo mass, computed using Eq. ( 14) for different redshifts and for a number of secondary halo properties: concentration C vir , spin λ B , ellipticity E h , Mach number M 5 , local dark matter overdensity ∆ dm , and tidal anisotropy T A .The dark-shaded stripe denotes a strip of ∆ 4−1 p−s ± 2 and the horizontal line marks ∆ 4−1 p−s = 0. ∆ 4−1 p−s > 0 (< 0) implies that upper (lower) quartiles display higher clustering strength than lower (upper) quartiles.

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., Montero-Dorta et al. 2020;Sato-Polito et al. 2019;Salcedo et al. 2022, for a similar analysis, but using the correlation function computed from this type of subsets).For simplicity, we concentrate here in 6 representative secondary properties, namely, halo spin, Article number, page 20 of 37 Author: Balaguera-Antolínez, Montero-Dorta & Favole  take place (when more redshifts are shown).For example, notice how tracers with low halo concentration at z ∼ 2 are more biased than low concentrated halos at all masses, a trend that remains at z = 0 but only for very high halo mass.In general, this plot illustrates a well-known result: virtually all intrinsic and external halo properties display a certain level of secondary halo bias in some mass ranges.
We have verified that the mean effective bias measured in quartiles of primary properties such as the maximum circular velocity or the velocity dispersion is statistically compat-Article number, page 22 of 37 ible with that from the full sample, i.e, these type (namely, primary) of properties do not display significant secondary bias.In other words, a halo property displays secondary bias as long as it has a weak connection with halo mass (or peak height).
As we aim to measure signals of secondary bias for a 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 θ 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 vir [h −1 M ⊙ ] = 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 et al. 2019;Johnson et al. 2019) is not observed in UNIT-Sim 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 ν, 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 where ⟨b h |θ p ⟩ (s) i is the mean sample bias measured in bins of a primary property θ p in the i-th quartile of secondary property θ s .Similarly, σ s,i denotes the standard error in the mean bias obtained in the ith quartile of the secondary property.For a pair of primary and secondary properties, a value |∆ i− j p−s | ≫ 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, that is, 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.( 14) 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 ∆ i− j p−s = ±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: -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 ν, this crossover moves toward 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.
Article number, page 23 of 37 -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 vir [h −1 M ⊙ ] = 12.2 at z = 0).This behavior might be due to the presence of splashback halos (see e.g., Tucci et al. 2021) in the UNITSim.-Halo ellipticity.At fixed halo mass, spherical halos (lower quartile in E h ) are more clustered than nonspherical halos at z = 0.This trend appears to be inverted at z ≳ 2 at a mass scale of log M vir ∼ 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 toward the lower halo masses: at fixed halo mass, halos with higher M 5 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-M 5 halos are more clustered on large scales).In terms of cosmic-web types, again the presence of an inversion depends on environment (e.g., it is not observed clearly in sheets).-Neighbor 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 ν is partially observed (with the signals shifting toward lower masses).There are particular redshifts and cosmic-web environments for which small and low (statistical) -significant 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 vir [M ⊙ h −1 ] ∼ 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.
Article number, page 24 of 37 Author: Balaguera-Antolínez, Montero-Dorta & Favole -The redshift evolution of secondary bias is approximately (but not completely) captured by the implicit redshift dependence of ν 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 ν is chosen as a primary property (see also Fig. B.1 from the Appendix).

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.

The link between halo bias and halo properties
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 θ a (different from halo mass or peak height) as which implies that the statistical properties of halo bias as a function of θ 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 P(b h |M vir ) = P(b h |ν = ν 0 (M vir )), according to Eq.(3).Each of the conditional probability distributions under the integral in Eq. ( 15) is a marginalization over the halo properties; denoting with η = {η i } the remaining set of halo properties (i.e, not including mass or θ a ), this is expressed as where P(M vir |θ a , η 1 , • • • , η i ) represents a multivariate distribution of the halo mass conditional to the properties {η} (see e.g., Evrard et al. 2014).In general, the moments of the scaling relation P(M vir |θ a ) will be dependent of the correlation between viral mass and the set {η}.In particular, the mean and the scatter around it are a consequence of the propagation of the statistical properties of the links P(M vir |η i ).
We now evaluate the impact of partial marginalization over a reduced number of secondary properties in the scaling relations involved in Eq.(15).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 P(b h |M vir ) to the information provided by the Nbody simulation, such that the impact of secondary properties (including environmental properties) in the P(b h |θ a ) link will be addressed through the scaling relation P(M vir |θ a ).
We measure the scaling relation P(θ a |M vir , {ψ}), where {ψ} ∈ {η} is a subset of secondary properties12 .We use this measurement to reassign a new version of the property θ a (labeled θa ) to each tracer of the simulated volume, i.e θi a ↶ P(θ a |M vir = M i vir , ψ = ψ (i) ). (17) By construction, the new property θa obeys the same scaling relation with respect to the halo mass and the set {ψ}, while the correlations between θa and the rest of the halo properties (i.e, the set {η} − {ψ} ) 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.( 16).
We compare the properties of P(b h | θa ) (in particular, the mean) against those from the original scaling relation P(b h |θ a )13 .Figure 22 shows the result of reassigning the each of the properties θ s = {V max , C vir , λ B } using only the virial mass in Eq. ( 17).The resulting mean effective bias displays large deviations with respect to the true one, especially for halo concentration and spin.Interestingly, ∼ 20% deviations in the bias measured in bins of V max with respect to the original behavior leads to the conclusion that V 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 nonnegligible 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 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 {ψ} 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 of new halo catalogs with newly assigned properties using variations of the set ψ = {M 5 , ∆ dm , ∆ 5 , T/U, E h , ω 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 max , C vir , λ B } reassigned under different configurations (sets of {ψ}) 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), can be obtained without systematic bias when nonlocal 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 V and halo ellipticity E h can be key to retrieve the effective bias as a function of V 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 mentioned, again at z = 0. Left panels (top and bottom) in this plot show that the single mass-dependency used to sample the secondary property following Eq.( 17) cannot account for the signal of secondary bias, even if partial information is present in the form of scatter, as suggested by Eq. ( 16).Explicitly including more information such as local halo overdensity ∆ 5 or the cosmic-web classification can help to bring the set of reconstructed properties to the true signal of secondary bias, as shown in the different panels 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.

Toward applications to the construction of halo 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 so-called "bias assignment method" (Balaguera-Antolínez et al. 2019;Balaguera-Antolínez et al. 2020;Pellejero-Ibañez 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 to 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 ∆ 5 or the relative Mach number M 5 can help improve the accuracy of the reconstruction of halo fields.
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 Once bias is assigned, halo properties would then be assigned in a learn-and-map procedure.For example, virial mass can be assigned as after which V max can be assigned using 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 {η} = {η new } = {∆ 5 , ∆ dm , M 5 , D 5 , T A , ω 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).

Discussion and conclusions
For this work, we measured the signal of the secondary halo bias for the first time in the UNIT simulation (UNITSim), Article number, page 27 of 37 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 ∼ 6 to z = 0, in a comoving box of a 1 Gpch −1 side.The primary and secondary dependencies of the bias have been obtained using a new approach for the 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 with a 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 the halo bias opens the possibility to analyze the statistical behavior of the halo effective bias as a function of a number of halo properties classified as intrinsic, nonlocal, and environmental.
Intrinsic properties (mass, velocity dispersion, and halo spin) are derived from the halo finder producing the catalogs of the UNITSim.These can be divided into 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 directly probing the depth of the halo potential well, that is, the mass, velocity dispersion, and maximum circular velocity) and, on the other hand, on statistical grounds (a principal component analysis, shown in Appendix A, shows that these contain a higher level of variability, and therefore information, across time).
Nonlocal halo properties refer to small-scale statistical diagnostics aimed 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 Mpch −1 .Similarly, we have used the neighbor statistics, which condenses information of the distribution of pair separations with respect to each tracer 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 nonlocal (i.e., Article number, page 28 of 37 Author: Balaguera-Antolínez, Montero-Dorta & Favole tidal field).The computation of the tidal field of the darkmatter distribution allows for the exploration of the 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 the halo properties and effective bias as a function of cosmic-web environments and tidal anisotropy.
In the first part of this paper, we have reported the measurements of halo scaling relations in different cosmological environments.The 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 and maximum circular velocity) vary approximately along with expectations from, for example, the virial theorem, the scaling relations between the mass and secondary properties are more complex (e.g., these are not simple power laws) and, in general, they display a larger scatter around the mean relations.One particular example is the massconcentration relation, whose slope varies with mass and redshift, ranging from positive values at a high redshift (i..e, more massive halos are more concentrated, with a Spearman's rank correlation of ρ s ∼ 0.1, see Fig. 6), passing through a nearly mass-independent relation at z ∼ 2 (with ρ s ∼ 0.04) and ending with a negative correlation (i.e., more massive halos are less concentrated) toward a lower redshift (with ρ s ∼ −0.3).This in agreement with a number of findings in the literature (see e.  los is heavily influenced by its environment (see Figs. 7 and Fig. 8).
The second (and main) part of this paper has been devoted to the assessment of the statistical properties of the large-scale halo bias as a function of the properties described in the first section of this work.The conclusions can be summarized as follows: -The object-by-object halo assignment proposed by Paranjape et al. (2018) yields estimates for the largescale halo bias in agreement with estimates based on the halo power spectrum (see Fig. 10).This approach has the great advantage of no explicit need to compute two-point statistics for subsamples divided into 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 the 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 the 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 into terms of cosmic-web environments.-In general, the effective halo bias expressed as a function of halo properties is quite sensitive to the cosmic-web environment (see Fig. 16).-The signal of the secondary bias is present when exploring the bias as a function of the 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 preserve a hierarchy across cosmic time, that is, the difference between the bias in the higher quartile and that in the lower quartile does not change sign with redshift, though the trend might differ from property to property.For instance, high spin halos have a larger bias than low-spin halos for all masses and redshifts explored ( §4.5).-Among the secondary properties, the signal of the secondary bias due to the halo concentration displays a well-known behavior, inverting the hierarchy at a different mass scale in a different snapshot, a fact seen in previous works (see e.g., Wechsler et al. 2006;Sato-Polito et al. 2019).Such a 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 the clustering trend in the secondary bias induced by the halo spin, as previously reported (e.g., Sato-Polito et al. 2019;Tucci et al. 2021).However, this type of crossover is observed when the analysis is conditioned to halos located in high density regions (knots).-Nonlocal 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 the secondary bias is approximately (but not completely) captured by the implicit redshift dependence of ν 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 ν is chosen as a primary property.
We have demonstrated that the signal of the 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 associated with high (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 as to the physical origin of the secondary bias.Our results conceptually agree with those derived from the analysis of Wang et al. (2021) and (Wang et al. 2023), which implemented different estimators for the large-sale halo (or galaxy) bias, definitions of halo and galaxy properties, and a slightly different set-up for cosmological simulations, including the halo finder.
Evidently, more effort is required to fully understand the physical origin of the secondary bias.In that regard, we shall finish by providing some information on forthcoming research which extends some of the results presented here.These include the following: -The object-by-object estimator of Eq.( 10) can be generalized to provide estimates of individual redshift-space effective halo bias.Such a generalization can be used to test different subjects linked to redshift distortions, Article number, page 30 of 37 Author: Balaguera-Antolínez, Montero-Dorta & Favole such as the velocity bias, measurement of the growth of structures, and marked statistics in redshift space.-Extensions to redshift surveys can also be complemented by assigning 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 the underlying dark matter density field has been reconstructed using a number of techniques (see e.g., Hada & Eisenstein 2018;Jasche & Lavaux 2019;Kitaura et al. 2021).Assigning a relative bias can be important 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;Balaguera-Antolínez 2014;Clerc & Finoguenov 2023); we can assign the relative effective bias using a single measurement for the cluster power spectrum from the the full cluster population.-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 TNG 14 ) in the context of galaxy assembly bias (i.e., the secondary dependencies of galaxy clustering and occupancy at a 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 that have yet to be understood.-Finally, the full applicability of the procedures described in §5.2 is focused on the reconstruction of the halo bias and its dependencies with halo properties to generate halo mock catalogs.This application aims to extend the procedures shown in Balaguera-Antolínez et al. (2023) which attempts to assign halo properties to mock tracers obtained from learning-based techniques, which is dependent on a precise signal for large-scale clustering.large scatter in the bias-mass relation (see Fig. 11), defining a subsample using a cut in effective (or relative) can induce statistical incompleteness, especially on large separations (see e.g., Balaguera-Antolínez et al. 2011).To account for this, we randomly redistribute the ranked marks within the tracers in the subsample and measure the marked correlation functions.We repeat this procedure 100 times, generating a set from which we compute a mean and a variance, denoted ⟨M 0 (r)⟩ and σ r (r j ) respectively.We report the measurements of M(r) ≡ M(r)/⟨M 0 (r)⟩, such that the value M(r) = 1 marks a scale in which the bias distribution is characterized by the mean bias of the sample.The uncertainty in this ratio is computed using σ r (r j ) and the uncertainty in M(r), which is obtained from a one-time deleted Jackknife technique (see e.g., Tukey 1958;Norberg et al. 2009;Favole et al. 2021, and references therein) using 64 sub-volumes.These two are combined in the form of error propagation neglecting the covariance between these two components. 16n Fig. C.2 we show the measurements of the quantity M(r) probing separations in the range 10 < r < 250 Mpc h −1 at two different redshifts, using both the effective and relative halo bias.We present the measurements from the high-bias tracers (b h ≫ ⟨b h ⟩) and the extremely low biased tracers (b h ≪ ⟨b h ⟩).All these follow a particular trend, namely, the bias of high (low) biased tracers depicts a positive (negative) deviation form the mean, reaching M(r) = 1 at r ∼ 50 Mpch −1 .An oscillatory pattern follows, with a period of nearly 40 Mpc h −1 and a decaying amplitude.This is more evident at high redshift, where actually low and high bias signals tend to be in phase, for both absolute and relative bias.In any case, we can interpret such oscillatory behavior as a reminiscent of Gibbs phenomena linked to the low pass filter used to define the effective (or relative) halo bias.Article number, page 37 of 37

Fig. 2 .
Fig. 2. Probability distribution of the 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 λ th = 0 in the cosmicweb classification, we identified under-densities with voids and overdensities with knots.

Fig. 3 .
Fig. 3. Distribution of intrinsic {M vir , V max , R s , C vir , λ B , T/|U|, T h , E h }, nonlocal {M 5 , D 5 , ∆ 5 }, and environmental (T A ) halo properties in the UNITSim at z = 0.In all panels, the black solid line represents the distribution from the full population.We also present the property distribution in the cosmic-web types defined in §3.3.1..

Fig. 4 .
Fig. 4. Mean scaling relations between halo virial mass and a number of halo properties described in §3.1, at different redshifts.Shown intrinsic properties are V max (panel a), R s (panel b) C vir (panel c), λ B (panel d), T/|U| (panel e), T h (panel f) , E h (panel g).Nonlocal properties: local Mach number M 5 (panel h), ∆ 5 (panel i), D 5 (panel j).Environmental properties: local dark matter density (panel k) and tidal anisotropy T A (panel l).The error bars denote the error in the mean in each mass bin.

Fig. 5 .
Fig. 5. Volume filling fraction (VFF) and the fraction of tracers f tr in each cosmic-web type defined in §3.3.1 as a function of redshift.

Fig
Fig. 6.Spearman's rank correlation coefficient ρ s between the halo effective bias and a number of halo properties for two different redshifts.

Fig. 7 .
Fig. 7. Spearman's rank correlation coefficient ρ s between the halo effective bias and a number of halo properties at z = 0.The upper diagonal shows the correlation computed in knots, lower diagonal in voids.

Fig. 8 .
Fig. 8. Mean scaling relations between the halo virial mass and a number of halo properties described in §3.1, at z = 0 evaluated in different cosmic-web types.Shown intrinsic properties are V max (panel a), R s (panel b) C vir (panel c), λ B (panel d), T/|U| (panel e), T h (panel f), E h (panel g).Nonlocal properties: local Mach number M 5 (panel h), ∆ 5 (panel i), D 5 (panel j).Environmental properties: local dark matter density (panel k) and tidal anisotropy T A (panel l).The error bars denote the error in the mean in each mass bin.

Fig. 9 .
Fig. 9. Slice of thickness ∼ 80 Mpc h −1 of the UNITSim (at z = 0) showing the spatial distribution of the effective halo bias, as obtained form the interpolation of the bias assigned to each halo according to Eq.(10) on 512 3 mesh.

Fig. 10 .
Fig. 10.Halo effective bias b hm computed with Eq.(10) at z = 0. Left panel: Bias as a function of the halo peak height; the points denote the mean bias in different ν bins, and the error bars denote the standard error of the mean in each bin.The dotted line shows the prediction of Tinker et al. (2010).Right panel: Halo bias b hm as a function of the halo virial mass.In both cases, the contours indicate a region of an equal number of tracers N H .

Fig
Fig. 12.Mean halo bias as a function of the 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.

Fig. 13 .
Fig. 13.Second and third order moments (variance, skewness, and kurtosis) of the halo scaling relation between the effective bias and peak height P(b h |ν), computed for different redshifts.The right column shows the behavior in high density regions (knots).

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

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

Fig. 16 .
Fig. 16.Mean effective bias ⟨b h |θ, w i ⟩ measured in bins of halo properties θ in the different cosmic-web ω i at z = 0. Intrinsic halo properties are: virial mass M vir (panel a), halo concentration C vir (panel b), halo spin λ B (panel c), the ratio T/|U| (panel d), halo triaxiality T h (panel d), halo ellipticity E h (panel g).Nonlocal properties: relative Mach number M 5 (panel g), local overdensity ∆ 5 (panel h) and neighbor statistics D 5 (panel i).Panel (j) shows the peak height-halo bias relation.Environmental properties: local dark matter density (panel k) and tidal anisotropy T A (panel l).Article number, page 17 of 37

Fig. 20 .
Fig. 20.Secondary halo bias in different cosmic-web environments at z = 0. From top to bottom: halo concentration C vir , spin λ B , ellipticity E h , Mach number M 5 , local overdensity ∆ dm and tidal anisotropy T A ).

Fig. 22 .
Fig. 22.Comparison between the measurements of the effective halo bias as a function V max , C vir , and λ B as read from the original catalog against the signal measured used after reassigning the same set of properties using a scaling relation based solely on the virial mass.The bottom panels show the ratio between the bias in the bins of the "reassigned" to that measured in bins of the "original" property.

Fig. 23 .
Fig. 23.Same as Fig. 22 but focused on V max .In each panel, different dependencies (named in the legends) are used to reassign maximum circular velocity.

Fig. 24 .
Fig. 24.Same as Fig. 22 but using more information in the scaling relation to assign the property θs (halo concentration in this case) as shown in each panel.

Fig. 25 .
Fig. 25.Same as Fig. 22 but using more information in the scaling relation to assign the property θs (halo spin in this case) as shown in each panel.

Fig. 26 .
Fig. 26.Reconstruction of the secondary bias: mean effective halo bias measured in bins of the 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 reassigned using Eqs.(17) with different combinations of nonlocal or environmental properties, including the information of tidal field T through the cosmic-web classification.

Fig. 27 .
Fig. 27.Same as Fig. 26 but for halo spin as secondary property.
g.,Macciò et al. 2007;Knebe & Power 2008).-Thevariation in the halo spin with the halo mass at different redshifts is weak for low-mass halos (log M vir ∼ 11.5).In particular, at z = 0 (with ρ s ∼ −0.01),only the very massive halos (log M vir > 14) display a relevant (negative) correlation with the halo mass.The dependency is slightly stronger at a higher redshift (e.g., at z ∼ 5, ρ s ∼ −0.02).-The halo-scaling relations involving secondary properties vary across the cosmic-web environment.Although the definition adopted based on the Hessian of the gravitational potential is arbitrary (as it depends on the threshold λ th fixed to zero in this paper), the clear distinction of the scaling relations in over-and underdensities indicates how the evolution and growth of ha-Article number, page 29 of 37

Fig. 28 .
Fig.28.Effective halo bias measured in a synthetic halo catalog with bias and mass reassigned according to the learn-andmap approach ( §5.2).The assignment of halo bias used the set of environmental properties {∆ 5 , ∆ dm , M 5 , D 5 , T A , ω i }, while halo masses are drawn from the scaling relation P(M vir |b h , ∆ dm , ∆ 5 , ω i ).

Fig. C. 1 .
Fig. C.1.Slices of 50 Mpch −1 thickness in an area of ∼ 300 (Mpch −1 ) 2 of the UNITSim volume, showing the dark matter density, the halo Mach number, and the effective bias (at z = 0) interpolated on a 512 3 mesh.It is evident how halo properties, in this case represented by the Mach number, follows the filamentary structure imposed by the dark matter density field, while the effective bias field filters out scales below ∼ 60 Mpc h −1 .

Fig
Fig. C.2. Ranked mark correlation function M(r) using the ranks of the effective bias as a mark for two different redshifts in the UNITSim.The dots with error bars show both the result using tracers with effective bias much larger (or lower) than the mean bias of the sample.The gray lines denote the measurements of the ranked-mark correlation function after randomizing the marks among the tracers.