© The Authors 2026

1. Introduction

Cosmic inflation–a phase of accelerated expansion during the first moments of the Universe–was postulated as a solution to some of the problems of the Hot Big Bang theory associated with the initial conditions (ICs; Guth 1981; Linde 1982). It also explains the origin of the inhomogeneities giving rise to the large-scale structure of the Universe (LSS). Although inflation is assumed to be part of the standard model in most modern cosmological analyses, the exact physics driving it remains unknown, as a wide variety of physical models can produce inflation.

One observable that discriminates among families of models is primordial non-Gaussianity (PNG). We focus on the simplest PNG, local-f_NL (Komatsu & Spergel 2001), which introduces quadratic expansion of a Gaussian gravitational field, ϕ_G,

$$\begin{array}{c}\hfill \varphi (\mathbf{x})={\varphi }_G(\mathbf{x})+f_{\mathrm{NL}}({\varphi }_G{(\mathbf{x})}^2-⟨{\varphi }_G{(\mathbf{x})}^2⟩),\end{array}$$(1)

where ⟨ ⋅ ⟩ represents a spatial average. Typically, this parameter is predicted to be f_NL > 𝒪(1) for inflation driven by multiple fields (Byrnes & Choi 2010; Alvarez et al. 2014).

Although the primary effect of local-f_NL is to generate a primordial bispectrum, the non-linear evolution of the LSS allows us to study PNG through two-point statistics. Galaxies (g) or haloes (h) are biased tracers (t) of the matter (m) density. At large scales, this is given by a linear relation of overdensity (δ):

$$\begin{array}{c}\hfill {\delta }_t=b·{\delta }_m.\end{array}$$(2)

For Gaussian initial conditions, the large-scale bias is simply a constant, b₁, the linear bias. However, when including PNG, the bias becomes scale-dependent (Dalal et al. 2008; Slosar et al. 2008; Matarrese & Verde 2008):

$$\begin{array}{c}\hfill b(k,z)=b_1(z)+\frac{f_{\mathrm{NL}}{b}_{\varphi }(z)}{\alpha (k,z)}.\end{array}$$(3)

Here b_ϕ is a bias parameter associated with the gravitational field, ϕ, related to the overdensity field via

$$\begin{array}{c}\hfill \begin{array}{c}\delta (k,z)=\alpha (k,z)\varphi (k,z),\hfill \\ \mathrm{with}\alpha (k,z)=k^{2}T(k)·\frac{2D(z)}{3{\mathrm{\Omega }}_M}\frac{c^2}{H_0^2},\hfill \end{array}\end{array}$$(4)

where D(z) is the growth factor normalised to (1 + z)⁻¹ in the radiation era, T(k) is the transfer function (with T(k → 0)→1), Ω_M is the current matter density, H₀ the Hubble constant, and c the speed of light. The PNG universality relation is given by (Dalal et al. 2008)

$$\begin{array}{c}\hfill b_{\varphi }=2{\delta }_{\mathrm{crit}}(b_1-p)\mathrm{with}p=1,\end{array}$$(5)

where δ_crit is the critical density for spherical collapse (1.686).

Many galaxy surveys, such as DESI, Euclid, or SPHEREx, aim to constrain f_NL with Eq. (3). Identifying observables that maximise the PNG signal would be highly valuable in order to better understand inflation. In that context, Castorina et al. (2018) postulated that unclustered (b₁ = 0) tracers would be optimal in a regime with negligible shot noise (P_shot), with an error given by

$$\begin{array}{c}\hfill \sigma {(f_{\mathrm{NL}})}^{-2}\propto \frac{b_1^2{b}_{\varphi }^2}{{(b_1^2{P}_{\mathit{mm}}(k)+P_{\mathrm{shot}})}^2},\end{array}$$(6)

with P_mm(k) the matter power spectrum.

The way to find bias-zero tracers proposed in that work is by binning tracers according to their local density. They showed that unclustered tracers can be found in dark matter halo catalogues from cosmological N-body simulations with Gaussian ICs.

