11institutetext: Departamento de Física Teórica, Universidad Autónoma de
Madrid, 28049 Madrid, Spain 22institutetext: Centro de Investigación Avanzada
en Física Fundamental (CIAFF), Facultad de Ciencias, Universidad Autónoma de
Madrid, ES-28049 Madrid, Spain 33institutetext: Instituto de Física Teorica
UAM-CSIC, c/ Nicolás Cabrera 13-15, Cantoblanco, 28049 Madrid, Spain
44institutetext: Institut de física d’altes energies (IFAE) The Barcelona
Institute of Science and Technology campus UAB, 08193 Bellaterra Barcelona,
Spain 55institutetext: Institute for Astronomy, University of Edinburgh, Royal
Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, United Kingdom
66institutetext: Donostia International Physics Centre, Paseo Manuel de
Lardizabal 4, 20018 Donostia-San Sebastian, Spain. 77institutetext: Department
of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt
Lake City, UT 84112, USA 88institutetext: Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA
94305, USA. 99institutetext: Berkeley Center for Cosmological Physics,
University of California, Berkeley, CA 94720, USA 1010institutetext:
International Centre for Radio Astronomy Research, University of Western
Australia, 35 Stirling Highway, Crawley, Western Australia 6009, Australia

# PNG-UNITsims: the response of halo clustering to Primodial Non-Gaussianities
as a function of mass

Adrián Gutiérrez Adame<EMAIL_ADDRESS>Santiago Avila
<EMAIL_ADDRESS>Violeta Gonzalez-Perez 1122 Gustavo Yepes
<EMAIL_ADDRESS>Marcos Pellejero 5566 Chia-Hsun Chuang 7788 Yu Feng
99 Juan Garcia-Bellido 1133 Alexander Knebe 11221010 Mike S. Wang 55

(Received December 21, 2023)

Cosmological simulations play a key role in modern cosmology as they are
routinely used to analyse and interpret observational data. In this paper we
present the largest full N-body simulation to date with local primordial non-
Gaussianities (l-PNG), the PNG-UNIT. This simulation has been designed to
study the clustering of cosmological tracers in the presence of local PNG.
Tracks the evolution of $4096^{3}$ particles within a periodic box with
$L_{\rm box}=1\;h^{-1}\,{\rm Gpc}$, leading to a mass resolution of
$m_{p}=1.24\times 10^{9}\;h^{-1}\,M_{\odot}$. This is enough to resolve
galaxies targeted by stage-IV spectroscopic surveys. PNG-UNIT has Fixed
initial conditions whose phases are also Matched to the pre-existing UNIT
simulation. These two features in the simulations reduce our uncertainty
significantly so we use 100 FastPM mocks to estimate this reduction. The
amplitude of the non-Gaussianities used to set the initial conditions of this
new simulation is $f_{\rm NL}^{\rm local}=100$. In this first study, we use
mass selected dark matter haloes from the PNG-UNIT simulation to constrain the
local PNG parameters. PNG induce a scale dependent bias, parameterised through
$b_{\phi}$ or $p$, which might depend on the type of cosmological tracer.
Those cases when $p=1$ are referred to as the universality relation. We
measure $p$ as a function of the halo mass. Haloes with masses between
$1\times 10^{12}$ and $2\times 10^{13}\,h^{-1}M_{\odot}$ are well described by
the universality relation. For haloes with masses between $2\times 10^{10}$
and $1\times 10^{12}\,h^{-1}M_{\odot}$ we find that $p<1$ at $3\sigma$.
Combining all the mass bins, we find $p$ consistent with a value of $0.955\pm
0.013$, which is $3\sigma$ away from universality, as low mass haloes are more
numerous. We also study the effect of using priors on $p$ when constraining
$f_{\rm NL}$. Using the values we obtain for $b_{\phi}$ as priors, we forecast
that a DESI-like (stage-IV) survey will be able to constrain $f_{\rm NL}$
better than it the universality relation is assumed.

###### Key Words.:

Cosmology – large-scale structure of Universe – Methods: numerical

##  1 Introduction

Primordial Non-Gaussianities (PNG) provide a unique window into the dynamics
of the Universe during the inflationary era. One of the most studied models is
that of local Primordial non Gaussianities. The amplitude of the deviation
from the Gaussianity is typically parameterised by $f_{\rm NL}^{local}$
($f_{\rm NL}$ hereafter). Slow-roll single-field inflation predicts a nearly
vanishing $f_{\rm NL}$ (Maldacena, 2003; Creminelli & Zaldarriaga, 2004).
Therefore, detecting a significantly larger value ($\mathcal{O}(f_{\rm
NL})>1$) would point to a more complex scenario, such as multi-field inflation
(Lyth et al., 2003; Byrnes & Choi, 2010; Pajer et al., 2013). The current
tightest constraints on $f_{\rm NL}$ come from the bispectrum measured for the
anisotropies of the Cosmic Microwave Background, with $f_{\rm NL}=-0.9\pm
5.1\;(1\sigma)$ (Planck Collaboration et al., 2019). However, this measurement
is fundamentally limited by cosmic variance, and it is not expected to reach
the benchmark of $\sigma(f_{\rm NL})=1$, even with the addition of
polarisation data from the CMB (Baumann et al., 2009). Such an small error
will allow to either detect or rule out most models of multi-field inflation.

The most promising way to detect local-PNG is through galaxy surveys, using
the measurement of the so-called scale-dependent bias (Dalal et al., 2008;
Slosar et al., 2008; Matarrese & Verde, 2008). This method has been used
during the last few years, and has resulted in increasingly precise
constraints (Slosar et al., 2008; Ross et al., 2012; Giannantonio et al.,
2014; Ho et al., 2015; Leistedt et al., 2014; Castorina et al., 2019; Cabass
et al., 2022; Rezaie et al., 2023). The most precise measurement up to date
from spectroscopic surveys comes from SDSS-IV/eBOSS , which measured $f_{\rm
NL}=-12\pm 21\,(1\sigma)$ (Mueller et al., 2022). Current spectroscopic galaxy
surveys, such as DESI 111https://www.desi.lbl.gov/ and Euclid
222https://www.esa.int/Science_Exploration/Space_Science/Euclid_overview, are
expected to reach a level of precision of $\sigma(f_{\rm NL})\sim 5$ (Laureijs
et al., 2011; Giannantonio et al., 2012; Sartoris et al., 2016; DESI
Collaboration et al., 2016). The next generation of cosmological surveys are
expected to reach $\sigma(f_{\rm NL})<1$ by combining the information coming
from galaxy clustering with other probes, such as ${\rm HI}$ (LSST Science
Collaboration et al., 2009; Yamauchi et al., 2014; and David J. Bacon et al.,
2020).

Interpreting current and future cosmological surveys requires the use of large
simulations. They are routinely used for various purposes, from estimating
covariance matrices (Manera et al., 2012; Avila et al., 2018; Zhao et al.,
2021) to validating the analysis pipeline for the observational data (Avila et
al., 2020; Alam et al., 2021; Rossi et al., 2020) and studying the theoretical
models and observational aspects in a sandbox (e.g. Ross et al., 2017; Avila
et al., 2021; Chan et al., 2022; Riquelme et al., 2023). To date, there is no
adequate simulation to perform these studies in the presence of local-PNG for
Stage-IV spectroscopic surveys, given the limited mass resolution of the
existing ones (Grossi et al., 2009; Desjacques et al., 2009; Pillepich et al.,
2009; Hamaus et al., 2011; Scoccimarro et al., 2012; Wagner & Verde, 2012).

To study the scale-dependent bias with the precision required by current
galaxy surveys such as DESI ($\sigma(f_{\rm NL})\sim 5$), the required
simulation would need to have $>16\,000^{3}$ particles in a volume larger than
$(4\;{h^{-1}{\rm Gpc}})^{3}$. As this is extremely expensive to run, several
techniques have been developed in the recent years to increase the statistical
significance of simulations. This can be done either changing the way the
initial conditions are generated (for example the Fixed & Paired technique
proposed by Angulo & Pontzen, 2016) or by using a set of approximate methods
that are computationally less demanding/expensive (Chartier et al., 2021; Ding
et al., 2022; Kokron et al., 2022; DeRose et al., 2023).

It was not clear if the Fixed & Paired technique can be used in the presence
of Local-PNG as fixing the amplitudes of the Fourier modes breaks the
Gaussianity of the initial conditions. But in Avila & Adame (2023) we
demonstrated that this can be used also in the presence of local-PNG to have a
bigger effective volume for simulations.

In addition, in that work we introduced a technique to further suppress the
noise in the $f_{\rm NL}$ measurement. This is done by matching the phases
between two simulations with and without PNG, which removes a great part of
the noise in the measured statistics associated with the cosmic variance.

In this work we present the PNG-UNITsims, the largest simulation to date with
non-Gaussian initial conditions. This simulation evolves $4096^{3}$ dark
matter particles within a cosmological volume of $(1h^{-1}\,{\rm Gpc})^{3}$,
resulting in a mass resolution of $1.24\times 10^{9}\,h^{-1}M_{\odot}$. This
is enough to resolve the haloes hosting the emission-line galaxies (ELGs)
targeted by DESI and Euclid $(\sim 10^{11}\,h^{-1}M_{\odot})$ (Gonzalez-Perez
et al., 2018; Knebe et al., 2022; Yuan et al., 2023; Rocher et al., 2023).
Other DESI tracers (BGS, LRG, QSO) are expected to probe even larger masses,
so they are also well resolved by this simulation (Yu et al., 2023; Rocher et
al., 2023; Prada et al., 2023). The inital conditions of PNG-UNITsims are
generated with amplitudes fixed to their expectation values, and phases
matched to one of the original UNIT simulations with Gaussian initial
conditions.

As a first use of the PNG-UNITsims, the aim of this paper is to put
constraints on the non-Gaussian bias parameter $b_{\phi}$ for mass-selected
haloes. This parameter controls how the clustering of haloes, or any given
tracer of the matter density field (cosmological tracer hereafter), responds
to the presence of PNGs. However, $b_{\phi}$ is completely degenerate with the
amplitude of the non Gaussianities, measured by $f_{\rm NL}$ (Barreira,
2022a). It is crucial to have a good understanding of $b_{\phi}$ if we aim to
measure $f_{\rm NL}$ from data. Thus, understanding $b_{\phi}$ from
simulations will be essential.

From assuming a mass function for the haloes where the mass only depends on
the height of the density peak $\nu=\delta_{c}/\sigma(M)$, a theoretical
prediction for $b_{\phi}$ can be done as a function of the linear bias $b_{1}$
(Dalal et al., 2008). This is what is usually called universality relation.
Nevertheless, a number of works have shown that $b_{\phi}$ may not be
accurately described by that theoretical expectation (Slosar et al., 2008;
Grossi et al., 2009; Hamaus et al., 2011; Scoccimarro et al., 2012; Wagner &
Verde, 2012; Biagetti et al., 2017). Moreover, it has been shown that
$b_{\phi}$ may be different depending on the particular cosmological tracer
(Barreira et al., 2020; Barreira, 2022b).

In the literature, $b_{\phi}$ has been constrained using other techniques,
such as the separate universe technique (Biagetti et al., 2017; Barreira et
al., 2020; Barreira, 2022a; Lazeyras et al., 2023). Nevertheless, this
technique can be applied only to simulations, not observational data. And
although $b_{\phi}$ can be measured very precisely with this approach, these
simulations cannot be used to validate non-Gaussian analysis tools of, for
example, galaxy clustering, weak lensing statistics and cluster counts given
they assume Gaussian initial conditions.

In this paper, we divide the haloes into 12 pseudo-logarithmically spaced mass
bins and study their clustering by measuring their power spectrum. This
statistic is sensitive to the product $b_{\phi}f_{\rm NL}$ through the scale
dependent bias (Dalal et al., 2008; Slosar et al., 2008). Here we get an
estimation of the variance of the power spectrum for haloes in PNG-UNITsims
from a set of fast mocks generated with FastPM (Feng et al., 2016) with fixed
initial conditions. We further improve the precision measuring $b_{\phi}f_{\rm
NL}$ by matching the phases of the initial conditions. We leave for a future
work a study of how $b_{\phi}$ depends with the environment and the physics of
galaxy formation in a non-Gaussian simulation.

This paper is structured as follows. A brief review of the theory of galaxy
clustering in the presence of local-PNG is given in section 2, focusing on the
scale-dependent bias. In section 3, we present the PNG-UNITsims suit of
simulations developed for this study. In section 4, we describe the methods
employed for defining and relating the different mass bins and estimate the
power spectrum variance as well as how to apply the matching technique to
improve the constraints on the PNG response parameters. In section 5, we
present and discuss our main findings, including the constraints on the
$b_{\phi}$ parameter for mass selected haloes, together with the impact of
using those priors to measure $f_{\rm NL}$. In section 6 we study the
convergence of our analysis with different sets of simulations and mass
resolutions. The set of tests performed to ensure the robustness of our
results are described in section 7. In section 8, we discuss how the results
we have obtained can be used as a prior for a DESI-like survey and how it
would affect the constraints on $f_{\rm NL}$. And finally, we present the
conclusions in section 9.

##  2 Theory

Dark matter haloes trace the matter distribution, as these objects are assumed
to be formed from high density peaks. However, the gravitational collapse of
dark matter haloes is a non-linear and complex process. On the quasi-linear
scales, the complexities can be absorbed into a set of operators, whose
coefficients are called bias parameters (Bardeen et al., 1986; Desjacques et
al., 2018). In the context of perturbation theory, the density contrast of
haloes and also other cosmological tracers, such as galaxies, can be written
as a function of position and redshift:

$\delta_{h}(\@vec{x},z)=\sum_{O}b_{O}(z)O(\@vec{x},z)\;.$ (1)

The operators $O$ describe all the fields that may affect the matter density
field, and $b_{O}$ are the bias parameters associated with each of these
operators. If only local operators are considered, then the bias expansion can
only contain terms with exactly two spatial derivatives of the potential, that
is $O\propto\partial_{i}\partial_{j}\phi$ (see Desjacques et al., 2018, (for
an extensive review on the bias formalism)). All the terms proportional to the
matter density contrast $\delta_{m}$ (e.g. $\propto\delta_{m}{}^{n}$ ) are
allowed in the local bias expansion given it is related with the potential
through the Poisson equation $(\nabla^{2}\phi\propto\delta_{m})$. The leading
order of Equation 1 can be expressed as:

$\delta_{h}=b_{1}\delta_{m}\;,$ (2)

where $\delta_{h}$ is the density field of haloes, $\delta_{m}$ is the density
field of the matter distribution and $b_{1}$ is the linear bias. This
expression is only valid on the purely linear regime, which typically are
considered those below $k_{\rm max}\lesssim 0.1\,h\,{\rm Mpc}^{-1}$. By
including more terms in Equation 1, the clustering of the haloes can be
modelled even in the mildly non-linear regime. This is not necessary as the
effect of local-PNG is mostly encoded in the largest scales (see §2.2) and
hence, most part of our analysis will focus on the purely linear regime.

Using the linear bias, $b_{1}$, from Equation 2, the power spectrum of haloes
can be expressed in terms of the linear matter power spectrum, $P_{m,m}(k,z)$,
as:

$P_{h,h}(k,z)=\langle\delta_{h}\delta_{h}\rangle=b_{1}^{2}(z)P_{m,m}(k,z)\;.$
(3)

