Predictions for local PNG bias in the galaxy power spectrum and bispectrum and the consequences for $f_{\rm NL}$ constraints

We use hydrodynamical separate universe simulations with the IllustrisTNG model to predict the local primordial non-Gaussianity (PNG) bias parameters $b_{\phi}$ and $b_{\phi\delta}$, which enter at leading order in the galaxy power spectrum and bispectrum. This is the first time that $b_{\phi\delta}$ is measured from either gravity-only or galaxy formation simulations. For dark matter halos, the popular assumption of universality overpredicts the $b_{\phi\delta}(b_1)$ relation in the range $1 \lesssim b_1 \lesssim 3$ by up to $\Delta b_{\phi\delta} \sim 3$ ($b_1$ is the linear density bias). The adequacy of the universality relation is worse for the simulated galaxies, with the relations $b_{\phi}(b_1)$ and $b_{\phi\delta}(b_1)$ being generically redshift-dependent and very sensitive to how galaxies are selected (we test total, stellar and black hole mass, black hole mass accretion rate and color). The uncertainties on $b_{\phi}$ and $b_{\phi\delta}$ have a direct, often overlooked impact on the constraints of the local PNG parameter $f_{\rm NL}$, which we study and discuss. For a survey with $V = 100{\rm Gpc}^3/h^3$ at $z=1$, uncertainties $\Delta b_{\phi} \lesssim 1$ and $\Delta b_{\phi\delta} \lesssim 5$ around values close to the fiducial can yield relatively unbiased constraints on $f_{\rm NL}$ using power spectrum and bispectrum data. We also show why priors on galaxy bias are useful even in analyses that fit for products $f_{\rm NL} b_{\phi}$ and $f_{\rm NL} b_{\phi\delta}$. The strategies we discuss to deal with galaxy bias uncertainties can be straightforwardly implemented in existing $f_{\rm NL}$ constraint analyses (we provide fits for some of the bias relations). Our results motivate more works with galaxy formation simulations to refine our understanding of $b_{\phi}$ and $b_{\phi\delta}$ towards improved constraints on $f_{\rm NL}$.


