Light Fields during Inflation
from BOSS and Future Galaxy Surveys

Daniel Green,★ Yi Guo,★ Jiashu Han ★ and Benjamin Wallisch ♠,★,◆,♣

★ Department of Physics, University of California San Diego, La Jolla, CA
92093, USA

♠ Oskar Klein Centre, Department of Physics, Stockholm University, 10691
Stockholm, Sweden

◆ School of Natural Sciences, Institute for Advanced Study, Princeton, NJ
08540, USA

♣ Center for Cosmology and Astroparticle Physics, Weinberg Institute for
Theoretical Physics,
Department of Physics, The University of Texas at Austin, Austin, TX 78751,
USA

Abstract
Primordial non-Gaussianity generated by additional fields present during
inflation offers a compelling observational target for galaxy surveys. These
fields are of significant theoretical interest since they offer a window into
particle physics in the inflaton sector. They also violate the single-field
consistency conditions and induce a scale-dependent bias in the galaxy power
spectrum. In this paper, we explore this particular signal for light scalar
fields and study the prospects for measuring it with galaxy surveys. We find
that the sensitivities of current and future surveys are remarkably stable for
different configurations, including between spectroscopic and photometric
redshift measurements. This is even the case at non-zero masses where the
signal is not obviously localized on large scales. For realistic galaxy number
densities, we demonstrate that the redshift range and galaxy bias of the
sample have the largest impact on the sensitivity in the power spectrum. These
results additionally motivated us to explore the potentially enhanced
sensitivity of Vera Rubin Observatory’s LSST through multi-tracer analyses.
Finally, we apply this understanding to current data from the last data
release of the Baryon Oscillation Spectroscopic Survey (BOSS DR12) and place
new constraints on light fields coupled to the inflaton.

###### Contents

  1. 1 Introduction
  2. 2 Light Fields and Galaxies
    1. 2.1 Galaxy Power Spectrum
    2. 2.2 Scale-Dependent Bias
  3. 3 Information Content of the Galaxy Power Spectrum
    1. 3.1 Fisher Information Matrix
    2. 3.2 Information and Survey Design
  4. 4 Constraints and Data Analysis for General Scaling Exponents
    1. 4.1 Search Strategy
    2. 4.2 BOSS DR12 Analysis
  5. 5 Multiple Tracers and the Dependence on Astrophysics
    1. 5.1 Optimizing Spectroscopic Multi-Tracer Analyses
    2. 5.2 Multi-Tracer LSST and Comparison to SPHEREx
  6. 6 Conclusions and Outlook
  7. A Forecasting Details
    1. A.1 Galaxy Power Spectrum Model
    2. A.2 Survey Specifications
  8. B Measuring the Scaling Behavior

##  1 Introduction

The statistics of the primordial density perturbations offer a window into the
dynamics of the very early universe [1, 2, 3], in particular the inflationary
epoch [4, 5]. While current maps of the universe are consistent with purely
Gaussian fluctuations [6], non-Gaussian correlation functions encode the
particles and interactions relevant to the origin of structure [7, 8, 5].
Inflation predicts a lower bound on the amount of primordial non-Gaussianity
(PNG) due to gravity alone that is approximately two orders of magnitude below
current constraints [9]. On the other hand, signals that are large enough to
be detected in the next generation of surveys arise in many models, including
examples with only Planck-suppressed interactions [10].

Due to the limited number of modes remaining to be measured to the cosmic
variance limit in the cosmic microwave background (CMB), surveys of the large-
scale structure (LSS) of the universe present the best opportunity to improve
our understanding of the primordial statistics [5, 11, 12, 13, 2]. With the
benefit of three-dimensional information, even current surveys have the raw
statistical power to compete with the CMB [14]. Unfortunately, for generic
non-Gaussian correlation functions, late-time nonlinearities make many modes
inaccessible and present a serious obstacle to our progress [15, 16, 17, 18].
In some circumstances, however, the manifestation of the primordial signal is
robust to nonlinear physics and can be measured reliably [1].

The fluctuations of extra fields beyond the inflaton present such an
opportunity. In the absence of these fields, correlation functions of matter
and galaxies are subject to the constraints of the single-field consistency
conditions [19, 20, 21]. In practice, these conditions are similar to the
equivalence principle [22] and ensure that correlations are determined by the
curvature of spacetime (i.e. derivatives of the metric fluctuations). Light
fields break these conditions and can lead to enhanced long-distance
correlations [23, 24, 25, 26, 27, 28]. These apparent violations of the
equivalence principle cannot arise from nonlinear gravity, and can therefore
be distinguished from astrophysical and gravitational sources of non-
Gaussianity.

From the perspective of both data analysis and survey design, the unique
advantage of these non-Gaussian signals is that they can be observed in the
two-point statistics of galaxies. Violations of the single-field consistency
conditions introduce couplings between long- and short-wavelength modes. The
short-wavelength modes control the number density of galaxies, which are then
non-trivially correlated on large scales in the presence of this type of non-
Gaussianity. The resulting (non-local) change to the shape of the galaxy power
spectrum, shown in Figure 1,

Figure 1: Illustration of the effect of the scale-dependent bias
$b_{\mathrm{NG}}^{\Delta}$ from a non-zero non-Gaussian amplitude
$f_{\mathrm{NL}}^{\Delta}$ on the linear galaxy power spectrum. We see that
the enhancement of power shifts from large to small scales as we increase the
exponent $\Delta$ (see §2.2 for its definition and details). The horizontal
dashed lines indicate the effective shot noise level for BOSS, DESI and
SPHEREx after scaling them to the displayed redshift $z=0$ and linear bias
$b_{1}=1.6$. The gray-shaded regions on large scales indicate the wavenumbers
below the minimum wavenumber of these surveys, $k<k_{\mathrm{min}}$, computed
based on their entire spherical volume. The gray-shaded regions for large
wavenumbers indicate the regimes where the scales exceed the nonlinear scale,
$k>k_{\mathrm{NL}}$, for the maximum redshift of BOSS and DESI, respectively.

is known as scale-dependent bias [29] and has been searched for in many
existing dataset (see e.g. [30, 31, 32, 33, 34, 35, 36, 37, 38, 39]).

The case of local primordial non-Gaussianity [40], for which the Gaussian and
non-Gaussian Newtonian potentials $\varphi$ and $\Phi$ are related by

$\Phi(\vec{x})=\varphi(\vec{x})-f_{\mathrm{NL}}^{\mathrm{loc}}\varphi^{2}(\vec{x})+\ldots\,,$
(1.1)

has been very well studied in the literature, both theoretically and
observationally. Local PNG arises in the presence of multiple massless fields,
such that $\varphi$ is an isocurvature mode during inflation, but is converted
to an adiabatic mode at later times [23, 24, 25]. Future measurements are
expected to reach particularly interesting thresholds for the physics of
multi-field inflation [15, 41, 42, 5].

The contribution to primordial non-Gaussianity from massive particles is much
less studied observationally. It is however a classic signature of quasi-
single-field inflation [26, 27, 28], which is also known as cosmological
collider physics [43]. Importantly, these models are compelling targets for
current and future surveys. In spite of this, little has been known about the
optimal survey strategy to search for these relics from inflation up to now.

In this paper, we will examine how future surveys can best constrain these
non-Gaussian signals using galaxy power spectra. For local PNG (massless
fields), the signal is dominated by the largest scales and most observational
strategies are designed accordingly [44, 45]. On the other hand, the signals
of massive fields may arise at large or small scales (see Fig. 1), which means
that the characteristics which best constrain these models is less clear a
priori [46, 47, 48]. We therefore study this signal in detail, and investigate
how redshift coverage, biasing, multi-tracer techniques and target selection
impact forecasts over the full range of scaling laws induced by these
particles.

We find that target selection, specifically finding highly biased objects, is
the largest factor in driving current and future measurements of this effect
when assuming a fixed universality relation. While large volumes help to
increase the number of modes, especially in spectroscopic surveys, we observe
that the benefits of going to higher redshifts are more driven by the large
biases of high-redshift objects and only secondarily by the increased volume.
In photometric surveys, the larger number density of galaxies also enables a
significant increase in sensitivity through sample-variance cancellation. This
requires splitting the sample according to the bias, which we explore in the
context of Vera Rubin Observatory’s Legacy Survey of Space and Time (LSST).

Using these insights, we explore current constraints on the scale-dependent
bias of galaxies using the BOSS DR12 dataset [49] for the full range of
scaling behaviors of light inflationary fields. We measure the contributions
to the power spectrum for one scaling exponent (PNG shape) at a time and
compute a correlation matrix to extrapolate between the discrete choices of
scaling dimensions (masses). While our inferred constraints are less sensitive
to the PNG amplitude than those from Planck, the CMB constraints for non-zero
masses are driven by bispectra in equilateral configurations rather than the
scaling behavior in the squeezed limit probed by the galaxy power spectrum.
The only previous constraint from LSS data in this regime was made in [50]. At
the same time, we forecast that galaxy power spectrum measurements in future
surveys will exceed the sensitivity of Planck to these light fields for a
sizable range of their masses.

This paper is organized as follows: In Section 2, we review the relevant
theoretical background for our non-Gaussian signal. In Section 3, we present
forecasts for a wide range of model parameters and experimental
configurations. Our goal is to identify the choices that most directly impact
the sensitivity to PNG beyond the local type. In Section 4, we apply our
understanding of the signal to BOSS DR12 galaxy clustering data and present
new constraints on non-Gaussianity from light fields. In Section 5, we discuss
the role of multi-tracer analyses and astrophysical effects on the forecasts.
We conclude in Section 6. A set of appendices contains technical details on
the modeling of the galaxy power spectrum and the survey specifications
(Appendix A), and a discussion on our ability to measure the scaling behavior
of the non-Gaussian signature from the galaxy power spectrum (Appendix B).

##  2 Light Fields and Galaxies

Our goal in this paper is to explore the signal of additional fields that
manifest themselves in the galaxy power spectrum via the scale-dependent bias.
This section reviews the necessary background to understand the signal.
Readers who are familiar with the scale-dependent bias due to general forms of
primordial non-Gaussianity (i.e. beyond the local type) may proceed to Section
3.

###  2.1 Galaxy Power Spectrum

The formation of structure in the universe is driven by the growth of density
fluctuations in the dark matter [51]. The evolution of the density contrast
Fourier mode in redshift space,
$\delta_{m}(\vec{k},z)\equiv\delta\rho_{m}(\vec{k},z)/\bar{\rho}_{m}(z)$, can
be solved at linear order to give

$\delta_{m}(\vec{k},z)=\frac{2k^{2}T(k)D(z)}{3\Omega_{m}H_{0}^{2}}\Phi(\vec{k})\equiv
k^{2}\mathcal{T}(k,z)\Phi(\vec{k})\,,$ (2.1)

where $k=|\vec{k}|$, $T(k)$ is the transfer function defined such that $T(k\to
0)\to 1$, and $D(z)$ is the linear growth factor normalized as $D(z)=1/(1+z)$
during matter domination. The primordial Newtonian potential $\Phi(\vec{k})$
encodes the primordial density fluctuations generated during inflation which
are concretely of the form

$\langle\Phi(\vec{k})\Phi(\vec{k}^{\prime})\rangle=\frac{9}{25}\frac{A_{\mathrm{s}}}{k^{3}}k^{n_{\mathrm{s}}-1}(2\pi)^{3}\delta_{D}(\vec{k}+\vec{k}^{\prime})\,,$
(2.2)

where $\delta_{D}$ is the Dirac delta function and we used
$\Phi=-\frac{3}{5}\zeta$, with curvature fluctuation $\zeta$, for modes that
re-entered the horizon during matter domination. The linear matter power
spectrum is therefore given by

$P_{\mathrm{lin}}(k)=k^{4}\mathcal{T}(k,z)^{2}\frac{A_{\mathrm{s}}}{k^{3-(n_{\mathrm{s}}-1)}}=\frac{4}{25}\frac{A_{\mathrm{s}}D(z)^{2}}{\Omega_{m}^{2}H_{0}^{4}}T(k)^{2}k^{n_{\mathrm{s}}}\,.$
(2.3)

The nonlinear matter power spectrum $P_{m}(k)$ follows the linear power
spectrum on large scales, $P_{m}(k)\approx P_{\mathrm{lin}}(k)$ for $k\ll
k_{\mathrm{NL}}$, for some scale
$k_{\mathrm{NL}}\approx$0.1\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
at $z=0$. Gravitational evolution corrects this behavior on smaller scales,
$k\gtrsim k_{\mathrm{NL}}$.

We will generally work with galaxies in the regime $k<k_{\mathrm{NL}}$. The
nonlinear effects that are relevant in this regime can be described solely in
terms of the bias expansion [52, 53]. In its most general form, the bias
expansion is simply the assumption that galaxy formation is a local
process:111Technically, biasing is best understood as a local process in
Lagrangian space. The non-locality in Eulerian space is captured by the
Zel’dovich approximation [53].

$\delta_{g}(\vec{x})=\sum_{i}b_{\mathcal{O}_{i}}\mathcal{O}_{i}(\vec{x})\,,$
(2.4)

where the operators $\mathcal{O}_{i}(\vec{x})$ are any locally measurable
quantities. Since we typically assume that the matter controls the formation
of galaxies, these are usually local products of $\delta_{m}(\vec{x})$ or the
tidal tensor [54],
$s_{ij}=\left(\frac{\nabla_{i}\nabla_{j}}{\nabla^{2}}-\frac{1}{3}\delta_{ij}\right)\delta_{m}(\vec{x})$,

$\mathcal{O}\in\\{\delta_{m},\nabla^{2}\delta_{m},\delta_{m}^{2},s_{ij}^{2},\delta_{m}^{3},\delta
s_{ij}^{2},\mathrm{Tr}[(\Pi^{[1]})^{3}]|^{(3)},\mathrm{Tr}[\Pi^{[1]}\Pi^{[2]}]|^{(3)},\cdots\\}\,,$
(2.5)

where the tensors $\Pi^{[n]}_{ij}$ are defined in [52].222The superscript
$(3)$ indicates that the operator includes only terms up to third order in
perturbation theory. In principle, there are an infinite number of operators
to consider. At any fixed accuracy, we however understand this as an expansion
in powers of the small density contrast, $\delta_{m}^{n}\ll 1$, and gradients
$R^{2}\nabla^{2}\ll 1$, for some fixed scale $R$. To model the galaxy power
spectrum at one-loop order, we consider the following set of operators (up to
third order in $\delta_{m}$) in the bias expansion:

$\displaystyle\delta_{g}(\vec{x})=$ $\displaystyle\
b_{1}\delta_{m}(\vec{x})+b_{\nabla^{2}}R^{2}\nabla^{2}\delta_{m}(\vec{x})+b_{\delta^{2}}(\delta_{m}^{2}(\vec{x})-\sigma^{2})+b_{s^{2}}(s^{2}(\vec{x})-\langle
s^{2}\rangle)$
$\displaystyle+b_{\Pi\Pi^{[2]}}\mathrm{Tr}[\Pi^{[1]}\Pi^{[2]}]|^{(3)}(\vec{x})+\cdots\,,$
(2.6)

with $\sigma^{2}=\langle\delta_{m}^{2}\rangle$. The coefficients
$b_{\mathcal{O}_{i}}=b_{1},b_{\nabla^{2}},b_{\delta^{2}},b_{s^{2}},b_{\Pi\Pi^{[2]}},\cdots$,
which we refer to as bias parameters, are constants that are determined by the
details of galaxy formation and evolution. Put differently, we parameterize
the complex physics of galaxies by these coefficients.

On large scales and for initial conditions set by single-field inflation, this
is the complete list of operators for practical purposes on large scales.
Because the fluctuations are necessarily adiabatic throughout cosmic history,
the evolution of the universe is controlled by a single statistical quantity,
which we have chosen to represent in terms of $\Phi$ and $\delta_{m}$. On
smaller scales, the evolution of different components of the universe can
however lead to some new types of biasing terms.

The redshift evolution of the bias parameters, particularly $b_{1}(z)$, are
important for projecting the sensitivity of future surveys. If
(proto-)galaxies form at a high redshift and primarily evolve with the
expansion of the universe, then we expect that the large-scale comoving galaxy
power spectrum remains constant, $P_{gg}(k,z)\approx P_{gg}(k,z^{\prime})$. In
other words, the power spectrum is due to the inhomogeneities in the density
field at the time they are formed, as one would find in Lagrangian biasing.
This holds if

$b_{1}(z)=D(z)^{-1}\,b_{1}(z=0)\,,$ (2.7)

which is a common simplifying assumption in many forecasts. Changes to the
galaxy sample through mergers, for example, can however lead to behavior that
differs significantly from Lagrangian biasing. When possible, we will
therefore use the biases for surveys that were defined by target selection,
but the assumption of this type of time evolution is also important in the
context of multi-tracer analyses discussed in Section 5.

