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

HAL
open science

# Bayesian field-level inference of primordial non-Gaussianity using next-generation galaxy surveys

Adam Andrews, Jens Jasche, Guilhem Lavaux, Fabian Schmidt

► To cite this version:

Adam Andrews, Jens Jasche, Guilhem Lavaux, Fabian Schmidt. Bayesian field-level inference of primordial non-Gaussianity using next-generation galaxy surveys. Monthly Notices of the Royal Astronomical Society, 2023, 520 (4), pp.5746-5763. 10.1093/mnras/stad432. hal-03625865

HAL Id: hal-03625865

https://hal.science/hal-03625865

Submitted on 15 Feb 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bayesian field-level inference of primordial non-Gaussianity using next-generation galaxy surveys

Adam Andrews<sup>1,★</sup> Jens Jasche,<sup>1,2</sup> Guilhem Lavaux<sup>1,2</sup> and Fabian Schmidt<sup>3</sup>

Accepted 2023 January 19. Received 2022 December 16; in original form 2022 March 16

## ABSTRACT

Detecting and measuring a non-Gaussian signature of primordial origin in the density field is a major science goal of next-generation galaxy surveys. The signal will permit us to determine primordial-physics processes and constrain models of cosmic inflation. While traditional approaches use a limited set of statistical summaries of the galaxy distribution to constrain primordial non-Gaussianity, we present a field-level approach by Bayesian forward modelling the entire three-dimensional galaxy survey. Since our method includes the entire cosmic field in the analysis, it can naturally and fully self-consistently exploit all available information in the large-scale structure, to extract information on the local non-Gaussianity parameter, $f_{nl}$. Examples include higher order statistics through correlation functions, peculiar velocity fields through redshift-space distortions, and scale-dependent galaxy bias. To illustrate the feasibility of field-level primordial non-Gaussianity inference, we present our approach using a first-order Lagrangian perturbation theory model, approximating structure growth at sufficiently large scales. We demonstrate the performance of our approach through various tests with self-consistent mock galaxy data emulating relevant features of the SDSS-III/BOSS-like survey, and additional tests with a *Stage IV* mock data set. These tests reveal that the method infers unbiased values of $f_{nl}$ by accurately handling survey geometries, noise, and unknown galaxy biases. We demonstrate that our method can achieve constraints of $\sigma_{f_{nl}} \approx 8.78$ for SDSS-III/BOSS-like data, indicating potential improvements of a factor ~2.5 over current published constraints. We perform resolution studies on scales larger than ~$16h^{-1}$ Mpc showing the promise of significant constraints with next-generation surveys. Furthermore, the results demonstrate that our method can consistently marginalize all nuisance parameters of the data model. The method further provides an inference of the three-dimensional primordial density field, providing opportunities to explore additional signatures of primordial physics. This first demonstration of a field-level inference pipeline demonstrates a promising complementary path forward for analysing next-generation surveys.

**Key words:** galaxies: statistics – cosmological parameters – inflation – large-scale structure of Universe.

## 1 INTRODUCTION

Deviations of primordial curvature fluctuations from a Gaussian random field, known as primordial non-Gaussianity (PNG), can constrain primordial universe physics and shed light on the number of fields driving inflation (Maldacena 2003; Chen 2010; Chen & Wang 2010; Komatsu 2010; Senatore, Smith & Zaldarriaga 2010; Barnaby & Shandera 2012; Alvarez et al. 2014; Celoria & Matarrese 2018). The deviation from the initial Gaussian conditions is commonly parametrized by the non-linearity parameter, $f_{nl}^{\text{local}}$, here denoted $f_{nl}$. To date, the strongest constraint on $f_{nl}$ is based on observations of the cosmic microwave background (CMB). CMB observations obtained by the Planck satellite yielded $f_{nl} = -0.9 \pm 5.1$ (Planck Collaboration IX 2019a). While the CMB is a powerful cosmological probe, its information content at the largest scales has reached its cosmic variance limit, and constraints of order unity, $|f_{nl}| < 1$, are expected to be beyond its reach (Camera et al.

2013; Ballardini, Matthewson & Maartens 2019; Meerburg et al. 2019; Moradinezhad Dizgah & Keating 2019; Karagiannis, Slosar & Liguori 2020).