Introduction
One of the main open questions in cosmology today concerns the origin of the primordial density fluctuations generated in the early Universe during the epoch of inflation. The simplest theoretical explanation is that of single-field models, in which the density perturbations are due to quantum fluctuations of a single scalar degree of freedom that rolls slowly down its potential. These models crucially predict density fluctuations that are Gaussian distributed [1][2][3][4][5][6], and as a result, the detection of any deviation from Gaussianity would immediately rule them out and open the door for more elaborate, multifield models [7,8]. Local-type primordial non-Gaussianity (PNG) is the most popular way to describe departures from perfectly Gaussian-distributed fluctuations. In this case, the primordial gravitational (Bardeen) potential φ(x) is expanded as [9] φ where φ G is a Gaussian distributed random field and · · · indicates ensemble averaging. The parameter f nl is a constant that quantifies the leading-order departure from Gaussianity (this expansion can continue to include third and higher powers of φ G ). Analyses of the three-point function of the cosmic microwave background (CMB) radiation measured by the Planck satelllite currently constrain f nl = −0.9 ± 5.1 (1σ) [10], but next-generation large-scale structure surveys hope to improve upon this bound. As galaxies form out of the primordial fluctuations, the statistics of their distribution can in principle be used to probe PNG. In cosmologies with local PNG, the deterministic galaxy density contrast δ g can be written to leading order (LO) as [11][12][13] δ g (x, z) where δ m is the total matter density contrast and q is the initial Lagrangian coordinate associated with the final Eulerian coordinate x. In keeping with the effective field theory approach to galaxy statistics and the galaxy bias expansion (see Ref. [14] for a comprehensive review), the wavelength of the perturbations δ m (x, z) and φ(q) is implicitly assumed to be much larger than the spatial scales over which galaxy formation takes place. The coefficients b 1 and b φ are called galaxy bias parameters and they describe, respectively, the linear responses of the galaxy number density to total mass and primordial potential perturbations with f nl = 0. The bias parameters effectively encode all of the complicated dependence of galaxy formation on the large-scale environment, are functions not only of redshift but also of the properties of the galaxies considered, and are thus extremely challenging to predict theoretically. The bias parameters are therefore normally fitted alongside the cosmological ones in inference analyses using galaxy clustering data, but as degeneracies between the two sets of parameters arise, one often finds it necessary/beneficial to take theoretical priors on bias into account in order to tighten the constraints on cosmology. These degeneracies are especially critical in f nl constraints. For example, Ref. [15] showed that the galaxy power spectrum P gg (k) (the Fourier transform of the two-point correlation function) in cosmologies with local PNG gets a contribution ∝ b φ f nl /k 2 , that becomes important on large scales and that can be used to constrain f nl . Given the perfect degeneracy between f nl and b φ , however, it is simply impossible to use P gg to constrain f nl , unless some theoretical prior on b φ is assumed 1 .
In the large-scale structure literature, the most popular way to circumvent this problem in f nl constraints is by relating b φ to b 1 using the so-called universality relation [12,[15][16][17] which follows from assuming that the halo mass function is universal, and δ c = 1.686 is the threshold overdensity for spherical collapse. In this way, the contribution scales as ∝ (b 1 − 1) f nl /k 2 , and since b 1 can in principle be constrained using the smaller-scale part of the power spectrum (where f nl contributes weakly), it then becomes possible to constrain f nl . The universality relation is adopted by almost all existing galaxy data constraints on f nl [16,[18][19][20][21][22][23][24][25] (the current tightest bound is f nl = −12 ± 21 (1σ) [25]), as well as in forecast studies [7,[26][27][28][29][30][31][32][33][34][35][36] for next-generation surveys. Despite its widespread adoption, there is however no reason to expect the universality relation to hold for real-life galaxy samples, and in fact, studies using N -body simulations have been indicating this to be the case already. For example, gravity-only simulations have shown that the b φ (b 1 ) relation of dark matter halos underpredicts slightly the universality relation [37][38][39][40][41], Ref. [42] showed using galaxy formation simulations that the b φ (b 1 ) relation of stellar-mass selected galaxies overpredicts the universality relation, and Refs. [16,43] discussed how the b φ (b 1 ) relation of recently-merged objects differs also from the universality prediction; we will return to these findings below when we reproduce some with our numerical results. The impact that uncertainties on the b φ (b 1 ) relation have on the resulting f nl constraints is only now beginning to be explored [44,45], but it is crucial that these investigations are made robust and mature since they will ultimately directly impact the final constraints on f nl . Beyond leading-order, the contributions at next-to-leading-order (NLO) are 2 δ g (x, z) where K ij = ∂ i ∂ j /∇ 2 − δ ij /3 δ m is a long-wavelength tidal field, and b 2 , b K 2 and b φδ are additional galaxy bias parameters. These terms contribute to the galaxy bispectrum B ggg (k 1 , k 2 , k 3 ) (the Fourier transform of the three-point correlation function), which is also a very well-known probe of f nl [46][47][48][49][50], and which may prove crucial to combine with the galaxy power spectrum if next-generation galaxy surveys are to improve over the current CMB constraints. For halos, the b 2 (b 1 ) relation is very well understood [51] and works have been progressively improving our understanding of b K 2 (b 1 ) as well [52][53][54][55][56][57][58][59]; the corresponding relations for simulated galaxies have also began to be recently studied [58]. In contrast, the b φδ parameter has never been the focus of dedicated simulation work. Assuming 1 Formally, next-to-leading-order contributions to the galaxy power spectrum (e.g. 1-loop terms and beyond) can break the degeneracy, but they are unimportant given the current observationally allowed range for fnl. 2 There is an additional contribution δg ⊃ b φ 2 f 2 nl φ 2 that is quadratic in fnl, but which is unimportant given the current observational bounds, and so we skip writing it explicitly. Should this contribution have been important, then all our motivations to study b φδ in this paper would apply equally to b φ 2 . universality of the halo mass function, one finds [12,14] (1.5) but unlike the case for b φ , the performance of this relation has never been checked, even for the simpler case of dark matter halos in gravity-only simulations. There are currently no real-data constraints on f nl using the galaxy bispectrum, but existing forecast studies effectively always assume the validity of this relation. Although it is formally possible to constrain f nl using the bispectrum without any prior on b φδ , we will see below that adopting priors based on a b φδ (b 1 ) relation is necessary for competitive constraints. This motivates the main goal of this paper, which is to use separate universe N -body simulations to obtain predictions for the bias parameters b φ and b φδ ; to the best of our knowledge, this is the first time that b φδ is estimated from dedicated simulation data. We will do so not only for the case of dark matter halos in gravity-only simulations, but also for simulated galaxies in hydrodynamical simulations with the IllustrisTNG galaxy formation model. One of our main new results is that the universality relation of Eq. (1.5) is not a perfect description for either one of these large-scale structure tracers, with the size of the departures varying depending on the galaxy selection criteria adopted. Recently, Ref. [58] showed that the b 2 (b 1 ) and b K 2 (b 1 ) relations of halos are broadly preserved for galaxies selected by a variety of criteria (total mass, stellar mass, color and black hole mass accretion rate), which suggests that priors around these relations may be adopted relatively safely in real-data analyses. In contrast, our results below will show that the b φ (b 1 ) and b φδ (b 1 ) relations are markedly more sensitive to the selection criteria, which makes the design of theoretical priors more challenging. In a second part of this paper, we build upon the idealized forecast study of Ref. [44] to study the impact that uncertainties on b φ and b φδ can have on f nl constraints. We will see that these bias parameters can have a marked impact on the resulting f nl bounds. One of our main takeaway messages is that more simulation work is needed to our current level The rest of this paper is organized as follows. In Sec. 2 we describe the simulation data we use in this work, as well as our methods to estimate the galaxy bias parameters. Section 3 contains our main numerical results on galaxy bias, where we show first the b φδ parameter measured in gravity-only simulations, and then compare the sensitivity of the b φ (b 1 ) and b φδ (b 1 ) relations to different galaxy selection criteria in hydrodynamical simulations. In Sec. 4 we show and discuss the results of a simple forecast study for a fictitious survey aiming to illustrate the impact that bias uncertainties have on f nl constraints. We summarize and conclude in Sec. 5. In App. A, we show estimates of the bias parameter b 2 using the same method we use to estimate b φδ . In App. B, we collect the theory model expressions for the galaxy power spectrum and bispectrum that we use in our forecast study.

Methodology
In this section we describe the estimation of the two leading-order local PNG galaxy bias parameters b φ and b φδ , as well as the linear bias parameter b 1 , which recall, contribute to the bias expansion as As explained below, we estimate (i) b 1 using the large-scale limit of the galaxy-matter cross-power spectrum; (ii) b φ as the response of the galaxy number density to long-wavelength primordial gravitational potentials with f nl ; and (iii) b φδ using the response of b 1 to long-wavelength primordial gravitational potentials with f nl .

Simulation data specifications
The simulations we use in this work have been presented previously in Ref. [42], and they were run using the moving-mesh gravity+hydrodynamical N -body code AREPO [60,61] with the IllustrisTNG model of galaxy formation [62][63][64]. This model is an improved version of its precursor Illustris [65,66], and it is characterized by sub-grid prescriptions for the physics of gas cooling, star formation, stellar feedback, chemical enrichment, and black hole growth/feedback, which were calibrated to broadly reproduce a number of observations such as the galaxy stellar mass function at low redshift, the star formation rate history, galaxy sizes and the gas mass fractions of galaxies and galaxy groups (see Refs. [67][68][69][70][71][72][73] for the first results with IllustrisTNG). The initial conditions were generated at z i = 127 with N-GenIC code [74] using the Zel'dovich approximation, and a linear matter power spectrum calculated using the CAMB code [75]. The galaxy formation simulations were run on a cubic box with size L box = 205 Mpc/h, containing N p = 1250 3 dark matter mass elements and N p = 1250 3 initial gas elements. At this resolution, we have both full-physics hydrodynamical simulations, as well as gravity-only counterparts (we label these as "Gravity"); we follow the standard IllustrisTNG nomenclature and refer to this resolution as TNG300-2. In addition, we consider also a set of gravityonly simulations run also with AREPO with N p = 1250 3 dark matter mass elements, but on a bigger simulation box size L box = 560 Mpc/h ≈ 800 Mpc. The cosmological parameters of our Fiducial cosmology are: mean baryon density today Ω b0 = 0.0486, mean total matter density today Ω m0 = 0.3089, mean dark energy density today Ω Λ0 = 0.6911, dimensionless Hubble rate h = 0.6774, primordial scalar spectral index n s = 0.967, and primordial scalar power spectrum amplitude A s = 2.068 × 10 −9 (at k pivot = 0.05/Mpc, corresponding to σ 8 (z = 0) = 0.816). In order to measure the galaxy bias parameters b φ and b φδ using the separate universe method (see below), we consider also two additional cosmologies, dubbed HighA s and LowA s , which differ from the Fiducial only in the value of A s → A s [1 + δA s ], where δA s = +0.05 for HighA s and δA s = −0.05 for LowA s .
We will show measurements of the bias parameters for both halos and subhalos/galaxies. The halos are identified with a Friends-of-Friends (FoF) algorithm run on the dark matter elements with linking length b = 0.2 times the mean interparticle distance. In turn, the subhalos correspond to the gravitationally bound structures found by the SUBFIND algorithm [76] inside each halo. In the hydrodynamical simulations, we refer to the subhalos that contain any mass in stars (M * > 0) as galaxies, and we do not explicitly distinguish between main and satellite galaxies. We also only consider objects with at least 100 member star particles to ensure we deal with objects that are sufficiently well resolved in the simulations. Throughout, we will show results for objects selected by their total mass M t , stellar mass M * , black hole mass M BH , black hole mass accretion rateṀ BH , and (dust-uncorrected) g − r color. When quoting the value of these quantities for a given object (halo or subhalo), we always consider the summed contribution of all member particles to that quantity.
There are two points worth emphasizing about our numerical setup. One is that the TNG300-2 resolution is below that at which the IllustrisTNG model was callibrated at, and as a result, the predictions of our simulations are not in as good agreement with the above-mentioned observations. The galaxy bias predictions are however expected to be less affected by numerical resolution compared to quantities like the galaxy number density itself. In fact, Ref. [42] has shown that the b 1 and b φ predictions at TNG300-2 resolution are in very good agreement with those obtained at a higher resolution (N p = 2 × 1250 3 , L box = 75Mpc/h; called TNG100-1.5 there) that is closer to the nominal IllustrisTNG one (N p = 2 × 1820 3 , L box = 75Mpc/h). The second point is that (as in past works with separate universe simulations of galaxy formation [42,77,78]) we keep the parameters of the IllustrisTNG model fixed when we adjust the parameter A s in the HighA s and LowA s cosmologies. This is the appropriate choice in order to interpret our galaxy bias measurements as predictions of the IllustrisTNG model, i.e., the galaxy bias parameters are the response of galaxy formation to long-wavelength perturbations, at fixed galaxy formation prescription (that of IllustrisTNG).

The linear density and local PNG bias parameters b 1 and b φ
We estimate the linear bias parameter b 1 in the Fiducial simulations using the large-scale limit of the ratio of the galaxy-matter cross-power spectrum P gm (k, z) and matter power spectrum P mm (k, z) which follows from Eq. (2.1) for f nl = 0 (note we dropped the redshift z from the arguments to lighten the notation). For the relatively small volume of the TNG300-2 simulations, the scale-dependence of this ratio can still be nonnegligible on the largest scales probed. Thus, rather than simply fitting a constant to it, we account for the leading-order scale-dependence by fitting instead for b 1 + Ak 2 , and take the constant coefficient as our estimate of the bias parameter. In practice, we use all modes with k < 0.15h/Mpc, and we have checked that the values of b 1 estimated from the TNG300-2 box agree with those from the L box ≈ 800 Mpc box that is less affected by this systematic. Our error bars on b 1 are the error estimate from the least-squares fitting procedure.
On the other hand, we estimate the linear local PNG bias parameter b φ from the definition where the first equality follows from Eq. (2.1) and the second equality from the separate universe (or peak-background split) equivalence between (i) structure formation inside a long-wavelength potential perturbation f nl φ and (ii) structure formation without the perturbation, but with a modified value of the primordial scalar power spectrum amplitude A s [15,16]. In particular, it is possible to show (see e.g. Sec. 7.1.2. of Ref. [14]) that if φ L is the amplitude of a long-wavelength potential perturbation, then galaxies forming inside it form as they would form in a cosmology with A s rescaled as  5) and where N g represents the number of objects in some selection variable bin (total mass, stellar mass, black hole mass, etc.) and the superscripts indicate in which cosmology the number of galaxies is counted. 3 The choice of |δA s | = 0.05 is motivated by the compromise between having a sizeable and measurable impact of the change in A s on the galaxy abundance, while keeping negligible higher-order corrections to the first-order finite-difference result. Further, having just a single realization of the initial conditions for each cosmology/resolution, it is not possible to estimate our measurement errors in a statistical ensemble sense. As a compromise, we take the difference between the values of b HighAs φ and b LowAs φ (which should be the same up to numerical noise) as our estimate of the error on b φ ; see Ref. [42] for a discussion of why this yields trustworthy error estimates.
Using the same simulations and methodology, Refs. [42,82] presented an indepth study of the total-and stellar-mass dependence of b φ . Here, we will reproduce some of these past results (while showing also additional ones for other galaxy selection criteria) to better compare with the results for the second-order b φδ parameter.

The second-order local PNG bias parameter b φδ
The same separate universe simulations can be used to estimate the bias parameter b φδ via the response of P gm (k) to long-wavelength primordial potential perturbations f nl φ. Concretely, in all our cosmologies, the large-scale galaxy-matter cross-power spectrum is described by P gm = b 1 P mm . Defining its linear local PNG response function as R φ,gm = ∂lnP gm /∂(f nl φ), then it follows that (the power spectrum response functions can be defined in analogy to the galaxy bias parameters by treating the local power spectrum as a biased tracer [83]) where R φ,mm = ∂lnP mm /∂(f nl φ). The first term on the right can be worked out by using that to yield where we have used also Eq. (2.3). The power spectrum response functions R φ,gm and R φ,mm can be measured straightforwardly using the separate universe simulations analogously to how b φ is estimated in Eqs. (2.4) and (2.5). For example, using the Fiducial and HighA s cosmologies, we would estimate the response as R HighAs φ,gm = 4 P HighAs gm /P Fiducial gm − 1 /δA s . Thus, given the estimates of b 1 and b φ described in the last subsection, as well as measurements of the power spectrum responses R φ,gm and R φ,mm on large scales, we can use Eq. (2.8) to fit for b φδ . The linear power spectrum response functions effectively describe the mode-coupling structure of squeezed-limit bispectra [83][84][85], which are sensitive to second-order bias parameters to leading order, and is what allows us to fit for b φδ (as well as other second-order bias parameters; cf. App. A). The parameter b φδ can also be estimated from the response of the galaxy power spectrum P gg = b 2 1 P mm +1/n g , although in this case the contribution from shot noise lowers the signal-to-noise unnecessarily; further, it can be straightforwardly shown that the response of this galaxy power spectrum model to f nl φ agrees with the expressions shown in Ref. [86] obtained with squeezed-limit bispectra.
The results on b φδ that we will show below were obtained, however, with a different (yet related) strategy. Rather than using the separate universe simulations to differentiate P gm and P mm , we first fit instead for b 1 in the three cosmologies separately using Eq. (2.2), and then we differentiate b 1 using finite differences. With this estimate, together with the estimates of b 1 and b φ , b φδ is then simply given by We have explicitly checked that both strategies above give consistent results, but the latter yielded slightly better signal-to-noise estimates, which is why we adopt it as the default. Concretely, we evaluate the response of b 1 as and we estimate its error analogously to as for b φ . The final uncertainty on b φδ is worked out by standard propagation of uncertainty in Eq. (2.9) assuming uncorrelated errors, which is conservative since b φ , b 1 and its response are measured from the same simulations and so their uncertainties are correlated. In App. A we validate this strategy to estimate b φδ by applying it (with the appropriate using the values of b1 measured for the same halos, with b φ = 2δc(b1 − 1) and b2(b1) given by Eq. (3.1).
modifications described there) to estimate the second-order bias parameter b 2 , whose values can be compared with known results in the literature. We note before proceeding that the b φδ parameter can also be estimated by fitting perturbation theory models to the bispectrum measured from simulations with local PNG initial conditions [45]. This requires however very large simulation volumes (in fact still currently out of reach for selfconsistent galaxy formation simulations) in order to measure the bispectrum precisely on large-scales where the effects of f nl dominate. At fixed volume, separate universe simulations offer thus the ideal method to study the local PNG bias parameters. Within the separate universe approach, there is an alternative way to estimate b φδ using separate universe simulations that incorporate simultaneously the effects of total mass and primordial gravitational potential perturbations. This would require however additional simulations for higher amplitudes of the long-wavelength modes to be sensitive to the second-order terms in the galaxy bias expansion, similarly to how Ref. [51] uses separate universe simulations to estimate higher-order bias parameters such as b 2 and b 3 .

Galaxy bias results
In this section we show and discuss our main numerical results on the galaxy bias parameters b φ and b φδ . We begin with the b φδ parameter measured for halos and subhalos in the gravity-only simulations, and then compare the results of the b φ (b 1 ) and b φδ (b 1 ) relations for a number of different galaxy samples in the IllustrisTNG simulations; we show results for galaxies selected by total mass M t , stellar mass M * , black hole mass M BH , black hole mass accretion rateṀ BH and (g − r) color.

The b φδ parameter in gravity-only simulations
The total-mass dependence of b φδ is shown in Fig. 1 for dark matter halos at different redshifts, and for the gravity-only TNG300-2 and L box ≈ 800 Mpc simulations, as labeled. The result is compared and b2(b1) given by Eq. (3.1). The dotted line shows the same but with b φ = 0.85 × 2δc(b1 − 1) to account for the departures from universality on this relation. The dashed line shows the outcome of Eq. (2.9) obtained using the formulae from Tinker et al [87,88].
to the universality prediction of Eq. (1.5) (shown in black), evaluated using the values of b 1 estimated for the same mass bins, together with the universality relation for the linear local PNG parameter b φ = 2δ c (b 1 − 1) and b 2 (b 1 ) given by the fit obtained for halos with separate universe simulations in Ref. [51]: The figure shows that the b φδ values measured in the simulations depart from the universality expectation, with the difference being both redshift-and mass-dependent. Concretely, for M t 10 12 M /h, there is a slight trend for b φδ to overpredict (less negative) the universality prediction at z ≤ 0.5, but to underpredict it at z ≥ 2 (more negative). On the other hand, for M t 10 12 M /h, the universality relation prediction is visibly above the measured b φδ at z ≤ 1, but the two are consistent within the errors at z ≥ 2. As one would expect for its smaller volume, the TNG300-2 results appear noisier compared to L box ≈ 800 Mpc, but it is nonetheless possible to discern a good overall agreement between the two resolutions on the mass scales where they overlap. Figure 2 shows the b φδ (b 1 ) relation for total-mass selected halos and subhalos in the gravity-only TNG300-2 and L box ≈ 800 Mpc simulations at different redshifts, as labeled. The relation is seen to depend only weakly on redshift; this can be better appreciated in the higher signal-to-noise results from the L box ≈ 800 Mpc box, but the TNG300-2 results display a consistent picture. Further, in the range 1 b 1 3, the b φδ values of the simulations are systematically below (more negative/less positive) the universality prediction shown by the solid black line; the largest difference occurs at b 1 ≈ 2 and is ∆b φδ ∼ 3. For b 1 0.8, one can also discern a trend for the measured b φδ to overpredict the universality relation, which is a manifestation of the same trend shown in Fig. 1 at low redshift and lower masses. Note also that the results for halos and subhalos are consistent within the precision of our measurements, although some differences are in general to be expected (even if small) since the relation between b φδ (b 1 ) is nonlinear (see Sec. 3.3 of Ref. [58] for a discussion).
The dotted lines in Fig. 2 show the prediction of a variant of the universality relation of Eq. (1.5), in which instead of using the universality relation for b φ (cf. Eq. (1.3)), one replaces it by b φ (b 1 ) = 0.85 × 2δ c (b 1 − 1), which offers a more adequate approximation to the b φ (b 1 ) relation in gravity-only simulations [37][38][39][40][41][42] (see the dotted line in the top left panel of Fig. 3 below). This variant of the universality relation with a modified b φ (b 1 ) relation does get slightly closer to the simulation measurements over the range 1 b 1 3, but not sufficiently to bring the two results into agreement. In other words, the breakdown of the universality relation for the second-order parameter b φδ cannot be attributed solely to the breakdown of the universality relation for the linear parameter b φ . The dashed lines in Fig. 2 show the outcome of Eq. (2.9) evaluated with the halo mass function and halo bias b 1 formulae from Tinker et al [87,88]. Concretely, the ∂lnb 1 /∂(f nl φ) ≡ 4∂lnb 1 /∂(δA s ) term is evaluated by finite-differencing the b 1 formula w.r.t. A s , and b φ is obtained by finite-differencing the halo mass function also w.r.t A s as in Eq. 1)). This semi-analytical calculation agrees extremelly well with our numerical estimates, which represents a good cross-check of both results. Finally, we have also found the following best-fitting quadratic polynomial to the dark matter halo results (fitted up to b 1 < 4 using all redshifts and both simulation boxes): but which we skip showing to avoid crowing the figure.

