© The Authors 2024

1. Introduction

Primordial non-Gaussianities (PNG) provide a unique window into the Universe’s dynamics during the inflationary era. One of the most widely studied models is that of local primordial non-Gaussianities, where the amplitude of the deviation from Gaussianity is typically parameterised using $f_{\mathrm{NL}}^{\mathrm{local}}$ (Komatsu & Spergel 2001). In this paper, we only consider the local PNG, so we refer to this quantity simply as f_NL hereafter. Slow-roll single-field inflation predicts a nearly vanishing f_NL (Maldacena 2003; Creminelli & Zaldarriaga 2004). Therefore, detecting a significantly larger value (f_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). Currently, the tightest constraints on f_NL come from the bispectrum measured for the anisotropies of the cosmic microwave background, with f_NL = −0.9 ± 5.1 (68% c.l.) (Planck Collaboration IX 2020). However, this measurement is fundamentally limited by cosmic variance and it is not expected to reach the benchmark of σ(f_NL) = 1, even with the addition of polarisation data from the CMB (Baumann et al. 2009). Such a small error would allow us to detect or rule out most multi-field inflation models.

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 & Percival 2014; Ho et al. 2015; Leistedt et al. 2014; Castorina et al. 2019; Cabass et al. 2022; Rezaie et al. 2024). The most precise measurement to date from spectroscopic surveys comes from SDSS-IV/eBOSS, which measured f_NL = −12 ± 21 (68% c.l.) (Mueller et al. 2022). Current spectroscopic galaxy surveys, such as DESI¹ and Euclid², are expected to reach a level of precision of σ(f_NL)∼5 (Laureijs et al. 2011; Giannantonio et al. 2012; Sartoris et al. 2016; DESI Collaboration 2016). The next generation of cosmological surveys is expected to reach σ(f_NL) < 1 by combining the information coming from galaxy clustering with other probes, such as HI intensity mapping (LSST Science Collaboration 2009; Yamauchi et al. 2014; Bacon et al. 2020; Kopana et al. 2024).

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 controlled environment (e.g., Ross et al. 2017; Avila et al. 2021; Chan et al. 2022; Riquelme et al. 2023). To date, there has been 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 (σ(f_NL)∼5), the required simulation would need to have > 16 000³ particles in a volume larger than (4 h⁻¹ Gpc)³. As this is extremely expensive to run, several techniques have been developed in recent years to increase the statistical significance of simulations. This increase can be obtained either by changing the way the initial conditions are generated (for example the fixed-and-paired technique proposed by Angulo & Pontzen 2016) or by using a set of approximate methods that are computationally less demanding and/or less expensive (Chartier et al. 2021; Ding et al. 2022; Kokron et al. 2022; DeRose et al. 2023).

First, we needed to clarify whether the fixed-and-paired technique could be used in the presence of local PNG, as fixing the Fourier modes’ amplitudes breaks Gaussianity of the initial conditions. However, in Avila & Adame (2023), we demonstrated that this can also be used in the presence of local PNG to have a larger effective volume for simulations. In addition, in that work, we introduced a technique further to suppress the noise in the f_NL measurement. This suppression is achieved by matching the phases between two simulations with and without PNG, which removes much of the noise in the measured statistics associated with the cosmic variance. We use the results from Avila & Adame (2023) as a crucial ingredient in the following analysis.

This work presents the PNG-UNITSIMS suite, which includes the largest simulation (in terms of number of particles) to date with non-Gaussian initial conditions, the PNG-UNIT. This simulation evolves 4096³ dark matter particles within a cosmological volume of (1 h⁻¹ Gpc)³, resulting in a mass resolution of 1.24 × 10⁹ h⁻¹ M_⊙. Most of the previously developed simulations with local-PNG were run with fewer particles, that is, by over a factor of 18 (e.g., Nishimichi et al. 2010; Baldauf et al. 2016; Biagetti et al. 2017; Coulton et al. 2023). Due to this low resolution, it was not possible to use these simulations to study the scale-dependent bias for low-mass halos and galaxies (i.e., M_halo ≲ 1 × 10¹³ h⁻¹ M_⊙). After the submission of this paper, Hadzhiyska et al. (2024) presented a new PNG simulations with the same number of particles as the PNG-UNIT. However, these simulations have a volume of (2 h⁻¹ Gpc)³, resulting in a mass resolution ×8 times lower than the PNG-UNIT. Therefore, their objectives are different and complementary to the project presented in this paper.

The resolution of the PNG-UNIT simulation is enough to finely resolve (with more than ∼100 particles) the halos hosting the emission-line galaxies (ELGs) targeted by DESI and Euclid ( ∼ 10¹¹ h⁻¹ M_⊙) (Gonzalez-Perez et al. 2018; Yuan et al. 2024; Rocher et al. 2023; Reyes-Peraza et al. 2024). Other DESI tracers (BGS, LRG, QSO) are expected to probe even larger masses, which are also well resolved by this simulation (Yu et al. 2024; Rocher et al. 2023; Prada et al. 2023). The initial conditions of the PNG-UNIT are generated with amplitudes fixed to their expectation values and phases matched to one of the original UNIT simulations with Gaussian initial conditions (Chuang et al. 2019).

We first used the PNG-UNITSIMS suite to constrain the non-Gaussian bias parameter, b_ϕ, for mass-selected halos. This parameter controls how the clustering of halos, or any given tracer of the matter density field (cosmological tracer hereafter), responds to the presence of PNGs. However, b_ϕ is completely degenerate with the amplitude of the non-Gaussianities measured by f_NL (Barreira 2022a). It is, therefore, crucial to have a good understanding of b_ϕ if we aim to measure f_NL from data. Thus, understanding b_ϕ from simulations will be essential.

An assumption is made for the mass function for halos where the mass only depends on the height of the density peak, namely, ν = δ_c/σ(M), where δ_c is the critical overdensity and σ(M)² is the variance of the density field smoothed at mass M. On this basis, a theoretical prediction for b_ϕ can be made as a function of the linear bias b₁ (Dalal et al. 2008):

$$\begin{array}{c}\hfill b_{\varphi }=2{\delta }_{\mathrm{c}}(b_1-1).\end{array}$$(1)

This result is what is usually called the universality relation. Nevertheless, several works have shown that b_ϕ 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 demonstrated that b_ϕ may differ depending on the particular cosmological tracer (Barreira et al. 2020; Barreira 2022b).

In the literature, b_ϕ has been constrained using not only the scale-dependent bias in simulations with PNG but also with other methods, namely the separate universe technique (Biagetti et al. 2017; Barreira et al. 2020; Barreira 2022a; Lazeyras et al. 2023). Although b_ϕ can be measured very precisely with this approach, these simulations cannot validate non-Gaussian analysis tools of, for example, galaxy clustering, weak lensing statistics and cluster counts given they assume Gaussian initial conditions. Moreover, in this paper we find that the errors for the b_ϕ obtained with the separate universe method may be underestimated.

In this work, we divided the halos into 12 pseudo-logarithmically spaced mass bins and studied their clustering by measuring their power spectrum. This statistic is sensitive to the product b_ϕf_NL through the scale-dependent bias (Dalal et al. 2008; Slosar et al. 2008). Here, we estimate the variance of the power spectrum for halos in PNG-UNIT from a set of fast mocks generated with FASTPM (Feng et al. 2016) with fixed initial conditions. We further improved the precision measuring b_ϕf_NL by matching the phases of the initial conditions. In the future, we plan to carry out a study of how b_ϕ depends on 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 Sect. 2, focusing on the scale-dependent bias. In Sect. 3, we present the PNG-UNITSIMS suit of simulations developed for this study. In Sect. 4, we describe the methods employed for defining and relating the different mass bins, estimate the power spectrum variance, and apply the matching technique to improve the constraints on the PNG response parameters. In Sect. 5, we present and discuss our main findings, including the constraints on the b_ϕ parameter for mass-selected halos. In Sect. 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 Sect. 7. In Sect. 8, we discuss how our results can be used as a prior for a DESI-like survey and how it would affect the constraints on f_NL. Finally, we present our conclusions in Sect. 9.

2. Theory

Dark matter halos trace the matter distribution, as these objects are assumed to form from high-density peaks. However, the gravitational collapse of dark matter halos is a non-linear and complex process. On 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 halos and also other cosmological tracers, such as galaxies, can be written as a function of position and redshift:

$$\begin{array}{c}\hfill {\delta }_{\mathrm{h}}(\mathit{x},z)={\displaystyle \sum _O}{b}_{\mathrm{O}}(z)O(\mathit{x},z).\end{array}$$(2)

The operators, O, describe all the fields that may affect the matter density field, while b_O represents the bias parameters associated with each operator. 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 ∝ ∂_i∂_jϕ (see Desjacques et al. 2018, for an extensive review on the bias formalism). All the terms proportional to the matter density contrast δ_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 (∇²ϕ ∝ δ_m). The leading order of Eq. (2) can be expressed as:

$$\begin{array}{c}\hfill {\delta }_h=b_1{\delta }_m,\end{array}$$(3)

where δ_h is the density field of halos, δ_m is the density field of the matter distribution, and b₁ is the linear bias. This expression is only valid on the purely linear regime, which typically is considered those below $k_{\max }\lesssim 0.1h{\mathrm{Mpc}}^{-1}$. By including more terms in Eq. (2), the clustering of the halos can be modelled even in the mildly non-linear regime. This modelling is unnecessary as the effect of local PNG is mostly encoded in the largest scales (see Sect. 2.2). Hence, most of our analysis is focussed on the purely linear regime.

Using the linear bias, b₁, from Eq. (3), the power spectrum of halos can be expressed in terms of the linear matter power spectrum, P_m, m(k, z), as:

$$\begin{array}{c}\hfill P_{\mathrm{h},\mathrm{h}}(k,z)=⟨{\delta }_{\mathrm{h}}{\delta }_{\mathrm{h}}⟩=b_1^2(z)P_{\mathrm{m},\mathrm{m}}(k,z).\end{array}$$(4)