For current state-of-the-art techniques which constrain PNG with the Large-Scale Structure (LSS), the most informative probe measures the two-point correlation function of galaxy populations. The reason for this is the scale-dependent bias, which gives the most information from the largest scales (Dalal et al. 2008; Slosar et al. 2008). Previous work has been successful in constraining PNG in available galaxy redshift surveys (Slosar et al. 2008; Ross et al. 2012; Karagiannis, Shanks & Ross 2014; Leistedt, Peiris & Roth 2014; Uhlemann et al. 2018; Castorina et al. 2019; Mueller, Percival & Ruggeri 2019; Mueller et al. 2021; Cabass et al. 2022; D'Amico et al. 2022); the tightest large-scale structure measurement comes from the SDSS-IV (DR 16, quasar sample) data, reaching a constraining power of $|f_{nl}| < 21$ (Mueller et al. 2021).

Next-generation galaxy redshift surveys [Vera C. Rubin Observatory (LSST Science Collaboration 2009), Euclid (Amendola et al. 2018), and SPHEREx (Doré et al. 2014)] will study the Universe at new levels of precision. These missions will be dominated by their

<footnote n="1"><sup>1</sup>The Oskar Klein Centre, Department of Physics, Stockholm University, AlbaNova University Centre, SE 106 91 Stockholm, Sweden</footnote>

<footnote n="2"><sup>2</sup>CNRS & Sorbonne Université, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany</footnote>

<footnote>*E-mail: adam.andrews@fysik.su.se</footnote>
# Bayesian field-level inference of primordial non-Gaussianity using next-generation galaxy surveys

Adam Andrews<sup>1,★</sup> Jens Jasche,<sup>1,2</sup> Guilhem Lavaux<sup>1,2</sup> and Fabian Schmidt<sup>3</sup>

Accepted 2023 January 19. Received 2022 December 16; in original form 2022 March 16

## ABSTRACT

Detecting and measuring a non-Gaussian signature of primordial origin in the density field is a major science goal of next-generation galaxy surveys. The signal will permit us to determine primordial-physics processes and constrain models of cosmic inflation. While traditional approaches use a limited set of statistical summaries of the galaxy distribution to constrain primordial non-Gaussianity, we present a field-level approach by Bayesian forward modelling the entire three-dimensional galaxy survey. Since our method includes the entire cosmic field in the analysis, it can naturally and fully self-consistently exploit all available information in the large-scale structure, to extract information on the local non-Gaussianity parameter, $f_{nl}$. Examples include higher order statistics through correlation functions, peculiar velocity fields through redshift-space distortions, and scale-dependent galaxy bias. To illustrate the feasibility of field-level primordial non-Gaussianity inference, we present our approach using a first-order Lagrangian perturbation theory model, approximating structure growth at sufficiently large scales. We demonstrate the performance of our approach through various tests with self-consistent mock galaxy data emulating relevant features of the SDSS-III/BOSS-like survey, and additional tests with a *Stage IV* mock data set. These tests reveal that the method infers unbiased values of $f_{nl}$ by accurately handling survey geometries, noise, and unknown galaxy biases. We demonstrate that our method can achieve constraints of $\sigma_{f_{nl}} \approx 8.78$ for SDSS-III/BOSS-like data, indicating potential improvements of a factor ~2.5 over current published constraints. We perform resolution studies on scales larger than ~$16h^{-1}$ Mpc showing the promise of significant constraints with next-generation surveys. Furthermore, the results demonstrate that our method can consistently marginalize all nuisance parameters of the data model. The method further provides an inference of the three-dimensional primordial density field, providing opportunities to explore additional signatures of primordial physics. This first demonstration of a field-level inference pipeline demonstrates a promising complementary path forward for analysing next-generation surveys.

**Key words:** galaxies: statistics – cosmological parameters – inflation – large-scale structure of Universe.

## 1 INTRODUCTION

Deviations of primordial curvature fluctuations from a Gaussian random field, known as primordial non-Gaussianity (PNG), can constrain primordial universe physics and shed light on the number of fields driving inflation (Maldacena 2003; Chen 2010; Chen & Wang 2010; Komatsu 2010; Senatore, Smith & Zaldarriaga 2010; Barnaby & Shandera 2012; Alvarez et al. 2014; Celoria & Matarrese 2018). The deviation from the initial Gaussian conditions is commonly parametrized by the non-linearity parameter, $f_{nl}^{\text{local}}$, here denoted $f_{nl}$. To date, the strongest constraint on $f_{nl}$ is based on observations of the cosmic microwave background (CMB). CMB observations obtained by the Planck satellite yielded $f_{nl} = -0.9 \pm 5.1$ (Planck Collaboration IX 2019a). While the CMB is a powerful cosmological probe, its information content at the largest scales has reached its cosmic variance limit, and constraints of order unity, $|f_{nl}| < 1$, are expected to be beyond its reach (Camera et al.

2013; Ballardini, Matthewson & Maartens 2019; Meerburg et al. 2019; Moradinezhad Dizgah & Keating 2019; Karagiannis, Slosar & Liguori 2020).

For current state-of-the-art techniques which constrain PNG with the Large-Scale Structure (LSS), the most informative probe measures the two-point correlation function of galaxy populations. The reason for this is the scale-dependent bias, which gives the most information from the largest scales (Dalal et al. 2008; Slosar et al. 2008). Previous work has been successful in constraining PNG in available galaxy redshift surveys (Slosar et al. 2008; Ross et al. 2012; Karagiannis, Shanks & Ross 2014; Leistedt, Peiris & Roth 2014; Uhlemann et al. 2018; Castorina et al. 2019; Mueller, Percival & Ruggeri 2019; Mueller et al. 2021; Cabass et al. 2022; D'Amico et al. 2022); the tightest large-scale structure measurement comes from the SDSS-IV (DR 16, quasar sample) data, reaching a constraining power of $|f_{nl}| < 21$ (Mueller et al. 2021).

Next-generation galaxy redshift surveys [Vera C. Rubin Observatory (LSST Science Collaboration 2009), Euclid (Amendola et al. 2018), and SPHEREX (Doré et al. 2014)] will study the Universe at new levels of precision. These missions will be dominated by their

<footnote n="1"><sup>1</sup>The Oskar Klein Centre, Department of Physics, Stockholm University, AlbaNova University Centre, SE 106 91 Stockholm, Sweden</footnote>

<footnote n="2"><sup>2</sup>CNRS & Sorbonne Université, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany</footnote>

<footnote>*E-mail: adam.andrews@fysik.su.se</footnote>
systematic effects, rather than by their statistical power (Graham et al. 2018). Therefore, a high degree of mitigation and modelling of the survey effects, instrumentation noise, and astrophysical contamination is required. If not taken into account, these effects can bias the results on the largest scales (Huterer, Cunha & Fang 2013; Leistedt et al. 2014; Ho et al. 2015; Jasche & Lavaux 2017, 2019).

A promising approach to analyse large-scale structure data is the field-level inference method (Jasche et al. 2010; Wang et al. 2014; Seljak et al. 2017; Lavaux, Jasche & Leclercq 2019; Schmittfull et al. 2019; Nguyen et al. 2020; Porqueres et al. 2020; Schmidt et al. 2020). This method has the goal of going beyond individual correlation functions and extract the maximum amount of information to solve cosmological problems, e.g. constraining cosmological parameters (Ramanah et al. 2019; Leclercq & Heavens 2021). Since the full cosmological map is included into the analysis, a field-level inference can incorporate other informative probes into the analysis, e.g. peculiar velocity fields and higher order statistics in the cosmic density field (Biagetti 2019; Baumann & Green 2021), which have been shown to be independently sensitive to PNG (Catelan & Scherrer 1995; Schmidt 2010; Baldauf, Seljak & Senatore 2011; Lam, Nishimichi & Yoshida 2011; Ma, Taylor & Scott 2013; Tasinato et al. 2014).

In this paper, we present for the first time a proof-of-concept to infer PNG in galaxy redshift surveys, with the aim of addressing current open issues. Our method is based on the physical forward modelling of the three-dimensional galaxy distribution in a Bayesian hierarchical framework. This implies that the algorithm performs a field-level inference of the full cosmic density field (Jasche & Wandelt 2013). In this way, our method self-consistently and naturally accounts for the current phenomenology of PNG when analysing galaxy data, and uses all of the information available in the cosmic density field when exploring the solution space (Jasche & Wandelt 2013). The physics-informed inference algorithm provides a path to account for various uncertainties and observational and systematic effects which otherwise can bias the cosmological conclusions drawn from data. These effects include survey geometry, selection effects, instrumentation noise, galaxy biases, and foreground contamination of data (Jasche & Lavaux 2017; Porqueres et al. 2019). Furthermore, the algorithm naturally accounts for higher order statistics and filamentary cosmic structure associated with non-linear structure formation, and can incorporate any probe associated with PNG. As a result, our method constrains $f_{nl}$ beyond the current state-of-the-art methods and can place constraints as low as $|f_{nl}| < 5.70$, for mock data emulating features of next-generation galaxy redshift surveys (see Fig. 1), while marginalizing out systematic effects.

A major advantage of the forward modelling approach is that the full 3D density field is part of the analysis (Jasche & Wandelt 2013; Baumann & Green 2021). Therefore, BORG naturally and fully self-consistently accounts for a variety of probes when constraining $f_{nl}$:

(i) higher order statistics of primordial origin in the dark matter density field (Baldauf et al. 2011; Tasinato et al. 2014);

(ii) mass distributions of the large-scale structure (for example, statistical moments such as skewness and kurtosis of the galaxies) (Chodorowski & Bouchet 1996; Durrer et al. 2000; Dalal et al. 2008; Yokoyama et al. 2011; Friedrich et al. 2019);

(iii) peculiar velocity fields (Catelan & Scherrer 1995; Schmidt 2010; Lam et al. 2011; Ma et al. 2013);

(iv) scale-dependent galaxy bias (de Putter, Roland & Doré 2017; Castorina et al. 2019).

For more details on the listed probes, the interested reader is referred to reviews found in the literature (Biagetti 2019). Another

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Conditional posterior distribution for the local type primordial non-Gaussianity parameter $f_{nl}$ obtained from a BORG application to a *Stage IV* mock galaxy survey as demonstrated in this paper. This posterior incorporates all signatures of primordial non-Gaussianity in the cosmic large-scale structure, for the scales included in this project. Furthermore, our method correctly handles systematic errors and survey effects, as well as marginalizing out bias parameters. The yellow band indicates the region $|f_{nl}| \le 1$, the target region of next-generation surveys. To compare with the current best constraints, we plot a Gaussian centred on the fiducial value of $f_{nl}$, with a width corresponding to the Planck18 constraints (Planck Collaboration X 2019b).</figcaption></figure>

advantage is the fact that our method allows us to simultaneously analyse effects that are degenerate in other analysis methods. In our physics model (see Section 2), the input field is iteratively built upon to create a prediction of the model. Thus, different effects enter the field at different locations in the pipeline, meaning that these effects can be disentangled. For example, while both structure formation and the perturbation of the primordial gravitational potential affect the bispectrum (or three-point correlation function) of the density field, BORG is able to jointly explore these two effects in its framework, because BORG explicitly accounts for the non-linear gravitational evolution.

The paper is structured as follows. We describe the Bayesian Origin Reconstruction from Galaxies (BORG) algorithm together with the developments made to the physical model to constrain $f_{nl}$ in Section 2. In Section 3, we describe the artificial mock data set generation, detailing both the SDSS-III/BOSS-like mock data set and the *Stage IV* mock data set. The results are described in Section 4, where we evaluate the performance of the PNG inference framework within BORG. We provide concluding thoughts and a summary in Section 5.

## 2 PHYSICAL MODEL AND STATISTICAL METHOD

In this work, we rely on the BORG algorithm to run a Markov chain which produces Monte Carlo samples of the joint posterior distribution of the 3D initial conditions, the galaxy bias parameters, and constraints on PNG (as parametrized by the local form). This section describes the implementation of the Bayesian forward modelling approach to infer the PNG signal from galaxy data, starting with a condensed overview of the BORG algorithm. Then we describe the necessary modifications of its data model to account for the phenomenology of PNG. Specifically, we will include the generation of primordial matter fluctuation fields (equation 4) and a
systematic effects, rather than by their statistical power (Graham et al. 2018). Therefore, a high degree of mitigation and modelling of the survey effects, instrumentation noise, and astrophysical contamination is required. If not taken into account, these effects can bias the results on the largest scales (Huterer, Cunha & Fang 2013; Leistedt et al. 2014; Ho et al. 2015; Jasche & Lavaux 2017, 2019).

A promising approach to analyse large-scale structure data is the field-level inference method (Jasche et al. 2010; Wang et al. 2014; Seljak et al. 2017; Lavaux, Jasche & Leclercq 2019; Schmittfull et al. 2019; Nguyen et al. 2020; Porqueres et al. 2020; Schmidt et al. 2020). This method has the goal of going beyond individual correlation functions and extract the maximum amount of information to solve cosmological problems, e.g. constraining cosmological parameters (Ramanah et al. 2019; Leclercq & Heavens 2021). Since the full cosmological map is included into the analysis, a field-level inference can incorporate other informative probes into the analysis, e.g. peculiar velocity fields and higher order statistics in the cosmic density field (Biagetti 2019; Baumann & Green 2021), which have been shown to be independently sensitive to PNG (Catelan & Scherrer 1995; Schmidt 2010; Baldauf, Seljak & Senatore 2011; Lam, Nishimichi & Yoshida 2011; Ma, Taylor & Scott 2013; Tasinato et al. 2014).

In this paper, we present for the first time a proof-of-concept to infer PNG in galaxy redshift surveys, with the aim of addressing current open issues. Our method is based on the physical forward modelling of the three-dimensional galaxy distribution in a Bayesian hierarchical framework. This implies that the algorithm performs a field-level inference of the full cosmic density field (Jasche & Wandelt 2013). In this way, our method self-consistently and naturally accounts for the current phenomenology of PNG when analysing galaxy data, and uses all of the information available in the cosmic density field when exploring the solution space (Jasche & Wandelt 2013). The physics-informed inference algorithm provides a path to account for various uncertainties and observational and systematic effects which otherwise can bias the cosmological conclusions drawn from data. These effects include survey geometry, selection effects, instrumentation noise, galaxy biases, and foreground contamination of data (Jasche & Lavaux 2017; Porqueres et al. 2019). Furthermore, the algorithm naturally accounts for higher order statistics and filamentary cosmic structure associated with non-linear structure formation, and can incorporate any probe associated with PNG. As a result, our method constrains $f_{nl}$ beyond the current state-of-the-art methods and can place constraints as low as $|f_{nl}| < 5.70$, for mock data emulating features of next-generation galaxy redshift surveys (see Fig. 1), while marginalizing out systematic effects.

A major advantage of the forward modelling approach is that the full 3D density field is part of the analysis (Jasche & Wandelt 2013; Baumann & Green 2021). Therefore, BORG naturally and fully self-consistently accounts for a variety of probes when constraining $f_{nl}$:

(i) higher order statistics of primordial origin in the dark matter density field (Baldauf et al. 2011; Tasinato et al. 2014);

(ii) mass distributions of the large-scale structure (for example, statistical moments such as skewness and kurtosis of the galaxies) (Chodorowski & Bouchet 1996; Durrer et al. 2000; Dalal et al. 2008; Yokoyama et al. 2011; Friedrich et al. 2019);

(iii) peculiar velocity fields (Catelan & Scherrer 1995; Schmidt 2010; Lam et al. 2011; Ma et al. 2013);

(iv) scale-dependent galaxy bias (de Putter, Roland & Doré 2017; Castorina et al. 2019).

For more details on the listed probes, the interested reader is referred to reviews found in the literature (Biagetti 2019). Another

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Conditional posterior distribution for the local type primordial non-Gaussianity parameter $f_{nl}$ obtained from a BORG application to a *Stage IV* mock galaxy survey as demonstrated in this paper. This posterior incorporates all signatures of primordial non-Gaussianity in the cosmic large-scale structure, for the scales included in this project. Furthermore, our method correctly handles systematic errors and survey effects, as well as marginalizing out bias parameters. The yellow band indicates the region $|f_{nl}| \le 1$, the target region of next-generation surveys. To compare with the current best constraints, we plot a Gaussian centred on the fiducial value of $f_{nl}$, with a width corresponding to the Planck18 constraints (Planck Collaboration X 2019b).</figcaption></figure>

advantage is the fact that our method allows us to simultaneously analyse effects that are degenerate in other analysis methods. In our physics model (see Section 2), the input field is iteratively built upon to create a prediction of the model. Thus, different effects enter the field at different locations in the pipeline, meaning that these effects can be disentangled. For example, while both structure formation and the perturbation of the primordial gravitational potential affect the bispectrum (or three-point correlation function) of the density field, BORG is able to jointly explore these two effects in its framework, because BORG explicitly accounts for the non-linear gravitational evolution.

The paper is structured as follows. We describe the Bayesian Origin Reconstruction from Galaxies (BORG) algorithm together with the developments made to the physical model to constrain $f_{nl}$ in Section 2. In Section 3, we describe the artificial mock data set generation, detailing both the SDSS-III/BOSS-like mock data set and the *Stage IV* mock data set. The results are described in Section 4, where we evaluate the performance of the PNG inference framework within BORG. We provide concluding thoughts and a summary in Section 5.

## 2 PHYSICAL MODEL AND STATISTICAL METHOD

In this work, we rely on the BORG algorithm to run a Markov chain which produces Monte Carlo samples of the joint posterior distribution of the 3D initial conditions, the galaxy bias parameters, and constraints on PNG (as parametrized by the local form). This section describes the implementation of the Bayesian forward modelling approach to infer the PNG signal from galaxy data, starting with a condensed overview of the BORG algorithm. Then we describe the necessary modifications of its data model to account for the phenomenology of PNG. Specifically, we will include the generation of primordial matter fluctuation fields (equation 4) and a
<figure><img src="image_2.png" /><figcaption>**Figure 2.** These two flow charts describe the overall process of the BORG algorithm. On the left side, we have the Bayesian hierarchical framework, which illustrates the parameter sampling steps. For each sample generated, one full cycle is performed. For more information on the individual sampling schemes, see Section 2.4. For each parameter evaluation, the algorithm runs one full cycle of the forward model, which is depicted on the right side of the plot. This flow chart illustrates the process of bringing a set of initial conditions (the white noise field $\epsilon$) to a prediction for the expected galaxy field. The likelihood is then computed (as in equation 12). In this way, the maximum amount of information is conserved in the data, since the algorithm can iteratively fit the field to the data. The yellow dashed regions are the stochastic phases of the algorithm, which are probabilistically explored in the sampling steps, while the blue lined regions are the deterministic regions; here, the algorithm transforms the input field accordingly. The letter (A)->(H) are references to more detailed descriptions in the text.¹</figcaption></figure>

scale-dependent bias (equation 10) to relate observed galaxies to the underlying matter distribution.

## 2.1 The BORG algorithm

The BORG algorithm is a Bayesian hierarchical inference framework aiming to analyse the cosmic structure by physical forward modelling three-dimensional galaxy fields in cosmological surveys (Jasche & Wandelt 2013; Jasche, Leclercq & Wandelt 2015; Lavaux & Jasche 2016; Jasche & Lavaux 2019; Ramanah et al. 2019). It uses models of gravitational structure formation to link the observed distribution of galaxies to the initial three-dimensional density field. In this way, the algorithm re-formulates the inverse problem of inferring the large-scale structures into an initial conditions problem. The naive procedure would be the following: start by proposing a set of initial conditions which is physically forward modelled to an observable galaxy field. First, the gravitational evolution of the linear matter field is simulated by a structure formation process. This yields physical realizations of the underlying dark matter density field in the present Universe. Second, the evolved density field is populated through a galaxy data model. Through this process, one would generate the model predictions which are compared to observed galaxy redshift data. A simplified overview of the algorithm for this paper is given in Fig. 2. BORG’s ability to infer the primordial matter fluctuation field from galaxy surveys is key to this work.

The BORG algorithm uses a combination of Hamiltonian Markov chain Monte Carlo and slice sampling techniques to explore the joint posterior of the data (Jasche & Wandelt 2012; Jasche & Lavaux 2019). This includes the full three-dimensional initial density fields, cosmological parameters, and nuisance parameters of the data model, such as galaxy bias parameters. The algorithm performs a statistically rigorous analysis by executing an iterative Markov chain. The primordial initial conditions, including PNG, are combined with the gravitational structure formation model to generate model predictions. In this case, model predictions are the galaxy number counts in individual volume elements. These galaxy number counts are then compared to the observed galaxy distribution using a likelihood estimation, which accounts for the observational noise and the galaxy bias. For every model algorithm iteration, the $f_{nl}$

parameter, the bias parameters, and initial conditions are updated. The iteration of these processes yields a valid Markov chain that allows us to quantify the significance of inferred quantities (Jasche & Wandelt 2013).

As outlined, the BORG algorithm permits us to infer the three-dimensional initial conditions out of which the present structures formed. To achieve optimal extraction, we will use a new forward model to jointly infer the three-dimensional initial white noise field $\epsilon$ and the non-Gaussianity parameter $f_{nl}$, conditioned on the galaxy number counts data $N_g^o$, with $g$ indicating the galaxy catalogue index and $o$ that it is an observed quantity. The vector $\epsilon$ describes the phases of the cosmic structure, while $f_{nl}$ is introduced in equation (4). Formally, the posterior $\pi(\epsilon, f_{nl}, \{b_i\}|N_g^o)$ can be written as

$$ \pi(\epsilon, f_{nl}, \{b_i\}|N_g^o) = \pi(f_{nl}) \pi(\epsilon) \pi(\{b_i^g\}) \frac{\pi(N_g^o|\epsilon, f_{nl}, \{b_i\})}{\pi(N_g^o)}, \quad (1) $$

where $\pi(\epsilon)$ is the Gaussian white-noise prior with zero mean and unit standard variance, $\pi(N_g^o|\epsilon, f_{nl})$ is a likelihood distribution, and $\{b_i\}$ are the bias parameters. It is important to remark that in this formulation, the entire data modelling is expressed within the likelihood distribution, which includes the generation of initial conditions and the subsequent forward modelling steps.

## 2.2 Forward modelling large-scale structures with primordial gravitational potential

In this section, we will discuss the modifications required to the physical model used in BORG to account for PNG.

### 2.2.1 Generating non-Gaussian initial conditions

In this study, we focus on inferring the local PNG at the leading order, expressed through the $f_{nl}$ parameter, using a forward modelling approach.

To describe the fluctuations of $\phi_G$, we convolve the Gaussian white noise field $\epsilon$ with a first transfer function $T_G(k)$. That transfer function is designed such that $\phi_G$ has a correlation structure given by a primordial power spectrum. In this work, we use the classical
<figure><img src="image_4.png" /><figcaption>**Figure 2.** These two flow charts describe the overall process of the BORG algorithm. On the left side, we have the Bayesian hierarchical framework, which illustrates the parameter sampling steps. For each sample generated, one full cycle is performed. For more information on the individual sampling schemes, see Section 2.4. For each parameter evaluation, the algorithm runs one full cycle of the forward model, which is depicted on the right side of the plot. This flow chart illustrates the process of bringing a set of initial conditions (the white noise field $\epsilon$) to a prediction for the expected galaxy field. The likelihood is then computed (as in equation 12). In this way, the maximum amount of information is conserved in the data, since the algorithm can iteratively fit the field to the data. The yellow dashed regions are the stochastic phases of the algorithm, which are probabilistically explored in the sampling steps, while the blue lined regions are the deterministic regions; here, the algorithm transforms the input field accordingly. The letter (A)->(H) are references to more detailed descriptions in the text.¹</figcaption></figure>

scale-dependent bias (equation 10) to relate observed galaxies to the underlying matter distribution.

## 2.1 The BORG algorithm

The BORG algorithm is a Bayesian hierarchical inference framework aiming to analyse the cosmic structure by physical forward modelling three-dimensional galaxy fields in cosmological surveys (Jasche & Wandelt 2013; Jasche, Leclercq & Wandelt 2015; Lavaux & Jasche 2016; Jasche & Lavaux 2019; Ramanah et al. 2019). It uses models of gravitational structure formation to link the observed distribution of galaxies to the initial three-dimensional density field. In this way, the algorithm re-formulates the inverse problem of inferring the large-scale structures into an initial conditions problem. The naive procedure would be the following: start by proposing a set of initial conditions which is physically forward modelled to an observable galaxy field. First, the gravitational evolution of the linear matter field is simulated by a structure formation process. This yields physical realizations of the underlying dark matter density field in the present Universe. Second, the evolved density field is populated through a galaxy data model. Through this process, one would generate the model predictions which are compared to observed galaxy redshift data. A simplified overview of the algorithm for this paper is given in Fig. 2. BORG’s ability to infer the primordial matter fluctuation field from galaxy surveys is key to this work.

The BORG algorithm uses a combination of Hamiltonian Markov chain Monte Carlo and slice sampling techniques to explore the joint posterior of the data (Jasche & Wandelt 2012; Jasche & Lavaux 2019). This includes the full three-dimensional initial density fields, cosmological parameters, and nuisance parameters of the data model, such as galaxy bias parameters. The algorithm performs a statistically rigorous analysis by executing an iterative Markov chain. The primordial initial conditions, including PNG, are combined with the gravitational structure formation model to generate model predictions. In this case, model predictions are the galaxy number counts in individual volume elements. These galaxy number counts are then compared to the observed galaxy distribution using a likelihood estimation, which accounts for the observational noise and the galaxy bias. For every model algorithm iteration, the $f_{nl}$

parameter, the bias parameters, and initial conditions are updated. The iteration of these processes yields a valid Markov chain that allows us to quantify the significance of inferred quantities (Jasche & Wandelt 2013).

As outlined, the BORG algorithm permits us to infer the three-dimensional initial conditions out of which the present structures formed. To achieve optimal extraction, we will use a new forward model to jointly infer the three-dimensional initial white noise field $\epsilon$ and the non-Gaussianity parameter $f_{nl}$, conditioned on the galaxy number counts data $N_g^o$, with $g$ indicating the galaxy catalogue index and $o$ that it is an observed quantity. The vector $\epsilon$ describes the phases of the cosmic structure, while $f_{nl}$ is introduced in equation (4). Formally, the posterior $\pi(\epsilon, f_{nl}, \{b_i\}|N_g^o)$ can be written as

$$ \pi(\epsilon, f_{nl}, \{b_i\}|N_g^o) = \pi(f_{nl}) \pi(\epsilon) \pi(\{b_i^g\}) \frac{\pi(N_g^o|\epsilon, f_{nl}, \{b_i\})}{\pi(N_g^o)}, \quad (1) $$

where $\pi(\epsilon)$ is the Gaussian white-noise prior with zero mean and unit standard variance, $\pi(N_g^o|\epsilon, f_{nl})$ is a likelihood distribution, and $\{b_i\}$ are the bias parameters. It is important to remark that in this formulation, the entire data modelling is expressed within the likelihood distribution, which includes the generation of initial conditions and the subsequent forward modelling steps.

## 2.2 Forward modelling large-scale structures with primordial gravitational potential

In this section, we will discuss the modifications required to the physical model used in BORG to account for PNG.

### 2.2.1 Generating non-Gaussian initial conditions

In this study, we focus on inferring the local PNG at the leading order, expressed through the $f_{nl}$ parameter, using a forward modelling approach.

To describe the fluctuations of $\phi_G$, we convolve the Gaussian white noise field $\epsilon$ with a first transfer function $T_G(k)$. That transfer function is designed such that $\phi_G$ has a correlation structure given by a primordial power spectrum. In this work, we use the classical
parametrization of the primordial power spectrum in terms of $A_s$ and $n_s$. In the discrete case, the covariance of the field $\phi_G$ becomes

$$ \langle \hat{\phi}_{G,a} \hat{\phi}_{G,b}^* \rangle = V \delta_{a,-b}^K A_s \frac{2\pi^2}{k_a^2} \left( \frac{k_a}{k_{\text{pivot}}} \right)^{n_s-1} \equiv \delta_{a,-b}^K T_G(k_a)^2, \quad (2) $$

with V the volume of the box holding the forward-modelled universe, **a** and **b** mesh indices, and $\delta_{a,b}^K$ the Kronecker delta. $A_s$ provides the amplitude of adiabatic, scalar, and potential fluctuations at the limit $f_{nl} = 0$. Thereafter, we generate the real space primordial gravitational potential as

$$ \phi_G(\mathbf{x}) = \frac{1}{V} \sum_{\mathbf{a}} \exp(i\mathbf{k}_{\mathbf{a}} \cdot \mathbf{x}) \hat{\phi}_{G,\mathbf{a}}, \quad (3) $$

with **a** a mesh index that runs through the 3D modes of the box, $\hat{\phi}_{G,a}$ the Fourier amplitude of the modes of the Gaussian random field $\phi_G(x)$.

In this model, a non-Gaussian primordial Bardeen potential $\Phi_{\text{NG}}(\mathbf{x})$ is generated from a given Gaussian field $\phi_G(x)$ through the following Taylor expansion, as highlighted in the box (C) in Fig. 2:

$$ \Phi_{\text{NG}}(\mathbf{x}) = \phi_G(\mathbf{x}) + f_{\text{nl}} (\phi_G(\mathbf{x})^2 - \langle \phi_G(\mathbf{x}) \rangle^2) \quad (4) $$

(Komatsu & Spergel 2001). The parameters $A_s$, $n_s$, and $f_{nl}$ can be related to the inflation potential (Salopek & Bond 1990; Verde et al. 2000).

Once the Bardeen potential is calculated, one can follow the usual modelling techniques to set up the initial conditions to start evolving large-scale structures (Peebles 1980; Bond & Efstathiou 1984; Crocce, Pueblas & Scoccimarro 2006). As evolution at high redshift is linear, the relation between matter density $\delta_{\text{L}}(\mathbf{k})$ and primordial potential is as highlighted by the box (D) in Fig. 2:

$$ \hat{\delta}_{\text{L}}(\mathbf{k}, z) = T(k, z) \Phi_{\text{NG}}(\mathbf{k}), \quad (5) $$

where $z$ is taken at an arbitrary reference epoch. The transfer function $T(k, z)$, which contains the linear growth factor, is provided by CLASS (Lesgourgues & Tram 2014) or by analytic computation (Eisenstein & Hu 1998, 1999; Takada, Komatsu & Futamase 2006).

### 2.2.2 Simulating structure growth

For this proof-of-concept, we construct the matter density field at present time by relying on the first-order Lagrangian perturbation theory (LPT), commonly known as Zel'dovich approximation. Although LPT is sufficient to describe large-scale features, small-scale features suffer from approximation issues. In this work, we focus on scales of $\ge 16 h^{-1}$ Mpc, where LPT is a good approximation. It should be noted that the implementation described here allows us to further improve and modify the forward model. As such, a more accurate structure formation model would allow the method to go beyond the current resolution limit. In particular, other structure growth models exist which can provide us with results including more physical effects, e.g. the particle mesh model (Jasche & Lavaux 2019) or second- (2LPT; Jasche & Lavaux 2015) and higher order LPT (Schmidt 2021). Thus, while in this proof of concept we rely on LPT to describe structure formation, we can improve the structure formation model to provide more realistic non-linear density field realizations, which we will investigate in future projects.

To compute the present-day density field, we populate the Lagrangian coordinates **q** with a set of particles that are displaced according to

$$ \Psi(\mathbf{q}, z) = -\nabla \Delta^{-1} \delta_{\text{L}}(\mathbf{q}, z), \quad (6) $$

as highlighted by the box (E) in Fig. 2. In this work, we evaluate at $z = 0$.

As a final step, we apply redshift-space distortion effects to the particles, by transforming them from their rest frame to the comoving frame, as highlighted by the box (F) in Fig. 2. Handling redshift-space distortion effects allows us to include the peculiar velocities when constraining $f_{nl}$ with our field-based approach (Tellarini et al. 2016; Di Dio et al. 2017; Karagiannis et al. 2018; Bharadwaj, Mazumdar & Sarkar 2020).

The redshift-space distortions are applied to the final matter field $\delta_m(\vec{x})$ with the help of the vector $\vec{s}$:

$$ \vec{s} = \vec{x} + \frac{a}{H(a)} \frac{\vec{v} \cdot \vec{x}}{|\mathbf{x}|^2} \vec{x}. \quad (7) $$

The observed final matter field $\delta_m(\vec{s})$ is then evaluated via a cloud-in-cell approach on the particles in $\vec{s}$-space. We note that we account only for linear and mildly non-linear redshift distortions in this fashion, which is adequate for the scales considered in this work. A detailed treatment of non-linear redshift distortions, notably the Finger of God effect, within the hierarchical forward modelling approach will be subject to an upcoming publication.

## 2.3 Scale-dependent galaxy biasing

The next step of the forward modelling is to model the relationship between the predicted matter field and the observation of galaxy number counts. In the field of large-scale structure cosmology, galaxies are considered tracers of the underlying gravitational potential. Specifically, galaxies are biased tracers, displaying clustering properties approximately aligned with the underlying dark matter (Kaiser 1984). The uncertainty of this model systematic effect is currently one of the most important unresolved issues in LSS cosmology, hindering the non-linear analysis of galaxy surveys (Neyrinck et al. 2014; Desjacques, Jeong & Schmidt 2018; Wechsler & Tinker 2018). For a review of the galaxy bias problem, the interested reader is referred to the literature (Desjacques et al. 2018; Schmidt et al. 2019).

The validity of the galaxy bias formalism at non-linear scales is currently being investigated in current research. Linear bias relationships have currently been employed in PNG studies to constrain $f_{nl}$ in galaxy redshift surveys, where the analyses use scale cuts, e.g. $k_{\text{max}} = 0.12 h \text{ Mpc}^{-1}$ (Mueller et al. 2021) and $k_{\text{max}} = 0.3 h \text{ Mpc}^{-1}$ (Castorina et al. 2019). To go beyond these scales into the non-linear regime, more advanced galaxy bias relationships are required. There are different approaches to do this; examples include applying Effective Field Theory of Large-Scale Structures [EFToLSS (Desjacques et al. 2018; Schmidt et al. 2019)], or applying machine learning techniques (Charnock et al. 2020). None the less, at the scales covered in this project, a linear bias relationship should be sufficient to model the formation of galaxies in the large-scale structure.

At linear order, the relationship is the linear bias function: $\rho_g = \langle N_g^o \rangle (1 + b_1^g \delta_m)$, where $\rho_g$ is the resulting galaxy field, $\langle N_g^o \rangle$ is the mean number of galaxies, $b_1^g$ is the linear bias, and $\delta_m$ is the final density field, in the present-time universe. $g$ is the index of the galaxy catalogue. In addition to a linear bias, non-Gaussian perturbations imprint a scale-dependent bias effect in dark matter halo formation. This effective parameter describes the effect of the PNG-induced coupling between the long-wavelength potential perturbations and the short-wavelength modes of the density field on galaxy formation (Dalal et al. 2008; Slosar et al. 2008; Desjacques, Seljak & Iliev 2009). Although the scale dependence of this bias is known, its amplitude $b_\phi$ depends on the details of the given galaxy sample
parametrization of the primordial power spectrum in terms of $A_s$ and $n_s$. In the discrete case, the covariance of the field $\phi_G$ becomes

$$ \langle \hat{\phi}_{G,a} \hat{\phi}_{G,b}^* \rangle = V \delta_{a,-b}^K A_s \frac{2\pi^2}{k_a^2} \left( \frac{k_a}{k_{\text{pivot}}} \right)^{n_s-1} = \delta_{a,-b}^K T_G(k_a)^2, \quad (2) $$

with V the volume of the box holding the forward-modelled universe, a and b mesh indices, and $\delta_{a,b}^K$ the Kronecker delta. $A_s$ provides the amplitude of adiabatic, scalar, and potential fluctuations at the limit $f_{nl} = 0$. Thereafter, we generate the real space primordial gravitational potential as

$$ \phi_G(\mathbf{x}) = \frac{1}{V} \sum_{\mathbf{a}} \exp(i\mathbf{k}_{\mathbf{a}} \cdot \mathbf{x}) \hat{\phi}_{G,a}, \quad (3) $$

with a a mesh index that runs through the 3D modes of the box, $\hat{\phi}_{G,a}$ the Fourier amplitude of the modes of the Gaussian random field $\phi_G(x)$.

In this model, a non-Gaussian primordial Bardeen potential $\Phi_{\text{NG}}(\mathbf{x})$ is generated from a given Gaussian field $\phi_G(x)$ through the following Taylor expansion, as highlighted in the box (C) in Fig. 2:

$$ \Phi_{\text{NG}}(\mathbf{x}) = \phi_G(\mathbf{x}) + f_{\text{nl}} (\phi_G(\mathbf{x})^2 - \langle \phi_G(\mathbf{x})^2 \rangle) \quad (4) $$

(Komatsu & Spergel 2001). The parameters $A_s$, $n_s$, and $f_{nl}$ can be related to the inflation potential (Salopek & Bond 1990; Verde et al. 2000).

Once the Bardeen potential is calculated, one can follow the usual modelling techniques to set up the initial conditions to start evolving large-scale structures (Peebles 1980; Bond & Efstathiou 1984; Crocce, Pueblas & Scoccimarro 2006). As evolution at high redshift is linear, the relation between matter density $\delta_{\text{L}}(\mathbf{k})$ and primordial potential is as highlighted by the box (D) in Fig. 2:

$$ \hat{\delta}_{\text{L}}(\mathbf{k}, z) = T(k, z) \Phi_{\text{NG}}(\mathbf{k}), \quad (5) $$

where $z$ is taken at an arbitrary reference epoch. The transfer function $T(k, z)$, which contains the linear growth factor, is provided by CLASS (Lesgourgues & Tram 2014) or by analytic computation (Eisenstein & Hu 1998, 1999; Takada, Komatsu & Futamase 2006).

### 2.2.2 Simulating structure growth

For this proof-of-concept, we construct the matter density field at present time by relying on the first-order Lagrangian perturbation theory (LPT), commonly known as Zel'dovich approximation. Although LPT is sufficient to describe large-scale features, small-scale features suffer from approximation issues. In this work, we focus on scales of $\ge 16 h^{-1}$ Mpc, where LPT is a good approximation. It should be noted that the implementation described here allows us to further improve and modify the forward model. As such, a more accurate structure formation model would allow the method to go beyond the current resolution limit. In particular, other structure growth models exist which can provide us with results including more physical effects, e.g. the particle mesh model (Jasche & Lavaux 2019) or second- (2LPT; Jasche & Lavaux 2015) and higher order LPT (Schmidt 2021). Thus, while in this proof of concept we rely on LPT to describe structure formation, we can improve the structure formation model to provide more realistic non-linear density field realizations, which we will investigate in future projects.

To compute the present-day density field, we populate the Lagrangian coordinates $\mathbf{q}$ with a set of particles that are displaced according to

$$ \Psi(\mathbf{q}, z) = -\nabla \Delta^{-1} \delta_{\text{L}}(\mathbf{q}, z), \quad (6) $$

as highlighted by the box (E) in Fig. 2. In this work, we evaluate at $z = 0$.

As a final step, we apply redshift-space distortion effects to the particles, by transforming them from their rest frame to the comoving frame, as highlighted by the box (F) in Fig. 2. Handling redshift-space distortion effects allows us to include the peculiar velocities when constraining $f_{nl}$ with our field-based approach (Tellarini et al. 2016; Di Dio et al. 2017; Karagiannis et al. 2018; Bharadwaj, Mazumdar & Sarkar 2020).

The redshift-space distortions are applied to the final matter field $\delta_m(\vec{x})$ with the help of the vector $\vec{s}$:

$$ \vec{s} = \vec{x} + \frac{a}{H(a)} \frac{\vec{v} \cdot \vec{x}}{|\mathbf{x}|^2} \vec{x}. \quad (7) $$

The observed final matter field $\delta_m(\vec{s})$ is then evaluated via a cloud-in-cell approach on the particles in $\vec{s}$-space. We note that we account only for linear and mildly non-linear redshift distortions in this fashion, which is adequate for the scales considered in this work. A detailed treatment of non-linear redshift distortions, notably the Finger of God effect, within the hierarchical forward modelling approach will be subject to an upcoming publication.

## 2.3 Scale-dependent galaxy biasing

The next step of the forward modelling is to model the relationship between the predicted matter field and the observation of galaxy number counts. In the field of large-scale structure cosmology, galaxies are considered tracers of the underlying gravitational potential. Specifically, galaxies are biased tracers, displaying clustering properties approximately aligned with the underlying dark matter (Kaiser 1984). The uncertainty of this model systematic effect is currently one of the most important unresolved issues in LSS cosmology, hindering the non-linear analysis of galaxy surveys (Neyrinck et al. 2014; Desjacques, Jeong & Schmidt 2018; Wechsler & Tinker 2018). For a review of the galaxy bias problem, the interested reader is referred to the literature (Desjacques et al. 2018; Schmidt et al. 2019).

The validity of the galaxy bias formalism at non-linear scales is currently being investigated in current research. Linear bias relationships have currently been employed in PNG studies to constrain $f_{nl}$ in galaxy redshift surveys, where the analyses use scale cuts, e.g. $k_{\text{max}} = 0.12 h \text{ Mpc}^{-1}$ (Mueller et al. 2021) and $k_{\text{max}} = 0.3 h \text{ Mpc}^{-1}$ (Castorina et al. 2019). To go beyond these scales into the non-linear regime, more advanced galaxy bias relationships are required. There are different approaches to do this; examples include applying Effective Field Theory of Large-Scale Structures [EFToLSS (Desjacques et al. 2018; Schmidt et al. 2019)], or applying machine learning techniques (Charnock et al. 2020). None the less, at the scales covered in this project, a linear bias relationship should be sufficient to model the formation of galaxies in the large-scale structure.

At linear order, the relationship is the linear bias function: $\rho_g = \langle N_g^o \rangle (1 + b_1^g \delta_m)$, where $\rho_g$ is the resulting galaxy field, $\langle N_g^o \rangle$ is the mean number of galaxies, $b_1^g$ is the linear bias, and $\delta_m$ is the final density field, in the present-time universe. $g$ is the index of the galaxy catalogue. In addition to a linear bias, non-Gaussian perturbations imprint a scale-dependent bias effect in dark matter halo formation. This effective parameter describes the effect of the PNG-induced coupling between the long-wavelength potential perturbations and the short-wavelength modes of the density field on galaxy formation (Dalal et al. 2008; Slosar et al. 2008; Desjacques, Seljak & Iliev 2009). Although the scale dependence of this bias is known, its amplitude $b_\phi$ depends on the details of the given galaxy sample
and cannot be predicted in general; see Pillepich, Porciani & Hahn (2009) and Biagetti et al. (2017) for measurements of $b_\phi$ on simulated dark matter haloes and Barreira et al. (2020) for a measurement on simulated galaxies.

Here, we will assume the universality relation for $b_\phi$ (Dalal et al. 2008; Slosar et al. 2008), which relates the amplitude of the scale-dependent bias to the linear bias $b_1^g$.

In terms of the bias contribution:

$$ \rho_g = \langle N_g^o \rangle \{1 + [b_1^g + \Delta b(k, f_{nl})] \delta_m \}, $$

the scale-dependent bias $\Delta b(k, f_{nl})$ then takes the following shape:

$$ \Delta b(k, f_{nl}) = b_{\phi} \alpha(k, z) = b_{\phi} f_{nl} \frac{3\Omega_m H_0^2}{2c^2 k^2 T(k, z)}. \quad (8) $$

Under the universal mass approximation (Press & Schechter 1974; Sheth & Tormen 1999; Pillepich et al. 2009), $b_\phi$ takes the form of (Matarrese & Verde 2008; Slosar et al. 2008; Schmidt, Jeong & Desjacques 2013):

$$ b_{\phi} = 2\delta_c(b_1^g - 1), \quad (9) $$

where $\delta_c$ is the spherical collapse threshold, $\delta_c \approx 1.686$ (Gunn & Gott 1972).

Thus, the galaxy bias model is formulated with linear bias $b_1^g$ and with the scale-dependent bias effect $b_\phi(k, f_{nl})$, which translates the dark matter density contrast field $\delta_m$ into the galaxy field $\rho_g$ as highlighted by the box (G) in Fig. 2:

$$ \rho_g(\delta_m, \langle N_g \rangle, b_1^g, f_{nl}) = \langle N_g \rangle \{1 + [b_1^g + \Delta b(k, f_{nl})] \delta_m \}. \quad (10) $$

This is the expression for the galaxy field in a universe with local PNG, to the leading order. We leave the incorporation of higher order bias terms of local PNG to future work (see e.g. Assassi, Baumann & Schmidt 2015; Barreira 2020, 2021a, b; Barreira et al. 2020; Moradinezhad Dizgah et al. 2021, for a description of the said higher order terms). Notice also that the scale-dependent bias should strictly involve the primordial potential at the Lagrangian (initial) position, rather than the late-time Eulerian position (Assassi et al. 2015). The difference is, however, only a second-order effect.

Finally, when constraining $f_{nl}$ from galaxy clustering, it is also in general important to include the so-called relativistic effects (see Jeong & Schmidt 2015, for a review). We leave the incorporation of these effects into the forward model for future work.

### 2.3.1 A likelihood model of galaxy data

With a prediction of $\rho_g$, the next step is to relate this field to $N_{g,p}$, the expected number of observed galaxies at a given voxel $p$. Since the galaxies are split into several catalogues, each with different survey geometries and nuisance parameters, we denote the catalogue for each quantity with $g$. In this way, by applying $R_{g,p}$, the linear response operator of the survey, we obtain as highlighted by the box (H) in Fig. 2:

$$ N_{g,p} = R_{g,p}[\rho_{g,p}(\delta_m, \langle N_g \rangle, b_1^g, f_{nl})], \quad (11) $$

where $\langle N_g \rangle$ is the mean number of galaxies in the catalogue $g$. The purpose of $R_{g,p}$ is to account for survey geometry and selection effects. We emphasize that the galaxy data model accounts for the fact that galaxies are biased tracers of the underlying density field, as well as for the discrete nature of galaxy formation (Layzer 1956; Peebles 1980).

As for the noise of the galaxy counts, in this work, we are at sufficiently large scales to describe it as a Gaussian likelihood

(Kitaura & Enßlin 2008; Jasche et al. 2010; Jasche & Lavaux 2017). For a galaxy catalogue $g$, the log-likelihood is given as

$$ \ln[\pi(N_{g,p}^o | N_{g,p}, \sigma_g)] = -P \ln[\sqrt{2\pi} \sigma_g] \\ - \frac{1}{2} \sum_{p=0}^{P-1} \left( \frac{N_{g,p}^o - N_{g,p}}{\sigma_g} \right)^2, \quad (12) $$

where $N_{g,p}^o$ is the observed galaxy number count in voxel $p$, and $N_{g,p}$ is the model prediction of the physical forward model, given $\epsilon$ and $f_{nl}$. $p$ is the voxel index, and $P$ is the total number of voxels. $\sigma_g$ is the noise of the galaxy catalogue; in this work, this is fixed to the mean number of galaxies in catalogue $g$:

$$ \sigma_g = \sqrt{\langle N_g^o \rangle} = \sqrt{\frac{\sum_p N_{g,p}}{\sum_p R_{g,p}}}, \quad (13) $$

where $R_{g,p}$ is the selection function for each voxel; summing over $R_{g,p}$ yields the total number of observed voxels. This is chosen to emulate Poisson noise. A comment on the likelihood model: While the uncertainty on galaxy counts is usually assumed to be Poisson-distributed (Layzer 1956; Peebles 1980; Ramanah et al. 2019), in this study we assume that it is Gaussian-distributed. The justification for this is that we have a low resolution in the voxels, resulting in a relatively high number of galaxies per voxel. As the Poisson distribution approaches a Gaussian distribution in the limit of large values, we argue that this choice has little effect on the outcome of the analysis (Nguyen et al. 2021).

In this way, we obtain a likelihood distribution that contains the statistical process of generating galaxy observations given a sampled set of initial conditions. Hence, due to the deterministic nature of structure formation, the expected number of galaxies per voxel $N_{g,p}$ is related to the white noise field $\epsilon$ through the physical forward model (described in Section 2.2) and the galaxy bias model (discussed in Section 2.3).

## 2.4 Sampling of density field, bias parameters, and $f_{nl}$

Within the framework of BORG, the main objective is to explore the full joint posterior distribution of the cosmic density field, bias parameters and $f_{nl}$. However, sampling directly from the joint posterior is difficult, possibly numerically infeasible. Instead, BORG evaluates the conditional posterior distributions separately, to sample new proposals for the white-noise field $\epsilon$, galaxy bias parameter $\{b_1^g\}$ (where $g$ is the galaxy catalogue number, and $i$ is the index for the galaxy bias parameter), and $f_{nl}$. In this way, the conditional posteriors can be sampled sequentially. This block sampling strategy is outlined in the following way:

$$ \begin{align} (1) \quad \epsilon^{n+1} &\cap \pi(\epsilon | \{b_1^g\}^n, f_{nl}^n) && (\text{Hamiltonian Monte Carlo}), \nonumber \\ (2) \quad \{b_1^g\}^{n+1} &\cap \pi(\{b_1^g\} | f_{nl}^n, \epsilon^{n+1}) && (\text{Slice sampler}), \nonumber \\ (3) \quad f_{nl}^{n+1} &\cap \pi(f_{nl}^n | \{b_1^g\}^{n+1}, \epsilon^{n+1}) && (\text{Slice sampler}), \tag{14} \end{align} $$

where $\cap$ represents the process of drawing a value from the probability distribution, $n$ is the index of the sample, and the sampling algorithm is described in parentheses. By consecutive sampling one parameter of the physical model at a time, the conditional distributions form a sample of the full joint posterior distribution (Hastings 1970). See the left-hand side of Fig. 2 for a graphical representation of the sampling schemes.

For the field $\epsilon$, the main problem consists of sampling a high-dimensional parameter space, with each dimension consisting of the amplitude in each voxel. Utilizing standard Markov chain
and cannot be predicted in general; see Pillepich, Porciani & Hahn (2009) and Biagetti et al. (2017) for measurements of $b_\phi$ on simulated dark matter haloes and Barreira et al. (2020) for a measurement on simulated galaxies.

Here, we will assume the universality relation for $b_\phi$ (Dalal et al. 2008; Slosar et al. 2008), which relates the amplitude of the scale-dependent bias to the linear bias $b_1^g$.

In terms of the bias contribution:

$$ \rho_g = \langle N_g^o \rangle \{1 + [b_1^g + \Delta b(k, f_{nl})] \delta_m \}, $$

the scale-dependent bias $\Delta b(k, f_{nl})$ then takes the following shape:

$$ \Delta b(k, f_{nl}) = b_{\phi} \alpha(k, z) = b_{\phi} f_{nl} \frac{3\Omega_m H_0^2}{2c^2 k^2 T(k, z)}. \quad (8) $$

Under the universal mass approximation (Press & Schechter 1974; Sheth & Tormen 1999; Pillepich et al. 2009), $b_\phi$ takes the form of (Matarrese & Verde 2008; Slosar et al. 2008; Schmidt, Jeong & Desjacques 2013):

$$ b_{\phi} = 2\delta_c(b_1^g - 1), \quad (9) $$

where $\delta_c$ is the spherical collapse threshold, $\delta_c \approx 1.686$ (Gunn & Gott 1972).

Thus, the galaxy bias model is formulated with linear bias $b_1^g$ and with the scale-dependent bias effect $b_\phi(k, f_{nl})$, which translates the dark matter density contrast field $\delta_m$ into the galaxy field $\rho_g$ as highlighted by the box (G) in Fig. 2:

$$ \rho_g(\delta_m, \langle N_g \rangle, b_1^g, f_{nl}) = \langle N_g \rangle \{1 + [b_1^g + \Delta b(k, f_{nl})] \delta_m \}. \quad (10) $$

This is the expression for the galaxy field in a universe with local PNG, to the leading order. We leave the incorporation of higher order bias terms of local PNG to future work (see e.g. Assassi, Baumann & Schmidt 2015; Barreira 2020, 2021a, b; Barreira et al. 2020; Moradinezhad Dizgah et al. 2021, for a description of the said higher order terms). Notice also that the scale-dependent bias should strictly involve the primordial potential at the Lagrangian (initial) position, rather than the late-time Eulerian position (Assassi et al. 2015). The difference is, however, only a second-order effect.

Finally, when constraining $f_{nl}$ from galaxy clustering, it is also in general important to include the so-called relativistic effects (see Jeong & Schmidt 2015, for a review). We leave the incorporation of these effects into the forward model for future work.

### 2.3.1 A likelihood model of galaxy data

With a prediction of $\rho_g$, the next step is to relate this field to $N_{g,p}$, the expected number of observed galaxies at a given voxel p. Since the galaxies are split into several catalogues, each with different survey geometries and nuisance parameters, we denote the catalogue for each quantity with g. In this way, by applying $R_{g,p}$, the linear response operator of the survey, we obtain as highlighted by the box (H) in Fig. 2:

$$ N_{g,p} = R_{g,p}[\rho_{g,p}(\delta_m, \langle N_g \rangle, b_1^g, f_{nl})], \quad (11) $$

where $\langle N_g \rangle$ is the mean number of galaxies in the catalogue g. The purpose of $R_{g,p}$ is to account for survey geometry and selection effects. We emphasize that the galaxy data model accounts for the fact that galaxies are biased tracers of the underlying density field, as well as for the discrete nature of galaxy formation (Layzer 1956; Peebles 1980).

As for the noise of the galaxy counts, in this work, we are at sufficiently large scales to describe it as a Gaussian likelihood

(Kitaura & Enßlin 2008; Jasche et al. 2010; Jasche & Lavaux 2017). For a galaxy catalogue g, the log-likelihood is given as

$$ \ln[\pi(N_{g,p}^o | N_{g,p}, \sigma_g)] = -P \ln[\sqrt{2\pi} \sigma_g] - \frac{1}{2} \sum_{p=0}^{P-1} \left( \frac{N_{g,p}^o - N_{g,p}}{\sigma_g} \right)^2, \quad (12) $$

where $N_{g,p}^o$ is the observed galaxy number count in voxel p, and $N_{g,p}$ is the model prediction of the physical forward model, given $\epsilon$ and $f_{nl}$. p is the voxel index, and P is the total number of voxels. $\sigma_g^o$ is the noise of the galaxy catalogue; in this work, this is fixed to the mean number of galaxies in catalogue g:

$$ \sigma_g = \sqrt{\langle N_g^o \rangle} = \sqrt{\frac{\sum_p N_{g,p}}{\sum_p R_{g,p}}}, \quad (13) $$

where $R_{g,p}$ is the selection function for each voxel; summing over $R_{g,p}$ yields the total number of observed voxels. This is chosen to emulate Poisson noise. A comment on the likelihood model: While the uncertainty on galaxy counts is usually assumed to be Poisson-distributed (Layzer 1956; Peebles 1980; Ramanah et al. 2019), in this study we assume that it is Gaussian-distributed. The justification for this is that we have a low resolution in the voxels, resulting in a relatively high number of galaxies per voxel. As the Poisson distribution approaches a Gaussian distribution in the limit of large values, we argue that this choice has little effect on the outcome of the analysis (Nguyen et al. 2021).

In this way, we obtain a likelihood distribution that contains the statistical process of generating galaxy observations given a sampled set of initial conditions. Hence, due to the deterministic nature of structure formation, the expected number of galaxies per voxel $N_{g,p}$ is related to the white noise field $\epsilon$ through the physical forward model (described in Section 2.2) and the galaxy bias model (discussed in Section 2.3).

## 2.4 Sampling of density field, bias parameters, and $f_{nl}$

Within the framework of BORG, the main objective is to explore the full joint posterior distribution of the cosmic density field, bias parameters and $f_{nl}$. However, sampling directly from the joint posterior is difficult, possibly numerically infeasible. Instead, BORG evaluates the conditional posterior distributions separately, to sample new proposals for the white-noise field $\epsilon$, galaxy bias parameter $\{b_1^g\}$ (where g is the galaxy catalogue number, and i is the index for the galaxy bias parameter), and $f_{nl}$. In this way, the conditional posteriors can be sampled sequentially. This block sampling strategy is outlined in the following way:

$$ \begin{align} (1) \quad \epsilon^{n+1} &\cap \pi(\epsilon | \{b_1^g\}^n, f_{nl}^n) && (\text{Hamiltonian Monte Carlo}), \nonumber \\ (2) \quad \{b_1^g\}^{n+1} &\cap \pi(\{b_1^g\} | f_{nl}^n, \epsilon^{n+1}) && (\text{Slice sampler}), \nonumber \\ (3) \quad f_{nl}^{n+1} &\cap \pi(f_{nl}^n | \{b_1^g\}^{n+1}, \epsilon^{n+1}) && (\text{Slice sampler}), \tag{14} \end{align} $$

where $\cap$ represents the process of drawing a value from the probability distribution, n is the index of the sample, and the sampling algorithm is described in parentheses. By consecutive sampling one parameter of the physical model at a time, the conditional distributions form a sample of the full joint posterior distribution (Hastings 1970). See the left-hand side of Fig. 2 for a graphical representation of the sampling schemes.

For the field $\epsilon$, the main problem consists of sampling a high-dimensional parameter space, with each dimension consisting of the amplitude in each voxel. Utilizing standard Markov chain
<table><caption>**Table 1.** An overview of the runs performed in this work.</caption><thead><tr><th>#</th><th>Survey mask</th><th>Box size h<sup>-1</sup> Mpc</th><th>Resolution h<sup>-1</sup> Mpc</th><th>k<sub>max</sub> h<sup>-1</sup> Mpc<sup>-1</sup></th><th>f<sub>nl</sub><sup>true</sup></th></tr></thead><tbody><tr><td>1</td><td>SDSS-III/BOSS</td><td>4000</td><td>62.5</td><td>0.1</td><td>5</td></tr><tr><td>2</td><td>SDSS-III/BOSS</td><td>4000</td><td>31.25</td><td>0.2</td><td>5</td></tr><tr><td>3</td><td>SDSS-III/BOSS</td><td>4000</td><td>15.625</td><td>0.4</td><td>5</td></tr><tr><td>4</td><td>Stage IV</td><td>8000</td><td>62.5</td><td>0.1</td><td>5</td></tr><tr><td>5</td><td>Stage IV</td><td>8000</td><td>31.25</td><td>0.2</td><td>5</td></tr><tr><td>6</td><td>Stage IV (Fixed bias)</td><td>8000</td><td>31.25</td><td>0.2</td><td>5</td></tr></tbody></table>

Note. The main run of this paper is run #5. The other runs serve as complementary runs, which have the purpose of investigating the impact of resolution, performance of the method on SDSS-III/BOSS-like data, to demonstrate the reliability of our method and the effect of fixing bias values.

<table><caption>**Table 2.** Configuration of the catalogues of the tracer and the bias parameters. $b_1$ is the linear bias parameter, $\bar{n}_{gal}$ is the number density, and $z$ is the redshift distance of the tracer bins. The noise $\sigma_g$, introduced in equation (12), is fixed to $\sqrt{\langle N_g^2 \rangle}$, the mean number of galaxies in the catalogue $g$.</caption><thead><tr><th>SDSS-III/BOSS</th><th>b<sub>1</sub><sup>g</sup></th><th>z<sub>min</sub>/z<sub>max</sub></th><th>\bar{n}_{gal} (10<sup>-4</sup> h<sup>3</sup> Mpc<sup>-3</sup>)</th><th>Stage IV</th><th>b<sub>1</sub><sup>g</sup></th><th>z<sub>min</sub>/z<sub>max</sub></th><th>\bar{n}_{gal} (10<sup>-4</sup> h<sup>3</sup> Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td>Cat. 1</td><td>3.22</td><td>0.4/0.75</td><td>3.01</td><td>Cat. 1</td><td>1.46</td><td>0.90/1.10</td><td>6.86</td></tr><tr><td>Cat. 2</td><td>2.48</td><td>0.4/0.75</td><td>2.99</td><td>Cat. 2</td><td>1.61</td><td>1.10/1.30</td><td>5.58</td></tr><tr><td>Cat. 3</td><td>1.78</td><td>0.2/0.5</td><td>6.31</td><td>Cat. 3</td><td>1.75</td><td>1.30/1.50</td><td>4.21</td></tr><tr><td>Cat. 4</td><td>0.87</td><td>0.2/0.5</td><td>10.4</td><td>Cat. 4</td><td>1.90</td><td>1.50/1.80</td><td>2.61</td></tr></tbody></table>

Monte Carlo methods yields insufficient results, mostly due to the numerical intensity of each forward model calculation, coupled with a high rejection rate. Instead, we apply a Hamiltonian Monte Carlo algorithm which guides the proposal procedure of new samples by utilizing dynamical physical symmetries (Duane et al. 1987; Betancourt 2017). This is possible since the physical model that is used to forward model the initial conditions is differentiable, so that the algorithm can utilize its gradient in the Hamiltonian Monte Carlo sampler. For a more detailed description of the sampling of the white-noise field $\epsilon$, see previous works on the subject (Jasche & Wandelt 2013; Lavaux et al. 2019; Ramanah et al. 2019).

In this study, the posterior distributions of $f_{nl}$ and galaxy bias parameters are probed using a slice sampling procedure (Neal 2000, 2003), which has a unit acceptance rate. For each realization of the density field, the algorithm samples the posterior distribution of $f_{nl}$ and accepts a new value for $f_{nl}$ for the next iteration. We use a prior with standard deviation set to 100. In the same fashion as for the $f_{nl}$ parameter, the galaxy bias parameters are sampled and updated by a slice sampling approach.

## 3 GENERATING ARTIFICIAL TEST DATA

In this section, we describe the mock data sets used to test our method. We are particularly interested in testing the method’s ability to constrain $f_{nl}$ when facing survey geometries and selection effects. Generated mock data sets will be designed to emulate the existing SDSS-III/BOSS survey (Alam et al. 2015) and the future Stage IV surveys (LSST Science Collaboration 2009; Doré et al. 2014; Amendola et al. 2018).

### 3.1 Creating mock galaxy surveys

To generate mock data, we will follow similar procedures as described in previous work, using the physical data model as outlined in Section 2 (Jasche & Kitaura 2010; Jasche & Wandelt 2013; Ramanah et al. 2019):

(i) For the sake of demonstration, we evaluate the physical forward model on a cubic Cartesian box of side length $L = 4000 h^{-1}$ Mpc (SDSS-III/BOSS-like) or $L = 8000 h^{-1}$ Mpc (Stage IV) and $N_{grid} = 128$ or 256. This yields grid resolutions in the range of $\Delta L \simeq 62.5 h^{-1}$ Mpc to $\Delta L \simeq 15.6 h^{-1}$ Mpc, well within the range of validity of LPT (Buchert, Melott & Weiss 1994).

(ii) A random three-dimensional field $\epsilon$, with zero mean and unit standard deviation, is generated. Given this white-noise field, we compute a primordial density field by applying the primordial power spectrum, perturbing it with the $f_{nl}$ parameter, and by applying the cosmological transfer function. This yields the linear matter field $\delta_L$, which is the starting point for the gravitational structure formation model. This step follows equations (3) to (5).

(iii) Although the BORG algorithm permits running full gravitational particle mesh simulations, at this coarse resolution LPT provides a viable approximation to the structure formation problem (Moutarde et al. 1991; Buchert et al. 1994; Bouchet et al. 1995; Scoccimarro 2000; Scoccimarro & Sheth 2002). To reduce the noise of the simulated particle distribution, we oversample the initial density by a factor of 2 yielding a total number of $(2N)^3$ simulation particles. Particles are then evolved to the present epoch, using LPT, and assigned to a three-dimensional Cartesian grid via the Cloud-In-Cell kernel to yield the evolved three-dimensional density field. In addition, we add the redshift-space distortions to transform the particles from rest frame to redshift frame.

(iv) To emulate a biased galaxy distribution, we applied the scale-dependent galaxy bias (described in Section 2.3) to the forward-modelled density field, to obtain the galaxy field.

(v) Finally, we apply the radial selection functions and survey geometry to the simulated galaxy field to emulate the observational effects of the respective surveys detailed below.

An overview of the runs and the data sets is organized in Tables 1 and 2, together with detailed parameter choices for the galaxy bias model. For a more detailed description of the SDSS-III/BOSS mock data, see Section 3.2. For a more detailed description of the Stage IV mock data, see Section 3.3.
<table><caption>**Table 1.** An overview of the runs performed in this work.</caption><thead><tr><th>#</th><th>Survey mask</th><th>Box size h<sup>-1</sup> Mpc</th><th>Resolution h<sup>-1</sup> Mpc</th><th>k<sub>max</sub> h<sup>-1</sup> Mpc<sup>-1</sup></th><th>f<sub>nl</sub><sup>true</sup></th></tr></thead><tbody><tr><td>1</td><td>SDSS-III/BOSS</td><td>4000</td><td>62.5</td><td>0.1</td><td>5</td></tr><tr><td>2</td><td>SDSS-III/BOSS</td><td>4000</td><td>31.25</td><td>0.2</td><td>5</td></tr><tr><td>3</td><td>SDSS-III/BOSS</td><td>4000</td><td>15.625</td><td>0.4</td><td>5</td></tr><tr><td>4</td><td>Stage IV</td><td>8000</td><td>62.5</td><td>0.1</td><td>5</td></tr><tr><td>5</td><td>Stage IV</td><td>8000</td><td>31.25</td><td>0.2</td><td>5</td></tr><tr><td>6</td><td>Stage IV (Fixed bias)</td><td>8000</td><td>31.25</td><td>0.2</td><td>5</td></tr></tbody></table>

Note. The main run of this paper is run #5. The other runs serve as complementary runs, which have the purpose of investigating the impact of resolution, performance of the method on SDSS-III/BOSS-like data, to demonstrate the reliability of our method and the effect of fixing bias values.

<table><caption>**Table 2.** Configuration of the catalogues of the tracer and the bias parameters. $b_1$ is the linear bias parameter, $\bar{n}_{gal}$ is the number density, and $z$ is the redshift distance of the tracer bins. The noise $\sigma_g$, introduced in equation (12), is fixed to $\sqrt{\langle N_g^2 \rangle}$, the mean number of galaxies in the catalogue $g$.</caption><thead><tr><th>SDSS-III/BOSS</th><th>b<sub>1</sub><sup>g</sup></th><th>z<sub>min</sub>/z<sub>max</sub></th><th>\bar{n}_{gal} (10<sup>-4</sup> h<sup>3</sup> Mpc<sup>-3</sup>)</th><th>Stage IV</th><th>b<sub>1</sub><sup>g</sup></th><th>z<sub>min</sub>/z<sub>max</sub></th><th>\bar{n}_{gal} (10<sup>-4</sup> h<sup>3</sup> Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td>Cat. 1</td><td>3.22</td><td>0.4/0.75</td><td>3.01</td><td>Cat. 1</td><td>1.46</td><td>0.90/1.10</td><td>6.86</td></tr><tr><td>Cat. 2</td><td>2.48</td><td>0.4/0.75</td><td>2.99</td><td>Cat. 2</td><td>1.61</td><td>1.10/1.30</td><td>5.58</td></tr><tr><td>Cat. 3</td><td>1.78</td><td>0.2/0.5</td><td>6.31</td><td>Cat. 3</td><td>1.75</td><td>1.30/1.50</td><td>4.21</td></tr><tr><td>Cat. 4</td><td>0.87</td><td>0.2/0.5</td><td>10.4</td><td>Cat. 4</td><td>1.90</td><td>1.50/1.80</td><td>2.61</td></tr></tbody></table>

Monte Carlo methods yields insufficient results, mostly due to the numerical intensity of each forward model calculation, coupled with a high rejection rate. Instead, we apply a Hamiltonian Monte Carlo algorithm which guides the proposal procedure of new samples by utilizing dynamical physical symmetries (Duane et al. 1987; Betancourt 2017). This is possible since the physical model that is used to forward model the initial conditions is differentiable, so that the algorithm can utilize its gradient in the Hamiltonian Monte Carlo sampler. For a more detailed description of the sampling of the white-noise field $\epsilon$, see previous works on the subject (Jasche & Wandelt 2013; Lavaux et al. 2019; Ramanah et al. 2019).

In this study, the posterior distributions of $f_{nl}$ and galaxy bias parameters are probed using a slice sampling procedure (Neal 2000, 2003), which has a unit acceptance rate. For each realization of the density field, the algorithm samples the posterior distribution of $f_{nl}$ and accepts a new value for $f_{nl}$ for the next iteration. We use a prior with standard deviation set to 100. In the same fashion as for the $f_{nl}$ parameter, the galaxy bias parameters are sampled and updated by a slice sampling approach.

## 3 GENERATING ARTIFICIAL TEST DATA

In this section, we describe the mock data sets used to test our method. We are particularly interested in testing the method’s ability to constrain $f_{nl}$ when facing survey geometries and selection effects. Generated mock data sets will be designed to emulate the existing SDSS-III/BOSS survey (Alam et al. 2015) and the future Stage IV surveys (LSST Science Collaboration 2009; Doré et al. 2014; Amendola et al. 2018).

### 3.1 Creating mock galaxy surveys

To generate mock data, we will follow similar procedures as described in previous work, using the physical data model as outlined in Section 2 (Jasche & Kitaura 2010; Jasche & Wandelt 2013; Ramanah et al. 2019):

(i) For the sake of demonstration, we evaluate the physical forward model on a cubic Cartesian box of side length $L = 4000 h^{-1}$ Mpc (SDSS-III/BOSS-like) or $L = 8000 h^{-1}$ Mpc (Stage IV) and $N_{grid} = 128$ or 256. This yields grid resolutions in the range of $\Delta L \simeq 62.5 h^{-1}$ Mpc to $\Delta L \simeq 15.6 h^{-1}$ Mpc, well within the range of validity of LPT (Buchert, Melott & Weiss 1994).

(ii) A random three-dimensional field $\epsilon$, with zero mean and unit standard deviation, is generated. Given this white-noise field, we compute a primordial density field by applying the primordial power spectrum, perturbing it with the $f_{nl}$ parameter, and by applying the cosmological transfer function. This yields the linear matter field $\delta_L$, which is the starting point for the gravitational structure formation model. This step follows equations (3) to (5).

(iii) Although the BORG algorithm permits running full gravitational particle mesh simulations, at this coarse resolution LPT provides a viable approximation to the structure formation problem (Moutarde et al. 1991; Buchert et al. 1994; Bouchet et al. 1995; Scoccimarro 2000; Scoccimarro & Sheth 2002). To reduce the noise of the simulated particle distribution, we oversample the initial density by a factor of 2 yielding a total number of $(2N)^3$ simulation particles. Particles are then evolved to the present epoch, using LPT, and assigned to a three-dimensional Cartesian grid via the Cloud-In-Cell kernel to yield the evolved three-dimensional density field. In addition, we add the redshift-space distortions to transform the particles from rest frame to redshift frame.

(iv) To emulate a biased galaxy distribution, we applied the scale-dependent galaxy bias (described in Section 2.3) to the forward-modelled density field, to obtain the galaxy field.

(v) Finally, we apply the radial selection functions and survey geometry to the simulated galaxy field to emulate the observational effects of the respective surveys detailed below.

An overview of the runs and the data sets is organized in Tables 1 and 2, together with detailed parameter choices for the galaxy bias model. For a more detailed description of the SDSS-III/BOSS mock data, see Section 3.2. For a more detailed description of the Stage IV mock data, see Section 3.3.
<figure><img src="image_3.png" /><figcaption>**Figure 3.** The sky maps displaying the observed and masked regions for the SDSS-III/BOSS-like runs of this study. In the left, we have illustrated the CMASS component, while in the right we have illustrated the LOWZ component. (Ross et al. 2017).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 4.** The radial selection functions for the SDSS-III/BOSS-like runs of this study. The solid and dashed lines depict the function for the two galaxy catalogues included in the CMASS sample, the North and the South Galactic Cap (NGC, SGC), while the dash-dotted and dotted lines describe the function for the two galaxy catalogues in the LOWZ component, NGC and SGC. (Ross et al. 2017).</figcaption></figure>

To calculate the cosmological power spectrum and transfer functions, we assume the following set of cosmological parameters: ($\Omega_m = 0.3111$, $\Omega_\Lambda = 0.6889$, $\Omega_b = 0.0490$, $h = 0.6766$, $\sigma_8 = 0.8102$, $n_s = 0.965$), which are taken to be the recent values of Planck + BAO (Planck Collaboration VI 2020). To calculate the required cosmological transfer function we use the prescription provided by Eisenstein & Hu (1998) and Eisenstein & Hu (1999).

The mock data sets are generated with $f_{\text{nl}}^{\text{fiducial}} = f_{\text{nl}}^{\text{true}} = 5$. To reiterate, the main aim of this study is to recover this $f_{\text{nl}}$ value and estimate the accuracy of our method.

## 3.2 Emulating SDSS-III/BOSS data

We seek to generate mock galaxy surveys emulating the major features of the combined CMASS and LOWZ galaxy samples of the SDSS-III/BOSS survey. We will account for the difference in sky completeness and selection effects by dividing the survey into four sub-catalogues: CMASS Northern Galactic Cap (NGC), CMASS Southern Galactic Cap (SGC), LOWZ NGC, and LOWZ SGC. These sky maps are illustrated in Fig. 3. The corresponding radial selection functions have been estimated numerically by binning the equivalent distribution of galaxies $N(d)$ in the actual CMASS and LOWZ surveys (Ross et al. 2017), where $d$ is the comoving distance from the observer. The results are plotted in Fig. 4. We note that by jointly

analysing four galaxy subcatalogues, we will perform a multitracer analysis. As has been previously demonstrated, using multitracer populations reduces the impact of cosmic variance when inferring $f_{\text{nl}}$ (Seljak 2009). However, since the galaxies are placed in separated subcatalogues, we neglect any correlated noise between the galaxy subcatalogues. The specific parameter settings for galaxy biases and the total number of galaxies for the four sub-catalogues are given in Table 2.

The value of expected galaxy number densities $\bar{n}_{\text{gal}}$ for each catalogue have been selected as table 1 in Alam et al. (2017). The expected number density is related to the observed number count as: $\bar{n}_{\text{gal}} = \langle N_g^0 \rangle / V_o$, where $V_o$ is the volume containing the observed galaxies. The linear biases $b_i^g$ for each tracer catalogue have been selected as per fig. 8 in Lavaux et al. (2019), with the light-cone effects adjusted.

As a final note, for the tracer populations, we are only including the galaxies of the SDSS-III/BOSS survey, and the observed quasar populations remain outside of the scope of our data cube. Hence, we are probing smaller volumes compared to other works (Leistedt et al. 2014; Castorina et al. 2019; Mueller et al. 2021).

## 3.3 Emulating Stage IV data

In this study, we apply our method to mock data, which are based on forecasts of next-generation galaxy redshift surveys. In the context of this paper, we dub these mock data sets as *Stage IV*. The purpose of these sets is to emulate typical features of next generation of surveys. In this way, by inferring $f_{\text{nl}}$ in the *Stage IV* mock data sets, we provide an estimate of how well our method could perform on data from the next generation of galaxy redshift surveys.

For the sky area of the *Stage IV* mock data, we design the mock galaxy survey to cover a total area of ~15 000 square degrees. This sky area is spread evenly over the north and south poles of the full sky, resulting in a uniform completeness for $|b| > 39^\circ$. In other words, we exclude observations in regions $\pm 39^\circ$ from the galactic plane. The resulting survey geometry is then represented as a full-sky HEALPIX map with $n_{\text{side}} = 512$ to be used for mock data generation (Gorski et al. 2005).

In the case of the radial selection function, we split the mock data into four galaxy catalogues, each fully covering a redshift bin as in Table 2, with a total redshift range of $0.9 < z < 1.8$. In other words, the radial selection function $S_g(z)$ for the gth galaxy catalogue
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The sky maps displaying the observed and masked regions for the SDSS-III/BOSS-like runs of this study. In the left, we have illustrated the CMASS component, while in the right we have illustrated the LOWZ component. (Ross et al. 2017).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The radial selection functions for the SDSS-III/BOSS-like runs of this study. The solid and dashed lines depict the function for the two galaxy catalogues included in the CMASS sample, the North and the South Galactic Cap (NGC, SGC), while the dash-dotted and dotted lines describe the function for the two galaxy catalogues in the LOWZ component, NGC and SGC. (Ross et al. 2017).</figcaption></figure>

To calculate the cosmological power spectrum and transfer functions, we assume the following set of cosmological parameters: ($\Omega_m = 0.3111$, $\Omega_\Lambda = 0.6889$, $\Omega_b = 0.0490$, $h = 0.6766$, $\sigma_8 = 0.8102$, $n_s = 0.965$), which are taken to be the recent values of Planck + BAO (Planck Collaboration VI 2020). To calculate the required cosmological transfer function we use the prescription provided by Eisenstein & Hu (1998) and Eisenstein & Hu (1999).

The mock data sets are generated with $f_{\text{nl}}^{\text{fiducial}} = f_{\text{nl}}^{\text{true}} = 5$. To reiterate, the main aim of this study is to recover this $f_{\text{nl}}$ value and estimate the accuracy of our method.

## 3.2 Emulating SDSS-III/BOSS data

We seek to generate mock galaxy surveys emulating the major features of the combined CMASS and LOWZ galaxy samples of the SDSS-III/BOSS survey. We will account for the difference in sky completeness and selection effects by dividing the survey into four sub-catalogues: CMASS Northern Galactic Cap (NGC), CMASS Southern Galactic Cap (SGC), LOWZ NGC, and LOWZ SGC. These sky maps are illustrated in Fig. 3. The corresponding radial selection functions have been estimated numerically by binning the equivalent distribution of galaxies $N(d)$ in the actual CMASS and LOWZ surveys (Ross et al. 2017), where $d$ is the comoving distance from the observer. The results are plotted in Fig. 4. We note that by jointly

analysing four galaxy subcatalogues, we will perform a multitracer analysis. As has been previously demonstrated, using multitracer populations reduces the impact of cosmic variance when inferring $f_{\text{nl}}$ (Seljak 2009). However, since the galaxies are placed in separated subcatalogues, we neglect any correlated noise between the galaxy subcatalogues. The specific parameter settings for galaxy biases and the total number of galaxies for the four sub-catalogues are given in Table 2.

The value of expected galaxy number densities $\bar{n}_{\text{gal}}$ for each catalogue have been selected as table 1 in Alam et al. (2017). The expected number density is related to the observed number count as: $\bar{n}_{\text{gal}} = \langle N_g^0 \rangle / V_o$, where $V_o$ is the volume containing the observed galaxies. The linear biases $b_l^g$ for each tracer catalogue have been selected as per fig. 8 in Lavaux et al. (2019), with the light-cone effects adjusted.

As a final note, for the tracer populations, we are only including the galaxies of the SDSS-III/BOSS survey, and the observed quasar populations remain outside of the scope of our data cube. Hence, we are probing smaller volumes compared to other works (Leistedt et al. 2014; Castorina et al. 2019; Mueller et al. 2021).

## 3.3 Emulating Stage IV data

In this study, we apply our method to mock data, which are based on forecasts of next-generation galaxy redshift surveys. In the context of this paper, we dub these mock data sets as *Stage IV*. The purpose of these sets is to emulate typical features of next generation of surveys. In this way, by inferring $f_{\text{nl}}$ in the *Stage IV* mock data sets, we provide an estimate of how well our method could perform on data from the next generation of galaxy redshift surveys.

For the sky area of the *Stage IV* mock data, we design the mock galaxy survey to cover a total area of ~15 000 square degrees. This sky area is spread evenly over the north and south poles of the full sky, resulting in a uniform completeness for $|b| > 39^\circ$. In other words, we exclude observations in regions $\pm 39^\circ$ from the galactic plane. The resulting survey geometry is then represented as a full-sky HEALPIX map with $n_{\text{side}} = 512$ to be used for mock data generation (Gorski et al. 2005).

In the case of the radial selection function, we split the mock data into four galaxy catalogues, each fully covering a redshift bin as in Table 2, with a total redshift range of $0.9 < z < 1.8$. In other words, the radial selection function $S_g(z)$ for the gth galaxy catalogue
<table><caption>**Table 3.** Summary of the measured $f_{nl}$ values for the different configurations. Our method is capable of consistently measuring the fiducial value $f_{nl}$, up to various degrees of uncertainty, for different survey masks, random mock data seeds, bias parameters, and resolutions (reso.). The main result of this paper is run #5.</caption><thead><tr><th>Run #</th><th>Description</th><th>Sampled f<sub>nl</sub></th><th>σ<sub>fnl</sub></th></tr></thead><tbody><tr><td>1</td><td>SDSS-III/BOSS, low reso.</td><td>21.55</td><td>17.5</td></tr><tr><td>2</td><td>SDSS-III/BOSS, mid reso.</td><td>10.11</td><td>11.9</td></tr><tr><td>3</td><td>SDSS-III/BOSS, high reso.</td><td>9.89</td><td>8.78</td></tr><tr><td>4</td><td><i>Stage IV</i>, low reso.</td><td>0.59</td><td>7.09</td></tr><tr><td>5</td><td><i>Stage IV</i>, high reso.</td><td>7.64</td><td>5.70</td></tr><tr><td>6</td><td><i>Stage IV</i>, fixed bias, high reso.</td><td>3.17</td><td>4.43</td></tr></tbody></table>

covering the redshift range $z_{min} \le z \le z_{max}$ can be written as

$$S_g(z) = \begin{cases} 1, & \text{if } z \in [z_{min}, z_{max}] \\ 0, & \text{otherwise.} \end{cases} \qquad (15)$$

For the choices of linear bias values and number densities, see Table 2. The purpose of these sets is to emulate typical features of next-generation surveys, and we take specifications similar to Euclid as a representative case from Euclid Collaboration (2020). For the configuration of resolution and box size, see Table 1. We use these specifications to generate a three-dimensional galaxy distribution that emulates relevant features and systematic effects of the Stage IV surveys. We illustrate the radial selection function in Fig. F4, and the sky maps in Fig. F3; these can be found in the appendix (Section C).

In addition, we also performed a fixed bias test run, using the configurations of the high-resolution Stage IV mock data. For this run, we keep the bias parameters of the MCMC-chain fixed to their fiducial values. This design choice is to test the scenario when we have a complete understanding of the galaxy biasing. In other words, with this simplified test case, we investigate how marginalizing out bias parameters affects the constraining power of $f_{nl}$.

# 4 RESULTS

The previous section outlined the generation of artificial galaxy data used to test our method. This section provides an analysis of the results of our method when applied to the described mock data. Specifically, we are interested in estimating the ability of the method to infer the $f_{nl}$ parameter from the data, which is subject to systematic effects such as survey geometries, selection effects, and noise. A summary of the constrained $f_{nl}$ values can be found in Table 3.

## 4.1 Testing the MCMC warm-up phase

The BORG algorithm executes a large-scale Markov chain Monte Carlo by exploring the joint posterior of initial density fields, the $f_{nl}$, and various nuisance parameters, such as bias parameters or unknown noise levels. To test the validity of the MCMC sampler, we initialize the Markov chain from an overdispersed state; BORG will start in a remote region of the probability space, and then coherently drift towards the region of highest probability and explore these probabilistically correctly. More specifically, we initialize our Markov chain with a randomly overdispersed Gaussian initial density field that is scaled down to have an amplitude the tenth of a realistic cosmological density field. The warm-up phase of the sampler can then be monitored by following the coherent drift of the power spectra estimated from sequential Markov samples. This coherent drift of the power at all Fourier modes is demonstrated in Fig. 5. It can be seen that the se-

<figure><img src="image_5.png" /><figcaption>**Figure 5.** The plot shows the systematic drift of present-time power spectrum amplitudes during burn-in phase of run #5. Different colours of lines denote measurements at different sequential sampling steps in the Markov chain as indicated by the colour bar. It can be seen that the Markov chain initially starts at a remote place in parameter space but approaches the target region in the parameter space after about ~7000 Markov chain transitions. Note that these power spectra have been computed during the post-production process; the forward model is independent of the matter power spectra. Thus, the retrieval of a correct power spectra implies that our method is able to retrieve the fiducial cosmology in the input data.</figcaption></figure>

quentially estimated power spectra approach the correct fiducial cosmological power spectrum after about 7000 Markov transitions.

After passing the warm-up phase, BORG correctly explores the parameter space of plausible large-scale structure realizations. This is reflected by the fact that power spectra correctly fluctuate around the fiducial power spectrum. The corresponding scatter reflects both cosmic variance and observational noise. Note that BORG recovers unbiased estimates of the primordial power spectrum throughout the entire range of Fourier modes considered in this work. This is because our approach correctly handles the systematic effects associated with survey geometries and selection effects. When not taken into account, these systematic effects will produce significant erroneous power at large scales, concealing the physical signal of PNG from $f_{nl}$. This successful test therefore demonstrates that our approach has the potential to study the primordial universe physics on the largest scales in galaxy surveys. We emphasize that the analysis pipeline adopts a pure forward model approach, as illustrated in Fig. 2; thus, by correctly sampling the fiducial power spectra, we
<table><caption>**Table 3.** Summary of the measured $f_{nl}$ values for the different configurations. Our method is capable of consistently measuring the fiducial value $f_{nl}$, up to various degrees of uncertainty, for different survey masks, random mock data seeds, bias parameters, and resolutions (reso.). The main result of this paper is run #5.</caption><thead><tr><th>Run #</th><th>Description</th><th>Sampled f<sub>nl</sub></th><th>σ<sub>fnl</sub></th></tr></thead><tbody><tr><td>1</td><td>SDSS-III/BOSS, low reso.</td><td>21.55</td><td>17.5</td></tr><tr><td>2</td><td>SDSS-III/BOSS, mid reso.</td><td>10.11</td><td>11.9</td></tr><tr><td>3</td><td>SDSS-III/BOSS, high reso.</td><td>9.89</td><td>8.78</td></tr><tr><td>4</td><td>Stage IV, low reso.</td><td>0.59</td><td>7.09</td></tr><tr><td>5</td><td>Stage IV, high reso.</td><td>7.64</td><td>5.70</td></tr><tr><td>6</td><td>Stage IV, fixed bias, high reso.</td><td>3.17</td><td>4.43</td></tr></tbody></table>

covering the redshift range $z_{min} \le z \le z_{max}$ can be written as

$$S_g(z) = \begin{cases} 1, & \text{if } z \in [z_{min}, z_{max}] \\ 0, & \text{otherwise.} \end{cases} \qquad (15)$$

For the choices of linear bias values and number densities, see Table 2. The purpose of these sets is to emulate typical features of next-generation surveys, and we take specifications similar to Euclid as a representative case from Euclid Collaboration (2020). For the configuration of resolution and box size, see Table 1. We use these specifications to generate a three-dimensional galaxy distribution that emulates relevant features and systematic effects of the Stage IV surveys. We illustrate the radial selection function in Fig. F4, and the sky maps in Fig. F3; these can be found in the appendix (Section C).

In addition, we also performed a fixed bias test run, using the configurations of the high-resolution Stage IV mock data. For this run, we keep the bias parameters of the MCMC-chain fixed to their fiducial values. This design choice is to test the scenario when we have a complete understanding of the galaxy biasing. In other words, with this simplified test case, we investigate how marginalizing out bias parameters affects the constraining power of $f_{nl}$.

# 4 RESULTS

The previous section outlined the generation of artificial galaxy data used to test our method. This section provides an analysis of the results of our method when applied to the described mock data. Specifically, we are interested in estimating the ability of the method to infer the $f_{nl}$ parameter from the data, which is subject to systematic effects such as survey geometries, selection effects, and noise. A summary of the constrained $f_{nl}$ values can be found in Table 3.

## 4.1 Testing the MCMC warm-up phase

The BORG algorithm executes a large-scale Markov chain Monte Carlo by exploring the joint posterior of initial density fields, the $f_{nl}$, and various nuisance parameters, such as bias parameters or unknown noise levels. To test the validity of the MCMC sampler, we initialize the Markov chain from an overdispersed state; BORG will start in a remote region of the probability space, and then coherently drift towards the region of highest probability and explore these probabilistically correctly. More specifically, we initialize our Markov chain with a randomly overdispersed Gaussian initial density field that is scaled down to have an amplitude the tenth of a realistic cosmological density field. The warm-up phase of the sampler can then be monitored by following the coherent drift of the power spectra estimated from sequential Markov samples. This coherent drift of the power at all Fourier modes is demonstrated in Fig. 5. It can be seen that the se-

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The plot shows the systematic drift of present-time power spectrum amplitudes during burn-in phase of run #5. Different colours of lines denote measurements at different sequential sampling steps in the Markov chain as indicated by the colour bar. It can be seen that the Markov chain initially starts at a remote place in parameter space but approaches the target region in the parameter space after about ~7000 Markov chain transitions. Note that these power spectra have been computed during the post-production process; the forward model is independent of the matter power spectra. Thus, the retrieval of a correct power spectra implies that our method is able to retrieve the fiducial cosmology in the input data.</figcaption></figure>

quentially estimated power spectra approach the correct fiducial cosmological power spectrum after about 7000 Markov transitions.

After passing the warm-up phase, BORG correctly explores the parameter space of plausible large-scale structure realizations. This is reflected by the fact that power spectra correctly fluctuate around the fiducial power spectrum. The corresponding scatter reflects both cosmic variance and observational noise. Note that BORG recovers unbiased estimates of the primordial power spectrum throughout the entire range of Fourier modes considered in this work. This is because our approach correctly handles the systematic effects associated with survey geometries and selection effects. When not taken into account, these systematic effects will produce significant erroneous power at large scales, concealing the physical signal of PNG from $f_{nl}$. This successful test therefore demonstrates that our approach has the potential to study the primordial universe physics on the largest scales in galaxy surveys. We emphasize that the analysis pipeline adopts a pure forward model approach, as illustrated in Fig. 2; thus, by correctly sampling the fiducial power spectra, we
<figure><img src="image_6.png" /><figcaption>**Figure 6.** The $f_{nl}$ values sampled of run #5: the high-resolution Stage IV data set. The figure displays the conditional posterior distribution of $f_{nl}$, for the given mock data. Notice how the fiducial $f_{nl}$ value is situated within $1\sigma$ of the sampled mean.</figcaption></figure>

demonstrate that our method is able to recover a key observable accurately.

We also perform the Gelman–Rubin convergence test of our chains, which they pass (Gelman & Rubin 1992). The outline of the test and the corresponding results can be found in Appendix D.

## 4.2 Inferring the non-Gaussianity parameter $f_{nl}$

In Fig. 6, the $f_{nl}$-posterior distribution of run #5 is plotted. Noticeably, the distribution displays little skewness, but is a near-Gaussian distribution, with the ground truth value $f_{nl}$ within one $\sigma_{f_{nl}}$.

For each run listed in Table 1, we have computed the uncertainty on the sampled $f_{nl}$ value. These results are summarized in Fig. 7. As can be seen, BORG is able to constrain $f_{nl}$ to $\sigma_{f_{nl}} \approx 9$ for a SDSS-III/BOSS-like data set, with a resolution at $256^3$. For a Stage IV galaxy survey, BORG can bring down the error estimate to $\sigma_{f_{nl}} \approx 6$. In addition, we note the impact of the resolution: by increasing the resolution by a factor of 2, we can bring down the error by a factor of $\sqrt{8}$, when changing the resolution from $128^3$ to $256^3$.

For the interested reader, a full listing of all posterior distributions can be found in the appendix; see Section B1.

## 4.3 Fixed bias test case

In addition to performing resolution studies, we also performed a simplified test case, as outlined in Section 3.3, representing the fixed bias scenario. Using high-resolution Stage IV mock data, we sample the density field and $f_{nl}$-parameter while keeping the bias parameters fixed. This is to test the scenario in which we have complete knowledge of the bias parameters. The improvement on the constrained $f_{nl}$-value can be seen in Table 3. For a graphical visualization, see Fig. 7. The outcome of the test shows that, for our simplified test case, having perfect information of the galaxy bias parameters can yield improvements by a factor of 0.78 on $\sigma_{f_{nl}}$, leading to $\sigma_{f_{nl}} \approx 4.43$.

## 4.4 Accessing the field of primordial matter fluctuations

The Bayesian forward modelling approach goes beyond the estimation of statistical summaries by fully incorporating the three-dimensional initial density field from which the observed structures

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Summary of the results of the runs performed. The vertical axis illustrates the estimated $f_{nl}$ uncertainty for each test run in this paper, as indicated by the horizontal axis. The blue dots represent the SDSS-III/BOSS-like runs, and the red dots represent the Stage IV runs. The corresponding $k_{max}$ for each run can be found in Table 1. Notice how the increase in resolution improves the error on $f_{nl}$, from left to right in each data set. For the last two points, notice the impact of running the analysis with fixed bias parameters compared to marginalizing them out.</figcaption></figure>

formed. In this work, the BORG algorithm infers the primordial gravitation potential from galaxy mock data, which emulates the SDSS-III/BOSS survey or the Stage IV surveys. As a result, the Markov chain provides a numerical representation of the posterior distribution of the initial three-dimensional density fields. Specifically, the Markov chain produces an ensemble of about $\approx 50000$ initial density field realizations that are conditioned on observations made at the present epoch. Given this ensemble, we can estimate any desired statistical summary of the Markov chain. For illustration, we plot the estimated mean and ensemble variance fields of the primordial initial conditions in Fig. 8. Notice how informative regions, i.e. data-filled regions, in the galaxy plot (left) correspond
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The $f_{nl}$ values sampled of run #5: the high-resolution Stage IV data set. The figure displays the conditional posterior distribution of $f_{nl}$, for the given mock data. Notice how the fiducial $f_{nl}$ value is situated within $1\sigma$ of the sampled mean.</figcaption></figure>

demonstrate that our method is able to recover a key observable accurately.

We also perform the Gelman–Rubin convergence test of our chains, which they pass (Gelman & Rubin 1992). The outline of the test and the corresponding results can be found in Appendix D.

## 4.2 Inferring the non-Gaussianity parameter $f_{nl}$

In Fig. 6, the $f_{nl}$-posterior distribution of run #5 is plotted. Noticeably, the distribution displays little skewness, but is a near-Gaussian distribution, with the ground truth value $f_{nl}$ within one $\sigma_{f_{nl}}$.

For each run listed in Table 1, we have computed the uncertainty on the sampled $f_{nl}$ value. These results are summarized in Fig. 7. As can be seen, BORG is able to constrain $f_{nl}$ to $\sigma_{f_{nl}} \approx 9$ for a SDSS-III/BOSS-like data set, with a resolution at $256^3$. For a Stage IV galaxy survey, BORG can bring down the error estimate to $\sigma_{f_{nl}} \approx 6$. In addition, we note the impact of the resolution: by increasing the resolution by a factor of 2, we can bring down the error by a factor of $\sqrt{8}$, when changing the resolution from $128^3$ to $256^3$.

For the interested reader, a full listing of all posterior distributions can be found in the appendix; see Section B1.

## 4.3 Fixed bias test case

In addition to performing resolution studies, we also performed a simplified test case, as outlined in Section 3.3, representing the fixed bias scenario. Using high-resolution Stage IV mock data, we sample the density field and $f_{nl}$-parameter while keeping the bias parameters fixed. This is to test the scenario in which we have complete knowledge of the bias parameters. The improvement on the constrained $f_{nl}$-value can be seen in Table 3. For a graphical visualization, see Fig. 7. The outcome of the test shows that, for our simplified test case, having perfect information of the galaxy bias parameters can yield improvements by a factor of 0.78 on $\sigma_{f_{nl}}$, leading to $\sigma_{f_{nl}} \approx 4.43$.

## 4.4 Accessing the field of primordial matter fluctuations

The Bayesian forward modelling approach goes beyond the estimation of statistical summaries by fully incorporating the three-dimensional initial density field from which the observed structures

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Summary of the results of the runs performed. The vertical axis illustrates the estimated $f_{nl}$ uncertainty for each test run in this paper, as indicated by the horizontal axis. The blue dots represent the SDSS-III/BOSS-like runs, and the red dots represent the Stage IV runs. The corresponding $k_{max}$ for each run can be found in Table 1. Notice how the increase in resolution improves the error on $f_{nl}$, from left to right in each data set. For the last two points, notice the impact of running the analysis with fixed bias parameters compared to marginalizing them out.</figcaption></figure>

formed. In this work, the BORG algorithm infers the primordial gravitation potential from galaxy mock data, which emulates the SDSS-III/BOSS survey or the Stage IV surveys. As a result, the Markov chain provides a numerical representation of the posterior distribution of the initial three-dimensional density fields. Specifically, the Markov chain produces an ensemble of about $\approx 50000$ initial density field realizations that are conditioned on observations made at the present epoch. Given this ensemble, we can estimate any desired statistical summary of the Markov chain. For illustration, we plot the estimated mean and ensemble variance fields of the primordial initial conditions in Fig. 8. Notice how informative regions, i.e. data-filled regions, in the galaxy plot (left) correspond
<figure><img src="image_8.png" /><figcaption>**Figure 8.** The galaxies of the SDSS-III/BOSS-like run ($N = 256$) (left panel), together with the mean (middle panel) and the standard deviation (right panel) over samples of the inferred linear density field, averaged over the z-axis. The galaxies represent the distribution of information in the data. In the mean linear density fields, regions containing galaxies are used by BORG to constrain the linear density field. In voxels absent of information, the density field is only constrained by the physical model and thus converges to the cosmic mean. In the standard deviation plot, the fluctuations are greatly reduced in galaxy-filled regions, as opposed to regions which lack observed data.</figcaption></figure>

to regions where we have a stronger signal in the mean linear density field ($\langle \delta_L \rangle$) (middle), and less variance (right). In other words, the observed regions allow BORG to constrain the underlying $\delta_L$-field, while in masked regions, the $\delta$-field is only constrained by prior coupled to the physical model.

By varying the $f_{nl}$ parameter for a given realization of the present-day density field, one can illustrate the PNG imprints on the cosmic LSS. An example of such a map can be found in Fig. 9. This is generated by subtracting two density fields originating from the same linear density field, which have been forward modelled using two different values of $f_{nl}$ (25, and 0). In this map, hot spots and cold spots represent positions in a universe where PNG has increased or decreased the amount of matter (in relation to a universe without PNG). In principle, these maps could be used by observational missions to further probe the signals of PNG, e.g. searching for small-scale effects of $f_{nl}$ or studying the properties of galaxies in these regions. In other words, our method can highlight regions where PNG is expected to be imprinted given the data (Kostić et al. 2022).

## 4.5 Tests of the scale-dependent bias model

As a concluding discussion, we outline the internal correlations of the model parameters. This is visualized in terms of marginalized posteriors of parameters, and contour plots. For the high-resolution mock data of the SDSS-III/BOSS run, the cross-correlations can be found in Fig. 10. The figure displays the marginalized posteriors of $f_{nl}$ and bias parameters, and the corresponding cross-contour plots. The fiducial values of the mock data are marked with blue lines. Notice how we are able to recover unbiased results of bias parameters (as the fiducial value falls within the $1\sigma$ range of the sampled values). More importantly, we are able to sample $f_{nl}$ while marginalizing the unknown bias parameters without any strong degeneracy within the physical model. This suggests that we can jointly and independently sample $f_{nl}$ together with the bias parameters, without introducing any degeneracies between the model parameters.

We also perform a variety of diagnostic tests; these have the purpose of, among others, investigating the state of the MCMC chains and the robustness of the galaxy bias model. These tests include the correlation matrix and the correlation lengths of the model parameters. The resulting plots, together with a longer discussion, can be found in the Appendix (Section B).

## 5 SUMMARY AND CONCLUSION

This work presents a novel approach to investigating the physics of the primordial universe using Bayesian physical forward modelling of the three-dimensional galaxy distribution in surveys. In contrast to traditional approaches (using only a limited set of statistical summaries), our method provides a full characterization of three-dimensional density and velocity fields including relevant observable signatures of PNG. In addition, our method can measure galaxy bias parameters and provide estimates on their uncertainties. Specifically, our model directly infers information on the cosmic initial conditions and PNG by fitting a physical structure formation model to observations. By modelling data at the field level, our Bayesian forward modelling approach permits us to naturally and fully self-consistently exploit the full phenomenology that PNG imprints on the cosmic large-scale structure, e.g. cluster and void abundances, velocity fields, higher order statistics, and scale-dependent galaxy bias. Moreover, BORG can simultaneously explore all relevant effects (e.g. both the primordial and structure-growth contributions to the bispectrum) and in this way break parameter degeneracies.

To achieve this goal, we have expanded the Bayesian inference algorithm BORG to account for physical models of PNG in the generation of initial conditions. Here, we illustrate our method in the case of PNG of the local type, as parametrized by the cosmological parameter $f_{nl}$. As discussed in Section 2, the inclusion of a physical model of PNG to BORG's forward model enters at two points in the large-scale structure posterior distribution. First, PNG affects the generation of the cosmic initial conditions with the amplitude of non-Gaussian contributions scaled by $f_{nl}$. When evaluated with BORG's
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The galaxies of the SDSS-III/BOSS-like run ($N = 256$) (left panel), together with the mean (middle panel) and the standard deviation (right panel) over samples of the inferred linear density field, averaged over the $z$-axis. The galaxies represent the distribution of information in the data. In the mean linear density fields, regions containing galaxies are used by BORG to constrain the linear density field. In voxels absent of information, the density field is only constrained by the physical model and thus converges to the cosmic mean. In the standard deviation plot, the fluctuations are greatly reduced in galaxy-filled regions, as opposed to regions which lack observed data.</figcaption></figure>

to regions where we have a stronger signal in the mean linear density field ($\langle \delta_L \rangle$) (middle), and less variance (right). In other words, the observed regions allow BORG to constrain the underlying $\delta_L$-field, while in masked regions, the $\delta$-field is only constrained by prior coupled to the physical model.

By varying the $f_{nl}$ parameter for a given realization of the present-day density field, one can illustrate the PNG imprints on the cosmic LSS. An example of such a map can be found in Fig. 9. This is generated by subtracting two density fields originating from the same linear density field, which have been forward modelled using two different values of $f_{nl}$ (25, and 0). In this map, hot spots and cold spots represent positions in a universe where PNG has increased or decreased the amount of matter (in relation to a universe without PNG). In principle, these maps could be used by observational missions to further probe the signals of PNG, e.g. searching for small-scale effects of $f_{nl}$ or studying the properties of galaxies in these regions. In other words, our method can highlight regions where PNG is expected to be imprinted given the data (Kostić et al. 2022).

## 4.5 Tests of the scale-dependent bias model

As a concluding discussion, we outline the internal correlations of the model parameters. This is visualized in terms of marginalized posteriors of parameters, and contour plots. For the high-resolution mock data of the SDSS-III/BOSS run, the cross-correlations can be found in Fig. 10. The figure displays the marginalized posteriors of $f_{nl}$ and bias parameters, and the corresponding cross-contour plots. The fiducial values of the mock data are marked with blue lines. Notice how we are able to recover unbiased results of bias parameters (as the fiducial value falls within the $1\sigma$ range of the sampled values). More importantly, we are able to sample $f_{nl}$ while marginalizing the unknown bias parameters without any strong degeneracy within the physical model. This suggests that we can jointly and independently sample $f_{nl}$ together with the bias parameters, without introducing any degeneracies between the model parameters.

We also perform a variety of diagnostic tests; these have the purpose of, among others, investigating the state of the MCMC chains and the robustness of the galaxy bias model. These tests include the correlation matrix and the correlation lengths of the model parameters. The resulting plots, together with a longer discussion, can be found in the Appendix (Section B).

## 5 SUMMARY AND CONCLUSION

This work presents a novel approach to investigating the physics of the primordial universe using Bayesian physical forward modelling of the three-dimensional galaxy distribution in surveys. In contrast to traditional approaches (using only a limited set of statistical summaries), our method provides a full characterization of three-dimensional density and velocity fields including relevant observable signatures of PNG. In addition, our method can measure galaxy bias parameters and provide estimates on their uncertainties. Specifically, our model directly infers information on the cosmic initial conditions and PNG by fitting a physical structure formation model to observations. By modelling data at the field level, our Bayesian forward modelling approach permits us to naturally and fully self-consistently exploit the full phenomenology that PNG imprints on the cosmic large-scale structure, e.g. cluster and void abundances, velocity fields, higher order statistics, and scale-dependent galaxy bias. Moreover, BORG can simultaneously explore all relevant effects (e.g. both the primordial and structure-growth contributions to the bispectrum) and in this way break parameter degeneracies.

To achieve this goal, we have expanded the Bayesian inference algorithm BORG to account for physical models of PNG in the generation of initial conditions. Here, we illustrate our method in the case of PNG of the local type, as parametrized by the cosmological parameter $f_{nl}$. As discussed in Section 2, the inclusion of a physical model of PNG to BORG's forward model enters at two points in the large-scale structure posterior distribution. First, PNG affects the generation of the cosmic initial conditions with the amplitude of non-Gaussian contributions scaled by $f_{nl}$. When evaluated with BORG's
<figure><img src="image_9.png" /><figcaption>**Figure 9.** The two left panels each display two identical fields generated by the same seed, varying only the $f_{nl}$ value (25 and 0). The difference in the same-seed fields are plotted in the right panel, which depicts the regions of increased and decreased matter as a response to the changing $f_{nl}$ value. These maps illustrate the difference maps as functions of $f_{nl}$, which is a novel data product that our method can infer.</figcaption></figure>

structure growth model, non-Gaussian initial conditions will directly yield imprints on the non-linear matter and velocity distributions.

Secondly, PNG affects the formation of galaxies, which we account for using a model of scale-dependent galaxy bias in our data model. Therefore, the updated physics model naturally and fully self-consistently accounts for a large phenomenology of PNG that is traditionally analysed independently. Besides, our field-level forward modelling approach provides direct solutions to account for selection effects and survey geometry, including measuring galaxy bias parameters. At this point, it should be noted that our physics-informed inference framework can also handle even unknown survey and foreground contamination, a feature that will be particularly relevant for the inference of PNG in the next-generation surveys (Porqueres et al. 2019). A final note: since the perturbation of $f_{nl}$ and gravitational evolution occur in two different steps in the forward model, BORG can differentiate between these two sources of non-linearity in the data analysis procedure.

To test and illustrate the performance of the method, we applied it to two artificial galaxy surveys that emulated the main survey characteristics of the SDSS-III/BOSS and the upcoming Stage IV surveys (see Section 3.1). The tests accounted for the respective survey geometries, selection functions, and noise. We chose to infer the corresponding initial conditions on a Cartesian equidistant grid of side length $L = 4000 \ h^{-1} \ \text{Mpc}$ and $L = 8000 \ h^{-1} \ \text{Mpc}$. We also performed a resolution study by inferring from $N_{\text{grid}} = 64$, $N_{\text{grid}} = 128$, and $N_{\text{grid}} = 256$, resulting in grid resolutions ranging

from $\Delta L = \sim 15.6 \ h^{-1} \ \text{Mpc}$ to $\Delta L = \sim 62.5 \ h^{-1} \ \text{Mpc}$. We also configure a run with fixed bias. This indicates that BORG is capable of handling unknown galaxy biases in a robust framework. Given this coarse resolution, we chose to approximate the formation of the gravitational structure using an LPT model. The investigation of higher resolution with BORG's particle mesh simulation will be the subject of a forthcoming publication.

The results of these tests have been presented and discussed in Section 4. To investigate the Markov chain's initial warm-up phase, we monitored the sequence of inferred posterior matter power spectra. The test revealed that the Markov chain completed the expected coherent drift towards the relevant region of the parameter space that contains most of the probability weight. From these tests, we conclude that the initial warm-up phase takes about 7000 Markov transitions, after which we start recording samples for the analysis. The test further reveals that the Markov chain correctly accounts for the survey geometries and selection effects by inferring the correct shape of the cosmic power spectrum throughout all Fourier modes considered in this analyses. Note that the correct inference of large-scale power is important for the inference of PNG in general and $f_{nl}$ in particular.

One of the major outcomes of this work is the demonstration that the $f_{nl}$ parameter can be constrained with a field-level forward modelling approach. We show that BORG can constrain $f_{nl}$ to $\sigma_{f_{nl}} = 8.78$ for a SDSS-III/BOSS-like survey and $\sigma_{f_{nl}} = 5.7$ for a Stage IV survey. A comparison with other papers is summarized
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The two left panels each display two identical fields generated by the same seed, varying only the $f_{nl}$ value (25 and 0). The difference in the same-seed fields are plotted in the right panel, which depicts the regions of increased and decreased matter as a response to the changing $f_{nl}$ value. These maps illustrate the difference maps as functions of $f_{nl}$, which is a novel data product that our method can infer.</figcaption></figure>

structure growth model, non-Gaussian initial conditions will directly yield imprints on the non-linear matter and velocity distributions.

Secondly, PNG affects the formation of galaxies, which we account for using a model of scale-dependent galaxy bias in our data model. Therefore, the updated physics model naturally and fully self-consistently accounts for a large phenomenology of PNG that is traditionally analysed independently. Besides, our field-level forward modelling approach provides direct solutions to account for selection effects and survey geometry, including measuring galaxy bias parameters. At this point, it should be noted that our physics-informed inference framework can also handle even unknown survey and foreground contamination, a feature that will be particularly relevant for the inference of PNG in the next-generation surveys (Porqueres et al. 2019). A final note: since the perturbation of $f_{nl}$ and gravitational evolution occur in two different steps in the forward model, BORG can differentiate between these two sources of non-linearity in the data analysis procedure.

To test and illustrate the performance of the method, we applied it to two artificial galaxy surveys that emulated the main survey characteristics of the SDSS-III/BOSS and the upcoming Stage IV surveys (see Section 3.1). The tests accounted for the respective survey geometries, selection functions, and noise. We chose to infer the corresponding initial conditions on a Cartesian equidistant grid of side length $L = 4000 \ h^{-1} \ \text{Mpc}$ and $L = 8000 \ h^{-1} \ \text{Mpc}$. We also performed a resolution study by inferring from $N_{\text{grid}} = 64$, $N_{\text{grid}} = 128$, and $N_{\text{grid}} = 256$, resulting in grid resolutions ranging

from $\Delta L = \sim 15.6 \ h^{-1} \ \text{Mpc}$ to $\Delta L = \sim 62.5 \ h^{-1} \ \text{Mpc}$. We also configure a run with fixed bias. This indicates that BORG is capable of handling unknown galaxy biases in a robust framework. Given this coarse resolution, we chose to approximate the formation of the gravitational structure using an LPT model. The investigation of higher resolution with BORG's particle mesh simulation will be the subject of a forthcoming publication.

The results of these tests have been presented and discussed in Section 4. To investigate the Markov chain's initial warm-up phase, we monitored the sequence of inferred posterior matter power spectra. The test revealed that the Markov chain completed the expected coherent drift towards the relevant region of the parameter space that contains most of the probability weight. From these tests, we conclude that the initial warm-up phase takes about 7000 Markov transitions, after which we start recording samples for the analysis. The test further reveals that the Markov chain correctly accounts for the survey geometries and selection effects by inferring the correct shape of the cosmic power spectrum throughout all Fourier modes considered in this analyses. Note that the correct inference of large-scale power is important for the inference of PNG in general and $f_{nl}$ in particular.

One of the major outcomes of this work is the demonstration that the $f_{nl}$ parameter can be constrained with a field-level forward modelling approach. We show that BORG can constrain $f_{nl}$ to $\sigma_{f_{nl}} = 8.78$ for a SDSS-III/BOSS-like survey and $\sigma_{f_{nl}} = 5.7$ for a Stage IV survey. A comparison with other papers is summarized
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Cross-correlations and marginal posteriors for all hyperparameters and $f_{nl}$, for the high-resolution SDSS-III/BOSS-like run. The darker regions in the contour plots represent the $1\sigma$ interval, the lighter regions represent the $2\sigma$ interval. The blue lines mark the fiducial values injected into the mock data, and constitute the true values of the data. The results indicate that our method is able to consistently sample the correct hyperparameters jointly with the true $f_{nl}$ parameter, since the lines overlap well with the contour plots. This demonstrates that our method can jointly infer the bias parameters of the galaxy bias model, and marginalize these out to constrain the $f_{nl}$ parameter.</figcaption></figure>

in Table 4. We studied the inference at two different resolutions, showing that the extracted information on $f_{nl}$ increases when using higher resolution. This result corresponds to our expectations: we expect the number of usable modes for the analysis to scale as $k^3$. The results suggest that more information can be extracted by further increasing resolution, which will be investigated in upcoming studies. Another key point is the fact that our method automatically marginalizes out the nuisance parameters associated with the galaxy bias model. This is an essential feature of the model, since it allows for unbiased estimates of cosmological parameters. As illustrated by Fig. 10, BORG's Markov chain Monte Carlo accurately accounts for

<table><caption>**Table 4.** Comparison of our inferred $f_{nl}$ constraints with other papers measuring $f_{nl}$ in galaxy redshift surveys. Note that while the two last data sets are observed with the same instrument, they cover different regions in the universe and have detected different tracers.</caption><thead><tr><td>Author</td><td>σ<sub>fnl</sub></td><td>Survey</td></tr></thead><tbody><tr><td>This paper</td><td>9</td><td>SDSS-III/BOSS (mock)</td></tr><tr><td>Mueller et al. (2021)</td><td>21</td><td>SDSS-IV/eBOSS</td></tr><tr><td>D'Amico et al. (2022)</td><td>29</td><td>SDSS-III/BOSS</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Cross-correlations and marginal posteriors for all hyperparameters and $f_{nl}$, for the high-resolution SDSS-III/BOSS-like run. The darker regions in the contour plots represent the $1\sigma$ interval, the lighter regions represent the $2\sigma$ interval. The blue lines mark the fiducial values injected into the mock data, and constitute the true values of the data. The results indicate that our method is able to consistently sample the correct hyperparameters jointly with the true $f_{nl}$ parameter, since the lines overlap well with the contour plots. This demonstrates that our method can jointly infer the bias parameters of the galaxy bias model, and marginalize these out to constrain the $f_{nl}$ parameter.</figcaption></figure>

in Table 4. We studied the inference at two different resolutions, showing that the extracted information on $f_{nl}$ increases when using higher resolution. This result corresponds to our expectations: we expect the number of usable modes for the analysis to scale as $k^3$. The results suggest that more information can be extracted by further increasing resolution, which will be investigated in upcoming studies. Another key point is the fact that our method automatically marginalizes out the nuisance parameters associated with the galaxy bias model. This is an essential feature of the model, since it allows for unbiased estimates of cosmological parameters. As illustrated by Fig. 10, BORG's Markov chain Monte Carlo accurately accounts for

<table><caption>**Table 4.** Comparison of our inferred $f_{nl}$ constraints with other papers measuring $f_{nl}$ in galaxy redshift surveys. Note that while the two last data sets are observed with the same instrument, they cover different regions in the universe and have detected different tracers.</caption><thead><tr><td>Author</td><td>σ<sub>f<sub>nl</sub></sub></td><td>Survey</td></tr></thead><tbody><tr><td>This paper</td><td>9</td><td>SDSS-III/BOSS (mock)</td></tr><tr><td>Mueller et al. (2021)</td><td>21</td><td>SDSS-IV/eBOSS</td></tr><tr><td>D'Amico et al. (2022)</td><td>29</td><td>SDSS-III/BOSS</td></tr></tbody></table>
and corrects the correlations between the parameter of interest, $f_{nl}$, and other nuisance parameters. When having complete knowledge of the galaxy bias parameter for a *Stage IV* data analysis, we find that our method can constrain $f_{nl}$ to $\sigma_{f_{nl}} = 4.43$. For more tests and their results, the interested reader is referred to the Appendix (Section B). In other words, using all higher order statistics of galaxy clustering, we can utilize more of the data than traditional data analysis tests (McQuinn 2021).

A major advantage of our approach is its ability to characterize the full three-dimensional field of initial density fluctuations. In Section 4.4, we show inferred three-dimensional initial conditions. In contrast to traditional methods that report only cosmological parameters or compressed statistical summaries, our approach therefore provides the possibility of performing posterior predictive tests on complementary data sets. Besides opening the possibility to search for new physics via cross-correlation with complementary data, the possibility of performing detailed posterior predictive tests adds to the reliability and robustness of our approach. For instance, this opens the immediate opportunity to confirm or detect PNG by cross-correlating our inferences with next-generation CMB data to study the kinetic Synyaev–Zeldovich effect (see e.g. Nguyen et al. 2020). In principle, one can fold the kinetic Synyaev–Zeldovich observable to jointly constrain $f_{nl}$ in the same field-level forward modelling framework, e.g. BORG.

Finally, we touch on some aspects of this project that will be further investigated in future work. First, in this work we utilize LPT to simulate structure growth in the forward model. In order to capture higher order effects in the density field in future projects, we will perform runs at both higher resolution and with a more accurate structure formation model, e.g. 2LPT and particle mesh methods. Second, for the galaxy biasing relationship, we use a linear bias component with a scale-dependent bias addition, which is fixed to the universal mass function prediction. At the moment, it is unclear how well the universal mass function approximation holds for future redshift galaxy surveys, or how tracer dependent the relationship is (Barreira 2020; Barreira et al. 2020; Moradinezhad Dizgah et al. 2021). As the research community makes progress on these topics, and on the galaxy bias relationship overall, this will be further developed in future papers.

In summary, this work presents the first fully Bayesian forward modelling solution to investigate primordial universe physics with galaxy surveys. The results demonstrate the potential of jointly exploiting the entire phenomenology that PNG imprints on the three-dimensional density and velocity fields. Besides providing estimates of the $f_{nl}$ parameter, the method produces maps of three-dimensional density and velocity fields permitting posterior predictive cross-correlation studies with complementary data. By accurately accounting for major survey systematic effects and marginalizing out bias parameters, our method can infer $f_{nl}$ at the level of 5.70 on next-generation galaxy surveys. Future work will focus on increasing the resolution of the reconstructions, thereby resolving more information for constraining PNG, as well as preparing the method for real data applications. Our results demonstrate the promise of Bayesian forward modelling to study the physics of the origin of the Universe with next-generation galaxy surveys.

## ACKNOWLEDGEMENTS

We are grateful to many enlightening discussions with Steffen Hagstotz, Alexandre Barreira, Fabio Finelli, and Yashar Akrami. We also thank Eleni Tsaprazi, Nhat-Minh Nguyen, and Natalia Porqueres for providing feedback on the manuscript. AA acknowl-

edges the travel funding supplied by the University of Oxford, and by the Birger och Gurli Grundströms forskarstipendiefond. JJ acknowledges support by the Swedish Research Council (VR) under the project 2020-05143 – ‘Deciphering the Dynamics of Cosmic Structure’. GL acknowledges support by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. FS acknowledges support from the Starting Grant (ERC-2015-STG 678652) ‘GrInflaGal’ of the European Research Council. This work was supported by the Simons Collaboration on ‘Learning the Universe’.

The computation and data processing in this study were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Tetralith, partially funded by the Swedish Research Council through grant agreement no. 2020-05143. This research utilized the HPC facility supported by the Technical Division of the Department of Physics, Stockholm University. This work is done within the Aquila Consortium.<fnref n="2" /> We acknowledge the use of the following packages: NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), GETDIST (Lewis 2019), and HEALPIX (Gorski et al. 2005).

## DATA AVAILABILITY

The data underlying this article will be shared on the basis of a reasonable request to the corresponding author.

## REFERENCES

* Alam S. et al., 2015, ApJS, 219, 12

* Alam S. et al., 2017, MNRAS, 470, 2617

* Alvarez M. et al., 2014, preprint (arXiv:1412.4671)

* Amendola L. et al., 2018, Living Rev. Relativ., 21, 2

* Assassi V., Baumann D., Schmidt F., 2015, J. Cosmol. Astropart. Phys., 2015, 043

* Baldauf T., Seljak U., Senatore L., 2011, J. Cosmol. Astropart. Phys., 2011, 006

* Ballardini M., Matthewson W. L., Maartens R., 2019, MNRAS, 489, 1950

* Barnaby N., Shandera S., 2012, J. Cosmol. Astropart. Phys., 2012, 034

* Barreira A., 2020, J. Cosmol. Astropart. Phys., 2020, 031

* Barreira A., 2021a, J. Cosmol. Astropart. Phys., 2022, 057

* Barreira A., 2021b, J. Cosmol. Astropart. Phys., 2022, 033

* Barreira A., Cabass G., Schmidt F., Pillepich A., Nelson D., 2020, J. Cosmol. Astropart. Phys., 2020, 013

* Baumann D., Green D., 2021, J. Cosmol. Astropart. Phys., 2022, 061

* Betancourt M., 2017, A Conceptual Introduction to Hamiltonian Monte Carlo

* Bharadwaj S., Mazumdar A., Sarkar D., 2020, MNRAS, 493, 594

* Biagetti M., 2019, Galaxies, 7, 71

* Biagetti M., Lazeyras T., Baldauf T., Desjacques V., Schmidt F., 2017, MNRAS, 468, 3277

* Bond J. R., Efstathiou G., 1984, ApJ, 285, L45

* Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, A&A, 296, 575

* Buchert T., Melott A. L., Weiss A. G., 1994, A&A, 288, 349

* Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022, Phys. Rev. Lett., 129, 21301

* Camera S., Santos M. G., Ferreira P. G., Ferramacho L., 2013, Phys. Rev. Lett., 111, 171302

* Castorina E. et al., 2019, J. Cosmol. Astropart. Phys., 2019, 010

* Catelan P., Scherrer R. J., 1995, ApJ, 445, 1

* Celoria M., Matarrese S., 2018, preprint (arXiv:1812.08197)

* Charnock T., Lavaux G., Wandelt B. D., Boruah S. S., Jasche J., Hudson M. J., 2020, MNRAS, 494, 50

* Chen X., 2010, A&A, 50, 638979

* Chen X., Wang Y., 2010, Phys. Rev. D, 81, 63511

<footnote n="2">²https://www.aquila-consortium.org/</footnote>
and corrects the correlations between the parameter of interest, $f_{nl}$, and other nuisance parameters. When having complete knowledge of the galaxy bias parameter for a Stage IV data analysis, we find that our method can constrain $f_{nl}$ to $\sigma_{f_{nl}} = 4.43$. For more tests and their results, the interested reader is referred to the Appendix (Section B). In other words, using all higher order statistics of galaxy clustering, we can utilize more of the data than traditional data analysis tests (McQuinn 2021).

A major advantage of our approach is its ability to characterize the full three-dimensional field of initial density fluctuations. In Section 4.4, we show inferred three-dimensional initial conditions. In contrast to traditional methods that report only cosmological parameters or compressed statistical summaries, our approach therefore provides the possibility of performing posterior predictive tests on complementary data sets. Besides opening the possibility to search for new physics via cross-correlation with complementary data, the possibility of performing detailed posterior predictive tests adds to the reliability and robustness of our approach. For instance, this opens the immediate opportunity to confirm or detect PNG by cross-correlating our inferences with next-generation CMB data to study the kinetic Synyaev–Zeldovich effect (see e.g. Nguyen et al. 2020). In principle, one can fold the kinetic Synyaev–Zeldovich observable to jointly constrain $f_{nl}$ in the same field-level forward modelling framework, e.g. BORG.

Finally, we touch on some aspects of this project that will be further investigated in future work. First, in this work we utilize LPT to simulate structure growth in the forward model. In order to capture higher order effects in the density field in future projects, we will perform runs at both higher resolution and with a more accurate structure formation model, e.g. 2LPT and particle mesh methods. Second, for the galaxy biasing relationship, we use a linear bias component with a scale-dependent bias addition, which is fixed to the universal mass function prediction. At the moment, it is unclear how well the universal mass function approximation holds for future redshift galaxy surveys, or how tracer dependent the relationship is (Barreira 2020; Barreira et al. 2020; Moradinezhad Dizgah et al. 2021). As the research community makes progress on these topics, and on the galaxy bias relationship overall, this will be further developed in future papers.

In summary, this work presents the first fully Bayesian forward modelling solution to investigate primordial universe physics with galaxy surveys. The results demonstrate the potential of jointly exploiting the entire phenomenology that PNG imprints on the three-dimensional density and velocity fields. Besides providing estimates of the $f_{nl}$ parameter, the method produces maps of three-dimensional density and velocity fields permitting posterior predictive cross-correlation studies with complementary data. By accurately accounting for major survey systematic effects and marginalizing out bias parameters, our method can infer $f_{nl}$ at the level of 5.70 on next-generation galaxy surveys. Future work will focus on increasing the resolution of the reconstructions, thereby resolving more information for constraining PNG, as well as preparing the method for real data applications. Our results demonstrate the promise of Bayesian forward modelling to study the physics of the origin of the Universe with next-generation galaxy surveys.

## ACKNOWLEDGEMENTS

We are grateful to many enlightening discussions with Steffen Hagstotz, Alexandre Barreira, Fabio Finelli, and Yashar Akrami. We also thank Eleni Tsaprazi, Nhat-Minh Nguyen, and Natalia Porqueres for providing feedback on the manuscript. AA acknowl-

edges the travel funding supplied by the University of Oxford, and by the Birger och Gurli Grundströms forskarstipendiefond. JJ acknowledges support by the Swedish Research Council (VR) under the project 2020-05143 – ‘Deciphering the Dynamics of Cosmic Structure’. GL acknowledges support by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. FS acknowledges support from the Starting Grant (ERC-2015-STG 678652) ‘GrInflaGal’ of the European Research Council. This work was supported by the Simons Collaboration on ‘Learning the Universe’.

The computation and data processing in this study were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Tetralith, partially funded by the Swedish Research Council through grant agreement no. 2020-05143. This research utilized the HPC facility supported by the Technical Division of the Department of Physics, Stockholm University. This work is done within the Aquila Consortium.<fnref n="2" /> We acknowledge the use of the following packages: NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), GETDIST (Lewis 2019), and HEALPIX (Gorski et al. 2005).

## DATA AVAILABILITY

The data underlying this article will be shared on the basis of a reasonable request to the corresponding author.

## REFERENCES

* Alam S. et al., 2015, ApJS, 219, 12

* Alam S. et al., 2017, MNRAS, 470, 2617

* Alvarez M. et al., 2014, preprint (arXiv:1412.4671)

* Amendola L. et al., 2018, Living Rev. Relativ., 21, 2

* Assassi V., Baumann D., Schmidt F., 2015, J. Cosmol. Astropart. Phys., 2015, 043

* Baldauf T., Seljak U., Senatore L., 2011, J. Cosmol. Astropart. Phys., 2011, 006

* Ballardini M., Matthewson W. L., Maartens R., 2019, MNRAS, 489, 1950

* Barnaby N., Shandera S., 2012, J. Cosmol. Astropart. Phys., 2012, 034

* Barreira A., 2020, J. Cosmol. Astropart. Phys., 2020, 031

* Barreira A., 2021a, J. Cosmol. Astropart. Phys., 2022, 057

* Barreira A., 2021b, J. Cosmol. Astropart. Phys., 2022, 033

* Barreira A., Cabass G., Schmidt F., Pillepich A., Nelson D., 2020, J. Cosmol. Astropart. Phys., 2020, 013

* Baumann D., Green D., 2021, J. Cosmol. Astropart. Phys., 2022, 061

* Betancourt M., 2017, A Conceptual Introduction to Hamiltonian Monte Carlo

* Bharadwaj S., Mazumdar A., Sarkar D., 2020, MNRAS, 493, 594

* Biagetti M., 2019, Galaxies, 7, 71

* Biagetti M., Lazeyras T., Baldauf T., Desjacques V., Schmidt F., 2017, MNRAS, 468, 3277

* Bond J. R., Efstathiou G., 1984, ApJ, 285, L45

* Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, A&A, 296, 575

* Buchert T., Melott A. L., Weiss A. G., 1994, A&A, 288, 349

* Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022, Phys. Rev. Lett., 129, 21301

* Camera S., Santos M. G., Ferreira P. G., Ferramacho L., 2013, Phys. Rev. Lett., 111, 171302

* Castorina E. et al., 2019, J. Cosmol. Astropart. Phys., 2019, 010

* Catelan P., Scherrer R. J., 1995, ApJ, 445, 1

* Celoria M., Matarrese S., 2018, preprint (arXiv:1812.08197)

* Charnock T., Lavaux G., Wandelt B. D., Boruah S. S., Jasje J., Hudson M. J., 2020, MNRAS, 494, 50

* Chen X., 2010, A&A, 2010, 638979

* Chen X., Wang Y., 2010, Phys. Rev. D, 81, 63511

<footnote n="2">²https://www.aquila-consortium.org/</footnote>
Chodorowski M. J., Bouchet F. R., 1996, MNRAS, 279, 557

Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369

D'Amico G., Lewandowski M., Senatore L., Zhang P., 2022, preprint (arXiv:2201.11518)

Dalal N., Doré O., Huterer D., Shirokov A., 2008, Phys. Rev. D, 77, 123514

de Putter R., Doré O., 2017, Phys. Rev. D, 95, 123513

Desjacques V., Seljak U., Iliev I. T., 2009, MNRAS, 396, 85

Desjacques V., Jeong D., Schmidt F., 2018, Phys. Rep., 733, 1

Di Dio E., Perrier H., Durrer R., Marozzi G., Dizgah A. M., Noreña J., Riotto A., 2017, J. Cosmol. Astropart. Phys., 2017, 006

Doré O. et al., 2014, preprint (arXiv:1412.4872)

Duane S., Kennedy A. D., Pendleton B. J., Roweth D., 1987, Phys. Lett. B, 195, 216

Durrer R., Juszkiewicz R., Kunz M., Uzan J.-P., 2000, Phys. Rev. D, 62, 021301

Eisenstein D. J., Hu W., 1998, ApJ, 496, 605

Eisenstein D. J., Hu W., 1999, ApJ, 511, 5

Euclid Collaboration, 2020, A&A, 642, A191

Friedrich O., Uhlemann C., Villaescusa-Navarro F., Baldauf T., Manera M., Nishimichi T., 2020, MNRAS, 498, 464

Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

Gorski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelman M., 2005, ApJ, 622, 759

Graham M. L., Connolly A. J., Ivezić v., Schmidt S. J., Jones R. L., Jurić M., Daniel S. F., Yoachim P., 2018, AJ, 155, 1

Gunn J. E., Gott J., Richard I., 1972, ApJ, 176, 1

Harris C. R. et al., 2020, Nature, 585, 357

Hastings W. K., 1970, Biometrika, 57, 97

Ho S. et al., 2015, J. Cosmol. Astropart. Phys., 2015, 040

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Huterer D., Cunha C. E., Fang W., 2013, MNRAS, 432, 2945

Jasche J., Kitaura F. S., 2010, MNRAS, 407, 29

Jasche J., Lavaux G., 2015, MNRAS, 447, 1204

Jasche J., Lavaux G., 2017, A&A, 606, A37

Jasche J., Lavaux G., 2019, A&A, 625, A64

Jasche J., Wandelt B. D., 2012, MNRAS, 425, 1042

Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894

Jasche J., Kitaura F. S., Wandelt B. D., Enßlin T. A., 2010, MNRAS, 406, 60

Jasche J., Leclercq F., Wandelt B. D., 2015, J. Cosmol. Astropart. Phys., 2015, 036

Jeong D., Schmidt F., 2015, Class. Quantum Gravity, 32, 044001

Kaiser N., 1984, ApJ, 284, L9

Karagiannis D., Shanks T., Ross N. P., 2014, MNRAS, 441, 486

Karagiannis D., Lazanu A., Liguori M., Raccanelli A., Bartolo N., Verde L., 2018, MNRAS, 478, 1341

Karagiannis D., Slosar A., Liguori M., 2020, J. Cosmol. Astropart. Phys., 2020, 052

Kitaura F. S., Enßlin T. A., 2008, MNRAS, 389, 497

Komatsu E., 2010, Class. Quantum Gravity, 27, 124010

Komatsu E., Spergel D. N., 2001, Phys. Rev. D, 63, 063002

Kostić A., Jasche J., Ramanah D. K., Lavaux G., 2022, A&A, 657, L17

LSST Science Collaboration, 2009, preprint (arXiv:0912.0201)

Lam T. Y., Nishimichi T., Yoshida N., 2011, MNRAS, 414, 289

Lavaux G., Jasche J., 2016, MNRAS, 455, 3169

Lavaux G., Jasche J., Leclercq F., 2019, preprint (arXiv:1909.06396)

Layzer D., 1956, AJ, 61, 383

Leclercq F., Heavens A., 2021, MNRAS, 506, L85

Leistedt B., Peiris H. V., Roth N., 2014, Phys. Rev. Lett., 113, 221301

Lesgourgues J., Tram T., 2014, J. Cosmol. Astropart. Phys., 2014, 032

Lewis A., 2019, preprint (arXiv:1910.13970)

Ma Y.-Z., Taylor J. E., Scott D., 2013, MNRAS, 436, 2029

Maldacena J., 2003, J. High Energy Phys., 2003, 013

Matarrese S., Verde L., 2008, ApJ, 677, L77

McQuinn M., 2021, J. Cosmol. Astropart. Phys., 2021, 024

Meerburg P. D. et al., 2019, BAAS, 51, 107

Moradinezhad Dizgah A., Keating G. K., 2019, ApJ, 872, 126

Moradinezhad Dizgah A., Biagetti M., Sefusatti E., Desjacques V., Noreña J., 2021, J. Cosmol. Astropart. Phys., 2021, 015

Moutarde F., Alimi J.-M., Bouchet F. R., Pellat R., Ramani A., 1991, ApJ, 382, 377

Mueller E.-M., Percival W. J., Ruggeri R., 2019, MNRAS, 485, 4160

Mueller E.-M. et al., 2021, preprint (arXiv:2106.13725)

Neal R. M., 2000, J. Comput. Graph., 9, 249

Neal R. M., 2003, Ann. Stat., 31, 705

Neyrinck M. C., Aragón-Calvo M. A., Jeong D., Wang X., 2014, MNRAS, 441, 646

Nguyen N.-M., Jasche J., Lavaux G., Schmidt F., 2020, J. Cosmol. Astropart. Phys., 2020, 011

Nguyen N.-M., Schmidt F., Lavaux G., Jasche J., 2021, J. Cosmol. Astropart. Phys., 2021, 058

Peebles P. J. E., 1980, Princeton Univ. Press, Princeton, NJ

Pillepich A., Porciani C., Hahn O., 2009, MNRAS, 402, 191

Planck Collaboration IX, 2019a, A&A, 641, A9

Planck Collaboration X, 2019b, A&A, 641, A10

Planck Collaboration VI, 2020, A&A, 641, A6

Porqueres N., Kodi Ramanah D., Jasche J., Lavaux G., 2019, A&A, 624, A115

Porqueres N., Hahn O., Jasche J., Lavaux G., 2020, A&A, 642, A139

Press W. H., Schechter P., 1974, ApJ, 187, 425

Ramanah D. K., Lavaux G., Jasche J., Wandelt B. D., 2019, A&A, 621, A69

Ross A. J. et al., 2012, MNRAS, 428, 1116

Ross A. J. et al., 2017, MNRAS, 464, 1168

Salopek D. S., Bond J. R., 1990, Phys. Rev. D, 42, 3936

Schmidt F., 2010, Phys. Rev. D, 82, 063001

Schmidt F., 2021, J. Cosmol. Astropart. Phys., 2021, 033

Schmidt F., Jeong D., Desjacques V., 2013, Phys. Rev. D, 88, 023515

Schmidt F., Elsner F., Jasche J., Nguyen N. M., Lavaux G., 2019, J. Cosmol. Astropart. Phys., 2019, 042

Schmidt F., Cabass G., Jasche J., Lavaux G., 2020, J. Cosmol. Astropart. Phys., 2020, 008

Schmittfull M., Simonović M., Assassini V., Zaldarriaga M., 2019, Phys. Rev. D, 100, 043514

Scoccimarro R., 2000, ApJ, 544, 597

Scoccimarro R., Sheth R. K., 2002, MNRAS, 329, 629

Seljak U., 2009, Phys. Rev. Lett., 102, 021302

Seljak U., Aslanyan G., Feng Y., Modi C., 2017, J. Cosmol. Astropart. Phys., 2017, 009

Senatore L., Smith K. M., Zaldarriaga M., 2010, J. Cosmol. Astropart. Phys., 2010, 028

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Slosar A., Hirata C., Seljak U., Ho S., Padmanabhan N., 2008, J. Cosmol. Astropart. Phys., 2008, 031

Takada M., Komatsu E., Futamase T., 2006, Phys. Rev. D, 73, 083520

Tasinato G., Tellarini M., Ross A. J., Wands D., 2014, J. Cosmol. Astropart. Phys., 2014, 032

Tellarini M., Ross A. J., Tasinato G., Wands D., 2016, J. Cosmol. Astropart. Phys., 2016, 014

Uhlemann C., Pajer E., Pichon C., Nishimichi T., Codis S., Bernardeau F., 2018, MNRAS, 474, 2853

Verde L., Wang L., Heavens A. F., Kamionkowski M., 2000, MNRAS, 313, 141

Wang H., Mo H. J., Yang X., Jing Y. P., Lin W. P., 2014, ApJ, 794, 94

Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435

Yokoyama S., Sugiyama N., Zaroubi S., Silk J., 2011, MNRAS, 417, 1074

## APPENDIX A: TESTING THE INFERENCE IN THE LINEAR REGIME

In this section, we outline a linearized test of the method. Our aim is to test whether the method is capable of retrieving $f_{nl}$ given data generated at a scale that can be considered linear. The run specifications are as follows: $L = 4096 h^{-1}$ Mpc, $N = 64$,
Chodorowski M. J., Bouchet F. R., 1996, MNRAS, 279, 557

Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369

D'Amico G., Lewandowski M., Senatore L., Zhang P., 2022, preprint (arXiv:2201.11518)

Dalal N., Doré O., Huterer D., Shirokov A., 2008, Phys. Rev. D, 77, 123514

de Putter R., Doré O., 2017, Phys. Rev. D, 95, 123513

Desjacques V., Seljak U., Iliev I. T., 2009, MNRAS, 396, 85

Desjacques V., Jeong D., Schmidt F., 2018, Phys. Rep., 733, 1

Di Dio E., Perrier H., Durrer R., Marozzi G., Dizgah A. M., Noreña J., Riotto A., 2017, J. Cosmol. Astropart. Phys., 2017, 006

Doré O. et al., 2014, preprint (arXiv:1412.4872)

Duane S., Kennedy A. D., Pendleton B. J., Roweth D., 1987, Phys. Lett. B, 195, 216

Durrer R., Juszkiewicz R., Kunz M., Uzan J.-P., 2000, Phys. Rev. D, 62, 021301

Eisenstein D. J., Hu W., 1998, ApJ, 496, 605

Eisenstein D. J., Hu W., 1999, ApJ, 511, 5

Euclid Collaboration, 2020, A&A, 642, A191

Friedrich O., Uhlemann C., Villaescusa-Navarro F., Baldauf T., Manera M., Nishimichi T., 2020, MNRAS, 498, 464

Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

Gorski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelman M., 2005, ApJ, 622, 759

Graham M. L., Connolly A. J., Ivezić v., Schmidt S. J., Jones R. L., Jurić M., Daniel S. F., Yoachim P., 2018, AJ, 155, 1

Gunn J. E., Gott J., Richard I., 1972, ApJ, 176, 1

Harris C. R. et al., 2020, Nature, 585, 357

Hastings W. K., 1970, Biometrika, 57, 97

Ho S. et al., 2015, J. Cosmol. Astropart. Phys., 2015, 040

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Huterer D., Cunha C. E., Fang W., 2013, MNRAS, 432, 2945

Jasche J., Kitaura F. S., 2010, MNRAS, 407, 29

Jasche J., Lavaux G., 2015, MNRAS, 447, 1204

Jasche J., Lavaux G., 2017, A&A, 606, A37

Jasche J., Lavaux G., 2019, A&A, 625, A64

Jasche J., Wandelt B. D., 2012, MNRAS, 425, 1042

Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894

Jasche J., Kitaura F. S., Wandelt B. D., Enßlin T. A., 2010, MNRAS, 406, 60

Jasche J., Leclercq F., Wandelt B. D., 2015, J. Cosmol. Astropart. Phys., 2015, 036

Jeong D., Schmidt F., 2015, Class. Quantum Gravity, 32, 044001

Kaiser N., 1984, ApJ, 284, L9

Karagiannis D., Shanks T., Ross N. P., 2014, MNRAS, 441, 486

Karagiannis D., Lazanu A., Liguori M., Raccanelli A., Bartolo N., Verde L., 2018, MNRAS, 478, 1341

Karagiannis D., Slosar A., Liguori M., 2020, J. Cosmol. Astropart. Phys., 2020, 052

Kitaura F. S., Enßlin T. A., 2008, MNRAS, 389, 497

Komatsu E., 2010, Class. Quantum Gravity, 27, 124010

Komatsu E., Spergel D. N., 2001, Phys. Rev. D, 63, 063002

Kostić A., Jasche J., Ramanah D. K., Lavaux G., 2022, A&A, 657, L17

LSST Science Collaboration, 2009, preprint (arXiv:0912.0201)

Lam T. Y., Nishimichi T., Yoshida N., 2011, MNRAS, 414, 289

Lavaux G., Jasche J., 2016, MNRAS, 455, 3169

Lavaux G., Jasche J., Leclercq F., 2019, preprint (arXiv:1909.06396)

Layzer D., 1956, AJ, 61, 383

Leclercq F., Heavens A., 2021, MNRAS, 506, L85

Leistedt B., Peiris H. V., Roth N., 2014, Phys. Rev. Lett., 113, 221301

Lesgourgues J., Tram T., 2014, J. Cosmol. Astropart. Phys., 2014, 032

Lewis A., 2019, preprint (arXiv:1910.13970)

Ma Y.-Z., Taylor J. E., Scott D., 2013, MNRAS, 436, 2029

Maldacena J., 2003, J. High Energy Phys., 2003, 013

Matarrese S., Verde L., 2008, ApJ, 677, L77

McQuinn M., 2021, J. Cosmol. Astropart. Phys., 2021, 024

Meerburg P. D. et al., 2019, BAAS, 51, 107

Moradinezhad Dizgah A., Keating G. K., 2019, ApJ, 872, 126

Moradinezhad Dizgah A., Biagetti M., Sefusatti E., Desjacques V., Noreña J., 2021, J. Cosmol. Astropart. Phys., 2021, 015

Moutarde F., Alimi J.-M., Bouchet F. R., Pellat R., Ramani A., 1991, ApJ, 382, 377

Mueller E.-M., Percival W. J., Ruggeri R., 2019, MNRAS, 485, 4160

Mueller E.-M. et al., 2021, preprint (arXiv:2106.13725)

Neal R. M., 2000, J. Comput. Graph., 9, 249

Neal R. M., 2003, Ann. Stat., 31, 705

Neyrinck M. C., Aragón-Calvo M. A., Jeong D., Wang X., 2014, MNRAS, 441, 646

Nguyen N.-M., Jasche J., Lavaux G., Schmidt F., 2020, J. Cosmol. Astropart. Phys., 2020, 011

Nguyen N.-M., Schmidt F., Lavaux G., Jasche J., 2021, J. Cosmol. Astropart. Phys., 2021, 058

Peebles P. J. E., 1980, Princeton Univ. Press, Princeton, NJ

Pillepich A., Porciani C., Hahn O., 2009, MNRAS, 402, 191

Planck Collaboration IX, 2019a, A&A, 641, A9

Planck Collaboration X, 2019b, A&A, 641, A10

Planck Collaboration VI, 2020, A&A, 641, A6

Porqueres N., Kodi Ramanah D., Jasche J., Lavaux G., 2019, A&A, 624, A115

Porqueres N., Hahn O., Jasche J., Lavaux G., 2020, A&A, 642, A139

Press W. H., Schechter P., 1974, ApJ, 187, 425

Ramanah D. K., Lavaux G., Jasche J., Wandelt B. D., 2019, A&A, 621, A69

Ross A. J. et al., 2012, MNRAS, 428, 1116

Ross A. J. et al., 2017, MNRAS, 464, 1168

Salopek D. S., Bond J. R., 1990, Phys. Rev. D, 42, 3936

Schmidt F., 2010, Phys. Rev. D, 82, 063001

Schmidt F., 2021, J. Cosmol. Astropart. Phys., 2021, 033

Schmidt F., Jeong D., Desjacques V., 2013, Phys. Rev. D, 88, 023515

Schmidt F., Elsner F., Jasche J., Nguyen N. M., Lavaux G., 2019, J. Cosmol. Astropart. Phys., 2019, 042

Schmidt F., Cabass G., Jasche J., Lavaux G., 2020, J. Cosmol. Astropart. Phys., 2020, 008

Schmittfull M., Simonović M., Assassini V., Zaldarriaga M., 2019, Phys. Rev. D, 100, 043514

Scoccimarro R., 2000, ApJ, 544, 597

Scoccimarro R., Sheth R. K., 2002, MNRAS, 329, 629

Seljak U., 2009, Phys. Rev. Lett., 102, 021302

Seljak U., Aslanyan G., Feng Y., Modi C., 2017, J. Cosmol. Astropart. Phys., 2017, 009

Senatore L., Smith K. M., Zaldarriaga M., 2010, J. Cosmol. Astropart. Phys., 2010, 028

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Slosar A., Hirata C., Seljak U., Ho S., Padmanabhan N., 2008, J. Cosmol. Astropart. Phys., 2008, 031

Takada M., Komatsu E., Futamase T., 2006, Phys. Rev. D, 73, 083520

Tasinato G., Tellarini M., Ross A. J., Wands D., 2014, J. Cosmol. Astropart. Phys., 2014, 032

Tellarini M., Ross A. J., Tasinato G., Wands D., 2016, J. Cosmol. Astropart. Phys., 2016, 014

Uhlemann C., Pajer E., Pichon C., Nishimichi T., Codis S., Bernardeau F., 2018, MNRAS, 474, 2853

Verde L., Wang L., Heavens A. F., Kamionkowski M., 2000, MNRAS, 313, 141

Wang H., Mo H. J., Yang X., Jing Y. P., Lin W. P., 2014, ApJ, 794, 94

Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435

Yokoyama S., Sugiyama N., Zaroubi S., Silk J., 2011, MNRAS, 417, 1074

## APPENDIX A: TESTING THE INFERENCE IN THE LINEAR REGIME

In this section, we outline a linearized test of the method. Our aim is to test whether the method is capable of retrieving $f_{nl}$ given data generated at a scale that can be considered linear. The run specifications are as follows: $L = 4096 h^{-1}$ Mpc, $N = 64$,
<figure><img src="image_11.png" /><figcaption>**Figure A1.** Cross-correlations and marginal posteriors of all hyperparameters and $f_{nl}$ for the linearized run. The darker regions in the contour plots represent the $1\sigma$ interval, and the lighter regions represent the $2\sigma$ interval. The blue lines mark the fiducial values injected into the mock data, and constitute the true values of the data. The results indicate that our method is able to infer $f_{nl}$ under a simplified, linearized density field. The test is mainly intended to be a consistency check.</figcaption></figure>

$\Delta L = 64h^{-1} \text{Mpc}$, $k_{\text{max}} \approx 0.1$, $\bar{n}_{\text{gal}} = 3.04 \times 10^{-4}h^3\text{Mpc}^{-3}$. The forward model and cosmology are the same as for the runs (see Sections 2 and 3. Galaxies are generated with $b_1 = 2.5$. For simplicity, neither a sky survey mask nor a radial selection has been applied. The resulting inference is presented in Fig. A1, in the form of a corner plot containing bias parameters and $f_{nl}$. As can be seen, the method is successful in inferring the correct fiducial value of $f_{nl} = 5$.

## APPENDIX B: ADDITIONAL RESULTS

### B1 $f_{nl}$ posterior distributions of runs

In Fig. B1, we have illustrated the posterior distributions of $f_{nl}$ for the secondary runs of this paper. Although their uncertainties have been provided in Table 3 and Fig. 7, here we provide the full inferred posterior distributions. Note that these distributions include any and all survey effects, e.g. survey geometry, instrumentation noise, and galaxy biases, except for run #6. The posterior distributions have their centres within $1\sigma$ of the fiducial values ($f_{nl}^{\text{fid.}} = 5$), indicating that our method can consistently infer the correct $f_{nl}$ distribution.

### B2 Correlation matrix

The correlation matrix of the high-resolution *Stage IV* mock data run can be found in Fig. B2. The figure displays the correlations and anticorrelations between model parameters throughout the samples in the chain. As can be seen, there is a very weak correlation between $f_{nl}$ and the bias parameters of the different tracer catalogues. This is a very promising fact, since it indicates that there is little-to-no relationship directly between $f_{nl}$ and the hyperparameters of the galaxy bias model. This implies that the cosmological parameter

<figure><img src="image_12.png" /><figcaption>(a) The marginalised posterior distribution of $f_{nl}$ for run #1; the low resolution SDSS-III/BOSS-like run.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) The marginalised posterior distribution of $f_{nl}$ for run #2; the mid resolution SDSS-III/BOSS-like run.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) The marginalised posterior distribution of $f_{nl}$ for run #3; the high resolution SDSS-III/BOSS-like run.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) The marginalised posterior distribution of $f_{nl}$ for run #4; the low resolution *Stage IV* run.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(e) The marginalised posterior distribution of $f_{nl}$ for run #6; the high resolution *Stage IV* run, with fixed bias sampling.</figcaption></figure>

<caption>**Figure B1.** The inferred posterior distributions of $f_{nl}$ for all other runs. These distributions include handling and marginalizing out all other effects, e.g. survey geometry, instrumentation noise, and galaxy biases, with the exception of run #6.</caption>

is being sampled and detected without strong dependencies on the bias parameters, which strengthens the validity of our method. Moreover, these results further suggest, together with the contour plot figure (Fig. 10), that the model has very few degeneracies between the model parameters.
<figure><img src="image_13.png" /><figcaption>**Figure A1.** Cross-correlations and marginal posteriors of all hyperparameters and $f_{nl}$ for the linearized run. The darker regions in the contour plots represent the $1\sigma$ interval, and the lighter regions represent the $2\sigma$ interval. The blue lines mark the fiducial values injected into the mock data, and constitute the true values of the data. The results indicate that our method is able to infer $f_{nl}$ under a simplified, linearized density field. The test is mainly intended to be a consistency check.</figcaption></figure>

$\Delta L = 64h^{-1} \text{Mpc}$, $k_{\text{max}} \approx 0.1$, $\bar{n}_{\text{gal}} = 3.04 \times 10^{-4}h^3\text{Mpc}^{-3}$. The forward model and cosmology are the same as for the runs (see Sections 2 and 3. Galaxies are generated with $b_1 = 2.5$. For simplicity, neither a sky survey mask nor a radial selection has been applied. The resulting inference is presented in Fig. A1, in the form of a corner plot containing bias parameters and $f_{nl}$. As can be seen, the method is successful in inferring the correct fiducial value of $f_{nl} = 5$.

## APPENDIX B: ADDITIONAL RESULTS

### B1 $f_{nl}$ posterior distributions of runs

In Fig. B1, we have illustrated the posterior distributions of $f_{nl}$ for the secondary runs of this paper. Although their uncertainties have been provided in Table 3 and Fig. 7, here we provide the full inferred posterior distributions. Note that these distributions include any and all survey effects, e.g. survey geometry, instrumentation noise, and galaxy biases, except for run #6. The posterior distributions have their centres within $1\sigma$ of the fiducial values ($f_{nl}^{\text{fid.}} = 5$), indicating that our method can consistently infer the correct $f_{nl}$ distribution.

### B2 Correlation matrix

The correlation matrix of the high-resolution *Stage IV* mock data run can be found in Fig. B2. The figure displays the correlations and anticorrelations between model parameters throughout the samples in the chain. As can be seen, there is a very weak correlation between $f_{nl}$ and the bias parameters of the different tracer catalogues. This is a very promising fact, since it indicates that there is little-to-no relationship directly between $f_{nl}$ and the hyperparameters of the galaxy bias model. This implies that the cosmological parameter

<figure><img src="image_14.png" /><figcaption>(a) The marginalised posterior distribution of $f_{nl}$ for run #1; the low resolution SDSS-III/BOSS-like run.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) The marginalised posterior distribution of $f_{nl}$ for run #2; the mid resolution SDSS-III/BOSS-like run.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) The marginalised posterior distribution of $f_{nl}$ for run #3; the high resolution SDSS-III/BOSS-like run.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) The marginalised posterior distribution of $f_{nl}$ for run #4; the low resolution *Stage IV* run.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(e) The marginalised posterior distribution of $f_{nl}$ for run #6; the high resolution *Stage IV* run, with fixed bias sampling.</figcaption></figure>

<caption>**Figure B1.** The inferred posterior distributions of $f_{nl}$ for all other runs. These distributions include handling and marginalizing out all other effects, e.g. survey geometry, instrumentation noise, and galaxy biases, with the exception of run #6.</caption>

is being sampled and detected without strong dependencies on the bias parameters, which strengthens the validity of our method. Moreover, these results further suggest, together with the contour plot figure (Fig. 10), that the model has very few degeneracies between the model parameters.
<figure><img src="image_17.png" /><figcaption>Figure B2. The plot shows the correlation matrix of run #3. The correlation matrix displays any internal correlations between model parameters. Most importantly, any correlations between $f_{nl}$ and the galaxy bias parameters are undesired, since they would indicate that the inferred $f_{nl}$ may be biased by the model parameters. The absence of these features indicate that we can robustly infer the correct $f_{nl}$ value with $1\sigma$ confidence.</figcaption></figure>

## B3 Correlation lengths

The correlation lengths of the high-resolution Stage IV runs can be found in Figs F2 and F1. This plot shows the correlation of the $f_{nl}$ parameter between the chains' samples. This can be described as

$$C(f_{nl})_n = \frac{1}{N-n} \sum_{i=0}^{N-n} \left( \frac{f_{nl}^i - \langle f_{nl} \rangle}{\sigma_{f_{nl}}} \frac{f_{nl}^{i+n} - \langle f_{nl} \rangle}{\sigma_{f_{nl}}} \right), \quad (B1)$$

where $n$ is the number of transition steps, $\langle f_{nl} \rangle$ and $\sigma_{f_{nl}}$ are the mean and standard deviation of the sampled $f_{nl}$ value, and $N$ is the total number of samples in the run. The plot shows that BORG is capable of sampling statistically independent samples roughly after $\approx 1000$ samples. These results demonstrate that the statistical efficiency of our $f_{nl}$ sampler is robust, since the number of steps between independent samples is much smaller than the total length of the chain.

## APPENDIX C: SURVEY CONFIGURATIONS FOR THE STAGE IV SURVEY

In this section, we provide plots of the radial selection functions and the sky map for the Stage IV survey. The design choices for these specifications are outlined in Section 3.3. The radial selection functions can be seen in Fig. F4. The sky map can be seen in Fig. F3.

## APPENDIX D: GELMAN-RUBIN TEST

In this section, we describe the Gelman-Rubin test, which evaluates the convergence of a MCMC chain (Gelman & Rubin 1992). The process in which it does so is the following: First, the chain is divided into $M$ different subchains, each containing $n$ samples, with $d$ inbetween (which are ignored in the analysis). Then, we compute the

potential scale reduction factor:

$$\text{PSRF} = \sqrt{\frac{V}{W}}, \quad (D1)$$

where

(i) $V = \frac{n-1}{n}W + \frac{M+1}{nM}$,

(ii) $W = \frac{1}{M} \sum_{m=1}^{M} M \sigma_m^2$,

(iii) $B = \frac{n}{M-1} \sum_{m=1}^{M} (\langle f_{nl} \rangle_m - \langle f_{nl} \rangle)^2$.

with $\sigma_m^2$ being the variance of $f_{nl}$ for the $m$th chain. The chain is considered to be converged as $\text{PSRF} \to 1$.

For run #5, we achieve the value of $V = 1.000021$, when $M = 2$, $n = 25000$, and $d = 6100$. Thus, this value indicates that the chain has converged (Gelman & Rubin 1992). Similar findings are made for the other runs, as well.

## APPENDIX E: ADJOINT GRADIENT CALCULATION OF THE $f_{NL}$ PERTURBATION

In this section, we outline the adjoint gradient of the $f_{nl}$ perturbation module, which is used by the algorithm to evaluate the sensitivity of model parameters. In this case, we want to evaluate the change in the non-Gaussian primordial density field $\delta^{\text{NG}}$, as the input Gaussian density field $\delta^{\text{G}}$ is altered. First, we drop the $m$-index for the matter density field $\delta_m$, such that $\delta_m \equiv \delta$, for the sake of brevity. Then, we start in the general case:

$$\frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_j^{\text{G}}} = \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}}. \quad (E1)$$

Next, we evaluate $\frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}}$. But in order to do so, we must first decompose $\delta_i^{\text{NG}}$:

$$\begin{align}
\delta_i^{\text{NG}} &= \sum_a \bar{F}_{i,a} \Phi_a, \nonumber \\
\Phi_a &= \phi_a + f_{nl} \phi_a^2, \nonumber \\
\phi_a &= \sum_b \mathcal{F}_{a,b} \delta_b^{\text{G}}, \tag{E2}
\end{align}$$

where $\mathcal{F}$ is the Fourier transform, $\Phi$ is the perturbed gravitational potential, and $\phi$ is the non-perturbed gravitational potential, both potentials transformed to redshift $z = 1000$, and acted under the correct transfer functions.

Now, we can evaluate $\frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}}$:

$$\frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}} = \frac{\partial \delta_i^{\text{NG}}}{\partial \Phi_k} \frac{\partial \Phi_k}{\partial \phi_l} \frac{\partial \phi_l}{\partial \delta_j^{\text{G}}}. \quad (E3)$$

The expressions for these partial derivatives are as follows:

$$\begin{align}
\frac{\partial \delta_i^{\text{NG}}}{\partial \Phi_k} &= \sum_b \bar{F}_{i,b} \delta_{b,k}^K, \nonumber \\
\frac{\partial \Phi_k}{\partial \phi_l} &= \delta_{k,l}^K + 2 \delta_{k,l}^K f_{nl} \phi_k, \nonumber \\
\frac{\partial \phi_l}{\partial \delta_j^{\text{G}}} &= \sum_q \mathcal{F}_{l,q} \delta_{l,j}^K, \tag{E4}
\end{align}$$