The b φ and b φδ parameters in galaxy formation simulations
We turn our attention now to the b φ and b φδ parameters measured from the hydrodynamical, fullphysics IllustrisTNG runs with the TNG300-2 box; in this subsection we focus only on the results for galaxies (i.e., subhalos with mass in stars). Figure 3 shows the b φ (b 1 ) and b φδ (b 1 ) relations at different redshifts for galaxies selected by their total mass M t , stellar mass M * , black hole mass M BH , black hole mass accretion rateṀ BH , and dust-uncorrected (g −r) color, as labeled. Concretely, we consider 6 M t bins log-spaced between 5 × 10 10 ; 8 × 10 12 M /h. For M * , M BH andṀ BH we consider 4 bins logspaced between 1 × 10 10 ; 3 × 10 11 M /h, 5 × 10 6 ; 5 × 10 8 M /h and 1 × 10 3 ; 1 × 10 9 M /Gyr, respectively. And for (g − r) we consider 4 bins linearly spaced within [0, 0.8]. In Fig. 4, we show the dependence of the three galaxy bias parameters b 1 , b φ and b φδ on the galaxy selection variables themselves for z = 0, z = 1 and z = 3. For the case of b φ , the total-mass and stellar-mass results in Fig. 3 (top two panels on the left) have been already discussed in detail in Ref. [42] using the same simulations as in this paper. Namely, for M t -selection, the figure shows the well-known result that the b φ (b 1 ) relation of the simulated objects is below the universality expectation for b 1 1.5, and that the modified relation b φ (b 1 ) = 0.85 × 2δ c (b 1 − 1) (shown by the dotted line) provides a more adequate description. For M * -selection, the b φ (b 1 ) relation of the galaxies systematically overpredicts the universality relation, with b φ (b 1 ) = 2δ c (b 1 − 0.55) being a more faithful approximation of the simulation measurements (center of the grey color map). The lower three panels on the left of Fig. 3 display the result for galaxies selected by M BH ,Ṁ BH and (g − r), and they illustrate the strong sensitivity of the b φ (b 1 ) relation to the galaxy selection criterion adopted, as well as the fact that the universality relation continues to be an inadequate description. We highlight for example the case of theṀ BH -selected galaxies, for which the b φ (b 1 ) relation becomes especially steep, and does not even admit a regular, redshift-independent function of b 1 .
The dashed line in theṀ BH panel for b φ shows the prediction of the variant of the universality relation b φ (b 1 ) = 2δ c (b 1 −1.6), which was derived by Ref. [16] for dark matter halos that had undergone a recent major merger. The authors further hypothesized it could be a better description of the b φ (b 1 ) relation of quasars compared to the universality relation, as recent mergers may correlate with strong active galactic nuclei (AGN) activity/luminosity. The reason why it is interesting to make this comparison is because AGN luminosity is also thought to be proportional to the accretion rate of the supermassive black hole. The result of Fig. 3 shows, however, that at least in the IllustrisTNG model, the relation b φ (b 1 ) = 2δ c (b 1 − 1.6) is not a good description of the objects selected byṀ BH (or by proxy, selected by their AGN luminosity). With just a single realization of the initial conditions of the simulations, and for a single galaxy formation model, our results do now allow us yet to conclude decisively on the b φ (b 1 ) relation of real-life quasars and AGN. However, in order to obtain competitive and unbiased constraints on f nl , theṀ BH results depicted in Fig. 3 do strongly motivate  more works with galaxy formation separate universe simulations to determine the precise bias relations for these objects. Note that the tightest constraints on f nl using large-scale structure to date were obtained precisely with quasar samples from eBOSS DR14 [24] and DR16 [25], who assumed b φ (b 1 ) = 2δ c (b 1 − 1.6) in parts of their analysis.
For the case of the b φδ results on the right of Fig. 3, a first point to note concerns the lower signal-to-noise of the measurements compared to b φ . This is not surprising since b φδ is a secondorder bias parameter, and our method to estimate it relies on estimating b 1 using the large-scale limit of the galaxy-matter cross-power spectrum (cf. Sec. 2.3), which can be somewhat uncertain on a L box = 205Mpc/h box; recall from Fig. 2 how the signal-to-noise improves substantially from the TNG300-2 to the bigger L box ≈ 800 Mpc box. There are nonetheless a few trends that one can discern. For example, similarly to the b φ (b 1 ) relation, there is also a visible trend for the b φδ values at fixed b 1 to be higher for stellar-mass relative to total-mass selection, and the b φδ (b 1 ) relation displays similar steep variations forṀ BH -selected galaxies, i.e., b φδ can vary rapidly in a narrow b 1 -interval. Also similarly to b φ , the b φδ (b 1 ) relations do not show evidence of being robust to changes in the galaxy selection criterion, although the poorer signal-to-noise here makes it more difficult to draw decisive conclusions.
The details of the b φ (b 1 ) and b φδ (b 1 ) relations can be understood by inspecting the corresponding dependencies of the bias parameters on the galaxy selection variables in Fig. 4. For example, the shape of the dependence of the bias parameters can vary quite significantly from one galaxy property to another, which is the reason behind the strong sensitivity of the b φ (b 1 ) and b φδ (b 1 ) relations to the galaxy selection criteria seen in Fig. 3. It is interesting to contrast this result for b φ (b 1 ) and b φδ (b 1 ), with the appreciably weaker sensitivity shown in Ref. [58] for the b 2 (b 1 ) and b K 2 (b 1 ) relations using also IllustrisTNG galaxies. As described in Refs. [42,58,82], this can be explained using the halo model and halo occupation distribution formalisms, but we leave a more detailed investigation along these lines to future work.
Finally, the grey color maps in the stellar-mass panels in Fig. 3 describe the shape of Gaussian priors on the b φ (b 1 ) and b φδ (b 1 ) relations that we will use in the next section to study the impact of galaxy bias uncertainties on f nl constraints. Concretely, for b φ and b φδ , our assumed priors are, respectively, where, we stress, the mean relations are expected to describe only the results for M * -selection, and the standard deviations are assumed b 1 -independent for simplicity and chosen to roughly match the uncertainty in our numerical results with IllustrisTNG.