The goal of this letter is to verify or disprove the feasibility of using bias-zero tracers selected by local overdensity to set strong constraints on PNGs. For that, we use state-of-the-art N-body simulations with PNGs from the PNG-UNITsim suites (Adame et al. 2024).

2. Simulations

The PNG-UNIT simulation is an N-Body simulation (Adame et al. 2024), which is a twin of one of the UNITsims (Chuang et al. 2019), with the same size (L = 1 h⁻¹ Gpc) and resolution (N_part = 4096³, m_p = 1.2 × 10⁹ h⁻¹ M_⊙), and that uses the same white noise for the ICs and ΛCDM parameters (based on Planck Collaboration XIII 2016); the only difference is the inclusion of a PNG of local-f_NL = 100. This makes it one of the PNG simulations with the highest mass resolution and one of the largest in terms of particles. For this work we used the z = 1 ROCKSTAR (Behroozi et al. 2013) halo catalogues. These simulations are publicly available¹.

The PNG-UNITsim-XL we used here is a L = 3 h⁻¹ Gpc simulation with N = 4096³ particles and initial conditions with f_NL = −20, which is a twin of a f_NL = 0 UNITsim-XL simulation. These simulations are fully described in Adame et al. (in prep.). Both these simulations have been populated with a halo occupation distribution (HOD) model (Avila et al. 2020) matching the galaxy bias of DESI’s Luminous Red Galaxies (LRG, using a snapshot at z = 0.74) and quasars (QSO, using z = 1.83). The HOD is first set to have a maximum of ∼1 for central galaxies and is later downsampled (by a factor of ∼2 and ∼40 for LRGs and QSOs, respectively) to match the number density of DESI. In this work we show the results for the full density samples, while we checked that results are consistent, but noisier, for the realistic-density sample.

3. Methodology

The potential of bias-zero tracers is illustrated in Fig. 1: an unclustered tracer in a Gaussian universe (b₁ = 0, black line) would cluster in PNG universes (green line) with a 1/k² pattern. In this section we describe the methodology to test that.

Fig. 1. Theoretical auto-power spectra for matter (b = 1, fNL = 0), a biased tracer in a Gaussian Universe (bg = 2, fNL = 0), a biased universal tracer with PNG (bg = 2, fNL = 100, p = 1), an unclustered but universal tracer (bg = 0, fNL = 100, p = 1), and an unclustered and non-PNG-responsive tracer (bg = 0, bϕ = 0).

3.1. Density classification

Unclustered tracers are expected to be found in regions that are neither very dense nor very empty. We follow the approach of Castorina et al. (2018) and classify tracers according to their density in spheres of radius R, δ_t, R, using by default R = 8 h⁻¹ Mpc. For this purpose, we use the SCIPY function CKDTREE. We considered different options for the parent samples:

• All haloes in the PNG-UNITsims above 20 particles;

• PNG-UNItsim haloes split into three mass bins: 10.5 < log M < 11, 11 < log M < 12, and log M > 12;

• LRGs and QSOs from a HOD applied to PNG-UNITsim-XL.

For this paper we focused on equal-number δ_t, R bins, starting from the parent catalogues. We used 20 bins for haloes and 10 bins for galaxies. For better comparison with Castorina et al. (2018), we also reproduced their method based on only splitting the parent sample into two local density subsamples. We then kept the lower one: ${\delta }_{t,R}<{\delta }_{t,R}^{\max }$. As we vary the parameter ${\delta }_{t,R}^{\max }$, we measure the linear bias of the sample, until we find b₁ ≈ 0.

3.2. Parameter fitting

For each subsample of interest, we compute the tracer auto-power spectrum P_tt(k) and the tracer-matter cross-power spectrum P_tm(k) with NBODYKIT (Hand et al. 2018). We additionally compute the matter power spectrum P_mm(k), so that we can estimate the bias parameters by fitting the equations