where $\delta^K$ is the Kronecker delta.
<figure><img src="image_19.png" /><figcaption>**Figure B2.** The plot shows the correlation matrix of run #3. The correlation matrix displays any internal correlations between model parameters. Most importantly, any correlations between $f_{nl}$ and the galaxy bias parameters are undesired, since they would indicate that the inferred $f_{nl}$ may be biased by the model parameters. The absence of these features indicate that we can robustly infer the correct $f_{nl}$ value with $1\sigma$ confidence.</figcaption></figure>

## B3 Correlation lengths

The correlation lengths of the high-resolution Stage IV runs can be found in Figs F2 and F1. This plot shows the correlation of the $f_{nl}$ parameter between the chains' samples. This can be described as

$$C(f_{nl})_n = \frac{1}{N-n} \sum_{i=0}^{N-n} \left( \frac{f_{nl}^i - \langle f_{nl} \rangle}{\sigma_{f_{nl}}} \frac{f_{nl}^{i+n} - \langle f_{nl} \rangle}{\sigma_{f_{nl}}} \right), \quad (B1)$$

where $n$ is the number of transition steps, $\langle f_{nl} \rangle$ and $\sigma_{f_{nl}}$ are the mean and standard deviation of the sampled $f_{nl}$ value, and $N$ is the total number of samples in the run. The plot shows that BORG is capable of sampling statistically independent samples roughly after $\approx 1000$ samples. These results demonstrate that the statistical efficiency of our $f_{nl}$ sampler is robust, since the number of steps between independent samples is much smaller than the total length of the chain.

