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

HAL
open science

# Impacts of the physical data model on the forward inference of initial conditions from biased tracers

Nhat-Minh Nguyen, Fabian Schmidt, Guilhem Lavaux, Jens Jasche

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

## To cite this version:

Nhat-Minh Nguyen, Fabian Schmidt, Guilhem Lavaux, Jens Jasche. Impacts of the physical data model on the forward inference of initial conditions from biased tracers. Journal of Cosmology and Astroparticle Physics, 2021, 03, pp.058. 10.1088/1475-7516/2021/03/058. hal-03047630

HAL Id: hal-03047630

https://hal.science/hal-03047630v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Impacts of the physical data model on the forward inference of initial conditions from biased tracers

To cite this article: Nhat-Minh Nguyen et al JCAP03(2021)058

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

You may also like

- **PATCHY BLAZAR HEATING; DIVERSIFYING THE THERMAL HISTORY OF THE INTERGALACTIC MEDIUM**
Astrid Lamberts, Philip Chang, Christoph Pfrommer et al.

- **The Bias from Hydrodynamic Simulations: Mapping Baryon Physics onto Dark Matter Fields**
Francesco Sinigaglia, Francisco-Shu Kitaura, A. Balaguera-Antolínez et al.

- **Inferring the drivers of language change using spatial models**
James Burridge and Tamsin Blaxter
# Impacts of the physical data model on the *forward* inference of initial conditions from biased tracers

Nhat-Minh Nguyen,ª Fabian Schmidt,ª Guilhem Lavauxᵇ and Jens Jascheᶜ

E-mail: minh@mpa-garching.mpg.de, fabians@mpa-garching.mpg.de, guilhem.lavaux@iap.fr, jens.jasche@fysik.su.se

Received November 25, 2020

Accepted February 3, 2021

Published March 17, 2021

**Abstract.** We investigate the impact of each ingredient in the employed physical data model on the Bayesian forward inference of initial conditions from biased tracers at the field level. Specifically, we use dark matter halos in a given cosmological simulation volume as tracers of the underlying matter density field. We study the effect of tracer density, grid resolution, gravity model, bias model and likelihood on the inferred initial conditions. We find that the cross-correlation coefficient between true and inferred phases reacts weakly to all ingredients above, and is well predicted by the theoretical expectation derived from a Gaussian model on a broad range of scales. The bias in the amplitude of the inferred initial conditions, on the other hand, depends strongly on the bias model and the likelihood. We conclude that the bias model and likelihood hold the key to an unbiased cosmological inference. Together they must keep the systematics — which arise from the sub-grid physics that are marginalized over — under control in order to obtain an unbiased inference.

**Keywords:** cosmological parameters from LSS, galaxy clustering, redshift surveys

**ArXiv ePrint:** 2011.06587

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

<footnote>ªMax-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany</footnote>

<footnote>ᵇSorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote>ᶜDepartment of Physics, Stockholm University, Albanova University Centre, SE-106 91 Stockholm, Sweden</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Inference setup</td><td>6</td></tr><tr><td>2.1</td><td>Gravity model</td><td>6</td></tr><tr><td>2.2</td><td>Bias model</td><td>6</td></tr><tr><td>2.2.1</td><td>Linear bias</td><td>6</td></tr><tr><td>2.2.2</td><td>Power-law bias</td><td>7</td></tr><tr><td>2.2.3</td><td>Broken power-law bias</td><td>7</td></tr><tr><td>2.3</td><td>Likelihood</td><td>8</td></tr><tr><td>2.3.1</td><td>Poisson likelihood</td><td>8</td></tr><tr><td>2.3.2</td><td>Gaussian likelihood</td><td>8</td></tr><tr><td>2.4</td><td>Grid resolution</td><td>9</td></tr><tr><td>3</td><td>N-body simulation setup</td><td>9</td></tr><tr><td>4</td><td>Results</td><td>9</td></tr><tr><td>4.1</td><td>Results with mock data</td><td>11</td></tr><tr><td>4.1.1</td><td>Tracer number density</td><td>12</td></tr><tr><td>4.1.2</td><td>Grid resolution</td><td>14</td></tr><tr><td>4.2</td><td>Results with N-body data</td><td>14</td></tr><tr><td>4.2.1</td><td>Gravity model</td><td>15</td></tr><tr><td>4.2.2</td><td>Bias model</td><td>16</td></tr><tr><td>4.2.3</td><td>Likelihood</td><td>18</td></tr><tr><td>5</td><td>Discussion</td><td>21</td></tr><tr><td>6</td><td>Conclusions</td><td>22</td></tr><tr><td>A</td><td>Gaussian expectation of the Fourier phase-correlation</td><td>23</td></tr><tr><td>B</td><td>Residual bias due to normalization using σ<sub>8</sub></td><td>23</td></tr><tr><td>C</td><td>Priors and thresholding</td><td>24</td></tr><tr><td>D</td><td>MCMC burn-in phase and chain thinning</td><td>25</td></tr></table>

## 1 Introduction

The distribution of large-scale structures (LSS) in our Universe, specifically galaxies — biased tracers of the underlying matter density field — is very far from the random distribution of initial quantum fluctuations that seed their formation. As complex as the process of galaxy formation is, on large and quasi-linear scales, the total matter distribution follows that of dark matter (DM) which only feels the effect of gravity. Going forward in time, the gravitational evolution of matter distribution is fully *deterministic*. It follows that, by forward-modeling this process, together with the relation between matter and galaxy distributions, one should
be able to retrieve a vast amount of information directly from the three-dimensional distributions of galaxies as probed by galaxy redshift surveys. This includes, but is not limited to the initial conditions of our observed patch of the Universe and parameters in our cosmological models.

Standard approaches in cosmological inference and reconstruction from LSS data, however, often rely on either modeling of *n*-point functions of galaxies [e.g. 1–4] or reversing the evolution of matter fluctuations by moving mass particles or galaxies backward (reconstruction hereafter) [e.g. 5–8]. Here, we are interested in the aforementioned more direct route: *forward modeling of the matter and galaxy clustering directly at the field-level* [9–24]. This alternative path not only circumvents the difficult and cumbersome problem of how to accurately predict and compute the higher *n*-point functions (plus their covariances) from theory [25–29] and data [1, 30], respectively, but also, in principle, accesses all information available in the three-dimensional data field above the smoothing scale, regardless of its Gaussianity. In fact, the latter statement is not guaranteed to apply to the *n*-point function approach if the observed galaxy field is non-Gaussian at the scales of interest [31].

There are four principal ingredients of the physical data model in such a field-level, forward inference approach [9, 18]:

1. **Prior for initial conditions:** this element encodes our prior knowledge of the statistical distribution of matter density fluctuations at the very early stage of the Universe that comes from the theory of initial fluctuations [32, 33] and observation of the Cosmic Microwave Background (CMB) [34, 35].

2. **Gravity model for matter:** this ingredient evolves matter distributions forward in time, from a specific set of initial conditions, under the effect of gravity.

3. **Bias model for deterministic matter-galaxy bias:** this constituent relates the forward-evolved matter field and the deterministic galaxy field by means of a deterministic relation (see [36] for a review).

4. **Likelihood for stochastic galaxy bias:** this final piece attempts to model the scatter in the deterministic matter-galaxy bias relation, i.e. the difference between the deterministic and observed galaxy fields, in the form of a conditional probability [18–21].

We summarize the connection between these ingredients of the physical data model in figure 1, where (and hereafter) we have generalized “galaxy” to “biased tracer”, partially due to the fact that our results are obtained using DM halos in N-body simulations, and also because most of our results and arguments apply to all generic biased tracers. We will discuss possible exceptions at the end of section 4.2.3.

How strongly does the quality of a forward inference depend on these ingredients? In particular, how sensitive are the inferred initial conditions to different choices for each ingredient? So far, this is an open question in the context of forward modeling (see [37] for a case study in the “backward” modeling approach). It is the aim of this paper to address that question. Specifically, inspired by previous findings in [22–24], our main goal is to *identify the key ingredients for a robust, unbiased cosmological forward inference from LSS*, and what should be the focus of future studies to realize this goal.

In this paper, we specifically measure: the cross-correlation coefficient between true and inferred phases (phase-correlation hereafter) $r(k)$ and the amplitude-bias in the corresponding power spectra (which is often referred to by the literature of reconstruction, and hence
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Simplified flowchart of the forward inference framework, as relevant for this work. The input data is illustrated by the top red panel. The target of the inference is represented by the blue diamond panel in the middle row. Different ingredients of the physical data model are displayed in the bottom row. We highlight the ingredients inspected in this work in red. Different options for each inspected ingredient are explicitly listed in the corresponding panels. Our only anchored assumption is the Gaussian prior on initial conditions, shown in the round-edge gray panel.</figcaption></figure>

hereafter as “transfer function”) $T_{\text{infer}}(k)$. Clearly, the former is important for studies making use of the inference for cross-correlations with other probes. The latter is important as it indicates systematic biases in cosmological parameter that are to be expected. Specifically, a transfer function that is biased in a *scale-independent* way implies a bias in the inferred primordial normalization of the power spectrum $A_s$, whereas a *scale-dependently biased* transfer function entails biases in not only $A_s$ but other cosmological parameters, e.g. $\Omega_m$, $n_s$, as well.

We use DM main halos identified in N-body simulations, for which the true initial conditions are known, as physical biased tracers. We further restrict our analysis to the halo rest frame, i.e. the case without redshift space distortion (RSD). We consider multiple inference setups corresponding to different choices for the ingredients two to four in the above list (bottom red panels in figure 1), as well as different options for halo number density and grid resolution, i.e. smoothing scale.

All forward inferences in this paper follow exactly the flowchart in figure 1, and are carried out using the BORG algorithm that was first introduced in [9]. As our parameter space actually consists of the whole three-dimensional field of initial conditions, augmented by the bias parameters, efficiently Monte-Carlo Markov-chain (MCMC) sampling this very high-dimensional space of the posterior is a truly daunting task. The task is even more
critical in the particular problem of forward inference, for each sample requires essentially a full simulation from the early-time Gaussian initial conditions to the late-time non-Gaussian observed data, as diagramed in figure 1, hence some non-negligible numerical expense. The BORG algorithm addresses this issue by utilizing the Hamiltonian Monte Carlo (HMC) sampling method [38–40] which combines geometrical information of the posterior, encoded in its gradient, and Hamiltonian dynamics. The latter ensures that the Metropolis-Hastings proposal always has a high acceptance rate — providing that errors in the numerical integration of the Hamiltonian equation of motion are kept under control [39, 41]. We encourage readers to see [9] for further conceptual and technical details specifically to BORG. We note that, in an application on real data from a galaxy redshift survey, BORG is also capable of routinely dealing with multiple sources of systematics, including survey geometries, selection functions, foreground contaminations and RSD effects [See, for example 17]. These steps are not included in figure 1 as they are not relevant for our analysis using rest-frame halos in N-body simulations.