Given the bias expansion for the galaxy overdensity, we can write down our
model for the galaxy power spectrum after including a few additional effects.
First, we work in Fourier space by transforming from real-space coordinates
$\vec{x}$ to $(k,\mu)$, where $\mu$ is the cosine between the wavevector
$\vec{k}$ and the line-of-sight direction. (Note that the composite operators
$\mathcal{O}_{i}(\vec{x})$ are simply related to
$\tilde{\mathcal{O}}_{i}(k,\mu)$ by a single Fourier integral.) Second, there
are redshift-space distortions which arise from the peculiar velocities of
galaxies. We account for the Kaiser effect by adding $\Delta
b(\mu,z)=f(z)\mu^{2}$, with the linear growth rate
$f\equiv\mathrm{d}\hskip-1.0pt\log{D}/\mathrm{d}a$, to the linear bias $b_{1}$
[55]. Collecting all terms in the bias expansion at linear order in
$\delta_{m}(\vec{k})$ therefore results in

$b(k,\mu,z)=b_{1}(z)+f(z)\mu^{2}+\sum_{n>0}b_{k^{2n}}(kR_{*})^{2n}\,,$ (2.8)

where the last term are the gradient biases in Fourier space with the comoving
Lagrangian radius $R_{*}$ of the halos of interest. Third, we describe
gravitational nonlinearities in the very mildly nonlinear regime via the one-
loop terms $P_{22}$ and $P_{13}$ in standard Eulerian perturbation theory.
Finally, the leading-order stochastic contribution is given by
$P_{\epsilon_{0}}$. Putting it all together, we therefore model the
theoretical anisotropic power spectrum for two galaxy samples $A$ and $B$ as
follows:

$\displaystyle P^{\mathrm{th}}_{AB}(k,\mu,z)=$ $\displaystyle\
b^{A}(k)b^{B}(k)\\!\left[P_{\mathrm{lin}}(k)+P_{22}(k)+P_{13}(k)\right]+P_{\epsilon_{0}}\delta_{AB}$
$\displaystyle+\left[b^{A}(k)b^{B}_{\delta^{2}}+b^{A}_{\delta^{2}}b^{B}(k)\right]\\!P_{\delta^{2}}(k)+\left[b^{A}(k)b^{B}_{s^{2}}+b^{A}_{s^{2}}b^{B}(k)\right]\\!P_{s^{2}}(k)$
$\displaystyle+b^{A}_{\delta^{2}}b^{B}_{\delta^{2}}\\!\left(P_{\delta^{2}\delta^{2}}(k)-2\sigma^{4}\right)+\left(b^{A}_{\delta^{2}}b^{B}_{s^{2}}+b^{A}_{s^{2}}b^{B}_{\delta^{2}}\right)\\!\left(\\!P_{\delta^{2}s^{2}}(k)-\tfrac{4}{3}\sigma^{4}\hskip-0.5pt\right)$
(2.9)
$\displaystyle+b^{A}_{s^{2}}b^{B}_{s^{2}}\\!\left(\\!P_{s^{2}s^{2}}(k)-\tfrac{8}{9}\sigma^{4}\hskip-0.5pt\right)+\tfrac{3}{2}\\!\left[b^{A}(k)b^{B}_{\Pi\Pi^{[2]}}+b^{A}_{\Pi\Pi^{[2]}}b^{B}(k)\right]\\!P_{\Pi\Pi^{[2]}}(k)\,,$

where we have suppressed the dependence on $\mu$ and $z$ on the right-hand
side, explicitly define all loop and bias contributions in Appendix A, and
$\delta_{AB}$ is the Kronecker delta. To illustrate these terms, we display
the individual contributions in Fig. 2.

Figure 2: Contributions to the galaxy power spectrum $P_{gg}^{(i)}(k)$ at
$z=0$ from the nonlinear, bias and non-Gaussian terms. The linear bias is
chosen to be $b_{1}=1.6$ and the gradient biases $b_{k^{2n}}$ (with $n\leq 2$)
are taken to be unity. The bias parameters $b_{s^{2}}$ and $b_{\Pi\Pi^{[2]}}$
are calculated in the Lagrangian local-in-matter-density model [52], and
$b_{\delta^{2}}$ is calculated using the halo simulation fit of [56]. For
comparison, we already include the scale-dependent bias from general PNG,
which we introduce in §2.2. The non-Gaussian parameter
$f_{\mathrm{NL}}^{\Delta}$ is set to unity for both cases of $\Delta=0$ and
$\Delta=2$. The dashed lines indicate negative contributions.

In practice, we compute the linear matter power spectrum
$P_{\mathrm{lin}}(k)$, which underlies all contributions of the theoretical
galaxy power spectrum $P^{\mathrm{th}}(k)$, with the Boltzmann solver CAMB
[57] using a fiducial $\Lambda$CDM cosmological model based on the Planck 2018
TT, TE, EE + lowE + lensing + BAO best-fit cosmology [58]. In Table 1,

Parameter | Description | Fiducial value
---|---|---
$\omega_{b}$ | Physical density of baryons $\omega_{b}\equiv\Omega_{b}h^{2}$ | $0.02242$
$\omega_{c}$ | Physical density of cold dark matter $\omega_{c}\equiv\Omega_{c}h^{2}$ | $0.11933$
$100\theta_{s}$ | Angular size of the sound horizon at recombination | $1.04119$
$\tau$ | Optical depth due to reionization | $0.0561$
$\ln{(${10}^{10}$A_{\mathrm{s}})}$ | Logarithm of the primordial scalar amplitude | $3.047$
$n_{\mathrm{s}}$ | Scalar spectral index | $0.9665$
$f_{\mathrm{NL}}^{\Delta}$ | Non-Gaussian amplitude | $0$
$b_{1}$ | Linear bias | $b_{1}(z)$
$b_{k^{2n}}$ | Gradient biases ($n\leq 2$) | $0$
$b_{\delta^{2}}$ | Quadratic bias | $0$
$b_{s^{2}}$ | Tidal bias | $0$
$b_{\Pi\Pi^{[2]}}$ | Evolution bias | $0$
Table 1: Parameters of the fiducial $\Lambda$CDM model, based on the Planck
2018 best-fit cosmology [58], with the sum of neutrino masses $\sum
m_{\nu}=$0.06\text{\,}\mathrm{e}\mathrm{V}$$, and the biasing model employed
in our forecasts. Except for the optical depth, we vary these parameters in
our forecasts unless stated otherwise.

we list the fiducial values of the $\Lambda$CDM, non-Gaussian and bias
parameters. We set the fiducial values of the three loop biases
$b_{\delta^{2}}$, $b_{s^{2}}$ and $b_{\Pi\Pi^{[2]}}$ to $0$ instead of using
the Lagrangian local-in-matter-density model [52] and fit results from halo
simulations [56], as we did in Fig. 2, since we cover a larger redshift range
than the fitted bias relations were derived from. This choice should have
minimal impact on the results of this paper because the forecasted constraints
on the loop biases are large enough for both choices of fiducial values to be
consistent with each other.

Finally, to relate the theoretical galaxy power spectrum
$P^{\mathrm{th}}_{AB}(k,\mu,z)$ to the observed galaxy power spectrum
$P^{\mathrm{obs}}_{AB}(k,\mu,z)$, we account for two observational effects
that are present in any real survey: redshift errors and the Alcock-Paczynski
effect. Redshift errors are observational uncertainties in the measurement of
galaxy redshifts, $\sigma_{z0}(1+z)$, which we model as an exponential
suppression of the power spectrum along the line of sight,

$f^{A}_{\sigma_{z}}(k,\mu,z)=\exp\hskip-1.0pt\left\\{-k^{2}\mu^{2}[\sigma^{A}_{z0}(1+z)]^{2}/H^{2}(z)\right\\}\,,$
(2.10)

where $\sigma^{A}_{z0}$ is the root-mean-square redshift error of sample $A$
at $z=0$ [59, 60]. (We neglect this effect for spectroscopic surveys by
setting $\sigma^{A}_{z0}\equiv 0$, i.e. $f^{A}_{\sigma_{z}}=1$.) The Alcock-
Paczynski effect arises when the true cosmology differs from the fiducial
cosmology that is used to convert the measured angular positions and redshifts
of the galaxies to comoving wavenumbers in Fourier space. The effect of this
mapping on the density contrast is captured by [61],

$f_{AP}(z)=\frac{D^{\mathrm{fid}}_{A}(z)^{2}}{D_{A}(z)^{2}}\frac{H(z)}{H_{\mathrm{fid}}(z)}\,,$
(2.11)

where $D_{A}(z)$ is the angular diameter distance to redshift $z$. To
summarize, our model for the observed galaxy power spectrum in redshift space
therefore is

$P^{\mathrm{obs}}_{AB}(k,\mu,z)=f_{AP}(z)\sqrt{f^{A}_{\sigma_{z}}(k,\mu,z)f^{B}_{\sigma_{z}}(k,\mu,z)}\,P^{\mathrm{th}}_{AB}(k,\mu,z)\,,$
(2.12)

where $P^{\mathrm{th}}_{AB}(k,\mu,z)$ is given by (2.9).

###  2.2 Scale-Dependent Bias

The conventional description of biasing applies when the initial conditions
are Gaussian. Gaussian initial conditions do not correlate modes of different
scales, i.e. the collapse of small-scale density fluctuations to form halos in
different regions is determined by locally observable quantities, such as the
matter density, that vary over cosmological distances due to the long-
wavelength fluctuations. On the other hand, mode coupling in the initial
conditions due to primordial non-Gaussianity can introduce long-range
correlations that are purely related to the statistics of the small-scale
fluctuations themselves. When interpreted through the bias expansion, this
mode coupling may appear to be non-local (scale-dependent) bias [29, 62, 63,
64] and/or long-range stochastic bias [65]. We will focus on the former since
the latter is typically a subdominant contribution to the signal.

The coupling of long and short modes is described by the squeezed or collapsed
limits of a non-Gaussian correlation function. In single-field inflation, the
bispectrum is constrained by the single-field consistency conditions to take
the following form in the squeezed limit [66, 19]:

$\lim_{\vec{k}_{1}\to
0}\big{\langle}\zeta(\vec{k}_{1})\zeta(\vec{k}_{2})\zeta(\vec{k}_{3})\big{\rangle}=-\left[(n_{\mathrm{s}}-1)+O(k_{1}^{2}/k_{3}^{2})\right]P_{\zeta}(k_{1})P_{\zeta}(k_{3})\,(2\pi)^{3}\delta_{D}\\!\left(\sum\nolimits_{i}\vec{k}_{i}\right).$
(2.13)

Furthermore, the leading term, which corresponds to
$f_{\mathrm{NL}}^{\mathrm{loc}}=-\frac{5}{12}(n_{\mathrm{s}}-1)$, is
unphysical [67, 68]. The leading physical term is suppressed by
$O(k_{1}^{2}/k_{3}^{2})$ and is due to the coupling of the short modes at
horizon crossing to the curvature of the universe due to long-wavelength modes
[69], which is typically proportional to the amplitude of equilateral non-
Gaussianity, $f_{\mathrm{NL}}^{\mathrm{eq}}$.

In the presence of additional fields, the short-scale power can depend on the
long-wavelength values of these fields and not just derivatives of the metric
fluctuations [70]. For an additional massless field, the late-time Newtonian
potential $\Phi$ may be nonlinearly related to a light (isocurvature) field
$\chi$ during inflation,
$\Phi(\vec{x})=\chi(\vec{x})-f_{\mathrm{NL}}^{\mathrm{loc}}\chi^{2}(\vec{x})$.
For a massive field, we have in the superhorizon limit

$\lim_{\vec{k}_{1}\to 0}\chi(\vec{k}_{1},t)\approx
c\frac{H}{k_{1}^{3/2}}\left(\frac{k_{1}}{aH}\right)^{\\!\Delta}\,,$ (2.14)

with a constant $c$. If $\chi$ alters the power spectrum of the short modes at
horizon crossing, $k_{2}=aH$, then we have a contribution to the bispectrum of
the form

$\lim_{\vec{k}_{1}\to
0}\big{\langle}\zeta(\vec{k}_{1})\zeta(\vec{k}_{2})\zeta(\vec{k}_{3})\big{\rangle}^{\prime}=f_{\mathrm{NL}}^{\Delta}\left(\frac{k_{1}}{k_{2}}\right)^{\\!\Delta}P(k_{1})P(k_{2})\,,$
(2.15)

where we defined $\langle\ldots\rangle=(2\pi)^{3}\hskip
0.5pt\delta_{D}(\sum_{i}\vec{k}_{i})\,\langle\ldots\rangle^{\prime}$. To have
such a contribution to the bispectrum, we also need $\chi$ to mix with
$\zeta$. Even without a mixing term, this same effect however appears in the
collapsed limit of the trispectrum. For light fields,
$\Delta=3/2-\sqrt{9/4-m^{2}/H^{2}}$ such that $m=0$ leads to local PNG
($\Delta=0$) as expected. For $m/H>3/2$, $\Delta=\frac{3}{2}\pm\mathrm{i}\nu$,
with $\nu=\sqrt{m^{2}/H^{2}-9/4}$, is complex in which case the complex
conjugate 333For heavy fields, there are two solutions of the long-wavelength
behavior which scale as $\Delta$ and its complex conjugate $\Delta^{\star}$.
We note that there is also a second solution $\Delta^{\prime}=(3-\Delta)$ for
light fields which is however suppressed relative to $\Delta$. also
contributes which results in the bispectrum being real valued, as required.
More generally, it is possible to construct models with a wider range of real
and complex values of $\Delta$ than those generated by a single massive field
[71, 72, 73, 74].

In single-field inflation, physical mode coupling must involve derivatives of
the metric fluctuations [67], whose leading behavior in the soft limit
corresponds to $\Delta=2$ [75]. One might naturally be surprised to find that
the $m^{2}\to\infty$ limit of $\Delta$ for a massive field does not yield this
value of single-field inflation. It was nicely explained in [43] that
$\Delta\approx 3/2+\mathrm{i}m/H$ is precisely what is expected from the
single-particle wavefunction of a massive particle in an expanding universe in
the high-mass limit. This contribution to the squeezed limit arises from the
physical production of these particles. This particle production is Boltzmann
suppressed in some simple models [76],
$f_{\mathrm{NL}}^{\Delta}\propto\mathrm{e}^{-\pi m/H}$, but can be enhanced
depending on the nature of the interaction [77]. In addition, very massive
fields will also produce a sub-leading contribution to the squeezed limit with
$\Delta=2$. This contribution arises from the virtual exchange of a massive
field and is equivalent to a purely local interaction (i.e. to integrating out
the massive field). It may be surprising that a local term will produce non-
local scale-dependent bias. When $T(k)\to 1$, $\Delta=2$ is indeed a local
term and it was long thought that $f_{\mathrm{NL}}^{\mathrm{eq}}$ could not be
measured via scale-dependent bias for this reason (see e.g. [64]). However,
because of the evolution of matter after horizon entry ($T(k)\neq 1$), local
interactions during inflation are distinguishable from local processes in
structure formation at late-times [48], i.e. equilateral PNG can be measured
in the power spectrum via scale-dependent bias.

In the conventional picture of galaxy biasing (see §2.1), the density contrast
of galaxies is determined by a long list of composite operators in the initial
density field. The origin of the unique signal of PNG in galaxies is that the
behavior of these composite operators becomes a proxy for the light fields
during inflation. For example, consider the impact of structure formation from
the local amplitudes of fluctuations on a scale $R$,

$\sigma_{R}(\vec{x})=\int\\!\mathrm{d}^{3}k\,\mathrm{e}^{-\mathrm{i}\vec{k}\cdot\vec{x}}\int\\!\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}\zeta(\vec{p})\zeta(\vec{k}-\vec{p})F(pR)\,,$
(2.16)

for some filter function $F(x)$. This is the local variance of the primordial
fluctuations which varies from place to place. Importantly, this variance is
correlated with the long-wavelength metric perturbations as follows:

$\langle\sigma_{R}(\vec{k})\zeta(\vec{k}^{\prime})\rangle^{\prime}\propto
f_{\mathrm{NL}}^{\Delta}(kR)^{\Delta}P_{\zeta}(k)\quad\Rightarrow\
\quad\langle\sigma(\vec{k})\delta_{m}(\vec{k}^{\prime})\rangle^{\prime}\propto
f_{\mathrm{NL}}^{\Delta}\frac{(kR)^{\Delta}}{k^{2}T(k)}P_{m}(k)\,.$ (2.17)