## APPENDIX C: SURVEY CONFIGURATIONS FOR THE STAGE IV SURVEY

In this section, we provide plots of the radial selection functions and the sky map for the Stage IV survey. The design choices for these specifications are outlined in Section 3.3. The radial selection functions can be seen in Fig. F4. The sky map can be seen in Fig. F3.

## APPENDIX D: GELMAN-RUBIN TEST

In this section, we describe the Gelman-Rubin test, which evaluates the convergence of a MCMC chain (Gelman & Rubin 1992). The process in which it does so is the following: First, the chain is divided into $M$ different subchains, each containing $n$ samples, with $d$ inbetween (which are ignored in the analysis). Then, we compute the

potential scale reduction factor:

$$\text{PSRF} = \sqrt{\frac{V}{W}}, \quad (D1)$$

where

(i) $V = \frac{n-1}{n}W + \frac{M+1}{nM}$,

(ii) $W = \frac{1}{M} \sum_{m=1}^{M} M \sigma_m^2$,

(iii) $B = \frac{n}{M-1} \sum_{m=1}^{M} (\langle f_{nl} \rangle_m - \langle f_{nl} \rangle)^2$.

with $\sigma_m^2$ being the variance of $f_{nl}$ for the $m$th chain. The chain is considered to be converged as $\text{PSRF} \to 1$.