Our results in section 4 are based on our analyses of the outputs of multiple BORG runs: each is an ensemble of MCMC samples that resembles the corresponding posterior density distribution. In our particular case, after marginalizing over bias parameters, the posterior parameter space is the three-dimensional field of initial conditions (middle blue diamond panel in figure 1), physically constrained by the three-dimensional input field of halos (top red square panel in figure 1) and our choices of models for each ingredient (bottom red rectangle panels).

The rest of this paper is organized as follows. In section 2, we describe in more details the three ingredients under investigation and other settings that together characterize the configuration of each inference. We further detail the N-body simulations and the halos that are used as the input data in section 3. In section 4, we first define the estimators used to assess the quality of the inference and later present our results. We then highlight the ingredients that demand better modeling in section 5, and conclude in section 6. The appendices provide additional information about the derivation of the Gaussian expectation in section 4, numerical implementation details and bias priors relevant for section 4.2, as well as MCMC chains analyzed throughout this paper.

## Convention and notation

Our Fourier notation follows that of [18, 36]:

$$ 
\begin{align} 
f(\boldsymbol{k}) &\equiv \int d^3\boldsymbol{x} f(\boldsymbol{x}) e^{-i\boldsymbol{k}\cdot\boldsymbol{x}} &\equiv \int_{\boldsymbol{x}} f(\boldsymbol{x}) e^{-i\boldsymbol{k}\cdot\boldsymbol{x}} \tag{1.1} \\ 
f(\boldsymbol{x}) &\equiv \int \frac{d^3\boldsymbol{k}}{(2\pi)^3} f(\boldsymbol{k}) e^{i\boldsymbol{k}\cdot\boldsymbol{x}} &\equiv \int_{\boldsymbol{k}} f(\boldsymbol{k}) e^{i\boldsymbol{k}\cdot\boldsymbol{x}}. 
\end{align} 
$$

For ensemble averages, we distinguish between $\langle O \rangle_k$, which denotes an average of $O$ over Fourier wavenumber $k$, and $\langle O \rangle_s$, which denotes an average of $O$ over MCMC samples. A spatial-average of $O$, on the other hand, is denoted by $\bar{O}$.

The matter and halo density contrast fields (on a fixed time slice) are defined as

$$ \delta_m(\boldsymbol{x}) \equiv \frac{\rho_m(\boldsymbol{x})}{\bar{\rho}_m} - 1, \quad \delta_h(\boldsymbol{x}) \equiv \frac{n_h(\boldsymbol{x})}{\bar{n}_h} - 1. \qquad (1.2) $$

Throughout the paper, we assume adiabatic, growing-mode initial conditions such that the initial fluctuations in the matter field can be described by only a single field $\delta_{m,ini}$.
Further, we adopt a flat $\Lambda$CDM cosmology with a primordial matter power spectrum of the form [42]:

$$P_R(k) = 2\pi^2 A_s k^{-3} \left(\frac{k}{k_{\text{pivot}}}\right)^{n_s-1}, \qquad (1.3)$$

such that the linear matter power spectrum is given by [42]:

$$P_L(k, a, \theta) = \frac{8\pi^2 A_s}{25 \Omega_m^2} D_+^2(a) T^2(k) \frac{k_s^n}{H_0^4 k_{\text{pivot}}^{n_s-1}}, \qquad (1.4)$$

where $\theta$ represents the set of cosmological parameters ($\Omega_m$, $A_s$, $n_s$, $H_0$), namely the total matter density parameter, the primordial spectrum normalization, the scalar spectral index and the Hubble constant, respectively; while $T(k)$ denotes the matter transfer function. We also adopt a pivot scale of $k_{\text{pivot}} = 0.05 h^{-1}\text{Mpc}$ following the Planck team's convention [43]. We further assume the prior on the distribution of $\delta_{m,\text{ini}} \equiv (\delta_{m,\text{ini}}(\mathbf{x}_1), \dots, \delta_{m,\text{ini}}(\mathbf{x}_n))$, i.e. ingredient one in the above list, to follow a zero-mean, multivariate Gaussian

$$\mathcal{P}_{\text{prior}}(\delta_{m,\text{ini}|\theta}) = \mathcal{N}(\delta_{m,\text{ini}|\mu=0}, C = \langle\delta_{m,\text{ini}}(\mathbf{x}_i)\delta_{m,\text{ini}}(\mathbf{x}_j)\rangle). \qquad (1.5)$$

In a Fourier-space representation of $\delta_{m,\text{ini}}$, the covariance matrix $C$ is diagonal, and the diagonal is set by the linear matter power spectrum in eq. (1.4), at $a = a_{\text{ini}}$.

For all inferences presented in this work, we proceed by fixing all cosmological parameters to the values used in the simulations, with one exception. Instead of specifying $A_s$, we specify the normalization of the present linear matter power spectrum $\sigma_8$, which is defined as

$$\sigma_8^2 = \sigma^2(R_{\text{TH}} = 8 h^{-1} \text{Mpc}, z=0) = \int_k P_L(k, z=0) W_{\text{TH}}^2(R_{\text{TH}} = 8 h^{-1} \text{Mpc}, k) \qquad (1.6)$$

where $W_{\text{TH}}(R_{\text{TH}}, k)$ is the Fourier transform of the top-hat filter, and fix it to the value in the corresponding simulation. This is due to the convention in the current numerical implementation of BORG. In practice, this subtlety might lead to a percent-level bias in the measured transfer function $T_{\text{infer}}(k)$, especially if the forward inference and the input simulation employ different matter transfer functions $T(k)$. In our case, BORG adopts the Eisenstein-Hu fitting formula [44, 45] while our N-body simulations take the matter transfer functions computed with CAMB or CLASS as inputs (cf. section 3). The implication of this difference is examined in more details in appendix B. In fact, the systematic shift this effect might induce on $T_{\text{infer}}(k)$ is negligible compared to the bias we observed throughout section 4.2.

The numerical implementation of gravity and bias models in the BORG framework discretizes the matter and halo density fields using a cloud-in-cell (CIC) projection onto a grid with cell size $l_{\text{grid}}$ such that

$$\delta_{h,m}^i \equiv \int d^3x W_{\text{CIC}}(l_{\text{grid}}, \mathbf{x}, \mathbf{x}^i) \delta_{h,m}(\mathbf{x}^i - \mathbf{x}) \qquad (1.7)$$

where the superscript $i$ denotes the cell index and $W_{\text{CIC}}(\mathbf{x}, \mathbf{x}^i)$ is the real-space, normalized CIC filter:

$$W_{\text{CIC}}(l_{\text{grid}}, \mathbf{x}, \mathbf{x}^i) = \frac{1}{l_{\text{grid}}} \begin{cases} 1 - |\mathbf{x}^i - \mathbf{x}|/l_{\text{grid}} & |\mathbf{x}^i - \mathbf{x}| < l_{\text{grid}} \\ 0 & \text{otherwise,} \end{cases} \qquad (1.8a)$$
whose Fourier-space representation is given by

$$W_{\text{CIC}}(l_{\text{grid}}, k) = \text{sinc}^2 \left( \frac{k l_{\text{grid}}}{2} \right). \qquad (1.8\text{b})$$

Let us briefly discuss implications of the shape of the smoothing filter here. A sharp-k filter, for example, corresponds to a particular scale cut-off $\Lambda$. This means only modes below the cut-off scale $\Lambda$ enter the analysis, and the higher this cut-off is, the more non-linear modes are included. The CIC filter adopted here, on the other hand, is not localized in Fourier space and thus does not completely remove modes with $k > 1/l_{\text{grid}}$. This might allow non-linear modes to slip into the inference, as further discussed in section 4.2.2 and section 4.2.3.

# 2 Inference setup

Below, for each inspected ingredient in the inference, we introduce and discuss each model or option considered within this study.

## 2.1 Gravity model

We choose the first- and second-order Lagrangian perturbation theories (LPT and 2LPT) [e.g. 46–48] as our gravity models for this comparison. We refer readers to [9] for specific details about the numerical implementation of these models in the BORG framework employed here.

Our two choices are the most common and relevant for forward inference of current galaxy survey data [see, for example 17]. A forward inference powered by any of the two approximations naturally includes optimal Baryon Acoustic Oscillation (BAO) reconstruction (to the degree that the forward model is accurate). Especially, with 2LPT, the final joint posterior should also include sub-leading effects of large-scale perturbation modes For a more detailed discussion on the relation with BAO reconstruction, see section 6 of [18].

## 2.2 Bias model

In this paper, we consider one linear and two non-linear, non-perturbative bias models. All are Eulerian, local-in-matter-density-field (LIMD) [36] bias relations, i.e. the biased tracer density at one location only depends on the matter density at that same location in the Eulerian frame. These models are the ones currently implemented in BORG and recently employed for the analyses in [8, 16, 17, 49, 50].

Our actual observable is the tracer count at a given cell, $n_h^i$. The deterministic or mean-field prediction for this quantity, $n_{h,\text{det}}^i$, is related to the quantity $\delta_{h,\text{det}}^i$ returned by a bias model for the same cell through a nuisance parameter $n_0$ such that

$$n_{h,\text{det}}^i \equiv n_0 (1 + \delta_{h,\text{det}}^i). \qquad (2.1)$$

Note that, in practice, the inferred value of $n_0$ can be different from the true mean number of halos in a grid cell $\langle n_h^i \rangle \equiv N_h/(N_{\text{grid}}^3)$.

### 2.2.1 Linear bias

The linear bias model links the forward-evolved matter density $\delta_{m,\text{fwd}}^i$ to the deterministic halo density $\delta_{h,\text{det}}^i$, via a single bias coefficient $b_1$:

$$\delta_{h,\text{det}}^i = b_1 \delta_{m,\text{fwd}}^i. \qquad (2.2)$$
It is worth noting that the linear bias coefficient $b_1$ in eq. (2.2) is different from the more familiar large-scale linear bias coefficient measured from n-point functions [51, 52]. In fact, all bias coefficients appearing in the field-level approach are moment biases (cf. section 4.2 of [36]), which explicitly depend on both the specific shape of the smoothing filter and the smoothing scale. In particular, on linear and quasi-linear scales where the perturbation theory (PT) for matter and biased tracers is guaranteed to converge to the correct result (when carried out to sufficiently high order), the moment bias $b_1^m$ is related to the large-scale $b_1^l$ as follows [36]