Since the number density of galaxies at a location $\vec{x}$ is related to
$\sigma_{R}(\vec{x})$, this introduces a long-distance correlation between
$\delta_{g}(\vec{x})$ and $\delta_{m}$ that is non-local. The resulting non-
local modification of the bias (2.8) is given by [65] 444It is important to
note that, in principle, the scale-dependent bias is distinguishable from the
expansion in $(kR_{\star})$ even for $\Delta=2n$ for any positive integer $n$.
At $k\ll k_{\mathrm{eq}}$, $\mathcal{T}(k)\to\mathrm{const}$ and the scale-
dependent bias scales as $k^{2n-2}$ which is degenerate with $b_{1}$ and
$b_{k^{2(n-1)}}$. However, since $k_{\mathrm{eq}}R_{\star}\ll 1$, the transfer
function introduces a non-trivial scale dependence in the regime
$k_{\mathrm{eq}}\lesssim k\ll R_{\star}^{-1}$ which is not captured in the
local bias expansion.

$b(k,\mu,z)=b_{1}(z)+f(z)\mu^{2}+\sum_{n>0}b_{k^{2n}}(kR_{*})^{2n}+Af_{\mathrm{NL}}\frac{b_{\phi}(z)}{k^{2}\mathcal{T}(k,z)}(kR_{*})^{\Delta}\,,$
(2.18)

where $b_{\phi}$ is the non-Gaussian bias parameter and we took $R=R_{*}$. The
constant $A$ normalizes the non-Gaussian parameter $f_{\mathrm{NL}}^{\Delta}$
in the conventional way for PNG in the squeezed limit of the
bispectrum,555Note that for the special bispectrum configurations, when
expressed in terms of the general exponent $\Delta$, the equivalent
$f_{\mathrm{NL}}^{\Delta}$ definitions are related by
$f_{\mathrm{NL}}^{\mathrm{loc}}=3f_{\mathrm{NL}}^{\Delta=0}$ and
$f_{\mathrm{NL}}^{\mathrm{eq}}=f_{\mathrm{NL}}^{\Delta=2}$. so that the last
term in (2.18) is [63, 78, 79, 80, 81]

$\displaystyle b_{\mathrm{NG}}^{\mathrm{loc}}(k,z)$
$\displaystyle=f_{\mathrm{NL}}^{\mathrm{loc}}\frac{b_{\phi}(z)}{k^{2}\mathcal{T}(k,z)}$
$\displaystyle\text{(local)},$ (2.19) $\displaystyle
b_{\mathrm{NG}}^{\mathrm{eq}}(k,z)$
$\displaystyle=3f_{\mathrm{NL}}^{\mathrm{eq}}\frac{b_{\phi}(z)}{k^{2}\mathcal{T}(k,z)}(kR_{*})^{2}$
$\displaystyle\text{(equilateral)},$ (2.20) $\displaystyle
b_{\mathrm{NG}}^{\Delta}(k,z)$
$\displaystyle=3f_{\mathrm{NL}}^{\Delta}\frac{b_{\phi}(z)}{k^{2}\mathcal{T}(k,z)}(kR_{*})^{\Delta}$
$\displaystyle\text{(general\leavevmode\nobreak\ exponent\leavevmode\nobreak\
$\Delta\in[0,2]$)}.$ (2.21)

Note that we do not include the physical orthogonal shape, as it is produced
in single-field inflation and, therefore, has the same squeezed limit as the
equilateral shape, i.e. corresponds to $\Delta=2$. Previous work on the scale-
dependent bias used an orthogonal template which has $\Delta=1$ (e.g. [47]),
which is however unphysical and has not been used in recent LSS analyses.

Without any additional theoretical input, $b_{\phi}(z)$ and
$f_{\mathrm{NL}}^{\Delta}$ are both unknown quantities and are completely
degenerate. If that was the case, we could only measure the combination
$f_{\mathrm{NL}}^{\Delta}b_{\phi}(z)$ in each redshift bin which would make it
challenging to infer limits on $f_{\mathrm{NL}}^{\Delta}$ alone. More
optimistically, we can fix $b_{\phi}(z)$ in terms of the linear bias
$b_{1}(z)$ using the universality relations found in phenomenological models
[52],

$b_{\phi}(z)=2\delta_{c}(b_{1}(z)-p)\,,$ (2.22)

where $\delta_{c}=1.686$ is the critical overdensity for spherical collapse
and $p=1$ if we assume a universal halo mass function. We will adopt this
relation in our forecasts, but refer to e.g. [82, 83, 84, 85] for a recent
discussion of the validity and risks of this choice.

To conclude, we note that we incorporate the general scale-dependent bias
(2.18) in the galaxy power spectrum as described in (2.9) and (2.12), i.e.
$b^{A}=b^{A}_{\Delta}(k,\mu,z)$. In the rest of the paper, we follow [48] in
assuming a minimum halo mass of $${10}^{13}$M_{\odot}$, which corresponds to
$R_{*}\approx$2.66\text{\,}{\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Mpc}$$,
and truncating the non-local gradient bias expansion at $n=2$. We illustrate
the effect of $b_{\mathrm{NG}}^{\Delta}$ on the linear galaxy power spectrum
in Fig. 1 and compare these non-Gaussian contributions to the nonlinear galaxy
power spectrum in Fig. 2 at redshift $z=0$ for linear bias $b_{1}=1.6$. We see
from the latter figure that the $\Delta=0$ contribution does not mimic the
behavior of any other bias or nonlinear terms at low $k$. On the other hand,
the $\Delta=2$ contribution is both much smaller at $k<k_{\mathrm{NL}}$ and
similar in behavior to those Gaussian contributions.

##  3 Information Content of the Galaxy Power Spectrum

Our goal in this section is to understand how to design a galaxy survey to
best measure $f_{\mathrm{NL}}^{\Delta}$ for $\Delta\in[0,2]$. To do so, we
will explore the Fisher forecasts for $f_{\mathrm{NL}}^{\Delta}$ as a function
of $\Delta$ for a variety of survey configurations and biasing models.666In
Appendix B, we also discuss our ability to measure the scaling behavior of the
non-Gaussian bias, where we also include $\Delta$ as an additional free
parameter in the Fisher matrix. We will compare these results with analytic
estimates to understand what aspects of the surveys drive major improvements
in sensitivity. While our forecasts will be rooted in specific survey
configurations, our goal is to isolate the features of these surveys that
impact the forecasted constraining power the most.

###  3.1 Fisher Information Matrix

In the following, we describe the well-known Fisher matrix forecasting
techniques for the galaxy power spectrum that we use. We review the basics
here and refer to Appendix A for the experimental specifications of the
surveys.

The Fisher matrix for a set of galaxy samples, labeled by $\\{A,B,\ldots\\}$,
is defined by

$F_{\alpha\beta}=\sum_{z_{i}}V(z_{i})\int_{-1}^{1}\\!\frac{\mathrm{d}\mu}{2}\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\frac{\mathrm{d}k\,k^{2}}{2\pi^{2}}\,\frac{1}{2}\mathrm{Tr}\left[\mathbf{C}_{,\alpha}\mathbf{C}^{-1}\mathbf{C}_{,\beta}\mathbf{C}^{-1}\right],$
(3.1)

where $\theta_{\alpha,\beta}$ are the model parameters, the trace is over the
samples and we defined
$\mathbf{C}_{,\alpha}\equiv\partial\mathbf{C}/\partial\theta_{\alpha}$. The
integrals are discretized, with the integration over wavenumbers $k$ being
limited by the minimum wavenumber $k_{\mathrm{min}}$ and the maximum
wavenumber $k_{\mathrm{max}}$, and bins of width $\Delta k=k_{\mathrm{min}}$.
We note that we vary all parameters listed in Table 1 throughout this paper,
unless specified otherwise. The elements of the matrix $\mathbf{C}$ are given
by

$\mathbf{C}_{AB}(k,\mu,z)\equiv
P^{\mathrm{obs}}_{AB}(k,\mu,z)+\delta_{AB}N_{A}(k,\mu,z)\,,$ (3.2)

with the shot noise term $N_{A}(k,\mu,z)=\bar{n}_{A}(z)^{-1}$, where
$\bar{n}_{A}(z)$ is the average comoving number density of tracers at redshift
$z$. Shot noise arises because we observe a finite number of discrete objects
in galaxy surveys. (Note that we absorbed the stochastic term $P_{\epsilon 0}$
of (2.9) into $N_{A}$.) We treat the redshift bins as spherical shells so that
the volume of the $i$-th redshift bin is

$V(z_{i})=\frac{4\pi}{3}{f_{\mathrm{sky}}}\left[d_{c}(z_{i}+\Delta
z_{i}/2)^{3}-d_{c}(z_{i}-\Delta z_{i}/2)^{3}\right],$ (3.3)

where ${f_{\mathrm{sky}}}$ is the sky fraction observed by the survey,
$d_{c}(z)$ is the comoving distance to redshift $z$ and $\Delta z_{i}$ is the
width of the $i^{\mathrm{th}}$ redshift bin.

To ensure that only modes under perturbative control in Eulerian perturbation
theory and the bias expansion enter our Fisher matrix calculation, and to
exclude long-wavelength modes inaccessible to a given survey due to its finite
volume, we restrict the integral in (3.1) to a maximum wavenumber
$k_{\mathrm{max}}$ and a minimum wavenumber $k_{\mathrm{min}}$, respectively.
We make a conservative choice for $k_{\mathrm{max}}$ by taking the smaller of
the two scales indicating the range of validity of perturbation theory and the
bias expansion. To be precise, we take it to be the minimum of the nonlinear
scale 777We use $k_{\mathrm{NL}}=\frac{\pi}{2R_{\mathrm{NL}}}$, with the
radius $R_{\mathrm{NL}}$ at which the variance of linear fluctuations is
$\sigma_{R}=1/2$ at a given $z$. $k_{\mathrm{NL}}$ and the wavenumber
$k_{\mathrm{halo}}\approx$0.19\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
associated with the Lagrangian scale corresponding to the minimum halo mass of
the surveyed population which we assume to be
$R_{*}\approx$2.66\text{\,}{\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Mpc}$$
(cf. §2.2), i.e.
$k_{\mathrm{max}}(z)=\min\\{k_{\mathrm{NL}}(z),k_{\mathrm{halo}}\\}$. For
$k_{\mathrm{min}}$, we assume a spherical survey geometry per redshift bin
which results in a conservative minimum wavenumber
$k_{\mathrm{min},i}=2\pi\hskip 1.0pt[3V(z_{i})/(4\pi)]^{1/3}$.

For most of the forecasts in Sections 3 and 5, we use a futuristic
spectroscopic survey with a billion objects divided into 10 redshift bins with
$z\leq 5$, which we will hereafter refer to as the billion-object survey. We
split the galaxy sample of this survey into two populations based on their
linear bias at $z=0$, chosen as $b_{1}^{(1)}(z=0)=2.0$ and
$b_{1}^{(2)}(z=0)=1.2$. We model their redshift evolution according to (2.7)
as previously discussed in §2.1. In addition, we assume a constant number
density across all redshift bins and both galaxy samples,
$\bar{n}^{(1)}_{g}(z)=\bar{n}^{(2)}_{g}(z)=\mathrm{const}$. The details of
this and other surveys used in our forecasts, i.e. the Baryon Oscillation
Spectroscopic Survey (BOSS) [49], the Dark Energy Spectroscopic Instrument
(DESI) [86], Euclid [87, 88], Vera Rubin Observatory’s LSST [89], the Spectro-
Photometer for the History of the Universe, Epoch of Reionization and Ices
Explorer (SPHEREx) [41], MegaMapper [90, 12] and the billion-object survey are
provided in Appendix A. Except for SPHEREx and the billion-object survey, we
treat all surveys as single-tracer surveys by combining different groups of
tracers (if available) into a single effective number density and bias. The
impact of this should not significantly impact our forecasts as we explain in
Appendix A.

Unless mentioned otherwise, we marginalize over all bias parameters and
include CMB information on the $\Lambda$CDM parameters (but not on PNG) from
Planck throughout the paper.888When marginalizing over the full biasing model
for the billion-object survey, for instance, the difference between using the
$\Lambda$CDM covariance from Planck or CMB-S4 [91, 92] for $\Delta\lesssim
1.5$ is less than two percent. In fact, this CMB information decreases
$\sigma(f_{\mathrm{NL}}^{\Delta})$ only by a few percent for $\Delta\lesssim
1.3$, and by a maximum of 20% (15%) with Planck (CMB-S4) at $\Delta=2$. We
achieve this by computing the $\Lambda$CDM Fisher matrix directly following
the Fisher methodology and the experimental specifications of [93] (see also
[94] for the latter), marginalized over the optical depth $\tau$. We then
combine this marginalized Planck Fisher matrix with the respective LSS Fisher
matrix for the five $\Lambda$CDM parameters, the non-Gaussian parameter
$f_{\mathrm{NL}}^{\Delta}$ and the galaxy bias parameters
$b_{\mathcal{O}_{i}}$. We treat each redshift bin with mean redshift
$\bar{z}_{i}$ as an independent survey and independently marginalize over the
biases within each redshift bin.

###  3.2 Information and Survey Design

Using Fisher forecasts, our aim is to understand the possible strategies for
optimizing the measurement of $f_{\mathrm{NL}}^{\Delta}$ for shapes including
but not limited to local PNG. The parameters that we most directly control in
the design of a survey are the sky fraction (${f_{\mathrm{sky}}}$), the
redshift range ($z\in[z_{\mathrm{min}},z_{\mathrm{max}}]$) and the number
density of tracers [$\bar{n}_{g}(z)$]. However, the actual statistical power
of the survey is controlled by the smallest ($k_{\mathrm{min}}$) and largest
($k_{\mathrm{max}}$) wavenumbers that can be reliably measured, the linear
bias of the sample [$b_{1}(z)$], etc. These factors are influenced by the
details of the survey, but can also be affected by systematics and
astrophysics.

We show the overall landscape of future power spectrum measurements of
$f_{\mathrm{NL}}^{\Delta}$ for $\Delta\in[0,2]$ in Fig. 3.

Figure 3: Forecasted constraints on $f_{\mathrm{NL}}^{\Delta}$ for general
scaling $\Delta\in[0,2]$ of the scale-dependent bias, for current and future
surveys, compared to the Planck bispectrum constraints for local and
equilateral non-Gaussianity. We also marginalize over the five $\Lambda$CDM
parameters, use Planck results as priors for the $\Lambda$CDM parameters, and
either fix the bias parameters (left) or marginalize over them (right). Note
that the SPHEREx forecast uses the five-tracer sample and the billion-object
survey forecast uses the double-tracer sample as provided in Appendix A. For
all the other surveys, we use their combined single-tracer sample.

Qualitatively, the forecasts show the overall behavior that we might expect:
(i) constraints at $\Delta=0$ are generally much stronger than at $\Delta=2$,
(ii) higher-order biasing affects larger $\Delta$ more than smaller $\Delta$,
(iii) larger surveys have more constraining power, and (iv) future LSS surveys
will be able to improve over Planck bispectrum constraints from the galaxy
power spectrum for $\Delta=0$, but not for $\Delta=2$ without additional
information on bias parameters. In detail, the forecasts however have elements
that are harder to understand without further investigation. First, the
forecasts show clear features at values of $\Delta\sim 1.4$ that depend on the
marginalization over biasing parameters (see also Fig. 6 below). This suggests
a qualitative change in where the constraining power is coming from as we vary
$\Delta$. Second, we see that SPHEREx and LSST produce very similar forecasts
despite having very different strengths and weaknesses. Furthermore, both
significantly exceed the constraining power of Euclid and DESI which have the
advantage of being three-dimensional (spectroscopic) surveys. In order to make
sense of these forecasts, we will break them down according to the number
density (noise), scales ($k_{\mathrm{min}}$/$k_{\mathrm{max}}$), the fiducial
bias [$b_{1}(z)$], the survey geometry (${f_{\mathrm{sky}}}$,
$z_{\mathrm{max}}$) and redshift errors ($\sigma_{z0}$) in the following.

#### Dependence on Number Density

The most basic parameter in any experiment is the signal to noise. In a galaxy
survey, shot noise is the dominant (irreducible) noise source. At fixed
volume, increasing the number of objects increases the range of wavenumbers
$k$ where $P_{gg}(k)$ is measured with at least signal to noise of order one,
$P_{gg}(k)/N(k)\gtrsim O(1)$.

For most cosmological parameters, sample variance presents an additional
irreducible source of noise. However, for scale-dependent bias, the parameter
$f_{\mathrm{NL}}^{\Delta}$ can, in principle, be measured without sample
variance when the signal to noise is large for multiple tracers [95].
Concretely, the relationship between $\delta_{g}(\vec{k})$ and
$\delta_{m}(\vec{k})$ in (2.6) is completely deterministic and is therefore
not limited by sample variance. However, this only works if we can measure the
same Fourier mode with high signal to noise for tracers with different values
of $b_{\phi}$. The noise introduced by sample variance is proportional to
$P(k)$, while shot noise is independent of $k$. As a result, sample-variance
cancellation dramatically alters how different scales contribute to the
measurement of $f_{\mathrm{NL}}^{\Delta}$ and, therefore, the qualitative
understanding of our forecasts. This is why it is important to distinguish
from the outset to what degree multi-tracer sample-variance cancellation is
relevant in our forecasts.