For run #5, we achieve the value of $V = 1.000021$, when $M = 2$, $n = 25000$, and $d = 6100$. Thus, this value indicates that the chain has converged (Gelman & Rubin 1992). Similar findings are made for the other runs, as well.

## APPENDIX E: ADJOINT GRADIENT CALCULATION OF THE $f_{NL}$ PERTURBATION

In this section, we outline the adjoint gradient of the $f_{nl}$ perturbation module, which is used by the algorithm to evaluate the sensitivity of model parameters. In this case, we want to evaluate the change in the non-Gaussian primordial density field $\delta^{\text{NG}}$, as the input Gaussian density field $\delta^{\text{G}}$ is altered. First, we drop the $m$-index for the matter density field $\delta_m$, such that $\delta_m \equiv \delta$, for the sake of brevity. Then, we start in the general case:

$$\frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_j^{\text{G}}} = \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}}. \quad (E1)$$

Next, we evaluate $\frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}}$. But in order to do so, we must first decompose $\delta_i^{\text{NG}}$:

$$\begin{aligned} \delta_i^{\text{NG}} &= \sum_a \bar{F}_{i,a} \Phi_a, \\ \Phi_a &= \phi_a + f_{nl} \phi_a^2, \\ \phi_a &= \sum_b \mathcal{F}_{a,b} \delta_b^{\text{G}}, \end{aligned} \quad (E2)$$