Below, we describe how PNG may arise from non-standard inflationary models
(§2.1) and how it induces a scale dependence in the linear galaxy/halo bias
(§2.2).

###  2.1 Primordial non-Gaussianities

The simplest inflationary models, such as single field slow-roll inflation,
predict a nearly Gaussian spectrum of perturbations (Maldacena, 2003;
Creminelli & Zaldarriaga, 2004). However, many models deviate from these
asumptions, which may lead to large primordial non-Gaussianities (PNG) (eg.
Byrnes & Choi, 2010; Ezquiaga et al., 2023). At first order, the amplitude of
these deviations from Gaussianity is typically measured in terms of the first-
order parameter $f_{\rm NL}$. In the case of Gaussian fields, all the
information is contained in the two-point function, and the higher-order
statistics, such as the bispectrum and the trispectrum, are determined by the
Wick theorem (Desjacques et al., 2018). However, in the presence of PNGs, this
relation is modified, and non-trivial contributions to the higher-order
statistics may arise.

At first order, this $f_{\rm NL}$ term induces a non-vanishing component in
the bispectrum. Specifically, by splitting the primordial potential into a
Gaussian and a non-Gaussian part as $\phi=\phi_{G}+\phi_{NG}$, it can be shown
that the leading order contribution to the bispectrum is given by:

$\langle\phi_{G},\phi_{G},\phi_{NG}\rangle=(2\pi)^{3}\delta_{D}(\@vec{k}_{1}+\@vec{k}_{2}+\@vec{k}_{3})\mathcal{B}_{\phi}(\@vec{k}_{1},\@vec{k}_{2},\@vec{k}_{3})$
(4)

where $\delta_{D}$ is the Dirac delta in 3D and
$\mathcal{B}_{\phi}(\@vec{k}_{1},\@vec{k}_{2},\@vec{k}_{3})$ is the primordial
bispectrum. Depending on the shape induced in the primordial bispectrum, 3
types of PNGs have been traditionally considered: local, equilateral and
orthogonal. Nevertheless, the landscape of possible non-Gaussianities is much
more rich and complex.

In this work, we focus only on the local type, which is one of the simplest
forms of PNGs. It can be described by the following parameterisation of the
primordial gravitational potential (Salopek & Bond, 1990; Komatsu & Spergel,
2001):

$\phi(\@vec{x})=\phi_{G}(\@vec{x})+f_{\rm
NL}\left(\phi_{G}(\@vec{x})^{2}-\langle\phi_{G}(\@vec{x})^{2}\rangle\right)\;,$
(5)

where $\phi_{G}$ is the Gaussian primordial gravitational potential. The total
primordial potential, $\phi$, is modified at quadratic order by a term
proportional to $f_{\rm NL}\phi_{G}^{2}$ and now is non-Gaussian.

Taking the expression for the local PNGs given in Equation 5, a non-vanishing
component to the three-point function is induced. In this case, the bispectrum
peaks in the squeezed limit ($k_{3}\ll k_{1},k_{2}$) and is given by:

$\mathcal{B}_{\phi}(\@vec{k}_{1},\@vec{k}_{2},\@vec{k}_{3})=2f_{\rm
NL}P_{\phi}(k_{1})P_{\phi}(k_{2})+(2\;{\rm cyc.\;perm}),$ (6)

where $P_{\phi}(k)$ is the power spectrum of the primordial potential $\phi$,
and $(2\;{\rm cyc.\;perm})$ are two cyclic permutations of the term
$P_{\phi}(k_{1})P_{\phi}(k_{2})$.

The perturbations in the potential lead to a density contrast in the matter
field through the Poisson equation. Then, the non-Gaussianities in the
primordial potential are induced also in the matter density field. These
fields can be related in the Fourier space through:

$\delta(k,z)=\alpha(k,z)\,\phi(k,z)\;,$ (7)

where $\alpha(k,z)$ is defined as:

$\alpha(k,z)=\frac{2D(z)}{3\Omega_{\rm
m,\;0}}\frac{c^{2}}{H_{0}^{2}}\frac{g(0)}{g(z_{\rm rad})}k^{2}T(k)\;.$ (8)

$\Omega_{\rm m,\;0}$ represents the matter density parameter today. $T(k)$ is
the transfer function normalised such that $T(k\rightarrow 0)=1$. The term
$D(z)$ is the growth factor, normalised to $D(z=0)=1$ following the Large
Scale-Structure (LSS) convention. For this reason, the factor $\frac{g(z_{\rm
rad})}{g(0)}$ appears in Equation 8, converting from the LSS convention to the
CMB one. Here, $g(z)$ is defined as $g(z)=(1+z)D(z)$. For the cosmology
described in Planck Collaboration et al. (2016) and we assume throughout the
paper, this factor turns out to be $\frac{g(z_{\rm rad})}{g(0)}=1.275$.

###  2.2 Scale-dependent bias

In the presence of local PNG, a new term proportional to $f_{\rm NL}\phi$ is
allowed in the bias expansion Equation 1. This results in the leading order
being:

$\centering\delta_{h}(k,z)=b_{1}(z)\delta_{m}(k,z)+b_{\phi}f_{\rm
NL}\phi=\left[b_{1}(z)+\frac{f_{\rm
NL}b_{\phi}(z)}{\alpha(k,z)}\right]\delta_{m}(k,z).\@add@centering$ (9)

Here, $b_{1}$ is the linear bias, $\alpha(k,z)$ is given in Equation 8.
$b_{\phi}$ is the bias parameter associated with the $\phi$ operator. This
term only appears in the presence of PNG, as it is associated to the
primordial bispectrum Equation 6 which vanishes in case of Gaussian initial
conditions. The term $\alpha$ introduces the dependence of $\delta_{h}$ with
the scale $k$, proportional to $1/k^{2}$.

Using this new bias expansion, the power spectrum for dark matter haloes, or
any other biased tracer of the total matter density, can be computed similarly
to Equation 3. The resulting halo power spectrum is given by:

$P_{h,h}(k,z)=\left[b_{1}(z)+\frac{f_{\rm
NL}b_{\phi}(z)}{\alpha(k,z)}\right]^{2}P_{m,m}(k,z).$ (10)

The relationship between $P_{h,h}(k,z)$ and $P_{m,m}(k,z)$ now depends on the
scale due to the term $\alpha(k,z)$, where the leading order effect is
$\propto 1/k^{2}$. This is what is commonly called ”scale-dependent bias”.

Figure 1: Theoretical power spectrum obtained from Equation 10. The blue line
shows the linear matter-matter power spectrum for the UNIT cosmology at $z=1$.
The $b_{1}$ parameter is the same for all cases and is chosen to be similar to
the one expected for the DESI emission-line galaxies (ELGs) at that redshift
(DESI Collaboration et al., 2016) . The orange line corresponds to the case of
a biased tracer with Gaussian initial conditions. The green line shows the
effect of a non-Gaussian contribution with $f_{\rm NL}=100$ and a PNG-response
parameter $b_{\phi}$ given by Equation 12 with $p=1$. The red line displays
how a different value of $f_{\rm NL}$ can lead to the same signal if the
$b_{\phi}$ parameter does not follow the universality relation. Red and green
lines are on top of each other, since the product $b_{\phi}f_{\rm NL}$ is the
same. Finally, the purple line shows the kind of signal we expect in case
$b_{\phi}f_{\rm NL}<0$.

However, different values of $f_{\rm NL}$ may produce the same signal
depending on the value of the bias parameter $b_{\phi}$ as they are perfectly
degenerate. This issue has already been discussed in the literature (Barreira,
2022a).

The parameter $b_{\phi}$, by definition, represents the response of the
abundance of haloes to the presence of perturbations in the primordial
potential $\phi$. As local-PNG induces a mixing between the long-wavelength
and short-wavelength perturbations (Equation 5), in the context of peak-
background split theory (Kaiser, 1984; Bardeen et al., 1986), it can be
demonstrated that $b_{\phi}$ can be thought also as the response of the number
density of haloes to a rescaling of the amplitude of the primordial power
spectrum (Desjacques et al., 2018; Barreira, 2022c):

$b_{\phi}=\frac{1}{\bar{n}_{h}}\frac{\partial\bar{n}_{h}}{\partial\log(f_{\rm
NL}\phi)}=\frac{1}{\bar{n}_{h}}\frac{\partial\bar{n}_{h}}{\partial\log\mathcal{A}_{s}}\;.$
(11)

Where $\bar{n}_{h}$ is the mean number of haloes and $\mathcal{A}_{s}$ is the
amplitude of the primordial scalar perturbations. Then, if an universal mass
function is assumed, then one can derive a theoretical prediction for the
value of $b_{\phi}$ (Dalal et al., 2008; Slosar et al., 2008; Desjacques et
al., 2018):

$b_{\phi}=2\delta_{c}\left(b_{1}-p\right)\;.$ (12)

Here, $\delta_{c}$ is the density threshold for the spherical collapse and it
is given by $\delta_{c}=1.686$ and $p$ may take different values depending on
the population of objects under consideration. By assuming an universal mass
function, then $p=1$ for all tracers of the matter distribution (Dalal et al.,
2008; Slosar et al., 2008). This is what is usually called universality
relation.

Recent studies have shown that this relation does not always accurately
describe the scale-dependent bias. For instance, Slosar et al. (2008) argued
that for haloes that have recently undergone a major merger process,
$b_{\phi}$ is better described by $2\delta_{c}(b_{1}-1.6)$, while Barreira et
al. (2020) showed that for galaxies selected by their stellar mass, the
relation with $b_{1}$ is better described by $2\delta_{c}(b_{1}-0.5)$.
Moreover, $b_{\phi}$ may depend not only on the mass of the halo (through
$b_{1}$) but also on other properties, such as the halo concentration
(Lazeyras et al., 2023; Fondi et al., 2023).

In Figure 1, we present the expected power spectrum of dark matter haloes at
linear order in perturbation theory for different values of the non-Gaussian
parameter, $f_{\rm NL}$. The value of $b_{1}$ has been chosen to be similar to
the one expected for the emission-line galaxies (ELGs), sample targeted by the
current galaxy surveys (Yu et al., 2023). For the Gaussian case, the halo
power spectrum is just a shifted version of the matter-matter power spectrum
in logarithmic scale, with the shift determined by the linear bias parameter
$b_{1}$, as given by Equation 3. At large scales, the term proportional to
$b_{\phi}f_{\rm NL}$, which appears only for the case with local-PNG, induces
and enhancement for $b_{\phi}f_{\rm NL}>0$ (red and green lines), and a
suppression if $b_{\phi}f_{\rm NL}<0$ (purple line) of the power spectrum
compared with the Gaussian case. When going to even larger scales, then the
term proportional to $(b_{\phi}f_{\rm NL})^{2}$ dominates and an enhancement
of the power spectrum is shown for both cases, $b_{\phi}f_{\rm NL}>0$ and
$b_{\phi}f_{\rm NL}<0$.

We remark that the scale-dependent bias is only sensitive to the product
$b_{\phi}f_{\rm NL}$. In this way, $f_{\rm NL}$ cannot be extracted from this
unless we have a prior knowledge of $b_{\phi}$ (Barreira, 2022a), as
illustrated in Figure 1, where the green and red lines show exactly the same
signal although the $f_{\rm NL}$ value is different.

##  3 Simulations

The UNITsims333http://www.unitsims.org/ simulation suite (Chuang et al., 2019)
aims to model the Large Scale Structure of the Universe to validate
theoretical models and analysis pipelines for extracting information about
cosmological parameters in the current galaxy surveys. For this purpose,
simulations need to have an effective volume equivalent or larger to these
surveys and a high mass resolution, enough to resolve the dark matter haloes
hosting the targeted galaxies/tracers.

$\Omega_{\rm m}$ | 0.3089
---|---
$\Omega_{\rm b}$ | 0.0486
$\Omega_{\Lambda}$ | 0.6911
$h\equiv H_{0}$/(100 km s-1 Mpc${}^{-1})$ | 0.6774
$\sigma_{8}$ | 0.8147
$n_{\rm s}$ | 0.9667
$f_{\rm NL}$ | 100
Table 1: The cosmological parameters of all UNITsim and PNG-UNITsims
simulation suite (Planck Collaboration et al., 2016).

Most state-of-the-art simulations, including the UNIT simulation (§3.1),
assume Gaussian cosmology and cannot be used directly to study the effect that
PNG has on the distribution of Large Scale Structure at later cosmic times. To
address this, we have developed the PNG-UNIT (§3.2), together with two lower-
resolution full N-body simulations (LR-UNIT §3.2.2) and a set of $100$ fast
N-body FastPM mocks (§3.3). In order to perform further robustness tests on
our analysis, we developed another set of smaller simulations, as we will
discuss in §3.4.

| UNIT | PNG-UNIT | LR-UNIT |  FastPM mocks
---|---|---|---|---
Code | L-Gadget 2 | L-Gadget 2 | L-Gadget 2 | FastPM
$N_{\rm particles}$ | $4096^{3}$ | $4096^{3}$ | $2048^{3}$ | $2048^{3}$
$m_{p}\;(h^{-1}\,M_{\odot})$ | $1.24\times 10^{9}$ | $1.24\times 10^{9}$ | $9.97\times 10^{9}$ | $9.97\times 10^{9}$
$f_{\rm NL}$ | 0 | 100 | 0 and 100 | 0 and 100
Halo Finder | Rockstar | Rockstar | Rockstar | FoF
Softening length $(h^{-1}{\rm kpc})$ | 6.1 | 6.1 | 12.2 | 244.14*
Table 2: Characteristics of the main simulations used throughout the present
paper. From top to bottom: the gravity solver, the number of dark matter
particles, the particle mass, he values of $f_{\rm NL}$ set in the initial
conditions, the halo finder and the softening length. For the FastPM mocks we
report the grid size used to compute the forces. All the simulations have a
box size of $L_{\rm box}=1\,{h^{-1}\,{\rm Gpc}}$. The cosmological parameters,
as well as the input power spectrum are exactly the same (with the exception
of $f_{\rm NL}$) between all the simulations. * FastPM is a particle-mesh code
so the force resolution is given by the size of the mesh.

All simulations presented in this paper share the same cosmology (Planck
Collaboration et al., 2016), which is summarised in Table 1. The
characteristics of the simulations are outlined in Table 2. Moreover, in this
work we focus our analysis in the snapshots at redshift $z=1.032$.

###  3.1 The UNITsims

The UNITsims are a set of 2 pairs of state-of-the-art, dark matter-only, full
N-body simulations with Fixed & Paired initial conditions (Angulo & Pontzen,
2016; Chuang et al., 2019). Each of them, tracks the evolution of $4096^{3}$
particles in a periodic box with a side length of $L=1\,h^{-1}\,{\rm Gpc}$
from redshift $z=99$ up to $z=0$. The mass resolution is $m_{p}=1.24\times
10^{9}\;h^{-1}\,M_{\odot}$, so the least massive haloes expected to host
galaxies targeted by current spectroscopic surveys, such as Euclid and DESI,
are well resolved in these simulations (Cochrane et al., 2017; Gonzalez-Perez
et al., 2018; Chuang et al., 2019; Knebe et al., 2022). In addition to the new
simulations we run, we are using also in this work one of these original
UNITsims, for reference the one labeled as FixAmpInvPhase_001.

The initial conditions were set using the second order Lagrangian Perturbation
Theory (2LPT), implementated in FastPM 444https://github.com/fastpm/fastpm
(Feng et al., 2016). The required initial power spectrum was computed using
the Boltzmann solver CAMB555https://camb.info/ (Lewis et al., 2000).