The impact of shot noise in both single- and multi-tracer scenarios is
illustrated in Fig. 4.

Figure 4: Dependence of $\sigma(f_{\mathrm{NL}}^{\Delta})$ on the shot noise,
i.e. the inverse survey number density $\bar{n}_{g}^{-1}$, for the billion-
object survey for four representative values of $\Delta$. The solid lines show
the single-tracer forecasts and the dashed lines display the double-tracer
results. The different colors indicate which bias parameters are marginalized
when calculating $\sigma(f_{\mathrm{NL}}^{\Delta})$. The vertical lines mark
the averaged shot noise of various surveys which are calculated by taking the
average of the number densities of the middle 50% of redshift bins of each
survey, i.e. they should be taken as an approximate illustration of the noise
level in these surveys.

Two points should be very clear from these curves: (1) sample-variance
cancellation has the potential to dramatically improve the sensitivity to
$f_{\mathrm{NL}}^{\Delta}$ for all values of $\Delta$, and (2) realistic
current and future surveys are very far from being in the regime where sample-
variance cancellation has a large impact. In fact, most near future surveys
are, at best, at the boundary between single and multiple tracers.999Note that
the curves in Fig. 4 are computed using the spectroscopic billion-object
survey (cf. Table 7). This means that the photometric LSST is not a direct
comparison and the line labeled as LSST only actually indicates a
spectroscopic follow-up of LSST. The potential sensitivity in a multi-tracer
analysis of LSST itself will be discussed in §5.2. As a result, sample-
variance cancellation is important for highly accurate forecasts, but the
qualitative behavior is consistent with single-tracer forecasts. We therefore
use a single-tracer emphasis for the rest of this section. We will revisit the
advantages offered by multi-tracer techniques and how astrophysics affects the
constraints on $f_{\mathrm{NL}}^{\Delta}$ in Section 5.

#### Dependence on Scales

The defining feature of scale-dependent bias is that it is enhanced at small
wavenumbers/large distances. For $\Delta<2$ and $f_{\mathrm{NL}}^{\Delta}\neq
0$, we have $P_{gg}(k)\gg b_{1}^{2}P_{m}(k)$ as $k\to 0$. Since the noise is
also a function of $k$, it is however not a given that the information resides
at small $k$ for all $\Delta<2$. We now investigate this analytically in the
Fisher matrix and numerically in our forecasts.

For a single-tracer analysis where we know the bias parameters exactly, i.e.
$b_{\mathcal{O}_{i}}$ are held fixed, the Fisher information for a fiducial
$f_{\mathrm{NL}}^{\Delta}=0$ is given by

$F_{f_{\mathrm{NL}}^{\Delta}f_{\mathrm{NL}}^{\Delta}}=\sum_{z_{i}}V(z_{i})\int\\!\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\frac{[6b_{1}b_{\phi}(z_{i})(kR_{*})^{\Delta}]^{2}P_{m}(k,\mu,z_{i})^{2}}{2k^{4}\mathcal{T}(k,z_{i})^{2}\,[b_{1}^{2}P_{m}(k,\mu,z_{i})+N(k,\mu,z_{i})]^{2}}\,.$
(3.4)

In the high signal-to-noise regime, $b_{1}^{2}P_{m}(k)\gg N$, this becomes

$F_{f_{\mathrm{NL}}^{\Delta}f_{\mathrm{NL}}^{\Delta}}\approx\frac{9}{\pi^{2}}\sum_{z_{i}}\frac{b_{\phi}^{2}(z_{i})}{b_{1}^{2}(z_{i})}R_{*}^{2\Delta}V(z_{i})\int\\!\mathrm{d}k\,\frac{k^{2\Delta-2}}{\mathcal{T}(k,z_{i})^{2}}\,.$
(3.5)

For $2\Delta<1$, the integral over $k$ is dominated by $k_{\mathrm{min}}$ and
the Fisher information arises from the smallest scales. However, due to the
limits placed by noise and volume, realistic surveys are far from the
asymptotic $k\to 0$ regime. For modes with the largest signal to noise,
$k\approx
k_{\mathrm{eq}}\approx${10}^{-2}\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$,
the matter power spectrum is flat, $P_{m}(k\sim k_{\mathrm{eq}})\propto k^{0}$
and, therefore, $\mathcal{T}(k,z)^{2}\sim k^{-1}$. In this regime, repeating
our analysis of (3.5) shows that the integral is dominated by large $k$ for
$\Delta\gtrsim 0$. In this regard, we expect that the transition from low to
high $k$ is gradual as we change $\Delta$.

This behavior is seen in forecasts by varying $k_{\mathrm{min}}$ as
illustrated in Fig. 5.

Figure 5: The information density with respect to $\log k_{\mathrm{min}}$ in
the billion-object survey, with fixed $\Lambda$CDM cosmology, and the same
$k_{\mathrm{min}}$, $\Delta
k=$0.003\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ and
$k_{\mathrm{max}}=$0.2\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
in each redshift bin. The different colors indicate which bias parameters are
marginalized when calculating $\sigma(f_{\mathrm{NL}}^{\Delta})$ using the
double-tracer version of the survey. A negative (positive) slope of the curves
indicates that there is more (less) information on smaller wavenumbers $k$. We
observe that the information density is dominated by small (large) wavenumbers
for small (large) scaling exponents.

Since increasing the minimum wavenumber $k_{\mathrm{min}}$ removes
information, the constraints $\sigma(f_{\mathrm{NL}}^{\Delta})$ always
decrease. We therefore show
$\mathrm{d}\sigma^{-2}(f_{\mathrm{NL}}^{\Delta})/\mathrm{d}\hskip-1.0pt\log
k_{\mathrm{min}}$ because this quantity provides a more quantitative
illustration of how much information is being lost as we change
$k_{\mathrm{min}}$.101010We consider the logarithmic derivative in order to
numerically extract the relevant scaling that we analytically discussed in the
previous paragraph. We additionally note that $\log k_{\mathrm{min}}$ is
implicitly normalized by
$1\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$. For $\Delta=0$,
the largest impact is for the smallest values of $k_{\mathrm{min}}$,
indicating that the lowest wavenumber has the most information. As we increase
$\Delta>0$, we see that this trend changes: the curves are mostly flat at
small $k_{\mathrm{min}}$ for $\Delta\lesssim 0.5$. This is a reflection of the
fact that there is no sharp transition in $\Delta$ due to the impact of
$\mathcal{T}(k,z)$ in (3.5).

As we move from large to small scales (or small to large $k$), we are
increasingly sensitive to the marginalization over the bias parameters. The
information density in Fig. 5 at large wavenumbers becomes highly suppressed
as we marginalize so that the true constraining power remains at small
wavenumbers. The oscillatory behavior in these figures highlights that the
baryon acoustic oscillations contained in the transfer function are imprinted
in the non-Gaussian signal and are not absorbed into the bias expansion.

At the same time, this short-distance information remains however crucial for
breaking degeneracies between bias and cosmological parameters. This is shown
in Fig. 6,

Figure 6: Dependence of $\sigma(f_{\mathrm{NL}}^{\Delta})$ on the maximum
wavenumber where we take $k_{\mathrm{max}}=k_{\mathrm{halo}}$ and
$k_{\mathrm{max}}=k_{\mathrm{NL}}(z)$, with the latter being less conservative
at higher redshifts. Here, we combine the middle two redshift bins ($2\leq
z\leq 3$) of the double-tracer billion-object survey into a single comoving
box with a total of about $1.3\text{\times}{10}^{8}$ objects. We show this
dependence of the constraining power for four biasing models. The choice of
$k_{\mathrm{max}}$ barely impacts $\sigma(f_{\mathrm{NL}}^{\Delta})$ at small
$\Delta$, but becomes significant for large $\Delta$. We explain the origin of
the feature around $\Delta=1$ in the main text.

which illustrates how the maximum wavenumber $k_{\mathrm{max}}$ impacts the
constraining power. We see that the constraint on $f_{\mathrm{NL}}^{\Delta}$
improves by roughly a factor of five to nine (depending on the biasing model)
for $\Delta=2$ when we choose the less conservative
$k_{\mathrm{max}}=k_{\mathrm{NL}}$. This improvement is consistent with the
number of additional modes available between $k_{\mathrm{halo}}$ and
$k_{\mathrm{NL}}$, and the fact that our signal depends on wavenumbers. This
difference becomes even more pronounced at higher redshifts since
$k_{\mathrm{NL}}(z)$ increases significantly, while $k_{\mathrm{halo}}$ is
independent of redshifts. For all of our other forecasts, we however use the
more conservative choice of
$k_{\mathrm{max}}(z)=\min\\{k_{\mathrm{halo}},k_{\mathrm{NL}}(z)\\}$ as noted
in §3.1.

When marginalizing over $b_{1}$, we notice a feature around $\Delta\approx 1$
in Fig. 6 and not around $\Delta=0.5$ as we might expect from (3.5). This is a
consequence of the information on $f_{\mathrm{NL}}^{\Delta}$ moving to larger
$k$ where $T(k)\neq 1$. Specifically, if we assume $P_{m}(k)\propto k^{n}$ for
large $k$, then $k^{2}T(k)\propto k^{n/2+3/2}$, which implies

$P_{gg}(k)\approx\left(b_{1}+2b_{\phi}f_{\mathrm{NL}}^{\Delta}C(kR_{*})^{\Delta-\frac{n}{2}-\frac{3}{2}}+\ldots\right)^{\\!2}P_{m}(k)\,,$
(3.6)

for some constant $C$. Using $n\approx-1$ for
$k\sim$0.1\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ [96], we
notice that $\Delta-\frac{3}{2}-\frac{n}{2}\approx\Delta-1$. At large
wavenumbers, the “scale-dependent bias” for $\Delta\approx 1$ is therefore
degenerate with the linear bias $b_{1}$, resulting in the observed feature.

#### Dependence on Fiducial Biases

All else being equal, a single-tracer survey with the largest possible
(absolute) value of $b_{\phi}\propto(b_{1}-1)$ will yield the most sensitive
measurement of $f_{\mathrm{NL}}$. It may therefore seem self-evident that
selecting highly biased targets is a central tool in the search for PNG.
Having said that, in practice, target selection involves numerous factors
which lie beyond the scope of these simple forecasts. Yet, when it comes to
understanding the performance of a given survey, the biases of the objects in
their sample will strongly influence the overall sensitivity. As a result, it
is important to separately understand the role that the fiducial biases play
from the aspects of the survey that we can control more directly.

For a single-tracer analysis, we only measure
$b_{\phi}f_{\mathrm{NL}}^{\Delta}$. At fixed redshift and assuming
$b_{\phi}\propto(b_{1}-1)$ from (2.22), we therefore have

$\sigma(f_{\mathrm{NL}}^{\Delta})_{\mathrm{single}}^{\text{fixed-$z$}}\propto\frac{1}{b_{1}-1}\,.$
(3.7)

As a result, we can expect a large improvement in a single-tracer analysis
from choosing a sample with $b_{1}\gtrsim 2$. Figure 7

Figure 7: Dependence of $\sigma(f_{\mathrm{NL}}^{\Delta})$ on the fiducial
value of the linear bias for both the single-tracer (solid) and double-tracer
(dashed) versions of the billion-object survey for $\Delta=0$ (left) and
$\Delta=1$ (right). We show this dependence as a function of the linear bias
$b_{1}(z=0)$ of the single-tracer survey in which case we see the improvement
in the constraining power of the survey for larger $b_{1}$. In the double-
tracer scenario, for which the galaxy sample is split into two populations of
equal number density with $b_{1}^{(A,B)}(z=0)=b_{1}(z=0)\pm 0.4$, the effects
of sample-variance cancellation can be observed. (Note that the $y$-axis is on
a linear and not on a logarithmic scale like in the previous figures.)

shows how the value of $b_{1}(z=0)$ influences the overall constraining power
of our billion-object survey in both single- and multi-tracer analyses. We see
that $\sigma(f_{\mathrm{NL}})$ increases sharply when the bias approaches one.
It does not diverge because of the redshift dependence of $b_{1}(z)$ which
ensures that some of the redshift bins have $b_{\phi}\neq 0$. On the other
hand, sample-variance cancellation in the multi-tracer configuration results
in a much weaker dependence of $\sigma(f_{\mathrm{NL}}^{\Delta})$ on the
linear bias and, in fact, benefits from having one of the tracers with
$b_{1}^{B}(z=0)<1$. At the same time, we however see that taking larger values
of $b_{1}(z)$ in a single-tracer analysis closes the gap to the double-tracer
case for $b_{1}(z=0)\gtrsim 3$.

Technically speaking, we cannot choose the linear bias $b_{1}(z)$ and the
number densities $\bar{n}_{g}(z)$ independently, in particular since more
highly biased objects are more rare. In practice, traditional spectroscopic
surveys are typically limited by the number of spectra that can be measured
and not the number of objects in the universe. However, other types of
surveying techniques, such as line intensity mapping (see e.g. [97, 98, 99,
100, 101, 102] in the context of PNG), may get closer to these fundamental
limits. We refer to [48] for a discussion of the sensitivities of surveys in
the regime where the galaxy power spectrum analysis is limited by the number
of available objects.

#### Dependence on the Volume and Redshift

Given a single redshift bin with volume $V$ and a fixed choice of targets, it
is easy to see the impact of this choice on $k_{\mathrm{min}}$ and
$k_{\mathrm{max}}$. However, when we design a survey, the choice of redshift
range, volume and targets is not fixed, but is part of the survey design. The
redshift range and volume affect $k_{\mathrm{min}}$ and $k_{\mathrm{max}}$,
but they also change the biases of the targets and the number densities of
objects in each redshift bin. For the purpose of designing a survey, we
therefore want to understand the optimal choices for the sky fraction
${f_{\mathrm{sky}}}$ and the redshift range as given by $z_{\mathrm{max}}$.

It is important to compare survey strategies assuming “constant effort”. In
practice, it is the time for acquiring spectra that limits the sensitivity of
a spectroscopic survey. This is why we vary ${f_{\mathrm{sky}}}$ for a given
redshift range while holding the total number of objects fixed. As a result,
there is a trade-off between increasing the volume and increasing the shot
noise, which needs to be optimized. When it comes to varying the redshift
range, the time required to observe each object changes with redshift which
means that the total number of objects in the sample decreases with larger
$z_{\mathrm{max}}$. To hold the observational effort fixed, we use a
simplified model where it takes twice the amount of time to obtain spectra for
galaxies at $z>3$ than for galaxies at $z<3$ (cf. [11]). We further assume
that it is desirable to maintain a constant level of shot noise over the
survey volume or that we obtain spectra evenly across the sky and along the
radial direction.

Higher-redshifts objects are generally expected to yield larger values of
$b_{1}(z)$ [see the discussion around (2.7)]. Since
$b_{\phi}(z)\propto(b_{1}(z)-1)$ from the universality relation (2.22), the
signal can be enhanced by a large amount by extending the survey to redshifts
where $b_{1}(z)>2$ instead of $b_{1}(z)\approx 1$. Assuming the bias evolution
model of (2.7), we see in Fig. 8

Figure 8: Dependence of $\sigma(f_{\mathrm{NL}}^{\Delta=1})$ on the sky
fraction ${f_{\mathrm{sky}}}$ and the maximum redshift $z_{\mathrm{max}}$ at
fixed observational effort for the double-tracer billion-object survey. The
left panel shows the fiducial case, $b_{1}(z)\propto D(z)^{-1}$ [cf. (2.7)]
and the right panel displays the constant bias case,
$b_{1}(z)=\mathrm{const}$. While we only show the results for $\Delta=1$, the
displayed trends are similar for all scaling exponents $\Delta$.

that we maximize our sensitivity to $f_{\mathrm{NL}}^{\Delta}$ even at
constant effort by maximizing $z_{\mathrm{max}}$. While we only show the
results for $\Delta=1$, the qualitative features in the figure are independent
of the scaling exponent $\Delta$. This implies that increasing the redshift
range at fixed observational effort increases the constraining power for the
entire range of $\Delta$.

Increasing the redshift range also increases the number of modes. If we
however hold $b_{1}(z)=b_{1}$ fixed and increase $z_{\mathrm{max}}$, we see in
the right panel of Fig. 8 that there is a more complex relationship between
$z_{\mathrm{max}}$ and the sensitivity to $f_{\mathrm{NL}}^{\Delta}$. For
instance, keeping $z_{\mathrm{max}}=3$ and increasing ${f_{\mathrm{sky}}}$
also increases the number of modes without the added observing time per object
needed at high redshifts. Having said that, the benefit of a large sky
fraction is much weaker than the impact of larger $b_{1}$ at high $z$ which
can be seen by the lower overall sensitivity to $f_{\mathrm{NL}}^{\Delta}$ in
the right panel.