where $\mathcal{F}$ is the Fourier transform, $\Phi$ is the perturbed gravitational potential, and $\phi$ is the non-perturbed gravitational potential, both potentials transformed to redshift $z = 1000$, and acted under the correct transfer functions.

Now, we can evaluate $\frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}}$:

$$\frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^{\text{G}}} = \frac{\partial \delta_i^{\text{NG}}}{\partial \Phi_k} \frac{\partial \Phi_k}{\partial \phi_l} \frac{\partial \phi_l}{\partial \delta_j^{\text{G}}}. \quad (E3)$$

The expressions for these partial derivatives are as follows:

$$\begin{aligned} \frac{\partial \delta_i^{\text{NG}}}{\partial \Phi_k} &= \sum_b \bar{F}_{i,b} \delta_{b,k}^K, \\ \frac{\partial \Phi_k}{\partial \phi_l} &= \delta_{k,l}^K + 2 \delta_{k,l}^K f_{nl} \phi_k, \\ \frac{\partial \phi_l}{\partial \delta_j^{\text{G}}} &= \sum_q \mathcal{F}_{l,q} \delta_{l,j}^K, \end{aligned} \quad (E4)$$

where $\delta^K$ is the Kronecker delta.
Now, combining yields the full expression of the adjoint gradient. The log-posterior distribution is therefore