$$\begin{array}{c}\hfill \begin{array}{c}{P}_{\mathit{tm}}(k)=(b_1+f_{\mathrm{NL}}{b}_{\varphi }\alpha {(k)}^{-1})P_{\mathit{mm}}(k),\hfill \\ P_{\mathit{tt}}(k)={(b_1+f_{\mathrm{NL}}{b}_{\varphi }\alpha {(k)}^{-1})}^2{P}_{\mathit{mm}}(k)+\frac{A_{\mathrm{sn}}}{n},\hfill \end{array}\end{array}$$(7)

where we fit the free parameters b₁, b_ϕ, and A_sn. Here A_sn represents the shot noise deviation from the Poisson prediction (1/n).

We use a Gaussian diagonal covariance (Smith 2009)

$$\begin{array}{c}\hfill \begin{array}{c}{\sigma }_{\mathit{tt}}{(k)}^2={(P_{\mathit{tt}}(k)+\frac{A_{\mathit{sn}}}{n})}^2·\frac{2{(2\pi )}^3}{4\pi k^2\mathrm{\Delta }kV}\mathrm{and}\hfill \\ {\sigma }_{\mathit{tm}}{(k)}^2=[(P_{\mathit{tt}}(k)+\frac{A_{\mathit{sn}}}{n})P_{\mathit{mm}}(k)+P_{\mathit{tm}}{(k)}^2]\frac{{(2\pi )}^3}{4\pi k^2\mathrm{\Delta }kV},\hfill \end{array}\end{array}$$(8)

where V = L³ is the volume and Δk is the k-bin width, which we set to the size of the fundamental mode k_f = 2π/L.

We fit the parameters using a Markov chain Monte Carlo (MCMC) algorithm with the EMCEE library from the cross-power spectrum up to k_max = 0.08 h Mpc⁻¹. We validate our methodology in Appendix A.

4. Results

4.1. Halo density bins

The linear bias of haloes grows monotonically with their local density δ_h, R. This is shown for the UNIT simulation (f_NL = 0) haloes with at least 20 particles (M_h > 2 × 10¹⁰) in Fig. 2. For this measurement, we fit the linear bias (b₁) and shot-noise parameter (A_sn), while fixing f_NL = 0. Figure 2 also shows the agreement between the results using the auto- and cross-power spectra. These results are in line with Castorina et al. (2018), and allow us to find the b₁ ≈ 0 tracers.

Fig. 2. Measured linear halo bias b1 (68% c.l.) as a function of halo overdensity δh, R in spheres of R = 8 h−1 Mpc for the UNIT simulation (fNL = 0). Similar results are found for fNL = 100.

For the PNG-UNIT simulation (f_NL = 100) haloes, we additionally fit for b_ϕ, while fixing f_NL = 100. The b₁(δ_h, R) relation measured from the cross-correlation is nearly identical to the f_NL = 0 case. However, for the auto-correlation, the fits near b₁ ≈ 0 are more numerically unstable, with b₁ absorbing the shot noise. For simplicity, we focus here on the cross-correlation results.

At linear scales, δ_h = b₁δ_M. Hence, one would expect that we measure something compatible with b₁ = 0 whenever δ_h = 0. However, this is not exactly what we see in Fig. 2, where the curve does not cross {b₁ = 0, δ_h, R = 0}. This arises because R = 8 h⁻¹ Mpc is mildly non-linear; we verified that increasing R does not change our conclusions, while approaching the {b₁ = 0, δ_h, R = 0} crossing.

Another interesting feature is that the best fit for A_sn typically ranges in 20 − 40, with A_sn = 27 for the bin closest to b₁ = 0. Castorina et al. (2018) already showed that these tracers were more stochastic than Poisson, but the factor found here is even higher, which might be due to the increased number density (and reduced baseline shot noise). This contribution of a factor of ∼30 in the increase in shot noise already suggests that using these tracers as optimal for PNG will be challenging, according to Eq. (6). We observe that A_sn becomes higher for larger R, suggesting this effect comes from the exclusion of tracers.