By comparing the two panels in Fig. 8, we can see that the optimal
configuration is still at the largest maximum redshift and sky fraction,
$z_{\mathrm{max}}=5$ and ${f_{\mathrm{sky}}}=0.7$, for either biasing model.
For smaller ${f_{\mathrm{sky}}}$, the large advantage offered by larger
$z_{\mathrm{max}}$ however disappears when going from evolving to fixed bias
due to the lower overall number density of galaxies at higher redshifts. In
this case, more emphasis is placed on high-redshift targets without the
corresponding boost in the signal. (The drop-off in $\sigma(f_{\mathrm{NL}})$
between $z_{\mathrm{max}}=3.0\text{ and }3.5$ is the result of our simplified
choice of dividing the observational effort at $z=3$.)

To summarize, at fixed observational effort, increasing the redshift range via
a larger $z_{\mathrm{max}}$ has the largest influence on the sensitivity of a
survey. The reason is that $z_{\mathrm{max}}$ increases both the size of the
signal and the number of modes. The former is particularly important since it
offsets the increased noise associated with a large-volume survey. We however
note that this conclusion does depend on the redshift dependence of the bias
and the observing time needed for acquiring spectra of these high-redshift
objects.

#### Photometric versus Spectroscopic Surveys

The benefit of a three-dimensional survey can simply be estimated by counting
modes for many cosmological parameters, including $f_{\mathrm{NL}}$ from the
bispectrum. The number of modes in a three-dimensional survey scales as
$k_{\mathrm{max}}^{3}$, while it scales as $k_{\mathrm{max}}^{2}$ for a two-
dimensional survey. The benefits of a three-dimensional survey are however
less clear for scale-dependent bias since the information for this signal
manifests itself at low $k$, in particular for $\Delta\lesssim 0.5$. In
reality, the trade-off is of course not between a two- and a three-dimensional
survey, but between surveys which measure redshifts spectroscopically and
photometrically. While photometric redshifts are less precise, they may be
good enough to provide three-dimensional information on the scales needed for
scale-dependent bias for a significant range of $\Delta$ which we investigate
now.111111Spectroscopic surveys will have additional benefits for controlling
systematics, such as projecting out large-scale effects, but typically at the
cost of larger shot noise. Quantifying the importance of spectra in this
context is a delicate issue which is beyond the scope of this work.

The impact of redshift errors on the constraining power of the billion-object
survey is shown in Fig. 9.

Figure 9: Comparison of $\sigma(f_{\mathrm{NL}}^{\Delta})$ for a range of
photometric redshift errors $\sigma_{z0}\in[0.001,0.2]$ for the billion-object
survey. We observe that photometric redshift errors only mildly degrade the
constraining power of a survey for $\Delta\lesssim 1.3$. Photometric surveys
therefore seem to be able to access most of the information about the scale-
dependent bias with the additional benefit of generally lower noise levels in
the galaxy power spectrum.

We see, somewhat surprisingly, that spectroscopic surveys do not provide a
significant improvement to the power spectrum measurement of
$f_{\mathrm{NL}}^{\Delta}$ for $\Delta\lesssim 1-1.3$. (The specific value of
this transition in $\Delta$ depends on the employed biasing model.) In
practice, photometric redshifts therefore appear to be good enough to access a
lot of the available information about the scale-dependent bias. This
conclusion is supported by the similarity of the forecasts for SPHEREx and
LSST. While the two surveys pursue very different strategies for obtaining
redshifts, the resulting forecasts are quite similar. At first sight, this
similarity is still surprising because it does not reflect the impact of the
difference in number density shown in Fig. 4. The reason for this, however, is
just that multi-tracer analyses are not implemented in conventional LSST
forecasts. This naturally raises the question whether a multi-tracer analysis
would lead to an advantage of LSST over some spectroscopic surveys given its
very high number density (cf. Fig. 3). We will explore this in detail in
Section 5.

Importantly, we have not discussed the impact of systematic effects [103] on
photometric and spectroscopic surveys. Given that the signals are the largest
on large angular scales (especially for $\Delta\lesssim 0.5$), a variety of
atmospheric, detector and astrophysical effects can contaminate our signal.
Having said that, we expect that many of the conclusions drawn from these
forecasts will be robust over a large parameter range since the scales on
which the signals dominate vary with the scaling exponent $\Delta$.

##  4 Constraints and Data Analysis for General Scaling Exponents

In the previous section, we focused on forecasting current and future
constraints on $f_{\mathrm{NL}}^{\Delta}$, and on how to optimize a survey to
measure PNG from the galaxy power spectrum for $\Delta\in[0,2]$. In this
section, we will apply this knowledge to constrain $f_{\mathrm{NL}}^{\Delta}$
over this same range of scaling exponents $\Delta$ using the BOSS DR12 galaxy
power spectra. We will first derive and propose an effective and convenient
search strategy. We will then perform this measurement and compare it to a
direct analysis.

###  4.1 Search Strategy

The scaling exponent $\Delta$ is fixed in some inflationary models, while it
is a free parameter in others. For the standard local and equilateral shapes,
the value of $\Delta$ are fixed by their templates to $\Delta=0$ and
$\Delta=2$, respectively. The physics of these cases is clear since they arise
from multi- and single-field dynamics. In quasi-single-field inflation, we
have $\Delta=3/2-\sqrt{9/4-m^{2}/H^{2}}$ which implies that it is natural to
scan over values of $\Delta\in[0,3/2]$ in an analysis of those models. The
natural question we face in a practical analysis of current data therefore is
whether it is prudent to vary $f_{\mathrm{NL}}^{\Delta}$ and $\Delta\in[0,2]$
as free parameters, or whether we should hold the scaling exponent fixed to
some value as suggested by a given model, $\Delta=\Delta^{\prime}$, and only
vary $f_{\mathrm{NL}}^{\Delta^{\prime}}$.

The strategy we will follow is to hold $\Delta$ fixed to a few discrete
values, e.g. $\Delta=0,2$ or $\Delta=0.0,0.1,\ldots,1.9,2.0$, and infer the
associated constraints on $f_{\mathrm{NL}}^{\Delta}$ at each respective
$\Delta$. The reason why we prefer this strategy over scanning over both
parameters simultaneously is that small changes to $\Delta$ have no impact on
the signal to noise and, therefore, represent a degenerate direction in the
analysis. This becomes evident in the Fisher matrix formalism, where we have
for a fiducial value of $f_{\mathrm{NL}}^{\Delta}=0$:

$\left.\frac{\partial}{\partial\Delta}P_{gg}(k)\right|_{f_{\mathrm{NL}}^{\Delta}=0}=6f_{\mathrm{NL}}^{\Delta}b_{\phi}(k)\,k^{-2}\mathcal{T}(k)^{-1}\log(kR_{\star})(kR_{\star})^{\Delta}\,b(k)P_{\mathrm{lin}}(k)\Big{|}_{f_{\mathrm{NL}}^{\Delta}=0}=0\,.$
(4.1)

We therefore have no reason to expect marginalizing over $\Delta$ would lead
to more accurate results until our data are comfortably excluding
$f_{\mathrm{NL}}^{\Delta}=0$.

We can quantify to what degree the data can distinguish two values of the
scaling exponent, $\Delta_{1}$ and $\Delta_{2}$, through the “cosine” between
the signals,

$\cos{(f_{\mathrm{NL}}^{\Delta_{1}},f_{\mathrm{NL}}^{\Delta_{2}})}=\frac{F_{f_{\mathrm{NL}}^{\Delta_{1}}f_{\mathrm{NL}}^{\Delta_{2}}}}{\sqrt{F_{f_{\mathrm{NL}}^{\Delta_{1}}f_{\mathrm{NL}}^{\Delta_{1}}}F_{f_{\mathrm{NL}}^{\Delta_{2}}f_{\mathrm{NL}}^{\Delta_{2}}}}}\equiv\frac{F_{12}}{\sqrt{F_{11}F_{22}}}\,,$
(4.2)

where we defined $F_{ij}\equiv
F_{f_{\mathrm{NL}}^{\Delta_{i}},f_{\mathrm{NL}}^{\Delta_{j}}}$. This is
precisely the same definition of the cosine used to define the PNG shapes in a
bispectrum analysis [104]. We can also generalize this to the Fisher matrix
after marginalizing over the bias parameters by inverting the Fisher matrix,
truncating to $f_{\mathrm{NL}}^{\Delta_{1}}$ and
$f_{\mathrm{NL}}^{\Delta_{2}}$, and inverting the truncated matrix to get the
effective Fisher matrix for these two parameters. We therefore take $F_{ij}$
to be the marginalized, two-dimensional Fisher matrix for the parameters
$f_{\mathrm{NL}}^{\Delta_{1}}$ and $f_{\mathrm{NL}}^{\Delta_{2}}$ from now on.

The correlation (4.2) between the measurements of
$f_{\mathrm{NL}}^{\Delta_{1}}$ and $f_{\mathrm{NL}}^{\Delta_{2}}$ is shown in
Fig. 10

Figure 10: Correlation matrix for measurements of galaxy power spectra with
different values of the non-Gaussian exponent $\Delta$ as defined in (4.2).
This was computed for BOSS with fiducial values of
$f_{\mathrm{NL}}^{\Delta_{i}}=0$, fixed $\Lambda$CDM parameters and
marginalized over all biases. The orange arrows indicate the coverage of the
local ($\Delta=0$) and equilateral ($\Delta=2$) templates with a correlation
coefficient of larger than $0.9$.

using the experimental specifications for the BOSS survey of Table 3 which
have yielded forecasted results consistent with performed data analyses in the
past (see e.g. [14]).121212Note that for the calculation of the correlation
matrix for BOSS, we assumed fiducial values of $b_{s^{2}}$ and
$b_{\Pi\Pi^{[2]}}$ based on the Lagrangian local-in-matter-density biasing
model [52] and $b_{\delta^{2}}$ based on halo simulation fit in [56]. This
choice (rather than setting them to $0$ as in most other cases in this paper)
leads to minimal differences in the correlation coefficient. In addition, we
impose the same Gaussian priors on the loop biases as quoted in [105] which
barely affect the forecasted constraints for small $\Delta$ and up to a factor
of four at large $\Delta$, as expected. Finally, we employ a slightly larger
minimum wavenumber
$k_{\mathrm{min}}=$0.01\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$.
All these choices are guided by the data analysis performed in §4.2. The
cosine computed for other surveys is similar (see Appendix B). The key take-
away from this correlation is that we can effectively constrain almost the
entire range $\Delta\in[0,2]$ using measurements only at $\Delta=0\text{ and
}2$. Specifically, using the Cauchy-Schwartz inequality and the Cramér-Rao
(CR) bound, we have

$(F_{22})^{-1}\leq\frac{F_{11}}{F_{12}^{2}}\leq\frac{\sigma^{2}(f_{\mathrm{NL}}^{\Delta_{1}})}{\cos^{2}(f_{\mathrm{NL}}^{\Delta_{1}},f_{\mathrm{NL}}^{\Delta_{2}})}\frac{F_{11}}{F_{22}}\,,$
(4.3)

where $\sigma^{2}(f_{\mathrm{NL}}^{\Delta_{1}})$ is the variance of the
measurement of $f_{\mathrm{NL}}^{\Delta_{1}}$ while holding
$f_{\mathrm{NL}}^{\Delta_{2}}$ fixed. The best possible constraint we could
place on $f_{\mathrm{NL}}^{\Delta_{2}}$,
$\sigma^{2}(f_{\mathrm{NL}}^{\Delta_{2}})_{\mathrm{CR}}=F_{22}^{-1}$, is
therefore bounded by the measurement of $f_{\mathrm{NL}}^{\Delta_{1}}$ via

$\sigma^{2}(f_{\mathrm{NL}}^{\Delta_{2}})_{\mathrm{CR}}\leq\frac{\sigma^{2}(f_{\mathrm{NL}}^{\Delta_{1}})}{\cos^{2}(f_{\mathrm{NL}}^{\Delta_{1}},f_{\mathrm{NL}}^{\Delta_{2}})}\frac{F_{11}}{F_{22}}\,.$
(4.4)

As a result, we can place an approximate constraint on
$f_{\mathrm{NL}}^{\Delta_{2}}$ by rescaling the measurement of
$f_{\mathrm{NL}}^{\Delta_{1}}$ using the elements of the marginalized Fisher
matrix $F_{ij}$.

It is important to highlight that this strategy is only effective if our
parameter inference is consistent with $f_{\mathrm{NL}}^{\Delta}=0$. To be
more precise, all values of $\Delta$ look similar at low signal to noise for
$f_{\mathrm{NL}}^{\Delta}$. On the other hand, the data is no longer
independent of $\Delta$ if one or more of your $f_{\mathrm{NL}}^{\Delta}$
measurements at specific $\Delta$ show significant evidence for
$f_{\mathrm{NL}}^{\Delta}\neq 0$. In that case, a measurement at a different
value $\Delta^{\prime}$ could produce a much larger or much smaller signal to
noise. We refer to Appendix B for more details on this scenario.

###  4.2 BOSS DR12 Analysis

Based on the strategy laid out above, we inferred constraints on
$f_{\mathrm{NL}}^{\Delta}$ from the BOSS DR12 dataset following the general
setup of the analysis for local, orthogonal and equilateral PNG performed in
[106, 36] (see also [35]) with a few different (generally more conservative)
choices. In the following, we provide a brief overview of our data analysis,
describe these differences and report the limits on PNG that we inferred
directly and via the correlation matrix (4.2).

We analyzed the power spectrum from galaxy clustering data of the twelfth and
final release of BOSS, referred to as DR12 [49]. This dataset contains the
positions of about $1.2\text{\times}{10}^{6}$ galaxies between redshifts $0.2$
and $0.75$ in a cosmic volume of roughly
$5.8\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$ divided
into four subsets: the Northern and Southern galactic caps of the two non-
overlapping redshift bins $[0.2,0.5]$ and $[0.5,0.75]$ with effective mean
redshifts $z=0.38\text{ and }0.61$. We employ the galaxy power spectrum
multipoles $P_{gg}^{(\ell)}\hskip-1.0pt(k)$, with $\ell=0,2,4$, measured using
a quadratic window-function-free estimator [107]. The corresponding covariance
matrices were computed from MultiDark-Patchy mock catalogs [108]. This follows
the analyses in [106, 36] restricted to the power spectrum multipoles. We
model the nonlinear galaxy power spectrum in redshift space as implemented in
CLASS-PT [109] (based on CLASS [110]) which follows the analyses of [111, 112]
based on the effective field theory of large-scale structure (EFTofLSS) [113,
114, 115, 116] (see [17, 117] for recent reviews) and galaxy bias expansion
[52], including the Alcock-Paczynski effect. We additionally include the non-
Gaussian contributions from the scale-dependent bias (2.21) due to
$f_{\mathrm{NL}}^{\Delta}\neq 0$ similar to the analyses in [106, 36].131313In
this paper, we have not considered the imprints of light inflationary fields
in the bispectrum, but focused on the scale-dependent bias in the power
spectrum. To remain model agnostic, we also did not include the contribution
of a non-zero primordial bispectrum to the galaxy power spectrum.

We performed a Markov chain Monte Carlo (MCMC) analysis of the data using the
Metropolis-Hastings sampler of MontePython [118, 119], varying the theoretical
galaxy power spectrum multipoles in each step. We fix the $\Lambda$CDM
parameters to the Planck 2018 best-fit values [58] (with the sum of neutrino
masses $\sum m_{\nu}=$0.06\text{\,}\mathrm{e}\mathrm{V}$$) instead of
including Planck information on $\Lambda$CDM since we saw essentially no
change in the forecasted constraints on $f_{\mathrm{NL}}^{\Delta}$ of Section
3 when marginalizing over all of our galaxy bias parameters. Per redshift bin
and galactic cap, we however separately vary the bias parameters
$b_{1}^{(i)}$, $b_{2}^{(i)}$ and $b_{\mathcal{G}_{2}}^{(i)}$ as defined in
[120, 121],141414The bias expansions used in our forecasts and data analysis
are equivalent. and analytically marginalize over the other eight biases and
EFTofLSS counterterms which appear linearly in the theoretical galaxy power
spectrum [122]. We impose flat priors on the linear biases
$b_{1}^{(i)}\in[1,4]$ and the non-Gaussian amplitude
$f_{\mathrm{NL}}^{\Delta}$ (infinitely wide), and Gaussian priors on the other
nuisance parameters as described in [105]. Following [106], we also
marginalize over the non-Gaussian bias $b_{\phi}\to N_{b_{\phi}}b_{\phi}$
separately for each data subset and impose a wide Gaussian prior on its
normalization, $N_{b_{\phi}}\sim\mathcal{N}(1,5)$, which is motivated by the
peak-background split model [78].151515It might appear surprising that we
marginalize over $b_{\phi}$ with a Gaussian prior given the degeneracy with
$f_{\mathrm{NL}}^{\Delta}$. Note that the mean value of the Gaussian is still
set by the universality relation (2.22) with $p=1$, but its variance is very
wide. In addition, we marginalize over $b_{\phi}$ in each of the four data
subsets separately, as we do with the Gaussian bias parameters. This procedure
alleviates degeneracies between $f_{\mathrm{NL}}^{\Delta}$ and $b_{1}$ which
can be severe if we assume that the universality relation is exact, especially
for $\Delta\sim 1$. This procedure is an imperfect solution to address the
uncertainty in $b_{\phi}$ which we plan to address as part of [123]. On the
other hand, the non-Gaussian parameter $f_{\mathrm{NL}}^{\Delta}$ is commonly
varied for all four data subsets at fixed exponent $\Delta$. We (more)
conservatively limit the range of wavenumbers to
$k\in$0.010.13\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$\text{
and }$0.010.16\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ for
the low- and high-redshift bin, respectively, after initial tests. We
therefore employ a smaller maximum wavenumber than in [36, 106], who tested
the validity of the analysis for $\Delta=0\text{ and }2$ on mock catalogs, and
a slightly higher (same) minimum (maximum) wavenumber than the forecasts shown
in Fig. 3. Since $\Delta
k=$0.005\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$, this
choice implies that we employ 24 and 30 $k$-bins for the respective redshift
bins, each power spectrum multipole and both galactic caps. All chains
converged with a Gelman-Rubin criterion of $R-1<0.01$ (usually much smaller)
for each parameter.