The dark matter particles were evolved using the N-body code L-Gadget2
(Springel, 2005). This code uses a Tree-PM algorithm to compute the forces
between particles, by splitting the Newtonian potential into a long-range and
a short-range components. We have employed the non-public version of the code
L-Gadget-2 as it is highly optimised in terms of memory usage. This code has
been widely used for generating many large-volume and high resolution
simulations such as the Millennium run (Springel et al., 2005) and the
Multidark suite (Klypin et al., 2016). The Plummer softening length was set to
$\epsilon=6\,{\rm kpc}/h$ and then the particles were evolved from $z=99$ up
to $z=0$, storing a total of $129$ snapshots logarithmically spaced in scale
factor.

Dark matter haloes were identified using the publicly available code
Rockstar666https://bitbucket.org/gfcstanford/rockstar/src/main/ (Behroozi et
al., 2012). This code takes advantage of the $6-D$ phase space information of
the particles to identify the dark matter haloes. These structures are
identified as having at least 20 particles per halo. For tracking their merger
histories and build the merger trees, we used ConsistentTrees. More details
regarding these simulations can be found in (Chuang et al., 2019).

###  3.2 The PNG-UNITsims

For this paper, we have developed the new PNG-UNITsims simulation suite, with
exactly the same features as the original UNIT simulations but with non-
Gaussian initial conditions. Specifically we set $f_{\rm NL}=100$. The suite
consist of a high-resolution simulation (§3.2.1), a lower resolution version
of it (§3.2.2) and a set of $100$ fast mocks generated with FastPM, half of
them with $f_{\rm NL}$ $=0$ and the other half with $f_{\rm NL}$ $=100$
(§3.3).

To generate the initial conditions of all of these simulations we have used
the FastPM code. In first place, the initial conditions are generated for
$\phi$ in the same way as if the simulation were Gaussian. Indeed, at this
point, the amplitudes of the Fourier modes can be fixed to their expectation
value as in Angulo & Pontzen (2016). Then, to induce the local PNGs, the
primordial potential is computed from the density field using Equation 7 and
then the transformation given in Equation 5 is applied. This procedure is
detailed and validated in Avila & Adame (2023).

####  3.2.1 The PNG-UNIT

The PNG-UNIT simulation is the main product of this work as it is the largest
full N-body simulations in terms of number of particles with primordial non-
Gaussianites to date. It is a dark matter-only, full N-body simulation that
includes local primordial non-Gaussianities. As the original UNIT simulation,
it tracks the evolution of $4096^{3}$ particles in a periodic box with a side
length of $L=1\,h^{-1}\,{\rm Gpc}$ from redshift $z=99$ up to $z=0$. The
details are shown in Table 2. Given the large volume and the high mass
resolution, this simulation is a unique laboratory to validate models of
galaxy clustering in the presence of local-PNG.

The value of $f_{\rm NL}$ used for generating the initial conditions is set to
$f_{\rm NL}=100$. This large value of $f_{\rm NL}$ is already ruled out by
Planck at $\sim 20\sigma$ (Planck Collaboration et al., 2019). However, we
wanted a strong signal from the scale dependent bias in order to understand
better this effect. With a more realistic $f_{\rm NL}$, given our limitation
on scales to $k_{\rm min}=k_{f}=0.006\,h\,{\rm Mpc}^{-1}$ we would not be able
detect the signal from the PNGs in the halo power spectrum.

In addition, the initial conditions of the PNG-UNITsims were generated using
the Fixed & Matched technique described in Avila & Adame (2023). The
amplitudes were set to the expectation value as in Angulo & Pontzen (2016).
This can be safely applied also in simulations with local primordial non-
Gaussianities as demonstrated in Avila & Adame (2023). Then, the phases were
chosen to be exactly the same as in the reference Gaussian UNIT simulation.

The N-body code and the halo finder were run using exactly the same
configuration as for the UNITsims.

####  3.2.2 The LR-UNIT simulations

In addition to the high-resolution simulation presented above, for this paper
we have developed two new, lower-resolution versions of the previous ones
(with $f_{\rm NL}=0$ and $f_{\rm NL}=100$). The main purpose of these LR-UNIT
simulations is to better understand the effects from mass resolution (see §6).

For this reason, we use exactly the same configuration as the high-resolution
simulations. The only difference is the resolution (see Table 2).

As we are using the same code and seed for generating the initial conditions,
all the long wavelength modes were set to the same exact values as in the
high-resolution simulations. Thus, the same region of the Universe is modelled
in these low resolution simulations, and the same structures can be identified
in both sets of simulations. We focus our analysis in the snapshot at redshift
$z=1$.

###  3.3 FastPM mocks

We have also generated a set of $100$ fast mocks using FastPM (Feng et al.,
2016). These mocks were split into two sets of $50$: one set had Gaussian
initial conditions ($f_{\rm NL}=0$) while the other set had non-Gaussian
initial conditions ($f_{\rm NL}=100$).

FastPM is a code that uses a hybrid method combining Perturbation theory with
a Particle-Mesh approach to quickly generate cosmological simulations at the
cost of poor resolution on the short wavelength modes. This code is
particularly useful for generating a large number of simulations in a short
amount of time which can be used, for example, for estimating covariance
matrices at a low computational cost. In this work, we use FastPM mocks to
obtain an estimation of the variance of the power spectrum of the haloes in
simulations with fixed initial conditions as well as for obtaining the
correlation coefficients required by the matching technique as we discuss in
subsection 4.3.

Each simulation box contains $2048^{3}$ dark matter particles within a
cosmological volume of $1\;(h^{-1}\,{\rm Gpc})^{3}$, the same volume as the
UNIT simulation. The mass resolution is $m_{p}=9.97\times
10^{10}\;h^{-1}\,{\rm M_{\odot}}$, which is $8$ times coarser than the high-
resolution simulations (as well as the LR-UNIT simulations) but still
sufficient for the purposes of this work. The force resolution in the
Particle-Mesh approach is given by the size of the grid used to compute the
potential and the forces. We use a regular grid of $4096^{3}$ cells giving a
spatial resolution of $244.14\;h^{-1}{\rm kpc}$.

Following the same procedure as for the full N-body simulations, the initial
conditions were set at $z=99$ using 2nd order Lagrangian Perturbation Theory
with fixed amplitudes. The particles were then evolved up to $z=1$ using $100$
timesteps, linearly spaced in time.

To generate the halo catalogues, we use the Friends-of-Friends (FoF)
algorithm, which is implemented by default in FastPM. The linking length is
set to the standard value of $l_{l}=0.2\,l_{p}$, where $l_{p}$ is the mean
interparticle distance.

The difference between the halo finder used for the FastPM and the PNG-
UNITsims catalogues leads to a different mass definition for dark matter
haloes. This issue is addressed in subsection 4.1.

###  3.4 Supporting Simulations

In addition to the PNG-UNITsims and the FastPM mocks presented above, we have
developed an additional set of simulations specially aimed to check the
robustness of our results. We want to check if the results obtained from
fitting $b_{\phi}$ in Equation 10 are consistent with measuring $b_{\phi}$
using the separate universe technique (see §6.1).Then, the set of simulations
is split into 2 groups:

  * •

High resolution (HR) separate universe

  * •

Low Resolution (LR) separate universe

Each group of separate universe simulations consists of 3 simulations with
$512^{3}$ dark matter particles within a periodic cubic box of $L_{\rm
box,\;HR}=67.5\,h^{-1}{\rm Mpc}$ and $L_{\rm box,\;LR}=250\,h^{-1}{\rm Mpc}$
for the high resolution and low resolution groups respectively. This leads to
a mass resolution which is 8 times better and 8 times worse respectively
compared to the baseline UNITsims. Again, the fiducial cosmology is the same
as the one of the UNIT simulation, summarised in Table 1. Following a similar
procedure as shown in Barreira (2022a), we produced three simulation with the
same initial conditions, but varying the $\mathcal{A}_{s}$ (or equivalently
$\sigma_{8}$) by a $5\%$. This allow us to measure $b_{\phi}$ by computing the
derivative Equation 19 numerically, as we discuss in subsection 6.1.

For all these simulations we evolved the particles using the L-Gadget 2 code,
just like the high resolution simulations and setting the Plummer softening
length to $1/40$ times the mean interparticle distance, which varies depending
on the mass resolution of the simulations.

###  3.5 Compared Halo mass functions

Figure 2 shows the cumulative mass function of the haloes from the PNG-
UNITsims and LR-UNITsims.

The mass function has been calculated from $20$ times the mass of the particle
in PNG-UNITsim $2.48\times 10^{10}h^{-1}\,M_{\odot})$ up to $1\times
10^{15}h^{-1}\,M_{\odot})$ in 50 bins, spaced on a logarithmic scale.

We find that the low-resolution simulations converge to the high-resolution
ones at $<5\%$ for $M\sim 100\,m_{\rm part}^{\rm LR}$, with $m_{\rm part}^{\rm
LR}$ being the particle mass of the lower resolution simulation. The shaded
area represents the threshold of $20$ particles per halo below which the mass
of the haloes is unreliable in the low resolution simulations.

The first thing we remark, is that the mass function of the FastPM mocks is
very different to ones from the full N-body simulations. The reason is that
the mass definition is different, since for the FastPM mocks haloes were
identified with a FoF algorithm while we used Rockstar for the haloes in the
other simulations. This highlights the need for establishing a correspondence
between the full N-body haloes and the FastPM haloes. We devote subsection 4.1
to this task.

Moreover, we find an increment of the high-mass haloes in the $f_{\rm NL}=100$
simulations with respect to the Gaussian ones. This was expected as a
consequence of the change of the PDF of $\delta$ due to the presence of local-
PNG (Matarrese et al., 2000; LoVerde et al., 2008).

The Gaussian UNIT simulation has an excess of low mass haloes compared with
the other simulations due to an issue with the initial conditions affecting
the older version of the FastPM code used to produce that simulation. This
problem produced an excess power in the power spectrum of the initial
conditions at intermediate scales ($k_{\rm nyquist}/4$ to $k_{\rm nyquist}/2$)
of the order of $10\%$. This results in an overproduction of haloes of the
same order of $\sim 10\%$ and explains the differences between the UNIT and
the PNG-UNIT simulations at the mass function level. However, we find that
this does not significantly affect the halo clustering measurements compared
to the low-resolution simulations, as we explain in section 6.

Figure 2: Upper: Halo mass functions for the simulations described in section
3. Lower: Ratio of halo mass functions relative to the Gaussian UNIT
simulation as a reference. Shaded areas represent the $M=20m_{\rm part}^{\rm
LR}$ mass threshold for the low-resolution simulation. The minimum mass shown
here represents the 20 particle limit for the higher resolution simulations.

##  4 Covariance matrix of the Power Spectrum in simulations with Fixed
initial conditions

In this work, we study the PNG parameters ($b_{\phi}$ and $f_{\rm NL}$) as a
function of halo mass. We measure them using the scale-dependent bias induced
by the PNG in the halo power spectrum (§2.2). The dark matter auto power
spectra of the simulations is obtained by interpolating the particles into a
regular grid with $N=512^{3}$ cells using a Cloud-in-Cell (CIC) scheme
(Hockney & Eastwood, 1981) with a k-binning width equivalent to the
fundamental frequency of the box $(k_{f}=2\pi/L_{\rm box}=0.00628\,h\,{\rm
Mpc}^{-1})$. We have used the Python library
Nbodykit777https://nbodykit.readthedocs.io/en/latest/ (Hand et al., 2018) for
this task. The halo power spectrum is computed in an equivalent way but using
the halo positions instead of the dark matter particles. In addition, we use a
large set of fast mocks to estimate the variance in the power spectrum from
the full N-body simulations. This approach allow us to take into account the
variance suppression due to the initial conditions (see §4.2). However, there
are fundamental differences between the fast mocks and the PNG-UNITsims. In
particular: how the haloes are defined, the mass and the resolution of the
simulations. We now discuss how to connect the haloes (and hence the measured
power spectrum) from the full N-body simulation with the corresponding ones
from the fast mocks.

###  4.1 Adapting FastPM mass bins to UNITsim

Throughout this work, we only consider main haloes, which are those that are
not substructure of a bigger halo. For simulations where we run Rockstar (see
Table 2) the halo mass is estimated from an spherical overdensity up to the
threshold of $\Delta=200\rho_{c}$, where $\rho_{c}$ is the critical density
(Behroozi et al., 2012). For the FastPM mocks, the halo mass is just given as
the sum of the particle mass of its FoF components. This leads to some
differences when choosing the mass bins. The edges of the mass bins are chosen
with pseudo-logarithmic spacing for UNIT, as summarised in Table 3.

The mass definitions in the full N-body simulations are very different from
those for the fast mocks. To relate haloes between these two approaches, we
aim to reproduce simultaneously the clustering and the abundances of the UNIT
haloes with the FastPM mocks. Then, for each bin of the UNITsims, we consider
a maximum mass for the haloes in the FastPM mocks and select the next $N$ most
massive haloes to recover the UNIT number density. Using a least squares
algorithm, we find the optimal maximum mass for FastPM that reproduces the
UNIT power spectrum for $k<0.1\,h\,{\rm Mpc}^{-1}$. The same mass cuts are
applied to the $f_{\rm NL}=100$ simulations.

Figure 3 shows the power spectrum of FastPM compared to that from the PNG-
UNITsims. The lines are the power spectrum computed for the PNG-UNITsims for
the mass bins from Table 3. The shaded areas are the mean and the standard
deviation of the power spectra obtained from the FastPM mocks using the
procedure described above.

For the three lightest bins, there are no haloes in the fast mocks since the
mass resolution is not high enough to resolve them. In the case of the two
most massive bins, although we are able to recover the clustering, we can not
recover simultaneously the abundances. This is associated to the differences
in the mass and bias functions due to the distinct halo finders used (see
Figure 2). In addition, the power spectra for the 3 most massive bins in the
PNG-UNITsims are very noisy since the number density is very low
$(n<10^{-5}\,h^{3}\,{\rm Mpc}^{-3})$.

