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

Adam Andrews1, Jens Jasche1,2, Guilhem Lavaux2, Fabian Schmidt3
1The Oskar Klein Centre, Department of Physics, Stockholm University, AlbaNova
University Centre, SE 106 91 Stockholm, Sweden
2CNRS & Sorbonne Université, Institut d’Astrophysique de Paris (IAP), UMR
7095, 98 bis bd Arago, 75014 Paris, France
3Max–Planck–Institut für Astrophysik, Karl–Schwarzschild–Straße 1, 85748
Garching, Germany E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### 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 utilise 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. Our method naturally and
fully self-consistently exploits the entirety of the large-scale structure,
e.g., higher-order statistics, peculiar velocity fields, and scale-dependent
galaxy bias, to extract information on the local non-Gaussianity parameter,
${f_{\mathrm{nl}}}$. We demonstrate the performance of our approach through
various tests with 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_{\mathrm{nl}}}$ by accurately handling survey geometries, noise, and
unknown galaxy biases. We demonstrate that our method can achieve constraints
of $\sigma_{{f_{\mathrm{nl}}}}\approx 8.78$ for SDSS-III/BOSS-like data, an
improvement of a factor $\sim 2.5$ over currently published constraints. Tests
with next-generation mock data show that significant further improvements are
feasible with sufficiently high resolution. Furthermore, the results
demonstrate that our method can consistently marginalise 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.

###### keywords:

statistics – large-scale structure of Universe – galaxies – inflation –
cosmological parameters

††pubyear: 2022††pagerange: Bayesian field-level inference of primordial non-
Gaussianity using next-generation galaxy surveys–16