Impact of galaxy bias uncertainties on f nl constraints
In this section we show and discuss the results of an idealized forecast setup for a fictitious survey, to study the impact that uncertainties on b φ and b φδ have on f nl constraints. This analysis effectively extends that of Ref. [44] to include also uncertainties on the b φδ (b 1 ) relation.

Forecast setup
We work with a Gaussian likelihood function with a data vector consisting of hypothetical measurements of the real-space multitracer power spectrum P gg [89,90] and bispectrum B ggg where the superscripts A and B indicate two galaxy samples, e.g.,P AB gg denotes the cross-power spectrum of the two samples. We consider only the bispectrum of sample A for simplicity; a more involved analysis could include also the bispectrum of sample B and the cross-bispectra B AAB ggg and B ABB ggg , but this additional complication is not critical to our discussion. In our theory model, we evaluate the galaxy power spectrum and bispectrum at tree level using the following galaxy bias expansion where in addition to the deterministic terms that appeared already before, we include now also the relevant stochastic contributions , δ , φ . This expansion contains all terms that are needed to selfconsistently derive the leading-order galaxy power spectrum and bispectrum [13]. We display the corresponding expressions for the power spectrum, bispectrum and our treatment of the covariance matrix in App. B.
We consider a galaxy sample at redshift z = 1 covering a volume of V = 100Gpc 3 /h 3 . For the multitracer part of the data vector, we split this galaxy sample into a low-and a high-stellarmass subsamples, with M * ∈ 5 × 10 10 ; 2 × 10 11 M /h for subsample A and M * > 2 × 10 11 M /h for subsample B. For the stellar mass function of the IllustrisTNG simulations this corresponds to the number densitiesn A g = 1.74 × 10 −3 h 3 /Mpc 3 ,n B g = 1.07 × 10 −4 h 3 /Mpc 3 , and linear bias parameters b A 1 = 1.58, b B 1 = 2.37. We evaluate b 2 and b K 2 for these samples using b 2 (b 1 ) = 0.30 − 0.79b 1 + 0.20b 2 1 + 0.12b 3 1 and b K 2 (b 1 ) = 0.66 − 0.57b 1 , which are fits for IllustrisTNG galaxies obtained by Ref. [58]. The fiducial values of b φ and b φδ are given by the mean relations in Eqs. (3.3) and (3.4). We assume Poissonian statistics for the fiducial power spectrum and bispectrum of the noise terms, i.e., P = 1/n g , P δ = b 1 /(2n g ) and B = 1/n 2 g (but note that we sample these in our constraints too; see Ref. [45] for the importance of k-dependent corrections to the shot noise in f nl constraints).
We generate our data vector as a noiseless realization of our theory model for the bias parameters and noise terms listed above, and at a fiducial cosmology that has the same parameters as our Fiducial simulations (cf. Sec. 2.1), except for f nl = 5. The adoption of a noiseless data vector is naturally idealized, but which we note is the most adequate choice for our purpose here to isolate the impact of galaxy bias uncertainties; see also Ref. [45], who takes halo power spectra and bispectra measurements from simulations as the data vector, and so the assessment of the true impact of bias uncertainties gets complicated by the limitations of the theory model to describe the simulation measurements. We consider 42 k-bins between k min = π/V 1/3 s and k max = 0.2 h/Mpc (in intervals of k min up to 0.01h/Mpc and of 10 × k min beyond that).
Similarly to Ref. [44], we show results for two ways to deal with galaxy bias uncertainties: • Parametrization 1: direct priors on b φ and b φδ .
In this case, we sample the following 13 dimensional parameter space: With this parametrization we will study the impact that different priors on b φ (b 1 ) and b φδ (b 1 ) have on f nl constraints.
• Parametrization 2: fit for products f nl b φ , f nl b φδ .
In this case, rather than making assumptions on the bias parameters, we fit directly for the parameter combinations f nl b φ and f nl b φδ as they enter the galaxy power spectrum and bispectrum (cf. App. B). This approach makes it harder to constrain f nl directly, but note it can still be powerful as it can provide evidence for f nl = 0 in a way that is independent of galaxy bias assumptions. This parametrization prevents however direct comparisons of the constraining power of galaxy and CMB data, as well as the combination of galaxy-and CMB-based bounds to obtain a tighter combined bound on f nl . The parameter space is also 13 dimensional: In this part of the analysis, we will also consider cases in which we do not fit for [f nl b A φδ ], and replace it with priors on b φδ (b 1 ) to illustrate the breaking of parameter degeneracies.
The majority of our results are for the full data vector of Eq. (4.1), but we shall also display results for power-spectrum-only analyses (shown in dashed-black and labeled as P gg -only). In this case, in parametrizations 1 and 2 we need to consider only the first 7 and 6 parameters, respectively; we keep A s fixed at the fiducial value in these cases. Except for the b φ and b φδ parameters, we always assume wide uninformative linear priors when sampling the parameter space, which we do using the EMCEE Python implementation [91] of the affine-invariant Markov Chain Monte Carlo (MCMC) sampler in Ref. [92]. We use 32 walkers with a chain convergence criteria that (i) the size of the chain must be 100 times the autocorrelation time and (ii) the latter having varied less than 1% since the last calculation point, which is every few thousand samples. We stress that our idealized forecast setup serves primarily the purpose to provide a simple framework to visualize the impact that galaxy bias uncertainties can have on f nl constraints, and that it is not intended to be representative of any specific current or future survey.