$$b_1^m(l_{\text{grid}}) = b_1^l + \mathcal{O}(R_*^2 \sigma_1^2(l_{\text{grid}})), \qquad (2.3)$$

where $R_*$ and $\sigma_1^2 \equiv \int_k k^2 P_L(k) W_{\text{CIC}}^2(k)$ are the physically relevant scale for the tracer formation process and the first generalized spectral moment, respectively. Given the grid resolutions considered in this study (cf. section 2.4), they yield a leading order correction term of $\mathcal{O}(R_*^2 \sigma_1^2(l_{\text{grid}})) \simeq 0.005 - 0.053$ and $\simeq 0.050 - 1.327$ for the two N-body simulation setups (cf. section 3).

### 2.2.2 Power-law bias

In the power-law bias model, the deterministic halo density in a given cell is related to the matter density at that cell also by only one free parameter, the power-law index $\beta$ [49]

$$\delta_{h, \text{det}}^i = (1 + \delta_{m, \text{fwd}}^i)^\beta - 1. \qquad (2.4)$$

For small perturbations, $\delta_{m, \text{fwd}}^i \ll 1$, eq. (2.4) simplifies to

$$\delta_{h, \text{det}}^i = \beta \delta_{m, \text{fwd}}^i + \frac{\beta(\beta - 1)}{2} (\delta_{m, \text{fwd}}^i)^2 + \mathcal{O}(\delta_{m, \text{fwd}}^3). \qquad (2.5)$$

Thus, for cells where the evolved matter density fluctuations are still small, the role of the power-law index in eq. (2.4) reduces to that of the linear (moment) bias coefficient in eq. (2.2), i.e. $\beta \simeq b_1$ for $\delta_{m, \text{fwd}}^i \ll 1$. As shown in figure 2, however, there is a significant fraction of cells where fluctuations, evolved by either N-body or 2LPT, have become highly non-linear for the typical grid resolution we study in this paper. Thus we expect no convergence of eq. (2.5), but instead significant deviation of $\beta$ from $b_1$ (cf. right panel of figure 11). We will return to this point again in the discussion near the end of section 4.2.3, where we will see that the same issue poses a problem for the validity of the Poisson likelihood (for halos).

### 2.2.3 Broken power-law bias

The broken power-law bias model, as its name indicates, introduces an additional power-law cut-off into the power-law bias relation in eq. (2.4), such that

$$\delta_{h, \text{det}}^i = (1 + \delta_{m, \text{fwd}}^i)^\beta \exp [-\rho(1 + \delta_{m, \text{fwd}}^i)^{-\epsilon}] - 1. \qquad (2.6)$$

This model essentially behaves like the power-law bias model but accounts for the exponential suppression of halo clustering inside voids [53]. It is characterized by three bias parameters: the power-law index $\beta$ and two hyperparameters $(\rho, \epsilon)$ appearing in the exponential cut-off.

We emphasize again that both the power-law and broken power-law bias are LIMD bias models as the right-hand sides (r.h.s) of eq. (2.4) and eq. (2.6) both include only one single operator constructed out of the local matter density field, that is $\delta_{m, \text{fwd}}^i$.
<figure><img src="image_6.png" /><figcaption>Figure 2. Histograms of the forward-evolved matter density fields $\delta_{m,fwd}^i$ in our fiducial simulation box. The solid blue histogram denotes the field evolved by N-body simulation while the dotted-dashed red histogram represents the one evolved by 2LPT. The total number of grid cells is $N_{\text{grid}} = 128^3$.</figcaption></figure>

## 2.3 Likelihood

For this case study, we compare Poisson and Gaussian distributions for the likelihood which captures stochasticity in the actual biased tracer field. These two distributions are the most commonly studied and considered in the literature of halo stochasticity [e.g. 54, 55]. Together, they form the basis for all active likelihoods in BORG [16, 17, 49, 50].

The joint likelihood is a product of single-cell likelihoods evaluated at each individual cell, such that

$$ \ln \mathcal{P} \left( n_h \middle| \{\delta_{m,fwd}\} \right) = \sum_{i=1}^{N_g^3} \ln \mathcal{P}^{(1)} \left( n_h^i \middle| n_{h,\text{det}}^i \right), \quad (2.7) $$

where $\mathcal{P}^{(1)}$ is defined as in eq. (2.8) (for the Poisson likelihood) or eq. (2.9) (for the Gaussian likelihood). The assumption in eq. (2.7) is that the stochasticity is uncorrelated between grid cells.

### 2.3.1 Poisson likelihood

The single-cell Poisson likelihood is given by

$$ \mathcal{P}_{\text{Poisson}}^{(1)} (n_h^i | n_{h,\text{det}}^i) = \frac{\left(n_{h,\text{det}}^i\right)^{n_h^i} e^{-n_{h,\text{det}}^i}}{(n_h^i)!}, \quad (2.8) $$

where $n_{h,\text{det}}^i$ is defined as in eq. (2.1). The Poisson likelihood requires that $n_{h,\text{det}}^i \ge 0$.

### 2.3.2 Gaussian likelihood

The single-cell Gaussian likelihood can be written as

$$ \mathcal{P}_{\text{Gaussian}}^{(1)} (n_h^i | n_{h,\text{det}}^i) = \frac{1}{\sqrt{2\pi}\sigma} \exp \left[ -\frac{(n_h^i - n_{h,\text{det}}^i)^2}{2\sigma^2} \right], \quad (2.9) $$
in which the Gaussian standard deviation $\sigma$ is an additional nuisance parameter and assumed to be universal for all cells. We will return to this assumption in section 4.2.3.

## 2.4 Grid resolution

To examine the effect of the smoothing scale, we run our inferences at five different grid resolutions: $l_{\text{grid}} = 31.3, 20.8, 15.6 \ h^{-1}\text{Mpc}$ with Box-1 (see below) plus its mock counterparts; and $l_{\text{grid}} = 15.6, 7.8, 5.2 \ h^{-1}\text{Mpc}$ with Box-2 plus its mock counterparts. These correspond to a grid of $N_{\text{grid}} = 64^3, 96^3, 128^3$ and $N_{\text{grid}} = 32^3, 64^3, 96^3$, respectively.

## 3 N-body simulation setup

We use one of the GADGET-2 [56] N-body simulations first presented in [57] and later used in [18, 22, 23] as our main input data. We label this fiducial box as Box-1. Box-1 has a box size of $L_{\text{box}} = 2000 \ h^{-1}\text{Mpc}$ and a total number of $N_{\text{part}} = 1536^3$ DM particles, giving it a mass resolution of $M_{\text{part}} = 1.8 \times 10^{11} \ h^{-1}M_{\odot}$. The simulation adopted a flat $\Lambda$CDM cosmology with $\Omega_{\text{m}} = 0.30, n_s = 0.967, h = 0.7, \sigma_8 = 0.850$. The input matter transfer function was computed with CLASS<fnref n="1" /> [58]. Initial conditions for the N-body run were generated at redshift $z_{\text{ini}} = 99$ using the 2LPTic algorithm [59, 60]. DM halos were identified at redshift $z = 0$ as spherical over-densities (SO) [61–63] using the Amiga Halo Finder (AHF) algorithm<fnref n="2" /> [64, 65], where an over-density threshold of 200 times the background matter density is chosen. Note that we only include main halos below $M_h = 10^{14.5}h^{-1}M_{\odot}$ in our analysis. We further divide these halos into three mass bins of $\log_{10} M_h = [13.0, 13.5), [13.5–14.0), [14.0–14.5) \ h^{-1}M_{\odot}$ with the first mass bin being our fiducial sample. Their corresponding comoving number densities are listed in table 1.

Whenever a higher grid resolution inference is required to verify a trend observed in lower resolution inferences, we employ a smaller GADGET-2 simulation box, first presented in [52], as our input data. We label this as Box-2. Box-2 has a box length $L_{\text{box}} = 500 \ h^{-1}\text{Mpc}$ and $N_{\text{part}} = 512^3$ particles, yielding a mass resolution of $M_{\text{part}} = 7.0 \times 10^{10} \ h^{-1}M_{\odot}$. Box-2 also assumed a flat $\Lambda$CDM cosmology, albeit with a slightly different cosmology: $\Omega_{\text{m}} = 0.27, n_s = 0.950, h = 0.7, \sigma_8 = 0.831$. The input matter transfer function was computed with CAMB<fnref n="3" /> [66]. Halos were identified in the same fashion as done in Box-1. We consider only a single mass bin of $\log_{10} M_h = [12.55, 13.55) \ h^{-1}M_{\odot}$ in our analysis of Box-2.

We summarize the halo samples in table 1. It is worth emphasizing that all inferences take halos in their own rest frame as input, i.e. neglecting the RSD effect.

## 4 Results

The phase-correlation coefficient between two given fields $\delta_1, \delta_2$ is given by:

$$r_{12}(k) = \frac{\langle \delta_1 \delta_2 \rangle_k}{\sqrt{\langle \delta_1 \delta_1 \rangle_k \langle \delta_2 \delta_2 \rangle_k}} \quad (4.1)$$

<footnote n="1">¹http://class-code.net.</footnote>

<footnote n="2">²http://popia.ft.uam.es/AHF/.</footnote>

<footnote n="3">³http://camb.info.</footnote>
<table><caption>**Table 1.** The halo samples used in our inferences. Throughout, masses $M \equiv M_{200m}$ are spherical-overdensity masses with respect to 200 times the background matter density.</caption><thead><tr><th>Simulation</th><th>Mass range<br>log<sub>10</sub> M[h<sup>-1</sup>M<sub>☉</sub>]</th><th>N<sub>h</sub></th><th>n&#773;<sub>h</sub> [(h Mpc<sup>-1</sup>)<sup>3</sup>]</th></tr></thead><tbody><tr><td>Box-1</td><td>[13.0, 13.5)</td><td>2672481</td><td>3.34 &times; 10<sup>-4</sup></td></tr><tr><td>Box-1</td><td>[13.5, 14.0)</td><td>882123</td><td>1.10 &times; 10<sup>-4</sup></td></tr><tr><td>Box-1</td><td>[14.0, 14.5)</td><td>232654</td><td>2.91 &times; 10<sup>-5</sup></td></tr><tr><td>Box-2</td><td>[12.55, 13.55)</td><td>127722</td><td>1.02 &times; 10<sup>-3</sup></td></tr></tbody></table>