##  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; Alvarez et al., 2014; Chen, 2010; Komatsu, 2010; Celoria & Matarrese,
2018; Senatore et al., 2010; Barnaby & Shandera, 2012; Chen & Wang, 2010). The
deviation from Gaussian initial conditions is commonly parameterised by the
nonlinearity parameter, $f_{\mathrm{nl}}^{\mathrm{local}}$, denoted as
${f_{\mathrm{nl}}}$ from now on. To date, the strongest constraint on
${f_{\mathrm{nl}}}$ is based on observations of the Cosmic Microwave
Background (CMB). CMB observations obtained by the Planck satellite yielded
${f_{\mathrm{nl}}}=-0.9\pm 5.1$ (Planck Collaboration, 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_{\mathrm{nl}}}|<1$, are expected to be beyond its reach (Camera et al.,
2013; Ballardini et al., 2019; Moradinezhad Dizgah & Keating, 2019;
Karagiannis et al., 2020; Meerburg et al., 2019).

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 yields the most information from the largest scales
(Dalal et al., 2008; Slosar et al., 2008). Previous works have been successful
in constraining PNG in available galaxy redshift surveys (Slosar et al., 2008;
Ross et al., 2012; Leistedt et al., 2014; Uhlemann et al., 2018; Karagiannis
et al., 2014; Mueller et al., 2019; Castorina et al., 2019; Mueller et al.,
2021; Cabass et al., 2022; D’Amico et al., 2022); the tightest large-scale
structure measurements comes from the SDSS-IV (DR 16, quasar sample) data,
reaching a constraining power of $|{f_{\mathrm{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 systematic effects, instead of by their
statistical power (Graham et al., 2018). Thus, a high degree of mitigation and
modelling of the survey effects, instrumentation noise, and astrophysical
contamination are required. If unaccounted for, these effects can bias the
results on the largest scales, (Huterer et al., 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 et al., 2019; Schmittfull et al., 2019; Porqueres et al., 2020; Schmidt
et al., 2020; Nguyen 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).

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 nonlinear
structure formation, and can incorporate any probe associated with PNG. As a
result, our method constrains ${f_{\mathrm{nl}}}$ beyond the current state-of-
the-art methods, and can place constraints as low as
$|{f_{\mathrm{nl}}}|<5.70$, for mock data emulating features of next-
generation galaxy redshift surveys (see Fig. 1), while marginalising out
systematic effects.

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

  1. 1.

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

  2. 2.

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)

  3. 3.

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

  4. 4.

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

For more details on the listed probes, the interested reader is referred to
reviews found in the literature (Biagetti, 2019). Another 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 model prediction. 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 both 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 nonlinear 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_{\mathrm{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. Results
are outlined 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.

Figure 1: Conditional posterior distribution for the local type primordial
non-Gaussianity parameter ${f_{\mathrm{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
marginalising out bias parameters.The yellow band indicates the region
$|{f_{\mathrm{nl}}}|\leq 1$, the target region of next-generation surveys. To
compare with the current best constraints, we plot a Gaussian centered on the
fiducial value of ${f_{\mathrm{nl}}}$, with a width corresponding to the
Planck18 constraints (Planck Collaboration, 2019b). 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 11). 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.222For a description on
(A) to (D), see section 2.2.1; for a description on (E) to (F), see section
2.2.2; for a description on (G), see section 2.3; for a description on (H),
see section 2.3.1.

##  2 Physical model and statistical method

In the present 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 primordial
non-Gaussianity (as parameterised 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. Thereafter, 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
scale-dependent bias (equation 9) 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 et al.,
2015; Lavaux & Jasche, 2016; Ramanah et al., 2019; Jasche & Lavaux, 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 realisations of the underlying dark
matter density field in the present Universe. Secondly, 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 the present 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_{\mathrm{nl}}}$ parameter, the bias parameters, and initial conditions
are updated. The iteration of these processes yields a valid Markov chain
permitting 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_{\mathrm{nl}}}$, conditioned on the galaxy number counts data
$N_{\textrm{g}}^{\textrm{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_{\mathrm{nl}}}$ is introduced in
equation (4). Formally, the posterior
$\pi(\epsilon,{f_{\mathrm{nl}}},\\{b_{i}\\}|N_{\textrm{g}}^{\textrm{o}})$ can
be written as:

$\displaystyle\pi(\epsilon,{f_{\mathrm{nl}}},\\{b_{i}\\}|N_{\textrm{g}}^{\textrm{o}})=\pi({f_{\mathrm{nl}}})\,\pi(\epsilon)\,\pi(\\{b_{i}^{\textrm{g}}\\})\frac{\pi(N_{\textrm{g}}^{\textrm{o}}|\epsilon,{f_{\mathrm{nl}}},\\{b_{i}\\})}{\pi(N_{\textrm{g}}^{\textrm{o}})}\,,$
(1)

where $\pi(\epsilon)$ is the Gaussian white-noise prior with zero mean and
unit standard variance,
$\pi(N_{\textrm{g}}^{\textrm{o}}|\epsilon,{f_{\mathrm{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.

# |  | Survey
---
mask
| Box size
---
$h^{-1}\;\text{Mpc}$
| Resolution
---
$h^{-1}\;\text{Mpc}$
$f_{\mathrm{nl}}^{\mathrm{true}}$
1 | SDSS-III/BOSS | 4000 | 31.25 | 5
2 | SDSS-III/BOSS | 4000 | 15.625 | 5
3 | Stage IV | 8000 | 62.5 | 5
4 | Stage IV | 8000 | 31.25 | 5
5 | Stage IV (Fixed bias) | 8000 | 31.25 | 5

The main run of this paper is run #4. 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 the bias values.

Table 1: An overview of the runs performed in this work.

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

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)

In this section, we will discuss the required modifications to the physical
model used in borg to account for primordial non-Gaussianity.

####  2.2.1 Generating non-Gaussian initial conditions

In this study, we are focused on inferring the local primordial non-
Gaussianity at the leading order, expressed through the ${f_{\mathrm{nl}}}$
parameter, using a forward modelling approach.

To describe the fluctuations of $\phi_{\textrm{G}}$, we convolve the Gaussian
white noise field ${\mathbf{\epsilon}}$ with a first transfer function
$T_{G}(k)$. That transfer function is designed such that $\phi_{\mathrm{G}}$
has a correlation structure given by a primordial power spectrum. In this
work, we use the classical parameterisation of the primordial power spectrum
in terms of $A_{s}$ and $n_{s}$. In the discrete case, the covariance of the
field $\phi_{\mathrm{G}}$ becomes as follows:

$\displaystyle\langle\hat{\phi}_{G,{\mathbf{a}}}\hat{\phi}^{*}_{G,{\mathbf{b}}}\rangle=V\delta^{K}_{{\mathbf{a}},-{\mathbf{b}}}A_{\textrm{s}}\frac{2\pi^{2}}{k_{\mathbf{a}}^{3}}\left(\frac{k_{\mathbf{a}}}{k_{\mathrm{pivot}}}\right)^{n_{s}-1}\equiv\delta^{K}_{{\mathbf{a}},-{\mathbf{b}}}T_{G}(k_{\mathbf{a}})^{2},$
(2)

with $V$ the volume of the box holding the forward modelled universe,
${\mathbf{a}}$ and ${\mathbf{b}}$ mesh indices, and
$\delta^{K}_{{\mathbf{a}},{\mathbf{b}}}$ the Kronecker delta. $A_{\textrm{s}}$
provides the amplitude of adiabatic, scalar, fluctuations of the potential in
the limit ${f_{\mathrm{nl}}}=0$. Thereafter, we generate the real space
primordial gravitational potential as:

$\displaystyle\centering\phi_{G}({\mathbf{x}})=\frac{1}{V}\sum_{{\mathbf{a}}}\exp(i{\mathbf{k}}_{{\mathbf{a}}}\cdot{\mathbf{x}})\hat{\phi}_{G,{\mathbf{a}}},\@add@centering$
(3)

with ${\mathbf{a}}$ a mesh index running over the 3D modes of the box,
$\hat{\phi}_{G,{\mathbf{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_{\textrm{NG}}({\mathbf{x}})$ is generated from a given Gaussian field
$\phi_{\textrm{G}}({\mathbf{x}})$ through the following Taylor expansion, as
highlighted by the box (C) in Figure 2:

$\displaystyle\Phi_{\textrm{NG}}({\mathbf{x}})=\phi_{\textrm{G}}({\mathbf{x}})+{f_{\mathrm{nl}}}\left(\phi_{\textrm{G}}({\mathbf{x}})^{2}-\langle\phi_{\textrm{G}}({\mathbf{x}})^{2}\rangle\right)\,.$
(4)

The parameters $A_{s},n_{s}$, and ${f_{\mathrm{nl}}}$ can be related to the
inflaton potential (Salopek & Bond, 1990; Verde et al., 2000; Komatsu &
Spergel, 2001).

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 et al., 2006). As
the evolution at high redshift is linear, the relation between matter density
$\hat{\delta}_{\mathrm{L}}({\mathbf{k}})$ and the primordial potential is, as
highlighted by the box (D) in Figure 2:

$\displaystyle\hat{\delta}_{\mathrm{L}}({\mathbf{k}},z)=T(k,z)\Phi_{\textrm{NG}}({\mathbf{k}})\,,$
(5)

where $z$ is taken at an arbitrary reference epoch. The transfer function
$T(k,z)$, which contains the linear growth factor, is provided either by CLASS
(Lesgourgues & Tram, 2014) or through analytic computation (Eisenstein & Hu,
1998, 1999; Takada et al., 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. To compute the present-day density
field, we populate the Lagrangian coordinates ${\mathbf{q}}$ with a set of
particles that are displaced according to:

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

The particle at the position
${\mathbf{x}}={\mathbf{q}}+{\mathbf{\Psi}}({\mathbf{q}},z)$ is then assigned,
with the Cloud-In-Cell kernel, to a regular grid, as highlighted by the box
(E) in Figure 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 Figure 2. This incorporates both the Kaiser and
Fingers of God effects, to the extent that they are captured by first-order
LPT, which may impact the inference of ${f_{\mathrm{nl}}}$ (Di Dio et al.,
2017; Tellarini et al., 2016; Bharadwaj et al., 2020; Tellarini et al., 2016;
Karagiannis et al., 2018). Handling RSDs allows us to include the peculiar
velocities when constraining ${f_{\mathrm{nl}}}$ with our field-based
approach.

Although LPT is sufficient to describe large-scale features, small-scale
features suffer from approximation issues. In this work, we focus on scales of
$\geq 16~{}h^{-1}\;\text{Mpc}$, which is 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 are
relying on LPT to describe structure formation, we can improve the structure
formation model to provide more realistic, nonlinear density field
realisations. Moreover, this allows us to combine the ${f_{\mathrm{nl}}}$
perturbation model with other cosmological models, e.g., a dark energy model
(as done in Ramanah et al., 2019). More information about this modular
approach, of which the ${f_{\mathrm{nl}}}$ perturbation model is a part of,
can be found in the literature (see Ramanah et al., 2019).

###  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 align 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 nonlinear analysis of galaxy surveys. For a review of
the galaxy bias problem, the interested reader is referred to the literature
(Desjacques et al., 2018; Schmidt et al., 2019).

At linear order, the relationship is the linear bias function:
$\rho_{\textrm{g}}=\langle
N_{\textrm{g}}^{\mathrm{o}}\rangle(1+b_{1}^{\mathrm{g}}\delta_{\mathrm{m}})$,
where $\rho_{\textrm{g}}$ is the resulting galaxy field, $\langle
N_{\textrm{g}}^{\mathrm{o}}\rangle$ is mean number of galaxies,
$b_{1}^{\mathrm{g}}$ is the linear bias and $\delta_{\mathrm{m}}$ is the final
density field, in the present-time universe. $\mathrm{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 et al., 2009). While 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 et al. (2009);
Biagetti et al. (2017) for measurements of $b_{\phi}$ on simulated dark matter
halos, and Barreira et al. (2020) for a measurement on simulated galaxies.

Here, as in all published constraints on ${f_{\mathrm{nl}}}$ so far, 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}^{\mathrm{g}}$.

In terms of the bias contribution:

$\rho_{\textrm{g}}=\langle
N_{\textrm{g}}^{\mathrm{o}}\rangle\left\\{1+\left[b_{1}^{\mathrm{g}}+\Delta
b(k,{f_{\mathrm{nl}}})\right]\delta_{\mathrm{m}}\right\\}\,,$

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

$\displaystyle\Delta
b(k,{f_{\mathrm{nl}}})=b_{\phi}\alpha(k,z)=b_{\phi}{f_{\mathrm{nl}}}\frac{3\Omega_{\mathrm{m}}H^{2}_{0}}{2c^{2}k^{2}T(k,z)}\,.$
(7)

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 et al., 2013):

$\displaystyle b_{\phi}=2\delta_{\mathrm{c}}(b_{1}^{\mathrm{g}}-1)\,,$ (8)

where $\delta_{\mathrm{c}}$ is the spherical collapse threshold,
$\delta_{\mathrm{c}}\approx 1.686$ (Gunn & Gott, 1972).

SDSS-III/BOSS | $b_{1}^{\mathrm{g}}$ |  $z_{\mathrm{min}}$/$z_{\mathrm{max}}$ |  $\bar{n}_{\mathrm{gal}}$ ($10^{-4}h^{3}\text{Mpc}^{-3}$) | Stage IV | $b_{1}^{\mathrm{g}}$ |  $z_{\mathrm{min}}$/$z_{\mathrm{max}}$ |  $\bar{n}_{\mathrm{gal}}$ ($10^{-4}h^{3}\text{Mpc}^{-3}$)
---|---|---|---|---|---|---|---
Cat. 1 | 3.22 | 0.4/0.75 | $3.01$ | Cat. 1 | 1.46 | 0.90/1.10 | 6.86
Cat. 2 | 2.48 | 0.4/0.75 | $2.99$ | Cat. 2 | 1.61 | 1.10/1.30 | 5.58
Cat. 3 | 1.78 | 0.2/0.5 | $6.31$ | Cat. 3 | 1.75 | 1.30/1.50 | 4.21
Cat. 4 | 0.87 | 0.2/0.5 | $10.4$ | Cat. 4 | 1.90 | 1.50/1.80 | 2.61
Table 2: Configuration of tracer catalogues and bias parameters. $b_{1}$ is
the linear bias parameter, $\bar{n}_{\mathrm{gal}}$ is the number density, and
$z$ is the redshift distance of the tracer bins. The noise
$\sigma_{\textrm{g}}$, introduced in equation (11), is fixed to $\sqrt{\langle
N_{\textrm{g}}^{\mathrm{o}}\rangle}$, the mean number of galaxies in catalogue
$g$.

Thus, the galaxy bias model is formulated with the linear bias
$b_{1}^{\mathrm{g}}$ and with the scale-dependent bias effect
$b_{\phi}(k,{f_{\mathrm{nl}}})$, which translates the dark matter density
contrast field $\delta_{\mathrm{m}}$ to the galaxy field $\rho_{\textrm{g}}$
as highlighted by the box (G) in Figure 2:

$\displaystyle\rho_{\textrm{g}}\left(\delta_{\mathrm{m}},\langle
N_{\textrm{g}}\rangle,b_{1}^{\mathrm{g}},{f_{\mathrm{nl}}}\right)=\langle
N_{\textrm{g}}\rangle\left\\{1+\left[b_{1}^{\mathrm{g}}+\Delta
b\left(k,{f_{\mathrm{nl}}}\right)\right]\delta_{\mathrm{m}}\right\\}.$ (9)

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 et al., 2015; Moradinezhad Dizgah
et al., 2021; Barreira et al., 2020; Barreira, 2020, 2021a, 2021b, 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_{\mathrm{nl}}}$ from galaxy clustering, it is
in general important also 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 to future work.

####  2.3.1 A likelihood model of galaxy data

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)

With a prediction of $\rho_{\textrm{g}}$, the next step is to relate this
field to $N_{\textrm{g,p}}$, the expected number of observed galaxies at a
given voxel $\mathrm{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 $\mathrm{g}$. In this way, by
applying $R_{\textrm{g,p}}$, the linear response operator of the survey, we
obtain as highlighted by the box (H) in Figure 2:

$\displaystyle
N_{\textrm{g,p}}=R_{\textrm{g,p}}[\rho_{\textrm{g,p}}(\delta_{\mathrm{m}},\langle
N_{\mathrm{g}}\rangle,b_{1}^{\mathrm{g}},{f_{\mathrm{nl}}})]\,,$ (10)

where $\langle N_{\mathrm{g}}\rangle$ is the mean number of galaxies in
catalogue $\mathrm{g}$. The purpose of $R_{\textrm{g,p}}$ is to account for
survey geometry and selection effects. We emphasise 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:

$\mathrm{ln}\left[\pi\left(N_{\textrm{g,p}}^{\mathrm{o}}\,\big{|}\,N_{\textrm{g,p}}\,,\sigma_{\textrm{g}}\right)\right]=-\textrm{P}\,\mathrm{ln}\left[\sqrt{2\pi}\,\sigma_{\textrm{g}}\right]-\\\
\frac{1}{2}\sum_{\textrm{p}=0}^{\textrm{P-1}}\left(\frac{N_{\textrm{g,p}}^{\mathrm{o}}-N_{\textrm{g,p}}}{\sigma_{\textrm{g}}}\right)^{2}\,,$
(11)

where $N_{\textrm{g,p}}^{\mathrm{o}}$ is the observed galaxy number count in
voxel $\mathrm{p}$, and $N_{\textrm{g,p}}$ is the model prediction of the
physical forward model, given $\epsilon$ and ${f_{\mathrm{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$:

$\displaystyle\sigma_{\textrm{g}}=\sqrt{\langle
N_{\textrm{g}}^{\mathrm{o}}\rangle}=\sqrt{\frac{\Sigma_{\textrm{p}}N_{\textrm{g,p}}}{\Sigma_{\textrm{p}}R_{\textrm{g,p}}}},$
(12)

where $R_{\textrm{p}}$ is the selection function for each voxel; summing over
$R_{\textrm{p}}$ yields the total number of observed voxels. This is chosen
such as 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_{\textrm{g,p}}$ is related to the
white-noise field $\epsilon$ via 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_{\mathrm{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_{\mathrm{nl}}}$. However, directly sampling 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_{\mathrm{i}}^{\textrm{g}}\\}$ (where $g$ is the galaxy catalogue number,
and $i$ is the index for the galaxy bias parameter), and ${f_{\mathrm{nl}}}$.
In this way, the conditional posteriors can be sampled from in a sequential
manner. This block sampling strategy is outlined in the following way:

$\displaystyle(1)$
$\displaystyle\>\epsilon^{\mathrm{n}+1}\curvearrowleft\pi\left(\epsilon|\\{b_{\mathrm{i}}^{\textrm{g}}\\}^{\mathrm{n}},f_{\mathrm{nl}}^{\mathrm{n}}\right)\,,$
$\displaystyle(2)$
$\displaystyle\>\\{b_{\mathrm{i}}^{\textrm{g}}\\}^{\mathrm{n}+1}\curvearrowleft\pi\left(\\{b_{\mathrm{i}}^{\textrm{g}}\\}|f_{\mathrm{nl}}^{\mathrm{n}},\epsilon^{\mathrm{n}+1}\right)\,,$
(13) $\displaystyle(3)$
$\displaystyle\>f_{\mathrm{nl}}^{\mathrm{n}+1}\curvearrowleft\pi\left(f_{\mathrm{nl}}^{\mathrm{n}}|\\{b_{\mathrm{i}}^{\textrm{g}}\\}^{\mathrm{n}+1},\epsilon^{\mathrm{n}+1}\right)\,,$

where $\curvearrowleft$ represents the process of drawing a value from the
probability distribution, and $n$ is the index of the sample. 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 Figure 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. Utilising standard Markov chain 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 utilising dynamical physical symmetries. For a more detailed
description of the sampling of the white-noise field $\epsilon$, see previous
works on the subject (Jasche & Wandelt, 2013; Ramanah et al., 2019).

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

Figure 5: The ${f_{\mathrm{nl}}}$ values sampled of run #4: the high
resolution Stage IV data set. The figure displays the conditional posterior
distribution of ${f_{\mathrm{nl}}}$, for the given mock data. Notice how the
fiducial ${f_{\mathrm{nl}}}$ value is situated within 1$\sigma$ of the sampled
mean.

##  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_{\mathrm{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; Amendola et al., 2018; Doré et al., 2014).

###  3.1 Creating mock galaxy surveys

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

  1. 1.

For the sake of demonstration, we aim to evaluate the physics forward model on
a cubic Cartesian box of side length $L=4000\,h^{-1}\;\text{Mpc}$ (SDSS-
III/BOSS-like) or $L=8000\,h^{-1}\;\text{Mpc}$ (Stage IV) and
$N_{\text{grid}}=128$ or $256$. This yields grid resolutions in the range of
$\Delta L\simeq 62.5h^{-1}\;\text{Mpc}$ to $\Delta L\simeq 15.6$.

  2. 2.

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_{\mathrm{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).

  3. 3.

Although the borg algorithm permits running full gravitational particle mesh
simulations, at this coarse resolution Lagrangian perturbation theory provides
a viable approximation to the structure formation problem (Moutarde et al.,
1991; Buchert et al., 1993; Bouchet et al., 1994; Scoccimarro, 2000;
Scoccimarro & Sheth, 2002). To reduce the noise of the simulated particle
distribution, we oversample the initial density by a factor $2$ yielding a
total number of $(2N)^{3}$ simulation particles. Particles are then evolved to
the present epoch, using LPT, and are assigned to a three-dimensional
Cartesian grid via the Cloud-In-Cell (CIC) 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.

  4. 4.

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

  5. 5.

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 are organised in Table 1 and Table
2, together with the 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.

To calculate the cosmological power spectrum and transfer functions, we assume
the following set of cosmological parameters ($\Omega_{\mathrm{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, 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_{\mathrm{nl}}^{\mathrm{fiducial}}=f_{\mathrm{nl}}^{\mathrm{true}}=5$. To
reiterate, the main aim of this study is to recover this ${f_{\mathrm{nl}}}$
value and to estimate the accuracy of our method.

###  3.2 Emulating SDSS-III/BOSS data

Figure 6: The plot shows the systematic drift of present-time power spectrum
amplitudes during burn-in phase of run #4. 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 $\sim 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.

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. Corresponding radial selection functions have
been estimated numerically by binning the equivalent distribution of galaxies
$N(d)$ in the actual CMASS and LOWZ survey (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 sub-catalogues, we will perform a
multi-tracer analysis. As has been previously demonstrated, using multi-tracer
populations reduces the impact of cosmic variance when inferring
${f_{\mathrm{nl}}}$ (Seljak, 2009). However, since the galaxies are placed in
separated sub-catalogues, we neglect any correlated noise between the galaxy
sub-catalogues. 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}_{\mathrm{gal}}$ for
each catalogue have been selected as table 1 in Alam et al. (2017). The
expected galaxy number density is related to the observed number count as:
$\bar{n}_{\mathrm{gal}}=\langle
N^{\mathrm{o}}_{\mathrm{g}}\rangle/V_{\mathrm{o}}$, where $V_{\mathrm{o}}$ is
the volume containing observed galaxies. The linear biases
$b_{1}^{\mathrm{g}}$ for each tracer catalogue has been selected as per figure
8 in Lavaux et al. (2019), with light-cone effects adjusted for.

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_{\mathrm{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
such as to cover a total area of $\sim$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 4
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_{\textrm{g}}(z)$ for the $g$th galaxy catalogue covering the
redshift range $z_{\text{min}}\leq z\leq z_{\text{max}}$ can be written as:

$\displaystyle S_{\textrm{g}}(z)=\left\\{\begin{array}[]{rl}1,&\text{if
}z\in[z_{\text{min}},z_{\text{max}}]\\\
0,&\text{otherwise}.\end{array}\right.$ (16)

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 Euclid-like specifications 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 Figure
16, and the sky maps in Figure 15; these can be found in the appendix (Section
B).

In addition, we also perform 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 marginalising out bias parameters affect the constraining power of
${f_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{nl}}}$, and various nuisance parameters, such as bias parameters
or unknown noise levels. To test the validity of the MCMC sampler, we
initialise 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 initialise 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. 6. It can be
seen that the sequentially 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 realisations. 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 unaccounted for, these
systematic effects will yield significant erroneous power at large scales,
concealing the physical signal of PNG from ${f_{\mathrm{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 emphasise that the analysis pipeline adopts a pure forward model approach,
as illustrated in Figure 2; thus, by correctly sampling the fiducial power
spectra, we demonstrate that our method is able to recover a key observable
accurately.

We also perform Gelman-Rubin convergence test of our chains, which they pass.
The outline of the test and the corresponding results can be found in Appendix
C.

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

Run # | Description | Sampled ${f_{\mathrm{nl}}}$ | $\sigma_{{f_{\mathrm{nl}}}}$
---|---|---|---
1 | SDSS-III/BOSS, low reso. | 10.11 | 11.9
2 | SDSS-III/BOSS, high reso. | 9.89 | 8.78
3 |  Stage IV low reso. | 0.59 | 7.09
4 |  Stage IV, high reso. | 7.64 | 5.70
5 |  Stage IV, fixed bias, high reso. | 3.17 | 4.43
Table 3: Summary of the measured ${f_{\mathrm{nl}}}$ values for the different
configurations. Our method is able to consistently measure the fiducial
${f_{\mathrm{nl}}}$ value, 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 #4. Figure 7: Summary of the
results of the runs performed. The vertical axis illustrates the estimated
${f_{\mathrm{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. Notice how the increase in
resolution improves the error on ${f_{\mathrm{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 marginalising them out. 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.

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

For each run listed in Table 1, we have computed the uncertainty on the
sampled ${f_{\mathrm{nl}}}$ value. These results are summarised in Fig. 7. As
can be seen, borg is able to constrain ${f_{\mathrm{nl}}}$ to
$\sigma_{{f_{\mathrm{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_{\mathrm{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 A.1.

###  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_{\mathrm{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_{\mathrm{nl}}}$-value can
be seen in Table 3. For a graphical visualisation, 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_{\mathrm{nl}}}$, leading to $\sigma_{f_{\mathrm{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 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, borg’s Markov chain provides a numerical representation of the
posterior distribution of three-dimensional initial density fields.
Specifically, the Markov chain produces an ensemble of about $\approx 50000$
initial density field realisations 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 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_{\mathrm{nl}}}$ parameter for a given realisation 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_{\mathrm{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_{\mathrm{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 visualised in terms of marginalised 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 marginalised posteriors of ${f_{\mathrm{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_{\mathrm{nl}}}$ while marginalising the unknown bias parameters
without any strong degeneracy within the physical model. This suggests that we
can jointly and independently sample ${f_{\mathrm{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
correlation lengths of model parameters. The resulting plots, together with a
longer discussion, can be found in the Appendix (Section A).

##  5 Summary and Conclusion

Figure 9: The two left panels each display two identical fields generated by
the same seed, varying only the ${f_{\mathrm{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_{\mathrm{nl}}}$ value. These maps illustrate the difference maps
as a function of ${f_{\mathrm{nl}}}$, which is a novel data product that our
method can infer.

This work presents a novel approach to investigating the physics of the
primordial universe, by 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 characterisation 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.

Figure 10: Cross-correlations and marginal posteriors for all hyperparameters
and ${f_{\mathrm{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_{\mathrm{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 marginalise these out to constrain the
${f_{\mathrm{nl}}}$ parameter.

To achieve this goal, we have expanded the Bayesian inference algorithm borg
to account for physical models of primordial non-Gaussianity in the generation
of initial conditions. Here, we illustrated our method in the case of PNG of
the local type, as parameterised by the cosmological parameter
${f_{\mathrm{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. Firstly, PNG affects the generation of the
cosmic initial conditions with the amplitude of non-Gaussian contributions
scaled by ${f_{\mathrm{nl}}}$. When evaluated with borg’s structure growth
model, non-Gaussian initial conditions will directly yield imprints on the
nonlinear matter and velocity distributions.

Secondly, PNG affects the formation of galaxies, which we account for by using
a model of scale-dependent galaxy bias in our data model. The updated physics
model therefore 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 remarked 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_{\mathrm{nl}}}$ and gravitational evolution
occur in two different steps in the forward model, borg can differentiate
between these two sources of nonlinearity in the data analysis procedure.

To test and illustrate the performance of the method, we applied it to two
artificial galaxy surveys emulating 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=4000h^{-1}\;\text{Mpc}$ and
$L=8000h^{-1}\;\text{Mpc}$. We also performed a resolution study by performing
inferences with $N_{\text{grid}}=128$ and $N_{\text{grid}}=256$, resulting in
grid resolutions ranging from $\Delta L=\sim 15.6h^{-1}\;\text{Mpc}$ to
$\Delta L=\sim 62.5h^{-1}\;\text{Mpc}$. We also configure a run with fixed
bias. This indicates that borg is able to handle unknown galaxy biases in a
robust framework. Given this coarse resolution, we chose to approximate the
gravitational structure formation using a Lagrangian perturbation theory
model. The investigation of higher resolution with borg’s particle mesh
simulation will be the subject of a coming publication.

The outcomes 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 of
importance for the inference of PNG in general and ${f_{\mathrm{nl}}}$ in
particular.

One of the major outcomes of this work is the demonstration that the
${f_{\mathrm{nl}}}$ parameter can be constrained with a field-level forward
modelling approach. We demonstrate that borg can constrain ${f_{\mathrm{nl}}}$
to $\sigma_{{f_{\mathrm{nl}}}}=8.78$ for a SDSS-III/BOSS-like survey, and
$\sigma_{{f_{\mathrm{nl}}}}=5.7$ for a Stage IV survey. A comparison with
other papers has been summarised in Table 4. We studied the inference at two
different resolutions, showing that the extracted information on
${f_{\mathrm{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 increasing the resolution further, which will be investigated in
upcoming studies. Another key point is the fact that our method automatically
marginalises 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 and corrects the correlations
between the parameter of interest, ${f_{\mathrm{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_{\mathrm{nl}}}$ to $\sigma_{{f_{\mathrm{nl}}}}=4.43$. For more tests and
their results, the interested reader is referred to the appendix (Section A).
In other words, by using all higher order statistics of galaxy clustering, we
can utilise more of the data than traditional data analysis tests (McQuinn,
2021).

A major advantage of our approach is its ability to characterise 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 either 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). One can, in
principle, fold in the kinetic Synyaev-Zeldovich observable to jointly
constrain ${f_{\mathrm{nl}}}$ in the same field-level forward modelling
framework, e.g., borg.

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_{\mathrm{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
marginalising out bias parameters, our method can infer ${f_{\mathrm{nl}}}$ at
the level of $5.70$ on next-generation galaxy surveys. Future work will be
focused on increasing the resolution of the reconstructions, thereby resolving
more information for constraining PNG. Our results demonstrates the promise of
Bayesian forward modelling to study the physics of the origin of the Universe
with next-generation galaxy surveys.

Author | $\sigma_{{f_{\mathrm{nl}}}}$ | Survey
---|---|---
This paper | $9$ | SDSS-III/BOSS (mock)
Mueller et al. (2021) | $21$ | SDSS-IV/eBOSS
D’Amico et al. (2022) | $29$ | SDSS-III/BOSS
Table 4: Comparison of our inferred ${f_{\mathrm{nl}}}$ constraints with other
papers measuring ${f_{\mathrm{nl}}}$ in galaxy redshift surveys. Note, while
the two last data sets are observed with the same instrument, they cover
different regions in the Universe and have detected different tracers.

## Acknowledgements

We are grateful to many enlightening discussions with Steffen Hagstotz,
Alexandre Barreira, Fabio Finelli, and Yashar Akrami. We would also like to
thank Eleni Tsaprazi, Nhat-Minh Nguyen and Natalia Porqueres for providing
feedback on the manuscript. AA acknowledges 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 utilised the HPC facility supported by
the Technical Division at the Department of Physics, Stockholm University.
This work is done within the Aquila Consortium333https://www.aquila-
consortium.org/. 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 et al. (2015) Alam S., et al., 2015, ApJ, 219, 12
  * Alam et al. (2017) Alam S., et al., 2017, MNRAS, 470, 2617–2652
  * Alvarez et al. (2014) Alvarez M., et al., 2014, arXiv:1412.4671
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Rev Relativ, 21, 2
  * Assassi et al. (2015) Assassi V., Baumann D., Schmidt F., 2015, J. Cosmology Astropart. Phys., 2015, 043
  * Baldauf et al. (2011) Baldauf T., Seljak U., Senatore L., 2011, JCAP, 2011, 006–006
  * Ballardini et al. (2019) Ballardini M., Matthewson W. L., Maartens R., 2019, MNRAS, 489, 1950
  * Barnaby & Shandera (2012) Barnaby N., Shandera S., 2012, JCAP, 2012, 034–034
  * Barreira (2020) Barreira A., 2020, JCAP, 2020, 031–031
  * Barreira (2021a) Barreira A., 2021a, arXiv:2112.03253
  * Barreira (2021b) Barreira A., 2021b, arXiv:2107.06887
  * Barreira et al. (2020) Barreira A., Cabass G., Schmidt F., Pillepich A., Nelson D., 2020, JCAP, 2020, 013–013
  * Baumann & Green (2021) Baumann D., Green D., 2021, arXiv:2112.14645
  * Bharadwaj et al. (2020) Bharadwaj S., Mazumdar A., Sarkar D., 2020, MNRAS, 493, 594
  * Biagetti (2019) Biagetti M., 2019, arXiv:1906.12244
  * Biagetti et al. (2017) Biagetti M., Lazeyras T., Baldauf T., Desjacques V., Schmidt F., 2017, MNRAS, 468, 3277–3288
  * Bond & Efstathiou (1984) Bond J. R., Efstathiou G., 1984, ApJ, 285, L45
  * Bouchet et al. (1994) Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1994, arXiv:astro-ph/9406013
  * Buchert et al. (1993) Buchert T., Melott A. L., Weiss A. G., 1993, arXiv:astro-ph/9309056
  * Cabass et al. (2022) Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022, arXiv:2201.07238
  * Camera et al. (2013) Camera S., Santos M. G., Ferreira P. G., Ferramacho L., 2013, PRL, 111
  * Castorina et al. (2019) Castorina E., et al., 2019, JCAP, 2019, 010
  * Catelan & Scherrer (1995) Catelan P., Scherrer R. J., 1995, ApJ, 445, 1
  * Celoria & Matarrese (2018) Celoria M., Matarrese S., 2018, arXiv:1812.08197,
  * Chen (2010) Chen X., 2010, A&A, 2010, 638979
  * Chen & Wang (2010) Chen X., Wang Y., 2010, Physical Review D, 81
  * Chodorowski & Bouchet (1996) Chodorowski M. J., Bouchet F. R., 1996, MNRAS, 279, 557
  * Crocce et al. (2006) Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369–381
  * D’Amico et al. (2022) D’Amico G., Lewandowski M., Senatore L., Zhang P., 2022, arXiv:2201.11518
  * Dalal et al. (2008) Dalal N., Doré O., Huterer D., Shirokov A., 2008, Phys. Rev. D, 77, 123514
  * Desjacques et al. (2009) Desjacques V., Seljak U., Iliev I. T., 2009, MNRAS, 396, 85–96
  * Desjacques et al. (2018) Desjacques V., Jeong D., Schmidt F., 2018, Physics Reports, 733, 1
  * Di Dio et al. (2017) Di Dio E., Perrier H., Durrer R., Marozzi G., Dizgah A. M., Noreña J., Riotto A., 2017, JCAP, 2017, 006
  * Doré et al. (2014) Doré O., et al., 2014, arXiv:1412.4872
  * Durrer et al. (2000) Durrer R., Juszkiewicz R., Kunz M., Uzan J.-P., 2000, Phys. Rev. D, 62, 021301
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, ApJ, 496, 605
  * Eisenstein & Hu (1999) Eisenstein D. J., Hu W., 1999, ApJ, 511, 5
  * Euclid Collaboration (2020) Euclid Collaboration 2020, A&A, 642, A191
  * Friedrich et al. (2019) Friedrich O., Uhlemann C., Villaescusa-Navarro F., Baldauf T., Manera M., Nishimichi T., 2019, arXiv:1912.06621
  * Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457
  * Gorski et al. (2005) Gorski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelman M., 2005, ApJ, 622, 759
  * Graham et al. (2018) Graham M. L., Connolly A. J., Ivezić v., Schmidt S. J., Jones R. L., Jurić M., Daniel S. F., Yoachim P., 2018, The Astronomical Journal, 155, 1
  * Gunn & Gott (1972) Gunn J. E., Gott J. Richard I., 1972, ApJ, 176, 1
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hastings (1970) Hastings W. K., 1970, Biometrika, 57, 97
  * Ho et al. (2015) Ho S., et al., 2015, JCAP, 2015, 040–040
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Huterer et al. (2013) Huterer D., Cunha C. E., Fang W., 2013, MNRAS, 432, 2945–2961
  * Jasche & Kitaura (2010) Jasche J., Kitaura F. S., 2010, MNRAS, 407, 29
  * Jasche & Lavaux (2015) Jasche J., Lavaux G., 2015, MNRAS, 447, 1204
  * Jasche & Lavaux (2017) Jasche J., Lavaux G., 2017, A&A, 606, A37
  * Jasche & Lavaux (2019) Jasche J., Lavaux G., 2019, A&A, 625, A64
  * Jasche & Wandelt (2012) Jasche J., Wandelt B. D., 2012, MNRAS, 425, 1042
  * Jasche & Wandelt (2013) Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894
  * Jasche et al. (2010) Jasche J., Kitaura F. S., Wandelt B. D., Enßlin T. A., 2010, MNRAS, 406, 60
  * Jasche et al. (2015) Jasche J., Leclercq F., Wandelt B. D., 2015, JCAP, 1, 036
  * Jeong & Schmidt (2015) Jeong D., Schmidt F., 2015, Class. Quant. Grav., 32, 044001
  * Kaiser (1984) Kaiser N., 1984, ApJ, 284, L9
  * Karagiannis et al. (2014) Karagiannis D., Shanks T., Ross N. P., 2014, MNRAS, 441, 486
  * Karagiannis et al. (2018) Karagiannis D., Lazanu A., Liguori M., Raccanelli A., Bartolo N., Verde L., 2018, MNRAS, 478, 1341–1376
  * Karagiannis et al. (2020) Karagiannis D., Slosar A., Liguori M., 2020, JCAP, 2020, 052–052
  * Kitaura & Enßlin (2008) Kitaura F. S., Enßlin T. A., 2008, MNRAS, 389, 497
  * Komatsu (2010) Komatsu E., 2010, CGQ, 27, 124010
  * Komatsu & Spergel (2001) Komatsu E., Spergel D. N., 2001, Phys. Rev. D, 63, 063002
  * Kostić et al. (2022) Kostić A., Jasche J., Ramanah D. K., Lavaux G., 2022, Astronomy & Astrophysics, 657, L17
  * LSST Science Collaboration (2009) LSST Science Collaboration 2009, arXiv:0912.0201
  * Lam et al. (2011) Lam T. Y., Nishimichi T., Yoshida N., 2011, MNRAS, 414, 289
  * Lavaux & Jasche (2016) Lavaux G., Jasche J., 2016, MNRAS, 455, 3169
  * Lavaux et al. (2019) Lavaux G., Jasche J., Leclercq F., 2019, arXiv:1909.06396,
  * Layzer (1956) Layzer D., 1956, Astron. J., 61, 383
  * Leclercq & Heavens (2021) Leclercq F., Heavens A., 2021, MNRAS, 506, L85–L90
  * Leistedt et al. (2014) Leistedt B., Peiris H. V., Roth N., 2014, PRL, 113, 221301
  * Lesgourgues & Tram (2014) Lesgourgues J., Tram T., 2014, JCAP, 2014, 032–032
  * Lewis (2019) Lewis A., 2019, arXiv:1910.13970
  * Ma et al. (2013) Ma Y.-Z., Taylor J. E., Scott D., 2013, MNRAS, 436, 2029
  * Maldacena (2003) Maldacena J., 2003, JHEP, 2003, 013
  * Matarrese & Verde (2008) Matarrese S., Verde L., 2008, ApJ, 677, L77
  * McQuinn (2021) McQuinn M., 2021, JCAP, 2021, 024
  * Meerburg et al. (2019) Meerburg P. D., et al., 2019, arXiv:1903.04409
  * Moradinezhad Dizgah & Keating (2019) Moradinezhad Dizgah A., Keating G. K., 2019, ApJ, 872, 126
  * Moradinezhad Dizgah et al. (2021) Moradinezhad Dizgah A., Biagetti M., Sefusatti E., Desjacques V., Noreña J., 2021, JCAP, 2021, 015
  * Moutarde et al. (1991) Moutarde F., Alimi J.-M., Bouchet F. R., Pellat R., Ramani A., 1991, ApJ, 382, 377
  * Mueller et al. (2019) Mueller E.-M., Percival W. J., Ruggeri R., 2019, MNRAS, 485, 4160
  * Mueller et al. (2021) Mueller E.-M., et al., 2021, arXiv:2106.13725
  * Neal (2000) Neal R. M., 2000, JCGS, 9, 249
  * Neal (2003) Neal R. M., 2003, Ann. Statist., 31, 705
  * Nguyen et al. (2020) Nguyen N.-M., Jasche J., Lavaux G., Schmidt F., 2020, JCAP, 2020, 011
  * Nguyen et al. (2021) Nguyen N.-M., Schmidt F., Lavaux G., Jasche J., 2021, J. Cosmology Astropart. Phys., 2021, 058
  * Peebles (1980) Peebles P. J. E., 1980, Princeton University Press, 1980
  * Pillepich et al. (2009) Pillepich A., Porciani C., Hahn O., 2009, MNRAS, 402, 191–206
  * Planck Collaboration (2019a) Planck Collaboration 2019a, arXiv:1905.05697
  * Planck Collaboration (2019b) Planck Collaboration 2019b, arXiv:1807.06211
  * Planck Collaboration (2020) Planck Collaboration 2020, A&A, 641, A6
  * Porqueres et al. (2019) Porqueres N., Kodi Ramanah D., Jasche J., Lavaux G., 2019, A&A, 624, A115
  * Porqueres et al. (2020) Porqueres N., Hahn O., Jasche J., Lavaux G., 2020, A&A, 642, A139
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  * De Putter, Roland & Doré (2017) De Putter, Roland Doré O., 2017, Phys. Rev. D, 95, 123513
  * Ramanah et al. (2019) Ramanah D. K., Lavaux G., Jasche J., Wandelt B. D., 2019, A&A, 621, A69
  * Ross et al. (2012) Ross A. J., et al., 2012, MNRAS, 428, 1116–1127
  * Ross et al. (2017) Ross A. J., et al., 2017, MNRAS, 464, 1168
  * Salopek & Bond (1990) Salopek D. S., Bond J. R., 1990, Phys. Rev. D, 42, 3936
  * Schmidt (2010) Schmidt F., 2010, Phys. Rev. D, 82, 063001
  * Schmidt (2021) Schmidt F., 2021, JCAP, 04, 033
  * Schmidt et al. (2013) Schmidt F., Jeong D., Desjacques V., 2013, Physical Review D, 88
  * Schmidt et al. (2019) Schmidt F., Elsner F., Jasche J., Nguyen N. M., Lavaux G., 2019, JCAP, 2019, 042
  * Schmidt et al. (2020) Schmidt F., Cabass G., Jasche J., Lavaux G., 2020, JCAP, 11, 008
  * Schmittfull et al. (2019) Schmittfull M., Simonović M., Assassi V., Zaldarriaga M., 2019, Physical Review D, 100, 043514
  * Scoccimarro (2000) Scoccimarro R., 2000, ApJ, 544, 597
  * Scoccimarro & Sheth (2002) Scoccimarro R., Sheth R. K., 2002, MNRAS, 329, 629
  * Seljak (2009) Seljak U., 2009, PRL, 102, 021302
  * Seljak et al. (2017) Seljak U., Aslanyan G., Feng Y., Modi C., 2017, JCAP, 12, 009
  * Senatore et al. (2010) Senatore L., Smith K. M., Zaldarriaga M., 2010, JCAP, 2010, 028–028
  * Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119–126
  * Slosar et al. (2008) Slosar A., Hirata C., Seljak U., Ho S., Padmanabhan N., 2008, JCAP, 2008, 031
  * Takada et al. (2006) Takada M., Komatsu E., Futamase T., 2006, Phys. Rev. D, 73, 083520
  * Tasinato et al. (2014) Tasinato G., Tellarini M., Ross A. J., Wands D., 2014, JCAP, 2014, 032–032
  * Tellarini et al. (2016) Tellarini M., Ross A. J., Tasinato G., Wands D., 2016, JCAP, 2016, 014
  * Uhlemann et al. (2018) Uhlemann C., Pajer E., Pichon C., Nishimichi T., Codis S., Bernardeau F., 2018, MNRAS, 474, 2853
  * Verde et al. (2000) Verde L., Wang L., Heavens A. F., Kamionkowski M., 2000, MNRAS, 313, 141
  * Wang et al. (2014) Wang H., Mo H. J., Yang X., Jing Y. P., Lin W. P., 2014, ApJ, 794, 94
  * Yokoyama et al. (2011) Yokoyama S., Sugiyama N., Zaroubi S., Silk J., 2011, MNRAS, 417, 1074

##  Appendix A Additional Results

(a) The marginalised posterior distribution of ${f_{\mathrm{nl}}}$ for run #1;
the low resolution SDSS-III/BOSS-like run.  (b) The marginalised posterior
distribution of ${f_{\mathrm{nl}}}$ for run #2; the high resolution SDSS-
III/BOSS-like run.  (c) The marginalised posterior distribution of
${f_{\mathrm{nl}}}$ for run #3; the low resolution Stage IV run.  (d) The
marginalised posterior distribution of ${f_{\mathrm{nl}}}$ for run #5; the
high resolution Stage IV run, with fixed bias sampling.

Figure 11: The inferred posterior distributions of ${f_{\mathrm{nl}}}$ for all
other runs. These distributions include handling and marginalising out all
other effects, e.g., survey geometry, instrumentation noise, and galaxy
biases, with the exception of run #5. Figure 12: The plot shows the
correlation matrix of run #2. The correlation matrix displays any internal
correlations between model parameters. Most importantly, any correlations
between ${f_{\mathrm{nl}}}$ and the galaxy bias parameters are undesired,
since they would indicate that the inferred ${f_{\mathrm{nl}}}$ may be biased
by the model parameters. The absence of these features indicate that we can
robustly infer the correct ${f_{\mathrm{nl}}}$ value with $1\sigma$
confidence.

###  A.1 ${f_{\mathrm{nl}}}$ posterior distributions of runs

In Fig. 11, we have illustrated the posterior distributions of
${f_{\mathrm{nl}}}$ for the secondary runs of this paper. While their
uncertainties have been provided in Table 3 and in 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 #5. The posterior distributions have their
centres within 1$\sigma$ of the fiducial values
(${f_{\mathrm{nl}}}^{\mathrm{fid.}}=5$), which indicates that our method can
consistently infer the correct ${f_{\mathrm{nl}}}$-distribution.

###  A.2 Correlation matrix

The correlation matrix of the high resolution Stage IV mock data run can be
found in Fig. 12. The figure displays the correlations and anti-correlations
between model parameters throughout the samples in the chain. As can be seen,
there is a very weak correlation between ${f_{\mathrm{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_{\mathrm{nl}}}$ and the hyper-parameters of the galaxy bias model. This
implies that the cosmological parameter 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 figure
of contour plots (Fig. 10), that the model has very few degeneracies between
the model parameters.

###  A.3 Correlation lengths

The correlation lengths of the high resolution Stage IV runs can be found in
Fig. 14 and 13. This plot displays the correlation of the ${f_{\mathrm{nl}}}$
parameter between the samples of the chains. This can be described as:

$\displaystyle
C({f_{\mathrm{nl}}})_{\mathrm{n}}=\frac{1}{N-n}\sum^{\mathrm{N}-\mathrm{n}}_{\mathrm{i}=0}\left(\frac{f_{\mathrm{nl}}^{\mathrm{i}}-\langle{f_{\mathrm{nl}}}\rangle}{\sigma_{{f_{\mathrm{nl}}}}}\frac{f_{\mathrm{nl}}^{\mathrm{i}+\mathrm{n}}-\langle{f_{\mathrm{nl}}}\rangle}{\sigma_{{f_{\mathrm{nl}}}}}\right)\,,$
(17)

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

##  Appendix B 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
Figure 16. The sky map can be seen in Figure 15.

##  Appendix C Gelman-Rubin test

In this section, we outline 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 $\mathrm{M}$
different subchains, each containing $\mathrm{n}$ samples, with $\mathrm{d}$
in-between (which are ignored in the analysis). Then, we compute the potential
scale reduction factor:

$\displaystyle\mathrm{PSRF}=\sqrt{\frac{V}{W}}\,,$ (18)

where:

  1. 1.

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

  2. 2.

$W=\frac{1}{M}\sum_{m=1}{M}\sigma_{m}^{2}$,

  3. 3.

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

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

For run #4, 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 D Adjoint Gradient Calculation of the ${f_{\mathrm{nl}}}$
Perturbation

In this section, we outline the adjoint gradient of the ${f_{\mathrm{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^{\mathrm{NG}}$, as the
input Gaussian density field $\delta^{\mathrm{G}}$ is altered. First, we drop
the $m$-index for the matter density field $\delta_{\mathrm{m}}$, such that
$\delta_{\mathrm{m}}\equiv\delta$, for the sake of brevity. Then, we start in
the general case:

$\displaystyle\frac{\partial\mathrm{ln}(\pi(\delta^{\mathrm{NG}}))}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}=\sum_{\mathrm{i}}\frac{\partial\mathrm{ln}(\pi(\delta^{\mathrm{NG}}))}{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}\frac{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}\,.$
(19)

Next, we evaluate
$\frac{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}$.
But in order to do so, we must first decompose
$\delta_{\mathrm{i}}^{\mathrm{NG}}$:

$\displaystyle\delta_{\mathrm{i}}^{\mathrm{NG}}=\sum_{\mathrm{a}}\mathcal{\bar{F}}_{\mathrm{i},\mathrm{a}}\Phi_{\mathrm{a}}\,,$
$\displaystyle\Phi_{\mathrm{a}}=\phi_{\mathrm{a}}+{f_{\mathrm{nl}}}\phi_{\mathrm{a}}^{2}\,,$
(20)
$\displaystyle\phi_{\mathrm{a}}=\sum_{\mathrm{b}}\mathcal{F}_{\mathrm{a},\mathrm{b}}\delta^{\textrm{G}}_{\mathrm{b}}\,,$

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_{\mathrm{i}}^{\mathrm{NG}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}$:

$\displaystyle\frac{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}=\frac{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}{\partial\Phi_{\mathrm{k}}}\frac{\partial\Phi_{\mathrm{k}}}{\partial\phi_{\mathrm{l}}}\frac{\partial\phi_{\mathrm{l}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}\,.$
(21)

The expressions for these partial derivatives are:

$\displaystyle\frac{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}{\partial\Phi_{\mathrm{k}}}=\sum_{\mathrm{b}}\mathcal{\bar{F}}_{\mathrm{i},\mathrm{b}}\delta^{\mathrm{K}}_{\mathrm{b},\mathrm{k}}\,,$
$\displaystyle\frac{\partial\Phi_{\mathrm{k}}}{\partial\phi_{\mathrm{l}}}=\delta^{\mathrm{K}}_{\mathrm{k},\mathrm{l}}+2\delta^{\mathrm{K}}_{\mathrm{k},\mathrm{l}}{f_{\mathrm{nl}}}\phi_{\mathrm{k}}\,,$
(22)
$\displaystyle\frac{\partial\phi_{\mathrm{l}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}=\sum_{\mathrm{q}}\mathcal{F}_{\mathrm{l},\mathrm{q}}\delta^{\mathrm{K}}_{\mathrm{l},\mathrm{j}}\,,$

where $\delta^{\mathrm{K}}$ is the Kronecker delta.

Now, combining yields the full expression of the adjoint gradient. The log-
posterior distribution is therefore:

$\displaystyle\frac{\partial\mathrm{ln}(\pi(\delta^{\mathrm{NG}}))}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}=\sum_{\mathrm{i}}\frac{\partial\mathrm{ln}(\pi(\delta^{\mathrm{NG}}))}{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}\frac{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}$
$\displaystyle=$
$\displaystyle\sum_{\mathrm{i}}\frac{\partial\mathrm{ln}(\pi(\delta^{\mathrm{NG}}))}{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}\frac{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}{\partial\Phi_{\mathrm{k}}}\frac{\partial\Phi_{\mathrm{k}}}{\partial\phi_{\mathrm{l}}}\frac{\partial\phi_{\mathrm{l}}}{\partial\delta_{\mathrm{j}}^{\mathrm{G}}}$
$\displaystyle=$
$\displaystyle\sum_{\mathrm{i}}\frac{\partial\mathrm{ln}(\pi(\delta^{\mathrm{NG}}))}{\partial\delta_{\mathrm{i}}^{\mathrm{NG}}}[\sum_{\mathrm{b}}\mathcal{\bar{F}}_{\mathrm{i},\mathrm{b}}\delta^{\mathrm{K}}_{\mathrm{b},\mathrm{k}}][\delta_{\mathrm{k},\mathrm{l}}+\delta^{\mathrm{K}}_{\mathrm{k},\mathrm{l}}{f_{\mathrm{nl}}}\phi_{\mathrm{k}}][\sum_{\mathrm{q}}\mathcal{F}_{\mathrm{l},\mathrm{q}}\delta^{\mathrm{K}}_{\mathrm{l},\mathrm{j}}]\,.$
(23)

##  Appendix E Adjoint Gradient Computation of Bias Model

Figure 13: The plot shows the autocorrelation of sampled ${f_{\mathrm{nl}}}$
values as a function of sample step of run #4. 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. Figure 14:
The plot shows the autocorrelation of sampled ${f_{\mathrm{nl}}}$ values as a
function of sample step of run #5.

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_{\mathrm{m}}$, such that $\delta_{\mathrm{m}}\equiv\delta$, for the
sake of brevity. Then, we perform the derivative of the log likelihood with
respect to $\delta$:

$\displaystyle\frac{\partial\mathrm{ln}(\pi(\delta))}{\partial\delta_{\mathrm{q}}}=\sum_{\mathrm{i}}\frac{\partial\mathrm{ln}(\pi(\delta))}{\partial\rho_{\mathrm{i}}}\frac{\partial\rho_{\mathrm{i}}}{\partial\delta_{\mathrm{q}}}\,.$
(24)

Now we examine and expand the last term:

$\displaystyle\frac{\partial\rho_{\mathrm{i}}}{\partial\delta_{\mathrm{q}}}$
$\displaystyle=$
$\displaystyle\bar{n}_{\mathrm{gal}}\left(b_{1}+\frac{\partial\delta^{\prime}_{\mathrm{i}}}{\partial\delta_{\mathrm{q}}}\right),$
(25)

where:

$\displaystyle\delta^{\prime}_{\mathrm{i}}({\mathbf{k}},{f_{\mathrm{nl}}})=\Delta
b({\mathbf{k}},{f_{\mathrm{nl}}})\delta_{\mathrm{i}}({\mathbf{k}}),$ (26)

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

$\displaystyle\delta^{\prime}_{\mathrm{p}}({\mathbf{x}})=\sum_{\mathrm{b}}\mathcal{F}_{\mathrm{p},\mathrm{b}}\delta_{\mathrm{b}}({\mathbf{k}})\,,$
$\displaystyle\delta^{\prime}_{\mathrm{j}}({\mathbf{k}})=\Delta
b({\mathbf{k}},{f_{\mathrm{nl}}})\delta_{\mathrm{j}}({\mathbf{k}})\,,$ (27)
$\displaystyle\delta_{\mathrm{i}}({\mathbf{k}})=\sum_{\mathrm{a}}\mathcal{\bar{F}}_{\mathrm{i},\mathrm{a}}\delta_{\mathrm{a}}({\mathbf{x}})\,,$

Taking the partial derivatives with respect to the inputs (in the forward
case), yields:

$\displaystyle\frac{\partial\delta^{\prime}_{\mathrm{p}}({\mathbf{x}})}{\partial\delta^{\prime}_{\mathrm{j}}({\mathbf{k}})}=\sum_{\mathrm{b}}\mathcal{F}_{\mathrm{p},\mathrm{b}}\delta^{\mathrm{K}}_{\mathrm{b},\mathrm{j}}\,,$
$\displaystyle\frac{\partial\delta^{\prime}_{\mathrm{j}}({\mathbf{k}})}{\partial\delta_{\mathrm{i}}({\mathbf{k}})}=\Delta
b({\mathbf{k}},{f_{\mathrm{nl}}})\delta^{\mathrm{K}}_{\mathrm{i},\mathrm{j}}\,,$
(28)
$\displaystyle\frac{\partial\delta_{\mathrm{i}}({\mathbf{k}})}{\partial\delta_{\mathrm{q}}({\mathbf{x}})}=\sum_{\mathrm{a}}\mathcal{\bar{F}}_{\mathrm{i},\mathrm{a}}\delta^{\mathrm{K}}_{\mathrm{a},\mathrm{q}}\,,$

The final expression is:

$\displaystyle\frac{\partial\mathrm{ln}(\pi(\delta))}{\partial\delta_{\mathrm{q}}}=\sum_{\mathrm{i}}\frac{\partial\mathrm{ln}(\pi(\delta))}{\partial\rho_{\mathrm{i}}}\times$
$\displaystyle\times\,\bar{n}_{\mathrm{gal}}\left\\{b_{1}+[\sum_{\mathrm{b}}\mathcal{F}_{\mathrm{p},\mathrm{b}}\delta^{\mathrm{K}}_{\mathrm{b},\mathrm{j}}][\Delta
b({\mathbf{k}},{f_{\mathrm{nl}}})\delta^{\mathrm{K}}_{\mathrm{i},\mathrm{j}}][\sum_{\mathrm{a}}\mathcal{\bar{F}}_{\mathrm{i},\mathrm{a}}\delta^{\mathrm{K}}_{\mathrm{a},\mathrm{q}}]\right\\}$
(29)

Figure 15: The sky map displaying the observed and masked regions for the
Stage IV runs of this study. Figure 16: The radial selection functions for the
Stage IV runs of this study.