We now move to the central idea of the paper using the PNG-UNITsim haloes with f_NL = 100. In the top panel of Fig. 3 we represent the fitted b_ϕ against the linear bias for all haloes when classified in bins of δ_t, R as explained in Sect. 3.1.

Fig. 3. PNG-response (bϕ) as a function of the measured bias (b1) for the tracers binned in their local overdensity. Top: Results for all haloes (purple band) together. The black dot-dashed line shows the universality relation. We find the PNG response well below the universality relation and compatible with no response (bϕ ∼ 0). We also show (as markers with errorbars) the results using mass bins as parent catalogues. Bottom: Same results, but for the LRGs and QSOs from PNG-UNITsim-XL.

We find that b_ϕ strongly deviates from the universality relation (Eq. (5)) marked by the dot-dashed line in Fig. 3. Moreover, we find that b_ϕ is compatible with zero for most of the δ_h, R bins. In particular, the bins close to b₁ = 0 are clearly vanishing.

4.2. Mass and halo density bins

A more realistic scenario is that you have a galaxy sample that is typically associated with a mass range. In Castorina et al. (2018), they test three mass bins as parent catalogues, similar to the three bins used here and defined in Sect. 3.1. We also show the b_ϕ(b₁) results for the mass bins in the top panel of Fig. 3. Remarkably, for the lower mass bin, we find a very flat b_ϕ(b₁)∼0 curve (as in Sect. 4.1). On the other hand, the medium- and high-mass bins show a mild PNG response for |b₁|> 1, but are still typically below the universality relation. The excess of shot noise becomes more moderate for mass bins; in increasing order, we find A_sn= 13, 7, and 2.

4.3. Galaxies in density bins

We now use a more realistic scenario, where we sample the halo catalogue with an HOD, representing DESI QSOs and LRGs in the PNG-UNITsim-XL simulations. While we show the results for the high-density samples (see Sect. 2), these are consistent with those from the realistic scenario, but less noisy.

The results are shown at the bottom of Fig. 3. In this case, the results follow a trend that is more similar to the universality relation. Nevertheless, the points near b₁ = 0 are still very compatible with |b_ϕ|  = 0 and clearly below the universality relation. In this case, the excess of shot noise is more moderated with A_sn = 1.26 for LRGs and A_sn = 1.21 for QSOs.

4.4. Density splits for all parent catalogues

Finally, we use an approach that is matched to the original Castorina et al. (2018) proposal: splitting each of the parent haloes in just two density bins, according to a ${\delta }_{r,R}^{\max }$. For each of the parent catalogues, we identify the ${\delta }_{t,R}^{\max }$ that makes b₁ closest to zero, and plot the values of b_ϕ for each of them in Fig. 4, and compare the results to the equivalent results from the δ_t, R bins. We find all of these cases to also be compatible with b_ϕ = 0.

Fig. 4. PNG response (bϕ) and linear bias (b1) measurements for the closest case to b1 = 0 for the density bins (blue) and density split (red) samples. We consider different parent catalogues: all haloes, haloes in mass bins, and HODs (Sect. 3.1). For all cases, the PNG response is clearly smaller than predicted by universality (dot-dashed line) and compatible with bϕ = 0.

5. Discussion and conclusions

In this work we have used state-of-the-art PNG simulations from the PNG-UNITsim suite to evaluate the efficiency of unclustered tracers (b₁ = 0) for the study of local-PNGs as postulated by Castorina et al. (2018). We classified a parent tracer catalogue (all haloes, haloes in mass bins, or HOD-based galaxies) by their local density δ_t, R in spheres of R = 8 h⁻¹ Mpc to measure their linear bias (b₁) and the PNG response (b_ϕ) from the cross-power spectrum of the δ_t, R-selected tracers and the matter field.