where $\langle \rangle_k$ denotes the wavenumber average. Hereafter, we define the correlation coefficient eq. (4.1) between the inferred initial modes in the MCMC sample $s'$ and the true ones as

$$r_{\mathrm{ri}}^{s'}(k) = \frac{\langle \delta_{\mathrm{infer}}^{s'} \delta_{\mathrm{input}} \rangle_k}{\sqrt{\langle \langle \delta_{\mathrm{infer}}^s \rangle_s \langle \delta_{\mathrm{infer}}^s \rangle_s \rangle_k \langle \delta_{\mathrm{input}} \delta_{\mathrm{input}} \rangle_k}} \quad (4.2)$$

where $\langle \rangle_s$ denotes the BORG ensemble average, i.e.

$$\langle O \rangle_s \equiv \frac{1}{N_{\mathrm{samples}}} \sum_{s}^{N_{\mathrm{samples}}} O^s, \quad (4.3)$$

where $N_{\mathrm{samples}} = 1401$ is the number of MCMC samples considered for each inference in our analysis. Additionally, we have shortened the notation by writing $\delta_{m,\mathrm{ini}}^{\mathrm{infer}} \equiv \delta_{\mathrm{infer}}$ and $\delta_{m,\mathrm{ini}}^{\mathrm{input}} \equiv \delta_{\mathrm{input}}$. Our estimator of the correlation coefficient for each MCMC chain is then the BORG ensemble average of eq. (4.2):

$$r_{\mathrm{ri}}(k) = \langle r_{\mathrm{ri}}^{s}(k) \rangle_{s} \quad (4.4)$$

whose variance is estimated using the sample variance over BORG samples,

$$\sigma_r^2(k) = \frac{1}{N_{\mathrm{samples}}} \sum_{s'}^{N_{\mathrm{samples}}} \left[ r_{\mathrm{ri}}^{s'}(k) - \langle r_{\mathrm{ri}}^s(k) \rangle_s \right]^2 . \quad (4.5)$$

The inferred phases of the initial conditions are unbiased if $r_{\mathrm{ri}}(k)$ is consistent with 1 within the uncertainty. Note that we use the MCMC ensemble mean $\langle \langle \delta_{\mathrm{infer}} \rangle_s \rangle_k$ instead of the single-sample $\langle \delta_{\mathrm{infer}}^{s'} \rangle_k$ in the denominator of the estimator eq. (4.2). This ensures that, if the inferred initial conditions are drawn from a posterior centered on the true initial conditions, $r_{\mathrm{ri}}(k)$ indeed asymptotes to 1. This would not be the case if one used $\langle \delta_{\mathrm{infer}}^{s'} \rangle_k$ in the denominator, since the inference noise would then always lead to an $r_{\mathrm{ri}}(k)$ below 1.

We use the Gaussian expectation of $r_{\mathrm{ri}}$ as a reference. This large-scale (low-$k$) limit corresponds to the case wherein $\delta_m$, $\delta_h$, and $\varepsilon_h$ (halo stochasticity) fields are all Gaussian fields, where the tracer and matter fields are linearly related. As derived in appendix A, the Gaussian expectation is given by:

$$r_{\mathrm{ri,lin}}(k) = \frac{1}{\sqrt{1 + P_{\varepsilon}(k) / (b_1^2 P_{\mathrm{L}}(k))}}, \quad (4.6)$$
where $P_L(k)$ and $P_\epsilon(k)$ denote the linear matter power spectrum and the halo stochasticity power spectrum, respectively.

Eq. (4.6) clearly shows that there is an expected limit for the inference, set by halo stochasticity and bias, if the inference only has access to information of the linear evolution of LSS. Additionally, eq. (4.6) also suggests that one expects $r_{ri,lin}$ to asymptote to unity on large scales (low-$k$ values) and to decrease at progressively smaller scales (higher-$k$ values). This has nothing to do with non-linear evolution, but just the fact that $P_\epsilon(k)$ asymptotes to a constant on large scales while $P_L(k)$ is a decreasing function of $k$. Note that at very low $k$ (below the range shown here), $r_{ri,lin}$ actually shrinks again due to the turnover in $P_L(k)$. For comparison purpose, since the halo stochasticity power spectrum cannot be precisely determined, we approximate eq. (4.6) by the correlation coefficient between the actual halo and the input initial Fourier modes, $r_{hi}$, which is given by (cf. eq. (4.1)):

$$r_{hi}(k) = \frac{\langle \delta_h \delta_{input} \rangle_k}{\sqrt{\langle \delta_h \delta_h \rangle_k \langle \delta_{input} \delta_{input} \rangle_k}}. \quad (4.7)$$

Below, we will show eq. (4.7) as dotted lines in all of our figures of correlation coefficients to serve as a reference. Note that eq. (4.7) can deviate from eq. (4.6) at high $k$, where the fields are significantly non-Gaussian. However, at low $k$, for whatever weighting applied (on the halos), eq. (4.6) always reduces to eq. (4.7).

The transfer function can be defined for each MCMC sample $s$ as

$$T_{\text{infer}}^{s}(k) = \sqrt{\frac{P_{\text{infer}}^{s}(k)}{P_{\text{true}}(k)}}, \quad (4.8)$$

where, as usual, we define the matter power spectrum as $\langle\delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle_k \equiv (2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k)$. The BORG ensemble mean and variance are then given by:

$$T_{\text{infer}}(k) = \left\langle \sqrt{\frac{P_{\text{infer}}(k)}{P_{\text{true}}(k)}} \right\rangle_s, \quad (4.9)$$

$$\sigma_T^2(k) = \frac{1}{N_{\text{samples}}} \sum_s^{N_{\text{samples}}} [T_{\text{infer}}^s(k) - T_{\text{infer}}(k)]^2. \quad (4.10)$$

The mean inferred amplitudes of the initial conditions at different wavenumbers are unbiased compared to the truth if $T_{\text{infer}}(k)$ is consistent with 1 within the uncertainty.

Below we vary each element in section 2 that could potentially affect the quality of the inference, while keeping the rest of the inference setup fixed. We classify our results into two categories: those from inferences using mock data and those from inferences using N-body DM halo catalogs. For additional clarity, we explicitly specify the details of the setup at the beginning of each section.

## 4.1 Results with mock data

In addition to the N-body halo datasets described in section 3, we generate mock “halo” fields consistently following the same procedure of the forward inference:

1. Draw a Gaussian realization of $\delta_{m,ini}$ from eq. (1.5).

2. Forward evolve $\delta_{m,ini} \rightarrow \delta_{m,fwd}$, assuming a certain gravity model (LPT or 2LPT).
3. Transform $\delta_{m,fwd} \rightarrow \delta_{h,\det}$, assuming one of the bias models presented in section 2.2, plus some fiducial values for the parameters of that model.

4. Draw a random realization of $P(n_h^i | n_{h,\det}^i)$, assuming either the Poisson or Gaussian likelihood presented in section 2.3, plus some fiducial values for $n_0$ (and $\sigma$ in the case of the Gaussian).

We use these mock datasets in section 4.1.1 and section 4.1.2 below, each as input for an inference employing exactly the same choices of gravity model, bias model and likelihood. The goal is to isolate and study the effect of tracer density and grid resolution, respectively, on the inference.

### 4.1.1 Tracer number density

All results in this section are obtained with the following setup:

* 2LPT forward model.

* Power-law bias model.

* Poisson likelihood.

* $l_{\text{grid}} = 15.6 h^{-1} \text{Mpc}$.

We consider three mock “halo” catalogs, as described at the beginning of this section. We emphasize again that, as the biased field is generated from the exact models used in the inference, any difference between the qualities of different inferences should really arise from the tracer number density alone (the same applies for the grid resolution examined in section 4.1.2 below). We consider three cases with tracer (comoving) number densities of $\bar{n}_h = 3.34 \times 10^{-4}$, $1.10 \times 10^{-4}$, $2.91 \times 10^{-5} (h \text{Mpc}^{-1})^3$. In figure 3, we show results for these inferences.

As shown in the left panel, upper plot of figure 3, the phase-correlation coefficient is scale-dependent. At moderately high $k$, the inference from mock data performs slightly better than the corresponding Gaussian expectation, while it perfectly agrees with the latter at low $k \le 0.04 h \text{Mpc}^{-1}$. This trend indicates two things:

* The forward inference does have access to information encoded in the non-linear evolution of LSS and can recover the small-scale phases that have been processed by gravity in the mildly non-linear regime. The inference quality relative to the Gaussian expectation depends very weakly on the tracer density, at least on the scales considered here.

* The Gaussian expectation correctly describes the phase inference on large scales, even in the absence of modeling error.

Indeed, the bottom plot in the left panel of figure 3, which shows the ratios between the phase-correlation coefficients and their corresponding Gaussian expectations, confirms that the Gaussian expectation is an excellent indicator of the cross-correlation between true and inferred initial conditions at low $k$ for a wide range of tracer densities.

The right panel of figure 3 shows the amplitude-bias, or transfer function, in the corresponding inferences. The transfer functions are consistent with 1, which demonstrates the
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Varying tracer density — mock data. Left: the phase-correlation coefficient, estimated by eq. (4.4), between the inferred and the true initial conditions, for different tracer densities, plotted in solid lines. The corresponding shaded bands denote the $1-\sigma$ uncertainties. For reference, the large-scale, Gaussian expectations (cf. eqs. (4.6)–(4.7)) are plotted in dotted lines. Right: the amplitude-bias between inferred and true initial conditions, measured in terms of the transfer function $T_{\text{infer}}$, for the same cases shown on the left. The solid lines again represent the ensemble means (cf. eq. (4.9)) while the corresponding shaded bands designate the $1-\sigma$ uncertainty (cf. eq. (4.10)).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Varying grid resolution — mock data. Top: the phase-correlation coefficient $r(k)$ in inferences with $l_{\text{grid}} = 31.3, 20.8, 15.6 \ h^{-1}\text{Mpc}$ (left) and $l_{\text{grid}} = 15.6, 7.8, 5.2 \ h^{-1}\text{Mpc}$ (right). Dotted lines again represent the Gaussian expectations. Bottom: the transfer function $T_{\text{infer}}(k)$ in the corresponding inferences.</figcaption></figure>
consistency of the inference pipeline. The right panel of figure 3 shows a direct comparison between transfer functions in three inferences, which interestingly appear to be different only at $\lesssim 2\%$ level at the lowest $k$ values. Another important observation is that the uncertainties, introduced by marginalizing over the Fourier phases, is of order of $\lesssim 2\%$ for both the phase-correlation and amplitude-bias. This shows that there is still a lot of information and constraining power even when the phases are sampled, as done in the field-level, forward-modeling approach.