The results of our separate MCMC data analyses for 21 fixed values of the non-
Gaussian exponent $\Delta$ are shown in Fig. 11.

Figure 11: Comparison of the measured (diamonds) and derived (lines) values
for twice the standard deviation of $f_{\mathrm{NL}}^{\Delta}$,
$2\sigma(f_{\mathrm{NL}}^{\Delta})$, for the BOSS DR12 dataset. The light
(dark) blue lines show the approximate value of
$2\sigma(f_{\mathrm{NL}}^{\Delta})$ derived from the data with
$\Delta=2\leavevmode\nobreak\ (0)$, i.e. the equilateral (local) PNG templates
using (4.5), with the minimum of these being shown in blue. The diamonds
indicate the constraints directly inferred from BOSS DR12 for fixed values of
$\Delta$, while the gray pentagons are the corresponding mean values. Filled
(unfilled) symbols and solid (dotted) lines represent positive or upper
(negative or lower) values. We do not find evidence for
$f_{\mathrm{NL}}^{\Delta}\neq 0$ for any $\Delta\in[0,2]$, and generally good
agreement between the measured and derived limits on general PNG. This
demonstrates that our proposed search strategy works on data.

Since the posterior distributions for $f_{\mathrm{NL}}^{\Delta}$ are slightly
non-Gaussian, we display the mean values in gray pentagons, and the upper and
lower $2\sigma$ constraints in filled and unfilled black diamonds,
respectively. We do not find any evidence for any PNG shape with
$\Delta\in[0,2]$ and observe the same characteristic functional dependence of
$2\sigma(f_{\mathrm{NL}}^{\Delta})$ on $\Delta$ as we forecasted in Fig.
3.161616The non-Gaussian amplitude $f_{\mathrm{NL}}^{\Delta}$ is defined
relative to the scalar amplitude,
$A_{\mathrm{s}}^{1/2}\approx$5\text{\times}{10}^{-5}$$ [58], such that weak
PNG corresponds to $f_{\mathrm{NL}}A_{\mathrm{s}}^{1/2}\ll 1$. This means that
our bounds $2\sigma(f_{\mathrm{NL}}^{\Delta})\sim${10}^{4}$$ are approaching
the limits of the range of validity of our implicit expansion in linear order
in $f_{\mathrm{NL}}^{\Delta}$. This concern is however alleviated even with a
very conservative inclusion of bispectrum data [106, 35]. We note that the
difference in forecasted constraining power of Fig. 3 and the data analysis of
Fig. 11 can be mainly attributed to the use of Gaussian priors on the higher-
order galaxy bias parameters in Section 4. In fact, we find reasonable
consistency if we include the priors on the Gaussian biases in our Fisher
matrices as discussed in §4.1. Comparing to [106, 36, 35], our constraints are
weaker, as expected due to our more conservative approach and use of less
data, but consistent. This is particularly evident for the case of the
equilateral shape ($\Delta=2$) that dramatically benefits from bispectrum
information in various ways which is absent in our analysis of the power
spectrum alone. On the other hand, our limits on $f_{\mathrm{NL}}^{\Delta}$
are comparable to the constraints on $f_{\mathrm{NL}}^{\mathrm{eq}}$ when
including bispectrum information [106] for $\Delta\lesssim 0.75$. While the
bispectrum will always add information, we might however imagine there is a
more limited gain in this regime.

We can also compare our bounds on $f_{\mathrm{NL}}^{\Delta}$ to those obtained
in previous analyses of CMB and LSS data. As already anticipated and explained
in Section 3, the CMB limits on local and equilateral PNG from the Planck 2018
bispectra are significantly better [6]. This is also the case for the
constraints for $\Delta\in[0,3/2]$ derived from Planck 2013 data [124],
$|f_{\mathrm{NL}}^{\Delta}|\lesssim 50$, which is dominated by the equilateral
bispectrum while our bounds arise from the squeezed limit of the quasi-single-
field shapes. The only previous analysis constraining
$f_{\mathrm{NL}}^{\Delta}$ from LSS data was performed in [50] on SDSS-III DR8
data with more limited modeling and range of scales. It is therefore
unsurprising that our bounds are considerably more constraining.

In addition to directly inferring constraints on $f_{\mathrm{NL}}^{\Delta}$
from the data, there exists a second path to inferring the bounds at a given
value of $\Delta$ based on the correlation of different values of $\Delta$ in
the marginalized galaxy power spectrum as discussed in §4.1. We take
$\Delta=\Delta_{1}=0\text{ and }2$ as the fiducial values for this inference
since these scaling exponents correspond to the signal from local and
equilateral PNG, respectively, which are usually considered in data analyses.
Based on (4.4), we obtain these derived bounds as follows:

$\left[2\sigma(f_{\mathrm{NL}}^{\Delta})\right]_{\mathrm{MCMC}}=\left[\frac{1}{\cos(f_{\mathrm{NL}}^{\Delta},f_{\mathrm{NL}}^{\Delta=0,2})}\frac{\sigma(f_{\mathrm{NL}}^{\Delta})}{\sigma(f_{\mathrm{NL}}^{\Delta=0,2})}\right]_{\mathrm{Fisher}}\left[2\sigma(f_{\mathrm{NL}}^{\Delta=0,2})\right]_{\mathrm{MCMC}}\,,$
(4.5)

where the subscript ‘Fisher’ (‘MCMC’) indicates that these quantities are
calculated based on the Fisher matrix (from the MCMC analysis). The
correlation term $\cos(f_{\mathrm{NL}}^{\Delta},f_{\mathrm{NL}}^{\Delta=0,2})$
is displayed in Fig. 10, while the standard deviations
$\sigma(f_{\mathrm{NL}}^{\Delta})$ are computed from the marginalized one-
dimensional Fisher matrix as in Section 3, but including the same Gaussian
priors as in the correlation term and the data analysis. The constraints
derived from the measurement at $\Delta=0\text{ and }2$ are displayed in Fig.
11 as the dark and light blue curves. We see that the constraint derived from
$2\sigma(f_{\mathrm{NL}}^{\Delta=0})$ dominates for $\Delta\lesssim 1.1$,
while the constraint derived from $2\sigma(f_{\mathrm{NL}}^{\Delta=2})$
dominates for larger $\Delta$, which is consistent with our expectations from
Fig. 10.

Overall, we find generally reasonable agreement between the directly inferred
limits on $f_{\mathrm{NL}}^{\Delta}$, and those jointly derived from the local
and equilateral measurements. This is especially true for those values of
$\Delta$ with large correlation and negligible mean values. This implies that
we can use (4.5), with the correlation term calculated in a Fisher matrix
approach for a given survey, to derive $2\sigma(f_{\mathrm{NL}}^{\Delta\neq
0,2})$ in the absence of a detection of
$f_{\mathrm{NL}}^{\mathrm{loc}},f_{\mathrm{NL}}^{\mathrm{eq}}\neq 0$. In other
words, we have demonstrated that the search strategy that we proposed in §4.1
works on BOSS data and allows to dramatically reduce the number of analyses
needed to constrain general forms of primordial non-Gaussianity and light
inflationary fields as parametrized by $\Delta$ from measurements of local and
equilateral non-Gaussianity. The agreement between our analysis and forecasts,
together with the consistency with prior analyses, is an additional sign that
searches for light fields with the galaxy power spectrum will be competitive
with the CMB for near-term galaxy surveys.

##  5 Multiple Tracers and the Dependence on Astrophysics

Sample-variance cancellation plays an important part in improving the
measurement of primordial non-Gaussianity from the galaxy power spectrum.
Applying this technique however requires knowledge of the details of the
specific galaxy samples used in the multi-tracer forecast or analysis. In this
section, we will first investigate how an optimal survey design can maximize
the scientific return of multiple tracers in spectroscopic surveys. Then, we
will consider the untapped potential of a multi-tracer configuration of LSST
and compare it with to SPHEREx.

###  5.1 Optimizing Spectroscopic Multi-Tracer Analyses

We now explore the potential of multi-tracer analyses and how astrophysical
details of the galaxy samples affect the constraints on
$f_{\mathrm{NL}}^{\Delta}$. Given our biasing model (2.6) with (2.18), this
comes down to the description of the non-Gaussian bias parameter
$b_{\phi}(z)$, which depends on three parameters when using the universality
relation (2.22): $\delta_{c}$, $b_{1}(z)$ and $p$. Since the parameter
$\delta_{c}$ is completely degenerate with $f_{\mathrm{NL}}^{\Delta}$, we will
focus on exploring what combination of $b_{1}(z)$ and $p$ leads to the biggest
improvement on $\sigma(f_{\mathrm{NL}}^{\Delta})$.

For galaxy power spectrum forecasts to recover all the information in the
primordial statistics, we need to apply the multi-tracer technique and its
ability for sample-variance cancellation at very high number density. For
example, we know the information in the primordial bispectrum for $\Delta>0$
is dominated by $k\sim k_{\mathrm{max}}$ [104]. On the other hand, the non-
Gaussian information in the single-tracer power-spectrum forecasts are limited
by sample variance at small wavenumbers for $\Delta\lesssim 0.5$. Sample-
variance cancellation eliminates this artificial dependence on the noise at
low $k$, limiting a measurement instead by the shot noise associated with
galaxy formation [125], which corresponds to modes at $k_{\mathrm{max}}\sim
k_{\mathrm{halo}}$.

Realistic spectroscopic surveys are unfortunately far from being in the fully
multi-tracer regime as we have seen in Fig. 4. As a result, it is harder to
get reliable intuition for how to optimize the measurement of
$f_{\mathrm{NL}}^{\Delta}$. We have previously explored how to maximize the
science of a spectroscopic survey by changing the experimental configuration
in terms of ${f_{\mathrm{sky}}}$, $z_{\mathrm{max}}$, etc. Here, we will
assume a fixed survey geometry and investigate what additional information can
be extracted by taking advantage of sample-variance cancellation.

In the fully multi-tracer regime, the inverse variance depends on the linear
bias $b_{1}^{(i)}$ and non-Gaussian bias $b_{\phi}^{(i)}$ as [126]

$\sigma(f_{\mathrm{NL}}^{\Delta})^{-2}\propto\left|b_{1}^{(1)}b_{\phi}^{(2)}-b_{1}^{(2)}b_{\phi}^{(1)}\right|.$
(5.1)

On the other hand, we simply maximize $b_{1}b_{\phi}$ in the single-tracer
regime. The question therefore is how we can choose the samples to minimize
$\sigma(f_{\mathrm{NL}}^{\Delta})$ in realistic surveys. In Figure
12,171717Figure 12 only shows the forecasted constraints and improvements for
the scaling exponent $\Delta=0$, i.e. local PNG. When we marginalize over all
biases, the scaling behavior of the scale-dependent bias however has little
effect on the improvement factor and the marginalized biasing model only
changes the overall constraining power as indicated by the color bar in the
upper panels. In other words, the overall qualitative behavior conveyed in
Fig. 12 remains intact for different biasing models and scaling exponents.

Figure 12: Double-tracer forecast for $\sigma(f_{\mathrm{NL}}^{\Delta})$ with
$\Delta=0$ for the billion-object survey (top) and relative comparison to a
single-tracer forecast,
$\sigma(f_{\mathrm{NL}}^{\Delta})_{\mathrm{double}}/\sigma(f_{\mathrm{NL}}^{\Delta})_{\mathrm{single}}$
(bottom), with fixed $\Lambda$CDM parameters. In the left column, we vary the
linear biases $b_{1}^{(i)}$ of the two samples at fixed
$p^{(i)}_{\vphantom{1}}=1$, i.e. similar to our forecasts in Section 3. In the
right column, we use the fiducial biases of the samples [cf. (2.7)], with
$\bar{b}_{1}^{(1)}(z=0)=2.0>1.2=\bar{b}_{1}^{(2)}(z=0)$, and vary the
parameter $p^{(i)}_{\vphantom{1}}$ in the universality relation for
$b_{\phi}^{(i)}$ [cf. (2.22)]. While we only show these results for
$\Delta=0$, the qualitative features of these panels remain the same for other
values of the scaling exponent.

we show the relative comparison of $\sigma(f_{\mathrm{NL}}^{\Delta})$ between
a double- and single-tracer forecast of the billion-object survey while
holding the $\Lambda$CDM parameters fixed.181818When marginalizing over the
biases, the difference between the forecasts for fixed and marginalized
$\Lambda$CDM parameters (with a Planck prior) is at most a few percent. We
therefore fix the $\Lambda$CDM parameters to their fiducial values for all
forecasts in this section. We show this for $\Delta=0$ (see also e.g. [127])
since we expect this case to benefit the most from a multi-tracer analysis,
but explicitly verified that the forecasts with $\Delta>0$ show at most
comparable levels of improvement.

We vary the linear biases of the two subsamples of the survey in the left
column of Fig. 12, while fixing $p=1$, as in all of our forecasts in Section
3. Note that we only vary the bias at $z=0$, as labeled in the figure, with
the biases of the actual galaxy samples then being calculated based on the
bias evolution model of (2.7). The values along the diagonal
$b_{1}^{(1)}=b_{1}^{(2)}$ correspond to a single-tracer forecast, in which
$\sigma(f_{\mathrm{NL}}^{\Delta})$ simply decreases with increasing
$b_{1}^{(i)}$. We observe that the multi-tracer improvement is most pronounced
here when the difference between $b_{1}^{(1)}$ and $b_{1}^{(2)}$ is maximized.
We additionally see that the improvement of the double- over the single-tracer
approach is limited to about 30%.

So far, we have assumed a universal halo mass function which corresponds to
$p=1$ in the universality relation (2.22). Even though this has been widely
adopted in the literature, this assumption is not a perfect description of
halos in simulations and will likely not be generally true for actual galaxy
samples, cf. e.g. [128, 129, 130, 131, 132, 133, 134, 135, 136, 137]. We
therefore now explore the possibility of improving the constraints on
$f_{\mathrm{NL}}^{\Delta}$ by selecting tracers within the parametrization of
$b_{\phi}$ through $p$ (assuming such knowledge is already reliably known from
simulation results).

We therefore vary the parameters $p^{(i)}$ of the two galaxy samples while
using their fiducial biases in the right column of Fig. 12. Note that the fact
that the linear biases for the two tracers are different, with
$b_{1}^{(1)}(z=0)=2.0$ and $b_{1}^{(2)}(z=0)=1.2$, is the reason for the
asymmetry of these panels. We observe the largest improvement for small
$p^{(1)}$ and large $p^{(2)}$ which is exactly when the combination shown in
(5.1) is maximized. In this case, the double-tracer approach can improve the
constraints by as much as a factor of two over a single-tracer analysis which
is much larger than what we saw for fixed $p^{(i)}$ and varying $b_{1}^{(i)}$
in the left panels. Figure 13

Figure 13: Forecasted constraints on $f_{\mathrm{NL}}^{\Delta}$ for different
values of the parameter $p^{(i)}$ in the universality relation (2.22) of
$b_{\phi}$ for the double-tracer configuration of the billion-object survey
and fixed $\Lambda$CDM parameters. We observe that
$\sigma(f_{\mathrm{NL}}^{\Delta})$ improves for all $\Delta$ when
$|b_{1}^{(1)}b_{\phi}^{(2)}-b_{1}^{(2)}b_{\phi}^{(1)}|$ increases which is the
direction indicated by the arrow.