When considering all haloes at once or low-mass haloes (10.5 < log M < 11), this method leads to b_ϕ(b₁)≈0 for all δ_t, R bins (except below b₁ ∼ −1.5). This implies that these tracers will remain unclustered after the inclusion of PNG. More graphically, they will follow the black line in Fig. 1 and not the green line. In the case of the medium-mass and higher mass bins (11 < log M < 12, log M > 12), we recover a slight hint for b_ϕ ≠ 0, but only for |b₁|> 1. When considering a parent sample of LRGs and QSOs from a HOD model, matching the galaxy bias of DESI, we recover a b_ϕ(b₁) curve closer to universality (Eq. (5)) for b₁ > 1, but that still flattens to b_ϕ = 0 for b₁ ≈ 0.

In summary, all unclustered tracers (b₁ ∼ 0) found here remain unclustered with the inclusion of local-PNGs (Fig. 4), strongly deviating from the universality relation. In addition to the density bins, we also considered splitting the sample into two subsamples given a ${\delta }_{t,R}^{\max }$, and found qualitatively similar results. We also verified that the results mentioned above are robust to scale-cuts, the usage of auto- and cross-power spectra, and the variation in the smoothing scale R (from 4 to 30 h⁻¹ Mpc). A similar trend is also found in Appendix B when classifying haloes in matter density instead of halo or galaxy density. We note that for b_ϕ = 0 the Fisher information in f_NL (Eq. (6)) will vanish, leaving f_NL unconstrained. In addition, we find that these samples tend to have super-Poissonian shot noise, sometimes by a factor of ∼30, also complicating the optimality.

This b_ϕ(b₁ ≈ 0)≈0 relation reduces the effectiveness of constraining PNGs. As an example, if we do a Fisher forecast for the LRGs and QSOs modelled in PNG-UNITsim-XL for the DESI volume and number density, we obtain, respectively, σ_fNL = 4.1 and σ_fNL = 3.8. When considering their δ_h, R bin closest to b₁ = 0, even if neglecting shot noise, we obtain respectively, σ_fNL = 133 and σ_fNL = 102. Additionally, the b_ϕ(b₁ = 0) = 0 relation would also hinder a multi-tracer analysis, where σ(f_NL)⁻² ∝ |b₁^Ab_ϕ^B − b₁^Bb_ϕ^A| (Barreira & Krause 2023).

The lack of a PNG response appears in tension with the results reported by Castorina et al. (2018), where they found a PNG response of up to b_ϕ ≈ 2.5 for a b₁ = 0 tracer from their log M > 12.5 parent catalogue². This difference could possibly arise because that work used separate-universe simulations (with f_NL = 0) to estimate b_ϕ as a derivative of the number of tracers with respect to the amplitude of primordial fluctuations (A_s). Whereas the separate-universe approach has been tested in halo mass bins (Biagetti et al. 2017), the selection based on the local halo density inherently brings higher order terms (in density and bias) to the selection that could in turn have cosmology dependence (including A_s). This possibility is supported by Morawetz et al. (2025) who found that δ_h, R-selected tracers (from a parent catalogue of log M > 13.5 haloes) have a PNG response that deviates from the separate-universe prediction (see their Fig. 9).

Exploring the PNG signal from alternative methods and samples is a powerful way to maximise the information from current and future missions. This work remarks on the importance of developing in parallel realistic suite of simulations to validate and improve the modelling of these less standard statistics.

Acknowledgments

We thank discussions with E. Paillas, J. L. Bernal, R. Angulo and I. Sevilla. SA has been funded by MCIN/AEI/10.13039/501100011033 and FSE+ (Europe) under project PID2024-156844NA-C22 and the RYC2022-037311-I fellowship. VGP is supported by the CNS2024-154242 grant. Contributions: CM led the analysis and figure making. SA led the conceptual design of the project and writing of the paper. AGA contributed to the analysis and interpretation. AA led a first version of the analysis. AGA, SA and VGP led the simulation construction. JMR constructed the HOD catalogues.

References

Appendix A: The lack of PNG signal