We emphasize again, however, that there is no modeling error since the mock datasets employed in these inferences are sampled from the model assumed in the inference. Additionally, the three mock datasets were generated with the same phases, which explains the correlated fluctuations in the three different mock tracer cases at high $k$ in the right panel of figure 3.

### 4.1.2 Grid resolution

The general setup in this section is as follows:

* 2LPT forward model.

* Power-law bias model.

* Poisson likelihood.

* Mock “halos” with $\bar{n}_h = 3.34 \times 10^{-4} (h \, \text{Mpc}^{-1})^3$ in one $L_{\text{box}} = 2000 \, h^{-1}\text{Mpc}$ (mimicking Box-1) and one $L_{\text{box}} = 500 \, h^{-1}\text{Mpc}$ box (mimicking Box-2).

We run three inferences on each mock data set, at five different grid resolution settings: $l_{\text{grid}} = 31.3, 20.8, 15.6 \, h^{-1}\text{Mpc}$ for the bigger box, and $l_{\text{grid}} = 15.6, 7.8, 5.2 \, h^{-1}\text{Mpc}$ for the smaller one. Increasing the grid resolution certainly increases the amount of input information. However, as discussed in section 1, one must be careful not to include information from highly non-linear scales where our model is no longer valid. Fortunately, with mock data, we can just focus on the effect of the former, i.e. being able to access more small-scale information, on the inference.

The top panels of figure 4 compare the phase-correlation coefficients from six inferences with five different grid resolutions. Higher resolutions are observed to considerably improve $r(k)$ at high-$k$ values. Interestingly, the improvement when going from $l_{\text{grid}} = 7.8 \, h^{-1}\text{Mpc}$ to $l_{\text{grid}} = 5.2 \, h^{-1}\text{Mpc}$ is marginal. This is presumably because we are reaching the shot noise dominated regime, where most cells are empty for this tracer density. We show, in the bottom panels of figure 4, a comparison between transfer functions from the corresponding inferences. Similarly to the case of varying tracer density, the transfer function appears to be rather insensitive to grid resolution. Thus, for such cases where the physical data model perfectly matches the underlying mechanism generating the data, the inference is robustly unbiased.

## 4.2 Results with N-body data

Results shown in figure 3 and figure 4 are obtained with mock data, for which there is no mismatch between the mechanisms generating the data and the models employed in the inference. When N-body halos from Box-1 are instead taken as input data in each of the setup shown above, i.e. all inference settings being the same, we obtain the results shown
in figure 5. These results essentially confirm the trends observed with mock data with one interesting exception, to which we will turn shortly.

First, for the phase-correlation, we observe the same trends, namely the scale-dependence behavior and the sharp decrease after some characteristic wavenumber. Despite possible differences between the physical data model and the actual data, the inferences continue to outperform the Gaussian expectations. Second, for the transfer function, we find no strong dependence on either tracer density or grid resolution. Interestingly, we observe — consistently across all inferences — a clear scale-dependent bias which cannot be explained by the scale-dependent, percent-level bias introduced by different normalization conventions between CLASS (CAMB) and BORG (see appendix B for more details).

The origin of this bias in the measured $T_{\text{infer}}(k)$ is perhaps easier to understand in the linear regime, where

$$ \delta_{m,\text{fwd}}(\mathbf{k}) = D_+(a)\delta_{m,\text{ini}}(\mathbf{k}), \quad (4.11) $$

$$ \delta_h(\mathbf{k}) = b_1\delta_{m,\text{fwd}}(\mathbf{k}) = D_+(a)b_1\delta_{m,\text{ini}}(\mathbf{k}), \quad (4.12) $$

with $D_+(a)$ being the linear growth factor. One immediately sees that the linear bias $b_1$ and the amplitude of initial fluctuations are perfectly degenerate at linear level. A wrong inferred $b_1$ (or $\beta$ in the context of the power-law bias model) thus entails a bias in the amplitude of the inferred modes. To correctly recover both bias parameter(s) and the amplitude, the physical data model must be flexible enough while still keeping sub-grid systematics under rigorous control [24]. In particular, the fact that the bias in $\beta$ only shows up with N-body input further implies that this bias is of physical origin, and thus can only be addressed by improving either the gravity model, the bias model, the likelihood, or any combination of these three ingredients of the physical data model. Below, we repeatedly return to this point in section 4.2.1-section 4.2.3. We further show that it is highly likely a combination of the bias model and likelihood (and to some extent, the gravity model) that is in charge and therefore should be improved.

### 4.2.1 Gravity model

To obtain results in this section, we use the following setup:

* Power-law bias model.

* Poisson likelihood.

* Box-1 (fiducial sample) and Box-2.

* $l_{\text{grid}} = 15.6 h^{-1}\text{Mpc}$ (Box-1) and $l_{\text{grid}} = 5.2 h^{-1}\text{Mpc}$ (Box-2).

We run two inferences using LPT and 2LPT as the forward model. We note that LPT speeds up the HMC sampler by approximately a factor of 3, on average. The warm-up phases, during which the chains approach the stationary target distribution (see appendix D), are very similar in both cases in terms of number of MCMC samples required.

As shown in figure 6, the performances of inferences using LPT mirror those of inferences using 2LPT extremely well at $l_{\text{grid}} = 15.6 h^{-1}\text{Mpc}$, in terms of both metrics. The ensemble means of the inferred initial density fields using LPT and 2LPT are $\gtrsim 99\%$ spatially-correlated. We further verify that this trend is universal among inferences using different bias models and likelihoods. To test the limit of LPT, we increase the resolution to
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Varying tracer density & grid resolution — N-body data. Top: the phase-correlation coefficient in different inferences with varying tracer density (left) and grid resolution (right), when Box-1 is taken as input data of the inference. The Gaussian expectations, plotted in dotted curves, serve as a reference for the performance of each inference. Bottom: the transfer function in the corresponding inferences.</figcaption></figure>

$l_{\text{grid}} = 5.2 \ h^{-1}\text{Mpc}$ using Box-2. In this case we observe a 5-10% difference in performance, favoring 2LPT. Taking into consideration also the computational demand, these results jointly suggest that, for large-volume surveys where the grid resolution will be anyway limited (due to numerical restrictions on number of grid cells), LPT could be a sufficient choice. In fact, the validity of LPT as the forward model for an inference has already been demonstrated in the BOSS-SDSS3 inference in [17] (see section 4.4, especially figure 6). The limit $l_{\text{grid}} = 5.2 \ h^{-1}\text{Mpc}$ — where LPT starts to fall behind 2LPT — however, should be noted.

Going from LPT to 2LPT shifts the transfer functions towards unity, which is the right direction. The shifts are however not nearly enough to yield transfer functions consistent with 1. We therefore argue that the gravity model is not the key to this issue.

### 4.2.2 Bias model

In this section, we adopt the following general setup:

* 2LPT forward model.

* Poisson likelihood.

* Box-1 (fiducial sample) and Box-2.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Varying forward model — N-body data. Top: the phase-correlation coefficient in different inferences with LPT and 2LPT at the grid resolution of $l_{\text{grid}} = 15.6 \ h^{-1}\text{Mpc}$ using Box-1 (left) and $l_{\text{grid}} = 5.2 \ h^{-1}\text{Mpc}$ using Box-2 (right). The dark-gray dotted line shows the Gaussian expectation, which is the same for both cases. Bottom: the transfer function in the corresponding inferences.</figcaption></figure>

• $l_{\text{grid}} = 15.6 \ h^{-1}\text{Mpc}$ (Box-1) and $l_{\text{grid}} = 5.2 \ h^{-1}\text{Mpc}$ (Box-2).

We consider three deterministic bias relations: the linear, the power-law, and the broken power-law models (cf. section 2.2). For clarity, when combining linear bias model (cf. eq. (2.2)) and Poisson likelihood (cf. eq. (2.8)), we need to implement a soft-thresholder which truncates $\delta_{h,\text{det}}$ to ensure that $n_{h,\text{det}} \ge 0$. Details about the thresholding procedure and priors on our bias parameters can be found in appendix C.

Figure 7 shows results of inferences using the three bias models described above, for Box-1 at $l_{\text{grid}} = 15.6 \ h^{-1}\text{Mpc}$ (top panels) and Box-2 at $l_{\text{grid}} = 5.2 \ h^{-1}\text{Mpc}$ (bottom panels). For the former, all bias models considered perform comparably in both metrics: the differences in $r(k)$ are all within the $1-\sigma$ uncertainties; the differences in transfer functions favor the linear bias which exhibits a bias of only roughly 6-10% less than the others. More interestingly, results for the higher-resolution case clearly support this trend: both metrics favor the linear bias model. The differences in the transfer function now increase significantly. Still, none of the cases show a transfer function consistent with unity.

The fact that the differences between bias models are amplified with a higher grid resolution strongly indicates that it is the cells whose fluctuations have become non-linear that dominate the r.h.s. of eq. (2.5), and the expansion itself fails to converge at these resolutions.
<figure><img src="image_11.png" /></figure>

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

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

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

**Figure 7.** Varying bias model — N-body data. Top: the phase-correlation coefficient in different inferences with different bias models described in section 2.2. The same Gaussian expectation for all three cases is illustrated by the dark-gray dotted line. Bottom: the transfer function in the same inferences.

Note that the linear bias model, as implemented with the Poisson likelihood here, is not truly linear. That is, the soft-thresholder, described in appendix C, which ensures that $n_{h,\text{det}}^i$ is always positive (as required for a Poisson variable), naturally deals with voids and extremely low density regions by setting $n_{h,\text{det}}^i = 0$ for those cells. A fairer comparison between the linear and power-law bias models requires, for example, a Gaussian likelihood. We return to this subtlety in figure 9 where it will become clear that the thresholding procedure applied in this section is highly unlikely to be the main reason for the superior performance of the linear bias model.

### 4.2.3 Likelihood

In this section, we adopt the following general setup:

* 2LPT forward model.

* Power-law bias model.

* Box-1 (fiducial sample) and Box-2

* $l_{\text{grid}} = 15.6 \ h^{-1}\text{Mpc}$ (Box-1) and $l_{\text{grid}} = 5.2 \ h^{-1}\text{Mpc}$ (Box-2).
We consider two forms of likelihoods: Poisson and Gaussian distributions (cf. section 2.3). In the case of the Gaussian likelihood, we find that letting the Gaussian variance $\sigma^2$ (cf. eq. (2.9)) free is generally problematic for cases where the average Poisson variance of tracer counts-in-cell is small, that is