displays the constraining power as a function of the scaling exponent $\Delta$
for a few pairs of tracers characterized by different values of
$p^{(i)}$.191919We also tested the possibility of using galaxy samples for
which the parameter $p^{(i)}$ is redshift dependent. The results shown here
are robust to these changes. This is further confirmation of our findings and
we explicitly see that maximizing the combination of biases given in (5.1)
leads to the same relative improvement in $\sigma(f_{\mathrm{NL}}^{\Delta})$
for all $\Delta$.

The results shown in the previous figures suggest that it may be possible to
improve $\sigma(f_{\mathrm{NL}}^{\Delta})$ by more than a factor of two if we
had knowledge about the parameter $p^{(i)}$ of the tracers. As previously
discussed, however, most forecasts for current or near-future galaxy surveys
assume the scale-dependent bias takes the form predicted by a universal halo
mass function. In particular, these surveys do not select targets based on the
value of $p$ (or $b_{\phi}$) which implies that there is some level of
uncertainty in the actual value of $p$ for the employed galaxy samples. Since
the parameter $p$ is degenerate with $f_{\mathrm{NL}}^{\Delta}$, this
uncertainty in $p$ could complicate the measurement of
$f_{\mathrm{NL}}^{\Delta}$. While complete ignorance of $p$ (i.e.
marginalizing over it as a free parameter in each redshift bin) would make
constraining $f_{\mathrm{NL}}^{\Delta}$ impossible, it would be helpful to
know the required level of precision so that the universality assumption does
not significantly affect the inference of $f_{\mathrm{NL}}^{\Delta}$.

We investigate this question by imposing different Gaussian priors on the
parameter $p$ for a fiducial value of $p=1$. Figure 14 shows that an
uncertainty in $p$ mostly affects $\Delta\lesssim 0.5$ when marginalization
over the biases. In addition, the degradation in the constraining power on
$f_{\mathrm{NL}}^{\Delta}$ is limited to within a factor of less than three
(two) for an extremely conservative prior with width $\sigma(p)=10\text{
(5)}$. This implies that a lack of knowledge of the precise value of $p$ for
the given galaxy samples does not significantly affect the forecasts. On the
other hand, a careful consideration of the values of $p$ in target selection
could potentially lead to important improvements in
$\sigma(f_{\mathrm{NL}}^{\Delta})$.

Figure 14: Forecast for $\sigma(f_{\mathrm{NL}}^{\Delta})$ for different
Gaussian priors $\sigma(p)$ on the parameter $p$ in the universality relation
(2.22) for the billion-object survey (at fixed $\Lambda$CDM parameters and
with fiducial value $\bar{f}_{\mathrm{NL}}^{\Delta}=1$ so that $p$ can be
independently varied). While marginalizing over $p$ with wide priors
significantly impacts the constraints for small $\Delta$ in the case of fixed
bias parameters (left), the realistic case of marginalizing over the biases
(right) shows a relatively small degradation even for very wide priors on $p$.

###  5.2 Multi-Tracer LSST and Comparison to SPHEREx

So far, we have seen that the multi-tracer forecasts offer somewhat limited
potential for improvements in the context of both current and future,
realistic spectroscopic surveys. Acquiring spectra takes more observational
time than imaging alone which makes it harder for spectroscopic surveys to
reach the high number densities needed for sample-variance cancellation to be
most effective. In principle, photometric surveys, such as LSST, are therefore
not limited in the same way, as indicated in Fig. 4, and have the potential to
benefit significantly from a multi-tracer analysis. Following the LSST science
book [89], our LSST forecasts however use a single tracer as shown in Table 8
and, consequently, do not capture the full potential of this survey which we
will now investigate.

In order to explore the potential for an LSST multi-tracer analysis, we
consider a simplified case in which we split the LSST sample of Table 8 into
two subsamples of equal number density but different linear biases
$b_{1}^{(i)}$.202020For simplicity, we split the LSST sample in the same way
as the sample for the billion-object survey, i.e. we take
$b_{1}^{(1,2)}(z=0)=b_{1}^{\mathrm{single}}(z=0)\pm 0.4$ with equal number
densities, $\bar{n}_{g}^{(1)}=\bar{n}_{g}^{(2)}$, to retain the single-tracer
LSST sample when combined. We refer to [138] for similar forecasts for an LSST
sample split into “blue” and “red” galaxies. In Figure 15, we compare the
constraints $\sigma(f_{\mathrm{NL}}^{\Delta})$ from single- and double-tracer
configurations of LSST and the billion-object survey. As for all of our
spectroscopic surveys with their given number densities, there is a relatively
small change in the constraining power for the billion-object survey.212121In
this context, it is important to remember the large maximum redshift for this
survey ($z_{\mathrm{max}}=5$) which results in a number density more or less
comparable to DESI/Euclid for many redshift bins. The same number of objects
in a smaller volume, as would for example be the case for a spectroscopic
follow-up survey to LSST with $z_{\mathrm{max}}=3$, would not have the same
limitation, but its effectiveness for most other science targets would be
reduced [2]. On the other hand, LSST shows a remarkable improvement which is
consistent with the factor of roughly three seen in Fig. 4. This increase in
sensitivity however depends of course on the biases of the samples and would
be diminished to less than a factor of two if we had instead taken
$b_{1}^{(1,2)}(z=0)=b_{1}^{\mathrm{single}}(z=0)\pm 0.2$, for instance. We
additionally note that the much larger number density of the LSST samples
appears to outweigh the larger number of (linear) modes in the billion-object
survey.

Figure 15: Comparison of the constraining power of single- and double-tracer
versions of LSST and the billion-object survey (with fixed $\Lambda$CDM
parameters) for $\Delta=0$ (left) and $2$ (right). We can clearly see the
significant improvement that a double-tracer analysis of LSST could provide
due to its large number density for the entire range of scaling exponents
$\Delta$.

In this regard, LSST has much more potential for testing inflationary physics
than is typically seen in forecasts. Of course, spectroscopic information may
prove essential in eliminating systematics or reliably splitting the LSST
sample. The improvements that we see in Fig. 15 however strongly encourage a
future investigation that accounts for systematic effects. In addition, we
have not included the information from cosmic-shear observations, which are
sensitive to primordial non-Gaussianity in their own right (see e.g. [139])
and may be useful in conjunction with the galaxy power spectra considered
here. Cross-correlations with other observables, such as CMB lensing and the
kinetic Sunyaev-Zel’dovich effect, would additionally provide complementary
constraining power [140, 141, 142].

Given the potential of LSST as a multi-tracer survey, we should also revisit
how it compares to other near-term surveys. We noted in §3.2 that LSST and
SPHEREx result in similar forecasts despite having different strengths and
weaknesses in their designs. SPHEREx is a spectro-photometric survey with a
larger sky coverage, while LSST is a photometric survey with a larger number
density. The comparison in §3.2 is additionally complicated by the fact that
we treated LSST as a single-tracer survey following [89] while multi-tracer
forecasts would maximize the benefit of the higher number density.

We might have naively imagined that the slight advantage of SPHEREx was due to
the smaller photometric redshift error of most of its samples. (Specifically,
the redshift uncertainty of LSST is assumed to be $\sigma_{z0}=0.05$, while
the SPHEREx samples are binned with maximum errors of
$\sigma_{z0}=\\{0.003,0.01,0.03,0.1,0.2\\}$.) This is however not what is
found in our forecasts. In Figure 9, we see that the potential impact of
photometric redshift errors on $f_{\mathrm{NL}}^{\Delta}$ for $\Delta\lesssim
1.2$ is at most $50\%$. At the same time, we observed in Fig. 15 that the
potential improvements in LSST from splitting the single sample into two is
potentially a factor of three or more. This suggests that the number density
is a larger effect than the quality of the photometric redshifts.222222While
the observed volume of SPHEREx is in principle more than twice as large (see
Appendix A), its number densities are much smaller, especially at high
redshifts, which results in an effective volume that is smaller by a factor of
almost four compared to LSST. At the same time, there are however many other
factors that differ between these surveys and we would like to identify the
most important individual design factors that set the two surveys apart (at
least for the purpose of constraining $f_{\mathrm{NL}}^{\Delta}$ from the
scale-dependent bias in the power spectrum).

In order to isolate the differences and control the individual factors in
survey design, we compare LSST and SPHEREx by their biases and number
densities in Fig. 16,

Figure 16: Constraining power of SPHEREx with its fiducial biases
$b_{1}^{(i)}$ and number densities $\bar{n}_{g}^{(i)}$ compared to a survey
with either of these observational design factors replaced by LSST-like
numbers for fixed (left) and marginalized (right) Gaussian bias parameters.
The solid/dashed lines show the forecasts for a single-/multi-tracer version
of these surveys (five samples for SPHEREx and two samples when considering
LSST-like design factors) for which we ignore the photometric redshift errors
in order to directly infer the sizable impact of these design choices.

while ignoring the photometric redshift errors for both surveys,
$\sigma_{z0}=0$. To be precise, we take SPHEREx as the fiducial survey and
substitute either the linear biases or the number densities from a survey
equivalent to LSST (see Appendix A for details). When we keep the number
densities of SPHEREx and use the LSST-like biases, i.e. two samples with equal
number densities, which add to the total number density of SPHEREx, and the
same biases as in our earlier LSST forecasts just over the larger redshift
range of SPHEREx with the same scaling proportional to $D(z)^{-1}$, SPHEREx
still outperforms in the single-tracer case, but underperforms in the double-
tracer scenario. In this sense, SPHEREx has a built-in advantage from the
fiducial biases used in our forecasts, while our relatively large splitting in
the linear biases of the LSST-like samples leads to a considerable improvement
in the regime of sample-variance cancellation. If we instead keep the biases
of SPHEREx and use the number densities of LSST, i.e. we include objects with
LSST-like number densities according to (A.13) beyond its $z_{\mathrm{max}}$,
the latter becomes more sensitive. In fact, this advantage becomes larger when
we additionally allow for a double-tracer configuration of LSST, which takes
full advantage of the difference in number densities. As a result, we should
understand that the slightly better sensitivity of SPHEREx indicated in Fig. 3
is primarily due to the single-tracer nature of the LSST forecast and the
larger biases of the SPHEREx samples.

Before concluding, we have to however state that the largest challenge for
both LSST and SPHEREx will of course be identifying and mitigating the impact
of large-angle (and other potential) systematics. In this regard, the survey
which will be most sensitive to $f_{\mathrm{NL}}^{\Delta}$ cannot be easily
anticipated ahead of time. Either way, both surveys have the potential to in
principle place constraints on general PNG from the galaxy power spectrum
alone (i.e. without the galaxy bispectrum) that exceed those inferred from the
Planck bispectrum for $\Delta\lesssim 1$.

##  6 Conclusions and Outlook

Primordial non-Gaussianity from scale-dependent bias offers one of the best
avenues to test the physics of inflation in the coming decade and beyond. The
signal imprinted by light particles coupled to the inflaton is localized on
relatively large scales where analyses are less influenced by gravitational
nonlinearities, baryonic physics and modeling errors. Equally importantly, the
scaling behavior of this signal encodes the mass of these light particles,
offering a unique window into the spectrum and interactions of particles
during the inflationary epoch.

Pragmatically, the scale-dependent bias offers a method to test for new
physics within the power spectra measured in large-scale structure. One might
expect that the best measurements of primordial physics should ultimately come
from higher-$N$-point statistics or even from analyses performed directly at
the level of the observed maps.232323In some cases, it has already been shown
that the signals of extra fields are best understood as map-level features
[143, 144, 145, 146, 147, 148] that may carry more signal to noise than low-
order statistics. Technical obstacles have however slowed down progress and
obscured the ultimate sensitivity of LSS surveys. In fact, analyses of the
bispectrum only recently became possible for single-field inflation
(equilateral PNG) and remain far from the demonstrated sensitivity of the CMB.
The established theoretical and observational control of LSS power spectra
therefore offer a reliable approach for near-term surveys to impact our
understanding of inflation.

Previous work has largely focused on local PNG for which the signal to noise
is dominated by only the largest angular scales accessible in a given survey.
These measurements are exciting since they have the potential to exceed the
sensitivity of the CMB and explore the interesting theoretical threshold of
$f_{\mathrm{NL}}^{\mathrm{loc}}\sim O(1)$. The much richer phenomenology of
inflation however raises the question if optimizing a survey for
$f_{\mathrm{NL}}^{\mathrm{loc}}$ is the same as optimizing for understanding
inflation more broadly. For example, degeneracies [149, 150, 151] and large-
angle (systematic) effects [103, 152, 153, 154, 155, 156] that limit the local
PNG measurement may not impact other changes to the statistics. Similarly,
some variations in the observational strategy or target selection may not
significantly affect the constraining power on
$f_{\mathrm{NL}}^{\mathrm{loc}}$, but could yield a vast increase in the
insights that we can gain on other inflationary observables.

In order to extend the reach of current and future surveys to these compelling
inflationary targets, we studied primordial non-Gaussianity associated with
light fields during inflation over a wide range of their masses. We found that
surveys optimized for $f_{\mathrm{NL}}^{\mathrm{loc}}$ tend to perform equally
well over a broader mass range beyond the massless (local) limit. Furthermore,
we placed new constraints on these types of PNG by analyzing the BOSS DR12
dataset in two different and consistent ways. These results are consistent
with previous analyses in the case of massless fields. In addition, they
extend to heavier fields and other types of non-Gaussianity. They also show
that bounds on these general non-Gaussian shapes can be simply inferred from
scale-dependent bias measurements of $f_{\mathrm{NL}}^{\mathrm{loc}}$ and
$f_{\mathrm{NL}}^{\mathrm{eq}}$ if they are consistent with zero. Our measured
constraints from BOSS data are also consistent with our BOSS forecasts which
suggests that our understanding of the constraining power of these surveys,
which we described in detail, is reflected in real data. These insights are
likely relevant to the analysis strategies for current or near-term
observations and the design of future surveys. For instance, we demonstrated
the significant improvements in sensitivity that may be available from a
multi-tracer analysis of LSST data.

More broadly speaking, additional fields present during inflation are
interesting in their own right, far beyond their potential impact on the
galaxy power spectrum. While we focused on models leading to power-law
corrections to the power spectrum, oscillatory contributions to inflationary
spectra can arise from (more) massive fields [76, 43, 157] or chemical
potentials [158, 159], for instance. This vast range of phenomenology
therefore indicates that there are many opportunities in the analysis of
galaxy surveys [160, 138] that remain under-explored and, as we also show in
this paper, minor additional work could result in much broader and deeper
insights into the primordial universe.

##### Acknowledgments

We are grateful to Alexandre Barreira, Daniel Baumann, Matteo Biagetti, Kyle
Dawson, Olivier Doré, Simone Ferraro, Marilena LoVerde, Azadeh Moradinezhad
Dizgah, and Charuhas Shiveshwarkar for helpful discussions. The authors were
supported by the US Department of Energy under Grant DE-SC0009919. B. W. also
acknowledges support from the Swedish Research Council (Contract No.
638-2013-8993). This work is based on observations obtained by the Sloan
Digital Sky Survey III (SDSS-III, http://www.sdss3.org/). Funding for SDSS-
III has been provided by the Alfred P. Sloan Foundation, the Participating
Institutions, the National Science Foundation and the US Department of Energy
Office of Science. We acknowledge the use of CAMB [57], CLASS [110], CLASS-PT
[109], FAST-PT [161], IPython [162] and MontePython [118, 119], and the Python
packages Matplotlib [163], NumPy [164] and SciPy [165].

##  Appendix A Forecasting Details

In this appendix, we collect additional information on our LSS Fisher
forecasts. We first provide the full set of definitions underlying our model
for the galaxy power spectrum and then describe the experimental
specifications of all galaxy surveys employed in the main text.

###  A.1 Galaxy Power Spectrum Model

We model the galaxy power spectrum at one-loop order in standard Eulerian
perturbation theory and in the bias expansion where the galaxy overdensity
$\delta_{g}$ is a function of terms up to third order in the linear matter
density contrast $\delta_{m}$ [cf. (2.6)]. The theoretical galaxy power
spectrum for two biased tracers $A$ and $B$ is then given by (2.9). The
various loop contributions in that equation are defined as follows:

$\displaystyle\sigma^{4}$
$\displaystyle=\int_{\vec{q}}P_{\mathrm{lin}}^{2}(q)\,,$ (A.1) $\displaystyle
P_{22}(k)$
$\displaystyle=2\int_{\vec{q}}\left[F^{(s)}_{2}(\vec{q},\vec{k}-\vec{q})\right]^{2}P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(|\vec{k}-\vec{q}|)\,,$
(A.2) $\displaystyle P_{13}(k)$
$\displaystyle=6P_{\mathrm{lin}}(k)\int_{\vec{q}}F^{(s)}_{3}(\vec{k},\vec{q},-\vec{q})P_{\mathrm{lin}}(q)\,,$
(A.3) $\displaystyle P_{\delta^{2}}(k)$
$\displaystyle=2\int_{\vec{q}}F^{(s)}_{2}(\vec{q},\vec{k}-\vec{q})P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(|\vec{k}-\vec{q}|)\,,$
(A.4) $\displaystyle P_{s^{2}}(k)$
$\displaystyle=2\int_{\vec{q}}F^{(s)}_{2}(\vec{q},\vec{k}-\vec{q})\left(\mu_{-}^{2}-\tfrac{1}{3}\right)P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(|\vec{k}-\vec{q}|)\,,$
(A.5) $\displaystyle P_{\delta^{2}\delta^{2}}(k)$
$\displaystyle=2\int_{\vec{q}}P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(|\vec{k}-\vec{q}|)\,,$
(A.6) $\displaystyle P_{\delta^{2}s^{2}}(k)$
$\displaystyle=2\int_{\vec{q}}\left(\mu_{-}^{2}-\tfrac{1}{3}\right)P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(|\vec{k}-\vec{q}|)\,,$
(A.7) $\displaystyle P_{s^{2}s^{2}}(k)$
$\displaystyle=2\int_{\vec{q}}\left(\mu_{-}^{2}-\tfrac{1}{3}\right)^{\\!2}P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(|\vec{k}-\vec{q}|)\,,$
(A.8) $\displaystyle P_{\Pi\Pi^{[2]}}(k)$
$\displaystyle=2P_{\mathrm{lin}}(k)\int_{\vec{q}}\left\\{\frac{2}{7}\left[\\!\left(\frac{\vec{k}\cdot\vec{q}}{k\hskip
1.0ptq}\right)^{\\!2}-1\right]\frac{2}{3}\mu_{-}^{2}+\frac{8}{63}\right\\}P_{\textrm{lin}}(q)\,,$
(A.9)

with $\int_{\vec{q}}\equiv\int\hskip-3.0pt\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}$
and $\mu_{-}\equiv\frac{\vec{q}\cdot(\vec{k}-\vec{q})}{q|\vec{k}-\vec{q}|}$.
The functions $F_{2}^{(s)}(\vec{k}_{1},\vec{k}_{2})$ and
$F_{3}^{(s)}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3})$ are the symmetric, second-
and third-order kernels of $\delta_{m}$ in standard perturbation theory [51],