Below, we describe how PNG may arise from non-standard inflationary models (Sect. 2.1) and how it induces a scale dependence in the linear galaxy and halo bias (Sect. 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 assumptions, which may lead to large primordial non-Gaussianities (PNG) (e.g., 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_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, this relation is modified in the presence of PNGs, and non-trivial contributions to the higher-order statistics may arise.

At first order, this f_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 ϕ = ϕ_G + ϕ_NG, it can be shown that the leading order contribution to the bispectrum is given by:

$$\begin{array}{c}\hfill ⟨{\varphi }_{\mathrm{G}}{\varphi }_{\mathrm{G}}{\varphi }_{\mathrm{NG}}⟩={(2\pi )}^3{\delta }_{\mathrm{D}}({\mathit{k}}_1+{\mathit{k}}_2+{\mathit{k}}_3){\mathcal{B}}_{\varphi }({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3)\end{array}$$(5)

where δ_D is the Dirac delta in 3D and ℬ_ϕ(k₁, k₂, k₃) is the primordial bispectrum. Depending on the shape induced in the primordial bispectrum, three types of PNGs have been traditionally considered: local, equilateral and orthogonal. Nevertheless, the landscape of possible non-Gaussianities is much richer and more complex.

This work focuses only on the local type, 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):

$$\begin{array}{c}\hfill \varphi (\mathit{x})={\varphi }_{\mathrm{G}}(\mathit{x})+f_{\mathrm{NL}}({\varphi }_{\mathrm{G}}{(\mathit{x})}^2-⟨{\varphi }_{\mathrm{G}}{(\mathit{x})}^2⟩),\end{array}$$(6)

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

One finds that a non-vanishing component to the three-point function is induced from the expression for the local PNGs given in Eq. (6). In this case, the bispectrum peaks in the squeezed limit (k₃ ≪ k₁, k₂) and is given by

$$\begin{array}{c}\hfill {\mathcal{B}}_{\varphi }({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3)=2f_{\mathrm{NL}}{P}_{\varphi }(k_1)P_{\varphi }(k_2)+(2\mathrm{cyc}.\mathrm{perm}),\end{array}$$(7)

where P_ϕ(k) is the power spectrum of the primordial potential ϕ, and (2 cyc. perm) are two cyclic permutations of the term P_ϕ(k₁)P_ϕ(k₂).

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 via:

$$\begin{array}{c}\hfill \delta (k,z)=\alpha (k,z)\varphi (k,z),\end{array}$$(8)

where α(k, z) is defined as

$$\begin{array}{c}\hfill \alpha (k,z)=\frac{2D(z)}{3{\mathrm{\Omega }}_{\mathrm{m},0}}\frac{c^2}{H_0^2}\frac{g(0)}{g(z_{\mathrm{rad}})}{k}^{2}T(k).\end{array}$$(9)

The parameter Ω_m,  0 is the matter density parameter today, c is the speed of light, and H₀ is the Hubble parameter today. The transfer function T(k) is normalised such that T(k → 0) = 1. The growth factor D(z) is also normalised to D(z = 0) = 1. The term $\frac{g(z_{\mathrm{rad}})}{g(0)}$ (with g(z) = (1 + z)D(z)) arises from the normalisation we assume for D(z) (Mueller et al. 2018). If the growth factor is normalised to (1 + z)⁻¹ during matter domination the term $\frac{g(z_{\mathrm{rad}})}{g(0)}$ can be omitted. For the cosmology described in Planck Collaboration XIII (2016) and we assume throughout the paper, we find that $\frac{g(z_{\mathrm{rad}})}{g(0)}=1.275$.

2.2. Scale-dependent bias

In the presence of local PNG, the leading-order term proportional to f_NLϕ in the bias expansion (Eq. (2)) gives

$$\begin{array}{c}\hfill {\delta }_{\mathrm{h}}(k,z)=b_1(z){\delta }_{\mathrm{m}}(k,z)+b_{\varphi }{f}_{\mathrm{NL}}\varphi =[b_1(z)+\frac{f_{\mathrm{NL}}{b}_{\varphi }(z)}{\alpha (k,z)}]{\delta }_{\mathrm{m}}(k,z).\end{array}$$(10)

Here b₁ is the linear bias, and α(k, z) is given in Eq. (9). b_ϕ is the bias parameter associated with the ϕ operator. This term only appears in the presence of PNG, as it is associated with the primordial bispectrum (Eq. (7)), which vanishes in the case of Gaussian initial conditions. The term α introduces the 1/k² dependence of δ_h on the scale k.

Using this new bias expansion, the power spectrum for dark matter halos, or any other biased tracer of the total matter density, can be computed similarly to Eq. (4). The resulting halo power spectrum is given by

$$\begin{array}{c}\hfill P_{\mathrm{h},\mathrm{h}}(k,z)={[b_1(z)+\frac{f_{\mathrm{NL}}{b}_{\varphi }(z)}{\alpha (k,z)}]}^2{P}_{\mathrm{m},\mathrm{m}}(k,z).\end{array}$$(11)

The relationship between P_h, h(k, z) and P_m, m(k, z) now depends on the scale due to the term α(k, z), where the leading-order effect is proportional to 1/k². This proportionality is what is commonly called the ‘scale-dependent bias’.

However, different values of f_NL may produce the same signal depending on the value of the bias parameter b_ϕ as they are perfectly degenerate. This issue has already been discussed in the literature (Barreira 2022a).

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

$$\begin{array}{c}\hfill b_{\varphi }=\frac{1}{{\overline{n}}_{\mathrm{h}}}\frac{\partial {\overline{n}}_{\mathrm{h}}}{\partial log(f_{\mathrm{NL}}\varphi )}=\frac{1}{{\overline{n}}_{\mathrm{h}}}\frac{\partial {\overline{n}}_{\mathrm{h}}}{\partial log{\mathcal{A}}_{\mathrm{s}}},\end{array}$$(12)

where ${\overline{n}}_{\mathrm{h}}$ is the mean number of halos and 𝒜_s is the amplitude of the primordial scalar perturbations. Then, if a universal mass function is assumed, one can derive a theoretical prediction for the value of b_ϕ (Dalal et al. 2008; Slosar et al. 2008; Desjacques et al. 2018):

$$\begin{array}{c}\hfill b_{\varphi }=2{\delta }_{\mathrm{c}}(b_1-p).\end{array}$$(13)

Here, δ_c = 1.686 is the density threshold for the spherical collapse, and p may take different values depending on the population of objects under consideration. By assuming a universal mass function, then p = 1 for all tracers of the matter distribution (Dalal et al. 2008; Slosar et al. 2008). This expression with p = 1 is usually called the universality relation for b_ϕ (noted as the universality relation hereafter).

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 halos that have recently undergone a major merger process, b_ϕ is better described by 2δ_c(b₁ − 1.6), while Barreira et al. (2020) showed that for galaxies selected by their stellar mass, the relation with b₁ is better described by 2δ_c(b₁ − 0.5). Moreover, b_ϕ may depend not only on the halo mass (through b₁) but also on other properties such as the halo concentration (Lazeyras et al. 2023; Fondi et al. 2024).

In Fig. 1, we present the expected power spectrum of dark matter halos at linear order in perturbation theory for different values of the non-Gaussian parameter, f_NL. We chose the value of b₁ to be similar to the one expected for the emission-line galaxies (ELGs), a tracer sample targeted by the current galaxy surveys (Yu et al. 2024). 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₁, as given by Eq. (4). At large scales, the term proportional to b_ϕf_NL, which appears only for the case with local PNG, induces an enhancement for b_ϕf_NL >  0 (red and green lines), and suppression if b_ϕf_NL < 0 (purple line) of the power spectrum compared with the Gaussian case. When going to even larger scales, the term proportional to (b_ϕf_NL)² dominates, and an enhancement of the power spectrum is shown for both cases, b_ϕf_NL >  0 and b_ϕf_NL < 0.

Fig. 1. Theoretical power spectrum obtained from Eq. (11). The blue line shows the linear matter-matter power spectrum for the UNIT cosmology at z = 1. The b1 parameter is the same for all cases, and we chose it to be similar to the one expected for the DESI emission-line galaxies (ELGs) at that redshift (DESI Collaboration 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 fNL = 100 and a PNG-response parameter bϕ given by Eq. (13) with p = 1. The red line displays how a different value of fNL can lead to the same signal if the bϕ parameter does not follow the universality relation. Red and green lines are on top of each other since the product bϕfNL is the same. Finally, the purple line shows the kind of signal we expect in case bϕfNL < 0.

The scale-dependent bias is only sensitive to the product b_ϕf_NL. Because of this, f_NL cannot be extracted from this unless we have prior knowledge of b_ϕ (Barreira 2022a), as illustrated in Fig. 1, where the green and red lines show the same signal although the f_NL value is different.

3. Simulations

The UNITSIMS³ 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 an effective volume equivalent to or larger than these surveys and a mass resolution high enough to resolve the dark matter halos hosting the targeted galaxies and tracers.

Most state-of-the-art simulations, including the UNIT simulation (Sect. 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 (Sect. 3.2), together with two lower-resolution full N-body simulations (LR-UNIT, Sect. 3.2.2) and a set of 100 fast N-body FASTPM mocks (Sect. 3.3). We developed another set of smaller simulations to perform further robustness tests on our analysis, which we discuss in Sect. 3.4.

All simulations presented in this paper share the same cosmology (Planck Collaboration XIII 2016), as we summarise in Table 1. We outline the characteristics of the simulations in Table 2. Moreover, we focus our analysis on the snapshots at redshift z = 1.032.

Throughout the paper, we refer to the suite as PNG-UNITSIMS/UNITSIMS and to the 4096³ particle simulation as PNG-UNIT/UNIT.

3.1. The UNITSIMS

The UNITSIMS are a set of two pairs of state-of-the-art, dark matter-only, full N-body simulations with fixed-and-paired initial conditions (Angulo & Pontzen 2016; Chuang et al. 2019). Each of them tracks the evolution of 4096³ particles in a periodic box with a side length of L = 1 h⁻¹ Gpc from redshift z = 99 down to z = 0. The mass resolution is m_p = 1.24 × 10⁹ h⁻¹ M_⊙, so the least massive halos 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 also using in this work one of the original UNITSIMS, for reference, the one labelled FixAmpInvPhase_001.

The initial conditions were set using the second-order Lagrangian Perturbation Theory (2LPT), implemented in FASTPM⁴ (Feng et al. 2016). The required initial power spectrum was computed using the Boltzmann solver CAMB⁵ (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 long-range and short-range components. We employed the non-public version of the code L-GADGET-2 as it is highly optimised regarding 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 ϵ = 6 kpc h⁻¹ and then the particles were evolved from z = 99 down to z = 0, stored in a total of 129 snapshots logarithmically spaced in the scale factor.

Dark matter halos were identified using the publicly available code ROCKSTAR⁶ (Behroozi et al. 2012). This code uses the particles’ 6D phase space information to identify the dark matter halos. These structures are identified as having at least 20 particles per halo. We used CONSISTENTTREES to track their merger histories and build the merger trees. More details regarding these simulations can be found in Chuang et al. (2019).

3.2. The PNG-UNITSIMS

We developed the new PNG-UNITSIMS suite for this paper with the same features as the original UNITSIM simulations but with non-Gaussian initial conditions. Specifically we set f_NL = 100. The suite consists of a high-resolution simulation (Sect. 3.2.1), a lower resolution version of it (Sect. 3.2.2) and a set of 100 fast mocks generated with FASTPM, half of them with f_NL = 0 and the other half with f_NL = 100 (Sect. 3.3).

To generate the initial conditions of all of these simulations, we have used the FASTPM code. In the first place, the initial conditions are generated for ϕ 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 Eq. (8) and then the transformation given in Eq. (6) 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 simulation in terms of number of particles with primordial non-Gaussianites to date. It is a dark matter-only, full N-body simulation with local primordial non-Gaussianities. As the original UNIT simulation, it tracks the evolution of 4096³ particles in a periodic box with a side length of L = 1 h⁻¹ Gpc from redshift z = 99 down 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_NL used for generating the initial conditions is set to f_NL = 100. This large value of f_NL is already ruled out by Planck at ∼20σ (Planck Collaboration IX 2020). However, we wanted a strong signal from the scale-dependent bias to understand this effect better. With a more realistic f_NL, given our limitation in scales with minimum wavenumber k_min = k_f = 0.006 h Mpc⁻¹ (where k_f is the fundamental wavenumber), we would not be able to detect the PNG signal in the halo power spectrum.

In addition, the initial conditions of the PNG-UNIT were generated using the fixed-and-matched technique described in Avila & Adame (2023). The amplitudes were set to the expectation value as in Angulo & Pontzen (2016). This method 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 the same as in the reference Gaussian UNIT simulation. The N-body code and the halo finder were run using the same configuration as for the UNITSIMS. We focus our analysis on the snapshot at redshift z = 1.

3.2.2. The LR-UNIT simulations

In addition to the high-resolution simulation presented above, for the purposes of this paper, we have developed two new lower-resolution versions of the previous ones (with f_NL = 0 and f_NL = 100). The primary purpose of these LR-UNIT simulations is to understand the effects of mass resolution better (see Sect. 6.1). For this reason, we use the same configuration as the high-resolution simulations. The only difference is the resolution (see Table 2).

As we used the same code and seed to generate the initial conditions, all the long-wavelength modes were set to the same 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.

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_NL = 0) while the other set had non-Gaussian initial conditions (f_NL = 100).

FASTPM is a code that combines perturbation theory with a particle–mesh (PM) approach to quickly generate cosmological simulations at the cost of poor resolution of the short-wavelength modes. This code is particularly useful for generating many simulations quickly, which can be used, for example, for estimating covariance matrices at a low computational cost. In this work, we use FASTPM mocks to estimate the variance of the power spectrum of the halos in simulations with fixed initial conditions as well as for obtaining the correlation coefficients required by the matching technique as we discuss in Sect. 4.3.

Each simulation box contains 2048³ dark matter particles within a cosmological volume of 1 (h⁻¹ Gpc)³, the same volume as the UNIT simulation. The mass resolution is m_p = 9.97 × 10¹⁰ h⁻¹ M_⊙, which is 8 times coarser than the high-resolution simulations (same resolution as LR-UNIT simulations) but still sufficient for this work. The force resolution in the PM approach is given by the grid size used to compute the potential and the forces. We use a regular grid of 4096³ cells, giving a spatial resolution of 244.14 h⁻¹ kpc.

Following the same procedure as the full N-body simulations, the initial conditions were set at z = 99 using second-order Lagrangian perturbation theory with fixed amplitudes. The particles were then evolved to z = 1 using 100 linearly spaced timesteps.

We use the friends-of-friends (FoF) algorithm to generate the halo catalogues, 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-UNIT catalogues leads to a different mass definition for dark matter halos. This issue is addressed in Sect. 4.1.

3.4. Supporting simulations

In addition to the PNG-UNIT and the FASTPM mocks presented above, we developed an additional set of simulations mainly aimed at checking the robustness of our results. We want to check if the results obtained from fitting b_ϕ in Eq. (11) are consistent with measuring b_ϕ using the separate universe technique (see Sect. 6.2). Hence, the set of simulations is split into two groups: high-resolution (HR) and low-resolution (LR) separate universe simulations.

Each group consists of 3 simulations with 512³ dark matter particles within a periodic cubic box of size L_box,  HR = 67.5 h⁻¹ Mpc and L_box,  LR = 250 h⁻¹ Mpc for the high-resolution and low-resolution groups respectively. This configuration leads to a mass resolution which is eight times better and eight times worse, respectively, compared to the baseline UNITSIMS. Again, the fiducial cosmology is the same as the one of the UNIT simulations, summarised in Table 1. Following a similar procedure as shown in Barreira (2022a), we produced three simulations with the same initial conditions but varied the 𝒜_s (or equivalently σ₈) by 5%. This variation allows us to numerically measure b_ϕ by computing the derivative (Eq. (22)), as we discuss in Sect. 6.2.

For all these simulations, we evolved the particles using the L-GADGET 2 code, as in the case of the high-resolution simulations and we set 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 halos from the PNG-UNIT and LR-UNIT simulations. The mass function has been calculated from 2.48 × 10¹⁰ h⁻¹ M_⊙ (20 times the particle mass in PNG-UNIT) up to 1 ×  10¹⁵ h⁻¹ M_⊙ in 50 logarithmically spaced bins. We find that the low-resolution simulations converge to the high-resolution ones at < 5% for $M\sim 100m_{\mathrm{part}}^{\mathrm{LR}}$, with $m_{\mathrm{part}}^{\mathrm{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 halos is unreliable in the low-resolution simulations.

Fig. 2. Halo mass functions for the simulations described in Section 3 (top). Ratio of halo mass functions relative to the Gaussian UNIT simulation as the reference (bottom). Shaded areas represent the $M=20m_{\mathrm{part}}^{\mathrm{LR}}$ mass threshold for the low-resolution simulation. The minimum mass shown here represents the 20 particle limit for the higher resolution simulations.

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

We checked that the differences in the mass function are mainly due to the different halo-finding algorithms by running Rockstar on one of FASTPM mocks. Running both halo finders on the same FASTPM mock, we find that for high masses, the mass function matches that of the full N-body simulations, while we find a deficit in the number of low-mass halos; this is probably due to differences in the resolution of the forces between the two approaches.

Moreover, we find an increment of the high-mass halos in the f_NL = 100 simulations with respect to the Gaussian ones. This increment was expected due to the change of the probability density function (PDF) of δ 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 halos 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 scales between k_nyquist/4 and k_nyquist/2 of the order of 10%. This results in an overproduction of halos of the same order of ∼10% and explains the differences between the UNIT and the PNG-UNIT simulations at the mass function level. However, this does not significantly affect the halo clustering measurements compared to the low-resolution simulations, as we explain in Sect. 6.1.

4. Uncertainty of the Power Spectrum in simulations with fixed-and-matched initial conditions

In this work, we study the PNG parameters b_ϕ and f_NL as a function of halo mass. We measure them using the scale-dependent bias induced by PNG in the halo power spectrum (Sect. 2.2). The dark matter auto-power spectra of the simulations are obtained by interpolating the particles onto a regular grid with N = 512³ cells using a cloud-in-cell (CIC) scheme (Hockney & Eastwood 1981) with a k-bin width equivalent to the fundamental wavenumber of the box (k_f = 2π/L_box = 0.00628 h Mpc⁻¹). We have used the PYTHON library NBODYKIT⁷ (Hand et al. 2018) for this task. The halo power spectrum is computed similarly 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 allows us to consider the variance suppression due to the initial conditions (see Sect. 4.2). However, fundamental differences exist between the fast mocks and the PNG-UNIT. The main differences are how the halos are defined and the mass resolution of the simulations. We now discuss how to connect the halos (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 halos, which are not substructures of a larger halo. For simulations where we run ROCKSTAR (see Table 2), the halo mass is estimated from a spherical overdensity up to the threshold of Δ = 200ρ_c, where ρ_c is the critical density (Behroozi et al. 2012). For the FASTPM mocks, the halo mass is the sum of the particle mass of its FoF components. This mass definition 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 differ significantly from those for the fast mocks. To relate halos between these two approaches, we aim to reproduce simultaneously the clustering and the abundances of the UNIT halos with the FASTPM mocks. To this end, for each mass bin of the UNITSIMS, we consider a maximum mass for the halos in the FASTPM mocks and select the following N most massive halos 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 Mpc⁻¹. We take the minimum and the maximum mass from the resulting bins and apply these cuts for the f_NL = 100 simulations.

Figure 3 shows the power spectrum of FASTPM compared to that from the PNG-UNIT. The lines are the power spectrum computed for the PNG-UNIT 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 above procedure.

Fig. 3. Solid lines represent the power spectrum of PNG-UNIT halos within the different halo mass bins (see Table 3 for the mass ranges). Shaded regions show the standard deviation around the mean power spectrum of the FASTPM mocks, generated using halos in mass bins that reproduce the clustering and the abundance of the N-body simulation (see Sect. 4.1). There are no halos in the two lowest-mass bins (1 and 2, shown in dark blues) in the FASTPM mocks since the mass resolution is not high enough to resolve them. The top panel shows the fNL = 0 case, and the bottom panel shows the fNL = 100 one.

There are no halos in the fast mocks for the three lowest-mass bins since the mass resolution is not high enough to resolve them. In the case of the two most massive bins, although we can recover the clustering, we cannot simultaneously recover the abundance. This problem is associated with the mass and bias function differences due to the distinct halo finders used (see Fig. 2). In addition, the power spectra for the three most massive bins in the PNG-UNIT are very noisy since the number density is very low (n < 10⁻⁵ h³ Mpc⁻³).

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):

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{S}}^2(k)=\frac{4{\pi }^2}{V{k}^2\mathrm{\Delta }k}{(P(k)+\frac{1}{\overline{n}})}^2,\end{array}$$(14)

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

As the initial conditions of the simulations were fixed to their expectation value (Sect. 3), a considerable suppression is expected with respect to Eq. (14), 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 σ(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 a problem obtaining this quantity for bins beyond the resolution limit of the FASTPM mocks.

Therefore, we propose fitting the ratio σ(P(k))/P(k) with a smooth function in order to get a less noisy estimate of the errors. For that, we measure σ(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-UNIT) and fit it to the following ansatz that we find convenient after several trials of different functional forms:

$$\begin{array}{c}\hfill \frac{\sigma ((P(k))}{P(k)}=f(k;a,b,c)=a(1+\frac{1}{1+bexp(-ck)}).\end{array}$$(15)

In Fig. 4, we show the fits of the FASTPM mocks to the previous expression in solid lines. To extrapolate σ(P(k))/P(k) beyond the mass resolution of the fast mocks, we fitted the three parameters a, b and c from Eq. (15) we obtained as a function of the mass. For parameter a, we use an exponential function; for parameters b and c, we use a polynomial of order two as a function of log₁₀M_halo. The extrapolations are shown in Fig. 4 as dashed lines.

Fig. 4. Variance over mean power spectrum for each halo mass bin (indicated by colour with the mass ranges defined in Table 3). Markers with error bars show the data taken directly from the FASTPM mocks, while solid lines show results obtained by fitting these curves to Eq. (15). The extrapolation to each mass bin as described in Sect. 4.2 is shown in dashed lines.

We will now adopt by default the errors coming from the fit to Eq. (15) for bins 3–10 and the extrapolation for bins 0–2 and 11–12. Nevertheless, we checked the impact of this assumption, as detailed in Sect. 7.1.

4.3. Reducing the errors using fixed-and-matched simulations

When two simulations share the same set of random phases in their initial conditions, their cosmic variance is highly correlated even if they have different cosmologies. Here, we discuss how we take advantage of the correlated noise from two simulations with different f_NL to constrain the PNG bias parameters with high precision. This noise reduction technique is what we call ‘matching the simulations’.

The initial conditions for an N-body simulation are a realisation of a Gaussian random field generated from the power spectrum of primordial perturbations. There are two degrees of freedom: the amplitude of the perturbations, whose expected value is determined by the power spectrum, and the phase of those perturbations, which are randomly uniformly distributed between 0 and 2π. This stochasticity in the generation of the initial conditions (ICs), usually called the cosmic variance, leads to some degree of noise when measuring quantities from simulations in the same way as we expect it to be in the real Universe. The level 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, we could cancel out much of the noise affecting measurements by matching that simulation with an existing one with the same ICs but a different cosmology.

We applied this in Avila & Adame (2023) for measuring f_NL, where only three ingredients were needed for the measurement: the value of f_NL obtained from the Gaussian simulation, the same quantity from the non-Gaussian simulation and the correlation coefficient ρ between both simulations. In that study, we used a set of full N-body simulations with matched ICs (i.e., the same stochastic part) to measure ρ.

In this study, we adopt an alternative strategy by leveraging only two simulations, UNIT and PNG-UNIT, along with a suite of computationally efficient mocks to estimate ρ. This approach circumvents the need for an extensive series of full N-body mocks, which would be impractically demanding on computational resources at this mass resolution. More precisely, we determine b_ϕf_NL for each bin within each FASTPM realisation and then calculate the Pearson correlation coefficient to obtain ρ.

We find that the measurements of b_ϕf_NL correlate more when the mass of the halos is lower. This finding may be because the lower the mean mass of the bins, the higher the number of halos included, and therefore, the less noisy the power spectrum measurement. We performed a linear fit of this ρ as a function of the mass to reduce the noise associated with the measurement of this quantity⁸. The results of the linear fit are reported in Table 3. However, studying the evolution of the correlation with mass in depth is beyond the scope of this paper.

Thus, we are left with the question of obtaining the correlation coefficients for those bins where we do not have FASTPM mocks. One possibility is to use the values obtained from the linear regression as a function of bin mass. This regression would give us a correlation coefficient of ρ ∼ 0.9 for the less massive halo bins (i.e., bins 0–2). While this does not affect the measurement of b_ϕ, it is critical for the error bars. However, we have no guarantee that this relation of ρ with the mass would saturate at some point. Therefore, we have chosen a more conservative approach (as it results in larger reported uncertainties) and assume for the bins 0–2 the average ρ that we have calculated from the FASTPM mocks (i.e., bins 3–10), giving for those bins ρ = 0.703.

In addition, we propose to apply this to constrain directly the product b_ϕf_NL (instead of just f_NL as done in Avila & Adame 2023), which is the quantity we are sensitive to. To this end, we construct this new estimator for b_ϕf_NL, which takes into account the matching between both simulations:

$$\begin{array}{c}\hfill \mathrm{\Delta }{b}_{\varphi }{f}_{\mathrm{NL}}={[b_{\varphi }{f}_{\mathrm{NL}}]}_{100}-{[b_{\varphi }{f}_{\mathrm{NL}}]}_0,\end{array}$$(16)

where [b_ϕf_NL]₁₀₀ is the quantity measured in the non-Gaussian simulation with f_NL = 100 and [b_ϕf_NL]₀ corresponds to the Gaussian one. Given that the noise in [b_ϕf_NL]₁₀₀ and [b_ϕf_NL]₀ is highly correlated, it is expected that the average Δb_ϕf_NL would converge to the mean much faster than the average of [b_ϕf_NL]₁₀₀ or [b_ϕf_NL]₀ alone, since we are cancelling some part of the noise. The variance of Δb_ϕf_NL is given by

$$\begin{array}{cc}\hfill {\sigma }^2(\mathrm{\Delta }{b}_{\varphi }{f}_{\mathrm{NL}})& ={\sigma }^2({[b_{\varphi }{f}_{\mathrm{NL}}]}_{100})+{\sigma }^2({[b_{\varphi }{f}_{\mathrm{NL}}]}_0)\hfill \\ \hfill & -2\rho \sigma ({[b_{\varphi }{f}_{\mathrm{NL}}]}_{100})\sigma ({[b_{\varphi }{f}_{\mathrm{NL}}]}_0),\hfill \end{array}$$(17)

where σ([b_ϕf_NL]₁₀₀) is the error of the quantity measured in the non-Gaussian simulation and σ([b_ϕf_NL]₀) corresponds to the Gaussian one, and ρ is the correlation coefficient between the two measured using the FASTPM mocks.

While, in principle, this method can be applied to any two values of f_NL, the advantage of choosing one value to be f_NL = 0 is that ⟨Δb_ϕf_NL⟩=⟨[b_ϕf_NL]₁₀₀⟩. From now on, when we apply this matching, we label the errors as σ_S + M if we use Eq. (14) for the variance of the power spectrum and σ_F + M if we use the method described in Sect. 4.2.

4.4. Variance estimation methods

In the previous sections, we discuss how to reduce the noise in the measurements of b_ϕf_NL by taking advantage of the fixed-and-matched initial conditions of the PNG-UNIT. We compare our results with those derived when no noise reduction techniques are employed to ensure that we are not introducing any bias when using these techniques.

We defined the following errors for the power spectrum, σ(P(k)):

• Standard errors, σ_S, come from Eq. (14) with only one simulation used to measure b_ϕf_NL. This error is not correlated with the other f_NL simulation. This method is expected to provide the largest estimate of the error.

• Fixed errors, σ_F, are obtained from the best-fit σ(P(k)) (Eq. (15)) measured from the FASTPM mocks for bins 3–10 in halo mass. For bins 0–2 and 10–12, we do not recover the clustering nor the abundances with FASTPM mocks simultaneously, so we need to perform an extrapolation for those mass bins (Sect. 4.3). To obtain these errors, we do not correlate f_NL = 0 with f_NL = 100.

In addition, we have errors that take into account the noise reduction due to the fixed-and-matched initial conditions (Eq. (17)):

• ‘standard + matched’ error, σ_S + M, where the standard error σ_S for the power spectrum from the simulations with f_NL = 0 is correlated with the ones with f_NL = 100 using the method described in Sect. 4.3;

• ‘fixed + matched’ errors, σ_F + M, where we assume the σ_F errors for the power spectrum and correlate matched simulations with different f_NL. We consider these errors for the baseline analysis.

In Sect. 7.2, we compare how these assumptions affect the constraints on b_ϕ and p.

5. Constraining the PNG bias parameters for mass-selected halos

Given the input value of f_NL used to generate the initial conditions for the PNG-UNIT, we can measure directly b_ϕ from the scale-dependent bias (Eq. (11)). In this section, we discuss how we constrain b_ϕ for mass-selected halos using the suppressed variance from the matching technique (Sect. 4.3).

First, we check how well Eq. (11) describes the clustering of the dark matter halos in the PNG-UNIT under the standard assumption that b_ϕ is given by the universality relation (i.e., p = 1). We measure the non-Gaussian signal b_ϕf_NL in our simulations (Sect. 5.2), and by assuming the universality relation, we report our results for f_NL (Sect. 5.3). Finally, using the input values of f_NL, we directly measure p (Sect. 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 parameter estimation common to all our analyses. We study the impact of varying these assumptions in Sect. 7.

To obtain the PNG parameters, we performed Markov chain Monte Carlo (MCMC) using the Python package EMCEE⁹ (Foreman-Mackey et al. 2013). All the cosmological parameters are fixed except b₁ and b_ϕf_NL for which we set flat priors over a wide range ([0, 10] for b₁ and [ − 10 000, 10 000] for b_ϕf_NL). We assume a Gaussian likelihood with χ² given by

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=0}^{N_k}}\frac{{(P_{\mathrm{halos}}(k_i)-P_{\mathrm{model}}(k_i;b_1;b_{\varphi }{f}_{\mathrm{NL}}))}^2}{{\sigma }^2(P(k_i))}.\end{array}$$(18)

Here, P_model(k) is obtained using Eq. (11) by inputting the matter power spectrum P_m, m(k) computed directly from the simulation. In the baseline analysis, we use the fixed+matched errors σ_F + M defined in Sect. 4.4 for σ(P(k)). Moreover, in the subsequent figures, we use solid circles when we are using the best fit to Eq. (15) and the mocks, and we use empty circles when we performed extrapolation (i.e. for bins 0–2). We then obtain σ(b_ϕf_NL) from the posterior distribution after marginalising over b₁. The model described by Eq. (11) is only accurate in the purely linear regime.

Due to the limitation of our model to the linear regime, we need to drop the modes associated with higher-order nonlinear terms in the power spectrum and consider only the linear modes by applying a scale cut. We choose to set k_max = 0.1 h Mpc⁻¹, as the next-to-leading order terms become important beyond this scale. In Sect. 7, we discuss this issue and study its impact on the b_ϕ and f_NL constraints, where we find that variations of 50% in k_max compared to our fiducial choice have a negligible effect.

5.2. Constraints on b_ϕf_NL

Following the pipeline described above, we fix all the cosmological parameters for our analysis, except b₁ and b_ϕf_NL in Eq. (11), and then search for the best-fit parameters.

We show the results in the third column of Table 3 and in Fig. 5, where the blue and orange symbols represent the best-fit 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 68% confidence interval obtained from our fits for b_ϕf_NL. For the blue and orange data points, we have used the variances obtained from the fits to the FASTPM mocks (see Sect. 4.2), while for the green symbols, we have used the fixed-and-matched technique (see Sect. 4.3). We expected to obtain exactly zero for the Gaussian simulation, while for the non-Gaussian one, the reference we compare to is the universality relation prediction (i.e., p = 1). The grey lines show these expectations. For halos with b₁ >  1, we clearly measure b_ϕf_NL. When b₁ ∼ 1, for the most popular assumption of p = 1 in Eq. (13), b_ϕ is nearly zero and the non-Gaussian signal is masked. In the same figure, we also observe fluctuations for the f_NL = 0 simulation around zero, which is expected due to the intrinsic cosmic variance of the simulation.

Fig. 5. Constraints on bϕfNL as a function of b1. 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 Sect. 4.2). By subtracting both, we obtain the green points. Thanks to the so-called matching technique, these points have smaller uncertainties (see Eqs. (16) and (17) in Sect. 4.3). Grey lines indicate the expected value by assuming the universality relation, with the dotted line for fNL = 0 and the dashed line for fNL = 100. In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method. Here, mass bins 0–9 are shown.

Since both simulations originate from the same random phases, we expect the measurements of b_ϕf_NL to be highly correlated (Sect. 4.3). Therefore, if we find perturbation in the Gaussian simulation (f_NL = 0) that lowers the measured b_ϕf_NL, an equivalent perturbation is expected to appear in the non-Gaussian (f_NL = 100) simulation, as we see in Fig. 5 indeed.

We apply Eqs. (16) and (17) to compute the matched data points from both simulations, obtaining the green points in Fig. 5. We expect that these results have less noise than the individual cases of the simulations because of the correlations in the initial conditions. We find the green points to be closer to the prediction from the universality relation (p = 1). While this is true for M_halo >  1 × 10¹² h⁻¹ M_⊙, we find that the measured p is below that of the universality relation for masses with M_halo < 1 × 10¹² h⁻¹ M_⊙, with a significance that reaches 3σ for bin 2.

Regarding the mass range from M_halo = 5 × 10¹³ h⁻¹ M_⊙ up to 5 × 10¹⁴ h⁻¹ M_⊙ (bins 10–12 in Table 3, not shown from Fig. 5 on), their power spectra exhibit a very low signal-to-noise ratio (S/N < 5) due to their low number densities (n_halos ≪ 10⁻⁴ h³ Mpc⁻³). As we discussed in Sect. 4.2, to estimate the covariance for the power spectra, we considered two cases: applying Eq. (14) or the extrapolation of Eq. (15). When considering σ_S(P(k)), a more detailed analysis of these bins reveals that considering a super-Poisson contribution to the shotnoise (∝(1 + α)/n_g) is necessary to have an accurate estimation of the variance of the power spectrum, σ(P(k)). At the same time, this extra noise, parameterised by α, is consistent with α = 0 for the rest of the bins (0–9), which is our assumed baseline. We leave for future analysis the usage of non-Poisson errors as our analysis baseline. Thus, in what follows, we only analyse halos with masses 2 × 10¹⁰ < M_halo(h⁻¹ M_⊙)≤5 × 10¹³, which correspond to bins 0–9 in Table 3.

5.3. Measuring f_NL from b_ϕf_NL assuming the universality relation

We cannot constrain both b_ϕ and f_NL simultaneously since they are perfectly degenerate. We can do so, however, by assuming the universality relation with p = 1. In this way, we extract f_NL from our simulation results. If this relation does not accurately describe the b_ϕ parameter, we expect to recover biased values of f_NL. Hence, we have tested the validity of the universality relation.

In Fig. 6, we compare the values of f_NL obtained in this way with the ones set in the initial conditions. We find that the perturbations that favour a positive value of f_NL in the Gaussian simulation are also present in the non-Gaussian one. The results are shown in blue symbols for the f_NL = 0 simulation and orange symbols for the f_NL = 100 simulation. We reduce the noise in f_NL by applying the matching technique. This involves the estimation of b_ϕf_NL using Eq. (16) and its variance through Eq. (17). In the Gaussian simulation, the expected value for b_ϕf_NL is 0. This means that the expected value of Δb_ϕf_NL is the same as the b_ϕf_NL of the simulation with f_NL = 100. Then, we divide by b_ϕ computed assuming p = 1. Using this correlation, we determine the combined measurement of f_NL, shown in green symbols in Fig. 6.

Fig. 6. Value of fNL for each mass bin, assuming bϕ is described by Eq. (13) with p = 1. Blue and orange symbols show results for the Gaussian and non-Gaussian simulations, respectively. In contrast, green symbols show the matching (difference) between them, where we used Eq. (16) for obtaining bϕfNL and Eq. (17) for its uncertainty. Dashed lines show the input fNL values for each simulation. For the mass bins 2 and 3, we have that b1 − 1 ≤ 0.1 and hence bϕ ∼ 0. This value of bϕ explains the large error bars compared to the other mass bins. In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method.

For masses between 5 × 10¹⁰ and 2 × 10¹¹ h⁻¹ M_⊙, the error bars are enlarged because b₁ ∼ 1 for those bins, so using the universality relation results in b_ϕ ∼ 0. Moreover, a slight deviation from the expected value in b_ϕf_NL assuming universality in these mass bins leads to an amplification of this deviation from the input value in f_NL, precisely because b_ϕ is close to 0. This is what we observe for bins 2 and 3. Nevertheless, considering also the growth of the errorbars, for halos with masses between M_halo = 1 × 10¹¹ h⁻¹ M_⊙ and M_halo = 1 × 10¹² h⁻¹ M_⊙ we find deviations from the input value, with a significance up to 4σ in the case of bin 2. These differences, which we have already seen in b_ϕf_NL in Fig. 5, can be interpreted as meaning that the population of halos in that mass bin have b_ϕ that does not follow the universality relation, as we assumed p = 1 for deriving these measurements of f_NL. Finally, we recover the input value of f_NL for both simulations for halos with M_halo >  1 × 10¹² h⁻¹ M_⊙ (see Table 3).

5.4. p as a function of mass

Since we know the input value of f_NL, we can use the deviations of our measurements from that value to find how far the assumption p = 1 is from the p measured for the simulation halos.

We fix the value of f_NL and derive b_ϕ from the previous constraints on b_ϕf_NL we obtained in Sect. 5.2. We then compute p using Eq. (13). We have also checked that we get the same results by fitting p directly and fixing f_NL to its input value in a similar way as done in Sect. 5.3.

The results are shown in Fig. 7, where we compare the different methods of estimating the variance of the power spectrum and reducing the noise introduced in Sect. 4.4. Our methodology is consistent as we recover the same results regardless of the method used for estimating the variance of the power spectrum with or without the matching technique. The clustering of the halos with 1 × 10¹² < M_halo < 2 × 10¹³ h⁻¹ M_⊙ are aptly described through the universality relation (p = 1), given the size of our uncertainties. Concerning the mass bins with M_halo < 5 × 10¹¹, we find a slight deviation from p = 1, with a significance that reach > 3σ for halos with masses between 1 × 10¹¹ and 2 × 10¹¹ h⁻¹ M_⊙. These lower values of p are related to the overprediction of the product b_ϕf_NL (Sect. 5.2).

Fig. 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 (Eq. (14)) and fixed variances, respectively. Green and red symbols show the same after applying the matching with the Gaussian simulation. In solid circles, the fixed variances is estimated by the best fit to the FASTPM variances with Eq. (15); in rings, we use the extrapolation method. The shaded areas represent the 1σ, 2σ and the 3σ regions (from darker to lighter respectively) around the value of p we find in Eq. (19).

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 zero at 1σ, suggesting that a constant p can indeed be assumed for the halos. By an average weighted by the inverse of the variance, we get:

$$\begin{array}{c}\hfill p=0.955\pm 0.013\end{array}$$(19)

with χ²/d.o.f. = 1.31 when we take into account all the bins (0–9), while the χ²/d.o.f. for p = 1 is 2.64. This best-fit value of p with the 1,2 and 3σ regions are shown as the shaded region in Fig. 7. If we neglect those bins where we used extrapolation (i.e., by removing the mass bins 0–2 from the fits and considering just 3–9), we get:

$$\begin{array}{c}\hfill p=0.950\pm 0.025\end{array}$$(20)

with a χ²/d.o.f. of 1.33.

6. Convergence tests

This section shows 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 (Sect. 6.1). We also demonstrate the robustness of our results by comparing them with an alternative method to measure b_ϕ, using the separate universe approach (Sect. 6.2). Finally, we investigate the impact of considering different mass definitions (Sect. 6.3).

6.1. Convergence with LR-UNIT simulations

First, we investigate the mass convergence of our results by comparing the PNG-UNIT results with those of the lower-resolution LR-UNIT simulations. Our methodology is the same as in Sect. 5.2, and we assume σ_F + M errors.

The upper panel of Fig. 8 compares the values of b_ϕ derived from the 4096³-particle simulations with those from the LR-UNIT simulations. The b_ϕ values agree for bins with M_halo ≳ 20m_part, LR. Below this threshold, the mass function from the LR-UNIT simulations drops as we cannot resolve such light halos (Fig. 2). This drop in the abundance of halos also affects the clustering, as demonstrated by the sudden change in b_ϕ. Thus, we prove that the limit of 20 particles is an adequate minimum threshold for the halo mass considered for our simulations.

Fig. 8. Non-Gaussian bias parameter bϕ as a function of the halo mass for the different sets of simulations. The colour of the symbols and the dashed lines represent the mass resolution. Red: mp = 1.5 × 108 [h−1 M⊙]. Green: mp = 1.2 × 109 [h−1 M⊙]. Blue-Dark blue: mp = 1.0 × 1010 [h−1 M⊙]. The dashed vertical lines represent the threshold of Mhalo = 20mpart for each simulation. The pink solid line is the universality relation p = 1 for the b1 values measured in the PNG-UNIT. The dashed pink line extrapolates the universality relation toward lower masses, using b1 coming from Tinker et al. (2010). In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method. Upper panel: measurements of bϕ in non-Gaussian simulations by fitting Eq. (11) (see Sect. 6.1). Lower panel: measurements of bϕ using the separate universe technique (see Sect. 6.2).

The universality relation (p = 1) is shown as a solid pink line in Fig. 8, assuming the b₁ values from the fits of the PNG-UNIT to the corresponding mass resolution. To extend the relation beyond the mass resolution of the N-body simulations, we have also used the theoretical values for b₁(M) derived using Tinker et al. (2010) for each bin in halo mass, given the simulation cosmology. This extrapolation is shown as dashed pink lines in the upper panel of Fig. 8. These two lines are consistent with each other.

The derived b_ϕ values agree in both resolutions for halos with more than 20 particles. Although not shown here, we have checked that the results from running LR-UNIT with 2LPTIC (Crocce et al. 2006, 2012) initial conditions are consistent with those from LR-UNIT initialised with FASTPM. In a forthcoming paper, we study the impact of using different codes for the initial conditions in simulations with primordial non-Gaussianities.

Finally, we measure p from the LR-UNIT simulations by fixing b₁ to Tinker et al. (2010) in Eq. (13). We obtain similar results for p using b₁ directly obtained from the simulations. Combining all the mass bins with M_halo >  20m_part, we obtain p = 0.959 ± 0.014 for the PNG-UNIT and p = 0.972 ± 0.017 for the LR-UNIT. These values are consistent with Eq. (19).

6.2. Comparison with separate universe simulations

As a further test of the robustness of our results, we compare two alternative methods for measuring b_ϕ: (i) using the scale-dependent bias as we have done throughout this paper and (ii) using the so-called separate universe technique. This latter technique exploits the fact that b_ϕ is defined as the response of the abundance of halos to the presence of a large-scale perturbation (Eq. (12)).

The separate universe technique consists of computing b_ϕ as the numerical derivative of two simulations with the same initial conditions but with a slightly different 𝒜_s, or equivalently, σ₈. The initial conditions for these simulations are Gaussian; thus, these simulations do not show the effect of the scale-dependent bias. Following a similar procedure as done in Barreira et al. (2020, 2022c), b_ϕ is obtained as follows:

$$\begin{array}{c}\hfill b_{\varphi }=\frac{b_{\varphi }^{\mathrm{high}}+b_{\varphi }^{\mathrm{low}}}{2},\end{array}$$(21)

with

$$\begin{array}{cc}& b_{\varphi }^{\mathrm{high}}=\frac{4}{\delta {\mathcal{A}}_s{N}_{\mathrm{h}}^{\mathrm{fid}}}(N_{\mathrm{h}}^{\mathrm{high}}-N_{\mathrm{h}}^{\mathrm{fid}}),\hfill \\ \hfill & b_{\varphi }^{\mathrm{low}}=\frac{4}{\delta {\mathcal{A}}_s{N}_{\mathrm{h}}^{\mathrm{fid}}}(N_{\mathrm{h}}^{\mathrm{fid}}-N_{\mathrm{h}}^{\mathrm{low}}).\hfill \end{array}$$(22)

In Eq. (22), δ𝒜_s is the variation in the amplitude of the initial power spectrum (which we set to 5%), and N_h is the number of halos within a particular mass bin measured in the simulations. For the separate universe technique, we estimate the error in b_ϕ, σ_SU(b_ϕ), following the same consideration as Barreira et al. (2020) and Barreira (2022c):

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{SU}}(b_{\varphi })\sim |b_{\varphi }^{\mathrm{high}}-b_{\varphi }^{\mathrm{low}}|.\end{array}$$(23)

The lower panel of Fig. 8 compares the b_ϕ values obtained from using the separate universe technique (Eq. (22)) for two sets of simulations with different mass resolutions. Comparing this method for measuring b_ϕ with respect to the one we have been considering in Sect. 5.2 is a way to check the consistency across the two different methods and whether they hold at low masses.

The universality relation is displayed as pink lines in Fig. 8. We used both the linear bias, b₁, coming from the measurements on the PNG-UNIT simulation (continuous line) and b₁ coming from Tinker et al. (2010) (dashed line). We do not measure b₁ from the separate universe simulations as the box size of the simulations is relatively small (L_box = 67.5 h⁻¹ Mpc and L_box = 250 h⁻¹ Mpc for the high-resolution and low-resolution simulations respectively) and direct measurement of the linear bias through the halo power spectrum may be significantly contaminated by the non-linearities of the small scales (k_f, HR ∼ 0.09 h Mpc⁻¹ and k_f, LR ∼ 0.02 h Mpc⁻¹).

The error bars obtained for the separate universe method might not be robust. On the one hand, we find that the high-resolution separate universe simulations do not agree with the low-resolution simulations, as shown in the lower panel of Fig. 8. The real error bars should, at least, reflect the difference between both. On the other hand, the scatter around the universality relation is much bigger than the error bar reported by this method. This is specially significant for M_halo between 1 × 10¹¹ h⁻¹ M_⊙ and 5 × 10¹¹ h⁻¹ M_⊙.

We estimate a global error in b_ϕ by combining the measurements from both the high and low-resolution simulations, taking the mean value and the absolute difference as the error, similar to what we did to compute b_ϕ from the high-As and the low-As simulations in Eqs. (21) and (23). With this procedure, we expect to obtain more reliable errors than those from a single set of simulations (i.e., high or low-resolution).

Finally, we measure p from the separate universe. We fit b_ϕ to Eq. (13) using a least-squares method, with b₁ coming from Tinker et al. (2010) for each mass bin. Utilising the error in b_ϕ for the combination of the high- and low-resolution separate universe simulations we have described, we get p = 0.8307 ± 0.0022. This result differs from the universality relation and the p obtained for the PNG-UNIT, Eq. (19). We argue that this suggests a shortcoming in the method used to obtain σ_SU(b_ϕ). The mass bins with the smallest errors drive the separate universe measurement of p. Moreover, χ²/d.o.f. = 519/20 for the p obtained with the separate universe technique is rather large.

Ignoring the separate universe error bars, we find p = 1.047 ± 0.057 by not weighting the bins by their σ(b_ϕ) in the least-squares fitting procedure. If we apply this method to the results coming from PNG-UNIT simulation, we get p = 0.984 ± 0.031. In both cases, the value is consistent with the universality relation, and it is less than 2σ away from the results shown in Eq. (19).

Overall, the results in this section suggest that one must be very careful when considering the errors in b_ϕ arising from the separate universe technique and consider a larger set of simulations to deal with them statistically. Given the setup considered in this paper, the error bars reported are more robust for PNG-UNITSIMS, and the underlying true uncertainty seems smaller for PNG-UNITSIMS.

6.3. Dependency on the mass definition

If b_ϕ depends on the halo mass only by its dependence on b₁, then the relationship between these two quantities should be independent of the mass definition used for halos.

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 ρ_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_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 ρ_c for M_500c and 200 times the background density ρ_b for M_200b. For the virial mass M_vir, the spherical overdensity is determined by the virial radius r_vir (Bryan & Norman 1998).

We have M_200c ≥ M_500c for a given halo. Thus, if we apply the same mass cuts for defining the bin, we expect that b₁(M_200c) < b₁(M_500c), as if we compare the M_200c of the halos selected in each case, we are taking more massive halos when we apply the cut in M_500c. However, although we expect that b_ϕ(M_200c)≠b_ϕ(M_500c) if b_ϕ only depends on the mass through b₁, we would expect that p does not vary with the mass definition.

We present our findings in Fig. 9 which compares four mass definitions outputted by ROCKSTAR: our fiducial M_200c, M_500c, M_200b and M_vir. There are several things to point out. We find that bins defined according to M_200b and M_vir are practically equal between them and very similar to M_200c. Moreover, they have almost the same b₁, although this is not shown here. For these three cases, the results are consistent with our findings in Sect. 5.4.

Fig. 9. PNG bias parameter p as a function of Mhalo for different mass definitions. The green symbols represent our choice for the fiducial analysis (M200c), while the grey dashed line represents the universality relation (p = 1). The shaded areas represent the 1σ, 2σ and the 3σ regions (from darker to lighter respectively) around the value of p we find in Eq. (19). Masses are shifted for visualisation purposes.

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 error bars cover p = 1. This result suggests that b_ϕ for dark matter halos may not depend solely on the mass through b₁ but also on other parameters. The b_ϕ parameter is defined as how the abundance of the dark-matter halos responds to a large-scale perturbation. However, PNG may not only affect the amplitude of the primordial perturbations (and hence the abundance of the high-mass halos) but also could affect other properties of the halo as suggested by recent works (Sullivan et al. 2023; Lazeyras et al. 2023; Lucie-Smith et al. 2023; Fondi et al. 2024). Investigating this dependence is beyond the scope of this paper and is left for future research.

We also tested the impact on b_ϕ if we consider the main halos and the subhalos identified by ROCKSTAR. This consideration might be necessary if one runs subhalo abundance matching (SHAM) algorithms to populate the simulation with galaxies. When the substructure is included in the analysis, we find that this induces variations in the linear bias b₁ of the mass bins and in p. We leave a detailed analysis of this for future work.

7. Robustness tests

Here, we test the impact of the assumptions we have made in our analysis to obtain b_ϕ or p. In particular, we study the impact of the scale cut, k_max (Sect. 7.1), and the variance estimation methods described in Sect. 4.4.

7.1. Choice of k_max

Equation (11) is only valid in the linear regime, thus, we must discard all modes with k greater than a certain k_max to avoid introducing non-linear modes our model does not properly describe. As perturbation theory breaks beyond k_NL ∼ 0.3 h Mpc⁻¹ (Desjacques et al. 2018), k_max should be less than this value. Even before reaching this threshold, we are at the risk of incorporating some mildly non-linear modes that our model does not adequately describe.

To determine the optimal k_max, we must first examine the limits of the linear theory using only the Gaussian simulation. In this case, Eq. (11) 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², where A is a nuisance parameter (Lazeyras & Schmidt 2018; Lazeyras et al. 2023). To establish k_max, we compare the fits obtained using only b₁ with those obtained using b₁ + Ak² after marginalising over the nuisance parameter A as a function of k_max. We find that the differences in the posterior distributions for b₁ obtained using both methods agree for k_max ≤ 0.1 h Mpc⁻¹ for all mass bins.

We proceed similarly with the non-Gaussian simulation. Starting from our fiducial choice of k_max = 0.1 h Mpc⁻¹, we test both the lower and higher values. Our results show that varying the k_max has a minimal impact on p, with all central values remaining within 1σ when k_max is varied from 0.07 to 0.15 h Mpc⁻¹. This is shown in panel a of Fig. 10. The error bars shrink by approximately ∼20% throughout this range when considering a larger k_max. Furthermore, we performed a linear fit for p as a function of mass (as described in Sect. 5.4) and we find that the slope is consistent with zero. If we compute the weighted average of bins 0–9, we find that p(k_max = 0.07 h Mpc⁻¹) = 0.933 ± 0.013 and p(k_max = 0.15 h Mpc⁻¹) = 0.947 ± 0.012, in agreement with our previous findings for our fiducial choice of k_max.

Fig. 10. PNG bias parameter p as a function of Mhalo for different analysis choices. The green symbols represent the fiducial choices, while the grey dashed line represents the universality relation (p = 1). In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method. The shaded areas represent the 1σ, 2σ and the 3σ regions (from darker to lighter respectively) around the value of p we find in Eq. (19). Masses are shifted for visualisation purposes. Panel a: variation in p with changes kmax (see Sect. 7.1). Panel b: constraints on p by considering different methods to compute σ(P(k)), namely: standard σS(P(k)) (blue), direct fit to FASTPM mocks (orange), best-fit to Eq. (15) (green) and extrapolation of Eq. (15) to the mass bins (red) (see Sect. 7.2).

Going beyond the purely linear modes may bias the constraints on b₁ (since those modes are not properly modelled) and hence the constraints on b_ϕ and p. We therefore set the cutoff at k_max = 0.1 h Mpc⁻¹. On the one hand, we have shown that the results are consistent even if this value is varied by ∼50%. On the other hand, this is a standard choice, and it is more conservative than continuing to push 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, σ(P(k)), is critical when we perform parameter estimation, given that this quantity enters the computation of the likelihood. However, we know that the standard approach of Eq. (14) overestimates the variance of our simulations as they have fixed initial conditions. In Sect. 4.2, we discuss how we estimated σ(P(k)) from the FASTPM mocks, which take into account the effect of the fixed ICs. We show how, by fitting σ(P(k))/P(k) to Eq. (15), we can even extrapolate to masses where we are not able to generate the ‘equivalent’ mass bins to the ones from the N-body simulation.

Here, we check the impact of considering different ways of computing the power spectrum variance on p. In particular, we consider the four different methods (see Sect. 4.2). First of all we consider the standard variances, σ_S(P(k)), computed from Eq. (14) (blue symbols in panel b of Fig. 10). Then we consider the direct FASTPM variances, sigma_F(P(k)), obtained by computing directly the standard deviation of the power spectrum of the fast simulations (orange symbols). We also consider the best-fit of FASTPM variances using Eq. (15) (green symbols). And finally, we consider the extrapolation of Eq. (15) to the mean halo mass of the bins (red symbols).

Our baseline analysis adopts the third method for bins 3–9, where we use the best fit of Eq. (15) for each one, and the fourth method for bins 0–2, which are beyond the resolution limits of the FASTPM mocks and hence, we cannot apply the third method directly.

We present the results of our comparison in panel b of Fig. 10, which shows the PNG bias parameter p as a function of M_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 σ_F(P(k)) are shown. Regardless of the method used to estimate the variance, we recover the same results within 1σ. Moreover, the FASTPM approach is more accurate than the standard approach as it considers the fixed initial conditions in our simulations. Additionally, the best fit of σ_F(P(k))/P(k) provides almost the same results compared to using the extrapolation method. These results are consistent with using the variances directly obtained from FASTPM mocks, even considering their small error bars. 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σ with the standard approach for all mass bins.

We now consider the fit to mass bins 0–9. We performed a linear fit on p and found that it is consistent with being a constant, confirming the results of Sect. 5.4. Furthermore, the weighted average yields p(σ_S + M) = 0.950 ± 0.056 and p(direct FASTPM σ_F+M) = 0.937 ± 0.021, which agrees within 1σ with our fiducial choice.

8. Constraining f_NL with priors on p

In this section, we study how setting a prior in p (or equivalently, on b_ϕ) affects the constraints on f_NL. This is because, with the scale-dependent bias, we are only sensitive to the product b_ϕf_NL.

Here, we assume a Gaussian prior for p centred on p_prior and with a width σ(p). To implement the priors in our MCMC procedure (see Sect. 5.1 for details), we modify the Gaussian likelihood by introducing this a second term in χ²:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=0}^{N_k}}\frac{{(P_{\mathrm{halos}}(k_i)-P_{\mathrm{model}}(k_i,b_{\varphi }{f}_{\mathrm{NL}},b_1))}^2}{{\sigma }^2(P(k_i))}+\frac{{(p-p_{\mathrm{prior}})}^2}{{\sigma }^2(p)}.\end{array}$$(24)

In Sect. 8.1, we demonstrate how increasing the prior width σ(p) deteriorates the constraints in f_NL. In Sect. 8.2, we show how using these priors can affect the measurement of f_NL in future cosmological surveys.

8.1. Effect of the prior size on σ(f_NL)

This section examines how the width of the prior for p affects the determination of f_NL. Nine mock halo power spectra are generated based on Eq. (11). CAMB is used to compute the matter power spectrum in each case. Three f_NL values are considered: 0, −12, and 50, with b₁ set to 1.4 for all cases. The choice of f_{NL,  fid 1} = 0 is motivated by the current tightest constraint for f_NL by Planck Collaboration IX (2020); f_{NL,  fid 1} = −12 and f_{NL,  fid 3} = 50 are motivated by two recent constraints coming from galaxy clustering analysis (Mueller et al. 2022; Rezaie et al. 2024). The values of b_ϕ are derived in each case from Eq. (13), assuming that p is 0.5, 1, and 1.6. MCMC analyses are conducted to measure f_NL and its uncertainty σ(f_NL) for each of the nine mocks.

We analyse the data in two ways. First, we fix p to its fiducial value, assuming perfect knowledge of this quantity. Second, we allow p to vary, introducing a Gaussian prior centred on the fiducial value of p with a varying width. The idealised scenario of having a fixed p, has only two free parameters: b₁ and f_NL. This scenario yields the most stringent constraints. Allowing p to vary requires a prior on this parameter to counter the degeneracy between p and f_NL. This is further discussed below.

In Fig. 11, we show the uncertainty σ(f_NL) for the setup described above as a function of the width of the prior in p, σ(p), for the case of f_NL = 0. The analysis for the three fiducial p values are shown in different colours. The values of σ(f_NL) from assuming fixed values for p are shown as horizontal dashed lines. Those assuming a prior in p are shown as solid lines.

Fig. 11. Statistical error in fNL as a function of the size of the prior on p, for three values for p, shown in different colours, as indicated in the legend. The mock data vector for this figure is generated using the linear power spectrum computed with CAMB in Eq. (11) and fixing b1 = 1.4 and fNL,  fid = 0. The horizontal dashed lines show the case where we performed the fits fixing p in the model to its corresponding fiducial value. The dashed vertical line displays the value of the error in p, obtained for halos with mass between 5 × 1011 and 1 × 1012 h−1 M⊙ (bin 4 in Table 3).

Figure 11 shows that constraints on f_NL are similar when considering either a fixed value of p or a prior, as long as this prior is narrow enough, σ(p)≲0.1. Beyond this threshold, the constraints weaken significantly, increasing in σ(f_NL). As we can see in Fig. 11, this threshold depends on the fiducial value of p. In fact, it also depends on b₁ and b_ϕ. For a fixed b₁, the lower the p (and hence the higher the b_ϕ), the tighter the constraints on f_NL. In Fig. 11 we only show results for f_NL = 0, but similar trends are found for f_NL = −12 and 50.

8.2. Forecast for a DESI-like survey

We go on to study the impact on the determination of f_NL in the context of current galaxy surveys when assuming either the universality relation or considering priors for p (or b_ϕ). We base our forecast on the clustering of the typical halos hosting cosmological tracers.

For this forecast, we assume survey properties similar to DESI, whose footprint will cover an area in the sky of ∼14 000 deg² (DESI Collaboration 2023). We consider the three cosmological tracers targeted by the DESI dark-time programme: luminous red galaxies (LRGs), emission-line galaxies (ELGs) and quasars (QSOs). For our analysis, we use the results given in Yu et al. (2024) for the DESI early data release (EDR), as they are based the authors’ analysis on the Gaussian UNIT simulation. Yu et al. (2024) found that DESI EDR LRGs, ELGs, and QSOs populate halos with typical masses of log₁₀(M_{halo,  LRG}) = 13.16, log₁₀(M_{halo,  ELG}) = 11.90, and log₁₀(M_{halo,  QSO}) = 12.66, respectively. These masses correspond to our halo mass bins 8, 4, and 6, respectively (see Table 3). Additionally, each tracer is restricted to a certain range of redshifts. The redshift range for LRGs is from z = 0.4 to z = 1.1, with an effective redshift of z_eff LRG = 0.814. For ELGs, the range goes from z = 0.8 to z = 1.6 with an effective redshift of z_eff ELG = 1.202. And for QSOs, we have considered the range from z = 0.8 to z = 2.5 with an effective redshift of z_eff QSO = 1.741. We computed each tracer’s total volume, assuming a spherical shell between the redshift ranges and the entire DESI footprint. In this way, we get V_{tot,  LRG} = 19.45 h⁻¹ Gpc, V_{tot,  ELG} = 34.37 h⁻¹ Gpc, and V_{tot,  QSO} = 64.94 h⁻¹ Gpc. The number densities Yu et al. (2024) found for the considered redshift range are n_LRG = 5.50 × 10⁻⁴ h³ Mpc⁻³, n_ELG = 7.26 × 10⁻⁴ h³ Mpc⁻³ and n_QSO = 2.4 × 10⁻⁵ h³ Mpc⁻³. For the linear bias we assume b_1,  LRG = 1.2/D(z), b_1,  ELG = 0.84/D(z), and b_1,  QSO = 1.2/D(z), where D(z) is the growth factor (DESI Collaboration 2016, 2023). At the effective redshifts of interest, this leads to b_1,  ELG = 1.503, b_1,  LRG = 2.568, and b_1,  QSO = 2.628.

We generated the mock halo power spectrum using CAMB in Eq. (11). As described in the previous section, here we consider three fiducial values for f_NL: f_{NL,  fid 1} = 0, f_{NL,  fid 2} = −12, and f_{NL,  fid 2} = 50. Regarding b_ϕ, we assume it follows the form given by Eq. (13). We use the values of p measured in this work for the mass bin corresponding to each tracer¹⁰: p_LRG = 1.11 ± 0.27, p_ELG = 0.895 ± 0.055 and p_QSO = 0.973 ± 0.077. The minimum wavenumber considered for the power spectrum is given by the fundamental wavenumber, $k_{\min ,\mathrm{tracer}}=2\pi /V_{\mathrm{tot},\mathrm{tracer}}^{1/3}$. The corresponding covariance is calculated using Eq. (14), considering the survey parameters described above for each tracer.

We used an MCMC fitting procedure to determine the impact on the measurement of f_NL from assuming a prior for p or not. We employed the same methodology as for our main results. In this analysis, all the cosmological parameters are fixed, except for b₁, f_NL, and b_ϕ. When fixing b_ϕ to the universality relation value, we defined the value of χ² using Eq. (18). When considering p as a free parameter, we incorporated a Gaussian prior centred on p_LRG, p_ELG and p_QSO for each case value and a width given by our errors. In this case, χ² is defined by Eq. (24).

In Table 4, we summarise our findings. Fixing p to the universality relation induce a bias of 1.9σ for f_NL fid = 50 and 0.9σ for f_NL fid = −12 for ELGs. These biases arise because the halos hosting ELGs have a value for p below one by ∼10%. For LRGs, we also find a 0.8σ difference for the case with f_NL fid = 50. Finally, for the QSO hosts, we correctly recovered the fiducial values of f_NL given that the value of p for this sample is consistent with p = 1.

When we fix p = p_fid, we recovered unbiased results for f_NL as expected. Nevertheless, the errors in f_NL increase or decrease if p_fid is higher or lower compared to the case with p = 1. These are the cases of the ELG hosts and LRG hosts, respectively. This variation in σ(f_NL) is because the larger b_ϕ is, the tighter the constraints on f_NL are, as we are only sensitive to the product of both quantities. We found negligible differences for the QSO hosts.

Finally, we recovered unbiased results for the three tracers when using our priors while, at the same time, the constraining power is comparable to fixing p to its fiducial value. In the case of the ELG hosts, our constraints are equal for $f_{\mathrm{NL}}^{\mathrm{fid}}=0$, deteriorating by 7% and a 27% for $f_{\mathrm{NL}}^{\mathrm{fid}}=-12$, and $f_{\mathrm{NL}}^{\mathrm{fid}}=50$, respectively. For the LRG hosts, our priors perform worse, deteriorating the constraints by a 9%, a 34% and a 85% for $f_{\mathrm{NL}}^{\mathrm{fid}}=0$, $f_{\mathrm{NL}}^{\mathrm{fid}}=-12$, and $f_{\mathrm{NL}}^{\mathrm{fid}}=50$, respectively. Nevertheless, we are still able to constrain f_NL without recovering the degeneration between this parameter and p. Finally, our priors for the QSO hosts are comparable to fixing p = 1, providing even tighter constraints for $f_{\mathrm{NL}}^{\mathrm{fid}}=-12$.

9. Summary and conclusions

Detecting primordial non-Gaussianities (PNG) would have a major impact on how we understand inflation. The key to measuring or accurately constraining PNG in the near future is to understand the scale-dependent bias parameter, b_ϕ, and how it varies for different cosmological tracers.

In this paper, we study the scale-dependent bias b_ϕ for dark matter halos with masses from 2 × 10¹⁰ h⁻¹ M_⊙ to 5 × 10¹⁴ h⁻¹ M_⊙. To achieve this, we have developed PNG-UNIT, a new full N-body simulation with local non-Gaussian initial conditions (f_NL = 100). PNG-UNIT is the largest simulation to date that incorporates non-Gaussian initial conditions. The simulation has 4096³ dark matter particles in a 1 (h⁻¹ Gpc)³ box, which was run using the L-GADGET2 code. The cosmology of the simulation is based on Planck Collaboration XIII (2016) (see Table 1). PNG-UNIT assumes local PNG with an amplitude of f_NL = 100. We have chosen a large value of f_NL to increase the signal from PNG while reducing the errors on b_ϕ. We significantly reduced the error on b_ϕ by using the method proposed in Avila & Adame (2023) for simulations with fixed-and-matched initial conditions. 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).

In addition to the main simulation, PNG-UNIT, we also developed supporting simulations: a set of full N-body simulations, either with lower resolution or smaller volumes than PNG-UNIT, along with 100 fast mocks, half of which have f_NL = 100 and the other half with f_NL = 0, generated with the approximated method FASTPM. We refer to the whole set of simulations as PNG-UNITSIMS.

The PNG-UNITSIMS suite (see Sect. 3.2) is a unique laboratory for studying galaxy clustering models in the presence of local PNGs. The simulations in this suite have effective volumes comparable to current spectroscopic galaxy surveys. The reference simulation, PNG-UNIT, has a mass resolution high enough to accurately resolve the dark matter halos expected to host the galaxies currently targeted by surveys such as DESI and Euclid. Thus, the PNG-UNITSIMS suite is a powerful tool to develop and test analysis pipelines to measure f_NL with current spectroscopic galaxy surveys.

In this first study of the PNG-UNITSIMS suite, we were able to constrain the non-Gaussian bias parameter b_ϕ using mass-selected samples of dark matter halos. We inferred the scale-dependent bias from the halo power spectrum measured for the Gaussian and the non-Gaussian simulations. From the power spectra, we inferred b₁ and b_ϕf_NL and correlated the measurements in both simulations Sect. 5. Moreover, given that we know the input values of f_NL, we can derive b_ϕ or p. Our findings can be summarised as follows:

• We measured b_ϕf_NL using the fix-and-match technique (Sect. 4.3) by combining the information of the Gaussian and the non-Gaussian simulation. This technique allows us to cancel out much of the cosmic variance (Sect. 5.2).

• By assuming the universality relation, we recovered the input values of f_NL for our simulations within 1σ for halos between 1 × 10¹² h⁻¹ M_⊙ and 1 × 10¹⁴ h⁻¹ M_⊙. The power spectra are very noisy for higher halo masses due to their low number densities. As the variance estimates are not robust for these massive halos, we have not included them in the rest of our analyses.

• For halos with masses between 5 × 10¹⁰ h⁻¹ M_⊙ and 1 × 10¹² h⁻¹ M_⊙ we find deviations from the input value of f_NL = 100, which reaches a significance of 4σ for halos between 1 × 10¹¹ h⁻¹ M_⊙ and 2 × 10¹¹ h⁻¹ M_⊙, suggesting that the value of p assumed (p = 1) does not accurately described the clustering of these halos.

• We measured p by fixing f_NL to the input values. For halos with masses higher than 1 × 10¹² h⁻¹ M_⊙ and between 2 × 10¹⁰ h⁻¹ M_⊙ and 5 × 10¹⁰ h⁻¹ M_⊙, we find p to be consistent with the universality relation. For halos with masses between 5 × 10¹⁰ h⁻¹ M_⊙ and 1 × 10¹² h⁻¹ M_⊙, we find p to be ∼10% below the universality relation, with a significance between 1.5σ and 3.1σ.

• Combining the information of b_ϕf_NL from mass bins 0–9 (M_halo between 2 × 10¹⁰ h⁻¹ M_⊙ and 5 × 10¹³ h⁻¹ M_⊙) we find a preferred value of p = 0.955 ± 0.013 (68% c.l.). This result is more than 3σ away from the universality relation (i.e., p = 1).

We assured the validity of our results for M_halo ≳ 2 × 10¹⁰ h⁻¹ M_⊙ by checking that the convergence with the LR-UNIT simulations (2048³) occur at M_halo ∼ 20m_part. Moreover, we also made a comparison with the separate universe technique, but the priors reported by the PNG-UNIT simulation are more robust. For the set-up considered here, we find that the separate universe tends to underestimate the uncertainty reported on b_ϕ. By looking at the scatter and lack of convergence of the separate universe results, we also conclude that the uncertainty associated with those b_ϕ measurements would be larger than the one found by PNG-UNITSIMS.

We explored how this affects variations on the mass definition used for the mass cuts to the measurements of p, finding consistency for all mass definitions except for M_500c (Sect. 6.3). We argue that b_ϕ might not depend only on the halo mass through b₁ but also on other parameters. We leave a more detailed exploration of these dependencies for future work.

We have also tested the robustness of the above results to variations in the choice of the scale cuts (Sect. 7.1) and the way we compute the variances of the power spectrum (Sect. 7.2). None of the above is shown to have a significant impact.

Finally, we construct priors on the scale-dependent bias for DESI galaxies based on the simplified assumption that those tracers can be identified directly by the mass of their host halo at z = 1. We show that these priors are:

1. required to get unbiased constraints on f_NL (compared to the universality relation that may bias the results, up to 2σ for the setup considered);

2. sufficiently tight to allow us to constrain f_NL with DESI at a comparable precision with respect to cases where p is fixed to the fiducial value.

All these results show the capability of the PNG-UNITSIMS suite as a valuable tool for interpreting observations with DESI and Euclid. Specifically, we address the issue of putting priors on b_ϕ. While we have analysed how the halo power spectrum can be used to constrain the scale-dependent bias in this work, the PNG-UNITSIMS suite is a unique tool to test and develop alternative methods for constraining PNG with cosmological surveys. In the future, we plan to include galaxies in our simulations and consider other clustering statistics such as the configuration space correlation functions, higher order statistics including the bispectrum or trispectrum, and one-point statistics such as counts-in-cells, density peaks, and galaxy cluster counts.

Acknowledgments

We thank V. Yankelevich, O. Hahn, M. Manera, and R. Scoccimarro for the valuable discussions on how to solve some early problems with initial conditions. We thank the Benasque science center for fostering useful discussions regarding this work. 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). AK further thanks Embrace for ‘gravity’. 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. 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, RES-AECT-1-0007, RES-AECT-2022-3-0030, RES-AECT-2023-1-0014, and RES-AECT-2023-2-0006.). Data availability: The halo catalogues used in this paper are available through the UNITSIM website (http://www.unitsims.org/) The other data products will be shared under email request to the corresponding author.

References

All Tables

All Figures

Fig. 1. Theoretical power spectrum obtained from Eq. (11). The blue line shows the linear matter-matter power spectrum for the UNIT cosmology at z = 1. The b1 parameter is the same for all cases, and we chose it to be similar to the one expected for the DESI emission-line galaxies (ELGs) at that redshift (DESI Collaboration 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 fNL = 100 and a PNG-response parameter bϕ given by Eq. (13) with p = 1. The red line displays how a different value of fNL can lead to the same signal if the bϕ parameter does not follow the universality relation. Red and green lines are on top of each other since the product bϕfNL is the same. Finally, the purple line shows the kind of signal we expect in case bϕfNL < 0.

In the text

	Fig. 2. Halo mass functions for the simulations described in Section 3 (top). Ratio of halo mass functions relative to the Gaussian UNIT simulation as the reference (bottom). Shaded areas represent the $M=20m_{\mathrm{part}}^{\mathrm{LR}}$ mass threshold for the low-resolution simulation. The minimum mass shown here represents the 20 particle limit for the higher resolution simulations.
In the text

Fig. 3. Solid lines represent the power spectrum of PNG-UNIT halos within the different halo mass bins (see Table 3 for the mass ranges). Shaded regions show the standard deviation around the mean power spectrum of the FASTPM mocks, generated using halos in mass bins that reproduce the clustering and the abundance of the N-body simulation (see Sect. 4.1). There are no halos in the two lowest-mass bins (1 and 2, shown in dark blues) in the FASTPM mocks since the mass resolution is not high enough to resolve them. The top panel shows the fNL = 0 case, and the bottom panel shows the fNL = 100 one.

In the text

	Fig. 4. Variance over mean power spectrum for each halo mass bin (indicated by colour with the mass ranges defined in Table 3). Markers with error bars show the data taken directly from the FASTPM mocks, while solid lines show results obtained by fitting these curves to Eq. (15). The extrapolation to each mass bin as described in Sect. 4.2 is shown in dashed lines.
In the text

Fig. 5. Constraints on bϕfNL as a function of b1. 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 Sect. 4.2). By subtracting both, we obtain the green points. Thanks to the so-called matching technique, these points have smaller uncertainties (see Eqs. (16) and (17) in Sect. 4.3). Grey lines indicate the expected value by assuming the universality relation, with the dotted line for fNL = 0 and the dashed line for fNL = 100. In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method. Here, mass bins 0–9 are shown.

In the text

Fig. 6. Value of fNL for each mass bin, assuming bϕ is described by Eq. (13) with p = 1. Blue and orange symbols show results for the Gaussian and non-Gaussian simulations, respectively. In contrast, green symbols show the matching (difference) between them, where we used Eq. (16) for obtaining bϕfNL and Eq. (17) for its uncertainty. Dashed lines show the input fNL values for each simulation. For the mass bins 2 and 3, we have that b1 − 1 ≤ 0.1 and hence bϕ ∼ 0. This value of bϕ explains the large error bars compared to the other mass bins. In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method.

In the text

Fig. 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 (Eq. (14)) and fixed variances, respectively. Green and red symbols show the same after applying the matching with the Gaussian simulation. In solid circles, the fixed variances is estimated by the best fit to the FASTPM variances with Eq. (15); in rings, we use the extrapolation method. The shaded areas represent the 1σ, 2σ and the 3σ regions (from darker to lighter respectively) around the value of p we find in Eq. (19).

In the text

Fig. 8. Non-Gaussian bias parameter bϕ as a function of the halo mass for the different sets of simulations. The colour of the symbols and the dashed lines represent the mass resolution. Red: mp = 1.5 × 108 [h−1 M⊙]. Green: mp = 1.2 × 109 [h−1 M⊙]. Blue-Dark blue: mp = 1.0 × 1010 [h−1 M⊙]. The dashed vertical lines represent the threshold of Mhalo = 20mpart for each simulation. The pink solid line is the universality relation p = 1 for the b1 values measured in the PNG-UNIT. The dashed pink line extrapolates the universality relation toward lower masses, using b1 coming from Tinker et al. (2010). In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method. Upper panel: measurements of bϕ in non-Gaussian simulations by fitting Eq. (11) (see Sect. 6.1). Lower panel: measurements of bϕ using the separate universe technique (see Sect. 6.2).

In the text

Fig. 9. PNG bias parameter p as a function of Mhalo for different mass definitions. The green symbols represent our choice for the fiducial analysis (M200c), while the grey dashed line represents the universality relation (p = 1). The shaded areas represent the 1σ, 2σ and the 3σ regions (from darker to lighter respectively) around the value of p we find in Eq. (19). Masses are shifted for visualisation purposes.

In the text

Fig. 10. PNG bias parameter p as a function of Mhalo for different analysis choices. The green symbols represent the fiducial choices, while the grey dashed line represents the universality relation (p = 1). In solid circles, we use the best fit to Eq. (15) for estimating the power spectrum variance; in rings, we use the extrapolation method. The shaded areas represent the 1σ, 2σ and the 3σ regions (from darker to lighter respectively) around the value of p we find in Eq. (19). Masses are shifted for visualisation purposes. Panel a: variation in p with changes kmax (see Sect. 7.1). Panel b: constraints on p by considering different methods to compute σ(P(k)), namely: standard σS(P(k)) (blue), direct fit to FASTPM mocks (orange), best-fit to Eq. (15) (green) and extrapolation of Eq. (15) to the mass bins (red) (see Sect. 7.2).

In the text

Fig. 11. Statistical error in fNL as a function of the size of the prior on p, for three values for p, shown in different colours, as indicated in the legend. The mock data vector for this figure is generated using the linear power spectrum computed with CAMB in Eq. (11) and fixing b1 = 1.4 and fNL,  fid = 0. The horizontal dashed lines show the case where we performed the fits fixing p in the model to its corresponding fiducial value. The dashed vertical line displays the value of the error in p, obtained for halos with mass between 5 × 1011 and 1 × 1012 h−1 M⊙ (bin 4 in Table 3).

In the text