Results from parametrization 1: direct priors on b φ and b φδ
The marginalized 1σ constraints on f nl are shown in the lower panel of Fig. 5 for parametrization 1. The upper panels show the two-dimensional 1σ marginalized constraints on the b A φ − f nl and b A φδ − f nl planes; we do not show the contours for all of the 13 parameters for brevity (the interested reader can find some of these triangle plots in Ref. [44]). The different colors show the result for different assumed priors on the b φ (b 1 ) and b φδ (b 1 ) relations, as labeled; the "Gaussian prior" result in orange corresponds to the use of Eqs. Assuming perfect knowledge of the bias relations, our fictitious galaxy sample would be able to constrain f nl with a 1σ uncertainty of σ fnl = 1.7 (solid black error bar). When we adopt our IllustrisTNG-inspired priors for the bias relations of M * -selected galaxies, which recall are centered around the fiducial relations and assume an uncertainty on b φ (b 1 ) and b φδ (b 1 ) of 1 and 5, respectively (cf. Eqs. (3.3) and (3.4)), then the f nl constraint becomes f nl = 4.1 +2.2 −1.6 (1σ) (orange error bar), i.e., the total 1σ interval increases by about 10% and the mean distribution value gets shifted from the truth by ≈ 0.5σ fnl . This shift is a manifestation of projection effects that can be understood as follows.
In the parts of the parameter space that are close to f nl = 0 (which our chains still explore), the theory predictions become insensitive to both b φ and b φδ , which can take on any value allowed by the assumed priors. The wider the priors, the greater the volume of the parameter space that is near the f nl = 0 direction, and so marginalizing over the poorly constrained b φ and b φδ will progressively center the marginalized constraints around f nl = 0. Indeed, the green and blue error bars show the result for 2× and 10× wider priors, which becomes more biased, as expected. On the other hand, halving the width of the priors (magenta), tightens the allowed range of b φ and b φδ sufficiently to the point where the constraints become effectively the same as the "perfect knowledge" case. A cautionary tale here is that contrary to what one might have naively expected, wide priors on the bias parameters are not conservative and should be interpreted carefuly in light of projection effects like these. 4 In real-life applications, however, the assumed priors will likely not be exactly centered around the bias of the observed galaxies, which introduces additional shifts in the best-fitting f nl and its inferred error bar. This is illustrated in Fig. 6, which shows the same as the lower panel of Fig. 5, but for offset priors on the bias parameters centered around the universality relations, instead of the fiducial relations used to generate the data vector. Assuming zero uncertainty on the wrong bias relations (black error bar), our fictitious survey would constrain f nl = 7.5 ± 2.5 (1σ), i.e., the inferred value would be 1σ away from the fiducial f nl = 5 (the upward shift in the constraint and the increase of the size of the error bar compared to Fig. 5 is as expected since the universality relations underpredict the fiducial bias values). If one is interested only on the overall detection of f NL = 0, then this represents a ≈ 3σ detection, which is the same significance as the constraint f nl = 5.0±1.7 (1σ) in Fig. 5 assuming perfect knowledge of the fiducial bias relations. Thus, if the true value is f nl = 0, then from the point of view of detection significance, the assumptions on the bias parameters are not as critical. Note, however, that different assumptions on the bias relations still directly impact the inferred precision of the constraint σ fnl , which can lead to misinterpretations about the true constraining power of the data. For detection significance purposes, a much cleaner approach is thus that based on parametrization 2 discussed below. Further, as the width of the prior increases in Fig. 6, the constraint progressively approaches f nl = 0 due to the projection effects discussed above. Interestingly, the result in orange shows the case for widths of ∆b φ = 1 and ∆b φδ = 5, for which the projection effects balance the shifts induced by centering the prior around the wrong b φ (b 1 ) and b φδ (b 1 ) relations.
Although not shown, we have also repeated the analyses behind Figs. 5 and 6, but for a fiducial value of f nl = 0, instead of f nl = 5. In this case, for both the cases with priors centered around and offset from the fiducial bias relations, the constraints on f nl are always unbiased, but too wide priors on b φ (b 1 ) and b φδ (b 1 ) eventually also artificially shrink the marginalized error bar due to the projection effects. For the case of assuming perfect knowledge of the bias relations and zero uncertainty around the wrong ones, our fictitious survey would constrain |f nl | < 1.2 (1σ) and |f nl | < 2.0 (1σ), respectively, further illustrating how wrong assumptions about the bias relations can directly impact the apparent constraining power of the data.
We emphasize that the absolute values of the constraints in Figs. 5 and 6 correspond strictly to our fictitious galaxy survey, but the relative impact of the galaxy bias uncertainties on those constraints is a more trustworthy measure of what to expect for real surveys. Taken at face value, our results in Fig. 5 suggest that priors on the b φ (b 1 ) and b φδ (b 1 ) relations centered around the truth with widths of order 0.5 − 1 and 2.5 − 5, respectively (i.e., in between the results shown in orange and magenta in Fig. 5), may be necessary to guarantee unbiased constraints on f nl . Within the IllustrisTNG model, Fig. 3 suggests that this target might be achievable for the case of M * -selected galaxies with more simulations to beat the statistical errors. However, the figure shows that the b φ (b 1 ) and b φδ (b 1 ) relations for other selection criteria differ by more than these target widths, and it is still currently unclear how sensitive the bias relations are to changes in the star formation and baryonic feedback models in galaxy formation simulations. Further, recall also the lesson from Fig. 6 that considerations about the width of the priors on the bias relations cannot be disentangled from discussions about the assumed center values. This all strongly motivates more work to study the bias relations in galaxy formation models beyond IllustrisTNG, as well as for more refined criteria to select the objects, in particular, criteria that resemble as closely as possible those adopted for real galaxies. We highlight that progress along these lines may prove critical to the success of future galaxy surveys to improve over the CMB constraints on local PNG.