$$ \sigma_{\text{Poisson}}^2 = \left\langle \left( n_h^i - n_{h,\text{det}}^i \right)^2 \right\rangle_{\text{Poisson}} = N_h/N_{\text{grid}}^3 \ll 1. \quad (4.13) $$

Typically for such cases, we find that the variance parameter $\sigma^2$ collapses to a very small value, i.e. $\sigma^2 \ll N_h/N_{\text{grid}}^3$, and the chain only explores a narrow region around this value in posterior space. In order avoid this issue, we fix $\sigma^2 = N_h/N_{\text{grid}}^3$ in all inferences employing a Gaussian likelihood. It would be worth trying to alleviate this issue in the future by, for example, relaxing the universal variance assumption in eq. (2.9) and allowing for the variance to be density-dependent, as discussed in [24].<fnref n="4" /> This way, the likelihood could essentially downweight regions with low signal-to-noise while upweighting those with higher signal-to-noise.<fnref n="5" />

As shown in the top panels of figure 8, the phase-correlation metric marginally favors the Poisson likelihood. Interestingly, however, the transfer function strongly prefers the Gaussian case. Perhaps it is worth keeping in mind that the amplitude-bias is robustly observed across different halo mass ranges, grid resolutions, forward models and bias models. We therefore conclude the Poisson likelihood is incompatible with halo stochasticity. This is perhaps not entirely surprising as previous literature has already pointed out that halo stochasticity can be either sub-Poisson, due to halo exclusion effect [55, 69, 70], or super-Poisson in high-density regions [54]. Intuitively, we only expect halo stochasticity to asymptote the Poisson limit at very high wavenumber $k$, however, precisely in this regime, our bias expansions no longer converge due to non-linearities in the matter density field (cf. figure 2).

To summarize, our results demonstrate that using a Poisson likelihood to model halo stochasticity typically leads to a biased transfer function, i.e. a bias in the amplitudes of the inferred modes. This is an important implication for future cosmological inference using the field-level, forward-modeling approach.

Note, however, that (main) halo and galaxy stochasticities might significantly differ: a likelihood, or a bias model for that matter, optimized on main halos might not be compatible with galaxy stochasticity or bias. Rather than a recommendation to always use the Gaussian rather than the Poisson likelihood, or the linear bias than a more sophisticated one, for example the broken power-law bias model, our results firmly suggest that the likelihood and the bias model should be rigorously derived from, for example, a perturbative approach, such that sub-grid systematics can be kept under control [e.g. 18–20, 22–24].

We now revisit the case of the linear bias model. Figure 9 shows a direct comparison between the Poisson and the Gaussian likelihoods, when combined with the linear bias model, at the highest resolution we have probed. The Gaussian likelihood does not involve any thresholding procedure. Interestingly, each metric now favors a different setup: the Poisson case seems to outperform the Gaussian slightly above $k \simeq 0.2 h \, \text{Mpc}^{-1}$ in terms of phase-correlation in the left panel of figure 9, whereas the Gaussian case results in a less biased transfer function up to $k \simeq 0.2 h \, \text{Mpc}^{-1}$. The fact that the bias in $T_{\text{infer}}(k)$ is further reduced with the Gaussian likelihood — and without the soft-thresholder (cf. appendix C)— suggests that the thresholding procedure plays no significant role in our results previously

<footnote n="4">⁴This is the case for the Poisson likelihood which specifically has a variance of $[\sigma^2]^i = n_{h,\text{det}}^i \equiv n_0(1+\delta_{h,\text{det}}^i)$.</footnote>

<footnote n="5">⁵The 2-point function analogy of which is the marked power spectrum [e.g. 67, 68].</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 8.** Varying likelihood — N-body data. Top: the phase-correlation coefficient in different inferences with Poisson and Gaussian likelihoods as described in section 2.3. The common Gaussian expectation is shown in dark-gray dotted line. Bottom: the transfer function in the corresponding inferences.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 9.** Revisiting the linear bias model — N-body data. Left: the phase-correlation coefficient in two inferences with Poisson and Gaussian likelihoods using the linear bias model. In the Gaussian case, we do not apply any thresholding procedure (cf. appendix C) as it is unnecessary. We show the linear Gaussian expectation in dark-gray dotted line. Right: the transfer function in the corresponding inferences.</figcaption></figure>
shown in figure 7. We additionally note that the case shown in figure 9 exhibits the most
significant difference in the phase correlations across all of our comparisons, going up to 5%
at $k = 0.2 h \text{ Mpc}^{-1}$.

# 5 Discussion

In this paper, we have examined the forward modeling approach of biased-tracer clustering in terms of inferring the initial conditions in the observed cosmological volume. Specifically, employing the BORG framework and N-body simulations, we inspect how different ingredients in the physical data model affect the performance of the inference of the initial conditions. These include: tracer density, grid resolution, gravity model, bias model and likelihood. We quantify the performance of each inference by two metrics:

* The cross correlation between inferred and input initial modes, quantified by the phase-correlation coefficient;

* The amplitude bias between the inferred and input amplitudes of these initial modes, quantified by the transfer function.

## Tracer density and grid resolution

First, in section 4.1.1, using mock tracers, we observe that — in the absence of modeling error — enhancing tracer density generally improves the inference, even at scales where the clustering signal dominates over shot noise. Similarly, increasing grid resolution increases the small-scale information the inference can access and consequently boosts the quality of the inference, as found in section 4.1.2. Both effects described above can be perfectly captured by the linear Gaussian expectation we derived in eq. (4.6). At low *k*, the Gaussian model is an *unbiased estimator* of the cross correlation between the inferred phases and the underlying truth. At high *k*, this relationship is biased but depends on neither tracer density nor grid resolution.

## Gravity model

Second, in section 4.2, using DM halos from N-body simulations as physical tracers, we have investigated the robustness of the inference when modeling error is present. In particular, the choice of gravity model considered in section 4.2.1, be it LPT and 2LPT, does not significantly affect the performance of the inference up to a moderate grid resolution of $l_{\text{grid}} = 15.6 \ h^{-1}\text{Mpc}$. This has a strong implication for cosmological inference using forward modeling approach with data from future large-volume surveys. We note however that 2LPT outperforms LPT at higher grid resolutions.

## Bias model and likelihood

Last, but perhaps most interesting, we observe, for all inferences with N-body data, that the Poisson likelihood results in a significant bias on the amplitude of the inferred modes. Further investigation is needed to determine exactly the origin of this effect, though as discussed in section 4.2.2 and section 4.2.3, our results suggest that this is an indication of un-controlled systematics from highly non-linear fluctuations (cf. figure 2) and their couplings (see also results in [24], who follow an approach based on effective field theory (EFT) which adopts a sharp-$k$ filter in place of the CIC filter employed here).

Let us summarize our findings below:
* The phase-correlation is sensitive to a change in any ingredient. Fortunately, the shifts are generally ≲ 5% in any direction and can be well predicted by the Gaussian expectation in all cases. This implies that analyses that cross-correlate the inference with other datasets [e.g. 71], will be robust to details of the physical data model.

* The transfer function is found to be only susceptible to the choices of bias model and likelihood. However, the differences between different likelihoods and bias models are significant. At this point, unfortunately, it is unclear how to predict these from first principles.

# 6 Conclusions

To our knowledge, our results represent the most stringent test for Bayesian forward modeling and inference of biased tracer clustering, within the BORG framework, so far. Our input data are main halos identified in N-body simulations, unlike the mock tracers generated from the same physical data model employed by previous analyses [see, e.g. 9, 50]. More realistic and rigorous tests are still awaiting ahead in order to achieve the goal of unbiased inference of initial conditions and cosmological parameters using this approach. Specifically, extending the input biased tracers to include galaxies — either populated inside DM halos using some halo occupation distribution model, or directly identified in hydrodynamic simulations — would be such a test. An even more realistic case must also consider RSD, survey geometries and selection effects applied to the galaxy field. We defer such studies to future work as well.

Clearly, a transfer function $T_{\text{infer}}(k)$ with a scale-dependent bias up to $\sim 30\%$ has a significant implication for cosmological inference, especially when future analyses are aiming to reach percent-level constraints on cosmological parameters. In particular, one would expect biases in not only $A_s$ but also $\Omega_m$, $n_s$ and probably the distance scale as well. We thus view this issue as a call for attention to studies and developments of accurate bias models and likelihoods. As already pointed out above, BORG and the work in this paper present a physically-motivated and technically-ready statistical framework for such studies.

Specifically, an EFT-based approach such as the one presented in [18, 22, 24] is certainly preferable. In fact, [24] showed that for more conservative scale cut-offs, one can recover an unbiased primordial power spectrum normalization within percent-level for the same fiducial halo sample studied in this work. The analysis in [24] however fixed the phases. It would be interesting to investigate how the EFT approach performs, specifically in terms of the level of systematic bias, when the phases are inferred from the data and marginalized over, as done in this paper. We leave this investigation for an upcoming work.

# Acknowledgments

We would like to thank Eiichiro Komatsu, Alex Barreira, Giovani Cabass, Ariel Sanchez as well as all Aquila Consortium members for helpful discussions. We especially thank Titouan Lazeyras for sharing the access to the N-body simulations. We extend our thanks to Elisa Ferreira and Giovanni Aricò for useful comments that helped improving the plots. MN and FS acknowledge support from the Starting Grant (ERC-2015-STG 678652) “GrInflaGal” of the European Research Council. GL acknowledges financial support from the ILP LABEX (under reference ANR-10-LABX-63) which is financed by French state funds managed by the
ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02. This work was supported by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. This work is done within the Aquila Consortium.<fnref n="6" />

## A Gaussian expectation of the Fourier phase-correlation

In this appendix we derive the Gaussian expectation of the phase-correlation coefficient, as shown in eq. (4.6).

As discussed in section 4, within this limit, we can assume that all fields (including the noise field) are Gaussian, and assume a linear bias relation of the following form:

$$ \delta_h(\mathbf{k}) = b_1 \delta_{m,\text{lin}}^{\text{input}}(\mathbf{k}) + \varepsilon_h(\mathbf{k}), \qquad (A.1) $$

where $\delta_{m,\text{lin}}^{\text{input}}$ here denotes the input matter density field linearly extrapolated to redshift zero. Since the halo stochasticity $\varepsilon_h$ has zero mean and is assumed to be uncorrelated with the matter density field, the optimal estimator for the reconstructed density field is simply given by