$\displaystyle F_{2}^{(s)}(\vec{k}_{1},\vec{k}_{2})=$ $\displaystyle\
\frac{5}{7}+\frac{2}{7}\left(\frac{\vec{k}_{1}\cdot\vec{k}_{2}}{k_{1}k_{2}}\right)^{\\!2}+\frac{1}{2}\frac{\vec{k}_{1}\cdot\vec{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right)\,,$
(A.10) $\displaystyle F_{3}^{(s)}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3})=$
$\displaystyle\
\frac{k^{2}}{27}\left[\frac{\vec{k}_{1}\cdot(\vec{k}_{2}+\vec{k}_{3})}{k_{1}^{2}(\vec{k}_{2}+\vec{k}_{3})^{2}}G_{2}^{(s)}(\vec{k}_{2},\vec{k}_{3})+2\leavevmode\nobreak\
\text{cyclic}\right]$
$\displaystyle+\frac{7}{54}\vec{k}\cdot\left[\frac{\vec{k}_{1}+\vec{k}_{2}}{(\vec{k}_{1}+\vec{k}_{2})^{2}}G_{2}^{(s)}(\vec{k}_{1},\vec{k}_{2})+2\leavevmode\nobreak\
\text{cyclic}\right]$ (A.11)
$\displaystyle+\frac{7}{54}\vec{k}\cdot\left[\frac{\vec{k}_{1}}{k_{1}^{2}}F_{2}^{(s)}(\vec{k}_{2},\vec{k}_{3})+\text{2
cyclic}\right],$

where $\vec{k}\equiv\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3}$ and
$G_{2}^{(s)}(\vec{k}_{1},\vec{k}_{2})$ is the symmetric, second-order kernel
of the velocity divergence $\theta_{m}$,

$G_{2}^{(2)}(\vec{k}_{1},\vec{k}_{2})=\frac{3}{7}+\frac{4}{7}\left(\frac{\vec{k}_{1}\cdot\vec{k}_{2}}{k_{1}k_{2}}\right)^{\\!2}+\frac{1}{2}\frac{\vec{k}_{1}\cdot\vec{k}_{2}}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}}+\frac{k_{2}}{k_{1}}\right).$
(A.12)

###  A.2 Survey Specifications

In the following, we provide detailed information about the spectroscopic and
photometric redshift surveys that we employ in our Fisher forecasts, including
the assumed redshift distribution in bias and number density.

#### Spectroscopic Surveys

The experimental specifications of the spectroscopic redshift surveys BOSS,
DESI, Euclid, SPHEREx, MegaMapper and the billion-object survey are provided
in Tables 3 to 7.

$\bar{z}$ | $0.35$ | $0.625$
---|---|---
$b_{1}$ | $1.634$ | $1.877$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $0.275$ | $0.142$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $2.18$ | $4.15$

Table 2: Basic specifications for a BOSS-like survey [14] (inspired by [166]
as detailed in [93]), covering a sky area of
$10\,252\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$ with a total of about
$1.2\text{\times}{10}^{6}$ objects in a volume of roughly
$6.3\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$.

$\bar{z}$ | $0.05$ | $0.15$ | $0.25$ | $0.35$ | $0.45$ | $0.65$ | $0.75$ | $0.85$ | $0.95$
---|---|---|---|---|---|---|---|---|---
$b_{1}$ | $1.40$ | $1.48$ | $1.55$ | $1.61$ | $1.67$ | $2.05$ | $1.71$ | $1.71$ | $1.53$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $38.8$ | $15.7$ | $3.96$ | $0.883$ | $0.0992$ | $0.591$ | $1.31$ | $0.920$ | $0.779$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $0.0356$ | $0.229$ | $0.560$ | $0.979$ | $1.44$ | $2.39$ | $2.83$ | $3.24$ | $3.61$

$\bar{z}$ | $1.05$ | $1.15$ | $1.25$ | $1.35$ | $1.45$ | $1.55$ | $1.65$ | $1.75$ | $1.85$
---|---|---|---|---|---|---|---|---|---
$b_{1}$ | $1.45$ | $1.48$ | $1.47$ | $1.47$ | $1.69$ | $1.68$ | $2.27$ | $2.45$ | $2.47$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $0.466$ | $0.398$ | $0.387$ | $0.180$ | $0.133$ | $0.110$ | $0.0387$ | $0.0197$ | $0.0208$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $3.94$ | $4.24$ | $4.49$ | $4.71$ | $4.90$ | $5.05$ | $5.18$ | $5.29$ | $5.37$

Table 3: Basic specifications for DESI [14] (derived from [86] as explained in
[93]), covering a sky area of
$14\,000\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$ with a total of about
$2.7\text{\times}{10}^{7}$ objects in a volume of roughly
$58\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$.

We provide the linear bias $b_{1}$ and number density $\bar{n}_{g}$ for each
sample and redshift bin with mean redshift $\bar{z}$ and spherical volume $V$.
For DESI, different types of tracers were combined into one galaxy sample with
a single effective number density and bias following [93, 14]. Since the
number density of DESI is not in the regime where sample-variance cancellation
is most effective (see Fig. 4), our results should be relatively insensitive
to this choice. For all surveys, we treat each redshift bin as being
independent so that the Fisher matrix for the entire survey is the sum of the
Fisher matrices associated with each separate redshift bin.

$\bar{z}$ | $0.65$ | $0.75$ | $0.85$ | $0.95$ | $1.05$ | $1.15$ | $1.25$ | $1.35$
---|---|---|---|---|---|---|---|---
$b_{1}$ | $1.06$ | $1.11$ | $1.16$ | $1.21$ | $1.27$ | $1.33$ | $1.38$ | $1.44$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $0.637$ | $1.46$ | $1.63$ | $1.50$ | $1.33$ | $1.14$ | $1.00$ | $0.837$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $2.56$ | $3.03$ | $3.47$ | $3.87$ | $4.23$ | $4.54$ | $4.81$ | $5.05$

$\bar{z}$ | $1.45$ | $1.55$ | $1.65$ | $1.75$ | $1.85$ | $1.95$ | $2.05$
---|---|---|---|---|---|---|---
$b_{1}$ | $1.51$ | $1.54$ | $1.63$ | $1.70$ | $1.85$ | $1.90$ | $1.26$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $0.652$ | $0.512$ | $0.357$ | $0.246$ | $0.149$ | $0.0904$ | $0.0721$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $5.25$ | $5.41$ | $5.55$ | $5.67$ | $5.76$ | $5.83$ | $5.88$

Table 4: Basic specifications for the Euclid survey [93] (derived from [167]), covering a sky area of $15\,000\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$ with a total of about $4.9\text{\times}{10}^{7}$ objects in a volume of roughly $71\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$. $\bar{z}$ | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 1.3 | 1.9 | 2.5 | 3.1 | 3.7 | 4.3
---|---|---|---|---|---|---|---|---|---|---|---
$V$ | 0.584 | 3.40 | 7.43 | 11.5 | 15.1 | 60.4 | 71.4 | 73.6 | 71.8 | 68.3 | 64.2
$b_{1}^{(1)}$ | 1.3 | 1.5 | 1.8 | 2.3 | 2.1 | 2.7 | 3.6 | 2.3 | 3.2 | 2.7 | 3.8
$${10}^{5}$\bar{n}_{g}^{(1)}$ | 997 | 411 | 50.1 | 7.05 | 3.16 | 1.64 | 0.359 | 0.0807 | 0.184 | 0.150 | 0.113
$b_{1}^{(2)}$ | 1.2 | 1.4 | 1.6 | 1.9 | 2.3 | 2.6 | 3.4 | 4.2 | 4.3 | 3.7 | 4.6
$${10}^{5}$\bar{n}_{g}^{(2)}$ | 1230 | 856 | 282 | 93.7 | 43.0 | 5.00 | 0.803 | 0.383 | 0.328 | 0.107 | 0.0679
$b_{1}^{(3)}$ | 1.0 | 1.3 | 1.5 | 1.7 | 1.9 | 2.6 | 3.0 | 3.2 | 3.5 | 4.1 | 5.0
$${10}^{5}$\bar{n}_{g}^{(3)}$ | 1340 | 857 | 362 | 294 | 204 | 21.2 | 0.697 | 0.202 | 0.143 | 0.193 | 0.0679
$b_{1}^{(4)}$ | 0.98 | 1.3 | 1.4 | 1.5 | 1.7 | 2.2 | 3.6 | 3.7 | 2.7 | 2.9 | 5.0
$${10}^{5}$\bar{n}_{g}^{(4)}$ | 2290 | 1290 | 535 | 495 | 415 | 79.6 | 7.75 | 0.787 | 0.246 | 0.193 | 0.136
$b_{1}^{(5)}$ | 0.83 | 1.2 | 1.3 | 1.4 | 1.6 | 2.1 | 3.2 | 4.2 | 4.1 | 4.5 | 5.0
$${10}^{5}$\bar{n}_{g}^{(5)}$ | 1490 | 752 | 327 | 250 | 183 | 73.4 | 25.3 | 5.41 | 2.99 | 0.941 | 0.204
Table 5: Basic specifications for SPHEREx [168], covering a sky fraction ${f_{\mathrm{sky}}}=0.75$ with a total of about $7.0\text{\times}{10}^{8}$ objects in a volume of roughly $450\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$. SPHEREx is a spectro-photometric survey and the observed objects are divided into five samples based on their photometric redshift uncertainty bin, with respective maximum error $\sigma_{z0}=\\{0.003,0.01,0.03,0.1,0.2\\}$. The volume $V$ and the galaxy number density $\bar{n}_{g}$ are given in units of ${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$ and ${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$, respectively. $\bar{z}$ | $2.0$ | $3.0$ | $4.0$ | $5.0$
---|---|---|---|---
$b_{1}$ | $2.5$ | $4.0$ | $3.5$ | $5.5$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $0.98$ | $0.12$ | $0.10$ | $0.040$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $54.0$ | $54.4$ | $50.0$ | $44.8$
Table 6: Basic specifications for MegaMapper [90], covering a sky area of $14\,000\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$ with a total of about $6.6\text{\times}{10}^{7}$ objects in a volume of roughly $200\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$. $\bar{z}$ | $0.25$ | $0.75$ | $1.25$ | $1.75$ | $2.25$ | $2.75$ | $3.25$ | $3.75$ | $4.25$ | $4.75$
---|---|---|---|---|---|---|---|---|---|---
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $4.78$ | $20.6$ | $32.8$ | $38.8$ | $40.8$ | $40.7$ | $39.4$ | $37.7$ | $35.9$ | $34.0$
$b_{1}^{(1)}$ | $2.28$ | $2.94$ | $3.66$ | $4.41$ | $5.17$ | $5.94$ | $6.72$ | $7.51$ | $8.29$ | $9.07$
$${10}^{3}$\bar{n}_{g}^{(1)}\,[${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$
$b_{1}^{(2)}$ | $1.37$ | $1.76$ | $2.19$ | $2.64$ | $3.10$ | $3.57$ | $4.03$ | $4.50$ | $4.97$ | $5.45$
$${10}^{3}$\bar{n}_{g}^{(2)}\,[${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$ | $1.54$
Table 7: Basic specifications for our (futuristic) spectroscopic billion-
object survey, covering a sky fraction ${f_{\mathrm{sky}}}=0.5$ with a total
of ${10}^{9}$ objects in a volume of roughly
$330\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$. The
observed objects are divided into two samples based on their linear bias,
$b_{1}(z=0)=2.0\text{ and }1.2$.

#### Photometric Surveys

To forecast Vera Rubin Observatory’s LSST, which is a photometric redshift
survey, we use the “gold” sample defined in their science book [89], which
includes more than four billion galaxies over
$20\,000\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$. The number of objects in
each redshift bin is calculated according to

$N_{g}(z_{\mathrm{min}},z_{\mathrm{max}})=\left.\frac{N_{\mathrm{tot}}\Omega_{\mathrm{sky}}}{2z_{0}^{2}}\left(2z_{0}^{2}+2z_{0}z+z^{2}\right)\mathrm{e}^{-z/z_{0}}\right|^{z=z_{\mathrm{max}}}_{z=z_{\mathrm{min}}},$
(A.13)

where
$N_{\mathrm{tot}}=$0.55\text{\,}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{m}\mathrm{i}\mathrm{n}^{-2}$$,
$\Omega_{\mathrm{sky}}=$20\,000\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$$
and $z_{0}=0.3$. The bias is taken to be $b_{1}(z)=0.95/D(z)$ and we assume a
photometric redshift error of $\sigma_{z0}=0.05$. Table 8

$\bar{z}$ | $0.125$ | $0.375$ | $0.625$ | $0.875$ | $1.125$ | $1.375$
---|---|---|---|---|---|---
$b_{1}$ | $1.01$ | $1.16$ | $1.31$ | $1.48$ | $1.65$ | $1.82$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $292$ | $183$ | $108$ | $63.6$ | $36.8$ | $21.0$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $0.709$ | $3.93$ | $8.10$ | $11.9$ | $14.9$ | $17.0$

$\bar{z}$ | $1.625$ | $1.875$ | $2.125$ | $2.375$ | $2.625$ | $2.875$
---|---|---|---|---|---|---
$b_{1}$ | $2.00$ | $2.18$ | $2.36$ | $2.54$ | $2.73$ | $2.91$
$${10}^{3}$\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $11.9$ | $6.46$ | $3.55$ | $1.94$ | $0.969$ | $0.646$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $18.4$ | $19.2$ | $19.7$ | $19.9$ | $19.8$ | $19.6$

Table 8: Basic specifications for Vera Rubin Observatory’s LSST [89], covering
a sky area of $20\,000\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$ with a
total of about $4\text{\times}{10}^{9}$ objects in a volume of roughly
$170\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$. The
photometric redshift error is $\sigma_{z0}=0.05$.

explicitly provides the derived experimental specifications for this survey
with $z_{\mathrm{max}}=3$ and twelve redshift bins with width $\Delta z=0.25$.

##  Appendix B Measuring the Scaling Behavior

We have assumed for all forecasts throughout the paper that the fiducial value
of the non-Gaussian amplitude $f_{\mathrm{NL}}^{\Delta}$ vanishes. This is of
course motivated by the fact that we have only been placing observational
bounds on PNG that are consistent with $f_{\mathrm{NL}}^{\Delta}=0$, including
in Section 4. We additionally always considered the scaling exponent $\Delta$
to be fixed. In this appendix, we will consider nonzero fiducial PNG
amplitudes and clarify to what degree measurements of the galaxy power
spectrum are sensitive to the precise value of $\Delta$.

We defined the correlation between the non-Gaussian signals of two different
values of $\Delta$ in (4.2) in terms of the respective Fisher matrix elements
(after marginalization) for a fiducial value of