Results from parametrization 2: fit for products
Our current limited understanding of the b φ (b 1 ) and b φδ (b 1 ) relations motivates the exploration of constraints on local PNG using parametrization 2, which constrains the products f nl b φ and f nl b φδ , and avoids the need to put priors on the local PNG bias parameters. For our forecast setup, the overall significance of the detection of local PNG is then determined by the combined significance of the detection of nonzero values for the parameters The results are displayed in Fig. 7 in red color. For f nl b A φ and f nl b B φ , the top left panel shows that our fictitious survey would comfortably detect these parameters at the 1σ level. However, as noted previously in Ref. [44], these constraints are dominated by the power spectrum part of the data vector, without much contribution from the bispectrum. This can be seen by comparing the solid red contour in the top left panel, which is for the combined power spectrum and bispectrum data, with the dashed black contour for the P gg -only analysis: the two are barely distinguishable, which shows that the bispectrum constrains the f nl b A φ and f nl b B φ parameters only weakly. In parametrization 2, the bispectrum thus contributes to the overall detection of local PNG only via a detection of the parameters f nl or f nl b φδ , but as the lower left panel of Fig. 7 shows, these are strongly degenerate and poorly constrained. Thus, while parametrization 2 can prove a powerful way to detect local PNG with power spectrum data without making assumptions on galaxy bias, the price to pay is that the bispectrum would not be able to further enhance the significance of the detection. This can be compared with the case of parametrization 1 in Fig. 5, in which the bispectrum helps visibly in improving the bounds on f nl : for the "perfect knowledge" case, the uncertainty on f nl shrinks by about 30% when the bispectrum information is considered (cf. the two error bars in black in the lower panel of Fig. 5).
As an attempt to "rescue" the constraining power of the bispectrum, we also show results where, in parametrization 2, we replace f nl b A φδ with b A φδ , for which we assume Gaussian priors again. The results are shown in orange, magenta and green in Fig. 7 for varying prior widths, as labeled. As φδ has no visible impact on the f nl b A φ , f nl b B φ constraints, as they are dominated by the power spectrum data where b φδ does not contribute. The main effect is thus the progressive tightening of the f nl bound for tighter b φδ priors, as shown on the right of Fig. 7. Specifically, when we fit for f nl b φδ , we find σ fnl ≈ 100, but for our default Gaussian prior on b φδ , a prior that is 2× narrower, and another prior that is 4× narrower, the uncertainties become, respectively, σ fnl ≈ 30, σ fnl ≈ 15 and σ fnl ≈ 5. That is, it is only after the prior width on b φδ becomes of order unity (green; 5/4 = 1.25 ∼ 1) that the constraints return σ fnl values comparable to our fiducial choice of f nl = 5, or in other words, that the bispectrum begins to add to the overall significance of the detection of local PNG.
Before concluding, we clarify the origin of the peculiar shape of the orange contour in the upper right panel of Fig. 7, namely its sudden widening in the f nl direction. It is possible to show that the ∝ 2b 3 1 f nl and ∝ b 2 1 f nl b φδ contributions to the galaxy bispectrum (cf. Eq. (B.6)) have a similar scale-dependence on the range of scales that have the most constraining power (this is especially the small-scale squeezed bispectrum); this is in fact the reason behind the strong degeneracy between f nl and f nl b φδ discussed above. If the degeneracy was perfect and if b φδ = −2b 1 , then the bispectrum would become independent of f nl in parametrization 2, and unable to constrain it. In reality, the degeneracy is not perfect, but as the chains approach b A φδ = −2b A 1 ≈ −3.16 (which is where the widening occurs), the sensitivity to f nl still gets significantly reduced, and the contours get wider.