$$ \delta_{m,\text{lin}}^{\text{infer}}(\mathbf{k}) = \frac{1}{b_1} \delta_h(\mathbf{k}). \qquad (A.2) $$

This together with eq. (A.1) implies that

$$ \langle \delta_{m,\text{lin}}^{\text{infer}}(\mathbf{k}) \delta_{m,\text{lin}}^{\text{input}}(\mathbf{k}) \rangle_k = \langle \delta_{m,\text{lin}}^{\text{input}}(\mathbf{k}) \delta_{m,\text{lin}}^{\text{input}}(\mathbf{k}) \rangle_k = P_L(k), \qquad (A.3) $$

where $P_L$ denotes the linear matter power spectrum at $z=0$. In addition, the power spectrum of the inferred modes can be expressed as (cf. eq. (A.2) and eq. (A.1))

$$ \begin{aligned} \langle \delta_{m,\text{lin}}^{\text{infer}}(\mathbf{k}) \delta_{m,\text{lin}}^{\text{infer}}(\mathbf{k}) \rangle_k &= \frac{1}{b_1^2} \langle \delta_h(\mathbf{k}) \delta_h(\mathbf{k}) \rangle_k \qquad (A.4) \\ &= P_L(k) + \frac{1}{b_1^2} P_\varepsilon(k). \end{aligned} $$

Combining eq. (A.3) with eq. (A.4) then yields the desired cross-correlation coefficient between $\delta_{m,\text{lin}}^{\text{infer}}$ and $\delta_{m,\text{lin}}^{\text{input}}$:

$$ r_{\text{ri,lin}}(k) = \frac{1}{\sqrt{1 + P_{\varepsilon}(k) / (b_1^2 P_{\text{L}}(k))}}. \qquad (A.5) $$

## B Residual bias due to normalization using $\sigma_8$

Our N-body simulations use input matter transfer functions computed by Boltzmann codes CLASS and CAMB, respectively. On the other hand, for performance reasons, BORG employs the physically-motivated Eisenstein-Hu fitting function [44, 45] to compute the same quantity. In either cases, the matter transfer functions were outputted at redshift zero, then back-scaled to the initial redshifts in the simulation and the forward model. Further, BORG cosmological setting requests one to specify the normalization of the linear matter power spectrum at redshift zero $\sigma_8$ instead of the normalization of the primordial power spectrum $A_s$. This convention, plus the difference between matter transfer functions computed by CLASS or

<footnote n="6">⁶https://aquila-consortium.org.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 10.** Ratios between the transfer functions at redshift zero, computed by CLASS (red) or CAMB (blue) and that determined from Eisenstein-Hu fitting formula.</figcaption></figure>

CAMB and those obtained from the Eisenstein-Hu fitting formula, together will result in a bias between the inferred and input initial conditions, and hence a deviation of $T_{\text{rec}}(k)$ from 1. In this appendix, we examine this effect by measuring the ratio between these transfer functions. Specifically, figure 10 shows that one can expect a scale-dependent bias in $T_{\text{rec}}$ up to 1.7% on the scales of interest in this paper. This clearly cannot explain the bias observed in section 4.2.

## C Priors and thresholding

For all inferences presented in this paper, we impose uniform priors and positivity constraints on all bias parameters, including the nuisance parameter $n_0$:

$$ \mathcal{P}_{\text{prior}}(n_0, \{b_O\}) = \Theta(n_0) \prod_O \Theta(b_O), \qquad (\text{C.1}) $$

in which $\Theta(x)$ is the Heaviside step function.

The Poisson likelihood in eq. (2.8) requires a non-negative $n_{h,\text{det}}$; in order to satisfy this condition in the particular case of the linear bias model, we additionally implement a soft-thresholder of the form:

$$ n_{h,\text{det}} = n_0 \begin{cases} (1 + \delta_{h,\text{det}}) & \delta_{h,\text{det}} > -0.9 \\ th(1 + \delta_{h,\text{det}}, c) & \text{otherwise,} \end{cases} \qquad (\text{C.2}) $$

with

$$ th(1 + \delta_{h,\text{det}}, c) = \ln(1 + \exp[c(1 + \delta_{h,\text{det}})] / c) \qquad (\text{C.3}) $$

being the softplus function [72]. We choose $c = 5$ for our implementation. We emphasize again that this thresholding procedure is only needed for, and thus applied to, inferences using both linear bias model and Poisson likelihood.
<figure><img src="image_18.png" /><figcaption>**Figure 11.** *Left:* one example of tracing the power spectrum posterior during and after warm-up phase. For visual clarity, we only plot 1-in-every-200 MCMC samples. The color scale indicates the MCMC identifier. The dotted line denotes the input power spectrum for the given cosmology. The low-k bias observed here corresponds to the bias in the transfer function seen in the bottom left panel of figure 5 (red line there). *Right:* trace plot of the bias parameters in the same inference. Again, each point in this plot represents 1-in-every-200 MCMC samples. Further, samples whose identifier is smaller than 100 were removed to enhance visibility of the y-axis.</figcaption></figure>

## D MCMC burn-in phase and chain thinning

This appendix describes in detail how we verify that our MCMC chains have indeed approached the target typical sets and how we select MCMC samples in our analysis.

All the Markov chains considered in this paper were initialized in an over-dispersed state, presumably remotely far from the stationary distribution region in the posterior space. We refer to this initial phase as the *warm-up* phase. For each chain, we ensure that samples included in our analysis do not come from this phase by carefully tracing the posterior power spectrum and bias parameters. In figure 11, we show two such examples. Both trace plots confirm that the chain has reached a stationary distribution after around 2000 samples.

To achieve a reasonable effective sample size while avoiding wasting disk storage, we only store and analyze 1-in-every-20 MCMC sample.

## References

[1] H. Gil-Marín, W.J. Percival, L. Verde, J.R. Brownstein, C.-H. Chuang, F.-S. Kitaura et al., *The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the power spectrum and bispectrum of the DR12 BOSS galaxies*, Mon. Not. Roy. Astron. Soc. 465 (2017) 1757 [arXiv:1606.00439].

[2] BOSS collaboration, *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample*, Mon. Not. Roy. Astron. Soc. 470 (2017) 2617 [arXiv:1607.03155].

[3] M.M. Ivanov, M. Simonović and M. Zaldarriaga, *Cosmological Parameters from the BOSS Galaxy Power Spectrum*, JCAP 05 (2020) 042 [arXiv:1909.05277].

[4] DES collaboration, *Dark Energy Survey Year 1 Results: Cosmological constraints from cluster abundances and weak lensing*, Phys. Rev. D 102 (2020) 023509 [arXiv:2002.11124].
[5] U. Frisch, S. Matarrese, R. Mohayaee and A. Sobolevski, *A Reconstruction of the initial conditions of the universe by optimal mass transportation*, *Nature* **417** (2002) 260 [astro-ph/0109483].

[6] N. Padmanabhan, X. Xu, D.J. Eisenstein, R. Scalzo, A.J. Cuesta, K.T. Mehta et al., *A 2 per cent distance to z = 0.35 by reconstructing baryon acoustic oscillations — I. Methods and application to the Sloan Digital Sky Survey*, *Mon. Not. Roy. Astron. Soc.* **427** (2012) 2132 [arXiv:1202.0090].

[7] M. Schmittfull, T. Baldauf and M. Zaldarriaga, *Iterative initial condition reconstruction*, *Phys. Rev. D* **96** (2017) 023505 [arXiv:1704.06634].

[8] Y. Wang, B. Li and M. Cautun, *Iterative removal of redshift space distortions from galaxy clustering*, *Mon. Not. Roy. Astron. Soc.* **497** (2020) 3451 [arXiv:1912.03392].

[9] J. Jasche and B.D. Wandelt, *Bayesian physical reconstruction of initial conditions from large scale structure surveys*, *Mon. Not. Roy. Astron. Soc.* **432** (2013) 894 [arXiv:1203.3639].

[10] H. Wang, H.J. Mo, X. Yang, Y.P. Jing and W.P. Lin, *ELUCID — Exploring the Local Universe with reConstructed Initial Density field I: Hamiltonian Markov Chain Monte Carlo Method with Particle Mesh Dynamics*, *Astrophys. J.* **794** (2014) 94 [arXiv:1407.3451].

[11] J. Jasche, F. Leclercq and B.D. Wandelt, *Past and present cosmic structure in the SDSS DR7 main sample*, *JCAP* **01** (2015) 036 [arXiv:1409.6308].

[12] M. Ata, F.-S. Kitaura and V. Müller, *Bayesian inference of cosmic density fields from non-linear, scale-dependent, and stochastic biased tracers*, *Mon. Not. Roy. Astron. Soc.* **446** (2015) 4250 [arXiv:1408.2566].

[13] H. Wang, H.J. Mo, X. Yang, Y. Zhang, J. Shi, Y.P. Jing et al., *ELUCID — Exploring the Local Universe with reConstructed Initial Density field III: Constrained Simulation in the SDSS Volume*, *Astrophys. J.* **831** (2016) 164 [arXiv:1608.01763].

[14] BOSS collaboration, *The Clustering of Galaxies in the Completed SDSS-III Baryon Oscillation Spectroscopic Survey: Cosmic Flows and Cosmic Web from Luminous Red Galaxies*, *Mon. Not. Roy. Astron. Soc.* **467** (2017) 3993 [arXiv:1605.09745].

[15] C. Modi, Y. Feng and U. Seljak, *Cosmological Reconstruction From Galaxy Light: Neural Network Based Light-Matter Connection*, *JCAP* **10** (2018) 028 [arXiv:1805.02247].

[16] J. Jasche and G. Lavaux, *Physical Bayesian modelling of the non-linear matter distribution: new insights into the Nearby Universe*, *Astron. Astrophys.* **625** (2019) A64 [arXiv:1806.11117].

[17] G. Lavaux, J. Jasche and F. Leclercq, *Systematic-free inference of the cosmic matter density field from SDSS3-BOSS data*, *arXiv:1909.06396*.

[18] F. Schmidt, F. Elsner, J. Jasche, N.M. Nguyen and G. Lavaux, *A rigorous EFT-based forward model for large-scale structure*, *JCAP* **01** (2019) 042 [arXiv:1808.02002].

[19] G. Cabass and F. Schmidt, *The EFT Likelihood for Large-Scale Structure*, *JCAP* **04** (2020) 042 [arXiv:1909.04022].

[20] G. Cabass and F. Schmidt, *The Likelihood for LSS: Stochasticity of Bias Coefficients at All Orders*, *JCAP* **07** (2020) 051 [arXiv:2004.00617].

