PAPER • OPEN ACCESS

# Taming assembly bias for primordial non-Gaussianity

To cite this article: Emanuele Fondi et al JCAP02(2024)048

View the [article online](#) for updates and enhancements.

## You may also like

- **DETECTION OF THE SPLASHBACK RADIUS AND HALO ASSEMBLY BIAS OF MASSIVE GALAXY CLUSTERS**
  Surhud More, Hironao Miyatake, Masahiro Takada et al.

- **Assembly bias in quadratic bias parameters of dark matter halos from forward modeling**
  Titouan Lazeyras, Alexandre Barreira and Fabian Schmidt

- **Toward Accurate Modeling of Galaxy Clustering on Small Scales: Halo Model Extensions and Lingering Tension**
  Gillian D. Beltz-Mohrmann, Adam O. Szewciw, Andreas A. Berlind et al.
RECEIVED: November 23, 2023

REVISED: January 23, 2024

ACCEPTED: January 24, 2024

PUBLISHED: February 26, 2024

# Taming assembly bias for primordial non-Gaussianity

Emanuele Fondi<sup>a</sup>, Licia Verde<sup>a,b</sup> Francisco Villaescusa-Navarro<sup>c,d</sup>

Marco Baldi<sup>e,f,g</sup> William R. Coulton,<sup>h,i,d</sup> Gabriel Jung<sup>j</sup> Dionysios Karagiannis<sup>k</sup>

Michele Liguori<sup>l,m</sup> Andrea Ravenni<sup>l,m</sup> and Benjamin D. Wandelt<sup>n,c</sup>

<figure><img src="image_1.png" /></figure>

<footnote><sup>a</sup>ICC, University of Barcelona, IEEC-UB, Martí i Franquès, 1, E08028 Barcelona, Spain</footnote>

<footnote><sup>b</sup>ICREA, Pg. Lluís Companys 23, Barcelona, 08010, Spain</footnote>

<footnote><sup>c</sup>Center for Computational Astrophysics, 162 5th Avenue, New York, NY, 10010, U.S.A.</footnote>

<footnote><sup>d</sup>Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, U.S.A.</footnote>

<footnote><sup>e</sup>Dipartimento di Fisica e Astronomia, Alma Mater Studiorum — University of Bologna, Via Piero Gobetti 93/2, 40129 Bologna BO, Italy</footnote>

<footnote><sup>f</sup>INAF — Osservatorio Astronomico di Bologna, Via Piero Gobetti 93/3, 40129 Bologna BO, Italy</footnote>

<footnote><sup>g</sup>INFN — Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna BO, Italy</footnote>

<footnote><sup>h</sup>Kavli Institute for Cosmology, Madingley Road, Cambridge, CB3 0HA, U.K.</footnote>

<footnote><sup>i</sup>DAMTP, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, U.K.</footnote>

<footnote><sup>j</sup>Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405, Orsay, France</footnote>

<footnote><sup>k</sup>Department of Physics & Astronomy, University of the Western Cape, Cape Town 7535, South Africa</footnote>

<footnote><sup>l</sup>INFN, Sezione di Padova, Via Marzolo 8, I-35131, Padova, Italy</footnote>

<footnote><sup>m</sup>Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova, Via Marzolo 8, I-35131, Padova, Italy</footnote>

<footnote><sup>n</sup>Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>
E-mail: emanuelefondi@icc.ub.edu, liciaverde@icc.ub.edu, villaescusa.francisco@gmail.com, marco.baldi5@unibo.it, william.coulton@gmail.com, gabriel.jung@protonmail.com, dakaragian@gmail.com, michele.liguori@unipd.it, andrea.ravenni@unipd.it, bwandelt@iap.fr

**ABSTRACT:** Primordial non-Gaussianity of the local type induces a strong scale-dependent bias on the clustering of halos in the late-time Universe. This signature is particularly promising to provide constraints on the non-Gaussianity parameter $f_{\text{NL}}$ from galaxy surveys, as the bias amplitude grows with scale and becomes important on large, linear scales. However, there is a well-known degeneracy between the real prize, the $f_{\text{NL}}$ parameter, and the (non-Gaussian) assembly bias i.e., the halo formation history-dependent contribution to the amplitude of the signal, which could seriously compromise the ability of large-scale structure surveys to constrain $f_{\text{NL}}$. We show how the assembly bias can be modeled and constrained, thus almost completely recovering the power of galaxy surveys to competitively constrain primordial non-Gaussianity. In particular, studying hydrodynamical simulations, we find that a proxy for the halo properties that determine assembly bias can be constructed from photometric properties of galaxies. Using a prior on the assembly bias guided by this proxy degrades the statistical errors on $f_{\text{NL}}$ only mildly compared to an ideal case where the assembly bias is perfectly known. The systematic error on $f_{\text{NL}}$ that the proxy induces can be safely kept under control.

**KEYWORDS:** cosmological parameters from LSS, cosmological simulations, galaxy surveys, inflation

ARXIV ePRINT: 2311.10088
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Theoretical background</td><td>3</td></tr><tr><td>2.1</td><td>b<sub>φ</sub> and the universality relation</td><td>5</td></tr><tr><td>2.2</td><td>Halo assembly bias and its relation with b<sub>φ</sub></td><td>6</td></tr><tr><td>3</td><td>N-body simulations and comparison with ePS-derived expressions</td><td>7</td></tr><tr><td>3.1</td><td>Describing halo formation history through ω<sub>f</sub> and P<sub>ωf</sub></td><td>9</td></tr><tr><td>3.2</td><td>Accuracy of P<sub>ωf</sub> modeling: results on Δb<sub>φ</sub> from non-Gaussian simulations</td><td>12</td></tr><tr><td>4</td><td>Galaxy proxy for halo assembly bias: results from IllustrisTNG simulations</td><td>15</td></tr><tr><td>4.1</td><td>Correlations between galaxy photometry and halo formation redshift</td><td>15</td></tr><tr><td>4.2</td><td>Optimal proxy for ω<sub>f</sub> from a combination of galaxy colors</td><td>16</td></tr><tr><td>5</td><td>Marginalizing over cosmology and astrophysics: tests on CAMELS-TNG</td><td>19</td></tr><tr><td>6</td><td>Forecasts on f<sub>NL</sub> measurement</td><td>21</td></tr><tr><td>6.1</td><td>Single tracer: statistical error</td><td>23</td></tr><tr><td>6.2</td><td>Single-tracer systematic error</td><td>24</td></tr><tr><td>6.3</td><td>Multi-tracer</td><td>26</td></tr><tr><td>7</td><td>Conclusions</td><td>28</td></tr><tr><td>A</td><td>Optimizing the number of galaxy colors combined in Y</td><td>31</td></tr></table>

## 1 Introduction

Inflation represents the standard paradigm for the description of the early evolution of our Universe. In its simplest form, it is determined by a single slow-roll scalar field and produces an almost Gaussian [1] distribution of primordial perturbations, which evolve gravitationally to form structure in the late Universe. Although current constraints from Cosmic Microwave Background (CMB) [2] and large-scale structure (LSS) [3] observables confirm the Gaussian scenario, the enhanced sensitivity of ongoing [4, 5] and forthcoming [6, 7] LSS experiments may have the potential to reveal the non-Gaussianity of the initial conditions, predicted by alternative inflationary models [8–13].
In particular, the presence of multiple fields during inflation would introduce non-linearities in the primordial Bardeen potential, generating what is known as Primordial non-Gaussianity (PNG) of the local type [14–16]. Local PNG can be expressed in its simplest form as a quadratic correction with amplitude $f_{\text{NL}}$ [17, 18]:

$$ \Phi(\mathbf{x}) = \phi(\mathbf{x}) + f_{\text{NL}} (\phi(\mathbf{x})^2 - \langle \phi(\mathbf{x})^2 \rangle) \quad (1.1) $$

where $\Phi(\mathbf{x})$ is the Bardeen potential and $\phi(\mathbf{x})$ is an auxiliary Gaussian random field.

The primordial bispectrum generated by local PNG couples small scale density perturbations with large scale potential modes, affecting dark matter halo clustering at late times. As a result, the halo bias takes a scale dependence on large scales, with an amplitude given by $b_{\phi}f_{\text{NL}}$, where $b_{\phi}$ is the local PNG bias parameter [19–23]. Detecting a nonzero $b_{\phi}f_{\text{NL}}$ would be sufficient to rule out the standard single-field, slow-roll inflationary model [24]. However, in order to constrain the value of $f_{\text{NL}}$ and shed light on the physics of inflation, discriminating between different models, a prior knowledge of $b_{\phi}$ is required.

Conventionally, predictions for $b_{\phi}$ rely on the universality relation, derived under the assumption of a halo mass function (HMF) which is universal, i.e., dependent on halo mass and redshift only through a single variable [25–27]. However, a number of recent results from N-body simulations highlighted the breakdown of this relation, due to assembly bias in $b_{\phi}$, i.e., its dependence on properties beyond the halo mass. Specifically, at fixed mass, halos with different concentration have significantly different $b_{\phi}$ [28–30]. Similar conclusions have been reported for samples of galaxies selected by their stellar mass, color and other galaxy properties [31–33].

Recognizing these departures from universality paved the way for further studies in various directions. In particular, the theoretical uncertainties called for new efforts to provide better priors on $b_{\phi}$ with the aid of machine learning techniques [30, 34], or to break the $b_{\phi}f_{\text{NL}}$ degeneracy introducing new summary statistics [35]. Furthermore, since samples with different $b_{\phi}$ represent an opportunity to employ the multi-tracer technique [36] and provide better $f_{\text{NL}}$ constraints, some analyses focused on optimizing sample selection strategies [29, 30, 37].

Although recently highlighted in many works, these issues were already pointed out by the seminal papers in the field [21, 36, 38, 39]. In particular, the attention was focused on how details of the halo formation history affect the value of $b_{\phi}$. Specifically, ref. [21] showed that at fixed mass, halos which assembled their mass recently have a $b_{\phi}$ lower than predicted by the universality relation. This result was generalized by ref. [39], where the authors effectively found an analytical relation between $b_{\phi}$ and halo formation time, by employing the extended Press-Schechter (ePS) formalism [40–42]. In the same way as the Press-Schechter formalism can be used to compute a universal HMF, its extension gives the conditional mass function (CMF) that describes the halo mass accretion history. The CMF can be used to estimate the departure of $b_{\phi}$ from the universality relation.

In this paper we use different suites of N-body simulations [43–46] to test the validity of the analytical ePS prediction, across various halo mass ranges and redshifts. To take into account inconsistencies, we suggest a 1-parameter extension of the ePS, calibrated on simulations. We then resort to hydrodynamical simulations [47–49] which also model galaxy formation to further connect the modeling of $b_{\phi}$ and assembly bias to observational
properties of galaxies. In particular, we study photometric properties of simulated galaxies to find an observational proxy for halo assembly bias and assess its robustness to changes in cosmological and astrophysical parameters. Finally, we report $f_{\text{NL}}$ error forecasts for single and multi-tracer analyses, providing insights on optimal galaxy sample selection strategies.

The rest of the paper is structured as follows. In section 2 we review the theoretical background of the non-Gaussian halo bias and the ePS prediction for $b_\phi$. In section 3 we test the ePS predictions on N-body simulations and calibrate a simple extension of the analytical result, which shows a remarkable good fit across different halo mass and redshift ranges. In section 4 we consider different galaxy samples and identify the optimal proxy for halo assembly bias from combinations of photometric galaxy colors. The robustness of this method is assessed in section 5 by using a different suite of simulations, which consider a broad range in the astrophysical parameters. Focusing on ELG-like and LRG-like simulated samples, we report Fisher forecasts for the precision of the $f_{\text{NL}}$ measurement from single and multi-tracer analyses in section 6. We finally summarize our main results and conclude in section 7.

## 2 Theoretical background

In this section we summarize the relevant results from the literature. While no novel results are presented, the section serves to lay out useful equations and formulae and to provide a unified notation.

Predictions for the distribution and clustering of dark matter halos in the presence of PNG can be obtained using conceptually distinct but mutually consistent approaches [19–22]. In the peak-background split formalism [27], the auxiliary Gaussian field $\phi$ in eq. 1.1 is decomposed into statistically independent short- and long-wavelength modes $\phi = \phi_l + \phi_s$.

The gravitational potential fluctuations described by eq. 1.1 are translated into (linear) matter density fluctuations at redshift $z$ through the Poisson equation

$$ \delta_{\text{m}}(\mathbf{k}, z) \equiv \mathcal{M}(k, z)\Phi(\mathbf{k}) = \frac{2c^2k^2T(k)}{3\Omega_mH_0^2}D(z)\Phi(\mathbf{k}) \quad (2.1) $$

where $D(z)$ denotes the linear growth rate of perturbations normalized to be $(1+z)^{-1}$ during matter-domination and $T(k)$ being the matter transfer function normalized to unity on large scales.

In the absence of PNG, small- and large-scale linear density fluctuations would evolve independently under the Poisson equation, while in the presence of a non-zero $f_{\text{NL}}$ in eq. 1.1, the short wavelengths in the density field become:

$$ \delta_s = \mathcal{M}(k, z)[(1 + 2f_{\text{NL}}\phi_l)\phi_s + f_{\text{NL}}\phi_s^2]. \quad (2.2) $$

Thus small scale density fluctuations are modulated by large scale potential modes, with a coupling given by $f_{\text{NL}}$. In turn, the number of collapsed halos will be affected, determined both by density and potential fluctuations. Writing the halo overdensity $\delta_h$ as a linear bias expansion, we can therefore express it in terms of the two fields $\delta_m$ and $\Phi$ as [23]

$$ \begin{aligned} \delta_h(\mathbf{k}|M, z) &= b_1(M, z)\delta_m(\mathbf{k}, z) + b_\phi(M, z)f_{\text{NL}}\Phi(\mathbf{k}) + \epsilon(\mathbf{k}, z) \\ &= \left[ b_1(M, z) + \frac{b_\phi(M, z)f_{\text{NL}}}{\mathcal{M}(k, z)} \right] \delta_m(\mathbf{k}, z) + \epsilon(\mathbf{k}, z) \end{aligned} \quad (2.3) $$
where $\epsilon$ denotes a shot noise (or stochastic) term and $b_1$ and $b_\phi$ are bias parameters. The first bias term in the bracket is scale-independent, while the second one, which depends on $k$ via $\mathcal{M}(k, z)$ represents the so-called scale-dependent bias. Both bias parameters can be related to the halo mass function $n(M, z)$, i.e. the number density of halos in a mass interval $dM$ around mass $M$ and redshift interval $dz$ around redshift $z$.

In the Gaussian case, the linear Lagrangian bias reduces to the scale-independent term:

$$b_1^G(M, z) = \frac{\partial \ln n^G(M, z)}{\partial \delta_l} = -\frac{\partial \ln n^G(M, z)}{\partial \delta_c} \qquad (2.4)$$

where $n^G$ denotes the Gaussian halo mass function. The second equality arises because the effect of modulating the density field by a long wavelength mode $\delta_l$ in some patch of the Universe can be seen as an additive change in the critical density for collapse $\delta_c$ in that region [50].

In the presence of non-Gaussianity there are two main corrections. The first one regards the scale-independent bias parameter $b_1$ [51, 52]:

$$b_1(f_{\text{NL}}, M, z) = b_1^G(M, z) + \Delta b_I(f_{\text{NL}}, M, z) = \frac{\partial \ln n^{\text{NG}}}{\partial \delta_l} = -\frac{\partial \ln n^G}{\partial \delta_c} - \frac{\partial \ln R^{\text{NG}}}{\partial \delta_c} \qquad (2.5)$$

where $R^{\text{NG}}(f_{\text{NL}}, M, z) = \frac{n^{\text{NG}}(f_{\text{NL}}, M, z)}{n^G(M, z)}$ denotes the ratio between the non-Gaussian mass function and the Gaussian one, that can be modeled in different ways [53–56]. The scale independent correction $\Delta b_I$ is due to the fact that the modulation induced by $f_{\text{NL}}$ changes the mean number density of collapsed objects $\bar{n}$ compared to the Gaussian case.

The second correction regards the scale-dependent bias term mentioned above. It can be computed within the peak-background split formalism [21] by making use of eq. (2.2), from which

$$\Delta b(k, f_{\text{NL}}, M, z) = 2 f_{\text{NL}} \frac{d\phi_l}{d\delta_l} \frac{\partial \ln n^G}{\partial(1 + 2f_{\text{NL}}\phi_l)} = \frac{2f_{\text{NL}}}{\mathcal{M}(k, z)} \frac{\partial \ln n^G}{\partial \ln \sigma_8^{\text{local}}} \qquad (2.6)$$

where we have used the fact that locally, the effect of non-Gaussianity is a rescaling of the small scale matter fluctuations $\sigma_8^{\text{local}} = (1 + 2f_{\text{NL}}\phi_l)\sigma_8$ and thus $d\sigma_8^{\text{local}} = \sigma_8 d(1 + 2f_{\text{NL}}\phi_l)$. This is known as separate universes argument.

This correction becomes dominant at large scales, due to the fact that $T(k \to 0) = 1$ and thus $\mathcal{M}(k, z) \propto k^2 D(z)$. Its amplitude is determined by the bias coefficient $b_\phi$ which, as we will show below, depends on the properties of the population of tracers considered. By comparing eq. (2.3) with (2.6) we can identify $b_\phi$ as

$$b_{\phi} = 2 \frac{\partial \ln n}{\partial \ln \sigma_8}. \qquad (2.7)$$

where we have dropped the G superscript and the local label for brevity. As pointed out by [21], we note that although $b_\phi$ drives the PNG-dependent contribution to the bias, it is expressed entirely in terms of Gaussian quantities, except for a cosmology with a suitably rescaled $\sigma_8$.
## 2.1 $b_\phi$ and the universality relation

As shown above, in order to provide predictions for $b_\phi$, it is sufficient to specify the halo mass function and measure its response to changes in $\sigma_8$. The standard result is obtained assuming a universal mass function [21, 25–27]

$$n(M, z) = \frac{\bar{\rho}}{M^2} \nu f(\nu) \frac{d \ln \nu}{d \ln M} \qquad (2.8)$$

where $\bar{\rho}$ is the background density, $\nu = \frac{\partial_c(z)}{\sigma(M)}$ is the ratio between the redshift-dependent collapse threshold $\mathfrak{D}_c(z) = \frac{\delta_c}{D(z)}$ and $\sigma(M)$, the amplitude of fluctuations in the linear density field, given by

$$S(M) \equiv \sigma^2(M) = \frac{1}{2\pi^2} \int_0^\infty P_{\text{lin}}(k) W_M^2(k) k^2 dk \qquad (2.9)$$

where $P_{\text{lin}}(k)$ is the power spectrum of the linear density field, which is smoothed with a top-hat filter $W_M(k)$ on a scale enclosing the mass $M$. In the Press-Schechter (PS) formalism, the number density of collapsed objects $n(M, z)$ is stated in terms of a first-crossing problem: at each spatial point, the trajectory $\delta(M)$ of the smoothed linear matter density as a function of the smoothing scale $M$ undergoes a random walk; the upcrossing of the barrier $\mathfrak{D}_c(z)$ corresponds then to the collapse of matter into halos. The solution of this problem is given by the first-crossing probability distribution $f(\nu)$. A number of different functional forms of $f(\nu)$ have been proposed in the literature [26, 27, 57, 58], depending on the choice of the filter and the barrier. However, its exact form does not need to be specified to compute $b_\phi$, as long as it is universal, i.e., only depends on mass and redshift via the variable $\nu$. For such a mass function, the prediction for $b_\phi$ (using the u superscript to indicate the assumption of a universal mass function) is [21]

$$b_{\phi}^{u} = 2\delta_{c}(b_{1} - 1) \qquad (2.10)$$

and agrees with previous predictions derived independently [19, 20].

As remarked in a number of recent works [28, 31–33] and already pointed out in the seminal papers in the field [21, 36, 38, 39], this result suffers from two main problems.

First, the value of $\delta_c$ depends on the assumptions made about the formation (collapse) of halos. For example, assuming a spherical collapse, as considered in the PS formalism [26], $\delta_c = 1.686$. This assumption, however, might not be robust and motivated the development of non-spherical collapse models, characterized by a non-constant barrier [57, 59–61]. Furthermore, as these predictions are tested against N-body simulations, the actual definition of halos also depends on the halo finder algorithm considered (as FoF, SO), which affects the halo number counts and thus their bias. Therefore, in order to include all these effects in a “fudge factor” $q$, a modification $\delta_c \to q\delta_c$ is often considered [51, 52, 62–64].

Moreover, and crucially in view of comparisons with observations, the result is based on the assumption of having a fair sample of all the halos in a given mass range. Equation (2.10) does not hold, for example, for a sample of objects populating only recently merged halos [21]. Therefore, it has been argued that it should be rather modified as

$$b_{\phi}^{s} = 2\delta_{c}(b_{1} - p) \qquad (2.11)$$
where we explicitly indicate that $b_\phi$ is specific to the sample $s$, and the sample characteristics define the value of $p$, with $p = 1$ for a fair sample of halos and $p = 1.6$ for recent mergers [21, 36].

## 2.2 Halo assembly bias and its relation with $b_\phi$

Dark matter halo tracers (e.g., galaxy) properties are well known to depend primarily on the host halo mass, however there are secondary parameters modulating galaxy formation, such as environment and halo formation history. If the selection of a sample of tracers depends crucially on the halo formation history (e.g., quasars, that are triggered by recent mergers), then eq. (2.10) will not hold.

Following [21, 39] we explicitly recognize that in the extended Press-Schechter (ePS) approach [40–42] it is possible to use the conditional mass function to make explicit the dependence of $b_\phi$ on the halo formation time. This can be done by considering the number density of halos of mass $M_o$, observed at redshift $z_o$, that had accreted<fnref n="1" /> a fraction $f$ of their final mass by formation time $z_f$. The $b_\phi$ associated to such halos is given by

$$b_{\phi}(M_o, z_o, z_f, f) = 2 \frac{\partial \ln n(M_o, z_o)}{\partial \ln \sigma_8} + 2 \frac{\partial \ln P_{z_f}(fM_o, z_f|M_o, z_o)}{\partial \ln \sigma_8} = b_{\phi}^u + \Delta b_{\phi}, \quad (2.12)$$

which follows from eq. (2.7). $P_{z_f}$ denotes conditional mass function: the probability that a halo with mass $M_o$ at $z_o$ had a mass $fM_o$ in a redshift interval $dz_f$ around a higher redshift $z_f$.

The first term, under the assumption of a universal mass function, reduces to eq. (2.10), hence if the generic parameterization of $b_\phi$ in terms of $p$ is considered, as in eq. (2.11), we can identify the second term with $\Delta b_\phi = 2\delta_c(1-p)$. Thus the conditional mass function, i.e., the merger history of the objects selected, can be used to put a prior on the value of $p$, or more generally on $b_\phi$. For a fair sample of halos of mass $M$, $\Delta b_\phi = 0$ and $p = 1$, otherwise it is a potentially very important correction [28, 31–33, 39].

As mentioned above, in the PS formalism the halo mass function is expressed in terms of the distribution $f(\nu)$ of first crossings of a barrier $\mathfrak{D}_c(z)$. Following [40–42], the conditional mass function can be computed in the same way, in the context of a diffusion problem between two barriers $\mathfrak{D}_c(z_o)$ and $\mathfrak{D}_c(z_f)$. Since for random walks the crossing probability does not depend on the path, the solution of this problem has the same form as $f(\nu)$, with a transformation of the origin of coordinates. In the specific case of the PS solution [26], expressed in terms of $S_f = S(fM_o)$ and $S_o = S(M_o)$ as defined by eq. (2.9), we have

$$\frac{df}{dS_o}(S_f, z_f|S_o, z_o) = \frac{1}{\sqrt{2\pi}} \frac{\mathfrak{D}_c(z_f) - \mathfrak{D}_c(z_o)}{(S_f - S_o)^{3/2}} e^{-\frac{|\mathfrak{D}_c(z_f) - \mathfrak{D}_c(z_o)|^2}{2(S_f - S_o)}} \quad (2.13)$$

which gives the conditional probability for objects of a given mass $M_o$. Then, we use eq. (2.8) to convert mass fractions to number of halos and define

$$\omega_f \equiv \frac{\mathfrak{D}_c(z_f) - \mathfrak{D}(z_o)}{\sqrt{\sigma^2(fM_o) - \sigma^2(M_o)}}, \quad (2.14)$$

which is the variable that we use to parameterize assembly bias and is proportional to $z_f - z_o$ at fixed $M_o$, as described in detail in section 3.1. After some calculations, we obtain the

<footnote n="1">¹Here we do not distinguish between major/minor mergers or continuous accretion, but see e.g., [65].</footnote>
conditional probability distribution

$$P_{\omega_f}(\omega_f) = \frac{1}{\sqrt{2\pi}} \int_0^1 \frac{M_o}{M(\tilde{S})} \left( \frac{\omega_f^2}{\tilde{S}^{5/2}} - \frac{1}{\tilde{S}^{3/2}} \right) e^{-\frac{\omega_f^2}{2\tilde{S}}} dS, \quad (2.15)$$

with $\tilde{S} = (S_f - S_o)S + S_o$.

Implementing this change of variable in eq. (2.12), the prediction for $\Delta b_\phi$ can be also expressed in terms of a derivative with respect to $\omega_f$:

$$\Delta b_{\phi} = 2 \frac{\partial \ln P_{z_f}(fM_o, z_f|M_o, z_o)}{\partial \ln \sigma_8} = -2 \left( 1 + \frac{d \ln P_{\omega_f}}{d \ln \omega_f} \right). \quad (2.16)$$

Consequently, if one considers the parameterization of $b_\phi$ through the variable $p$, as in eq. (2.11), then the expression for $p$ becomes

$$p = 1 + \frac{1}{\delta_c} \left( 1 + \frac{d \ln P_{\omega_f}}{d \ln \omega_f} \right). \quad (2.17)$$

Therefore, the correction to the universality relation, computed by including the effect of assembly bias, depends on redshifts ($z_f$ and $z_o$) and mass ($M_o$ and $fM_o$) only through the single variable $\omega_f$.

The approach considered here and the analytical expression of $P_{\omega_f}$ in eq. (2.15) provide a remarkably good description at high halo masses ($M_o \gtrsim 5 \times 10^{12} M_\odot$), as already tested in previous works [39, 42]. This result is further corroborated by comparisons with different N-body simulations below. However, for lower halo masses, the analytical prediction breaks down, due to $P_{\omega_f}$ taking a dependence on $M_o$ and $z_o$.

We highlight here that the form of Equations (2.16) and (2.17) does not depend on the assumptions made about the analytical form of $P_{\omega_f}$. In fact, they could directly follow from eq. (2.12) by adopting the change of variable defined by eq. (2.14). Therefore, a prediction for $p$, or equivalently for $\Delta b_\phi$, can be obtained from $P_{\omega_f}$ even if its expression is calibrated on dark matter simulations and not derived analytically. Building upon this argument, we provide below a simple modification of the ePS formula (which we refer to as e$^2$PS) for $P_{\omega_f}$, which fits also the simulated low mass halos well.

## 3 N-body simulations and comparison with ePS-derived expressions

In this work, we consider different sets of simulations with a double aim: on one hand, we make use of N-body simulations to evaluate the accuracy of the theoretical predictions described above; on the other hand, we study through hydrodynamical simulations how galaxy properties might be used as a proxy for halo assembly bias. To this end, we consider both the halo catalogs and halo merger trees available in these sets, which allow us to study the $b_\phi$ dependence on the halo formation history.

One set is the Quijote-PNG [44] suite of gravity-only simulations, which includes primordial non-Gaussianity of the local type in the initial conditions. This suite of N-body simulations consists in 1000 different realizations, each of them having a different value of local $f_{\text{NL}} \in [-300, 300]$ and different seed of the initial conditions. The cosmological
<table><caption>**Table 1.** Summary of the cosmological parameters and properties of the simulations used in the analysis.</caption><thead><tr><th></th><th>Ω<sub>m</sub></th><th>Ω<sub>Λ</sub></th><th>Ω<sub>b</sub></th><th>σ<sub>8</sub></th><th>h</th><th>n<sub>s</sub></th><th>L<br>(Mpc/h)</th><th>m<sub>p</sub><br>(M<sub>☉</sub>/h)</th></tr></thead><tbody><tr><td>Quijote-PNG</td><td>0.3175</td><td>0.6825</td><td>0.049</td><td>0.834</td><td>0.6711</td><td>0.9624</td><td>1000</td><td>6.6 × 10<sup>11</sup></td></tr><tr><td>IllustrisTNG</td><td>0.3089</td><td>0.6911</td><td>0.049</td><td>0.816</td><td>0.6774</td><td>0.9667</td><td>205</td><td>7 × 10<sup>7</sup></td></tr></tbody></table>

parameters are fixed to match those used in the fiducial Gaussian suite Quijote [43], which we also use in our work. The size of the simulation box is $L = 1 h^{-1} \text{Gpc}$ and the particle mass $m_p = 6.6 \times 10^{11} h^{-1} M_{\odot}$.

Another set is the IllustrisTNG [47] suite of high-resolution magneto-hydrodynamical simulations with Gaussian initial conditions. These simulations use a model for galaxy formation and evolution, calibrated to reproduce a number of galaxy properties such as the star formation rate, photometry and metallicities. Specifically, we focus on the TNG300 simulation, which uses a box with size $L = 205 h^{-1} \text{Mpc}$ and has a mass resolution of $m_p = 7 \times 10^7 h^{-1} M_{\odot}$. In this section we make use only of the dark matter aspect of the IllustrisTNG, the baryonic (galaxy properties) part is considered in section 4.

Table 1 summarizes the main properties of the simulations, including their cosmological parameters. The mass resolution introduces a mass limit under which the halos are not well resolved and their formation redshift is not well defined. This is due to the fact that for small $M_o$, there may not exist any progenitor with mass $fM_o$ (from which $z_f$ is defined) across the merger tree. Therefore, in the following we will consider only halos above a certain mass threshold: $M_o > 10^{14} h^{-1} M_{\odot}$ for Quijote-PNG and $M_o > 5 \times 10^{11} h^{-1} M_{\odot}$ for IllustrisTNG. These thresholds guarantee that $z_f$ is well defined for all the halos at $z_o = 0$ and are sufficient for $z_o > 0$ halo samples. Due to the different mass resolution, the two simulations can be used in a complementary way to test the ePS predictions on different mass ranges.

In order to study the predictions on the halo merger history, we make use of Quijote-PNG (IllustrisTNG) halo catalogs constructed with the Rockstar (FoF) halo finder and merger trees generated with the ConsistentTrees (Sublink) algorithm. For each halo identified in the catalog at $z_o$, its mass $M(z)$ is tracked back along the main (most massive) progenitor branch. The mass history is then reconstructed by interpolating $M(z)$ across the available snapshots and the formation redshift $z_f$ is defined as the redshift at which the halo mass is $M(z_f) = fM_o = M_o/2$, with $M_o$ being the mass at $z_o$. Hence each halo identified at $z_o$ is associated a mass $M_o$ and a formation redshift $z_f$. Then, the variable $\omega_f$ is computed for each halo by using eq. (2.14).

Our choice of $f = 0.5$ follows the original prescription of the ePS presented in ref. [40]. Although different values, as for example $f = 0.75$, have been chosen in the literature to describe cluster abundances [66], ref. [39] showed that, in the context of the non-Gaussian halo assembly bias, the two choices give consistent results. More generally, the description of assembly bias may be improved by considering other, integrated, quantities beyond the formation redshift $z_f$, which is subject to the variability of the halo mass accretion history, and thus might be noisy. In this particular context, especially where we consider large halo
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Relation between $\omega_f$ and $z_f - z_o$ for IllustrisTNG ($\bar{M} = 2 \times 10^{12} M_{\odot} h^{-1}$) and Quijote ($\bar{M} = 2 \times 10^{14} M_{\odot} h^{-1}$) halos at $z_o = 1$. The lines and the bands represent the mean and the standard deviation $\omega_f$ computed in $z_f - z_o$ bins of full halo sample. We show both the approximated $\omega_f$, computed through eq. (3.1), and the actual one, computed using eq. (2.14). The values of $n$ are obtained by fitting the approximation against the actual values of $\omega_f$ in the two samples.</figcaption></figure>

samples, the noise gets effectively averaged out. The parameterization chosen here benefits from having a transparent relation with the analytical prediction given by the ePS. Hence the adoption of less noisy parameterizations is left to future work.

## 3.1 Describing halo formation history through $\omega_f$ and $P_{\omega_f}$

In the context of the ePS formalism, the variable $\omega_f$ was introduced in order to express the conditional mass function in terms of a single variable, as $P_{\omega_f}(\omega_f)$. Strictly speaking, this is only valid under the approximation of an Einstein de Sitter cosmology with a power law linear matter power spectrum, as the one adopted in [40, 41].

In fact, we see from eq. (2.15) that the ePS prediction also depends on the ratio $\frac{M_o}{M(S)}$ which introduces a mass dependence, although it is very mild. Only by knowing the relation $S(M)$ analytically, this ratio can be integrated out and $P_{\omega_f}$ depends on $\omega_f$ only. Specifically, refs. [40, 41] consider an EdS universe, where $D(z) = (1 + z)^{-1}$ and assume a power-law linear matter power spectrum $P_{\text{lin}}(k) \propto k^n$. In this case, the variance can be analytically computed from eq. (2.9) and written as $S(M) = \delta_c^2 \left(\frac{M_*}{M}\right)^{\frac{n+3}{3}}$, where $M_*$ is a characteristic mass such that $S(M_*) = \delta_c^2$.

In this approximation, we can better understand the relation between $\omega_f$ and $M_o, z_o, z_f$. By substituting these approximations in eq. (2.14), we obtain

$$ \omega_f = \left( f^{-\frac{n+3}{3}} - 1 \right)^{-1/2} \left( \frac{M_o}{M_*} \right)^{\frac{n+3}{6}} (z_f - z_o) \quad (3.1) $$

Therefore at fixed mass $M_o, \omega_f$ is a linear rescaling of $(z_f - z_o)$, with a slope dependent on the power spectrum exponent $n$. The power-law approximation is valid in a sufficiently small
<figure><img src="image_3.png" /><figcaption>**Figure 2.** $z_f - z_o$ as a function of halo mass for the full sample of IllustrisTNG halos at $z_o=1$. The halos are sorted by their $\omega_f$ and divided in subsamples. Each subsample, represented by a different colour, contains 10% of the full halo sample. The black crosses indicate the mean mass and mean $z_f - z_o$ of the subsamples. The mean value of $\omega_f$ in each subsample is reported in the legend.</figcaption></figure>

range of $k$, which translates into a range of halo masses through the top-hat filter smoothing. Therefore, the power-law index $n$ will vary according to the mean mass $\bar{M}_o$ of the halo sample considered. In figure 1 we compare this approximation with the original expression in eq. (2.14) for both IllustrisTNG ($\bar{M}_o = 2 \times 10^{12} M_{\odot} h^{-1}$) and (Gaussian) Quijote ($\bar{M}_o = 2 \times 10^{14} M_{\odot} h^{-1}$) halos at $z_o = 1$. For these two different mass ranges, the approximation works well for a power-law power spectrum with $n = -1.9$ and $n = -1$, respectively. The scatter around the mean linear relation is due to the fact that the halo mass $M_o$ varies inside the samples considered. The maximum $z_f - z_o$ in the Quijote sample is lower than the IllustrisTNG one, reflecting the fact that heavier halos tend to form later than lighter ones, as can be noticed by looking at figure 2. In what follows, unless specified, we will always use the eq. 2.14 to compute $\omega_f$, but this linear relation between $\omega_f$ and $z_f - z_o$ indicates that, in small mass bins, selecting a certain fraction of the oldest (youngest) halos is equivalent to selecting the same fraction of halos with largest (smallest) $\omega_f$, as assumed by ref. [39]. However, $z_f$- and $\omega_f$- subsamples selected from the full sample of halos do not exactly match. In figure 2 we illustrate how quantiles of halos selected by their $\omega_f$ value, each one containing 10% of the full halo sample, project onto the $(M_o, z_f - z_o)$ plane. Generally, the same $\omega_f$ quantile contains both low mass old halos and higher mass younger halos. This “tradeoff” can be understood by looking at the approximation (3.1) for $\omega_f$ fixed. Incidentally, due to this tradeoff, the different $\omega_f$ quantiles share roughly the same mean halo mass (differing at most by 3%).

For the approximation described above, and for $f = 0.5$, we can solve the integral in
eq. (2.15) analytically and obtain the following form for $P_{\omega_f}$:

$$P_{\omega_f} = 2\omega_f \text{erfc} \left( \frac{\omega_f}{\sqrt{2}} \right) \qquad (3.2)$$

where erfc($\omega_f$) = 1 - erf($\omega_f$) is the complementary error function. The function in eq. (3.2) only depends on $\omega_f$ and as such it should be universal, i.e., valid regardless the mass and the redshift of the halos considered. Although this is true at large halo masses [39], a comparison with the $P_{\omega_f}$ constructed from the simulations considered in the present work shows a residual mass and redshift dependence. In particular, as we show in figure 3, as the mean mass and the redshift of the halo sample decrease, the ePS prediction does no longer accurately describe the $P_{\omega_f}$ obtained from simulations.

However, the shape of the distribution is preserved, although broadened and with a mean shifted towards larger $\omega_f$ values (i.e., earlier formation redshift). In order to account for this dependence, we propose a simple modification of eq. (3.2), shifting the argument of the erfc by a constant $\omega_0$:

$$P_{\omega_f} = P_{\omega_f}(\omega_f|\omega_0) = N(\omega_0)\omega_f \text{erfc} \left( \frac{\omega_f - \omega_0}{\sqrt{2}} \right), \qquad (3.3)$$

where we fit for the value of $\omega_0$, dependent on the mass and redshift of the sample of halos considered. The factor $N(\omega_0)$ ensures that the probability distribution is correctly normalized (and is such that $N(\omega_0 = 0) = 2$). The functional form of eq. (3.3), which we will refer to as e$^2$PS for the rest of the paper, can be considered as a phenomenological model, or a fitting formula which includes the ePS prediction as a particular case with $\omega_0 = 0$. The e$^2$PS formula shows a remarkably good fit with simulations across different redshifts and mass ranges, as shown in figure 3.

The best-fit value of $\omega_0$ approaches $\omega_0 = 0$ as the mass and the redshift increase and correspondingly, the e$^2$PS starts to coincide with the ePS, as illustrated in figure 3.

In the right panel, we report the prediction for $\Delta b_\phi$ as computed from eq. (2.16) by taking the derivative of the conditional mass function. We show both the ePS prediction and the one obtained from simulations, by taking the numerical derivative of the histograms reported on the left. Moreover, we also report the prediction for $\Delta b_\phi$ as derived from the best-fit e$^2$PS curve, calibrated on the corresponding $P_{\omega_f}$. The prediction for large values of $\Delta b_\phi$ is less precise, as it is related to the tail of the distribution $P_{\omega_f}$, where the noise increases due to small number statistics.

At redshifts snapshots $z = 0$ and $z = 0.7$ there is no evidence for any significant bias induced in $b_\phi$ by the e$^2$PS predictions keeping in mind that extremely large $\omega_f$ and $\Delta b_\phi$ values apply only to extremely rare objects with large statistical errors. At $z = 1$ there may be a hint of a possible mis-match for $\Delta b_\phi > 5$, $\omega_f > 2$, but statistics are poor; we return to this point in figure 4.

The mass and redshift dependence of $P_{\omega_f}$ is a result of the breakdown of the universality of the mass function, i.e., the Press-Schechter solution (2.13). This effect has already been discussed in the literature [63, 67] and modifications of the assumptions made in the PS approach might explain this discrepancy. This is beyond the scope of this work, but some insights on the theoretical modeling can be gleaned by looking at the form of the
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison between results from N-body simulations and the predictions from ePS and e$^2$PS. The grid on the left shows the results for $P_{\omega_f}$. The mass $M_o$ of the halos considered increases from left to right, while $z_o$ increases from bottom to top. The first two columns show the distributions of IllustrisTNG halos, while the third one for Quijote. In each panel, the histograms represent the results from simulations, the grey curve corresponds to the ePS prediction, while the coloured one is the best-fit e$^2$PS curve. On the right, we show the $\Delta b_\phi$ computed from eq. (2.16). The solid gray (coloured) line correspond to the ePS (e$^2$PS) model, the dotted lines correspond to the simulations. The coloured band around the simulation curve quantifies the error propagated from the bin count errors in the corresponding histogram on the grid. For brevity, we only report the results related to the first column of the grid. The dashed black horizontal line corresponds to the universality relation prediction, $\Delta b_\phi = 0$.</figcaption></figure>

e$^2$PS formula (3.3). In particular, the presence of a correction $\omega_0 = \omega_0(M, z)$ suggests the modification $\partial_c(z) \rightarrow \partial_c(M, z)$ in eq. (2.14). Such a modification would require the determination of the upcrossing statistics for what is known in the literature as a *moving barrier* [55, 57, 59–61, 68].

## 3.2 Accuracy of $P_{\omega_f}$ modeling: results on $\Delta b_\phi$ from non-Gaussian simulations

In the previous section we have shown how the variable $\omega_f$ and its distribution $P_{\omega_f}$ can be used to study halo assembly bias, reporting results from Gaussian simulations. Once $P_{\omega_f}$ is known, the deviation from the universality relation, $\Delta b_\phi$, can be obtained from eq. (2.16). Beyond being measured from simulations, these quantities have been tested against the ePS prediction and a fitting function has been provided when this prediction breaks down. In this section we validate these results on the Quijote-PNG non-Gaussian simulations, fitting $\Delta b_\phi$ by measuring the scale dependent bias induced by local PNG.
Of the Quijote-PNG suite we consider 500 simulations with $|f_{NL}| > 150$. This selection serves to amplify the signal and obtain better constraints on $b_\phi$.

For each simulation, we study the two snapshots at $z_o = 0, 1$. After having identified each halo in the snapshot, its $z_f$ is computed by interpolating across the merger tree and the pair of variables ($M_o, z_f$), where $M_o$ is the halo mass, are transformed into the variable $\omega_f$ by using eq. (2.14).

Our aim is to divide the full halo sample in each snapshot in 10 $\omega_f$ bins so that each bin has the same number of halos, as those in figure 2. As it is evident from figure 2 even though these subsamples have similar mean masses, it is not guaranteed they would have matching mass functions. As we can see from eq. (2.12), this is required to isolate the universality relation contribution to $b_\phi$ and fit $\Delta b_\phi$, as also discussed in [39]. Therefore, following [39], we proceed as follows.

We first divide the full sample in 30 small mass bins, each one containing the same number of halos (and having larger width for larger masses, accordingly to the mass function). Each mass bin is then divided into 10 $\omega_f$ quantiles. Finally, halos belonging to the same $\omega_f$ quantile in different mass bins are stacked together. By construction, we end up having 10 samples with mass functions matching that of the full sample.

From these 10 samples, we construct 10%, 20%, ..., 90% quantiles of the oldest (large $\omega_f$) and youngest (small $\omega_f$) halos. We test the predictions discussed above by analyzing the clustering properties of these quantiles, following the same approach as [39]. For each quantile, we construct the dark matter density field $\delta_m$ from the dark matter simulations snapshots and the halo number overdensity $\delta_h$ from the corresponding halo catalog.

The scale-dependent bias affects the largest scales, where the number of power spectrum modes (when averaged over $k$ bins) is relatively low. Therefore, rather than computing a power spectrum, we choose to fit eq. (2.3) mode by mode, up to $|k| < 0.03$ Mpc $h^{-1}$ and assume to have Poisson shot noise, $\langle \epsilon(\mathbf{k})\epsilon^*(\mathbf{k}) \rangle = \bar{n}^{-1}$, with $\bar{n}$ mean number density of the sample.

The variance about the model is

$$ \left\langle \left( \delta_h^* \delta_m - \left( b_1 + b_\phi \frac{f_{\mathrm{NL}}}{\mathcal{M}(k, z_o)} \right) \delta_m^* \delta_m \right)^2 \right\rangle = \frac{P_{mm}(k)}{\bar{n}}, \quad (3.4) $$

where the average is computed over the Poisson noise, assumed to be uncorrelated with $\delta_m$. In order to perform the fit, we fit for the values of $b_1$ and $b_\phi$ that minimize the $\chi^2$, computed as follows:

$$ \chi^2 = \sum_{\mathbf{k}} \frac{\left( \mathrm{Re}[\delta_h^* \delta_m] - \left( b_1 + b_\phi \frac{f_{\mathrm{NL}}}{\mathcal{M}(k, z_o)} \right) \delta_m^* \delta_m \right)^2}{P_{mm}(k)/(2\bar{n})} \quad (3.5) $$

We only consider the real component of $\delta_h^* \delta_m$, hence the noise is a factor of 2 smaller than in eq. (3.4).

For each simulation with different $f_{NL}$ and for each $z_o$, the fit procedure is applied to the full sample and the 9 cumulative subsamples of both the oldest and youngest halos. Then, $\Delta b_\phi$ is estimated from the difference between the $b_\phi$ fitted from the full sample and the one fitted from the subsample:

$$ \Delta b_{\phi} = b_{\phi}^{\mathrm{full}} - b_{\phi}^{\mathrm{sub}}. \quad (3.6) $$
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Average $\Delta b_{\phi}$ measured from the cumulative subsamples of the youngest and oldest halos. ePS predictions for both fractions are also shown. On the left, we show results for the population of halos at $z_o = 0$. On the right, for $z_o = 1$. The average is computed over 500 Quijote-PNG simulations with different $f_{\text{NL}}$ (and the errors denote the error on the mean). The halo mass threshold chosen to avoid unresolved halos is $M_o > 10^{14} M_{\odot}/h$.</figcaption></figure>

Finally, we take the average of $\Delta b_{\phi}$ over the different simulations, at same redshift $z_o$.

In figure 4 we report the results obtained at $z_o = 0$ and $z_o = 1$ for $\langle \Delta b_{\phi} \rangle$, averaged over the 500 realizations with different $f_{\text{NL}}$. Each of the quantile we consider in the analysis extends across a range of $\omega_f$ values, [$\omega_f^{\text{min}}$, $\omega_f^{\text{max}}$]. The theoretical prediction for the $\Delta b_{\phi}$ of each quantile is then found by averaging eq. (2.16),

$$ \langle \Delta b_{\phi} \rangle = -2 \left\langle 1 + \frac{d \ln P_{\omega_f}}{d \ln \omega_f} \right\rangle, \qquad (3.7) $$

where the average is computed over the range [$\omega_f^{\text{min}}$, $\omega_f^{\text{max}}$], using the probability distribution $P_{\omega_f}$. In this case we use the ePS prediction for $P_{\omega_f}$, which, as discussed in the previous section, which is shown to be accurate for the range of halo masses of Quijote-PNG. The main result is that old halos have a larger $b_{\phi}$, while young ones have a lower $b_{\phi}$ with respect to that of the full sample. The effect is asymmetric between the two, as $\Delta b_{\phi}$ increases rapidly at large $\omega_f$, while gradually at small $\omega_f$, as illustrated in figure 3. This asymmetry is in accordance with similar results which parameterized assembly bias through the halo concentration instead of the halo formation redshift [28–30].

If we consider the parameterization in terms of $p$, the 10% of the oldest halos have a $\Delta b_{\phi}$ which corresponds to $p \simeq -1.4$, while the 10% of youngest ones, $p \simeq 1.6$, as already found in [21] for recent mergers.

The ePS prediction shows a remarkable accuracy across the two different redshifts $z_o = 0, 1$, especially for the oldest halos; there are some small discrepancies for small $\lesssim 0.2$ fractions of the youngest halos but the agreement is good for fractions $\gtrsim 0.3$. The errorbars are larger for $z_o = 1$ as a result of the increased shot noise contribution with respect to $z_o = 0$. These results are consistent with those of [39], where they used different non-Gaussian simulations and range of masses.
In summary, as figure 4 shows, there is no evidence of significant systematic effect on $\Delta b_{\phi}$ introduced by the $e^2$PS modeling above the statistical errors corresponding to a volume of 500 (Gpc/h)$^3$.

# 4 Galaxy proxy for halo assembly bias: results from IllustrisTNG simulations

So far we have presented how to describe the impact of assembly bias on the local PNG bias parameter $b_{\phi}$, focusing on simulated halo catalogs and merger trees. By knowing ($M_o, z_f$) of a sample of halos, and computing $\omega_f$, we have shown how to obtain a prediction for $b_{\phi}$. However, halo properties such as $M_o, z_f$ are not directly observable and their inference is tied to the galaxy-halo connection, which description includes modeling halo occupation distribution (HOD), halo mass accretion history, as well as galaxy formation (see ref. [69] for a detailed review). Quantities which can actually be observed in LSS surveys are only (some) galaxy properties. Therefore, our aim is to find a galaxy-related quantity which can be connected to the host halo formation redshift $z_f$, or even better, to the $\omega_f$ of the corresponding host halo. In other terms, the goal is to find an observable proxy of halo assembly bias.

In the literature, different galaxy proxies have been used, most of which are also proxies for the age of the galaxy, which is expected to be statistically correlated to the age of the halo. Examples of properties used as proxies are stellar mass $M_*$ [70], specific star formation rate (sSFR) [71, 72] and the galaxy color $g-r$ [73]. However, some of these properties are to some extent model-dependent: for example $M_*$ and sSFR of galaxies are indirect measurements inferred from spectroscopic and photometric data [74, 75].

Our approach here is to remain as close to observations as possible, focusing on assembly bias proxies constructed from (combination of) photometric bands beyond the optical ones $g, r$. While the proxy could certainly be improved by adding extra information e.g., a detailed spectroscopic analysis, photometry is more readily available and we focus on this here.

## 4.1 Correlations between galaxy photometry and halo formation redshift

We consider the sample of IllustrisTNG halos at $z_o = 1$ used for the results in section 3. For each halo in the sample we select its central galaxy, identified by the largest subhalo within the parent halo. All the galaxy properties, including $M_*$, SFR and photometric quantities are taken directly from the Subfind subhalo catalogs released by IllustrisTNG. We only consider star-forming (SFR> 0) galaxies with stellar mass $M_* > 10^8 M_{\odot}h^{-1}$.

The full sample contains galaxies with different properties, spanning over a large range of halo and stellar masses. Real galaxy samples used for observational purposes, instead, always rely on some selection criteria such as magnitude cuts, which reduce the range of variation of some properties. Therefore, besides the full sample of galaxies we also consider 8 subsamples, through which we can assess how the optimal proxy depends on the sample selection criterion. These 8 subsamples are the following ones:

* 3 subsamples selected considering different halo mass ranges $M_h \in [5 \times 10^{11}, 2 \times 10^{12}], [2 \times 10^{12}, 10^{13}], [10^{13}, 10^{14}] M_{\odot}h^{-1}$ which we refer to as $M_h^1, M_h^2, M_h^3$ respectively;
* 3 subsamples obtained with stellar mass cuts $M_* \in [10^8, 2 \times 10^9], [2 \times 10^9, 3 \times 10^{10}], [3 \times 10^{10}, 10^{12}]$ $M_{\odot}h^{-1}$ which we refer to as $M_s^1, M_s^2, M_s^3$ respectively;

* 2 subsamples identified by specific star formation rate cuts, $\log[sSFR] > -9.23$ and $\log[sSFR] < -9.23$, respectively mimicking the populations of emission-line galaxies (ELG) and luminous red galaxies (LRG).

We only consider the heaviest galaxies within the ELG- and LRG- like subsamples, as to match the DESI target number densities for ELG, $n = 7 \times 10^{-4}\text{Mpc}^{-3}h^3$ and LRG, $n = 2 \times 10^{-4}\text{Mpc}^{-3}h^3$ [76, 77]. These cuts follow the sample definitions adopted by [30], based on previous studies about the link between IllustrisTNG and DESI target galaxies [78, 79]. For each galaxy in the sample considered, we focus on the set of 8 photometric bands $X = \{U, B, V, K, g, r, i, z\}$, with AB magnitudes based on the summed-up luminosities of all the stars within the galaxy. Taking the differences between these magnitudes, we form a set of 28 colors $C = X_i - X_j$, with $X_{i,j} \in X$.

As a preliminary step, figure 5 illustrates the correlations between these galaxy colors and the $\omega_f$ of their host halo, for both the full sample and the 8 subsamples defined above.

The correlation between a given galaxy color and $\omega_f$ depends on the subsample, but some common trends can be identified. In particular, colors involving the U band appear to be the least correlated with $\omega_f$, while colors involving the K band show the largest correlations. This can be understood by taking into account that this band, centered in the infrared, can be used to track the age of a galaxy [10, 80], which is in turn related to the formation time of the host halo.

Galaxy colors in the halo mass-selected subsamples manifest larger correlations, reaching values beyond $r = 0.6$, so such samples would be ideal targets to estimate $\omega_f$ from photometric quantities. However, their construction is based on the knowledge of the halo mass, which needs to be evaluated indirectly, as for example by measuring the linear bias $b_1$ from the power spectrum [81, 82].

Subsamples more directly related to single-object observational properties are instead the stellar mass selected $M_s^i$ and the ELG/LRG-like ones. Galaxy colors in the $M_s^i$ samples roughly share the same behavior, although we can observe that for $M_s^1$, the correlation varies on a wider range, compared to higher $M_s$ subsamples as $M_s^3$. In other words, for low stellar mass galaxies, we can clearly distinguish between colors that are notably good proxies for $\omega_f$ and colors that do not give much information about the halo formation time. On the other hand, for more massive galaxies there are not sharp differences among the different colors. This features can be also identified within the $M_h^i$ samples.

As it is evident, the LRG-like sample shares common correlation properties with the $M_h^3$ sample, since LRGs are tipically hosted by massive halos. The colors of ELGs, on the other hand do not appear to track the halo formation time.

## 4.2 Optimal proxy for $\omega_f$ from a combination of galaxy colors

In order to select the optimal galaxy proxy of halo assembly bias, we choose a simple *ansatz* and assume that $\omega_f$ can be recovered from a linear combination of 3 galaxy colors $C_i$ for $i = 1, 2, 3$, $Y = a_1C_1 + a_2C_2 + a_3C_3$. We have verified that including more than 3 colors does
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Correlation between the 28 galaxy colors and $\omega_f$. Each pixel $(i, j)$ identifies a galaxy color $C = X_i - X_j$ and is color-coded depending on the correlation $r(C, \omega_f)$, ranging from dark purple for low $r$ to light yellow for high $r$. The different panels illustrate the results for different subsamples, as defined in the text.</figcaption></figure>

not improve the results further, as discussed in the appendix A. For each possible choice of $C_1, C_2, C_3$, the coefficients $a_1, a_2, a_3$ are optimized as to maximize the Pearson correlation $r(Y, \omega_f)$. Among the 3276 possible combinations of 3 colors that can be extracted from the full set of 28 colors, the one which gives the highest $r$ is selected as optimal proxy. This procedure is repeated for both the full sample and the 8 subsamples mentioned above, thus we obtain a different optimal proxy $Y$ for each sample.

Once the optimal proxy has been determined for each sample, we test the proxy-recovered $\Delta b_{\phi}^{(r)}$ against the true $\Delta b_{\phi}^{(t)}$ (which assumes an exact knowledge of $\omega_f$, through $M_o, z_f$ of the host halo). To this end, for each sample we construct $P_{\omega_f}$, fit for $\omega_0$ in eq. (3.3) and
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Recovered $\Delta b_{\phi}^{(r)}$ of subsamples of galaxies selected by using the optimal galaxy proxy $Y$, as a function of their true $\Delta b_{\phi}^{(t)}$. Each point identifies a subsample, representing the 10% quantile of the sample of galaxies (IllustrisTNG galaxies at $z_o = 1$). For both the recovered and the true ones, the mean value and the 16th-84th percentiles of the $\Delta b_{\phi}$ distribution within each quantile are reported. In each panel, the correlation $r(Y, \omega_f)$ as well as the RMSE of the recovered vs true $\Delta b_{\phi}$ are shown, together with the identity relation indicated by the orange line. The different panels illustrate the results for different subsamples, as defined in the text.</figcaption></figure>

then use eq. (2.16) to predict $\Delta b_{\phi}$ for each galaxy. Then, on one hand galaxies are ranked ordered according to their $\omega_f$; on the other hand, according to their proxy $Y$. These two samples, containing the same galaxies with different ordering, are divided into 10% quantiles. In the ideal case of maximum correlation $r(Y, \omega_f) = 1$, the orderings of the two samples coincide, while in general there is a mismatch.

Consequently, within each Y-ranked quantile the distribution of $\Delta b_{\phi}^{(r)}$ will differ from that of $\Delta b_{\phi}^{(t)}$, related to the corresponding $\omega_f$-ranked quantile. We evaluate this mismatch by considering the root mean square error RMSE = $\sqrt{\frac{\sum_{n}^{N} (\Delta b_{\phi}^{(r)} - \Delta b_{\phi}^{(t)})^2}{N}}$, with N the size of the galaxy sample, that quantifies the scatter around the identity relation.

In figure 6 we report the mean value and 16th-84th percentiles of the $\Delta b_{\phi}^{(r)}$ and $\Delta b_{\phi}^{(t)}$ distributions within the 10 quantiles, one as a function of the other, for all the samples considered in the analysis.

Combining galaxy colors together allows one to get larger correlations with $\omega_f$, up to $r = 0.65$. The different “performance” of the various samples (i.e. the maximum correlation
between galaxy colors and $\omega_f$) mirrors the results in figure 5, with the halo mass selected samples manifesting largest correlations and ELG-like one, the smallest ones. Generally, larger $r$ does not guarantee to have a lower RMSE. In other terms, a good prediction of $\omega_f$ through the proxy $Y$ does not ensure to recover $\Delta b_\phi$ optimally. This is due to the properties of the $\omega_f$ distribution within each sample and the nonlinear relation between $\omega_f$ and $\Delta b_\phi$.

More specifically, as one can see in figure 3, $\Delta b_\phi$ varies rapidly for large $\omega_f$, while it stays roughly constant for small $\omega_f$. This propagates into a larger dispersion in both $\Delta b_\phi^{(t)}$ and $\Delta b_\phi^{(r)}$ within the high $\omega_f$ quantiles. As a consequence, samples containing old halos are affected, as the LRG-like one, which has a relatively large RMSE, although $r = 0.58$.

In general, although most of the proxies have good correlations with $\omega_f$, our simple model is not flexible enough to fully predict $\omega_f$ and recover $\Delta b_\phi^{(t)}$ completely. Specifically, for low (high) $\omega_f$ quantiles, $\Delta b_\phi^{(r)}$ is always over-(under-)estimated with respect to $\Delta b_\phi^{(t)}$. We will return to this systematic shift in section 6.2. On the one hand, predictions may be improved by considering nonlinear models or using machine learning techniques, possibly adding further information beyond galaxy colors; this is left for future work. On the other hand, simple models have the advantage of being more robust to changes in simulations settings and parameters, as will be shown in section 5.

# 5 Marginalizing over cosmology and astrophysics: tests on CAMELS-TNG

Galaxy properties are determined by complex physical aspects, often referred to as baryonic effects, which also affect the underlying dark matter distribution through feedback processes. These include feedback from two main sources: active galactic nuclei (AGNs) and starbursts, the latter determined by supernova explosions and stellar winds ejected from newborn stars.

Describing these processes in cosmological simulations is very challenging: the scales involved are beyond the resolution of large cosmological boxes and phenomenological models are usually adopted, with parameters chosen as to match observations. The details of their implementation in IllustrisTNG can be found in [83]. Specific simulations make specific implementation choices, but there is an intrinsic theoretical uncertainty associated. In section 4 the assembly bias proxy has been calibrated for the specific implementation of IllustrisTNG; calibration on a different phenomenological prescription might result in a different proxy.

In order to account for the theoretical uncertainty associated to these phenomenological models, an approach is to marginalize over the parameters which define them. This is one of the purposes of the suite of CAMELS simulations [49], which scan different subgrid models for galaxy formation and various choices of astrophysical and cosmological parameters. Among the different models considered in CAMELS, we focus here on the subset of CAMELS-TNG simulations, run with the same subgrid prescription of IllustrisTNG. Although tests on other subsets as CAMELS-SIMBA and CAMELS-ASTRID [84] would provide a broader comparison among different models, galaxy photometry is not yet available for these sets.

The set of CAMELS-TNG simulations includes 1000 different simulations run with different initial random seed and cosmological and astrophysical parameters ranging in a Latin Hypercube with $\Omega_m \in [0.1, 0.5]$, $\sigma_8 \in [0.6, 1.0]$, $A_{\text{SN1}} \in [0.25, 4.0]$, $A_{\text{SN2}} \in [0.5, 2.0]$,
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Recovered $\Delta b_{\phi}^{(r)}$ of subsamples of galaxies selected by using the optimal galaxy proxy $Y$, as a function of their true $\Delta b_{\phi}^{(t)}$, similarly to figure 6. In the top panels, we report the results of the calibration of the proxy on CAMELS. Both the full sample and samples where either an $\Omega_m$ prior or a prior on astrophysical parameters are considered. In the bottom panels, the CAMELS-calibrated proxy is applied to IllustrisTNG full sample, in comparison with the first panel of figure 6.</figcaption></figure>

$A_{\text{AGN1}} \in [0.25, 4.0]$, $A_{\text{AGN2}} \in [0.5, 2.0]$. The 4 astrophysical parameters regulate starbursts and AGN feedback and are set to 1 in IllustrisTNG. CAMELS implements a resolution comparable with that of the original TNG300, while simulating a smaller box with $L = 25$ Mpc $h^{-1}$. We refer the reader to [49] for further details on the simulations.

Our goal is to use CAMELS-TNG galaxies to assess the robustness of the proxy calibration of section 4, which is tuned to the specific choice of astrophysical (and cosmological) parameters in IllustrisTNG. We analyse the snapshots at $z_0 = 0.95$ of each simulation. For each simulation, we consider the photometric data of the central galaxies and compute $\omega_f$ of their host halo, which depends on cosmology via $D(z)$ and $\sigma(M)$ appearing in eq. (2.14). All the galaxies belonging to different simulations (i.e., formed under different astrophysical models depending on the values of the baryonic feedback parameters) are put together in what we refer to as “full sample”. By analysing this sample we are effectively marginalizing over cosmology, (IllustrisTNG) astrophysics and including the effect of cosmic variance.

The different simulation settings in CAMELS give rise to differences with the IllustrisTNG galaxy sample used in section 4. In particular, the smaller size of the simulation box has an impact on the most massive galaxies which can be simulated, making the selection of an LRG-like sample impossible. Therefore we restrict our analysis to the full sample.

By using the same method described in section 4, we calibrate the optimal proxy on CAMELS-TNG full sample and then test it on IllustrisTNG full sample. The results are reported in the left column of figure 7, where the top panel is related to the calibration on CAMELS-TNG sample and the bottom panel to the test on IllustrisTNG sample.

Not unexpectedly, the proxy performance in the calibration sample degrades, effect that is
driven in part by the variation in cosmologies and in part by the variation of the astrophysical parameters. In particular, marginalizing over $\Omega_m$ has a large impact on $P_{\omega_f}$. Indeed, at fixed formation redshift $z_f$, extreme values of this parameter can produce very different values of $\omega_f$ (see eq. (2.14)), broadening the distribution $P_{\omega_f}$. This larger scatter in $\omega_f$ is propagated to $\Delta b_\phi$ differently for low and high $\omega_f$ values, as we can notice from figure 3. In particular, as $\Delta b_\phi$ grows rapidly for older halos, the $\Delta b_\phi^{(t)}$ distribution of the last 10% quantile gets skewed towards larger values with respect to the IllustrisTNG counterpart, illustrated in the bottom left panel. The distribution of the recovered $\Delta b_\phi^{(r)}$ is also affected by the marginalization over the astrophysical parameters, which mixes the relations between $\omega_f$ and the galaxy colors.

However, the test on IllustrisTNG of the CAMELS-calibrated proxy shows good results if compared with the IllustrisTNG-calibrated one (illustrated in the first panel of figure 6). The correlation with $\omega_f$ only lowers from $r = 0.55$ to $r = 0.49$ and the scatter around equality is practically unchanged, from RMSE = 3.75 to RMSE = 3.77. Therefore, despite the large range of variation of the CAMELS galaxy features, due to the marginalization over cosmological and astrophysical parameters, the performance of a CAMELS-calibrated proxy is marginally affected.

To partly disentangle the two effects (marginalization over cosmology vs marginalization over astrophysics) we repeat the analysis by imposing either a cosmology prior or priors on the astrophysical parameters.

The cosmology prior is imposed by considering the subset of 500 simulations that exclude extreme values of $\Omega_m$, which has been tested to have a larger impact on $P_{\omega_f}$ with respect to $\sigma_8$. Specifically, our prior reduces its range of variation to $\Omega_m \in [0.2, 0.4]$, symmetrically around the fiducial value of 0.3. The range of variation of the astrophysical parameters is instead left unchanged. The results about the calibration and the test are shown in the middle column of figure 7.

The astrophysical parameters prior is imposed by considering the subset of 500 simulations that are closer to the IllustrisTNG fiducial values. The ranges of variation are $A_{\text{SN1}} \in [0.3, 3.0]$, $A_{\text{SN2}} \in [0.58, 1.7]$, $A_{\text{AGN1}} \in [0.3, 3.0]$, $A_{\text{AGN2}} \in [0.58, 1.7]$ and the results are in figure 7, on the right.

Within the calibration sample, imposing a cosmology prior improves the correlation coefficient and reduces the scatter, while imposing an astrophysical parameters prior does not. As regards the test on IllustrisTNG, both proxies have a performance roughly as good as the one calibrated on the full sample.

In table 2 we summarize the combinations which define the optimal proxies, as obtained from the analyses in sections 4, 5.

# 6 Forecasts on $f_{\text{NL}}$ measurement

The observational proxies calibrated in sections 4 and 5 provide an inference for $b_\phi$ for each galaxy sample, which can be used as a prior when inferring constraints on $f_{\text{NL}}$. Although a prior knowledge of $b_\phi$ is not needed to detect local PNG, it is required to constrain the value of $f_{\text{NL}}$. The precision and accuracy of the $b_\phi$ prior will affect the error of the $f_{\text{NL}}$ measurement in two ways: statistical, $\sigma_{f_{\text{NL}}}$, and systematic $\sigma_{f_{\text{NL}}^{\text{(sys.)}}}$.
<table><caption>**Table 2.** Summary of the combinations of galaxy colors obtained for the optimal proxies calibrated on the respective sample, reported in the first column. The optimal proxy is reported in the second column and the third column reports the Pearson correlation coefficient between the proxy and the physical variable $\omega_f$ (see text for more details).</caption><thead><tr><th>calibration sample</th><th>proxy Y</th><th>r(Y, ω<sub>f</sub>)</th></tr></thead><tbody><tr><td>LRG-like</td><td>6.3(V − K) + 3.1(B − r) + 0.005(i − U)</td><td>0.58</td></tr><tr><td>ELG-like</td><td>4(V − r) + 2(i − g) + 0.8(g − z)</td><td>0.42</td></tr><tr><td>IllustrisTNG full</td><td>3.8(V − K) + 1.8(B − z) + 0.01(K − U)</td><td>0.55</td></tr><tr><td>CAMELS full</td><td>12.4(V − i) + 1.9(K − B) + (r − U)</td><td>0.41</td></tr><tr><td>M<sub>h</sub><sup>1</sup></td><td>4(i − r) + 2.5(V − z) + 0.2(z − U)</td><td>0.62</td></tr><tr><td>M<sub>h</sub><sup>2</sup></td><td>2(V − g) + (B − r) + 0.02(V − U)</td><td>0.65</td></tr><tr><td>M<sub>h</sub><sup>3</sup></td><td>2(V − g) + (B − r) + 0.014(z − U)</td><td>0.61</td></tr><tr><td>M<sub>s</sub><sup>1</sup></td><td>8.25(g − B) + 1.02(i − U) + (K − r)</td><td>0.51</td></tr><tr><td>M<sub>s</sub><sup>2</sup></td><td>10(V − r) + 3(V − B) + 0.2(K − z)</td><td>0.51</td></tr><tr><td>M<sub>s</sub><sup>3</sup></td><td>1.96(V − B) + 1.73(B − g) + (K − r)</td><td>0.54</td></tr></tbody></table>

depends on the bias $b_\phi$ of the sample because the noise in the data is fixed effectively, for small non-Gaussianity, by the power spectrum amplitude, the volume surveyed and the shot noise but the signal itself for a given $f_{\text{NL}}$ is proportional to $b_\phi$. Thus samples with larger $|b_\phi|$ yield better signal to noise in the data and thus smaller errors on $f_{\text{NL}}$. A sample that spans a range of $\omega_f$ will tend to have reduced $|b_\phi|$ compared to a sample with well selected $\omega_f$.

On the other hand, a systematic error in the estimate of $b_\phi$ of $\sigma_{b_\phi}^{(\text{sys.})}$ will introduce a systematic error in $f_{\text{NL}}$ of $\sigma_{f_{\text{NL}}}^{(\text{sys.})} = \sigma_{b_\phi}^{(\text{sys.})}b_\phi^{-1}f_{\text{NL}}$, biasing the $f_{\text{NL}}$ constraint. In what follows, we report the results of the analysis of these two different contributions to the forecasted error on $f_{\text{NL}}$.

The statistical error can be estimated from the Fisher information matrix as $\sigma_{f_{\text{NL}}} = 1/\sqrt{\mathcal{F}_{f_{\text{NL}}}}$, where $\mathcal{F}_{f_{\text{NL}}}$ reads

$$ \mathcal{F}_{f_{\text{NL}}} = \int_{k_{\text{min}}}^{k_{\text{max}}} dk \frac{\partial \mathbf{D}^T(k)}{\partial f_{\text{NL}}} \mathbf{Cov}^{-1}(k) \frac{\partial \mathbf{D}(k)}{\partial f_{\text{NL}}}, \quad (6.1) $$

with **D** being the data vector (i.e. the power spectrum) and **Cov** the covariance matrix of the data which we assume to be diagonal. The reason is that the scale-dependent bias plays a role on large scales, where non-linearities sourced by gravitational evolution are negligible and different Fourier modes evolve independently. Moreover we do not explore here possible degeneracies between $f_{\text{NL}}$ and other cosmological (or astrophysical) parameters.

In order to provide somewhat realistic forecasts, in this section we will consider the IllustrisTNG full sample, as well as the ELG- and LRG-like subsamples, as defined in section 4. For completeness, we also compare the results of the CAMELS- and the IllustrisTNG-calibrated proxies applied to the full sample. We provide both single tracer and multitracer results and refer the reader to [11, 33, 85] for the details on the derivation of the equations we report in what follows.
## 6.1 Single tracer: statistical error

Let us first consider a data vector $D(k) = P_{gg}(k)$ consisting in the power spectrum measurements of a single sample of galaxies at redshift $z_o$ with number density $\bar{n}$, where $P_{gg}(k)$ is given by

$$P_{gg}(k) = P_{\text{th}}(k) + P_{\epsilon} \equiv [b_1 + b_{\phi} f_{\text{NL}} \mathcal{M}^{-1}(k, z_o)]^2 P_{\text{lin}}(k) + \frac{1}{\bar{n}} \quad (6.2)$$

and we have assumed Poisson shot noise. By using eq. 6.1, the $f_{\text{NL}}$ error Fisher forecast becomes [29, 85]

$$\sigma_{f_{\text{NL}}}(k) = \sqrt{\frac{2(2\pi)^3}{V_S V_k} \frac{b_1 + b_{\phi} f_{\text{NL}} \mathcal{M}^{-1}(k, z_o)}{2b_{\phi} \mathcal{M}^{-1}(k, z_o)} \left(1 + \frac{1}{\bar{n} P_{\text{th}}(k)}\right)}, \quad (6.3)$$

where $V_S$ is the volume of the survey and $V_k = 4\pi k^2 \Delta k$ is the volume in Fourier space of the power spectrum shell at wavenumber $k$, with wavenumber bin size $\Delta k$. In order to provide quantitative results, we consider a realistic case with galaxy samples to have the same empirically calibrated DESI linear bias [4], $b_1 = 1.38, 1.79$ respectively for ELG- and LRG-like samples, computed at $z_o = 1$. The nominal number densities of the two samples are $\bar{n} = 7 \times 10^{-4} \text{Mpc}^{-3}h^3$ for ELG and $\bar{n} = 2 \times 10^{-4} \text{Mpc}^{-3}h^3$ for LRG. For the full sample, we use the same linear bias and number density as the ELG-like sample. We adopt the DESI $k_{\text{min}} = 0.0023 h \text{ Mpc}^{-1}$ as reported by [30] and the corresponding survey volume is assumed to be $V_S = (2\pi/k_{\text{min}})^3 = 20.4 \text{ Gpc}^3/h^3$.

From eq. (6.3), we can see that for a given number density $\bar{n}$, $\sigma_{f_{\text{NL}}}$ decreases with increasing $b_{\phi}$ due to both the denominator and the $P_{\text{th}}(k)$ in the bracket, which grows with $b_{\phi}^2$. By selecting a sample of galaxies with large $\omega_f$, hosted by the oldest halos, we can maximize $b_{\phi}$ and obtain the most precise constraints on $f_{\text{NL}}$. The smaller the sample, the larger the $b_{\phi}$. However, for small samples the number density $\bar{n}$ effectively introduces a tradeoff, due to the enhanced shot noise, which worsens the precision, increasing $\sigma_{f_{\text{NL}}}$. This is illustrated in figure 8, where we report $\sigma_{f_{\text{NL}}}$ as a function of the cumulative quantile of oldest halos, for both the ideal (where we exactly know $\omega_f$ for each halo) and the proxy selected samples. Selecting a subsample of galaxies with a large value of $b_{\phi}$ leads to better constraints than those obtained from the full sample for a volume of $20.4 \text{ Gpc}^3 h^3$. The improvement, in the ideal case ($\omega_f$ perfectly known), becomes almost a factor of 2 for LRG and a factor of 3 for ELG. This difference between the two is due to the fact that the $\omega_f$ distribution in the ELG sample is broader, so there is more room for improvement than in the LRG case. The same conclusion applies to the full sample, in which case the improvement is almost a factor of 3.

In the realistic case, in which 20-30% of galaxies most-likely hosted by older halos is selected by using a proxy as discussed above, one can enhance the precision on the $f_{\text{NL}}$ constraint by 20-25%. For the ELG-like sample, this is only a fraction of the potential improvement, due to the fact that the proxy is not optimal, as discussed in section 4. Therefore, within the context of a single-tracer analysis, the LRG-like sample appears to be a better target to provide $f_{\text{NL}}$ constraints, with the minimum $\sigma_{f_{\text{NL}}}$ being almost half of the ELG-like one. Certainly, an improved calibration of the proxy in the ELG-like sample would produce better results, but we leave this to future work. As regards the full sample,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Forecasted single-tracer error $\sigma_{f_{\text{NL}}}$ for the ELG-like, LRG-like subsamples and the full sample. In all panels the cumulative quantiles (on the x axis) are obtained with the true $\omega_f$ for the ideal case and the proxy-selected $\omega_f$ for the proxy case. Only cumulative quantiles of the oldest (largest $\omega_f$) halos are considered. For larger quantile fractions, increasingly younger halos are included in the sample. The green line in the right panel is obtained using the CAMELS-calibrated proxy applied to the IllustrisTNG galaxies.</figcaption></figure>

we observe that the CAMELS- and the IllustrisTNG-calibrated proxies have practically the same performance, as expected from the discussion in section 5. The improvement in this case is somewhat in-between the LRG- and the ELG-like sample related ones.

The upturn at low quantiles in figure 8 is related to the ratio between the shot noise $1/\bar{n}$ and $P_{th}$, contributing to $\sigma_{f_{\text{NL}}}$ in eq. (6.3). If the subsample selected is too small, i.e. $\lesssim 20\%$, the larger shot noise dominates the error $\sigma_{f_{\text{NL}}}$, which increases. This effect does not show up in the ideal case, where maximizing $b_\phi$ determines a larger $P_{th}$ which compensates the shot noise contribution.

As figure 8 shows, the use of the proxy only degrades the $f_{\text{NL}}$ statistical constraints by a factor 1.5 or less for LRG, by a factor 3 or less for ELG, and by a factor of 2 or less for the full sample, compared to the ideal case where the galaxy sample can be selected by $\omega_f$ perfectly.

## 6.2 Single-tracer systematic error

As figures 6 and 7 show, there is a residual systematic trend in $\Delta b_\phi^{(r)}$ with respect to $\Delta b_\phi^{(t)}$. These residuals, if uncorrected, would be responsible for introducing a systematic shift in the inferred $f_{\text{NL}}$.

In the top left panel of figure 9, we report the residuals as a function of the true $\Delta b_\phi$, for the relevant samples considered in this section. As shown, the recovered $\Delta b_\phi^{(r)}$ of younger halos is overestimated, while that of older halos is underestimated. In the top right panel, the residuals are shown as a function of the cumulative proxy-selected $\omega_f$ quantiles of older halos. The behaviour of the residuals is strikingly similar for all the cases considered: ELG-like, LRG-like and full samples, also when for the latter one the proxy is calibrated on CAMELS instead of IllustrisTNG. This indicates that the residuals could be modeled as to reduce the systematic error on $f_{\text{NL}}$ quite significantly. As an initial estimate of the maximum systematic contribution we assume it will not be corrected at all, and propagate the full extent of $\Delta b_\phi^{(r)} - \Delta b_\phi^{(t)}$. This is therefore a conservative estimate.

In the bottom left panel of figure 9 we report the relative (conservative) systematic error $\sigma_{f_{\text{NL}}^{(\text{sys.})}} f_{\text{NL}}^{-1}$, varying with the cumulative proxy-selected quantile of old halos, for all the samples considered in this section.
statistical error to be dominant for small $f_{NL}$, while the systematic one to dominate for large $f_{NL}$. This is illustrated in the bottom right panel of figure 9, where we show the systematic to statistical relative error for the LRG-like sample, for several representative values of $f_{NL}$. For the quantile where the statistical error is the smallest ($\simeq 20-30\%$), the systematic error is below the statistical one as long as $f_{NL} \lesssim 10$.

## 6.3 Multi-tracer

By selecting subsamples with different clustering properties from the full sample of galaxies, one can employ the multi-tracer approach [36] and enhance the precision on $f_{NL}$ constraints. The improvement is due to the suppression of cosmic variance, which has a strong contribution on large scales, where we have most of the local PNG signal. In this work we consider a 2-tracer approach (labeled hereafter A and B), in which case the data vector is composed by the autopower spectra of the two tracers and the cross power spectrum between them: $\mathbf{D}(k) = \{P_{gg}^A(k), P_{gg}^{A \times B}(k), P_{gg}^B(k)\}$. By using eq. (6.1), one can obtain [37]

$$ \mathcal{F}_{f_{\text{NL}}} = \sum_{A,B} \int dk \mathcal{N}^{-1}(k) \left\{ \begin{aligned} & 4 P_{\text{lin}}(k)^2 (b_{\phi}^B)^2 [P_{\epsilon}^A + P_{\epsilon}^B R(k)^2]^2 [b^B(k)]^2 \\ & + 8 P_{\text{lin}}(k)^2 P_{\epsilon}^B b_{\phi}^B [b_1^B b_{\phi}^A - b_1^A b_{\phi}^B] [P_{\epsilon}^A + P_{\epsilon}^B R(k)^2] R(k) b^B(k) \\ & + 2 P_{\text{lin}}(k)^2 [b_1^B b_{\phi}^A - b_1^A b_{\phi}^B]^2 [P_{\epsilon}^A P^B(k) + P_{\epsilon}^B \{P_{\epsilon}^B + P^B(k)\} R^2(k)] \end{aligned} \right\}, \quad (6.4) $$

where

$$ \mathcal{N}(k) \equiv \frac{2(2\pi)^3}{V_S V_k} \mathcal{M}(k, z_o)^2 \{P_{\epsilon}^A P^B(k) + P_{\epsilon}^B [P^B(k) - P_{\epsilon}^B] R(k)^2\}^2, \quad (6.5) $$

$$ R(k) \equiv \frac{b^A(k)}{b^B(k)} = \frac{b_1^A + b_\phi^A f_{\text{NL}} \mathcal{M}^{-1}(k, z_o)}{b_1^B + b_\phi^B f_{\text{NL}} \mathcal{M}^{-1}(k, z_o)}. \quad (6.6) $$

The key factor in this expression, encoding the dependence on the different bias parameters of the two traces is $|b_1^B b_\phi^A - b_1^A b_\phi^B|$ [29, 37], which can be maximized by selecting subsamples with different specific properties in order to minimize $\sigma_{f_{NL}}$. The proxy developed in section 4-5 can be used to select such suitable subsamples. Although the proxy is imperfect, as we show below, it performs very efficiently.

In particular, we use the proxy to select and combine large $b_\phi^A$ with small $b_\phi^B$ subsamples, as to maximize the mentioned factor. In figure 10 we report the results on $\sigma_{f_{NL}}$, as a function of the size of ELG- and LRG-like proxy-selected cumulative quantiles. The numerical values used for the computation of eq. (6.4) are the same as used in the single-tracer case. We refer to the subsample of LRG-like galaxies hosted by old halos as “old LRG” for brevity, and similarly for “young ELG”. We do not consider “old ELG” subsamples as their $\Delta b_\phi^{(r)}$ is smaller than in the LRG case.

The results in figure 10 illustrate the power of the multi-tracer technique. If we consider the LRG-like sample, as shown in figure 8, the error on $f_{NL}$ can be brought from $\sigma_{f_{NL}} \gtrsim 4$ to $\sigma_{f_{NL}} = 3.2$ by selecting a subsample of galaxies hosted by old halos, a marginal improvement. However, if this subsample is cross-correlated (via the multitracer approach) with a selection
<figure><img src="image_10.png" /><figcaption>Figure 10. Statistical error $\sigma_{f_{NL}}$ in the multi-tracer case. For the two samples ELG- and LRG-like, we consider the combination of cumulative quantiles of the oldest and youngest halos. The results are reported as a function of the size of the cumulative subsamples, with color coding depending on the value of $\sigma_{f_{NL}}$. On the left, we show the results of the ideal case in which $\omega_f$ is known. Results relative to proxy-selected subsamples are shown in the middle and on the right.</figcaption></figure>

of ELGs formed in young halos, as shown in figure 10, we can get an improvement of a factor of 4 in the statistical errors, obtaining $\sigma_{f_{NL}} \simeq 1$. Although this forecast may be further improved by adopting a better proxy, it is only $\simeq 30\%$ larger than the ideal case in which $\omega_f$ is known, as shown in the left panel of figure 10. Due to the broader $\omega_f$ distribution of the ELG-like sample, a subsample of young halos can be selected efficiently, while this does not apply for the LRG-like sample. Therefore, cross-correlating old and young halos hosting LRGs lead to suboptimal results, as illustrated in the right panel of figure 10.

In both cases, similarly to the single-tracer results, we can identify the optimal fractions which minimize $\sigma_{f_{NL}}$, roughly corresponding to selecting a third of the full sample. As already discussed in section 5.1, this is related to the tradeoff between signal and (shot) noise driven by the number density.

As regards the systematic error on $f_{NL}$, induced by systematic shifts in the estimated $b_\phi$ of the two multi-tracer samples, we refer to [29] and evaluate it by using the following expression

$$ \frac{\sigma_{f_{NL}}^{(sys.)}}{f_{NL}} = \frac{b_{\phi}^{A,t} b_{\phi}^{A,r} / \sigma_A^2 + b_{\phi}^{B,t} b_{\phi}^{B,r} / \sigma_B^2}{(b_{\phi}^{A,r})^2 / \sigma_A^2 + (b_{\phi}^{B,r})^2 / \sigma_B^2} - 1 \qquad (6.7) $$

where the subscript $t(r)$ indicates the true (recovered) $b_\phi$ of the sample and $\sigma_{A,B}$ stands for the single tracer's statistical error on $f_{NL}$.

Similarly to the single tracer case, the systematic error scales linearly with $f_{NL}$. Differently from the previous case, as discussed in [29], the multi-tracer technique should be more robust to systematic shifts in $b_\phi$. However, ref. [29] assumed both $b_\phi$ of the two samples to be shifted in the same direction, by the same amount, in which case their effects in eq. (6.7) cancel each other. In our case, instead, the residuals of the sample containing young halos have opposite sign with respect to those related to old halos, as shown in the top left panel of figure 9. Therefore, this cancellation does not apply here.

In figure 11 we show the ratio between the systematic error in $f_{NL}$ and the statistical one, as a function of the quantiles of the two samples, for different representative values of $f_{NL}$.

In the same way as illustrated in figure 9, as expected due to the linear scaling with $f_{NL}$ of the systematic error, for values as low as $f_{NL} \lesssim 5$, the systematic error is a fraction of the statistical one. For larger values of $f_{NL}$, the systematic error can be reduced to be
<figure><img src="image_11.png" /><figcaption>**Figure 11.** $\sigma_{f_{\text{NL}}}^{(\text{sys.})}$ in the multi-tracer case. The samples considered are the old LRG-like (large $b_{\phi}$) and the young ELG-like (small $b_{\phi}$) ones. The different panels illustrate the results for the representative values of $f_{\text{NL}} = 5, 10$, as a function of the size of the cumulative subsamples. In all the plots, the white color denotes the limit in which the systematic error is larger than the statistical one.</figcaption></figure>

below the statistical one by choosing a larger subsample size, for which the $b_{\phi}$ systematic shift is lower. Although this may be a suboptimal choice from the statistical error point of view, it is worth to highlight that if taken at face value, the multi-tracer systematic error is significantly lower than the single-tracer one.

Figure 12 is an attempt to summarize the main takeaways of this section. We report a selection of forecasts for two different representative values $f_{\text{NL}} = 5, 10$. We recall that statistical errors are for a volume of 20.4 Gpc$^3$/h$^3$.

We also reiterate that the systematic error estimate is also likely conservative, as we have taken the full residuals into account without modeling them.

# 7 Conclusions

The presence of local Primordial non-Gaussianity (PNG) in the initial conditions of the Universe affects dark matter clustering, inducing a scale dependence in the halo bias on large scales. This feature provides, in principle, a very competitive approach for constraining the amplitude of local PNG, parameterized by $f_{\text{NL}}$. However, the determination of $f_{\text{NL}}$ critically depends on the knowledge of the local PNG bias parameter $b_{\phi}$, currently the subject of ongoing discussions due to its sensitivity to assembly bias. Ref. [39] investigated the non-Gaussian halo assembly bias, explicitly showing the dependence of $b_{\phi}$ on halo formation time and providing theoretical predictions based on the extended Press-Schechter prescription. Our work serves as an extension of this prior research, highlighting how assembly bias can be modeled very accurately from first principles and priors on $b_{\phi}$ can be imposed by considering the conditional halo mass function. This is confirmed by tests on N-body simulations.

While theoretical predictions apply to dark matter halos, observations provide galaxy surveys. The poorly known halo-galaxy connection represents a limitation to the application of the theoretical modeling of $b_{\phi}$ to data and thus to obtain observational constraints on $f_{\text{NL}}$. This work provides a first attempt to bridge this gap.

We propose an observational proxy whereby photometric properties of galaxies can be used to infer statistical information on the assembly history of their host halos. The proxy
<figure><img src="image_12.png" /></figure>

**Figure 12.** Summary of the main $f_{\text{NL}}$ error forecast results reported throughout section 6. For each of the two representatives values $f_{\text{NL}} = 5, 10$, we show the forecasted error for the ideal and the proxy-selected quantiles, reporting both the statistical and the systematic error in the latter case. For the single tracer case (top two groups), results are shown for the 30 % quantile of oldest halos in the LRG and ELG-like samples. As for multi-tracer, we consider the “old LRG X young ELG” combination, with two different choices of quantile sizes. We recall that statistical errors are for a volume of 20.4 Gpc$^3/h^3$.

can then be used to identify assembly-bias-selected subsamples and provide an estimate of the subsamples' $b_\phi$. In doing so we also identify which types of galaxy samples represent optimal targets to provide more precise constraints on $f_{\text{NL}}$.

Our main results can be summarized as follows:

* At fixed mass, halos that accreted their mass earlier have a larger $b_\phi$ than predicted by the universality relation and viceversa for late-forming halos. The departure from universality can be related to the conditional halo mass function, independently of its specific functional form (this effect goes under the name of assembly bias).

* The ePS formalism allows us to derive an analytical prediction for the conditional mass function, which is independent of halo mass and redshift. However, its accuracy decreases when considering halos of mass $M \lesssim 10^{13} M_\odot h^{-1}$, while a simple 1-parameter extension shows a remarkably good fit to simulations across all relevant mass and redshift ranges. Hence the assembly bias and thus $b_\phi$ for a given halo sample depends through a simple (physically motivated and calibrated on N-body simulations) analytic function on a single variable $\omega_f$, which encodes a specific property of the halo sample assembly history.

* Using state-of-the art cosmological simulations that model both the dark matter clustering and galaxy properties (IllustrisTNG and the CAMELS-TNG suite) we calibrated an observational proxy for $\omega_f$, the quantity which determines halo assembly bias and thus $b_\phi$. The proxy is built (by design) from a linear combination of photometric bands (i.e.
colors) of the galaxies and tuned to minimize the difference between the true $b_\phi$ of a given halo sample and the one recovered from the photometry of the halos' central galaxies.

* The proxy is used to identify subsamples which are (statistically) rank ordered according to their $\omega_f$ and thus $b_\phi$ and to provide an estimate of their $b_\phi$. The proxy is of course not perfect: $\omega_f$ is recovered with some scatter and $b_\phi$ is recovered with some dispersion and a systematic shift which depends on the sample selection.

* Marginalization over cosmological and IllustrisTNG astrophysical parameters (performed with the CAMELS-TNG simulations suite) does not degrade the proxy performance significantly. However, further work should be done to assess the robustness on galaxy formation models beyond IllustrisTNG.

* The statistical error on the $b_\phi$ recovered from the proxy does not degrade significantly the forecasted errors on $f_{\text{NL}}$ compared to the ideal case when $b_\phi$ is perfectly known.

* The systematic error on the proxy-recovered $b_\phi$ is remarkably regular across different galaxy selection strategies. This indicates that the shift could be modeled and subtracted. The proxy performance could probably be improved by generalizing it to a non-linear combination of galaxies' colors. Even without these improvements, and taking the full uncorrected shift in $b_\phi$ as a source of systematic error on $f_{\text{NL}}$, we show that the resulting systematic error on $f_{\text{NL}}$ can be kept reasonably under control.

* The proxy can also be used to observationally select halo subsamples with $b_\phi$ as different as possible as to optimize the constraints on $f_{\text{NL}}$ with the multi-tracer approach. In this case, forecasted errors on $f_{\text{NL}}$ reach the level of $\pm\mathcal{O}(1)$ for volumes $\sim \mathcal{O}(20)$ Gpc$^3$.

We hope that the findings of this work will serve to demystify the impact of assembly bias on PNG constraints and will motivate improvements on the simple observational proxy proposed here. Even just with the naive, linear combination of central galaxy colors proxy, the power spectrum of biased tracers offer a competitive and viable window into PNG.

## Acknowledgments

EF acknowledges the support from “la Caixa” Foundation (ID 100010434, code LCF/BQ/DI21/11860061). Funding for this work was partially provided by project PGC2018-098866-B-I00 MCIN/AEI/10.13039/501100011033 y FEDER “Una manera de hacer Europa”, and the “Center of Excellence Maria de Maeztu 2020-2023” award to the ICCUB (CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033). LV acknowledges support of European Union’s Horizon 2020 research and innovation programme ERC (BePreSysE, grant agreement 725327). The CAMELS project is supported by the Simons Foundation and NSF grant AST 2108078. AR acknowledges support from PRIN-MIUR 2020 METE, under contract no. 2020KB33TP. DK is supported by the South African Radio Astronomy Observatory and the National Research Foundation (Grant No. 75415). GJ acknowledges support from the ANR LOCALIZATION project, grant ANR-21-CE31-0019 / 490702358 of the French Agence Nationale de la Recherche. The Center for Computational Astrophysics and the Flatiron Institute are supported by the Simons Foundation.
<figure><img src="image_13.png" /><figcaption>**Figure 13.** Pearson correlation coefficient $r(Y, \omega_f)$ between the observational proxy $Y$ and $\omega_f$, as a function of the number of colors included in the linear combination which defines $Y$. The results are reported for the two representative samples of ELG- and LRG-like galaxies in the IllustrisTNG simulations, at $z_o = 1$.</figcaption></figure>

## A Optimizing the number of galaxy colors combined in Y

As discussed in section 4, in order to estimate $\omega_f$ from observational galaxy properties, we make use of the information contained in IllustrisTNG galaxy colors. In this appendix we show that, as long as a simple linear combination $Y = \sum_{i=1}^{N} a_i C_i$ is considered, $N = 3$ colors are sufficient to maximize the Pearson correlation $r(Y, \omega_f)$. In particular, in figure 13 we report $r$ as a function of the number of colors $N$ for two representative galaxy samples, the ELG- and LRG-like in IllustrisTNG at $z_o = 1$.

As illustrated in figure 13, adding colors to the combination improves the correlation between the optimal proxy $Y$ and $\omega_f$, until $N = 3$. Then $r$ saturates and the addition of further colors becomes redundant. Specifically, for the LRG-like sample, $r$ only increases by 0.6% passing from a combination of 3 colors to a combination of 5, while the computational time needed to optimize the model is 30 times larger. Similar results apply for what regards the ELG-like sample.

Consequently, for all the galaxy samples considered in section 4, we restrict our optimization of the proxy to a linear combination of 3 galaxy colors.

## References

[1] J.M. Maldacena, *Non-Gaussian features of primordial fluctuations in single field inflationary models*, JHEP 05 (2003) 013 [astro-ph/0210603] [INSPIRE].

[2] PLANCK collaboration, *Planck 2018 results. IX. Constraints on primordial non-Gaussianity*, Astron. Astrophys. 641 (2020) A9 [arXiv:1905.05697] [INSPIRE].
[3] E.-M. Mueller et al., *Primordial non-Gaussianity from the completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey II: measurements in Fourier space with optimal weights*, Mon. Not. Roy. Astron. Soc. **514** (2022) 3396 [arXiv:2106.13725] [INSPIRE].

[4] DESI collaboration, *The DESI experiment part I: science, targeting, and survey design*, arXiv:1611.00036 [INSPIRE].

[5] L. Amendola et al., *Cosmology and fundamental physics with the Euclid satellite*, Living Rev. Rel. **21** (2018) 2 [arXiv:1606.00180] [INSPIRE].

[6] SPHEREx collaboration, *Cosmology with the SPHEREX all-sky spectral survey*, arXiv:1412.4872 [INSPIRE].

[7] D.J. Schlegel et al., *Astro2020 APC white paper. The MegaMapper: a $z > 2$ spectroscopic instrument for the study of inflation and dark energy*, Bull. Am. Astron. Soc. **51** (2019) 229 [arXiv:1907.11171] [INSPIRE].

[8] T. Giannantonio et al., *Constraining primordial non-Gaussianity with future galaxy surveys*, Mon. Not. Roy. Astron. Soc. **422** (2012) 2854 [arXiv:1109.0958] [INSPIRE].

[9] M. Alvarez et al., *Testing inflation with large scale structure: connecting hopes with reality*, arXiv:1412.4671 [INSPIRE].

[10] R. de Putter and O. Doré, *Designing an inflation galaxy survey: how to measure $\sigma(f_{\text{NL}}) \sim 1$ using scale-dependent galaxy bias*, Phys. Rev. D **95** (2017) 123513 [arXiv:1412.3854] [INSPIRE].

[11] D. Karagiannis et al., *Constraining primordial non-Gaussianity with bispectrum and power spectrum from upcoming optical and radio surveys*, Mon. Not. Roy. Astron. Soc. **478** (2018) 1341 [arXiv:1801.09280] [INSPIRE].

[12] S. Ferraro, N. Sailer, A. Slosar and M. White, *Snowmass2021 cosmic frontier white paper: cosmology and fundamental physics from the three-dimensional large scale structure*, arXiv:2203.07506 [INSPIRE].

[13] A. Achúcarro et al., *Inflation: theory and observations*, arXiv:2203.08128 [INSPIRE].

[14] A. Gangui, F. Lucchin, S. Matarrese and S. Mollerach, *The three point correlation function of the cosmic microwave background in inflationary models*, Astrophys. J. **430** (1994) 447 [astro-ph/9312033] [INSPIRE].

[15] L.-M. Wang and M. Kamionkowski, *The cosmic microwave background bispectrum and inflation*, Phys. Rev. D **61** (2000) 063504 [astro-ph/9907431] [INSPIRE].

[16] F. Bernardeau and J.-P. Uzan, *Non-Gaussianity in multifield inflation*, Phys. Rev. D **66** (2002) 103506 [hep-ph/0207295] [INSPIRE].

[17] L. Verde, L.-M. Wang, A. Heavens and M. Kamionkowski, *Large scale structure, the cosmic microwave background, and primordial non-Gaussianity*, Mon. Not. Roy. Astron. Soc. **313** (2000) L141 [astro-ph/9906301] [INSPIRE].

[18] E. Komatsu and D.N. Spergel, *Acoustic signatures in the primary microwave background bispectrum*, Phys. Rev. D **63** (2001) 063002 [astro-ph/0005036] [INSPIRE].

[19] N. Dalal, O. Doré, D. Huterer and A. Shirokov, *The imprints of primordial non-gaussianities on large-scale structure: scale dependent bias and abundance of virialized objects*, Phys. Rev. D **77** (2008) 123514 [arXiv:0710.4560] [INSPIRE].

[20] S. Matarrese and L. Verde, *The effect of primordial non-Gaussianity on halo bias*, Astrophys. J. Lett. **677** (2008) L77 [arXiv:0801.4826] [INSPIRE].
[21] A. Slosar et al., Constraints on local primordial non-Gaussianity from large scale structure, JCAP 08 (2008) 031 [arXiv:0805.3580] [INSPIRE].

[22] N. Afshordi and A.J. Tolley, Primordial non-Gaussianity, statistics of collapsed objects, and the Integrated Sachs-Wolfe effect, Phys. Rev. D 78 (2008) 123507 [arXiv:0806.1046] [INSPIRE].

[23] P. McDonald, Primordial non-Gaussianity: large-scale structure signature in the perturbative bias model, Phys. Rev. D 78 (2008) 123519 [arXiv:0806.1061] [INSPIRE].

[24] P. Creminelli and M. Zaldarriaga, Single field consistency relation for the 3-point function, JCAP 10 (2004) 006 [astro-ph/0407059] [INSPIRE].

[25] J.R. Bond, S. Cole, G. Efstathiou and N. Kaiser, Excursion set mass functions for hierarchical Gaussian fluctuations, Astrophys. J. 379 (1991) 440 [INSPIRE].

[26] W.H. Press and P. Schechter, Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation, Astrophys. J. 187 (1974) 425 [INSPIRE].

[27] R.K. Sheth and G. Tormen, Large scale bias and the peak background split, Mon. Not. Roy. Astron. Soc. 308 (1999) 119 [astro-ph/9901122] [INSPIRE].

[28] T. Lazeyras, A. Barreira, F. Schmidt and V. Desjacques, Assembly bias in the local PNG halo bias and its implication for $f_{NL}$ constraints, JCAP 01 (2023) 023 [arXiv:2209.07251] [INSPIRE].

[29] A. Barreira and E. Krause, Towards optimal and robust $f_{NL}$ constraints with multi-tracer analyses, JCAP 10 (2023) 044 [arXiv:2302.09066] [INSPIRE].

[30] J.M. Sullivan, T. Prijon and U. Seljak, Learning to concentrate: multi-tracer forecasts on local primordial non-Gaussianity with machine-learned bias, JCAP 08 (2023) 004 [arXiv:2303.08901] [INSPIRE].

[31] A. Barreira, On the impact of galaxy bias uncertainties on primordial non-Gaussianity constraints, JCAP 12 (2020) 031 [arXiv:2009.06622] [INSPIRE].

[32] A. Barreira et al., Galaxy bias and primordial non-Gaussianity: insights from galaxy formation simulations with IllustrisTNG, JCAP 12 (2020) 013 [arXiv:2006.09368] [INSPIRE].

[33] A. Barreira, Can we actually constrain $f_{NL}$ using the scale-dependent bias effect? An illustration of the impact of galaxy bias uncertainties using the BOSS DR12 galaxy power spectrum, JCAP 11 (2022) 013 [arXiv:2205.05673] [INSPIRE].

[34] L. Lucie-Smith, A. Barreira and F. Schmidt, Halo assembly bias from a deep learning model of halo formation, Mon. Not. Roy. Astron. Soc. 524 (2023) 1746 [arXiv:2304.09880] [INSPIRE].

[35] G. Jung et al., Quijote-PNG: the information content of the halo mass function, Astrophys. J. 957 (2023) 1 [arXiv:2305.10597] [INSPIRE].

[36] U. Seljak, Extracting primordial non-Gaussianity without cosmic variance, Phys. Rev. Lett. 102 (2009) 021302 [arXiv:0807.1770] [INSPIRE].

[37] D. Karagiannis et al., Multi-tracer power spectra and bispectra: formalism, arXiv:2305.04028 [INSPIRE].

[38] A. Slosar, Optimal dataset combining in $f_{NL}$ constraints from large scale structure, JCAP 03 (2009) 004 [arXiv:0808.0044] [INSPIRE].

[39] B.A. Reid et al., Non-Gaussian halo assembly bias, JCAP 07 (2010) 013 [arXiv:1004.1637] [INSPIRE].

[40] C. Lacey and S. Cole, Merger rates in hierarchical models of galaxy formation, Mon. Not. Roy. Astron. Soc. 262 (1993) 627 [INSPIRE].
[41] C.G. Lacey and S. Cole, *Merger rates in hierarchical models of galaxy formation*. 2. *Comparison with N body simulations*, Mon. Not. Roy. Astron. Soc. **271** (1994) 676 [astro-ph/9402069] [INSPIRE].

[42] F.C. Bosch, *The universal mass accretion history of CDM haloes*, Mon. Not. Roy. Astron. Soc. **331** (2002) 98 [astro-ph/0105158] [INSPIRE].

[43] F. Villaescusa-Navarro et al., *The Quijote simulations*, Astrophys. J. Suppl. **250** (2020) 2 [arXiv:1909.05273] [INSPIRE].

[44] W.R. Coulton et al., *Quijote-PNG: the information content of the halo power spectrum and bispectrum*, Astrophys. J. **943** (2023) 178 [arXiv:2206.15450] [INSPIRE].

[45] W.R. Coulton et al., *Quijote-PNG: simulations of primordial non-Gaussianity and the information content of the matter field power spectrum and bispectrum*, Astrophys. J. **943** (2023) 64 [arXiv:2206.01619] [INSPIRE].

[46] G. Jung et al., *Quijote-PNG: quasi-maximum likelihood estimation of primordial non-Gaussianity in the nonlinear dark matter density field*, Astrophys. J. **940** (2022) 71 [arXiv:2206.01624] [INSPIRE].

[47] D. Nelson et al., *The IllustrisTNG simulations: public data release*, arXiv:1812.05609 [INSPIRE].

[48] D. Nelson et al., *First results from the IllustrisTNG simulations: the galaxy colour bimodality*, Mon. Not. Roy. Astron. Soc. **475** (2018) 624 [arXiv:1707.03395] [INSPIRE].

[49] CAMELS collaboration, *The CAMELS project: Cosmology and Astrophysics with MachinE Learning Simulations*, Astrophys. J. **915** (2021) 71 [arXiv:2010.00619] [INSPIRE].

[50] S. Cole and N. Kaiser, *Biased clustering in the cold dark matter cosmogony*, Mon. Not. Roy. Astron. Soc. **237** (1989) 1127 [INSPIRE].

[51] V. Desjacques, U. Seljak and I. Iliev, *Scale-dependent bias induced by local non-Gaussianity: a comparison to N-body simulations*, Mon. Not. Roy. Astron. Soc. **396** (2009) 85 [arXiv:0811.2748] [INSPIRE].

[52] C. Wagner and L. Verde, *N-body simulations with generic non-Gaussian initial conditions*. Part II. Halo bias, JCAP **03** (2012) 002 [arXiv:1102.3229] [INSPIRE].

[53] S. Matarrese, L. Verde and R. Jimenez, *The abundance of high-redshift objects as a probe of non-Gaussian initial conditions*, Astrophys. J. **541** (2000) 10 [astro-ph/0001366] [INSPIRE].

[54] M. LoVerde, A. Miller, S. Shandera and L. Verde, *Effects of scale-dependent non-Gaussianity on cosmological structures*, JCAP **04** (2008) 014 [arXiv:0711.4126] [INSPIRE].

[55] M. Maggiore and A. Riotto, *The halo mass function from excursion set theory*. Part III. Non-Gaussian fluctuations, Astrophys. J. **717** (2010) 526 [arXiv:0903.1251] [INSPIRE].

[56] G. D'Amico, M. Musso, J. Noreña and A. Paranjape, *An improved calculation of the non-Gaussian halo mass function*, JCAP **02** (2011) 001 [arXiv:1005.1203] [INSPIRE].

[57] R.K. Sheth, H.J. Mo and G. Tormen, *Ellipsoidal collapse and an improved model for the number and spatial distribution of dark matter haloes*, Mon. Not. Roy. Astron. Soc. **323** (2001) 1 [astro-ph/9907024] [INSPIRE].

[58] A. Jenkins et al., *The mass function of dark matter halos*, Mon. Not. Roy. Astron. Soc. **321** (2001) 372 [astro-ph/0005260] [INSPIRE].

[59] A. Paranjape and R.K. Sheth, *Peaks theory and the excursion set approach*, Mon. Not. Roy. Astron. Soc. **426** (2012) 2789 [arXiv:1206.3506] [INSPIRE].
[60] A. Paranjape, R.K. Sheth and V. Desjacques, *Excursion set peaks: a self-consistent model of dark halo abundances and clustering*, Mon. Not. Roy. Astron. Soc. **431** (2013) 1503 [arXiv:1210.1483] [INSPIRE].

[61] E. Castorina, A. Paranjape, O. Hahn and R.K. Sheth, *Excursion set peaks: the role of shear*, arXiv:1611.03619 [INSPIRE].

[62] M. Grossi et al., *Large-scale non-Gaussian mass function and halo bias: tests on N-body simulations*, Mon. Not. Roy. Astron. Soc. **398** (2009) 321 [arXiv:0902.2013] [INSPIRE].

[63] A. Pillepich, C. Porciani and O. Hahn, *Universal halo mass function and scale-dependent bias from N-body simulations with non-Gaussian initial conditions*, Mon. Not. Roy. Astron. Soc. **402** (2010) 191 [arXiv:0811.4176] [INSPIRE].

[64] M. Biagetti et al., *Verifying the consistency relation for the scale-dependent bias from local primordial non-Gaussianity*, Mon. Not. Roy. Astron. Soc. **468** (2017) 3277 [arXiv:1611.04901] [INSPIRE].

[65] E. Alizadeh and B. Wandelt, *Dark matter halo merger and accretion probabilities in the excursion set formalism*, arXiv:0810.2301 [INSPIRE].

[66] P.T.P. Viana and A.R. Liddle, *The cluster abundance in flat and open cosmologies*, Mon. Not. Roy. Astron. Soc. **281** (1996) 323 [astro-ph/9511007] [INSPIRE].

[67] J.L. Tinker et al., *Toward a halo mass function for precision cosmology: the limits of universality*, Astrophys. J. **688** (2008) 709 [arXiv:0803.2706] [INSPIRE].

[68] T. Lazeyras, M. Musso and F. Schmidt, *Large-scale assembly bias of dark matter halos*, JCAP **03** (2017) 059 [arXiv:1612.04360] [INSPIRE].

[69] R.H. Wechsler and J.L. Tinker, *The connection between galaxies and their dark matter halos*, Ann. Rev. Astron. Astrophys. **56** (2018) 435 [arXiv:1804.03097] [INSPIRE].

[70] D.J. Croton, L. Gao and S.D.M. White, *Halo assembly bias and its effects on galaxy clustering*, Mon. Not. Roy. Astron. Soc. **374** (2007) 1303 [astro-ph/0605636] [INSPIRE].

[71] DES collaboration, *Predicting galaxy star formation rates via the co-evolution of galaxies and haloes*, Mon. Not. Roy. Astron. Soc. **446** (2015) 651 [arXiv:1403.1578] [INSPIRE].

[72] J.L. Tinker, C. Hahn, Y.-Y. Mao and A.R. Wetzel, *Halo histories versus galaxy properties at z = 0 — III. The properties of star-forming galaxies*, Mon. Not. Roy. Astron. Soc. **478** (2018) 4487.

[73] A.P. Hearin and D.F. Watson, *The dark side of galaxy color*, Mon. Not. Roy. Astron. Soc. **435** (2013) 1313 [arXiv:1304.5557] [INSPIRE].

[74] SDSS collaboration, *Star formation rate indicators in the Sloan Digital Sky Survey*, Astrophys. J. **599** (2003) 971 [astro-ph/0306621] [INSPIRE].

[75] P.D. Mitchell, C.G. Lacey, C.M. Baugh and S. Cole, *How well can we really estimate the stellar masses of galaxies from broadband photometry?*, Mon. Not. Roy. Astron. Soc. **435** (2013) 87 [arXiv:1303.7228] [INSPIRE].

[76] DESI collaboration, *Target selection and validation of DESI luminous red galaxies*, Astron. J. **165** (2023) 58 [arXiv:2208.08515] [INSPIRE].

[77] A. Raichoor et al., *Target selection and validation of DESI emission line galaxies*, Astron. J. **165** (2023) 126 [arXiv:2208.08513] [INSPIRE].
[78] B. Hadzhiyska, S. Tacchella, S. Bose and D.J. Eisenstein, *The galaxy-halo connection of emission-line galaxies in IllustrisTNG*, Mon. Not. Roy. Astron. Soc. **502** (2021) 3599 [arXiv:2011.05331] [INSPIRE].

[79] S. Yuan, B. Hadzhiyska, S. Bose and D.J. Eisenstein, *Illustrating galaxy-halo connection in the DESI era with illustrisTNG*, Mon. Not. Roy. Astron. Soc. **512** (2022) 5793 [arXiv:2202.12911] [INSPIRE].

[80] E.F. Bell, D.H. McIntosh, N. Katz and M.D. Weinberg, *The optical and near-infrared properties of galaxies*. 1. Luminosity and stellar mass functions, Astrophys. J. Suppl. **149** (2003) 289 [astro-ph/0302543] [INSPIRE].

[81] R. Scoccimarro, R.K. Sheth, L. Hui and B. Jain, *How many galaxies fit in a halo? Constraints on galaxy formation efficiency from spatial clustering*, Astrophys. J. **546** (2001) 20 [astro-ph/0006319] [INSPIRE].

[82] J.L. Tinker et al., *The large scale bias of dark matter halos: numerical calibration and model tests*, Astrophys. J. **724** (2010) 878 [arXiv:1001.3162] [INSPIRE].

[83] A. Pillepich et al., *Simulating galaxy formation with the IllustrisTNG model*, Mon. Not. Roy. Astron. Soc. **473** (2018) 4077 [arXiv:1703.02970] [INSPIRE].

[84] CAMELS collaboration, *The CAMELS project: expanding the galaxy formation model space with new ASTRID and 28-parameter TNG and SIMBA suites*, Astrophys. J. **959** (2023) 136 [arXiv:2304.02096] [INSPIRE].

[85] C. Carbone, L. Verde and S. Matarrese, *Non-Gaussian halo bias and future galaxy surveys*, Astrophys. J. Lett. **684** (2008) L1 [arXiv:0806.1950] [INSPIRE].