Summary and conclusions
Observational constraints on f nl using galaxy power spectrum and bispectrum data are determined largely from contributions ∝ b φ f nl and ∝ b φδ f nl , and as a result, competitive constraints on f nl using these data require priors on the bias parameters b φ and b φδ to be assumed. The most popular approach encountered in the literature (both in real-data constraints and forecasts) involves using the so-called universality relations of Eqs. (1.3) and (1.5) (or certain variants thereof) to relate these two parameters to the parameter b 1 that can be constrained using the parts of the data where f nl contributes only weakly. The problem is that these relations are derived for dark matter halos assuming that their mass function is universal, and thus, there is no reason to expect them to hold for real-life tracers of the large-scale structure like galaxies.
In fact, prior to this work, it was already known that the universality relation does not describe perfectly the b φ (b 1 ) relation of halos in gravity-only simulations, nor the relation for stellar-mass selected galaxies in hydrodynamical simulations. In this paper, we tested, for the first time with dedicated measurements from separate universe N -body simulations, the validity of the universality relation for the b φδ parameter as well, which enters at leading order in the galaxy bispectrum. We carried out our analysis for both dark matter halos in gravity-only simulations, as well as simulated galaxies in hydrodynamical simulations with the IllustrisTNG galaxy formation model. We studied in particular the sensitivity of the b φ (b 1 ) and b φδ (b 1 ) relations to different ways of selecting the galaxy samples, including in terms of total mass, stellar mass, black hole mass, black hole mass accretion rate and color. Using an idealized forecast setup with galaxy power spectrum (multitracer) and bispectrum data for a fictitious survey with V = 100Gpc 3 /h 3 at z = 1, and a fiducial value of f nl = 5, we then explored ways for how to incorporate uncertainties on the b φ (b 1 ) and b φδ (b 1 ) relations in f nl constraint analyses.
Our main results can be summarized as follows: • The b φδ (b 1 ) relation of dark matter halos is approximately redshift-independent, but it is not adequately described by the universality relation. The most notable difference is an overprediction by the latter in the range 1 b 1 3, by up to ∆b φδ ≈ 3 (cf. Fig. 2).
• For the simulated galaxies, the b φ (b 1 ) and b φδ (b 1 ) relations are both very sensitive to the galaxy selection criteria, are poorly described by the corresponding universality relations, and do not generically admit a regular, redshift-independent function of b 1 (cf. Fig. 3).
• Priors on b φ (b 1 ) and b φδ (b 1 ) centered close to the fiducial with widths of order 1 and 5, respectively, proved sufficient to yield relatively unbiased f nl constraints (cf. Figs. 5 and 6). This sets a rough target for how precisely these relations may need to be determined from simulations.
• Fitting for products of f nl b φ and f nl b φδ still allows to rule out f nl = 0 without assumptions on galaxy bias, but renders the bispectrum less useful in the constraints. For our fictitious survey analysis, the bispectrum only begins improving the overall detection of local PNG if priors on b φδ (b 1 ) of order 1 are assumed (cf. Fig. 7).
Overall, our current poor understanding of the b φ and b φδ parameters strongly motivates more works like this one to pin down the expected values of these parameters for real-life galaxy samples. For example, it is important to determine the extent to which the b φ (b 1 ) and b φδ (b 1 ) relations depend on the assumed galaxy formation physics. This can be done with simulations for variants of the IllustrisTNG model parameters, or for different self-consistent models of galaxy formation altogether. Given the sensitivity to the galaxy selection criteria that we encountered in this paper, it is also important that future works with galaxy formation simulations attempt also to mimic as closely as possible the selection strategy applied to the real-life galaxy samples considered. Further, beyond galaxies as tracers, it would be interesting to study also the b φ and b φδ parameters of the gas distribution [82], in particular, the distribution of neutral hydrogen mapped with 21cm line intensity mapping observations. It is important to acknowledge, however, that the field of cosmological hydrodynamical simulations of galaxy formation has still many associated uncertainties, and that it can be nontrivial to work out the selection functions of real-life galaxies and apply them on simulated ones. This suggests that obtaining very tight theoretical priors on the b φ (b 1 ) and b φδ (b 1 ) relations can prove challenging, but this is work that needs to be carried out anyway to help shape our approach to f nl constraints using large-scale structure data. For example, even if at the end of the day the conclusion is that our priors on these bias relations are not satisfactory, this will still be informative in that it will suggest abandoning approaches that focus on the numerical value of f nl (like in parametrization 1), and limiting ourselves to analyses based on parametrization 2 more focused on simply detecting f nl = 0.
We finish by highlighting also the importance for future studies on f nl to begin taking galaxy bias uncertainties into account. The strategies that we described in Sec. 4 are straightforward to implement in real-data constraint analyses like those of Refs. [24,25] using eBOSS survey data, as well as in the many forecast codes that exist in the literature for surveys like Euclid [93], SphereX [27] or SKA [94]. Doing so is important not only to guarantee robust and trustworthy bounds on f nl , but also to determine more precisely the theoretical precision requirements on the b φ and b φδ parameters for future galaxy surveys to reach their f nl targets.
A Validation of the method to estimate b φδ using the b 2 parameter In this appendix we estimate the bias parameter b 2 using a similar method to that used to estimate the b φδ parameter. This allows us to assess the performance of our method by comparing against known results for the b 2 (b 1 ) relation of dark matter halos in the literature. Up to a subtle caveat that we comment on below, the method described here is in all similar to that used also by Ref. [95].
Starting again from the galaxy-matter cross-power spectrum described by P gm = b 1 P mm , and defining its linear response to mass perturbations δ m as R 1,gm = ∂lnP gm /∂δ m , we have where R 1,mm = ∂lnP mm /∂δ m is the linear response of the matter power spectrum, and in the second equality we have used that b 2 = (∂ 2n g /∂δ 2 m )/n g . Thus, given measurements of the power spectrum responses from separate universe simulations and estimates of the values of b 1 , we can use the above equation to fit for b 2 . Alternatively, and as we did in our main results, we can estimate b 2 with higher signal-to-noise by computing the response of b 1 measured in the separate universe simulations, i.e., using directly where ∂lnb 1 ∂δ m = 1 2 simulations. This is the same as Fig. 2 in the main text, but for b2, instead of b φδ . In both panels, the dashed line shows the fit obtained by Ref. [51] for dark matter halos using separate universe simulations (cf. Eq. (3.1)).
In these equations, δ Highδm L (z) and δ Lowδm L (z) are the amplitudes of the matter perturbation in the separate universe cosmologies, and the superscripts indicate in which cosmology b 1 is evaluated using Eq. (2.2). We evaluate the above equations using the separate universe simulations presented in Refs. [42,77], which correspond to the same Fiducial cosmology used in the main body of the paper, and two separate universe cosmologies, Highδ m and Lowδ m , characterized by δ Highδm L (z = 0) = +0.05 and δ Lowδm L (z = 0) = −0.05, respectively. The caveat that is important to mention here concerns the choice of the box size in these separate universe simulations. Contrary to the changes in A s that characterize the HighA s and LowA s cosmologies in the main body of the paper, the values of δ L (z) modify the time-evolution of the scale-factor a(t), and consequently, the relationship between comoving and physical distances. The separate universe simulations in Refs. [42,77] are for fixed comoving size of the boxes at all epochs in units of Mpc. This is important for how to interpret the corresponding response measurements since the finite-differences evaluated with these simulations are then being taken at fixed comoving volume, but strictly, the full responses in the equations above include the effects from changing the volume. For the case of the power spectrum responses, this means that the measured response corresponds actually to the so-called growth-only piece. Explicitly, the power spectrum response functions can be decomposed as [85,96]  where G 1,mm and G 1,gm are the growth-only responses that one measures with fixed-comoving-volume separate universe simulations, and the remaining terms account for so-called reference density and dilation effects, that describe respectively, the fact that in the separate universe cosmology the power spectrum is measured w.r.t. a modified mean density, and that the physical scales get modified (or diluted) by the modified scale factor; the missing +1 term in R 1,gm has to do with the fact that there is one less instance of δ m in P gm and the galaxy overdensity is typically measured with respect to the observed local, mean galaxy density (and not the global one). Plugging these expressions in Eq. (A. 1) gives (noting that on large scales the two logarithmic derivative terms are the same) where the second equality defines the fixed-volume derivative of b 1 that we actually evaluate with our separate universe simulations using Eqs. We thus estimate b 2 using b 2 = ∂lnb 1 ∂δ m fixed vol.
The result is in Fig. 8, which shows the same as Fig. 2, but for b 2 instead of b φδ . Our estimated b 2 (b 1 ) relation recovers well the expected result shown by the dashed line, obtained by Ref. [51] for dark matter halos using separate universe simulations (cf. Eq. (3.1)). Likewise for b φδ , the results have higher signal-to-noise in the L box ≈ 800 Mpc, compared to the TNG300-2 box. Overall, this successful recovery of the expected b 2 (b 1 ) relation validates our estimates of b φδ (b 1 ) in the main body of the paper, which were obtained effectively in the same manner.
In our analysis, we consider k max = 0.2h/Mpc, on which corrections to our tree-level theory model do become important, but since the data vector is obtained from the theory model this does not impact our main conclusions on the impact of galaxy bias uncertainties. For simplicity, we also skip modeling the effects of redshift space distortions [98], so-called projection/relativistic effects [28,30,35,[99][100][101][102][103], and observational systematics [104]. Any reduction in constraining power from taking these complications into account reduces the overall importance of the uncertainties on galaxy bias, but only in the sense that there are more sources of uncertainty worsening the constraints on f nl . Our results in the main body of the paper thus assess how well we need to understand the b φ and b φδ parameters, assuming negligible sources of error from other aspects of observational constraints on local PNG.
The covariance matrix of the data vector can be written as where Cov P P , Cov BP and Cov BB indicate the covariance of the power spectrum part of the data vector, the cross-covariance of the bispectrum and power spectrum, and the covariance of the bispectrum part, respectively. For Cov P P we consider only the Gaussian+shot noise contribution, for Cov BB we consider in addition to the Gaussian contribution also contributions from leading-order non-Gaussian terms, and we set Cov BP = 0 for simplicity. We do not repeat our covariance expressions here, but the interested reader can find them in Ref. [44]; see also App. B there for a look into the impact of different covariance compositions (including the cross-covariance term Cov BP = 0) on the final f nl bounds, which is another aspect of these analyses that is often overlooked too.