[21] G. Cabass, *The EFT Likelihood for Large-Scale Structure in Redshift Space*, *JCAP* **01** (2021) 067 [arXiv:2007.14988].

[22] F. Elsner, F. Schmidt, J. Jasche, G. Lavaux and N.-M. Nguyen, *Cosmology inference from a biased density field using the EFT-based likelihood*, *JCAP* **01** (2020) 029 [arXiv:1906.07143].

[23] F. Schmidt, G. Cabass, J. Jasche and G. Lavaux, *Unbiased Cosmology Inference from Biased Tracers using the EFT Likelihood*, *JCAP* **11** (2020) 008 [arXiv:2004.06707].
[24] F. Schmidt, *Sigma-Eight at the Percent Level: The EFT Likelihood in Real Space*, arXiv:2009.14176.

[25] M. Takada and B. Jain, *The three-point correlation function in cosmology*, Mon. Not. Roy. Astron. Soc. **340** (2003) 580 [astro-ph/0209167].

[26] E. Sefusatti, M. Crocce, S. Pueblas and R. Scoccimarro, *Cosmology and the Bispectrum*, Phys. Rev. D **74** (2006) 023522 [astro-ph/0604505].

[27] A. Lazanu, *Bispectrum Modelling in Large Scale Structure — A Three Shape Model*, 9, 2017 [arXiv:1709.09425].

[28] I. Hashimoto, Y. Rasera and A. Taruya, *Precision cosmology with redshift-space bispectrum: a perturbation theory based model at one-loop order*, Phys. Rev. D **96** (2017) 043526 [arXiv:1705.02574].

[29] N.S. Sugiyama, S. Saito, F. Beutler and H.-J. Seo, *Perturbation theory approach to predict the covariance matrices of the galaxy power spectrum and bispectrum in redshift space*, Mon. Not. Roy. Astron. Soc. **497** (2020) 1684 [arXiv:1908.06234].

[30] Z. Slepian et al., *The large-scale three-point correlation function of the SDSS BOSS DR12 CMASS galaxies*, Mon. Not. Roy. Astron. Soc. **468** (2017) 1070 [arXiv:1512.02231].

[31] J. Carron, *On the incompleteness of the moment and correlation function hierarchy as probes of the lognormal field*, Astrophys. J. **738** (2011) 86 [arXiv:1105.4467].

[32] J.M. Bardeen, J.R. Bond, N. Kaiser and A.S. Szalay, *The Statistics of Peaks of Gaussian Random Fields*, Astrophys. J. **304** (1986) 15.

[33] S. Weinberg, *Cosmology*, Oxford University Press (2008).

[34] WMAP collaboration, *Wilkinson Microwave Anisotropy Probe (WMAP) three year results: implications for cosmology*, Astrophys. J. Suppl. **170** (2007) 377 [astro-ph/0603449].

[35] PLANCK collaboration, *Planck 2018 results. IX. Constraints on primordial non-Gaussianity*, Astron. Astrophys. **641** (2020) A9 [arXiv:1905.05697].

[36] V. Desjacques, D. Jeong and F. Schmidt, *Large-Scale Galaxy Bias*, Phys. Rept. **733** (2018) 1 [arXiv:1611.09787].

[37] J. Birkin, B. Li, M. Cautun and Y. Shi, *Reconstructing the baryon acoustic oscillations using biased tracers*, Mon. Not. Roy. Astron. Soc. **483** (2019) 5267 [arXiv:1809.08135].

[38] S. Duane, A.D. Kennedy, B.J. Pendleton and D. Roweth, *Hybrid Monte Carlo*, Phys. Lett. B **195** (1987) 216.

[39] R.M. Neal, *MCMC using Hamiltonian dynamics*, arXiv:1206.1901 (2012).

[40] J. Jasche and F.S. Kitaura, *Fast Hamiltonian sampling for large scale structure inference*, Mon. Not. Roy. Astron. Soc. **407** (2010) 29 [arXiv:0911.2496].

[41] M.J. Betancourt, S. Byrne, S. Livingstone and M. Girolami, *The Geometric Foundations of Hamiltonian Monte Carlo*, arXiv:1410.5110 (2014).

[42] S. Dodelson and F. Schmidt, *Modern Cosmology*, second edition, Academic Press (2021) [https://doi.org/10.1016/C2017-0-01943-2].

[43] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. **641** (2020) A6 [arXiv:1807.06209].

[44] D.J. Eisenstein and W. Hu, *Baryonic features in the matter transfer function*, Astrophys. J. **496** (1998) 605 [astro-ph/9709112].

[45] D.J. Eisenstein and W. Hu, *Power spectra for cold dark matter and its variants*, Astrophys. J. **511** (1999) 5 [astro-ph/9710252].
[46] F. Moutarde, J.M. Alimi, F.R. Bouchet, R. Pellat and A. Ramani, *Precollapse Scale Invariance in Gravitational Instability*, Astrophys. J. 382 (1991) 377.

[47] F.R. Bouchet, R. Juszkiewicz, S. Colombi and R. Pellat, *Weakly nonlinear gravitational instability for arbitrary Omega*, Astrophys. J. Lett. 394 (1992) L5.

[48] F.R. Bouchet, S. Colombi, E. Hivon and R. Juszkiewicz, *Perturbative Lagrangian approach to gravitational instability*, Astron. Astrophys. 296 (1995) 575 [astro-ph/9406013].

[49] G. Lavaux and J. Jasche, *Unmasking the Masked Universe: the 2M++ catalogue through Bayesian eyes*, Mon. Not. Roy. Astron. Soc. 455 (2016) 3169 [arXiv:1509.05040].

[50] D.K. Ramanah, G. Lavaux, J. Jasche and B.D. Wandelt, *Cosmological inference from Bayesian forward modelling of deep galaxy redshift surveys*, Astron. Astrophys. 621 (2019) A69 [arXiv:1808.07496].

[51] J.L. Tinker, A.V. Kravtsov, A. Klypin, K. Abazajian, M.S. Warren, G. Yepes et al., *Toward a halo mass function for precision cosmology: The Limits of universality*, Astrophys. J. 688 (2008) 709 [arXiv:0803.2706].

[52] T. Lazeyras, C. Wagner, T. Baldauf and F. Schmidt, *Precision measurement of the local bias of dark matter halos*, JCAP 02 (2016) 018 [arXiv:1511.01096].

[53] M.C. Neyrinck, M.A. Aragon-Calvo, D. Jeong and X. Wang, *A halo bias function measured deeply into voids without stochasticity*, Mon. Not. Roy. Astron. Soc. 441 (2014) 646 [arXiv:1309.6641].

[54] R. Casas-Miranda, H.J. Mo, R.K. Sheth and G. Boerner, *On the Distribution of Haloes, Galaxies and Mass*, Mon. Not. Roy. Astron. Soc. 333 (2002) 730 [astro-ph/0105008].

[55] T. Baldauf, U. Seljak, R.E. Smith, N. Hamaus and V. Desjacques, *Halo stochasticity from exclusion and nonlinear clustering*, Phys. Rev. D 88 (2013) 083507 [arXiv:1305.2917].

[56] V. Springel, *The Cosmological simulation code GADGET-2*, Mon. Not. Roy. Astron. Soc. 364 (2005) 1105 [astro-ph/0505010].

[57] M. Biagetti, T. Lazeyras, T. Baldauf, V. Desjacques and F. Schmidt, *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].

[58] J. Lesgourgues, *The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview*, arXiv:1104.2932.

[59] M. Crocce, S. Pueblas and R. Scoccimarro, *Transients from Initial Conditions in Cosmological Simulations*, Mon. Not. Roy. Astron. Soc. 373 (2006) 369 [astro-ph/0606505].

[60] R. Scoccimarro, L. Hui, M. Manera and K.C. Chan, *Large-scale Bias and Efficient Generation of Initial Conditions for Non-Local Primordial Non-Gaussianity*, Phys. Rev. D 85 (2012) 083002 [arXiv:1108.5512].

[61] J.E. Gunn and J.R. Gott, III, *On the Infall of Matter into Clusters of Galaxies and Some Effects on Their Evolution*, Astrophys. J. 176 (1972) 1.

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

[63] A. Kravtsov and S. Borgani, *Formation of Galaxy Clusters*, Ann. Rev. Astron. Astrophys. 50 (2012) 353 [arXiv:1205.5556].

[64] S.P.D. Gill, A. Knebe and B.K. Gibson, *The Evolution substructure 1: A New identification method*, Mon. Not. Roy. Astron. Soc. 351 (2004) 399 [astro-ph/0404258].

[65] S.R. Knollmann and A. Knebe, *Ahf: Amiga's Halo Finder*, Astrophys. J. Suppl. 182 (2009) 608 [arXiv:0904.3662].
[66] A. Lewis, A. Challinor and A. Lasenby, *Efficient computation of CMB anisotropies in closed FRW models*, *Astrophys. J.* **538** (2000) 473 [astro-ph/9911177].

[67] E. Massara, F. Villaescusa-Navarro, S. Ho, N. Dalal and D.N. Spergel, *Using the Marked Power Spectrum to Detect the Signature of Neutrinos in Large-Scale Structure*, *Phys. Rev. Lett.* **126** (2021) 011301 [arXiv:2001.11024].

[68] O.H.E. Philcox, E. Massara and D.N. Spergel, *What does the marked power spectrum measure? Insights from perturbation theory*, *Phys. Rev. D* **102** (2020) 043516 [arXiv:2006.10055].

[69] R.E. Smith, R. Scoccimarro and R.K. Sheth, *The Scale Dependence of Halo and Galaxy Bias: Effects in Real Space*, *Phys. Rev. D* **75** (2007) 063512 [astro-ph/0609547].

[70] N. Hamaus, U. Seljak, V. Desjacques, R.E. Smith and T. Baldauf, *Minimizing the Stochasticity of Halos in Large-Scale Structure Surveys*, *Phys. Rev. D* **82** (2010) 043515 [arXiv:1004.5377].

[71] N.-M. Nguyen, J. Jasche, G. Lavaux and F. Schmidt, *Taking measurements of the kinematic Sunyaev-Zel'dovich effect forward: including uncertainties from velocity reconstruction with forward modeling*, *JCAP* **12** (2020) 011 [arXiv:2007.13721].

[72] X. Glorot, A. Bordes and Y. Bengio, *Deep sparse rectifier neural networks*, in *Proceedings of Machine Learning Research*, Fort Lauderdale, FL, U.S.A., 11–13 Apr, 2011, vol. 15, pp. 315–323 [http://proceedings.mlr.press/v15/glorot11a.html].