Bin | Mass Range | $b_{\phi}f_{\rm NL}$ | $b_{\phi}f_{\rm NL}$ | $f_{\rm NL}\,(p=1)$ | $f_{\rm NL}\,(p=1)$ | $p$
---|---|---|---|---|---|---
| $(h^{-1}\,M_{\odot})$ | $\sigma_{S,\;{\rm UNIT}}$ | $\sigma_{F+M}$ | $\sigma_{S}$ | $\sigma_{F+M}$ | $\sigma_{F+M}$
0 | $(2,5]\times 10^{10}$ | $16\pm 55$ | $-37.9\pm 7.0$ | $30\pm 130$ | $94\pm 10$ | $0.993\pm 0.021$
1 | $(5,10]\times 10^{10}$ | $4\pm 54$ | $7.4\pm 9.1$ | $-90\pm 400$ | $-55\pm 45$ | $0.938\pm 0.027$
2 | $(1,2]\times 10^{11}$ | $2\pm 63$ | $48\pm 12$ | $490\pm 600$ | $409\pm 74$ | $0.892\pm 0.035$
3 | $(2,5]\times 10^{11}$ | $-8\pm 69$ | $67\pm 11$ | $140\pm 150$ | $138\pm 25$ | $0.946\pm 0.034$
4 | $(5,10]\times 10^{11}$ | $-15\pm 83$ | $128\pm 18$ | $157\pm 100$ | $138\pm 18$ | $0.895\pm 0.055$
5 | $(1,2]\times 10^{12}$ | $10\pm 100$ | $168\pm 29$ | $112\pm 68$ | $111\pm 15$ | $0.951\pm 0.086$
6 | $(2,5]\times 10^{12}$ | $-50\pm 110$ | $244\pm 26$ | $79\pm 52$ | $104\pm 11$ | $0.973\pm 0.077$
7 | $(5,10]\times 10^{12}$ | $-110\pm 140$ | $307\pm 43$ | $60\pm 46$ | $85\pm 12$ | $1.16\pm 0.13$
8 | $(1,2]\times 10^{13}$ | $-230\pm 160$ | $476\pm 90$ | $54\pm 39$ | $93\pm 13$ | $1.11\pm 0.27$
9 | $(2,5]\times 10^{13}$ | $50\pm 260$ | $630\pm 120$ | $89\pm 43$ | $88\pm 16$ | $1.26\pm 0.35$
10 | $(5,10]\times 10^{13}$ | $-2138\pm 940$ | $1020\pm 290$ | $-58\pm 57$ | $85\pm 26$ | $1.53\pm 0.93$
11 | $(1,2]\times 10^{14}$ | $-1796\pm 1800$ | $880\pm 600$ | $-70\pm 100$ | $48\pm 36$ | $3.8\pm 1.9$
12 | $(2,5]\times 10^{14}$ | $-548\pm 1800$ | $-4230\pm 740$ | $-15\pm 79$ | $-145\pm 27$ | $22.4\pm 2.4$
Table 3: PNG parameters derived for haloes in different mass bins. The first
column shows the ID of each bin. The second column presents the mass range for
the haloes as taken from the PNG-UNITsims (§4.1). In the third column, we show
the measurements we get for the product $b_{\phi}f_{\rm NL}$ using the
standard errors ($\sigma_{S}(P(k))$) for the Gaussian UNIT simulation. In the
fourth column we present the results of $b_{\phi}f_{\rm NL}$ after applying
the matching between the PNG-UNIT and the UNIT simulations (§5.2). The fifth
column provides the $f_{\rm NL}$ values obtained from the PNG-UNIT simulation
using standard errors and assuming the universality relation, $p=1$, in
Equation 12. In the sixth column, we present the results when applying the
reduced errors derived from the FastPM mocks and matching with the Gaussian
simulation (§5.3). The last column displays the values of $p$ obtained from
Equation 12 which are used to correct the value of $f_{\rm NL}$ in the fifth
column to the input value (§5.4). For the bins in gray, we find that we do not
recover $b_{\phi}f_{\rm NL}=0$ for the Gaussian UNIT simulation using the
standard errors (Equation 13) within $3\sigma$. This suggest that we may be
underestimating the $\sigma(P(k))$ in these bins, not only for this case but
also for the PNG-UNIT simulation.

Figure 3: Solid lines present the power spectrum of PNG-UNIT haloes within
different mass bins (see Table 3). Shaded regions show the mean $\pm 1\sigma$
from the FastPM mocks, generated using haloes in mass bins that reproduce the
clustering and the abundance of the N-body simulation as (see §4.1). For the
two lightest bins (1 and 2, shown in dark blues), there are no haloes in the
FastPM mocks since the mass resolution is not high enough to resolve them. The
top panel shows the $f_{\rm NL}=0$ case and the bottom one the$f_{\rm NL}=100$
one.

###  4.2 Variance estimation from fast mocks for Fixed simulations

We focus now on the computation of the covariance matrix for the power
spectrum. In the linear regime, the covariance matrix can be assumed to be
nearly diagonal. Therefore, the expected variance for a Gaussian density field
can be computed following the expression (Feldman et al., 1994):

$\sigma_{S}^{2}(k)=\frac{4\pi^{2}}{Vk^{2}\Delta
k}\left(P(k)+\frac{1}{\bar{n}}\right)^{2}\;,$ (13)

where $P(k)$ is the power spectrum, $\Delta k$ is the width of the power
spectrum bins, $\bar{n}$ is the mean number density of tracers, and $V$ is the
simulation volume. We refer to this quantity as $\sigma_{S}(P(k))$, the
standard error.

As the initial conditions of the simulations were fixed to their expectation
value (§3), a huge suppression is expected with respect to Equation 13,
particularly relevant for the largest scales (Angulo & Pontzen, 2016;
Villaescusa-Navarro et al., 2018; Avila & Adame, 2023; Maion et al., 2022).

Our approach is to estimate directly $\sigma(P(k))$ as the standard deviation
of the FastPM mocks and apply them for the parameter estimation using the full
N-body simulation. However, this quantity is still noisy due to the limited
number of mocks available to calculate it. Moreover, we still have the problem
of how to obtain this quantity for those bins that are beyond the resolution
limit of the FastPM mocks.

Therefore, we propose to fit the ratio $\sigma(P(k))/P(k)$ with a smooth
function in order to get a less noisy estimate of the errors. For that, we
measure $\sigma(P(k))/P(k)$ directly from the FastPM mocks (for the mass bins
3-10 from Table 3 whose power spectra are consistent with the ones coming from
the PNG-UNITsims) and we fit it to the following ansatz:

$\frac{\sigma((P(k))}{P(k)}=f(k;a,b,c)=a\left(1+\frac{1}{1+b\exp{(-ck})}\right)\;.$
(14)

In Figure 4, we show on solid lines the fits of the FastPM mocks to the
previous expression.

To extrapolate $\sigma(P(k))/P(k)$ beyond the mass resolution of the fast
mocks, we fitted the three parameters $a$, $b$ and $c$ from Equation 14 we
obtained as a function of the mass. For the parameter $a$, we use an
exponential function and for parameters $b$ and $c$ we use a polynomial of
order 2 as a function of $\log_{10}M_{\rm halo}$. The extrapolations are shown
in Figure 4 as dashed lines.

We will now adopt by the default the errors coming from the fit to Equation 14
for bins 3-9 and the extrapolation for bins 0-2 and 10-12. Nevertheless, we
will check in subsection 7.1 the impact of this assumption.

Figure 4: Variance over mean power spectrum for each mass bin (indicated by
color). Symbols show the data taken directly from the FastPM mocks, while
solid lines show results obtained by fitting these curves to Equation 14. The
extrapolation to each mass bin as described in subsection 4.2 is shown in
dashed lines.

###  4.3 Reducing the errors using Fixed & Matched simulations

When two simulations share the same set of random phases in their initial
conditions the cosmic variance is highly correlated between them, even in the
case they have different cosmologies. Here we discuss how we take advantage of
the correlated noise from two simulations with different $f_{\rm NL}$ for
constraining the PNG-bias parameters with a very high precision. This is what
we call matching the simulations.

The initial conditions for a N-body simulation are a realisation of a Gaussian
random field generated from the power spectrum of the primordial
perturbations. Then, two degrees of freedom remain, the amplitude of the
perturbations, whose expected value is determined by the power spectrum, and
the phase of those perturbations (which in principle are random uniformly
distributed between $0$ and $2\pi$). This stochasticity during the generation
of the ICs lead to some degree of noise when measuring quantities from the
simulation in the same way as we expect to be in the Universe. This is what is
usually called cosmic variance. The degree of noise can be suppressed, for
example, by fixing the amplitudes of the modes to their expectation value
(Angulo & Pontzen, 2016).

The idea presented in Avila & Adame (2023) is that if the initial conditions
of two simulations share the same stochastic part, that is the same random
phases and the same deviations in the amplitude from their corresponding
average (in case they do not have fixed ICs), then the cosmic variance is
expected to be highly correlated. Hence, one could cancel out a great part of
the noise affecting a measurement by matching that simulation with an existing
one with the same ICs but different cosmology.

In Avila & Adame (2023), we applied this for measuring $f_{\rm NL}$. For that
case, the only three ingredients are needed in the measurement: the value of
$f_{\rm NL}$ obtained from the Gaussian simulation, the same quantity from the
non-Gaussian simulation and the correlation coefficient $\rho$ between both
simulations. In that paper, we used a set of full N-body simulations with
matched (same stochastic part) ICs to measure $\rho$.

Here, we propose a slightly different approach by using just two simulations,
the UNIT and the PNG-UNIT, and a set of cheap fast mocks for estimating
$\rho$, instead of a large number of full N-body mocks, which for this mass
resolution would be computationally prohibitive. This is done by fitting
$b_{\phi}f_{\rm NL}$ for each bin and for each FastPM realisation, and then
the Pearson correlation coefficient is computed.

We find that the measurements are more correlated the lower the mass of the
haloes. This may be due to the fact that the lower the mean mass of the bins
we have defined, the higher the number of haloes that are included, and
therefore the less noisy the power spectrum measurement. However, studying
this evolution of the correlation with mass is beyond the scope of this paper.

Thus, we are left with the question of how to obtain the correlation
coefficients for those bins where we do not have the FastPM mocks. One
possibility is to apply a linear regression as a function of bin mass. This
would give us a correlation coefficient of $\rho\sim 0.9$ for the less massive
halo bins (ie. bins 0-2). While this does not affect the measurement of
$b_{\phi}$, it is critical for the error bars. However, we have no guarantee
that this relation of $\rho$ with the mass will saturate at some point. We
have therefore chosen to take a more conservative approach and assume for the
bins 0-2 the average of the $\rho$ that we have calculated from the FastPM
mocks (ie. bins 3-10).

In addition, we propose to apply this to constraint directly the product
$b_{\phi}f_{\rm NL}$ (instead of just $f_{\rm NL}$ as done in Avila & Adame
(2023)), which is the quantity we are sensitive to. Then, we construct this
new estimator for $b_{\phi}f_{\rm NL}$ which takes into account the matching
between both simulations:

$\Delta b_{\phi}f_{\rm NL}=[b_{\phi}f_{\rm NL}]_{100}-[b_{\phi}f_{\rm
NL}]_{0}\;,$ (15)

where $[b_{\phi}f_{\rm NL}]_{100}$ is the quantity measured in the non-
Gaussian simulation and $[b_{\phi}f_{\rm NL}]_{0}$ corresponding in the
Gaussian one. Given that the noise in $[b_{\phi}f_{\rm NL}]_{100}$ and
$[b_{\phi}f_{\rm NL}]_{0}$ is highly correlated, it is expected that $\Delta
b_{\phi}f_{\rm NL}$ will converge to the mean much faster than
$[b_{\phi}f_{\rm NL}]_{100}$ or $[b_{\phi}f_{\rm NL}]_{0}$ alone, given we are
cancelling some part of the noise. The error in $\Delta b_{\phi}f_{\rm NL}$ is
given by:

$\displaystyle\sigma^{2}(\Delta b_{\phi}f_{\rm NL})$
$\displaystyle=\sigma^{2}([b_{\phi}f_{\rm
NL}]_{100})+\sigma^{2}([b_{\phi}f_{\rm NL}]_{0})$ (16)
$\displaystyle-2\rho\sigma([b_{\phi}f_{\rm NL}]_{100})\sigma([b_{\phi}f_{\rm
NL}]_{0})\;,$

where $\sigma([b_{\phi}f_{\rm NL}]_{100})$ is the error of the quantity
measured in the non-Gaussian simulation and $\sigma([b_{\phi}f_{\rm NL}]_{0})$
corresponding to the Gaussian one. $\rho$ is the correlation coefficient
between both simulations, which is measured using the FastPM mocks.

While in principle this method can be applied for any two values of $f_{\rm
NL}$, the advantage of choosing one value to be $f_{\rm NL}=0$ is that
$\langle\Delta b_{\phi}f_{\rm NL}\rangle=\langle[b_{\phi}f_{\rm
NL}]_{100}\rangle$. From now on, when we apply this matching, we label the
errors as $\sigma_{S+M}$ if we use Equation 13 for the variance of the power
spectrum and $\sigma_{F+M}$ if we use the method described in subsection 4.2.

###  4.4 Variance estimation methods

In the previous sections we have discussed how we reduce the noise in the
measurements of $b_{\phi}f_{\rm NL}$ taking advantage of the fixed & matched
initial conditions of the PNG-UNITsims. To ensure that we are not introducing
any bias when using these techniques, we compare results with those derived
when no reduction of different reductions of the noise are assumed.

We define the following errors for the power spectrum, $\sigma(P(k))$ :

  * •

Standard errors, $\sigma_{S}$, come from Equation 13, using just one
simulation for measuring $b_{\phi}f_{\rm NL}$. Note that this error is not
correlated with the other $f_{\rm NL}$ simulation. This is expected to be the
largest estimation for the error.

  * •

Fixed errors, $\sigma_{F}$, are obtained from the best fit of $\sigma(P(k))$
measured from the FastPM mocks to Equation 14 for bins $3-10$ in halo mass.
For bins $0-2$ and $10-12$, we do not recover simultaneously the clustering
and the abundances with FastPM mocks, so we need to perform extrapolations for
those mass bins (subsection 4.3). We do not correlate $f_{\rm NL}=0$ with
$f_{\rm NL}=100$ to obtain these errors.

Errors that take into account the reduction in the noise due to the fixed &
matched initial conditions and used for propagating the errors using Equation
16:

  * •

Standard + Matched error, $\sigma_{S+M}$. We use the Standard error
$\sigma_{S}$ for the power spectrum but now we correlate the simulations with
$f_{\rm NL}=0$ with the ones with $f_{\rm NL}=100$ using the method described
in subsection 4.3.

  * •

Fixed + Matched errors, $\sigma_{F+M}$. In this case, we assume the
$\sigma_{F}$ errors for the power spectrum and we do correlate the simulations
with different $f_{\rm NL}$.

In subsection 7.2 we compare how these assumptions affect the constraints on
$b_{\phi}$ and $p$.

##  5 Constraining the PNG bias parameters for mass selected haloes

Given the input value of $f_{\rm NL}$ used to generate the initial conditions
for the PNG-UNITsims, we can measure directly $b_{\phi}$ from the scale-
dependent bias Equation 10. In this section, we discuss how we constraint
$b_{\phi}$ for mass selected haloes using the suppressed variance from the
matching technique (§4.3).

First, we check how well Equation 10 describes the clustering of the dark
matter haloes in the PNG-UNITsims under the standard assumption that
$b_{\phi}$ is given by the universality relation (i.e. p=1).

We measure the non-Gaussian signal $(b_{\phi}f_{\rm NL})$ in our simulations
(§5.2) and by assuming the universality relation we report our results for
$f_{\rm NL}$ (§5.3). Finally, using the input values of $f_{\rm NL}$, we
measure directly $p$ (§5.4).

In Table 3 we summarise the main results from this section.

###  5.1 Parameter estimation and analysis choices

Here we describe the pipeline applied for the parameter estimation, which is
common to all our analyses. We study the impact of varying these assumptions
in section 7.

To obtain the PNG parameters we run a MCMC using the Python package
emcee888https://emcee.readthedocs.io/en/stable/ (Foreman-Mackey et al., 2013).
All the cosmological parameters are fixed except $b_{1}$ and $b_{\phi}f_{\rm
NL}$ for which we set flat priors over a wide range ($[0,10]$ for $b_{1}$ and
$[-10000,10000]$ for $b_{\phi}f_{\rm NL}$). We assume a Gaussian likelihood
with $\chi^{2}$ given by:

$\chi^{2}=\sum^{N_{k}}_{i=0}\frac{\left(P_{\rm haloes}(k_{i})-P_{\rm
model}(k_{i};b_{1};b_{\phi}f_{\rm NL})\right)^{2}}{\sigma^{2}(P(k_{i}))}\;.$
(17)

Here, $P_{\rm model}(k)$ is computed using Equation 10 by inputting
$P_{m,m}(k)$ which computed directly from the simulation. In the baseline
analysis, for $\sigma(P(k))$, we use the fixed errors $\sigma_{F}$ and
fixed+matched errors $\sigma_{F+M}$ defined in subsection 4.4. Moreover, in
the following we will use solid circles when we are using the best fit to
Equation 14 and the mocks, and we will use empty circles when we use the
extrapolation (ie. for bins 0-2).

Then, $\sigma(b_{\phi}f_{\rm NL})$ is obtained from the posterior after
marginalising over $b_{1}$. The model described by Equation 10 is only
accurate in the purely linear regime.

Due to the limitation of our model to the linear regime, we need to drop out
the modes associated with those higher order terms. This is done by applying a
scale cut in the power spectrum. We choose to set $k_{\rm max}$ to
$0.1\;h\,{\rm Mpc}^{-1}$, as the next-to-leading order terms start to become
important beyond this scale. We discuss this issue and study its impact on the
constraints on $b_{\phi}$ and $f_{\rm NL}$ and we find that variations of the
$50\%$ in $k_{\rm max}$ with respect to our fiducial choice have a negligible
effect (see §7).

###  5.2 Constraints on $b_{\phi}f_{\rm NL}$

Following the pipeline described above, we fix all the cosmological parameters
for our analysis except $b_{1}$ and $b_{\phi}f_{\rm NL}$ in Equation 10. Then
we search for the best fit parameters.

Figure 5: Constraints on $b_{\phi}f_{\rm NL}$ as a function of $b_{1}$. Blue
and orange symbols show results for the Gaussian and the non-Gaussian
simulations, respectively, assuming that the $P(k)$ variance is computed from
the FastPM mocks (see §4.2). Gray lines indicate the expected value by
assuming the universality relation, the dotted line for $f_{\rm NL}=0$ and the
dashed line for $f_{\rm NL}=100$. Results from the matching technique (see
§4.3) are shown in green. In solid circles we use the best fit to Equation 14
and in rings we use the extrapolation. Bins 0-9 are shown.

The results are shown in the third column of Table 3 and in Figure 5, where
the blue and orange symbols represent the values for the Gaussian and non-
Gaussian simulations, respectively, and the dashed and dotted lines indicate
the expected values for each case. The error bars show the $1\sigma$
confidence interval obtained from our fits for $b_{\phi}f_{\rm NL}$. For the
blue and orange symbols, we have used the variances obtained from the fits to
the FastPM mocks (see §4.2), while for the green symbols we use the matching
technique (see §4.3). For the Gaussian simulation, we expect to get exactly 0
while for the non-Gaussian one, the reference we assume is the prediction of
the universality relation (ie. p=1). These expectations are shown by the grey
lines.

For haloes with $b_{1}>1$, we clearly measure $b_{\phi}f_{\rm NL}$. When
$b_{1}\sim 1$, for the most popular assumption of $p=1$ in Equation 12,
$b_{\phi}$ is nearly $0$ and then the non-Gaussian signal is masked.

In the same figure, we also observe fluctuations for the $f_{\rm NL}=0$
simulation with respect to zero, which is expected due to the intrinsic cosmic
variance of the simulation.

Since both simulations originate from the same random phases, we expect the
measurements of $b_{\phi}f_{\rm NL}$ to be highly correlated (§4.3).
Therefore, if a perturbation is found in the Gaussian ($f_{\rm NL}=0$)
simulation that lowers the measured $b_{\phi}f_{\rm NL}$, an equivalent
perturbation is expected to be in the non-Gaussian ($f_{\rm NL}=100$)
simulation, as we actually see in Figure 5.

We do this to compute the matched data points from both simulations using
Equation 15 and obtain the green points in Figure 5. These data points are
expected to have less noise than the individual cases of the simulations and
we find they are closer to the prediction from the universality relation
($p=1$). The error bars increase with the halo mass. This is mainly due to the
increase in the variance of the power spectra resulting from the decrease in
the number densities.

Although not shown in the figure, for the bin 10 we do not recover
$b_{\phi}f_{\rm NL}=0$ at more than $2\sigma$ even assuming the standard
errors, $\sigma_{S}(P(k))$, for the UNIT simulation (see column 2 in Table 3).
Given this simulation has Gaussian initial conditions, the expected value for
$f_{\rm NL}$ is exactly 0. It is possible that a statistical fluctuation may
appear that favours a different measure of $b_{\phi}f_{\rm NL}$. But in such a
case, the error bars should be wide enough to contain $b_{\phi}f_{\rm NL}=0$.
This suggest that we may be underestimating the errors in Equation 13,
probably because of the term associated with the Poissonian shotnoise. For the
sake of accuracy, we are not considering these bins in the following analyses.

For $M_{\rm halo}<1\times 10^{12}h^{-1}M_{\odot}$, we find that the
universality relation is below what we measure. By fitting all the mass bins
using the relation of $b_{\phi}$ as a function of $b_{1}$ and $p$, Equation
12, and fixing the input value of $f_{\rm NL}$, we find that they are
described by $p=0.955\pm 0.013$. If we not take into account those bins in
which we used the extrapolation for $\rho$ and $\sigma_{F}(P(k))$, just using
these quantities measured directly from the FastPM mocks, we get $p=0.950\pm
0.025$.

###  5.3 Measuring $f_{\rm NL}$ from $b_{\phi}f_{\rm NL}$ assuming the
universality relation

We cannot constrain both $b_{\phi}$ and $f_{\rm NL}$ simultaneously since they
are perfectly degenerate.

To measure $f_{\rm NL}$, we assume the universality relation, with $p=1$. In
this way, we extract $f_{\rm NL}$ from our simulation results. If this
relation does not describe accurately the $b_{\phi}$ parameter, then we expect
to recover biased values of $f_{\rm NL}$. Hence, we are testing the validity
of the universality relation.

Figure 6: Value of $f_{\rm NL}$ for each mass bin, assuming $b_{\phi}$ is
described by Equation 12 with $p=1$. Blue and orange symbols show results for
the Gaussian and non-Gaussian simulations, respectively, while green symbols
show the matching between them. Dashed lines show the input $f_{\rm NL}$
values for each simulation. The shaded area represents the mass bins that have
$b_{1}-1\leq 0.1$ and hence $b_{\phi}\sim 0$. In solid circles we use the best
fit to Equation 14 and in rings we use the extrapolation.

In Figure 6, we compare the values of $f_{\rm NL}$ obtained in this way with
the ones set in the initial conditions. For the three lightest bins the
errorbars are enhanced. This is because for those bins $b_{1}\sim 1$ and
hence, the universality relation results in $b_{\phi}\sim 0$.

When $b_{1}$ is different from $1$, $f_{\rm NL}$ can be adequately measured.
The perturbations that favours a positive value of $f_{\rm NL}$ in the
Gaussian simulation are also present in the non-Gaussian one. This is shown in
blue symbols for the $f_{\rm NL}=0$ simulation and in and orange symbols for
the $f_{\rm NL}=100$ simulation.

We reduce the variance in $f_{\rm NL}$ by using the matching technique (see
§4.3) and those results are shown as green symbols .

We recover the input value of $f_{\rm NL}$ for both simulations for haloes
with $M_{\rm halo}>1\times 10^{12}\;h^{-1}M_{\odot}$ (see Table 3). For the
measurements in each of the simulations, we find that they have some noise
that causes the measurements to be up to $1.2\sigma$ away from the fiducial
value. However, this noise goes in the same direction in both simulations, and
with the matching technique we are able to cancel much of it, and thus recover
the fiducial value of $\Delta f_{\rm NL}=100$ for haloes with $M_{\rm
halo}>1\times 10^{12}\;h^{-1}M_{\odot}$.

For haloes with masses between $M_{\rm halo}2\times 10^{11}\;h^{-1}M_{\odot}$
and $M_{\rm halo}1\times 10^{12}\;h^{-1}M_{\odot}$, we find some deviations
from the input value. These differences, which we already see at the level of
$b_{\phi}f_{\rm NL}$ in Figure 5, can be interpreted as meaning that the
population of haloes in that mass bin have a $b_{\phi}$ that does not follow
the universality relation.

###  5.4 p as a function of mass

Since we know the input value of $f_{\rm NL}$, we can use the deviations of
our measurements from that value to find how far is the assumption $p=1$ is
from the $p$ measured for the simulation haloes.

We fix the value of $f_{\rm NL}$, and derive $b_{\phi}$ from the previous
constraints on $b_{\phi}f_{\rm NL}$ we got in subsection 5.2. Then, $p$ is
computed using Equation 12. We also checked that we get the same results by
fitting $p$ directly and fixing $f_{\rm NL}$ to its input value, in a similar
way as done in subsection 5.3.

The results are shown in Figure 7, where we compare the different methods of
estimating the variance of the power spectrum and reducing the noise
introduced in subsection 4.4.

Figure 7: PNG-response parameter $p$ as a function of the halo mass, using
different methods to treat variances. Blue and orange symbols show results for
the non-Gaussian simulation, using standard variances and FastPM variances
respectively. Green and red symbols show the same, after applying the matching
with the Gaussian simulation. In solid circles we use the best fit to Equation
14 and in rings we use the extrapolation.

We find that our methodology is consistent as we recover the same results
regardless of the method used for estimating the variance of the power
spectrum or using or not the matching technique.

The clustering of the haloes with $5\times 10^{11}<M_{\rm halo}<2\times
10^{13}\;h^{-1}M_{\odot}$ are well described through the universality relation
$(p=1)$. With respect to the mass bins with $M_{\rm halo}<5\times 10^{11}$, we
find a slight deviation from $p=1$. These lower values of $p$ are related to
the overprediction of the product $b_{\phi}f_{\rm NL}$ (§5.2).

We explore the possibility of $p$ being a constant independent of the halo
mass. We do this by performing a linear fit. We obtain a slope consistent with
$0$ at $1\sigma$, suggesting that a constant $p$ can indeed be assumed for the
haloes. By a weighted average we get $p=0.955\pm 0.013$ (with $\chi^{2}/{\rm
dof}=1.31$) in case we take into account all the bins, and $p=0.950\pm 0.025$
(with $\chi^{2}/{\rm dof}=1.33$) if neglect those where we used the
extrapolations.

##  6 Convergence Tests

The purpose of this section is to show how far we can push our results to low
masses. To do so, we test the convergence of the PNG-UNIT against the LR-UNIT
simulations with a lower resolution and compare also with an alternative way
to measure $b_{\phi}$ using the separate universe approach (§6.1). We also
investigate what is the impact of considering different mass definitions
(§6.2).

Figure 8: The non-Gaussian bias parameter $b_{\phi}$ as a function of the halo
mass for the different sets of simulations. The color of the symbols and the
dashed lines represent the mass resolution. Red: $m_{p}=1.5\times
10^{8}\,[h^{-1}\,M_{\odot}]$. Green: $m_{p}=1.2\times
10^{9}\,[h^{-1}\,M_{\odot}]$. Blue/Dark blue: $m_{p}=1.0\times
10^{10}\,[h^{-1}\,M_{\odot}]$. The dashed vertical lines represent the
threshold of $M_{\rm halo}=20m_{\rm part}$ for each simulation. The pink solid
line is the universality relation $p=1$ for the $b_{1}$ values measured in the
PNG-UNIT. The dashed pink line is an extrapolation of the universality
relation toward lower masses, using $b_{1}$ coming from Tinker et al. (2010).
In solid circles we use the best fit to Equation 14 and in rings we use the
extrapolation. Upper panel: Measurements of $b_{\phi}$ in non-Gaussian
simulations by fitting Equation 10. Lower panel: Measurements of $b_{\phi}$
using the separated universes technique (see §6.1)

###  6.1 Convergence with Supporting simulations

The upper panel of Figure 8 compares the values of $b_{\phi}$ derived from the
$4096^{3}$-particle simulations with those from the LR-UNIT simulations. The
$b_{\phi}$ values agree for bins with $M_{\rm halo}\gtrsim 20m_{\rm part,LR}$.
Below this threshold, the mass function from the LR-UNIT simulations drops as
the we cannot resolve such small haloes (Figure 2). This drop in halo numbers
affects not only the abundances but also the clustering, as demonstrated by
the sudden change in $b_{\phi}$. Thus, we prove that $20$ particles is an
adequate minimum threshold for the halo mass considered for our simulations.

The universality relation ($p=1$) is shown as a solid pink line in Figure 8,
assuming the $b_{1}$ values from the fits of the PNG-UNITsims to the
corresponding mass resolution. This is required as there is not a direct
theoretical relation between $b_{\phi}$ and the halo mass. In order to extend
the relation beyond the mass resolution from the N-body simulations, we have
also used the theoretical values for $b_{1}(M)$ derived using Tinker et al.
(2010) for each bin in halo mass, given the simulation cosmology. This is
shown as dashed pink lines in Figure 8. The two universality relations are
consistent between them.

The derived $b_{\phi}$ values are in agreement with the universality relation,
for haloes with more than $20$ particles,

Although not shown here, we have tested that the results from running LR-UNIT
with 2LPTic initial conditions are consistent with those from LR-UNIT
initialised with FastPM. In a forthcoming paper, we will study more in detail
the impact of using different codes for the initial conditions in simulations
with Primordial non-Gaussianities.

As a further test of the robustness of our results, we compare two alternative
methods for measuring $b_{\phi}$: (i) using the scale-dependent bias as we
have done throuhgout this paper; and (ii) using the so called separated
universes technique. This last technique exploits the fact that $b_{\phi}$ is
defined as the response of the abundance of haloes to the presence of a large
scale perturbation (Equation 11).

The separated universes technique consists in computing $b_{\phi}$ as the
numerical derivative of two simulations with the same initial conditions but
with a slightly different $\mathcal{A}_{s}$, or equivalently, $\sigma_{8}$.
The initial conditions for these simulations are expected to be Gaussian and
thus, these simulations cannot show the effect of the scale-dependent bias.
Following a similar procedure as done in Barreira et al. (2020); Barreira
(2022c), $b_{\phi}$ is obtained as dollows:

$b_{\phi}=\frac{b_{\phi}^{\rm high}+b_{\phi}^{\rm low}}{2}\,,$ (18)

with:

$\displaystyle b_{\phi}^{\rm high}$
$\displaystyle=\frac{4}{\delta\mathcal{A}_{s}N_{h}^{\rm fid}}(N_{h}^{\rm
high}-N_{h}^{\rm fid}),$ (19) $\displaystyle b_{\phi}^{\rm low}$
$\displaystyle=\frac{4}{\delta\mathcal{A}_{s}N_{h}^{\rm fid}}(N_{h}^{\rm
fid}-N_{h}^{\rm low}).$

In Equation 19, $\delta\mathcal{A}_{s}$ is the variation in the amplitude of
the initial power spectrum (which we set to $5\%$); and $N_{h}$ is the number
of haloes within a certain mass bin measured in the simulations. For the
separated universes technique, we estimate the error in $b_{\phi}$,
$\sigma_{\rm SU}(b_{\phi})$, following the same consideration as Barreira et
al. (2020); Barreira (2022c):

$\sigma_{\rm SU}(b_{\phi})\sim|b_{\phi}^{\rm high}-b_{\phi}^{\rm low}|.$ (20)

The lower panel of Figure 8 compares the $b_{\phi}$ values obtained from the
perturbation theory definition for simulations with PNG (Equation 10) and from
using the separated universes technique (Equation 19). Comparing the two
methods for measuring $b_{\phi}$, not only tests the limitations due to mass
resolution effects, but this also serves as a consistency test.

In Figure 8, we find a good agreement between the $b_{\phi}$ values obtained
from applying the separated universes technique to either the low or high
resolution simulations.

The $b_{\phi}$ calculated with the separated universes technique are also
close to the universality relation, displayed as pink lines in Figure 8. These
relations is found from Equation 10, as the boxes of the separate universes
simulations are quite small and a direct measurement of the linear bias
through the halo power spectrum is very contaminated by the non-linearities of
the small scales ($k_{f,\,{\rm HR}}\sim 0.09\;h\,{\rm Mpc}^{-1}$ and
$k_{f,\,{\rm LR}}\sim 0.02\;h\,{\rm Mpc}^{-1}$).

Comparing both panels in Figure 8, we find values of $b_{\phi}$ that are
consistent for haloes with masses above $M_{\rm halo,\;min}\sim 20\,m_{\rm
part}$. The differences between the average $b_{\phi}$ values from the
separate universe technique are within $2\sigma$ of the corresponding values
from the PNG-UNITsims. This comparison is done bin by bin in halo mass. For
the separate universe technique, we estimate a global error in $b_{\phi}$, by
combining the measurements from both the high and low-resolution simulations,
using Equation 18 and Equation 20. With this procedure, we expect to obtain
more statistically reliable errors, compared to those coming from a single set
of simulations.

Finally, we measure $p$ from the separate universe and from the LR-UNIT
simulations. Information on the linear bias, $b_{1}$, is required to obtain
$p$ from Equation 12. For the simulations with initial conditions with PNG, we
obtain similar results for $p$ either using $b_{1}$ coming from Tinker et al.
(2010) or those directly obtained from the simulations. In particular, we
obtain $p=0.959\pm 0.014$ for the PNG-UNIT and $p=0.972\pm 0.017$ for the LR-
UNIT. These values are consistent with our previous results.

In the case of the separate universe technique, for each mass bin, we fit
$b_{\phi}$ to Equation 12 using a least squares-method, with $b_{1}$ coming
from Tinker et al. (2010) for each mass bin. Utilising the separate universe
error in $b_{\phi}$, $\sigma_{\rm SU}(b_{\phi})$ (Equation 20), we get
$p=0.8307\pm 0.0022$. This result is very different to the universality
relation and the $p$ obtained for the PNG-UNIT, $p=0.955\pm 0.013$. We argue
that this suggest a shortcoming in the method used to obtain $\sigma_{\rm
SU}(b_{\phi})$. The separate universe measurement of $p$ is driven by the mass
bins with the smallest errors. Despite only having three separate universe
simulations for each mass resolution, the error bars for some of the mass bins
are smaller than for PNG-UNITsims, where we are using several variance
reduction techniques (section 4). Moreover, $\chi^{2}/{\rm dof}=519/20$ for
the $p$ obtained with the separate universe technique.

Ignoring the separate universe error bars, by not weighting the bins by their
$\sigma(b_{\phi})$ in the least squares fitting procedure, we find $p=1.047\pm
0.057$. This value is consistent with the universality relation and it is less
than $2\sigma$ away the results shown in subsection 5.2 and subsection 5.4.

###  6.2 Dependency with mass definition

If $b_{\phi}$ depends on the halo mass only by its dependence on $b_{1}$, then
the relationship between these two quantities should be independent of the
mass definition used for haloes.

Throughout this work, the halo masses have been obtained by computing the
enclosed mass within a spherical overdensity up to $200$ times the critical
density $\rho_{c}$ at a given $z$, $M_{200c}$. Here, we redo the analysis to
obtain $p$, using other mass definitions, namely: $M_{500c}$, $M_{200b}$ and
$M_{\rm vir}$ . The radii used for these mass definitions are determined by
the threshold density achieved within the spherical overdensity, which reaches
$500$ the critical density $\rho_{c}$ for $M_{500c}$ and $200$ times the
background density $\rho_{b}$ for $M_{200b}$. For the virial mass $M_{\rm
vir}$ the spherical overdensity is determined by the virial radius $r_{\rm
vir}$ (Bryan & Norman, 1998).

For a given halo, $M_{200c}\geq M_{500c}$. Thus, if we apply the same mass
cuts for defining the bin, we expect that $b_{1}(M_{200c})<b_{1}(M_{500c})$,
as if we compare the $M_{200c}$ of the haloes selected in each case, we are
taking more massive haloes when we apply the cut in $M_{500c}$. However,
although we expect that $b_{\phi}(M_{200c})\neq b_{\phi}(M_{500c})$, if
$b_{\phi}$ only depends on the mass through $b_{1}$, we would expect that $p$
does not vary with the mass definition.

We present our findings in Figure 9 which compares five mass definitions
outputted by Rockstar: our fiducial $M_{200c}$, $M_{500c}$, $M_{200b}$ and
$M_{\rm vir}$. There are several things to point out. We find that bins
defined according to $M_{200b}$ and $M_{\rm vir}$ are practically equal
between them and very similar to $M_{200c}$. Moreover, they have almost the
same $b_{1}$ although this is not shown here. For these three cases the
results are consistent with our findings in subsection 5.4.

However, the most interesting result comes from the bins defined according to
$M_{500c}$. We find a systematic suppression of $p$ below the predicted value
from the universality relation and compared with the other mass definitions,
although the errorbars cover $p=1$. This result suggests that $b_{\phi}$ for
dark matter haloes may not depend solely on the mass through $b_{1}$, but also
on other parameters. The $b_{\phi}$ parameter is defined as how the abundance
of the dark-matter haloes responds to a large-scale perturbation. This may not
only affect the amplitude of the primordial perturbations (and hence the
abundance of the high-mass haloes) but also could affect to other properties
of the halo as suggested by recent works (Lazeyras et al., 2023; Lucie-Smith
et al., 2023; Fondi et al., 2023). However, investigating this dependence, is
beyond the scope of this paper and is left for future research.

Additionally, we also tested the impact on $b_{\phi}$ if we consider not only
the main haloes, but also the subhaloes identified by Rockstar. This might be
important if one is about to run subhalo abundance matching (SHAMs) algorithms
to populate the simulation with galaxies. What we find is that when the
substructure is included in the analysis, this induces variations not only in
the linear bias $b_{1}$ of the mass bins but also in $p$. We leave a detailed
analysis of this for a future work.

Figure 9: The PNG bias parameter $p$ as a function of $M_{\rm halo}$ for
different mass definitions. The green symbols represent our choice for the
fiducial analysis $(M_{200c})$, while he gray dashed line represents the
universality relation ($p=1$). Masses are shifted for visualisation purposes.

##  7 Robustness Tests

Here, we test the impact of the assumptions we have made in our analysis to
obtain $b_{\phi}$ or $p$. In particular, we study the impact of the scale cut,
$k_{\rm max}$ (§7.1), and the variance estimation methods described in
subsection 4.4.

###  7.1 Choice of $k_{\rm max}$

Equation 10 is only valid in the linear regime, and thus, we must discard all
modes with $k$ greater than a certain $k_{\rm max}$ to avoid introducing non-
linear modes that are not properly described by our model. As the perturbation
theory breaks beyond $k_{\rm NL}\sim 0.3\;h\,{\rm Mpc}^{-1}$ (Desjacques et
al., 2018), $k_{\rm max}$ should be less than this value. Even before this
threshold, we are at the risk of incorporating some mildly non-linear modes
which are not adequately described by our model.

To determine the optimal $k_{\rm max}$, we first examine the limits of the
linear theory using only the Gaussian simulation. In this case, Equation 10 is
reduced to $P_{h,h}(k)=b_{1}^{2}P_{m,m}(k)$. The leading correction of the
non-linear bias to $P_{h,h}(k)$ has a dependence with $Ak^{2}$, where $A$ is a
nuisance parameter (Lazeyras & Schmidt, 2018; Lazeyras et al., 2023). To
establish $k_{\rm max}$, we compare the fits obtained using only $b_{1}$ with
those obtained using $b_{1}+Ak^{2}$ after marginalising over the nuisance
parameter $A$, as a function of $k_{\rm max}$. We find that the differences in
the posteriors for $b_{1}$ obtained using both methods agreed for $k_{\rm
max}\leq 0.1\;h\,{\rm Mpc}^{-1}$ for all mass bins.

We proceed in a similar way with the non-Gaussian simulation. Starting from
our fiducial choice of $k_{\rm max}=0.1\;h\,{\rm Mpc}^{-1}$, we tested both
lower and higher values. Our results showed that varying the $k_{\rm max}$ has
a minimal impact on $p$, with all central values remaining within $1\sigma$
when $k_{\rm max}$ was varied from $0.07$ to $0.15\;h\,{\rm Mpc}^{-1}$. This
is shown in panel (a) of Figure 10. The error bars shrank by approximately
$\sim 20\%$ throughout this range when considering a larger $k_{\rm max}$.

Figure 10: The PNG bias parameter $p$ as a function of $M_{\rm halo}$ for
different analysis choices. The green symbols represent the fiducial choices,
while the gray dashed line represents the universality relation ($p=1$). In
solid circles we use the best fit to Equation 14 and in rings we use the
extrapolation. Masses are shifted for visualisation purposes. (a) Variation in
$p$ with changes $k_{\rm max}$ (see §7.1). (b) Constrains on $p$ by
considering different methods to compute $\sigma(P(k))$, namely: Standard
$\sigma_{S}(P(k))$ (blue), direct fit to FastPM mocks (orange), best fit to
Equation 14 (green) and extrapolation of Equation 14 to the mass bins (red)
(see §7.2).

Going beyond the purely linear modes may bias the constraints on $b_{1}$
(since those modes are not properly modelled) and hence in $b_{\phi}$ and $p$.
We have therefore set the cutoff at $k_{\rm max}=0.1\;h\,{\rm Mpc}^{-1}$. On
the one hand, we have shown that the results are consistent even if this value
is varied by $\sim 50\%$. On the other hand, this is a standard choice, and it
is more conservative than keeping on pushing to higher $k$ where we might have
to deal with non-linear terms.

###  7.2 Power spectrum variance from FastPM

An accurate estimation of the variance of the halo power spectrum,
$\sigma(P(k))$, is critical when we do the parameter estimation, given that
this quantity enters in the computation of the likelihood. However, we know
that the standard approach of Equation 13 overestimates the variance of our
simulations as they have fixed initial conditions. In subsection 4.2, we
discussed how we estimate $\sigma(P(k))$ from the FastPM mocks which take into
account the effect of the fixed ICs. Moreover, we show how by fitting
$\sigma(P(k))/P(k)$ to Equation 14 we can even extrapolate to masses where the
mocks do not reproduce accurately the N-body simulation.

Here, we check the impact of considering different ways of computing the power
spectrum variance have on $p$. In particular we consider the following methods
(see §4.2):

  * •

Standard variances, $\sigma_{S}(P(k))$, computed from Equation 13 (blue).

  * •

Direct FastPM variances, $\sigma_{F}(P(k))$, obtained by computing directly
the standard deviation of the power spectrum of the fast simulations (orange).

  * •

Best fit of FastPM variances using Equation 14 (green).

  * •

Extrapolation of Equation 14 to the mean halo mass of the bins (red).

The way we have computed the variance of the power spectrum for our baseline
analysis is the third case, where we use the best fit of Equation 14 for each
mass bin.

We present the results of our comparison in panel (b) of Figure 10, which
shows the PNG bias parameter $p$ as a function of $M_{\rm halo}$ for the
different methods of estimating the power spectrum variance described above.
For bins $0$-$3$ we do not have enough mass resolution with the FastPM mocks,
so only results assuming the Standard variances and the extrapolation of
$\sigma_{F}(P(k))$ are shown. We find that regardless of the method used to
estimate the variance, we recover the same results within $1\sigma$. Moreover,
the FastPM approach is more accurate than the standard approach as it takes
into account the fixed initial conditions in our simulations. Additionally,
the best fit of $\sigma_{F}(P(k))/P(k)$ provides almost the same results
compared to using the extrapolation. These results are also consistent with
using directly the variances obtained from FastPM even considering their small
errorbars.

While we cannot test and compare our fitting procedure outside the mass domain
of the FastPM mocks, we remark that our results are still consistent within
$1\sigma$ with the Standard approach.

##  8 Constraining $f_{\rm NL}$ with priors on $p$

We devote this section to study how setting a prior in $p$ (or equivalently,
on $b_{\phi}$) affects the constraints on $f_{\rm NL}$. In the case where
these two parameters are left completely free in the MCMC, the constraints are
very weak because with the scale-dependent bias we are only sensitive to the
product $b_{\phi}f_{\rm NL}$.

In subsection 8.1 we find that the constraints in $f_{\rm NL}$ deteriorate as
we increase the size of the prior $\sigma(p)$. In subsection 8.2 we show how
these results can affect the measurement of $f_{\rm NL}$ in upcoming
cosmological surveys.

Here, we assume a Gaussian prior for p, around $p_{\rm prior}$ with width
$\sigma(p)$. To implement this in the MCMC we modify the Gaussian likelihood
by introducing this new term in the $\chi^{2}$:

$\chi^{2}=\sum^{N_{k}}_{i=0}\frac{\left(P_{\rm haloes}(k_{i})-P_{\rm
model}(k_{i},b_{\phi}f_{\rm
NL},b_{1})\right)^{2}}{\sigma^{2}(P(k_{i}))}+\frac{p-p_{\rm
prior}}{\sigma(p)}\;.$ (21)

###  8.1 Effect of the prior size on $\sigma(f_{\rm NL})$

Figure 11: Statistical error in $f_{\rm NL}$ as a function of the size of the
prior on $p$. The mock data vector for this figure is generated using the
linear power spectrum computed with CAMB in Equation 10 and fixing $b_{1}=1.4$
and $f_{\rm NL,\;fid}=0$ while we consider 3 values for $p$, which are shown
in different colours, as indicated in the legend. The horizontal dashed lines
show the case when we perform the fits fixing $p$ in the model to its fiducial
value in each case. The dashed vertical line displays the value of the error
in the value of $p$ we obtain for halos with mass between $5\times 10^{11}$
and $1\times 10^{12}\,h^{-1}M_{\odot}$ (bin 4 in Table 3).

In this section, we investigate how the size of a prior on $p$ affects the
determination of $f_{\rm NL}$ compared to assuming a fixed value for $p$. To
this end, we generate 9 mock galaxy power spectra using Equation 10. The
matter power spectrum is computed with CAMB. We fix $b_{1}=1.4$ and consider
three values for $f_{\rm NL}$, $0$, $-12$ and $50$. $b_{\phi}$ is obtained
from Equation 12, where we consider also three different values for $p$,
namely $0.5$, $1$ and $1.6$. In this way, we run an MCMC to determine $f_{\rm
NL}$ and $\sigma(f_{\rm NL})$ for each case.

The fits are done in two ways, by fixing $p$ to the fiducial value and by
allowing $p$ to vary with a Gaussian prior around the fiducial value. When we
fix $p$, we are assuming that we know the correct $p$. In this ideal case, we
only have two free parameters, $b_{1}$ and $f_{\rm NL}$, and the constraints
will be the tightets.

When we allow $p$ to vary, we have to assume a prior, otherwise we will not be
able to constrain $p$ and $f_{\rm NL}$ simultaneously due to the degeneracy of
these parameters. Thus, we establish a Gaussian prior around $p_{\rm fid}$
with a width of $\sigma(p)$, and the free parameters of the MCMC will be
$b_{1}$, $f_{\rm NL}$ and $p$.

We find that the uncertainties in $f_{\rm NL}$ arise primarily from two
factors: the size of the prior on $b_{\phi}$ or $p$ and the specific value of
$b_{\phi}$. Increasing the size of the priors in $p$ leads to larger
$\sigma(f_{\rm NL})$ due to the perfect degeneracy between $f_{\rm NL}$ and
$b_{\phi}$. However, when the prior in $p$ becomes sufficiently narrow,
approaching the limit where $b_{\phi}$ is fixed, the uncertainties in $f_{\rm
NL}$ become similar in both cases. This behaviour is demonstrated in Figure
11.

For the fiducial values we have considered here we find that if
$\sigma(p)\lesssim 0.1$, the constraints on $f_{\rm NL}$ obtained by using the
priors are equivalent to fixing $b_{\phi}$ to its fiducial value. Going beyond
that threshold, the constraints deteriorate rapidly due to the degeneracy
between $f_{\rm NL}$ and $b_{\phi}$. This results in an increase of
$\sigma(f_{\rm NL})$ of several orders of magnitude.

At large $\sigma(p)$, Figure 11 shows that $\sigma(f_{\rm NL})$ saturates at a
large value. This is a consequence of the choice for the size of the parameter
space in $f_{\rm NL}$ that we allow to be explored when running the MCMC
chains. The larger the range, the larger the saturation value.

We note that the threshold at $\sigma(p)$ depends on the particular value of
$b_{\phi}$ (and hence $b_{1}$ and $p$) and $f_{\rm NL}$. In this case, we find
significant variations on $\sigma(f_{\rm NL})$ when $p$ is varied. As $b_{1}$
is fixed for the three cases, the lower $p$ is, the larger $b_{\phi}$ is. As
$b_{\phi}$ increases, the constraints on $f_{\rm NL}$ become tighter. If
$b_{\phi}$ gets close to zero, then the constraints on $f_{\rm NL}$ are
weaker. Meanwhile, this relation between the value of $b_{\phi}$ and the
possible constrain on $f_{\rm NL}$ is what matters for cosmological surveys.

Similar trends to those presented in Figure 11 are found for $f_{\rm NL}=-12$
and $f_{\rm NL}=50$.

###  8.2 Forecast for a DESI-like survey

| ELG hosts $(b_{\rm 1,\;ELG}=1.15\pm 0.02)$ | LRG hosts $(b_{\rm 1,\;LRG}=2.16\pm 0.02)$ | QSO hosts $(b_{\rm 1,\;QSO}=1.64\pm 0.12)$
---|---|---|---
| universality | prior (bin 4) | universality | prior (bin 8) | universality | prior (bin 6)
$f_{{\rm NL},{\rm fid}}=0$ | $-1\pm 20$ | $0\pm 12$ | $-0.5\pm 4.5$ | $-0.4\pm 5.6$ | $-2\pm 21$ | $-4\pm 20$
$f_{{\rm NL},{\rm fid}}=-12$ | $-21\pm 19$ | $-14\pm 12$ | $-11.2\pm 4.0$ | $-13.7\pm 6.2$ | $-15\pm 22$ | $-16\pm 22$
$f_{{\rm NL},{\rm fid}}=50$ | $84\pm 25$ | $51\pm 17$ | $44.3\pm 6.5$ | $54\pm 15$ | $51\pm 19$ | $48\pm 19$
Table 4: Forecasts on $f_{\rm NL}$ for haloes hosting the different tracers
targeted by DESI. In the universality columns, $b_{\phi}$ is fixed to Equation
12 with $p=1$. In the prior columns, we use our results for $p$ to conduct the
fits with: $p_{\rm LRG}=1.11\pm 0.27$, $p_{\rm ELG}=0.895\pm 0.055$ and
$p_{\rm QSO}=0.973\pm 0.077$ respectively.

We now study what is the impact on the determination of $f_{\rm NL}$ in the
context of the current galaxy surveys if we assume the universality relation
or consider instead a prior on $p$ or $b_{\phi}$ as the ones obtained in this
work.

For this forecast, we assume that the survey properties are similar to DESI,
whose footprint will cover an area in the sky of around $14.000\deg^{2}$ (DESI
Collaboration et al., 2023). We consider the three types of cosmological
tracers observed in this survey, that is Luminous Red Galaxies (LRGs),
Emission-line galaxies (ELGs) and Quasars (QSOs). Although the redshift
distribution of these tracers is very wide and different among them, we limit
our analysis only to the range between $z=0.8$ and $z=1.1$. We assume that our
simulation can describe the clustering of the haloes hosting these galaxies at
redshift $z=1$. The combination of this redshift range with the total survey
area results in an effective volume of $V_{\rm eff,\;DESI}\sim
11\,(h^{-1}\,{\rm Gpc})^{3}$. This yields $k_{\rm min}=2\pi/V_{\rm
DESI}^{1/3}\sim 2.83\times 10^{-3}\,h\,{\rm Mpc}^{-1}$. We consider modes up
to $k_{\rm max}=0.1\,h\,{\rm Mpc}^{-1}$.

The properties of DESI LRGs, ELGs and QSOs have been studied in detail for the
Early Data Release (Yu et al., 2023; Yuan et al., 2023; Rocher et al., 2023;
Prada et al., 2023). We focus on Yu et al. (2023) as they based their analysis
on the simulation we have used in this work, the Gaussian UNIT simulation. Yu
et al. (2023) find that DESI EDR LRGs, ELGs and QSOs populate mostly haloes
whose mass is $\log_{10}(M_{\rm halo,\,LRG})=13.16$, $\log_{10}(M_{\rm
halo,\,ELG})=11.90$ and $\log_{10}(M_{\rm halo,\,QSO})=12.66$ respectively.
This corresponds to our mass bins 8 , 4 and 6 respectively (see Table 3). In
addition, the number densities they find for the considered redshift range are
$n_{\rm LRG}=4.36\times 10^{-4}\;h^{3}\,{\rm Mpc}^{-3}$, $n_{\rm
ELG}=1.047\times 10^{-3}\;h^{3}\,{\rm Mpc}^{-3}$ and $n_{\rm QSO}=2.9\times
10^{-5}\;h^{3}\,{\rm Mpc}^{-3}$. For the linear bias we assume the values from
Yu et al. (2023), which for the considered redshift range and tracers they are
given by $b_{1,\;{\rm ELG}}=1.15\pm 0.02$ , $b_{1,\;{\rm LRG}}=2.16\pm 0.02$
and $b_{1,\;{\rm QSO}}=1.64\pm 0.12$. We emphasise that our forecast is based
on the clustering of the haloes hosting these tracers.

We generate the mock halo power spectrum using CAMB in Equation 10 and the
corresponding covariance is calculated using Equation 13, considering the
survey parameters described above for each tracer. We consider three choices
for the fiducial value of $f_{\rm NL}$: $f_{\rm NL,\,fid\,1}=0$, $f_{\rm
NL,\,fid\,2}=-12$ and $f_{\rm NL,\,fid\,2}=50$. The choice of $f_{\rm
NL,\,fid\,1}=0$ is motivated by the current tightest constraint for $f_{\rm
NL}$ by Planck Collaboration et al. (2019); $f_{\rm NL,\,fid\,1}=-12$ and
$f_{\rm NL,\,fid\,3}=50$ are motivated by being the two recent constraints
coming from galaxy clustering analysis (Mueller et al., 2022; Rezaie et al.,
2023). Regarding $b_{\phi}$, we assume it follows the form given Equation 12.
We use the values of $p$ measured in this work: $p_{\rm LRG}=1.11\pm 0.27$,
$p_{\rm ELG}=0.895\pm 0.055$ and $p_{\rm QSO}=0.973\pm 0.077$

We conduct an MCMC fitting procedure to find out the effect that assuming a
prior on $b_{\phi}$ or not has on the measurement of $f_{\rm NL}$. We employ
the same methodology as for our main results. In this analysis, we keep all
the cosmological parameters fixed, except for $b_{1}$, $f_{\rm NL}$ and
$b_{\phi}$. When fixing $b_{\phi}$ to the universality relation, we define the
$\chi^{2}$ using Equation 17 . When considering $p$ as a free parameter we
incorporate a Gaussian prior centred in $p_{\rm LRG}$, $p_{\rm ELG}$ and
$p_{\rm QSO}$ for each case value and a width given by our errors.

The summary of our findings are shown in Table 4. The constraints for haloes
hosting LRGs are stronger than for the equivalent for ELGs and QSOs. This is
mainly because they have the largest linear bias in the redshift range we are
considering compared to the other tracers, and therefore the scale-dependent
bias signal is larger. Assuming our priors, we recover the fiducial $f_{\rm
NL}$ value within $<0.2\sigma$. Fixing $b_{\phi}$ to the universality
relation, while not inducing any bias in the case of QSOs, does induce a
deviation of $0.47\sigma$ for $f_{\rm NL\;fid}=-12$ and $1.36\sigma$ for
$f_{\rm NL\;fid}=50$ in the case of ELGs, were the value of $p$ is a $\sim
10\%$ lower than $p=1$. For LRGs we find also a difference of $0.87\sigma$ for
the case with $f_{\rm NL\;fid}=50$.

Regarding $\sigma(f_{\rm NL})$, we observe three different behaviours. For the
halos hosting ELGs, using our prior provides tighter constraints compared to
fixing $p=1$ while for the halos hosting LRGs is precisely the opposite. We
claim that this is not due to the different size of the prior but the central
value it has ($p_{\rm LRG}=1.1$ vs $p_{\rm ELG}=0.89$).

In subsection 8.1 we have shown that if the width of the Gaussian prior on $p$
are small enough, the effect of the prior in $\sigma(f_{\rm NL})$ is
equivalent to fixing $p$ to that value. And if the prior is too wide, then we
recover the degeneracy between $p$ and $f_{\rm NL}$. The priors we have in $p$
allow us to break this degeneracy and then, the lower $p$ is, the larger
$b_{\phi}$ is and hence the constraints on $f_{\rm NL}$ become tighter.

##  9 Summary and Conclusions

Detecting Primordial non-Gaussianities (PNG) would have a major impact on how
we understand inflation. The key to being able to measure PNG in the near
future is to have a good understanding of the scale dependent bias parameter,
$b_{\phi}$, and how it varies for different cosmological tracers.

In this paper we study the scale dependent bias, $b_{\phi}$, for dark matter
haloes with a range of masses from $2\times 10^{10}$ to $5\times
10^{14}\,h^{-1}M_{\odot}$ . To achieve this we have developed a new full
N-body simulation with non-Gaussian initial conditions $(f_{\rm NL}=100)$,
PNG-UNITsims. This is the largest simulation to date with non-Gaussian initial
conditions. The simulation has $4096^{3}$ dark matter particles in a
$1\;(h^{-1}\,{\rm Gpc})^{3}$ box, which was run using the L-Gadget2 code. The
cosmology assumed is based on Planck Collaboration et al. (2016) (see Table
1). We have implemented local PNG with an amplitude, given by $f_{\rm
NL}=100$. Such a large value of $f_{\rm NL}$ has been chosen to have a large
signal from PNG and reduce the errors on $b_{\phi}$. In addition to the main
simulation, we also present a set of low-resolution UNIT simulations and a
suite of 100 fast mocks generated using FastPM, half of them with $f_{\rm
NL}=100$ and the other half with $f_{\rm NL}=0$.

This suite of non-Gaussian simulations are a unique laboratory to study galaxy
clustering models in the presence of local PNGs. The simulations in this suite
have effective volume comparable to current spectroscopic galaxy surveys. The
reference simulation, PNG-UNITsims, has a mass resolution high enough to
resolve the dark matter haloes expected to host the galaxies currently being
targeted by surveys such as DESI and Euclid. Thus, this suite is a powerful
tool to develop and test analysis pipelines to measure $f_{\rm NL}$ with
current spectroscopic galaxy surveys.

In this first study of the PNG-UNITsims project, we constrain the non-Gaussian
bias parameter, $b_{\phi}$, for mass selected dark matter haloes using these
simulations. We make use of the method proposed in Avila & Adame (2023) of
fixed & matched initial conditions to reduce significantly the error on
$b_{\phi}$. The amplitudes of the modes in the initial conditions are fixed to
their expectation values and the phases are matched to one of the realisations
of the Gaussian UNIT simulations (Chuang et al., 2019), following Avila &
Adame (2023).

As a novelty, we use fast mocks generated with FastPM instead of full N-body
simulations in order to compute the reduced variance of the power spectra due
to fixing the initial conditions and the correlation coefficientes needed for
the reduction of the error due to the matched phases.

We measure the scale dependent bias in the halo-halo power spectrum for mass
selected haloes for both, the Gaussian and the non-Gaussian simulation. From
that, we infer $b_{1}$ and $b_{\phi}f_{\rm NL}$ and we match the measurements
in both simulations. Moreover, given we know the input values of $f_{\rm NL}$,
we can derive $b_{\phi}$ or $p$. Our findings can be summarised in the
following:

  * •

We measure $b_{\phi}f_{\rm NL}$ with high precision using the fix and match
technique (§4.3) by combining the information of the Gaussian and the non-
Gaussian simulation. This allows us to cancel out a great part of the cosmic
variance (§5.2).

  * •

Combining the information of $b_{\phi}f_{\rm NL}$ from all the mass bins and
by fitting $b_{\phi}$ as a function of $b_{1}$ and $p$, we find a preferred
value of $p=0.955\pm 0.013\,(1\sigma)$.

  * •

By assuming the universality relation (ie. p=1), we recover the input values
of $f_{\rm NL}$ for our simulations within the $1\sigma$ errorbars for haloes
between $1\times 10^{12}$ and $1\times 10^{14}\,h^{-1}M_{\odot}$. For larger
halo masses we do not recover $f_{\rm NL}=0$ in the Gaussian UNIT simulation
at $3\sigma$, even using the standard variances for $P(k)$. This suggest that
we are underestimating the $\sigma(P(k))$ for those bins so we do not consider
them in the rest of the analysis. For haloes with masses between $5\times
10^{10}\,h^{-1}M_{\odot}$ and $1\times 10^{12}\,h^{-1}M_{\odot}$ we find
deviations at $>1\sigma$ from the input value.

  * •

We fixing $f_{\rm NL}$ to the input values we measure $p$. For haloes with
masses larger than $1\times 10^{12}\,h^{-1}M_{\odot}$ and between $2\times
10^{10}\,h^{-1}M_{\odot}$ and $5\times 10^{10}\,h^{-1}M_{\odot}$ we find $p$
to be consistent with the universality relation. For haloes with masses
between $5\times 10^{10}\,h^{-1}M_{\odot}$ and $1\times
10^{12}\,h^{-1}M_{\odot}$ we find $p$ a $10\%$ below the universality
prediction.

We assure the strength of our results for $M_{\rm halo}\gtrsim 20m_{\rm
part}\sim 2\times 10^{10}\,h^{-1}M_{\odot}$ by checking convergence with lower
resolution simulations and an alternative methods to constrain $b_{\phi}$,
such as the separate universe approach.

We have tested the robustness of the above results to variations in the choice
of the scale cuts (§7.1), the way we compute the variances of the power
spectrum (§7.2) and the mass definition used for the mass cuts (§6.2). None of
the above has a significant impact.

We also check how changes our results if we consider a different the mass
definition. We find that our measurements of $p$ for $M_{\rm 500c}$ tend to be
lower compared with the other mass definitions, although they are consistent
with the errorbars. We argue that when we switch from $M_{200c}$ to $M_{500c}$
while keeping the mass limits of the bins constant, other properties apart
from the mass of the haloes considered are changing. This shows the need of a
better understanding of the impact of those properties on $b_{\phi}$.

We show that the priors we build with the PNG-UNITsims are:

  1. 1.

Necessary to get unbiased constraints on $f_{\rm NL}$ (with respect to the
universality that may give biased results).

  2. 2.

Sufficiently tight to allow us to constrain $f_{\rm NL}$ with DESI without
increasing the error with respect to the universality assumption

This is done based on the assumption that DESI tracers can be identified
directly with our halo mass bins in our PNG-UNIT simulation snapshot at $z=1$.
In the future we plan to improve these priors by exploring the variation with
redshift and other halo properties. We also plan to populate our simulation
with galaxies and QSOs.

This shows the potential of the PNG-UNITsim project to solve one of the main
concerns on theoretical systematics to measure/constrain PNG with DESI or
Euclid. Whereas we have focused here in the galaxy/halo power spectrum and the
scale-dependent bias, this simulation is a unique tool to test/develop other
ways to constrain PNG with cosmological surveys. This may include
configuration space correlation functions, higher order statistics, such us
the bispectrum or trispectrum, statistics of the PDF such as counts-in-cells,
density peaks, galaxy cluster counts, etc.

## Data Availability

The data will be shared under email request to the corresponding author.

###### Acknowledgements.

We thank V. Yankelevich O. Hahn, M. Manera and R. Scoccimarro for useful
discussions on how to solve some early problems with initial conditions. This
work has been supported by Ministerio de Ciencia e Innovación (MICINN) under
the following research grants: PID2021-122603NB-C21 (AGA, VGP, GY and AK),
PID2021-123012NB-C41 (SA) and PID2021-123012NB-C43 (JGB). SA and VGP have been
or are supported by the Atracción de Talento Contract no. 2019-T1/TIC-12702
granted by the Comunidad de Madrid in Spain. JGB has also been supported by
the Centro de Excelencia Severo Ochoa Program CEX2020-001007-S at IFT. IFAE is
partially funded by the CERCA program of the Generalitat de Catalunya. AK is
further thanks Embrace for ‘gravity’. The PNG-UNIT simulation has been run
thanks to the computer resources at MareNostrum granted by the Red Española de
Supercomputación and the technical support provided by Barcelona
Supercomputing Center (RES-AECT-2021-3-0004 and RES-AECT-2022-1-0007).

## References

  * Alam et al. (2021) Alam, S., de Mattia, A., Tamone, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 504, 4667
  * and David J. Bacon et al. (2020) and David J. Bacon, Battye, R. A., Bull, P., et al. 2020, Publications of the Astronomical Society of Australia, 37
  * Angulo & Pontzen (2016) Angulo, R. E. & Pontzen, A. 2016, Monthly Notices of the Royal Astronomical Society: Letters, 462, L1
  * Avila & Adame (2023) Avila, S. & Adame, A. G. 2023, Monthly Notices of the Royal Astronomical Society, 519, 3706
  * Avila et al. (2018) Avila, S., Crocce, M., Ross, A. J., et al. 2018, Monthly Notices of the Royal Astronomical Society, 479, 94
  * Avila et al. (2020) Avila, S., Gonzalez-Perez, V., Mohammad, F. G., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 5486
  * Avila et al. (2021) Avila, S., Vos-Ginés, B., Cunnington, S., et al. 2021, Monthly Notices of the Royal Astronomical Society, 510, 292
  * Bardeen et al. (1986) Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15
  * Barreira (2022a) Barreira, A. 2022a, Journal of Cosmology and Astroparticle Physics, 2022, 013
  * Barreira (2022b) Barreira, A. 2022b, Journal of Cosmology and Astroparticle Physics, 2022, 057
  * Barreira (2022c) Barreira, A. 2022c, Journal of Cosmology and Astroparticle Physics, 2022, 033
  * Barreira et al. (2020) Barreira, A., Cabass, G., Schmidt, F., Pillepich, A., & Nelson, D. 2020, Journal of Cosmology and Astroparticle Physics, 2020, 013
  * Baumann et al. (2009) Baumann, D., Jackson, M. G., Adshead, P., et al. 2009, in AIP Conference Proceedings (AIP)
  * Behroozi et al. (2012) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2012, The Astrophysical Journal, 762, 109
  * Biagetti et al. (2017) Biagetti, M., Lazeyras, T., Baldauf, T., Desjacques, V., & Schmidt, F. 2017, Monthly Notices of the Royal Astronomical Society, 468, 3277
  * Bryan & Norman (1998) Bryan, G. L. & Norman, M. L. 1998, The Astrophysical Journal, 495, 80
  * Byrnes & Choi (2010) Byrnes, C. T. & Choi, K.-Y. 2010, Advances in Astronomy, 2010, 1
  * Cabass et al. (2022) Cabass, G., Ivanov, M. M., Philcox, O. H., Simonović, M., & Zaldarriaga, M. 2022, Physical Review D, 106
  * Castorina et al. (2019) Castorina, E., Hand, N., Seljak, U., et al. 2019, Journal of Cosmology and Astroparticle Physics, 2019, 010
  * Chan et al. (2022) Chan, K. C., Ferrero, I., Avila, S., et al. 2022, Monthly Notices of the Royal Astronomical Society, 511, 3965
  * Chartier et al. (2021) Chartier, N., Wandelt, B., Akrami, Y., & Villaescusa-Navarro, F. 2021, Monthly Notices of the Royal Astronomical Society, 503, 1897
  * Chuang et al. (2019) Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, Monthly Notices of the Royal Astronomical Society, 487, 48
  * Cochrane et al. (2017) Cochrane, R. K., Best, P. N., Sobral, D., et al. 2017, Monthly Notices of the Royal Astronomical Society, 469, 2913
  * Creminelli & Zaldarriaga (2004) Creminelli, P. & Zaldarriaga, M. 2004, Journal of Cosmology and Astroparticle Physics, 2004, 006
  * Dalal et al. (2008) Dalal, N., Doré , O., Huterer, D., & Shirokov, A. 2008, Physical Review D, 77
  * DeRose et al. (2023) DeRose, J., Chen, S.-F., Kokron, N., & White, M. 2023, Journal of Cosmology and Astroparticle Physics, 2023, 008
  * DESI Collaboration et al. (2023) DESI Collaboration, Adame, A. G., Aguilar, J., et al. 2023, arXiv e-prints, arXiv:2306.06308
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, The DESI Experiment Part I: Science,Targeting, and Survey Design
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Physics Reports, 733, 1
  * Desjacques et al. (2009) Desjacques, V., Seljak, U., & Iliev, I. T. 2009, Monthly Notices of the Royal Astronomical Society, 396, 85
  * Ding et al. (2022) Ding, Z., Chuang, C.-H., Yu, Y., et al. 2022, Monthly Notices of the Royal Astronomical Society, 514, 3308
  * Ezquiaga et al. (2023) Ezquiaga, J. M., García-Bellido, J., & Vennin, V. 2023, Phys. Rev. Lett., 130, 121003
  * Feldman et al. (1994) Feldman, H. A., Kaiser, N., & Peacock, J. A. 1994, The Astrophysical Journal, 426, 23
  * Feng et al. (2016) Feng, Y., Chu, M.-Y., Seljak, U., & McDonald, P. 2016, Monthly Notices of the Royal Astronomical Society, 463, 2273
  * Fondi et al. (2023) Fondi, E., Verde, L., Villaescusa-Navarro, F., et al. 2023, Taming assembly bias for primordial non-Gaussianity
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Giannantonio et al. (2012) Giannantonio, T., Porciani, C., Carron, J., Amara, A., & Pillepich, A. 2012, Monthly Notices of the Royal Astronomical Society, 422, 2854
  * Giannantonio et al. (2014) Giannantonio, T., Ross, A. J., Percival, W. J., et al. 2014, Phys. Rev. D, 89, 023511
  * Gonzalez-Perez et al. (2018) Gonzalez-Perez, V., Comparat, J., Norberg, P., et al. 2018, MNRAS, 474, 4024
  * Grossi et al. (2009) Grossi, M., Verde, L., Carbone, C., et al. 2009, Monthly Notices of the Royal Astronomical Society, 398, 321
  * Hamaus et al. (2011) Hamaus, N., Seljak, U., & Desjacques, V. 2011, Physical Review D, 84
  * Hand et al. (2018) Hand, N., Feng, Y., Beutler, F., et al. 2018, The Astronomical Journal, 156, 160
  * Ho et al. (2015) Ho, S., Agarwal, N., Myers, A. D., et al. 2015, Journal of Cosmology and Astroparticle Physics, 2015, 040
  * Hockney & Eastwood (1981) Hockney, R. W. & Eastwood, J. W. 1981, Computer Simulation Using Particles
  * Kaiser (1984) Kaiser, N. 1984, ApJ, 284, L9
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, Monthly Notices of the Royal Astronomical Society, 457, 4340
  * Knebe et al. (2022) Knebe, A., Lopez-Cano, D., Avila, S., et al. 2022, Monthly Notices of the Royal Astronomical Society, 510, 5392
  * Kokron et al. (2022) Kokron, N., Chen, S.-F., White, M., DeRose, J., & Maus, M. 2022, Journal of Cosmology and Astroparticle Physics, 2022, 059
  * Komatsu & Spergel (2001) Komatsu, E. & Spergel, D. N. 2001, Phys. Rev. D, 63, 063002
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193
  * Lazeyras et al. (2023) Lazeyras, T., Barreira, A., Schmidt, F., & Desjacques, V. 2023, Journal of Cosmology and Astroparticle Physics, 2023, 023
  * Lazeyras & Schmidt (2018) Lazeyras, T. & Schmidt, F. 2018, Journal of Cosmology and Astroparticle Physics, 2018, 008
  * Leistedt et al. (2014) Leistedt, B., Peiris, H. V., & Roth, N. 2014, Physical Review Letters, 113
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, The Astrophysical Journal, 538, 473
  * LoVerde et al. (2008) LoVerde, M., Miller, A., Shandera, S., & Verde, L. 2008, Journal of Cosmology and Astroparticle Physics, 2008, 014
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, LSST Science Book, Version 2.0
  * Lucie-Smith et al. (2023) Lucie-Smith, L., Barreira, A., & Schmidt, F. 2023, Halo assembly bias from a deep learning model of halo formation
  * Lyth et al. (2003) Lyth, D. H., Ungarelli, C., & Wands, D. 2003, Phys. Rev. D, 67, 023503
  * Maion et al. (2022) Maion, F., Angulo, R. E., & Zennaro, M. 2022, Journal of Cosmology and Astroparticle Physics, 2022, 036
  * Maldacena (2003) Maldacena, J. 2003, Journal of High Energy Physics, 2003, 013
  * Manera et al. (2012) Manera, M., Scoccimarro, R., Percival, W. J., et al. 2012, Monthly Notices of the Royal Astronomical Society, 428, 1036
  * Matarrese & Verde (2008) Matarrese, S. & Verde, L. 2008, The Astrophysical Journal, 677, L77
  * Matarrese et al. (2000) Matarrese, S., Verde, L., & Jimenez, R. 2000, The Astrophysical Journal, 541, 10
  * Mueller et al. (2022) Mueller, E.-M., Rezaie, M., Percival, W. J., et al. 2022, Monthly Notices of the Royal Astronomical Society, 514, 3396
  * Pajer et al. (2013) Pajer, E., Schmidt, F., & Zaldarriaga, M. 2013, Physical Review D, 88
  * Pillepich et al. (2009) Pillepich, A., Porciani, C., & Hahn, O. 2009, Monthly Notices of the Royal Astronomical Society, 402, 191
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, Astronomy & Astrophysics, 594, A13
  * Planck Collaboration et al. (2019) Planck Collaboration, Akrami, Y., Arroja, F., et al. 2019, Planck 2018 results. IX. Constraints on primordial non-Gaussianity
  * Prada et al. (2023) Prada, F., Ereza, J., Smith, A., et al. 2023, The DESI One-Percent Survey: Modelling the clustering and halo occupation of all four DESI tracers with Uchuu
  * Rezaie et al. (2023) Rezaie, M., Ross, A. J., Seo, H.-J., et al. 2023, arXiv e-prints, arXiv:2307.01753
  * Riquelme et al. (2023) Riquelme, W. et al. 2023, Mon. Not. Roy. Astron. Soc., 523, 603
  * Rocher et al. (2023) Rocher, A., Ruhlmann-Kleider, V., Burtin, E., et al. 2023, The DESI One-Percent survey: exploring the Halo Occupation Distribution of Emission Line Galaxies with AbacusSummit simulations
  * Ross et al. (2017) Ross, A. J., Banik, N., Avila, S., et al. 2017, Monthly Notices of the Royal Astronomical Society, 472, 4456
  * Ross et al. (2012) Ross, A. J., Percival, W. J., Carnero, A., et al. 2012, Monthly Notices of the Royal Astronomical Society, 428, 1116
  * Rossi et al. (2020) Rossi, G., Choi, P. D., Moon, J., et al. 2020, Monthly Notices of the Royal Astronomical Society
  * Salopek & Bond (1990) Salopek, D. S. & Bond, J. R. 1990, Phys. Rev. D, 42, 3936
  * Sartoris et al. (2016) Sartoris, B., Biviano, A., Fedeli, C., et al. 2016, Monthly Notices of the Royal Astronomical Society, 459, 1764
  * Scoccimarro et al. (2012) Scoccimarro, R., Hui, L., Manera, M., & Chan, K. C. 2012, Physical Review D, 85
  * Slosar et al. (2008) Slosar, A., Hirata, C., Seljak, U., Ho, S., & Padmanabhan, N. 2008, Journal of Cosmology and Astroparticle Physics, 2008, 031
  * Springel (2005) Springel, V. 2005, Monthly Notices of the Royal Astronomical Society, 364, 1105
  * Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629
  * Tinker et al. (2010) Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, The Astrophysical Journal, 724, 878
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Naess, S., Genel, S., et al. 2018, The Astrophysical Journal, 867, 137
  * Wagner & Verde (2012) Wagner, C. & Verde, L. 2012, Journal of Cosmology and Astroparticle Physics, 2012, 002
  * Yamauchi et al. (2014) Yamauchi, D., Takahashi, K., & Oguri, M. 2014, Physical Review D, 90
  * Yu et al. (2023) Yu, J., Zhao, C., Gonzalez-Perez, V., et al. 2023, The DESI One-Percent Survey: Exploring A Generalized SHAM for Multiple Tracers with the UNIT Simulation
  * Yuan et al. (2023) Yuan, S., Zhang, H., Ross, A. J., et al. 2023, The DESI One-Percent Survey: Exploring the Halo Occupation Distribution of Luminous Red Galaxies and Quasi-Stellar Objects with AbacusSummit
  * Zhao et al. (2021) Zhao, C., Chuang, C.-H., Bautista, J., et al. 2021, Monthly Notices of the Royal Astronomical Society, 503, 1149