$$ \begin{aligned} \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_j^G} &= \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^G} \\ &= \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \frac{\partial \delta_i^{\text{NG}}}{\partial \Phi_k} \frac{\partial \Phi_k}{\partial \phi_l} \frac{\partial \phi_l}{\partial \delta_j^G} \\ &= \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \left[ \sum_b \bar{\mathcal{F}}_{i,b} \delta_{b,k}^K \right] [\delta_{k,l} + \delta_{k,l}^K f_{nl} \phi_k] \\ & \quad \times \left[ \sum_q \mathcal{F}_{1,q} \delta_{1,j}^K \right]. \end{aligned} \quad (E5) $$

## APPENDIX F: ADJOINT GRADIENT COMPUTATION OF BIAS MODEL

As in the previous section, we explicitly describe the adjoint gradient of the bias model. First, we drop the *m*-index for the matter density field $\delta_m$, such that $\delta_m \equiv \delta$, for the sake of brevity. Then, we perform the derivative of the log likelihood with respect to $\delta$:

$$ \frac{\partial \ln(\pi(\delta))}{\partial \delta_q} = \sum_i \frac{\partial \ln(\pi(\delta))}{\partial \rho_i} \frac{\partial \rho_i}{\partial \delta_q}. \quad (F1) $$

Now we examine and expand the last term:

$$ \frac{\partial \rho_i}{\partial \delta_q} = \bar{n}_{\text{gal}} \left( b_1 + \frac{\partial \delta'_i}{\partial \delta_q} \right), \quad (F2) $$

where

$$ \delta'_i(\mathbf{k}, f_{nl}) = \Delta b(\mathbf{k}, f_{nl}) \delta_i(\mathbf{k}), \quad (F3) $$

with $\Delta b(\mathbf{k}, f_{nl})$ defined as in equation (8). Writing down each computation step, we get

$$ \begin{aligned} \delta'_{p}(\mathbf{x}) &= \sum_{b} \mathcal{F}_{p,b} \delta_{b}^{K}(\mathbf{k}), \\ \delta'_{j}(\mathbf{k}) &= \Delta b(\mathbf{k}, f_{nl}) \delta_{j}^{K}(\mathbf{k}), \\ \delta_{i}(\mathbf{k}) &= \sum_{a} \bar{\mathcal{F}}_{i,a} \delta_{a}^{K}(\mathbf{x}), \end{aligned} \quad (F4) $$

Taking the partial derivatives with respect to the inputs (in the forward case) yields the following.

$$ \begin{aligned} \frac{\partial \delta'_{p}(\mathbf{x})}{\partial \delta'_{j}(\mathbf{k})} &= \sum_{b} \mathcal{F}_{p,b} \delta_{b,j}^{K}, \\ \frac{\partial \delta'_{j}(\mathbf{k})}{\partial \delta_{i}(\mathbf{k})} &= \Delta b(\mathbf{k}, f_{nl}) \delta_{i,j}^{K}, \\ \frac{\partial \delta_{i}(\mathbf{k})}{\partial \delta_{q}(\mathbf{x})} &= \sum_{a} \bar{\mathcal{F}}_{i,a} \delta_{a,q}^{K}, \end{aligned} \quad (F5) $$

The final expression is as follows:

$$ \frac{\partial \ln(\pi(\delta))}{\partial \delta_q} = \sum_i \frac{\partial \ln(\pi(\delta))}{\partial \rho_i} \times \bar{n}_{\text{gal}} \left\{ b_1 + \left[ \sum_b \mathcal{F}_{p,b} \delta_{b,j}^K \right] [\Delta b(\mathbf{k}, f_{nl}) \delta_{i,j}^K] \times \left[ \sum_a \bar{\mathcal{F}}_{i,a} \delta_{a,q}^K \right] \right\} \quad (F6) $$

<figure><img src="image_18.png" /><figcaption>Figure F1. The plot shows the autocorrelation of sampled $f_{nl}$ values as a function of sample step of run #5. Note in the figure how the samples become uncorrelated early in the chain (by dropping from 1 to 0 quickly), and oscillates around 0. This means that the number of steps needed for the MCMC algorithm to produce statistically independent samples are relatively few, compared to the overall length of the full chain.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure F2. The plot shows the autocorrelation of sampled $f_{nl}$ values as a function of sample step of run #6.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure F3. The sky map displaying the observed and masked regions for the Stage IV runs of this study.</figcaption></figure>
Now, combining yields the full expression of the adjoint gradient. The log-posterior distribution is therefore

$$ \begin{aligned} \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_j^G} &= \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \frac{\partial \delta_i^{\text{NG}}}{\partial \delta_j^G} \\ &= \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \frac{\partial \delta_i^{\text{NG}}}{\partial \Phi_k} \frac{\partial \Phi_k}{\partial \phi_l} \frac{\partial \phi_l}{\partial \delta_j^G} \\ &= \sum_i \frac{\partial \ln(\pi(\delta^{\text{NG}}))}{\partial \delta_i^{\text{NG}}} \left[ \sum_b \bar{\mathcal{F}}_{i,b} \delta_{b,k}^K \right] [\delta_{k,l} + \delta_{k,l}^K f_{nl} \phi_k] \\ & \quad \times \left[ \sum_q \mathcal{F}_{1,q} \delta_{1,j}^K \right]. \end{aligned} \quad (E5) $$

## APPENDIX F: ADJOINT GRADIENT COMPUTATION OF BIAS MODEL

As in the previous section, we explicitly describe the adjoint gradient of the bias model. First, we drop the *m*-index for the matter density field $\delta_m$, such that $\delta_m \equiv \delta$, for the sake of brevity. Then, we perform the derivative of the log likelihood with respect to $\delta$:

$$ \frac{\partial \ln(\pi(\delta))}{\partial \delta_q} = \sum_i \frac{\partial \ln(\pi(\delta))}{\partial \rho_i} \frac{\partial \rho_i}{\partial \delta_q}. \quad (F1) $$

Now we examine and expand the last term:

$$ \frac{\partial \rho_i}{\partial \delta_q} = \bar{n}_{\text{gal}} \left( b_1 + \frac{\partial \delta'_i}{\partial \delta_q} \right), \quad (F2) $$

where

$$ \delta'_i(\mathbf{k}, f_{nl}) = \Delta b(\mathbf{k}, f_{nl}) \delta_i(\mathbf{k}), \quad (F3) $$

with $\Delta b(\mathbf{k}, f_{nl})$ defined as in equation (8). Writing down each computation step, we get

$$ \begin{aligned} \delta'_{p}(\mathbf{x}) &= \sum_{b} \mathcal{F}_{p,b} \delta_{b}^{K}(\mathbf{k}), \\ \delta'_{j}(\mathbf{k}) &= \Delta b(\mathbf{k}, f_{nl}) \delta_{j}^{K}(\mathbf{k}), \\ \delta_{i}(\mathbf{k}) &= \sum_{a} \bar{\mathcal{F}}_{i,a} \delta_{a}^{K}(\mathbf{x}), \end{aligned} \quad (F4) $$

Taking the partial derivatives with respect to the inputs (in the forward case) yields the following.

$$ \begin{aligned} \frac{\partial \delta'_{p}(\mathbf{x})}{\partial \delta'_{j}(\mathbf{k})} &= \sum_{b} \mathcal{F}_{p,b} \delta_{b,j}^{K}, \\ \frac{\partial \delta'_{j}(\mathbf{k})}{\partial \delta_{i}(\mathbf{k})} &= \Delta b(\mathbf{k}, f_{nl}) \delta_{i,j}^{K}, \\ \frac{\partial \delta_{i}(\mathbf{k})}{\partial \delta_{q}(\mathbf{x})} &= \sum_{a} \bar{\mathcal{F}}_{i,a} \delta_{a,q}^{K}, \end{aligned} \quad (F5) $$

The final expression is as follows:

$$ \frac{\partial \ln(\pi(\delta))}{\partial \delta_q} = \sum_i \frac{\partial \ln(\pi(\delta))}{\partial \rho_i} \times \bar{n}_{\text{gal}} \left\{ b_1 + \left[ \sum_b \mathcal{F}_{p,b} \delta_{b,j}^K \right] [\Delta b(\mathbf{k}, f_{nl}) \delta_{i,j}^K] \times \left[ \sum_a \bar{\mathcal{F}}_{i,a} \delta_{a,q}^K \right] \right\} \quad (F6) $$

<figure><img src="image_20.png" /><figcaption>**Figure F1.** The plot shows the autocorrelation of sampled $f_{nl}$ values as a function of sample step of run #5. Note in the figure how the samples become uncorrelated early in the chain (by dropping from 1 to 0 quickly), and oscillates around 0. This means that the number of steps needed for the MCMC algorithm to produce statistically independent samples are relatively few, compared to the overall length of the full chain.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure F2.** The plot shows the autocorrelation of sampled $f_{nl}$ values as a function of sample step of run #6.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure F3.** The sky map displaying the observed and masked regions for the Stage IV runs of this study.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure F4.** The radial selection functions for the *Stage IV* runs of this study.</figcaption></figure>

This paper has been typeset from a *TEX/LATEX* file prepared by the author.
<figure><img src="image_23.png" /><figcaption>**Figure F4.** The radial selection functions for the *Stage IV* runs of this study.</figcaption></figure>

This paper has been typeset from a *TEX/LATEX* file prepared by the author.