The flat curve for b_ϕ(b₁) = 0 at the top of Fig. 3 can be somewhat surprising at first. To discard errors in the analysis, we verified that for simple halo-mass bins we recover a b_ϕ(b₁) law which is very close to the universality relation (Eq. 5). In fact, we are able to reproduce the results from Adame et al. (2024) with the analysis tools used for this letter.

Another way to verify the tools used here is to select a bin in mass with bias b₁ ∼ 3, and compare it to one of the bins in δ_h, R with a similar bias value (using all haloes as parent catalogue). This is shown in Fig. A.1, where we observe a very clear PNG effect at low k for the mass bins, whereas we do not observe any PNG response for the density bins.

Fig. A.1. Power spectrum of a bin with b1 ∼ 3.2 selected by halo mass (Mh, left) and by local halo overdensity (δh, R, right). Whereas the effect of introducing PNG (fNL = 100) has a strong effect on mass-selected bins, the δh, R bins show near identical clustering for fNL = 0 and fNL = 100.

Appendix B: Classifying by dark matter overdensity

In the main body, we focused on using biased tracers (haloes or galaxies) for the estimation of local density, δ_t, R. This is motivated by assuming this is what we can do with galaxy surveys. One could think of other ways of estimating local density, such as lensing (either from CMB or from galaxy shear), which will be a more direct probe of the matter field, although typically projected to a 2D sky over a long radial kernel.

For this reason, in this appendix we test an ideal scenario where we would have access to the matter field, δ_m, R and we can classify haloes according to this density. From the resulting classification, we can compute the linear (b₁) and PNG (b_ϕ) bias as before, obtaining Fig. B.1. We find that there is a strong deviation from universality, with negative b_ϕ for 0.5 < b₁ < 1.5. Nevertheless, what remains unchanged with respect to other tracers is that for b₁ = 0 we still find a clear b_ϕ = 0.

Fig. B.1. PNG response bϕ as a function of linear bias b1 for haloes classified as a function of local matter density δm, R using all PNG-UNITsim (fNL = 100) haloes above 20 particles.

All Figures

	Fig. 1. Theoretical auto-power spectra for matter (b = 1, fNL = 0), a biased tracer in a Gaussian Universe (bg = 2, fNL = 0), a biased universal tracer with PNG (bg = 2, fNL = 100, p = 1), an unclustered but universal tracer (bg = 0, fNL = 100, p = 1), and an unclustered and non-PNG-responsive tracer (bg = 0, bϕ = 0).
In the text

	Fig. 2. Measured linear halo bias b1 (68% c.l.) as a function of halo overdensity δh, R in spheres of R = 8 h−1 Mpc for the UNIT simulation (fNL = 0). Similar results are found for fNL = 100.
In the text

Fig. 3. PNG-response (bϕ) as a function of the measured bias (b1) for the tracers binned in their local overdensity. Top: Results for all haloes (purple band) together. The black dot-dashed line shows the universality relation. We find the PNG response well below the universality relation and compatible with no response (bϕ ∼ 0). We also show (as markers with errorbars) the results using mass bins as parent catalogues. Bottom: Same results, but for the LRGs and QSOs from PNG-UNITsim-XL.

In the text

	Fig. 4. PNG response (bϕ) and linear bias (b1) measurements for the closest case to b1 = 0 for the density bins (blue) and density split (red) samples. We consider different parent catalogues: all haloes, haloes in mass bins, and HODs (Sect. 3.1). For all cases, the PNG response is clearly smaller than predicted by universality (dot-dashed line) and compatible with bϕ = 0.
In the text

	Fig. A.1. Power spectrum of a bin with b1 ∼ 3.2 selected by halo mass (Mh, left) and by local halo overdensity (δh, R, right). Whereas the effect of introducing PNG (fNL = 100) has a strong effect on mass-selected bins, the δh, R bins show near identical clustering for fNL = 0 and fNL = 100.
In the text

	Fig. B.1. PNG response bϕ as a function of linear bias b1 for haloes classified as a function of local matter density δm, R using all PNG-UNITsim (fNL = 100) haloes above 20 particles.
In the text