# Constraining general multi-field inflation using the SPHEREx all-sky survey

Charuhas Shiveshwarkar  Thejs Brinckmann  Marilena Loverde

###### Abstract

We investigate how well the SPHEREx all-sky survey can constrain local
primordial non-Gaussianity beyond the parameter $f_{\mathrm{NL}}$ using the
galaxy power spectrum. We forecast joint constraints on the parameters
$f_{\mathrm{NL}}$, $g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ obtained
assuming a simple two-field curvaton model of inflation. The parameters
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ characterise the squeezed limits of
the primordial bispectrum and trispectrum respectively, and lead to a
characteristic scale-dependence of the galaxy bias that increases out to
arbitrarily large scales. Values of the parameter
$\tau_{\mathrm{NL}}>(\frac{6}{5}f_{\mathrm{NL}})^{2}$ cause the galaxy power
spectrum to have a stochastic component which also increases out to
arbitrarily large scales. Our MCMC forecasts indicate that SPHEREx can provide
joint constraints on any two of the three parameters
$f_{\mathrm{NL}},g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$. Due to strong
degeneracies among these parameters, measurements of the galaxy power spectra
alone may not be sufficient to jointly constrain all three. Constraints on
$f_{\mathrm{NL}},\ g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ obtained from
galaxy power spectrum observations depend on the modelling of underlying
nuisance parameters. We study the robustness of our forecast constraints to
modelling choices and note that even with relatively weak assumptions, SPHEREx
galaxy power spectra can provide strong evidence of local non-Gaussianity,
even if the particular values of $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$
cannot be measured precisely.

##  1 Introduction

A major goal of cosmology is to determine the source of the initial seeds of
structure in our Universe [1, 2]. In the simplest scenario, the primordial
curvature perturbations were sourced by quantum fluctuations in a single
inflaton field, whose energy density is also responsible for driving the
accelerated expansion of the Universe during inflation [3, 4, 5]. In this
case, the primordial perturbations are nearly Gaussian, with departures from
Gaussianity only arising if the inflaton has strong self-interactions [6, 7,
8] or unusual scenarios in which the inflaton field is not on the slow-roll
attractor solution to the equations of motion [9, 10]. If additional light
($m\ll H$, where $H$ is the Hubble rate during inflation) fields are present
during inflation, then these fields will acquire a spectrum of primordial
fluctuations in addition to the fluctuations from the inflaton. Such
fluctuations can be converted into fluctuations in curvature through processes
that are generically nonlinear but local in configuration space. Examples of
such scenarios are the curvaton model and modulated reheating [11, 12, 13, 14,
15]. Curvature perturbations generated in this manner have distinct non-
Gaussian statistical distributions that can be used to infer their presence.

The non-Gaussian primordial curvature perturbation, $\zeta$, arising from
fluctuations in additional light fields can typically be written as,

$\zeta(\vec{x})=\mathcal{Z}(\vec{x})+\frac{3}{5}f_{\mathrm{NL}}(\mathcal{Z}(\vec{x})^{2}-\langle\mathcal{Z}^{2}\rangle)+\frac{9}{25}g_{\mathrm{NL}}(\mathcal{Z}(\vec{x})^{3}-3\langle\mathcal{Z}^{2}\rangle\mathcal{Z}(\vec{x}))+\dots$
(1.1)

where $\mathcal{Z}$ is a Gaussian random field and the coefficients
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ characterise nonlinear terms which are
subdominant in the limit
$f_{\mathrm{NL}}\langle\mathcal{Z}^{2}\rangle^{1/2}\,,g_{\mathrm{NL}}\langle\mathcal{Z}^{2}\rangle\ll
1$.111This definition of $g_{\mathrm{NL}}$ (with the
$-3\langle\mathcal{Z}^{2}\rangle\mathcal{Z}(\vec{x})$ subtracted) is chosen so
that the variance of $\mathcal{Z}$ equals the variance of $\zeta$ at lowest
order in $g_{\mathrm{NL}}$. The expansion in Eq. (1.1) is said to describe
Local Primordial non-Gaussianity (LPnG) because it is generated through a
mapping of a Gaussian field that is local in configuration space. The
primordial bispectrum of $\zeta$ from Eq. (1.1) is proportional to
$f_{\mathrm{NL}}$, while the primordial trispectrum contains two terms, one
proportional to $g_{\mathrm{NL}}$ and another proportional to
$f_{\mathrm{NL}}^{2}$. A characteristic of LPnG initial conditions is that the
higher-order correlation functions of $\zeta$ peak in the so-called squeezed-
limit, when one of the wavenumbers of the $\zeta$ fields is taken to zero.
This is a feature that is distinct from primordial curvature generated by
fluctuations from a single inflaton field [16]. The detection of LPnG would
therefore be of tremendous importance, and imply that the initial conditions
of our Universe were sourced by a new light field, over and above the inflaton
[17]. The best constraints on the parameters $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ currently come from the Planck cosmic microwave background
(CMB) datasets [18], according to which $f_{\mathrm{NL}}=-0.9\pm 5.1$ and
$g_{\mathrm{NL}}=(-5.8\pm 6.5)\times 10^{4}$ at $68\%$ confidence. Typical
multi-field inflationary models (including the aforementioned curvaton models)
predict $f_{\mathrm{NL}}\approx 1$ [1]. A robust detection of
$f_{\mathrm{NL}}\approx 1$ is therefore an important theoretical target and
would conclusively establish the existence of one or more additional light
fields during the epoch of inflation.

A generalisation of the non-Gaussian initial conditions in Eq. (1.1) is to
allow for an observable contribution to the primordial curvature $\zeta$ that
is Gaussian and sourced by the inflaton, in addition to that sourced the
curvaton (see, e.g, Refs. [19, 20]).222From this point forward, for simplicity
we refer to non-Gaussian initial curvature perturbations as in Eq. (1.1) as
sourced by the curvaton, though other inflationary models with additional
light fields may generate similar initial conditions. In this case, the
initial curvature perturbations can be written as a sum of terms coming from
inflaton, $\zeta_{I}$ and curvaton, $\zeta_{C}$,

$\displaystyle\zeta(\vec{x})=\zeta_{I}(\vec{x})+\zeta_{C}(\vec{x})\ .$ (1.2)

Above, $\zeta_{I}$ is a Gaussian random field and $\zeta_{C}$ is given by an
expansion such as Eq. (1.1). As we shall see, these scenarios have suppressed
non-Gaussianity, relative to a model where $\zeta_{I}$ fluctuations are
negligible, yet have an $\mathcal{O}(f_{\mathrm{NL}}^{2})$ trispectrum that is
enhanced relative to the bispectrum in single-source initial conditions from
Eq. (1.1). This trispectrum peaks in the collapsed limit and its amplitude is
characterised by a parameter $\tau_{\mathrm{NL}}$ that is bounded from below
as $\tau_{\mathrm{NL}}\geq(\frac{6}{5}f_{\mathrm{NL}})^{2}$ [21, 22]. Initial
conditions that exceed this bound can be thought of as having a collapsed
trispectrum that is not fully correlated with the square of the primordial
bispectrum. Current bounds on $\tau_{\mathrm{NL}}$ from Planck CMB data are
$\tau_{\mathrm{NL}}\leq 2800$ at $95\%$ confidence and therefore far above the
minimum threshold set by $f_{\mathrm{NL}}$ constraints.

Although at present the best constraints on primordial non-Gaussianity come
from CMB datasets, major advances are expected to come from large-scale
structure [1]. For scenarios with LPnG, the squeezed-limit $N+1$-point
primordial correlation functions modulate $N$-point functions by the large-
scale curvature field $\zeta(\vec{x})$. This in turn manifests as a modulation
of the galaxy distribution by $\zeta(\vec{x})$, or a scale-dependent galaxy
bias with respect to the density field [23]. The scale-dependent galaxy bias
is a key signal (in addition to the galaxy bispectrum) used by large scale
structure surveys to constrain $f_{\mathrm{NL}}$ [24, 25, 26, 27, 28, 1, 29,
30, 31]. Current constraints on $f_{\mathrm{NL}}$ from the BOSS galaxy survey
are $f_{\mathrm{NL}}=9^{+33}_{-35}$ from the galaxy power spectrum alone and
$f_{\mathrm{NL}}=-33\pm 28$ with the addition of the galaxy bispectrum, both
at $68\%$ confidence [30]. Yet, the scale-dependent bias is the combined
effect of the squeezed limits of arbitrary higher point primordial correlation
functions. In other words, all the non-linear coefficients $f_{\mathrm{NL}},\
g_{\mathrm{NL}}\,\dots$ in Eq. (1.1) each contribute to the scale-dependent
galaxy bias (see, e.g, Refs. [32, 33]). As these non-linear coefficients all
contribute to the galaxy power spectrum, there is a strong degeneracy among
them when constrained using the galaxy power spectrum alone. On the other
hand, as we will discuss, scale-dependent bias constraints can be interpreted
as constraints on all higher-order terms appearing in Eq. (1.1). The galaxy
power spectrum therefore provides broad constraining power on LPnG. The
tightest constraints on $g_{\mathrm{NL}}$ from scale-dependent bias come from
quasar power spectra, $-4<g_{\mathrm{NL}}\times 10^{-5}<4.7$ at $95\%$
confidence. Joint constraints on $f_{\mathrm{NL}}$, $g_{\mathrm{NL}}$ have
been performed in [26, 27]. Relatedly, for initial conditions of the form Eq.
(1.2), the halo bias also acquires a stochastic component that is scale-
dependent on horizon scales, and has been used to constrain the relative
contributions of inflaton and curvaton fluctuations [20]. Generalising future
galaxy power spectrum constraints on primordial non-Gaussianity to include
these terms beyond $f_{\mathrm{NL}}$ is important, as information about these
higher order statistics may be helpful in distinguishing between different
multi-field inflation models, thus helping to further constrain inflationary
physics.

The SPHEREx All-sky galaxy survey aims to achieve constraints on the parameter
$f_{\mathrm{NL}}$ at a level $\sigma_{f_{\mathrm{NL}}}\lesssim 0.5$ [28].
Reaching this threshold would represent an order-of-magnitude improvement over
CMB constraints and the first instance of reaching values of $f_{\mathrm{NL}}$
predicted by the simplest two-field inflation scenarios [11, 12]. SPHEREx will
achieve this target through joint analysis of galaxy auto- and cross-power
spectra in combination with galaxy bispectra. While the signatures of
$f_{\mathrm{NL}}$, $g_{\mathrm{NL}}$, and $\tau_{\mathrm{NL}}$ on the galaxy
bispectrum and trispectrum are distinct, owing to the different shapes of
their primordial correlation functions, these parameters have nearly common
signatures in the galaxy power spectrum. As we will review, both
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ contribute to scale-dependent galaxy
bias, while $\tau_{\mathrm{NL}}\geq(\frac{6}{5}f_{\mathrm{NL}})^{2}$ generates
a stochasticity between galaxies and matter that is also probed by galaxy
auto- and cross-power spectra. Therefore, the galaxy power spectra measured by
SPHEREx can be used to put broad constraints on primordial non-Gaussianity
beyond just those set by interpreting data in terms of the $f_{\mathrm{NL}}$
parameter. Our goal in this paper is to forecast constraints achievable by
SPHEREx from galaxy power spectra on LPnG where the primordial trispectrum,
characterised by $g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$, is significant.
In this, we will pay special attention to degeneracies among these parameters
and possible limitations to interpretation due to uncertainties in modelling
their signatures on the galaxy bias. While additional information on these
parameters will be available from galaxy bispectra and trispectra, in
particular information that could break degeneracies among them, the analysis
is considerably more complicated and model-dependent. Galaxy power spectra,
therefore, could play an important role in establishing the existence of LPnG,
that could subsequently be studied in more detail through higher-order galaxy
correlation functions.

The rest of this paper is organised as follows. In Section 2, we review the
derivation of scale-dependent bias from LPnG, allowing for LPnG that is
sourced by arbitrarily high orders in an expansion such as Eq. (1.1). In the
same section, we review the derivation of scale-dependent stochastic bias
sourced by multi-field initial conditions such as Eq. (1.2). Here, we also
allow the curvaton field to have polynomial terms of arbitrary high order,
before specialising to the case where terms after
$g_{\mathrm{NL}}\mathcal{Z}^{3}$ can be neglected. In Section 3, we present
our likelihood and forecasting pipeline, as well as models we use for the
SPHEREx galaxy power spectra. In Section 4, we present our forecast
constraints on $f_{\mathrm{NL}}$, $g_{\mathrm{NL}}$, and $\tau_{\mathrm{NL}}$,
along with a study of the robustness of these constraints to bias modelling
choices. Our conclusions are presented in Section 5.

##  2 Galaxy bias from local primordial non-Gaussianity

In this section, we review the derivation of the effect of LPnG on the galaxy
power spectrum measured by galaxy surveys [23], with particular focus on the
impact of higher-order primordial correlation functions, beyond the
bispectrum. Throughout our analysis, we assume a weak primordial non-
Gaussianity, such that the dimensionless cumulants
$\langle\zeta^{N}\rangle_{c}/\langle\zeta^{2}\rangle^{N/2}\ll 1$, for $N>2$
and where c refers to the connected part of the correlation function.333In
what follows we use $\langle...\rangle$ and $\langle...\rangle_{c}$
interchangeably to mean connected moments. We will, however, allow
dimensionless cumulants of different orders to be of comparable magnitude,
i.e. we will consider examples where
$\langle\zeta^{M}\rangle_{c}/\langle\zeta^{2}\rangle^{M/2}\sim\langle\zeta^{N}\rangle_{c}/\langle\zeta^{2}\rangle^{N/2}$,
for $M$, $N>2$ but $M\neq N$. For this to happen with the initial conditions
in Eq. (1.1), higher-order coefficients, such as $g_{\mathrm{NL}}$, must be
much larger in magnitude than lower-order coefficients, such as
$f_{\mathrm{NL}}$. For the bispectrum and trispectum specifically, one would
require that $g_{\mathrm{NL}}\sim f_{\mathrm{NL}}A_{S}^{-1/2}$ for them to
have comparable magnitude, where $A_{s}$ characterises the scalar power
spectrum of $\zeta$. Initial conditions such as this can be thought of as
having a primordial kurtosis that is small, but as large as the primordial
skewness, in comparison to the variance. Non-Gaussianity of this sort is
completely consistent with current datasets (and in fact explains the relative
magnitude of the current constraints on $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$).

In Section 2.1, we review how LPnG leads to a characteristic scale-dependence
in the galaxy bias that diverges quadratically at the largest scales. In the
following Section 2.2, we review how LPnG gives rise to a stochastic
contribution to the galaxy power spectrum in the case of a two-field curvaton
model of inflation wherein both the inflaton and the massless curvaton
contribute comparably to the primordial curvature perturbation [20, 34]. In
both Sections 2.1 and 2.2, we first present a more general result for scale-
dependent bias and stochasticity that holds for non-trivial squeezed limits of
arbitrary higher point primordial correlation functions and then provide
concrete expressions for primordial non-Gaussianity with only a primordial
bispectrum and trispectrum. We work exclusively in the comoving synchronous
gauge. This has the advantage that it is the unique gauge in which the linear
bias relation $(\delta_{g}=b_{g}\delta_{cb})$ remains valid (to linear order)
at large scales [35].444Unlike other gauges where there are additional terms
which become important at horizon scales [35].

###  2.1 Scale-dependent galaxy bias from local primordial non-Gaussianity

LPnG is characterised by the primordial curvature perturbation having a non-
trivial $N$-point correlation function ($N\geq 3$) that peaks in the squeezed
limit, i.e. when one of the external momenta becomes soft. It is well
understood that to lowest order, the squeezed-limit of a primordial $(N+1)$
correlation function encodes the response of the small-scale $N$-point
correlation function in the presence of a long wavelength mode of the field
whose momentum becomes soft [36, 37]. In particular, for the squeezed-limit of
the primordial correlation functions, one can show that for $N\geq 2$,

$\displaystyle\lim_{q\to
0}\frac{\langle\zeta(\vec{q})\zeta(\vec{k_{1}})\zeta(\vec{k_{2}})...\zeta(\vec{k_{N}})\rangle}{P_{\zeta}(q)}=\frac{\partial\langle\zeta(\vec{k_{1}})...\zeta(\vec{k_{N}})\rangle}{\partial\zeta_{L}}\
,$ (2.1)

where the right-hand side denotes the variation of the small-scale $N$ point
function w.r.t. a constant background curvature perturbation $\zeta_{L}$. The
above equation also holds for the linear Cold Dark Matter (CDM) and baryon
density contrast $(\delta^{lin}_{cb}$) evaluated at late times, i.e.

$\displaystyle\lim_{q\to
0}\frac{\langle\zeta(\vec{q})\delta^{lin}_{cb}(\vec{k_{1}})...\delta^{lin}_{cb}(\vec{k_{N}})\rangle}{P_{\zeta}(q)}=\frac{\partial\langle\delta^{lin}_{cb}(\vec{k_{1}})...\delta^{lin}_{cb}(\vec{k_{N}})\rangle}{\partial\zeta_{L}}\
,$ (2.2)

where the linear CDM+baryon density fields are evaluated at typical redshifts
probed by galaxy surveys. Under initial conditions with LPnG, the left-hand
side of the above equations approaches a non-trivial value implying that the
fluctuations of CDM+baryon densities acquire a dependence on an ambient large-
scale curvature perturbation $(\zeta_{L})$ that does not vanish even at the
largest scales. The equivalence principle forbids such a dependence from being
generated from nonlinear gravitational evolution [17].

To analyse the effect of LPnG on the halo bias, we adopt the peak-background-
split (PBS) approach [38]. In this approach, the matter density fluctuations
can be separated into a large-scale component that does not vary significantly
over scales as large as the size of a typical halo and a small scale component
which undergoes gravitational collapse to form dark matter haloes. The
comoving number density of dark matter tracers, e.g. galaxies or quasars,
$n_{g}$ is therefore a function of the connected moments of the linear matter
density field $\delta_{cb,R}$ smoothed over some small scale $R$ corresponding
to the typical Lagrangian radius of dark matter haloes.555While in what
follows, we write the tracer density $n_{g}$ as a function of the smoothed
fields $\delta_{cb,R}$, the methodology is more general and only requires a
separation of scales between the Fourier modes that produce our tracers and
long-wavelength variations in the number density of the tracers $n_{g}$.  For
Gaussian initial conditions, the only non-trivial such moment is the smoothed
variance of the linear density field
$\sigma^{2}_{R}=\langle\delta^{2}_{cb,R}\rangle$. For non-Gaussian initial
conditions, one can have non-trivial higher order moments
$S_{N}(R;z)=\langle\delta^{N}_{cb,R}\rangle$ given by:

$\displaystyle
S_{N}=\int\prod^{N}_{i=1}\frac{d^{3}\vec{k_{i}}}{(2\pi)^{3}}W(k_{i}R)T_{cb}(k_{i},z)\langle\zeta(\vec{k_{1}})...\zeta(\vec{k_{N}})\rangle\
,$ (2.3)

where $W(kR)$ is the window function used to calculate the smoothed matter
density field $\delta_{cb,R}$ and $T_{cb}(k,z)$ is the linear transfer
function of matter density fluctuations
$\delta_{cb}(k,z)=T_{cb}(k,z)\zeta(\vec{k})$. Due to the coupling between
small-scale and large-scale modes induced by LPnG following Eq. (2.1), the
moments of smoothed linear matter density fluctuations $S_{N}(R;z)$ acquire a
dependence on a background large scale curvature perturbation $\zeta_{L}$,
through the dependence of $S_{N}$ on the squeezed-limit of the primordial
$N+1$-point correlation function,666In comoving coordinates, the squeezed
limit can have unphysical terms which arise as a consequence of the slight
scale-dependence of primordial fluctuations [6, 17, 37]. These terms however
are a gauge artefact, arising from the choice of how physical separations are
separated into the product of a scale factor and comoving separations, and are
not observable [39, 40] and do not contribute to the scale-dependent bias
[41].

$\displaystyle\frac{\partial S_{N}}{\partial\zeta_{L}}$ $\displaystyle=$
$\displaystyle\int\prod^{N}_{i=1}\frac{d^{3}\vec{k_{i}}}{(2\pi)^{3}}W(k_{i}R)T_{cb}(k_{i},z)\frac{\partial\langle\zeta(\vec{k_{1}})...\zeta(\vec{k_{N}})\rangle}{\partial\zeta_{L}}$
(2.4) $\displaystyle=$
$\displaystyle\int\prod^{N}_{i=1}\frac{d^{3}\vec{k_{i}}}{(2\pi)^{3}}W(k_{i}R)T_{cb}(k_{i},z)\lim_{q\to
0}\frac{\langle\zeta(\vec{q})\zeta(\vec{k_{1}})\zeta(\vec{k_{2}})...\zeta(\vec{k_{N}})\rangle}{P_{\zeta}(q)}\
.$ (2.5)

The quantity $n_{g}$ will generically depend on the small-scale statistics of
the matter field such as its variance and the higher-order cumulants $S_{N}$.
Within a large-scale matter density perturbation $\delta_{L}$,
$n_{g}[\delta_{L};\\{S_{N}(\zeta_{L})\\}]$ acquires an additional dependence
on the large scale curvature perturbation $\zeta_{L}$. At linear order in
$\zeta_{L}$, we have,

$n_{g}[\delta_{L};\\{S_{i}(\zeta_{L})\\}]=n_{g}[0;S^{o}_{i}]+\frac{\partial
n_{g}}{\partial\delta_{L}}[0;S^{o}_{i}]\delta_{L}+\left(\sum^{\infty}_{N=2}\frac{\partial
n_{g}}{\partial S_{N}}\cdot\frac{\partial
S_{N}}{\partial\zeta_{L}}\right)\zeta_{L}\ ,$ (2.6)

where $S^{o}_{i}$ are the moments of the smoothed linear density field in the
absence of an ambient large-scale curvature perturbation. For typical models
with weak non-Gaussianity, only finitely many terms contribute to the
summation in Eq. (2.6). Relating $\zeta_{L}$ to $\delta_{L}$ using the
transfer function $T_{cb}(k,z)$, one can write a linear bias relation for
fluctuations in the abundance of matter density tracers

$\displaystyle\frac{\delta n_{g}}{n_{g}}=\left(\frac{\partial\log
n_{g}}{\partial\delta_{L}}+\left(\sum^{\infty}_{N=2}\frac{\partial\log
n_{g}}{\partial S_{N}}\cdot\frac{\partial
S_{N}/\partial\zeta_{L}}{T_{cb}(k,z)}\right)\right)\delta_{L}$ (2.7)

We define $b_{L}=\partial\log n_{g}/\partial\delta_{L}$ as the scale-
independent linear Lagrangian halo bias which acquires a scale-dependent
correction in the presence of LPnG,

$\displaystyle\Delta b_{NG}(k,z)=\left(\sum^{\infty}_{N=2}\frac{\partial\log
n_{g}}{\partial S_{N}}\cdot\frac{\partial
S_{N}}{\partial\zeta_{L}}\right)\cdot\frac{1}{T_{cb}(k,z)}\ .$ (2.8)

The CDM+Baryon transfer function $T_{cb}(k,z)$ can be obtained by inverting
the Poisson equation and is typically written as,

$\displaystyle T_{cb}(k,z)$ $\displaystyle=$
$\displaystyle\frac{3}{5}\alpha(k,z)\ ,$ (2.9) $\displaystyle\alpha(k,z)$
$\displaystyle=$
$\displaystyle\frac{2k^{2}\mathrm{T}(k)D(z)}{3H^{2}_{0}\Omega_{m}}\ ,$ (2.10)

where the function $T(k)$ approaches 1 as $k\rightarrow 0$777Note that
$-(3/5)T(k)$ is the transfer function for the gravitational potential deep in
the matter-dominated era. and $D(z)$ is the linear growth factor normalised to
be the scale factor during the matter-dominated era. It is convenient to write
the scale-dependent contribution to the bias (Eq. (2.8)) in the following form

$\displaystyle\Delta b_{NG}(k,z)=\frac{\mathcal{B}(z)}{\alpha(k,z)}\ ,$ (2.11)

where

$\displaystyle\mathcal{B}(z)=\frac{5}{3}\sum^{\infty}_{N=2}\frac{\partial\log
n_{g}}{\partial S_{N}}\cdot\frac{\partial S_{N}}{\partial\zeta_{L}}\ .$ (2.12)

Note that in Eq. (2.12), $\mathcal{B}(z)$ will depend on properties of the
tracer population $n_{g}$ (for instance, the mass or luminosity of galaxies),
but we have suppressed that argument here. In this way, we see that LPnG
induces a scale-dependent galaxy bias that diverges as $1/k^{2}$ at the
largest scales and that its magnitude is parameterised by a weighted sum over
the squeezed limits of primordial correlation functions (given by the factors
$\partial S_{N}/\partial\zeta_{L}$, see Eqs. (2.4) and (2.5)) weighted by the
marginal contribution of small-scale statistics of the linear density field to
the abundance of dark matter tracers (given by the factors $\partial\log
n_{g}/\partial S_{N}$). This fact makes scale-dependent bias a broad probe of
any type of local primordial non-Gaussianity. This is in stark contrast to
probes such as $N$-point functions of the CMB, which only constrain $N$-point
functions of $\zeta$ at the same order.

We now specialise to the simplest model of LPnG constrained by the CMB
datasets – in this model, the primordial curvature fluctuation can be
expressed as a polynomial (at most cubic) in a Gaussian random field as

$\zeta(\vec{x})=\mathcal{Z}(\vec{x})+\frac{3}{5}f_{\mathrm{NL}}(\mathcal{Z}(\vec{x})^{2}-\langle\mathcal{Z}^{2}\rangle)+\frac{9}{25}g_{\mathrm{NL}}(\mathcal{Z}(\vec{x})^{3}-3\langle\mathcal{Z}^{2}\rangle\mathcal{Z}(\vec{x}))\
.$ (2.13)

Current constraints on the parameters $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$
satisfy the condition for weak non-Gaussianity, i.e.
$g_{\mathrm{NL}}P_{\zeta},f_{\mathrm{NL}}P_{\zeta}^{1/2}\ll 1$, yet allow for
$g_{\mathrm{NL}}P_{\zeta}\sim f_{\mathrm{NL}}P_{\zeta}^{1/2}$. Under these
conditions, the only relevant primordial correlation functions in addition to
the power spectrum are the bispectrum and trispectrum, proportional to
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$,

$\displaystyle B_{\zeta}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3})$
$\displaystyle=$
$\displaystyle\frac{6}{5}f_{\mathrm{NL}}\left(P_{\zeta}(k_{1})P_{\zeta}(k_{2})+\dots
2\,\textrm{permutations}\right)\ ,$ (2.14) $\displaystyle
T_{\zeta}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3},\vec{k}_{4})$ $\displaystyle=$
$\displaystyle\frac{54}{25}g_{\mathrm{NL}}\left(P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(k_{3})+\dots
3\,\textrm{permutations}\right)$
$\displaystyle+\,\frac{36}{25}f_{\mathrm{NL}}^{2}\left(P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(|\vec{k}_{1}+\vec{k}_{3}|)+11\,\textrm{permutations}\right)\,.$

In our studies of $g_{\mathrm{NL}}$, we are working in the limit that
$g_{\mathrm{NL}}P_{\zeta}\sim f_{\mathrm{NL}}P_{\zeta}^{1/2}$ so that
$f_{\mathrm{NL}}\ll g_{\mathrm{NL}}$ and the second term in Eq. (2.1) can be
dropped. In the next section, we will generalise to two-field initial
conditions, which can enhance the $f_{\mathrm{NL}}^{2}$ term. The only
relevant non-trivial moments of the smoothed linear matter density field
($S_{N}$ defined in Eq. (2.3)) are then proportional to $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ and are as follows,

$\displaystyle S_{2}(R,z)$ $\displaystyle=$ $\displaystyle\sigma^{2}(R,z)\ ,$
(2.16) $\displaystyle S_{3}(R,z)$ $\displaystyle=$
$\displaystyle\frac{3}{5}f_{\mathrm{NL}}\Sigma_{3}(R,z)\ ,$ (2.17)
$\displaystyle S_{4}(R,z)$ $\displaystyle=$
$\displaystyle\frac{9}{25}g_{\mathrm{NL}}\Sigma_{4}(R,z)\ ,$ (2.18)

where the factors $\Sigma_{3}$ and $\Sigma_{4}$ can be obtained using Eq.
(2.3) and are given by:

$\Sigma_{3}=6\int\left(\prod_{i=1}^{2}\frac{d^{3}k_{i}}{(2\pi)^{3}}W(k_{i}R)T_{cb}(k_{i},z)P_{\zeta}(k_{i})\right)W(|k_{1}+k_{2}|R)T_{cb}(|k_{1}+k_{2}|,z)\
,$ (2.19)
$\Sigma_{4}=24\int\left(\prod_{i=1}^{3}\frac{d^{3}k_{i}}{(2\pi)^{3}}W(k_{i}R)T_{cb}(k_{i},z)P_{\zeta}(k_{i})\right)W(|k_{1}+k_{2}+k_{3}|R)T_{cb}(|k_{1}+k_{2}+k_{3}|,z)\
.$ (2.20)

Eqs. (2.4) and (2.5) show that the response of the moment $S_{N}$ to an
ambient long-wavelength curvature perturbation is proportional to the smoothed
squeezed limit of the $(N+1)$ linear matter correlation function. Since the
only non-trivial linear matter correlation functions (for the model given by
Eq. (2.13)) are the bispectrum and the trispectrum respectively, it follows
that the only non-trivial responses are those of $S_{2}$ and $S_{3}$ (i.e. of
the smoothed variance and skewness of the matter density field respectively).
From Eqs. (2.4), (2.16) and (2.17), we obtain:

$\displaystyle\frac{\partial S_{2}}{\partial\zeta_{L}}$ $\displaystyle=$
$\displaystyle\frac{12}{5}f_{\mathrm{NL}}S_{2}\ ,$ (2.21)
$\displaystyle\frac{\partial S_{3}}{\partial\zeta_{L}}$ $\displaystyle=$
$\displaystyle\frac{27}{25}g_{\mathrm{NL}}\Sigma_{3}\ ,$ (2.22)

where $\Sigma_{3}$ is given by Eq. (2.19). Substituting Eqs. (2.21) and (2.22)
in Eq. (2.12), we obtain the following equation for the scale-dependent bias
in this simple LPnG model:

$\left.\mathcal{B}(z)\right|_{f_{\mathrm{NL}},g_{\mathrm{NL}}}=4f_{\mathrm{NL}}S_{2}\cdot\frac{\partial\log
n_{g}}{\partial
S_{2}}+\frac{9}{5}g_{\mathrm{NL}}\Sigma_{3}\cdot\frac{\partial\log
n_{g}}{\partial S_{3}}\ .$ (2.23)

Using Eq. (2.17), we can express $\Sigma_{3}\left(\partial\log n_{g}/\partial
S_{3}\right)=(5/3)\times\partial\log n_{g}/\partial f_{\mathrm{NL}}$ and thus
rewrite the non-Gaussian contribution to the galaxy bias in its most often
quoted form:

$\displaystyle\Delta
b_{\text{NG}}=\frac{f_{\mathrm{NL}}\beta_{f}(z)+g_{\mathrm{NL}}\beta_{g}(z)}{\alpha(k,z)}\
,$ (2.24)

where the quantities $\beta_{f}(z)=2\partial\log n_{g}/\partial\log\sigma$ and
$\beta_{g}(z)=3\partial\log n_{g}/\partial f_{\mathrm{NL}}$ denote the
responses of the abundance of matter density tracers to the variance and
skewness of the smoothed matter density field, respectively. They are
dependent on the complicated physics of galaxy (or other tracer) formation and
may not be expressible analytically.

In a galaxy power spectrum forecast, we need to make modelling choices for the
parameters $\beta_{f}$ and $\beta_{g}$. A strong modelling choice is to assume
that the mean tracer density $n_{g}$ has the form of a universal halo mass
function – wherein $n_{g}$ is only a function of the ratio
$\delta_{c}/\sigma$, where $\delta_{c}$ is the threshold density for spherical
collapse in an Einstein-deSitter universe. In this case, one can show that
$\beta_{f}=2b_{L}\delta_{c}$. Modelling $\beta_{g}$ requires further
assumptions about how the tracer density $n_{g}$ depends on $f_{\mathrm{NL}}$.
For our SPHEREx forecasts (under strong modelling assumptions), we use the
following fitting formula derived in Ref. [33] (see their Eq. (51)):

$\displaystyle\beta_{g}=\kappa_{3}\left(-0.4(\nu-1)+1.5(\nu-1)^{2}+0.6(\nu-1)^{2}\right)\
,$ (2.25)

where $\kappa_{3}=0.000329(1+0.09z)b_{L}^{-0.09}$ and
$\nu=(1+0.5\beta_{f})^{1/2}$.

The universality assumption has been widely used to model the scale-dependent
galaxy bias induced by LPnG. However, hydrodynamical separate universe
simulations (see Refs. [42, 43]) have shown that the universality assumption
may not hold true for realistic galaxy populations. In particular, Ref. [42]
has shown that for galaxies selected by stellar mass, $\beta_{f}$ can be
modelled more accurately as $\beta_{f}=2\delta_{c}(b_{E}-p)$, where $0.4\leq
p\leq 0.7$, and $b_{E}$ is the (scale-independent) Eulerian bias. Inspired by
this result, we also obtain MCMC forecasts with a non-universal modelling
choice whereby $\beta_{f}=2\delta_{c}(b_{E}-0.5)$ and $\beta_{g}$ is given by
Eq. (2.25), with $\nu=(1+\delta_{c}(b_{E}-0.5))^{1/2}$. In Sec. 4.3, we study
the robustness of our forecast constraints to these choices by allowing $p$ to
be a free parameter.

If $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ are of comparable magnitude, then
our expressions for $\beta_{g}$ and $\beta_{f}$ show that the contribution of
$g_{\mathrm{NL}}$ to the scale-dependent bias signal is suppressed w.r.t.
contribution of $f_{\mathrm{NL}}$ by four orders of magnitude. In this
scenario, the primordial correlation functions scale hierarchically, e.g.
$S_{N+1}\sim S_{N}\sigma$, and scale-dependent bias is primarily probing the
amplitude of the largest non-trivial cumulant, $f_{\mathrm{NL}}$. In our
forecasts with non-zero $g_{\mathrm{NL}}$ we focus on the opposite limit, in
which $g_{\mathrm{NL}}\sim 10^{4}f_{\mathrm{NL}}$, leading to a primordial
kurtosis of comparable magnitude to the primordial skewness, $S_{3}\sim
S_{4}$, so that both terms contribute comparably to the scale-dependent bias.
This situation is consistent with current constraints on both parameters and
is also the only scenario in which $g_{\mathrm{NL}}$ can be detected in the
foreseeable future.

As $b_{E}$ is typically taken as a free parameter for each tracer population,
Eq. (2.24) shows that the parameters $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$
are degenerate with respect to their impact on the scale-dependent bias. In
fact, for a single redshift and single tracer, $f_{\mathrm{NL}}$ and
$(\beta_{g}/\beta_{f})g_{\mathrm{NL}}$ cannot be distinguished from each other
by galaxy power spectrum measurements alone. However, observations of multiple
tracer populations, which typically have different values of $\beta_{g}(z)$
and $\beta_{f}(z)$, may break this near-perfect degeneracy. In this paper, we
perform MCMC forecasts to quantify the possible degradation in the joint
constraints on both $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ that can be
obtained from the data collected by the SPHEREx all-sky survey [28].

###  2.2 Stochastic galaxy bias from local primordial non-Gaussianity

As stated before, any model of cosmic inflation that produces LPnG in the
primordial curvature perturbations necessarily needs to include at least one
light field in addition to the inflaton. A simple model that can produce non-
Gaussian primordial fluctuations compatible with current measurements of
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ incorporates a single massless field
called the curvaton [11, 12, 13] in addition to the inflaton, such that the
inflationary expansion of the universe is driven mainly by the dynamics of the
inflaton, whereas the primordial fluctuations are completely dominated by
fluctuations in the curvaton. If both the inflaton and the curvaton contribute
comparably to generating the primordial curvature perturbations, the
primordial statistics are changed. In the simplest scenario, the contribution
of the inflaton remains Gaussian, while the curvaton contribution remains non-
Gaussian and admits an expansion like Eq. (2.13), which is local in real space
in terms of a Gaussian random variable, i.e.:

$\displaystyle\zeta(\vec{x})$ $\displaystyle=$
$\displaystyle\zeta_{I}(\vec{x})+\zeta_{C}(\vec{x})\ ,$ (2.26)
$\displaystyle\zeta_{C}(\vec{x})$ $\displaystyle=$
$\displaystyle\mathcal{Z}_{C}(\vec{x})+\frac{3}{5}\overline{f_{\mathrm{NL}}}(\mathcal{Z}_{C}(\vec{x})^{2}-\langle\mathcal{Z}_{C}^{2}\rangle)+\frac{9}{25}\overline{g_{\mathrm{NL}}}(\mathcal{Z}_{C}(\vec{x})^{3}-3\langle\mathcal{Z}_{C}^{2}\rangle\mathcal{Z}(\vec{x}))+\dots\
,$ (2.27)

where $\zeta_{I}$ is the Gaussian contribution to the primordial curvature
perturbation that is sourced from the inflaton fluctuations, whereas the non-
Gaussian $\zeta_{C}$ is sourced from the fluctuations of the curvaton. The
quantities $\overline{f_{\mathrm{NL}}}$ and $\overline{g_{\mathrm{NL}}}$
characterise the amplitudes of the bispectrum and trispectrum of the
$\zeta_{C}$ field.

The correlation between the fields $\zeta_{I}$ and $\zeta_{C}$ is dependent on
the particular features of the inflationary model. For this work, we consider
a simple model in which $\zeta_{I}$ and $\zeta_{C}$ are uncorrelated and have
proportional power spectra. We define the contribution of the inflaton to the
primordial curvature perturbation w.r.t. the curvaton contribution by

$\xi^{2}=P_{\zeta_{I}}/P_{\zeta_{C}}\,,$ (2.28)

a scale-independent constant (for generalisations of this, see, e.g, Ref.
[19]). Assuming this simple model, the bispectrum and trispectrum can be
computed as,

$\displaystyle B_{\zeta}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3})$
$\displaystyle=$
$\displaystyle\frac{6}{5}\frac{\overline{f_{\mathrm{NL}}}}{(1+\xi^{2})^{2}}\left(P_{\zeta}(k_{1})P_{\zeta}(k_{2})+\dots
2\,\textrm{permutations}\right)\ ,$ (2.29) $\displaystyle
T_{\zeta}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3},\vec{k}_{4})$ $\displaystyle=$
$\displaystyle\frac{54}{25}\frac{\overline{g_{\mathrm{NL}}}}{(1+\xi^{2})^{3}}\left(P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(k_{3})+\dots
3\,\textrm{permutations}\right)$
$\displaystyle+\,\frac{36}{25}\frac{\overline{f_{\mathrm{NL}}}^{2}}{(1+\xi^{2})^{3}}\left(P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(|\vec{k}_{1}+\vec{k}_{3}|)+11\,\textrm{permutations}\right)\
.$

From above, we identify the parameters $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$
that characterize the bispectrum and trispectrum of $\zeta$ as,

$\displaystyle f_{\mathrm{NL}}$ $\displaystyle=$
$\displaystyle\frac{\overline{f_{\mathrm{NL}}}}{(1+\xi^{2})^{2}}\ ,$ (2.31)
$\displaystyle g_{\mathrm{NL}}$ $\displaystyle=$
$\displaystyle\frac{\overline{g_{\mathrm{NL}}}}{(1+\xi^{2})^{3}}\ .$ (2.32)

We see then that for a given level of non-Gaussianity in the $\zeta_{C}$
field, allowing $\zeta$ to have a contribution from Gaussian inflaton
fluctuations reduces the overall amplitude of non-Gaussianity by factors of
$P_{\zeta_{C}}/P_{\zeta}=1/(1+\xi^{2})$. On the other hand, the collapsed
($\mathcal{O}(f_{\mathrm{NL}}^{2})$) part of the trispectrum is enhanced
relative to the bispectrum by a factor of $(1+\xi^{2})$. Defining

$\displaystyle\tau_{\mathrm{NL}}$ $\displaystyle=$
$\displaystyle\frac{36}{25}\frac{\overline{f_{\mathrm{NL}}}^{2}}{(1+\xi^{2})^{3}}$
(2.33) $\displaystyle=$
$\displaystyle\left(\frac{6}{5}f_{\mathrm{NL}}\right)^{2}(1+\xi^{2})\ ,$
(2.34)

we see that $\tau_{\mathrm{NL}}$ is enhanced by contributions from the
inflaton [21]. The most recent Planck datasets [44] lead to the constraint
$\tau_{\mathrm{NL}}<2800$ at $95\%$ confidence – this constraint is consistent
with the assumption of weak primordial non-Gaussianity
($\langle\zeta^{N}\rangle/\langle\zeta^{2}\rangle^{N/2}\ll 1$). Under these
conditions, the only non-trivial higher point functions of the primordial
curvature perturbations (Eq. (2.26)) are the bispectrum and the trispectrum of
$\zeta_{C}$, which have the functional forms given by Eqs. (2.14) and (2.1),
respectively.

To analyse the impact of a primordial curvature perturbation of the kind given
by Eq. (2.26) on the galaxy number density fluctuations, we note that small-
scale correlation functions of the curvature perturbation $\zeta$ couples only
to an ambient long-wavelength mode of the non-Gaussian part of primordial
curvature perturbation, i.e. $\zeta_{C}$. In particular, Eq. (2.1) holds in a
slightly modified form where now only the non-Gaussian contribution to the
curvature, i.e. $\zeta_{C}$, modulates the small-scale correlation functions,

$\lim_{q\to
0}\frac{\langle\zeta_{C}(\vec{q})\zeta(\vec{k_{1}})...\zeta(\vec{k_{N}})\rangle}{P_{\zeta_{C}}(q)}=\frac{\partial\langle\zeta(\vec{k_{1}})...\zeta(\vec{k_{N}})\rangle}{\partial\zeta_{C,L}}\
,$ (2.35)

whereas for $N\geq 2$ the curvature perturbations sourced by the inflaton do
not modulate small-scale correlation functions,

$\displaystyle\frac{\partial\langle\zeta(\vec{k_{1}})..\zeta(\vec{k_{N}})\rangle}{\partial\zeta_{I,L}}$
$\displaystyle=$ $\displaystyle\lim_{q\to
0}\frac{\langle\zeta_{I}(\vec{q})\zeta(\vec{k_{1}})...\zeta(\vec{k_{N}})\rangle}{P_{\zeta_{I}}(q)}$
$\displaystyle=$ $\displaystyle 0\ .$

The moments of the smoothed small-scale linear matter density field
$S_{N}(R,z)$ (defined in Eq. (2.3)) in the presence of a long-wavelength
curvature perturbation $\zeta_{L}=\zeta_{I,L}+\zeta_{C,L}$ depend on
$\zeta_{C,L}$ similarly as Eq. (2.4),

$\frac{\partial
S_{N}}{\partial\zeta_{C,L}}=\int\prod^{N}_{i=1}\frac{d^{3}\vec{k_{i}}}{(2\pi)^{3}}W(k_{i}R)T_{cb}(k_{i},z)\lim_{q\to
0}\frac{\langle\zeta_{C}(\vec{q})\zeta(\vec{k_{1}})...\zeta(\vec{k_{N}})\rangle}{P_{\zeta_{C}}(q)}\
.$ (2.37)

From the arguments of the previous section, it follows therefore that the
galaxy number density fluctuations obey a relation similar to Eq. (2.7),
namely

$\displaystyle\frac{\delta n_{g}}{n_{g}}=\frac{\partial\log
n_{g}}{\partial\delta_{L}}\delta_{L}+\left(\sum^{\infty}_{N=2}\frac{\partial\log
n_{g}}{\partial S_{N}}\cdot\frac{\partial
S_{N}}{\partial\zeta_{C,L}}\right)\zeta_{C,L}\ ,$ (2.38)

where
$\delta_{L}=T_{cb}(k,z)\zeta_{L}(\vec{k})=T_{cb}(k,z)(\zeta_{I,L}+\zeta_{C,L})$
is the large scale (linear) density perturbation. Eq. (2.38) shows that the
galaxy number density contrast $\delta_{g}=\delta n_{g}/n_{g}$ is not fully
correlated with the ambient large scale matter density field $\delta_{L}$ – it
is also partially correlated with $\zeta_{C,L}$, which has fluctuations that
cannot be completely specified from the realisation of $\zeta$ or $\delta_{L}$
alone. This introduces stochasticity in the distribution of galaxies, relative
to the matter field, and is known as stochastic bias [20]. From Eq. (2.38),
the galaxy bias $b_{g}$ and the galaxy power spectrum $P_{gg}$ are given by

$\displaystyle b_{g}(k,z)$ $\displaystyle=$
$\displaystyle\frac{\langle\delta_{g}\delta_{L}\rangle}{\langle\delta_{L}\delta_{L}\rangle}$
(2.39) $\displaystyle=$ $\displaystyle\frac{\partial\log
n_{g}}{\partial\delta_{L}}+\frac{1}{1+\xi^{2}}\left(\sum^{\infty}_{N=2}\frac{\partial\log
n_{g}}{\partial S_{N}}\cdot\frac{\partial
S_{N}}{\partial\zeta_{C,L}}\right)\frac{1}{T_{cb}(k,z)}\ ,$ (2.40)

where, as before, $b_{L}=\partial\log n_{g}/\partial\delta_{L}$ is the scale-
independent, Lagrangian bias, while the non-Gaussian, scale-dependent
correction to the galaxy bias is

$\displaystyle\Delta b_{NG}=\frac{\mathcal{B}_{C}(z)}{\alpha(k,z)}\ ,$ (2.41)

where $\alpha(k,z)$ is defined in the previous section (Eq. (2.9)) and

$\displaystyle\mathcal{B}_{C}(z)=\frac{5/3}{(1+\xi^{2})}\sum^{\infty}_{N=2}\frac{\partial\log
n_{g}}{\partial S_{N}}\cdot\frac{\partial S_{N}}{\partial\zeta_{C,L}}\ .$
(2.42)

Equation (2.38) can now be written in a concise form to give the modified
linear bias relation for (Lagrangian) galaxy number density contrasts,

$\displaystyle\delta_{g}=b_{L}\delta_{L}+(1+\xi^{2})\Delta
b_{NG}(k,z)\delta_{C,L}\ ,$ (2.43)

where $\delta_{C,L}=T_{cb}(k,z)\zeta_{C,L}$ is the part of the density field
sourced by $\zeta_{C}$. The only non-trivial responses of the smoothed moments
of the density field (for the model given by Eqs. (2.26) and (2.27)) can be
obtained using Eq. (2.37) and are given by:

$\displaystyle\frac{\partial S_{2}}{\partial\zeta_{L}}$ $\displaystyle=$
$\displaystyle\frac{12}{5}f_{\mathrm{NL}}(1+\xi^{2})S_{2}\ ,$ (2.44)
$\displaystyle\frac{\partial S_{3}}{\partial\zeta_{L}}$ $\displaystyle=$
$\displaystyle\frac{27}{25}g_{\mathrm{NL}}(1+\xi^{2})\Sigma_{3}\ ,$ (2.45)

where $\Sigma_{3}$ is given by Eq. (2.19). Substituting Eqs. (2.44) and (2.45)
in Eq. (2.42), we see that the non-Gaussian correction to the galaxy bias
$\Delta b_{NG}$ in the two-field initial conditions given by Eq. (2.26) is the
same as Eq. (2.24). The galaxy power spectrum
$P_{gg}=\langle\delta_{g}\delta_{g}\rangle$ on the other hand can be shown to
be,888Note that
$\langle\delta_{C,L}\delta_{C,L}\rangle\propto\langle\delta_{L}\delta_{L}\rangle/(1+\xi^{2})$.

$\displaystyle P_{gg}(k,z)=\left(\left(b_{E}+\Delta
b_{NG}\right)^{2}+\xi^{2}\Delta b_{NG}^{2}\right)P_{cc}(k,z)\ ,$ (2.46)

where $b_{E}=1+b_{L}$ is the Eulerian galaxy bias, $P_{cc}(k,z)$ is the linear
CDM+Baryon power spectrum, and the reader is reminded that $\xi^{2}$
represents the ratio of inflaton-to-curvaton power, Eq. (2.28), so that
$\tau_{\mathrm{NL}}=(\frac{6}{5}f_{\mathrm{NL}})^{2}(1+\xi^{2})$. Equation
(2.46) is the main result of this section, which shows that the galaxy power
spectrum acquires a stochastic contribution proportional to the relative
contribution of the inflaton to the primordial curvature perturbation w.r.t.
the curvaton. We note here an interesting feature of this scenario is that the
different $k$-dependences of the $\Delta b_{NG}$ and $\Delta b_{NG}^{2}$ terms
in Eq. (2.46) allow the combination $(1+\xi^{2})$ to be extracted from
$P_{gg}$ even if the coefficient $\mathcal{B_{C}}$ in Eqs. (2.41) & (2.42)
cannot be determined, for instance due to the challenges associated with
modelling $n_{g}(S_{N})$.

For the stochastic bias signal to be large enough for a viable measurement of
$\tau_{\mathrm{NL}}$, we need the stochastic part of the galaxy power spectrum
to be comparable to its non-stochastic part. In other words, a robust
measurement of $\tau_{\mathrm{NL}}$ can be made using the stochastic halo bias
if (following Eq. (2.46)), $\xi^{2}\Delta b_{NG}/(b_{E}+\Delta
b_{NG})^{2}\sim\mathcal{O}(1)$ in some wavelength and redshift bin which has a
significant constraining power. Our choices of the fiducial values for
$\tau_{\mathrm{NL}}$ will be motivated by imposing this requirement at $z=0$
and at the largest scales accessible to the SPHEREx survey (i.e.
$k_{\text{min}}=0.001\ \text{h/Mpc}$).

In general, the parameter $\tau_{\mathrm{NL}}$, which parameterises the
collapsed limit of the primordial 4-point function, can be shown to satisfy
the Suyama-Yamaguchi inequality [21, 22]:
$\tau_{\mathrm{NL}}\geq(\frac{6}{5}f_{\mathrm{NL}})^{2}$. It can also be shown
(see, for example, Ref. [32]) that a significant stochasticity in the galaxy
bias arises whenever this inequality is strictly obeyed
$\tau_{\mathrm{NL}}>(\frac{6}{5}f_{\mathrm{NL}})^{2}$.999In general, it can be
shown that the collapsed limits of all higher point primordial correlation
function leads to stochasticity in the distribution of galaxies [32]; somewhat
analogous to how squeezed limits of all higher point primordial correlations
contribute to the scale-dependent bias. However, typical inflationary models
constrained by data suggest that the dominant contribution to halo-
stochasticity comes from the collapsed limit of the primordial 4-point
function whose magnitude is given by $\tau_{\mathrm{NL}}$ [32]. Reflecting
this fact, we choose to forecast measurements of $\tau_{\mathrm{NL}}$ and not
other collapsed-limit non-Gaussianity parameters in this paper. In this
section, we have shown how this is true for a simple model101010Note that even
a non-trivial $g_{\mathrm{NL}}$ can also lead to stochasticity in the galaxy
bias, but this stochasticity is proportional to $P_{\zeta}^{2}$ and is smaller
than the stochasticity due to $\tau_{\mathrm{NL}}$ by a factor of
$P_{\zeta}\sim 10^{-9}$ and can be safely ignored in the leading-order
calculation relevant for this paper. See, for example, Ref. [32] for a more
detailed analysis on how primordial non-Gaussianity can lead to stochastic
galaxy bias. (given by Eq. (2.26)), which we use to obtain forecasts for
$\tau_{\mathrm{NL}}$ for the SPHEREx all-sky survey.

##  3 Analysis Method

In this paper, we use the Bayesian likelihood pipeline developed in [45] to
forecast the sensitivity of the SPHEREx all-sky survey in jointly constraining
the LPnG parameters $f_{\mathrm{NL}},\ g_{\mathrm{NL}}\ \text{and}\
\tau_{\mathrm{NL}}$ around fiducial values consistent with state-of-the-art
constraints from CMB datasets. This is a standard Bayesian Markov Chain Monte
Carlo (MCMC) method used to obtain parameter forecasts for cosmological
surveys – in particular, we use a multitracer galaxy power spectrum likelihood
(developed in Ref. [45]) for a SPHEREx-like galaxy survey constructed within
the MontePython v3.4 MCMC sampling package [46, 47]. In this section, we give
a broad overview of our MCMC analysis method used for the forecasts obtained
in this paper and refer the reader to Ref. [45] for more details.

Our MCMC analysis involves the following steps:

  * •

Modelling the observables (in this case the galaxy auto- and cross-power
spectra at different scales and redshifts) for a given set of cosmological
parameters.

  * •

Obtaining the posterior likelihood for a given set of cosmological parameters
(using the modelling framework defined above) given a ‘mock data’ provided by
observables modelled in the fiducial cosmology.

  * •

Obtaining confidence intervals for the cosmological parameters by sampling the
posterior likelihood function using Markov Chain Monte Carlo (MCMC) methods.
Alternatively, one could obtain (potentially) less accurate confidence
intervals by computing a Fisher matrix. In this paper, unless otherwise
specified, we sample the posterior likelihood function using the Metropolis
Hastings MCMC sampling algorithm within the MontePython v3.4 MCMC sampling
package.

Obtained in this way, the mean cosmological parameters should be the same
(within numerical uncertainties) as the fiducial parameters, whereas the
confidence intervals provide the forecast sensitivity to the cosmological
parameters.

The SPHEREx all-sky survey [28] is designed to minimise systematics in
observations of galaxy clustering at large scales (going to
$k_{\text{min}}=0.001\ \text{h/Mpc}$). It will map a large cosmic volume and
measure spectroscopic redshifts of galaxies belonging to the all-sky
catalogues of the WISE [48], Pan-STARRS [49] and DES [50] galaxy surveys. It
is also a deep survey – observing galaxies in eleven redshift bins
$\left[z_{i}^{\text{min}},z_{i}^{\text{max}}\right]_{i=1}^{i=11}$ with mean
redshifts $\overline{z}_{i}$ going from $\overline{z}_{\text{min}}=0.1$ to
$\overline{z}_{\text{max}}=4.3$. Due to its low spectral resolution (as
compared to a true spectroscopic galaxy survey), SPHEREx is expected to
measure spectroscopic galaxy redshifts to varying degrees of precision and
will divide its observed galaxies into five galaxy populations according to
their redshift uncertainties. In this way, the data collected by SPHEREx will
lend itself to the use of multitracer techniques which can be used to obtain
significantly tighter constraints on primordial non-Gaussianity. The
observables to be obtained from SPHEREx are the galaxy power spectra and
galaxy bispectra in redshift space – in this paper, we study the constraining
power of the SPHEREx all-sky survey in constraining LPnG parameters using the
galaxy power spectra alone. As discussed, galaxy power spectra alone could
provide a first detection of any LPnG that could later be studied in more
detail through higher-order correlation functions. Note that since SPHEREx
will divide its observed galaxies into five galaxy populations, our
observables include the cross-power spectra as well the auto-power spectra of
different galaxy samples.

As the first step in our likelihood analysis, we need to model the galaxy
power spectra for a given set of cosmological parameters. To do so, we use the
linear bias relation ($\delta_{g}=b_{g}\delta_{cb}$) that is valid (at the
linear level) at late times and in the comoving synchronous gauge (see Eq.
(2.46)). According to the linear bias model, the theoretical galaxy power
spectrum in the presence of LPnG between the $i^{th}$ and $j^{th}$ SPHEREx
galaxy samples is given by

$\displaystyle P^{th}_{ij}(k,z)$ $\displaystyle=$ $\displaystyle
b_{g,i}b_{g,j}P_{cc}(k,z)+\xi^{2}\Delta b_{NG,i}\Delta b_{NG,j}P_{cc}(k,z)\ ,$
(3.1)

where $b_{g,i}$ and $\Delta b_{NG,i}$ are respectively the total galaxy bias
(including non-Gaussian corrections) and the non-Gaussian contribution (due to
both $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$) to the galaxy bias of the
$i^{th}$ SPHEREx galaxy population. As mentioned in the previous section,
$P_{cc}$ is the linear CDM+Baryon power spectrum. The parameter $\xi^{2}$ is
defined such that $\tau_{\mathrm{NL}}=(6/5f_{\mathrm{NL}})^{2}(1+\xi^{2})$ and
encodes the relative contribution of the inflaton w.r.t. the curvaton in the
primordial curvature fluctuations. As shown in Ref. [45] the sensitivity of
SPHEREx to $f_{\mathrm{NL}}$ is entirely driven by low-$k$ Fourier modes where
use of the linear power spectrum and linear bias factors are justified. Since
the $k$-dependence of the bias is identical for $f_{\mathrm{NL}}$,
$g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$, this is also justified here. We
use the Boltzmann code CLASS v3.2.0 [51, 52, 53] to compute the linear
CDM+Baryon power spectrum $P_{cc}(k,z)$, and model the non-Gaussian correction
$\Delta b_{NG}$ to the galaxy bias according to the analysis of Section 2 – in
particular following Eq. (2.24) for a cosmology with non-zero
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ and Eq. (2.46) for
$\tau_{\mathrm{NL}}\neq(\frac{6}{5}f_{\mathrm{NL}})^{2}$.

For the galaxy power spectrum forecasts obtained in this paper, we model the
nuisance parameter $\beta_{f}$ as

$\beta_{f}=2\delta_{c}(b_{E}-p)\,.$ (3.2)

Our first sets of forecasts will fix either $p=1$ and $p=0.5$ with $b_{E}$ a
free parameter for each galaxy sample and each redshift. We subsequently
explore consequences of treating $p$ as a free parameter for each SPHEREx
galaxy sample. The choice $p=1$ holds in the particular case when the mean
number density of observed galaxies has the same form as a universal halo mass
function. If this universality assumption does not hold true, the value of
$\beta_{f}$ (which in general is tracer-dependent) is difficult to estimate
from first principles as it encodes all the intractable details of the non-
linear physics of galaxy formation. However, the results of hydrodynamical
simulations within the separate universe framework conducted in Refs. [42, 43]
suggest that $\beta_{f}$ is related to the galaxy bias as
$\beta_{f}=(b_{E}-p)\delta_{c}$ with a tracer-dependent number $p$. Our choice
of $p=0.5$ is representative of modelling for $\beta_{f}$ following the
observations in Refs. [42, 43] that for galaxies selected by stellar mass,
$0.4\leq p\leq 0.7$ is a better fit to the results of hydrodynamical separate
universe simulations. SPHEREx galaxy samples are not selected by stellar mass,
so we forecast two choices of $p$ to illustrate the $p$-dependence of
constraints.

For modelling $\beta_{g}$, we use the fitting formula Eq. (2.25) derived in
Ref. [33] which relates $\beta_{g}$ to $\beta_{f}$ through
$\nu=(1+0.5\beta_{f})^{1/2}$, with $\beta_{f}$ as in Eq. (3.2). Note that our
first-round forecasts are obtained using these strong modelling choices that
assume a perfect knowledge of the functional dependence of the nuisance
parameters $\beta_{f}$ and $\beta_{g}$ on $b_{E}$. We also study forecasts
under the weaker modelling choice where $p$ a free, but redshift-independent,
parameter. If $p$ were completely free (i.e. allowed to vary independently at
different redshifts), then in the limit of large uncertainty on $p$,
$\beta_{f}$ and $\beta_{g}$ would become free-parameters independent of
$b_{E}$. In the case where there is large uncertainty on the form of
$\beta_{f}$ and $\beta_{g}$, $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ cannot be
distinguished from scale-dependence bias alone. For our free-$p$ analysis the
parameters become sufficiently degenerate that we restrict to forecasting
$f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ only. In this limit constraints on
$f_{\mathrm{NL}}$ can be thought of as representing constraints on an
effective amplitude that combines LPnG and tracer-dependent bias parameters,
such as the $\mathcal{B}$ parameter in Eq. (2.12), while
$(1+\xi^{2})=\tau_{\mathrm{NL}}/(\frac{6}{5}f_{\mathrm{NL}})^{2}$ should, in
principle, still be constrained even with completely free $\beta_{f}$.
Throughout our analysis, we use $\delta_{c}=1.42$ because it shows better
agreement between the above fitting formulas and N-body simulations [33].

The galaxy power spectrum observed by a given galaxy survey differs from the
theoretical galaxy power spectrum modelled according to Eq. (3.1) in several
ways. The peculiar velocities of galaxies along the line-of-sight break the
isotropy of the matter power spectrum and cause the observed galaxy power
spectrum to only be isotropic around the line of sight. This effect is known
as redshift-space distortion and can be modelled at the linear and mildly non-
linear scales relevant for the SPHEREx survey by making the replacement
$b_{g,i}\rightarrow b_{g,i}+f(z)\mu^{2}$ and $b_{g,j}\rightarrow
b_{g,j}+f(z)\mu^{2}$ in Eq. (3.1); where $f(z)$ is the linear growth rate of
small-scale matter density perturbations $f(z)=d\log\delta_{cb}/d\log a$, and
$\mu$ is the cosine of the angle between the wave vector $\vec{k}$ and the
line of sight (this is the well-known Kaiser term [54]). In addition to
redshift-space distortion effects, there are additional multiplicative factors
needed to map Eq. (3.1) to a model of the observed galaxy power spectra. Our
final expression for the observed galaxy power spectrum between the $i^{th}$
and $j^{th}$ galaxy samples of SPHEREx is (see Ref. [45])

$P_{ij}=\large(f_{BF}\times(f_{\sigma_{z,i}}f_{\sigma_{z,j}})^{1/2}\times
f_{AP}\large[(b_{g,i}+f\mu^{2})(b_{g,j}+f\mu^{2})P_{cc}+\xi^{2}\Delta
b_{NG,i}\Delta b_{NG,j}P_{cc}\large)\large]\ ,$ (3.3)

where $f_{AP}(z)$ denotes the contribution of the Alcock-Paczynski effect and
$f_{BF}(k,z)$ encodes the effect of non-linear bulk flows. The Alcock-
Paczynski term corrects for the error in the galaxy power spectrum measured
from data assuming a fiducial cosmology which is different from the true
cosmology. The term $f_{AP}(z)$ relates the comoving volume elements in the
fiducial and true cosmologies and if given by :

$\displaystyle
f_{AP}(z)=\frac{H(z)D^{\text{fid}}_{A}(z)^{2}}{H^{\text{fid}}(z)D_{A}(z)^{2}}\
.$ (3.4)

The term $f_{BF}(k,z)$ models the smearing of (small-scale) BAO features in
the galaxy power spectrum due to the effect of non-linear bulk flows [28].
Following Ref. [28],

$\displaystyle
f_{BF}(k,z)=\exp\left(-\frac{1}{2}k^{2}\Sigma_{\perp}^{2}-\frac{1}{2}k^{2}\mu^{2}(\Sigma_{||}^{2}-\Sigma_{\perp}^{2})\right)\
.$ (3.5)

The Lagrangian displacement fields $\Sigma_{\perp}$ and $\Sigma_{||}$ are
given by [28]

$\displaystyle\Sigma_{\perp}(z)$ $\displaystyle=$ $\displaystyle
c_{\text{rec}}D(z)\Sigma_{0}\ ,$ (3.6) $\displaystyle\Sigma_{||}(z)$
$\displaystyle=$ $\displaystyle c_{\text{rec}}D(z)(1+f(z))\Sigma_{0}\ ,$ (3.7)

where $D(z)$ is the linear growth factor and $f(z)$ is the linear growth rate.
We set the parameter $c_{\text{rec}}=0.5$ [28] and $\Sigma_{0}=11\
\text{h/Mpc}$ for a fiducial $\sigma_{8}=0.8$. The terms
$f_{\sigma_{z,i}}(k,z),\ f_{\sigma_{z,j}}(k,z)$ are exponential damping
factors that suppress the observed galaxy power spectrum signal at small
scales along the line of sight due to the uncertainty in the determination of
galaxy redshifts (namely $\sigma_{z,i}\ \text{and}\ \sigma_{z,j}$ which denote
redshift uncertainties in the $i^{th}$ and $j^{th}$ galaxy populations to be
observed by SPHEREx). For the $i^{th}$ SPHEREx galaxy sample,
$f_{\sigma_{z,i}}(k,z)=\exp(-k^{2}\mu^{2}\sigma_{z,i}^{2}/H(z)^{2})$, where
$H(z)$ is the Hubble rate at redshift $z$ and $\sigma_{z,i}$ is the redshift
uncertainty of the $i^{th}$ galaxy sample. We refer the reader to Refs. [45,
28] for a detailed discussion of these additional terms and their importance
in the context of the SPHEREx all-sky survey.

Proceeding further, we treat the observed linear galaxy number density
fluctuations at a given (fiducial) k-mode and redshift as Gaussian random
variables [45] with covariance (assuming Poisson shot noise) given by

$\displaystyle\langle\delta^{\dagger}_{g,i}\delta_{g,j}\rangle=P_{ij}(\vec{k},z)+\frac{\delta^{K}_{ij}}{\overline{n}_{g,i}(z)}\
,$ (3.8)

where $P_{ij}$ denotes the observed the cross-power spectrum of the $i^{th}$
and $j^{th}$ galaxy populations of SPHEREx (see Eq. (3.3)) and
$\overline{n}_{g,i}$ is the mean number density of observed galaxies in the
$i^{th}$ galaxy population of SPHEREx at redshift $z$. This translates into an
exponential likelihood for the galaxy power spectra at a given k-mode and
redshift given by

$\displaystyle\mathcal{L}_{\textbf{k},z}=\mathcal{J}\frac{1}{\det\mathcal{C}}\cdot\exp\left[-\text{Tr}\left(\mathcal{C}^{-1}\mathcal{D}\right)\right]\
,$ (3.9)

where the proportionality constant $\mathcal{J}$ is independent of the sampled
cosmological parameters. In Eq. (3.9), both $\mathcal{C}$ and $\mathcal{D}$
are the covariance matrices of galaxy number density fluctuations (at a given
k-mode and redshift) that respectively depend on the sampled, s, and fiducial,
fid, cosmological parameters in an MCMC run:

$\displaystyle\mathcal{C}_{ij}$ $\displaystyle=$ $\displaystyle
P^{s}_{ij}+\frac{\delta^{K}_{ij}}{\overline{n}_{g,i}}\ ,$ (3.10)
$\displaystyle\mathcal{D}_{ij}$ $\displaystyle=$ $\displaystyle
P^{fid}_{ij}+\frac{\delta^{K}_{ij}}{\overline{n}_{g,i}}\ .$ (3.11)

The entire SPHEREx power spectrum likelihood is the product of
$\mathcal{L}_{\textbf{k},z}$ over all redshifts bins and independent k-modes.
As stated before, the SPHEREx survey is a large-scale survey [28] with the
observed k-modes between $k_{\text{min}}=0.001\ \text{h/Mpc}$ and
$k_{\text{max}}=0.2\ \text{h/Mpc}$ and eleven redshift bins with mean
redshifts going from $\overline{z}_{\text{min}}=0.1$ to
$\overline{z}_{\text{max}}=4.3$. Note that for the purposes of obtaining
power-spectrum forecasts, we assume that different redshift bins and
independent k-modes are uncorrelated – this is a standard approximation used
in forecasting the sensitivity of cosmological surveys (see Refs. [28, 55,
45]). An MCMC run used to obtain cosmological forecasts typically maximises
the likelihood by minimising $\chi^{2}$ defined as

$\displaystyle\chi^{2}=-2\sum_{z_{i}}V(z_{i})\int_{k_{\text{min}}}^{k_{\text{max}}}\int_{-1}^{1}\frac{k^{2}dkd\mu}{2(2\pi^{2})}\log\mathcal{L}_{\vec{k},z}\
,$ (3.12)

where $V(z_{i})$ is the comoving volume surveyed in the $i^{th}$ redshift bin.

Having thus obtained the likelihood, we use the boosted MCMC sampler
MontePython v3.4 to conduct MCMC runs to fit for the following cosmological
parameters: {$\omega_{b}$ , $\omega_{cdm}$, $100\theta_{s}$, $A_{s}$, $n_{s}$,
$z_{reio}$, $M_{\nu}$}, as well as combinations of {$f_{\mathrm{NL}}$,
$g_{\mathrm{NL}}$, $\tau_{\mathrm{NL}}$}. In addition, we marginalise over the
scale-independent, Gaussian parts of the galaxy biases (given by
$b_{0,i}=b_{g,i}-\Delta b_{NG,i}$) as well as the bulk flow parameter
$\Sigma_{0}$ which parametrises the effect of non-linear bulk flows given by
$f_{BF}(k,z)$ in Eq. (3.3) and Eqs. (3.5)–(3.7) (see Refs. [28, 45] for more
details). During our MCMC analysis, we combine our power spectrum likelihood
(Eq. (3.12), (3.9)) with a mock CMB power spectrum likelihood named
‘fake_planck_realistic’ in MontePython – this likelihood is designed to mimic
the sensitivity of the Planck 2018 dataset in constraining cosmological
parameters (see Ref. [56] for more details). This mock CMB likelihood does not
include the CMB lensing power spectrum – this ensures that it is not
significantly correlated with our SPHEREx galaxy power spectrum likelihood.
Moreover, the mock CMB likelihood does not include the effect of primordial
non-Gaussianity and does not by itself provide any constraints on LPnG
parameters – this makes it equivalent to a Planck-based prior imposed in order
to resolve degeneracies between $\nu\Lambda CDM$ cosmological parameters (as
is done in Ref. [28]).

We use a $\nu\Lambda CDM$ cosmology with fiducial parameters
$\omega_{b}=0.02218$, $\omega_{\text{cdm}}=0.1205$,
$100\times\theta_{s}=1.04113$, $A_{s}=2.032692\times 10^{-9}$, $n_{s}=0.9613$,
$z_{reio}=7.68$, $M_{\nu}\equiv\Sigma_{i}m_{\nu i}=0.06$ eV. For simplicity
neutrinos are implemented as three degenerate species (see, e.g, Ref. [57]).
In this paper, we forecast constraints on non-Gaussian initial conditions with
primordial trispectra that contribute comparably in magnitude to galaxy power
spectra as primordial bispectra with $f_{\mathrm{NL}}\sim\mathcal{O}(1)$. To
achieve this scenario in which the primordial kurtosis is comparable to the
skewness we choose fiducial values of $f_{\mathrm{NL}}=1.0$ and
$g_{\mathrm{NL}}=1.0\times 10^{4}$. For the two-source initial conditions
described in Section 2.2 that produce a $\tau_{\mathrm{NL}}$ not specified by
$f_{\mathrm{NL}}$ alone, we choose to obtain forecasts around two fiducial
values of $\tau_{\mathrm{NL}}$, namely $\tau_{\mathrm{NL}}=1.3\times 10^{2}$
and $\tau_{\mathrm{NL}}=1.3\times 10^{3}$, well above our fiducial value of
$f_{\mathrm{NL}}=1$ and consistent with CMB datasets [44].

##  4 Results

$p$ | $\sigma(f_{\mathrm{NL}})$ | $\sigma(g_{\mathrm{NL}})$
---|---|---
1.0 | 1.06 | $0.92\times 10^{4}$
0.5 | 0.89 | $0.66\times 10^{4}$
Table 1: Single-parameter LPnG forecasts for $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ obtained from the SPHEREx multitracer likelihood for the
$p=1$ and $p=0.5$ modelling choices (see Eq. (2.25) and Eq. (3.2)).

We start by presenting single-parameter forecasts for $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ which parametrise the squeezed-limit primordial bispectrum
and trispectrum respectively. Table 1 shows the SPHEREx multitracer MCMC
forecasts for the individual parameters $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ around their respective fiducial values
$f_{\mathrm{NL}}=1.0$, and $g_{\mathrm{NL}}=1.0\times 10^{4}$ obtained within
the framework of the two-field inflationary model. Note that for the two-field
inflationary model considered in this paper, it doesn’t make sense to obtain
individual forecasts for $\tau_{\mathrm{NL}}$ because $\tau_{\mathrm{NL}}$
cannot be non-vanishing unless $f_{\mathrm{NL}}\neq 0$. The forecasts in Table
1 are obtained under strong modelling assumptions (see Section 3, Eqs. (2.24)
and (2.25) and surrounding text), where $\beta_{f}=\delta_{c}(b_{E}-p)$ and
$\nu=(1+\delta_{c}(b_{E}-p))^{1/2}$ with $p=0.5$ or $p=1$ for all galaxy
samples, yet $b_{E}$ is a free parameter for each galaxy population at each
redshift. As expected, the forecast constraints on $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ depend on the assumed value of $p$. Our individual forecasts
for $f_{\mathrm{NL}}$ or $g_{\mathrm{NL}}$ alone are also consistent with the
observation at the end of Section 2.1 that the contribution of
$g_{\mathrm{NL}}$ to the scale-dependent bias signal is suppressed w.r.t. that
of $f_{\mathrm{NL}}$ by four orders of magnitude, due to the relative
magnitudes of the dimensionless cumulants
$\langle\delta_{cb}^{3}\rangle/\langle\delta_{cb}^{2}\rangle\sim
f_{\mathrm{NL}}\sigma$ and
$\langle\delta_{cb}^{4}\rangle_{c}/\langle\delta_{cb}^{2}\rangle\sim
g_{\mathrm{NL}}\sigma^{2}$ for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ of
similar size.

Having established the forecast constraints on either $f_{\mathrm{NL}}$ or
$g_{\mathrm{NL}}$ alone, we now consider combined constraints. We conduct MCMC
runs to obtain joint forecasts for the pairs of parameters
$(f_{\mathrm{NL}},g_{\mathrm{NL}})$ and
$(f_{\mathrm{NL}},\tau_{\mathrm{NL}})$. Expanding this to include all three
parameters is a challenging and computationally intensive sampling problem due
to strong degeneracies between the parameters. We expect the sensitivity would
be washed out and barely more illuminating than the two-parameter case, since
we would only see notable constraining power for the parameter(s) with the
greatest contribution to the observable signal. The exception to this is the
somewhat fine-tuned case where all three parameters contribute comparably, in
which case we would expect decreased sensitivity due to the strong three-way
degeneracy. For the subcase with two parameters giving comparable
contributions to the scale-dependent bias, we expect the result would largely
asymptote to that of the two-parameter case, possibly with slightly decreased
sensitivity to the parameters. Given this expectation, we do not consider it
worth the significant computational resources required (and associated carbon
footprint) and omit the three-parameter case.

###  4.1 Joint forecasts on $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$

To obtain joint forecasts for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$, we
conduct MCMC runs with the galaxy bias modelled as Eq. ((2.24)), with
$\beta_{f}$ and $\beta_{g}$ as in Eq. (3.2) and Eq. (2.25), and simultaneously
fit for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$. Table 2 shows our joint MCMC
forecasts for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ as well as their
covariance around fiducial values $f_{\text{NL,fid}}=1.0$ and
$g_{\text{NL,fid}}=1\times 10^{4}$. Our fiducial values for $f_{\mathrm{NL}}$
and $g_{\mathrm{NL}}$ are consistent with the most recent CMB constraints [18]
and are such that the dimensionless skewness and kurtosis of the primordial
curvature ($\langle\zeta^{3}\rangle/\langle\zeta^{2}\rangle^{3/2}$ and
$\langle\zeta^{4}\rangle/\langle\zeta^{2}\rangle^{2}$) are of similar order of
magnitude. As mentioned in Section 2.1, both $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ contribute comparably to the observed signal of primordial
non-Gaussianity (i.e. the scale-dependent bias) in this regime. We again
consider two representative values of $p$, $p=1$ and $p=0.5$. Controlling for
the difference in the order of magnitude of $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$, Table 2 shows that we obtain nearly equal forecast
constraints for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$, analogous to the
individual forecasts in Table 1. However, these forecasts are significantly
weaker than their respective individual forecasts. This degradation in
forecasts for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ can be attributed to the
strong degeneracy between $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$, which is
evident from Eq. (2.24). Both terms produce the same $k$-dependence in the
bias, the terms differ only in how the amplitude of the $k$-dependent terms
vary with tracer population through $\beta_{f}$ and $\beta_{g}$, which are
themselves degenerate through the dependence on $b_{E}$. This degeneracy is
reflected in Figure 1, which shows the joint $68\%$ and $95\%$ confidence
regions for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ obtained from our MCMC
analysis. The degeneracy between $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ is
also reflected in their covariance (see Table 2), which is $\sim-0.9$ for both
the $p=1$ and $p=0.5$, with only slight variations in the covariance between
the two $p$-values (visible in Figure 1) due to the small changes in the
relationship between $\beta_{g}$ and $b_{E}$ in each case. This is consistent
with the Fisher matrix analysis done by Ref. [58], where it was found that
while multitracer analysis using scale-dependent halo bias can distinguish
between $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$, their covariance is always
close to $-1$ due to their degeneracy w.r.t. the scale-dependent bias.

Figure 1: Joint MCMC forecasts around fiducial $f_{\mathrm{NL}}=1.0$ and
fiducial $g_{\mathrm{NL}}=1\times 10^{4}$ obtained from the galaxy auto- and
cross-power spectra of all five SPHEREx galaxy samples. The two sets of
contours illustrate constraints for two modelling choices for the scale-
dependent bias, $p=1$ and $p=0.5$.

As shown in Table 1, our SPHEREx multitracer forecasts under non-Gaussian
initial conditions for only $f_{\mathrm{NL}}\neq 0$ or only
$g_{\mathrm{NL}}\neq 0$ are such that $\Delta
f_{\mathrm{NL}}/f_{\mathrm{NL}}\sim\Delta g_{\mathrm{NL}}/g_{\mathrm{NL}}\sim
1$ for fiducial values $f_{\mathrm{NL}}=1.0$ and $g_{\mathrm{NL}}=1.0\times
10^{4}$. This means that with systematic uncertainties eliminated and in the
scenario where either $f_{\mathrm{NL}}\neq 0$ or $g_{\mathrm{NL}}\neq 0$,
SPHEREx power spectrum measurements alone can potentially detect LPnG (at the
small amplitudes here) with at most $1\sigma$ significance (with the modelling
assumptions made here). It is worth noting that this continues to hold even
with both $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ being as large as their
fiducial values – in fact, the prospect of SPHEREx detecting LPnG with
$1\sigma$ significance is slightly improved when both are non-zero. For
forecasts obtained with $p=1$, the point $f_{\mathrm{NL}}=g_{\mathrm{NL}}=0$
is now disfavoured with $\Delta\chi^{2}=4.6$, which puts it far outside the
$68\%$ confidence region contour (as compared to the forecasts made for a
single parameter in which case the point $f_{\mathrm{NL}}=0$ is nearly at the
boundary of the 68% confidence interval). On the other hand, the point
$f_{\mathrm{NL}}=g_{\mathrm{NL}}=0$ is disfavoured with $\Delta\chi^{2}=6.5$
for forecasts obtained using the modelling assumption $p=0.5$ – suggesting
that SPHEREx power spectrum measurements alone can potentially detect LPnG at
the (small) amplitudes considered here with as much as $2\sigma$ significance
if SPHEREx galaxies are better characterized by
$p=0.5$.111111$\Delta\chi^{2}=2.3(6.17)$ for the boundary of the
$68.3\%(95.4\%)$ confidence region.

$p$ | $\sigma(f_{\mathrm{NL}})$ | $\sigma(g_{\mathrm{NL}})$ | $\text{Cov}(f_{\mathrm{NL}},g_{\mathrm{NL}})$
---|---|---|---
1.0 | 2.54 | $2.53\times 10^{4}$ | $-0.93$
0.5 | 2.86 | $2.10\times 10^{4}$ | $-0.95$
Table 2: Joint MCMC forecast for $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$
obtained from the SPHEREx multitracer likelihood for $p=1$ and $p=0.5$
modelling choices.

###  4.2 Joint forecasts on $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$

(a) Fiducial $f_{\mathrm{NL}}=1.0$ and fiducial $\tau_{\mathrm{NL}}=1.3\times
10^{2}$ (b) Fiducial $f_{\mathrm{NL}}=1.0$ and fiducial
$\tau_{\mathrm{NL}}=1.3\times 10^{3}$

Figure 2: Joint MCMC forecasts for $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$
(obtained around the fiducial values assumed in this paper) from the scale-
dependent bias of all five SPHEREx galaxy samples for $p=1$ and $p=0.5$
modelling choices.

We now consider joint constraints on $f_{\mathrm{NL}}$ and
$\tau_{\mathrm{NL}}$, with $g_{\mathrm{NL}}=0$. To obtain joint multitracer
forecasts for $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$, we conduct two MCMC
runs with fiducial values $f_{\mathrm{NL}}=1,\ \tau_{\mathrm{NL}}=1.3\times
10^{2}$ and $f_{\mathrm{NL}}=1,\ \tau_{\mathrm{NL}}=1.3\times 10^{3}$
respectively. As mentioned in Section 2, these fiducial values are well within
the current constraints, yet produce a large stochastic component of the
galaxy bias. Note that while conducting our MCMC runs, we need to impose the
condition $\tau_{\mathrm{NL}}\geq(6/5f_{\mathrm{NL}})^{2}$ following the
Suyama-Yamaguchi inequality [21].

Table 3 shows the joint MCMC forecast around fiducial values
$f_{\mathrm{NL}}=1,\ \tau_{\mathrm{NL}}=1\times 10^{2}$ and
$f_{\mathrm{NL}}=1,\ \tau_{\mathrm{NL}}=1\times 10^{3}$, respectively, for
$p=1.0$ and $p=0.5$. Figure 2 shows the joint $68\%$ and $95\%$ contours for
$f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ around their respective fiducial
values for examples with $p=1$ and $p=0.5$. In all cases, the covariance
between $\tau_{\mathrm{NL}}$ and $f_{\mathrm{NL}}$ is $\approx-0.6$ – this is
consistent with the fact that $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ are
only partly degenerate. This is because the contribution of
$\tau_{\mathrm{NL}}$ to the galaxy bias (unlike that of $g_{\mathrm{NL}}$)
goes as $\Delta b_{NG}^{2}$ ($\propto 1/k^{4}$ on large scales, see Eq.
(2.46)) and does not have the same scale-dependence as the contribution of
$f_{\mathrm{NL}}$ ($\propto 1/k^{2}$ on large scales). The contour plots in
Figure 2 as well as the covariances in Table 3 indeed show that while
$f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ are partly degenerate, they are not
as degenerate as $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$. We also see that the
choice of $p$ does not appear to effect the degeneracy between
$f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$. Comparing the left and right
panels of Figure 2 or the two sub-tables in Table 3, we see that larger
fiducial values of $\tau_{\mathrm{NL}}$ lead to tighter constraints on
$\tau_{\mathrm{NL}}$. This is, perhaps, unsurprising as a large
$\tau_{\mathrm{NL}}$ should lead to a larger signal. On the other hand, we see
that constraints on $f_{\mathrm{NL}}$ are weakened in the presence of larger
$\tau_{\mathrm{NL}}$. Our MCMC forecasts lead to the observation that in the
absence of systematics, SPHEREx can potentially constrain large values of
$\tau_{\mathrm{NL}}$ more precisely at the cost of degraded constraints on
$f_{\mathrm{NL}}$. In particular, in the scenario with
$\tau_{\mathrm{NL}}=1\times 10^{3}$, SPHEREx could detect non-trivial
$\tau_{\mathrm{NL}}$ at $\approx 5\sigma$ significance while not, at the same
time, having unambiguously detected $f_{\mathrm{NL}}\neq 0$. This means that
SPHEREx could infer the presence of an additional light field during inflation
(and thus indirectly rule out $f_{\mathrm{NL}}=0$) through a robust
measurement of scale-dependent stochastic bias, if the amplitude is
sufficiently large.

$p$ | $\sigma(f_{\mathrm{NL}})$ | $\sigma(\tau_{\mathrm{NL}})$ | $\text{Cov}(f_{\mathrm{NL}},\tau_{\mathrm{NL}})$
---|---|---|---
1.0 | 1.79 | $0.78\times 10^{2}$ | $-0.59$
0.5 | 1.64 | $0.67\times 10^{2}$ | $-0.61$
(a) Fiducial $f_{\mathrm{NL}}=1.0$ and fiducial $\tau_{\mathrm{NL}}=1.3\times
10^{2}$

$p$ | $\sigma(f_{\mathrm{NL}})$ | $\sigma(\tau_{\mathrm{NL}})$ | $\text{Cov}(f_{\mathrm{NL}},\tau_{\mathrm{NL}})$
---|---|---|---
1.0 | 2.42 | $0.24\times 10^{3}$ | $-0.56$
0.5 | 2.22 | $0.21\times 10^{3}$ | $-0.59$
(b) Fiducial $f_{\mathrm{NL}}=1.0$ and fiducial $\tau_{\mathrm{NL}}=1.3\times
10^{3}$

Table 3: Joint MCMC forecast for $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$
obtained from the SPHEREx multitracer likelihood. For each fiducial value of
$\tau_{\mathrm{NL}}$, we consider two example values of $p=1$ and $p=0.5$,
entering the model for $\beta_{f}$ as in Eq. (3.2).

###  4.3 Impact of modelling choices

Figure 3: Joint Fisher forecasts for $f_{\mathrm{NL}}$ and
$\tau_{\mathrm{NL}}$ around $f_{\mathrm{NL}}=1$ with
$\tau_{\mathrm{NL}}=1.3\times 10^{2}$ (left) and $\tau_{\mathrm{NL}}=1.3\times
10^{3}$ (right) and $\beta_{f}=\delta_{c}(b_{E}-p)$ as a function of priors on
$p$. The red plots show forecasts obtained with fiducial $p=0.5$ and the black
plots show forecasts obtained with fiducial $p=1.0$.

Our single-parameter, as well as joint, MCMC forecasts for $f_{\mathrm{NL}},\
g_{\mathrm{NL}},\ \text{and}\ \tau_{\mathrm{NL}}$ (see Tables 1, 2, 3) make it
amply clear that the constraining power of any galaxy survey in constraining
LPnG parameters using the scale-dependent galaxy bias varies with the
modelling of the nuisance parameters $\beta_{f}$ and $\beta_{g}$ (see
subsection 2.1). Given that $\beta_{f}$ and $\beta_{g}$ are strongly dependent
on tracer population, they are difficult to model from first principles and
present a big limiting factor to fully exploiting galaxy power spectrum data
to constrain LPnG. The non-trivial relation between $\beta_{f}$ and the
Eulerian galaxy bias with Gaussian initial conditions has been investigated
with hydrodynamical separate universe simulations (see Refs. [42, 43]), which
provide some evidence that supports modelling $\beta_{f}$ as
$\beta_{f}=\delta_{c}(b_{E}-p)$, with $p$ being a tracer-dependent number.
However, modelling $\beta_{g}$ in terms of other observable parameters, such
as $b_{E}$ is a much more intractable problem and represents a bigger
modelling ambiguity. The is due both to the complicated nature of modelling
$n_{g}$ from first principles, and the observation that even in the simplest
models, $\beta_{g}$ depends non-linearly on $b_{E}$, and the $\beta_{g}$ for
each sample will therefore depend on how tracers are selected and grouped
together in data analyses [33].

In this paper, we have so far obtained forecasts for the constraining power of
the SPHEREx galaxy surveys under relatively strong modelling choices, which
assume a complete knowledge of the relation between $\beta_{f},\ \beta_{g}$
and $b_{E}$. In particular, we assume (see Section 3)
$\beta_{f}=\delta_{c}(b_{E}-p)$ and $\beta_{g}$ is given by Eq. (2.25) with
$\nu=(1+0.5\beta_{f})^{1/2}$. Different tracer populations could have
different values of $p$, which can change our MCMC forecasts proportionally
(i.e. smaller p gives tighter constraints) given that the contributions of
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ are proportional to $\beta_{f}$ and
$\beta_{g}$ respectively. We have so far presented constraints for $p=1.0$ and
$p=0.5$ as representative values that reflect the universal halo mass function
assumption and some results of hydrodynamical separate universe simulations,
respectively. A more realistic approach is to make $\beta_{f}$ and $\beta_{g}$
free parameters and impose heuristic priors on them in our analysis. This can
potentially degrade our forecasts, but may give a more complete picture of how
well a given galaxy survey can constrain LPnG using the galaxy power spectrum.

If $\beta_{f}$ and $\beta_{g}$ are free parameters, Eq. (2.24) illustrates
that they are completely indistinguishable from scale-dependent bias
measurements alone. In other words, without predictions for the values of
$\beta_{f}$ and $\beta_{g}$, the galaxy power spectrum loses its utility as a
signal that can potentially distinguish between a primordial bispectrum and
trispectrum. This situation is not all bad, however: it is still true that in
the presence of squeezed-limits of higher-point primordial correlation
functions, SPHEREx could potentially detect a non-trivial LPnG at high
significance, yet with uncertainty in which squeezed-limit correlation
functions were contributing. This is because the scale-dependent bias
measurement is a measurement of a redshift-dependent weighted combination of
the parameters $f_{\mathrm{NL}}$, $g_{\mathrm{NL}}$ together with all the
higher order coefficients in the expansion (see Eq. (2.13)). For instance,
non-trivial squeezed limits of higher-point primordial correlation functions
can potentially contribute to increasing the observable signal of LPnG – our
results show that this happens around the fiducial values of $f_{\mathrm{NL}}$
and $g_{\mathrm{NL}}$ consistent with current constraints from CMB datasets
chosen in this paper.121212Note that in principle, it is possible for
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ to be of opposite signs and be tuned
such that their respective contributions to the scale-dependent bias are
(nearly) equal and of opposite signs – this could reduce the detectability of
LPnG through the scale-dependent bias, though the cancellation between the
parameters would vary with the values of $\beta_{f}$ and $\beta_{g}$, and
thereby vary with the tracer population. In this scenario, one should then
consider constraining an effective amplitude of LPnG determined by
$\mathcal{B}(z)$ in Eq. (2.12) which, if determined to be non-zero, would
already definitively prove the existence of LPnG and thereby additional light
fields during inflation – even if the individual coefficients in Eq. (1.1)
were undetermined. The challenge is that $\mathcal{B}(z)$ depends on both the
primordial non-Gaussian parameters and the tracer populations. An individual
measurement of $\mathcal{B}$ from one redshift and tracer population will have
much lower signal-to-noise than a measurement of LPnG that combines data from
all redshift bins and tracer populations. Indeed, Fisher matrix computations
for the most constraining galaxy sample of SPHEREx performed using our
likelihood pipeline indicate that a non-trivial LPnG (with
$f_{\mathrm{NL}}\sim 1$, $g_{\mathrm{NL}}\sim 10^{4}$) cannot be ruled out
with any significance by fitting for the individual $\mathcal{B}(z_{i})$ when
these parameters are free to vary independently at each redshift $z_{i}$.
Consequently, one needs a model for $\mathcal{B}(z)$ – for instance, a model
relating $\mathcal{B}(z)$ to the galaxy bias $b_{E}(z)$ or a parameterization
of the redshift dependence of $\mathcal{B}(z)$ – to be able to harness the
full power of the data to unambiguously detect LPnG at these threshold
amplitudes.

Despite the challenge of forecasting constraints on $g_{\mathrm{NL}}$ (and
other higher-order non-Gaussian parameters) with weak modelling assumptions,
one can move forward with the simpler $f_{\mathrm{NL}}$, $\tau_{\mathrm{NL}}$
scenario. In this case, we consider $g_{\mathrm{NL}}$ to be small in
comparison to $\sim 10^{4}$, so that we can neglect that contribution to
$\Delta b_{NG}$ and proceed with forecasting constraints on $f_{\mathrm{NL}}$
and $\tau_{\mathrm{NL}}$ alone with relaxed assumptions about the form of
$\beta_{f}$. To assess the impact of modelling ambiguities on our joint
forecasts of $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$, we obtained Fisher
forecasts131313Although an MCMC analysis is more precise than a Fisher matrix
analysis, we find (throughout our analysis) that the results of MCMC runs are
in good agreement with the results of Fisher forecasts. for the single LPnG
parameter $f_{\mathrm{NL}}$ as well as joint fisher forecasts for the pair of
parameters $(f_{\mathrm{NL}},\ \tau_{\mathrm{NL}})$ under a weaker modelling
assumption whereby $\beta_{f}=\delta_{c}(b_{E}-p)$ with $p$ now being a free
parameter for each galaxy sample of SPHEREx. This is still a more restrictive
choice than allowing $\beta_{f}$ (or nearly equivalently, $p$) for each galaxy
sample and redshift to be a free parameter, but may not be unrealistic and
introduces only five additional parameters (one for each galaxy sample of
SPHEREx). Note that this modelling of $\beta_{f}$ is supported by evidence
from separate universe hydrodynamical simulations [42, 43]. For simplicity, we
assumed the same fiducial values for each galaxy sample (namely $p=1$ and
$p=0.5$) and obtained joint Fisher forecasts for $f_{\mathrm{NL}}$ and
$\tau_{\mathrm{NL}}$ (around the fiducial values assumed in this paper) using
our SPHEREx likelihood with varying priors on $p$. While the fiducial values
of $p$ may in general different for the different galaxy samples, we note that
SPHEREx power spectrum constraints on $f_{\mathrm{NL}}$ are dominated by a
single galaxy sample [45]. This simple exercise should therefore help us
understand how the forecasts obtained with the modelling choices in this paper
could change under more conservative modelling assumptions.

Figure 3 shows the joint Fisher forecasts for $f_{\mathrm{NL}}$ and
$\tau_{\mathrm{NL}}$ as function of the prior uncertainty on $p$ around the
fiducial values $f_{\mathrm{NL}}=1.0,\ \tau_{\mathrm{NL}}=1.3\times 10^{2}$
(left panel of Figure 3) and $f_{\mathrm{NL}}=1.0\
\tau_{\mathrm{NL}}=1.3\times 10^{3}$ (right panel of Figure 3). Note that
these Fisher forecasts are obtained with an independent Gaussian prior on each
$p$ (i.e. one for each galaxy population of SPHEREx) whose width goes from
$\sigma(p)=0$ to $\sigma(p)=10$. From both Figure 3, we see that imposing more
and more liberal priors on $p$ (and thus making the modelling of $\beta_{f}$
more ambiguous) does indeed worsen Fisher forecasts. However, their difference
from Fisher forecasts obtained with a fixed value of $p$ (i.e. $\sigma(p)=0$)
is $\lesssim 20\%$. The point at which $\sigma_{f_{\mathrm{NL}}}$ ceases to
grow, and the value at which the degradation of constraints on
$\sigma_{f_{\mathrm{NL}}}$ saturates, is dependent on the tracer population
through fiducial values of $\beta_{f}$. In Figure 3 this is evident from the
slight differences between the shapes of curves with different values of $p$.
Since we observe a close agreement (within a few percent) between our MCMC
analysis and the corresponding Fisher forecasts (both for fixed $p$), Figure 3
shows that our MCMC results (presented in Tables 3 and Figure 2) give a good
quantitative idea of how well SPHEREx can jointly constrain $f_{\mathrm{NL}}$
and $\tau_{\mathrm{NL}}$ under more realistic (and consequently weaker)
modelling assumptions for the nuisance parameter $\beta_{f}$. Finally, we note
again that from Eq. (2.46) and Eq. (2.41), the ratio of $\tau_{\mathrm{NL}}$
to $f_{\mathrm{NL}}^{2}$ is independent of $\beta_{f}$ (or $\mathcal{B}_{c}$)
and should therefore, in principle, be measurable even with weak (or no)
assumptions about those parameters. While a detection of $\xi^{2}$ could only
be achieved if $\xi^{2}$ were very large (we estimate $\mathcal{O}(10^{3})$)
it would, independently, represent a detection of additional light degrees of
freedom during inflation.

##  5 Discussion

The simplest inflationary models that generate initial conditions with LPnG
have a primordial curvature perturbation which has non-Gaussianity
predominantly characterised by a bispectrum (parametrised by
$f_{\mathrm{NL}}$) and which is sourced wholly by the fluctuations of a light
field separate from the inflaton. Upcoming galaxy surveys such as SPHEREx show
good promise in being able to constrain such models from observations of the
galaxy power spectrum (and bispectrum) (see Ref. [28]). In this paper, we
analyse how well the upcoming SPHEREx all-sky survey can constrain LPnG models
beyond this simplest class of models. In particular, we consider a two-field
inflation model in which the primordial curvature perturbations have a non-
trivial trispectrum whose squeezed limit is parameterised by $g_{\mathrm{NL}}$
(in addition to a squeezed-limit bispectrum parametrised by
$f_{\mathrm{NL}}$). We also consider the scenario in which primordial
curvature fluctuations are sourced from fluctuations of both the inflaton and
the curvaton. In this case, the contribution of inflaton fluctuations to the
primordial curvature perturbations is Gaussian whereas the contribution of
curvaton fluctuations shows LPnG. A substantial inflaton contribution to the
primordial curvature perturbations decorrelates the collapsed limit of the
primordial trispectrum from the square of the primordial bispectrum – thus
introducing an additional parameter in the model (called $\tau_{\mathrm{NL}}$)
which parametrises the collapsed limit of the primordial trispectrum. Thus,
the non-Gaussianity in this model is characterised by three parameters namely
$f_{\mathrm{NL}}$, $g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$. In this paper,
we present joint MCMC forecasts for the LPnG parameters $f_{\mathrm{NL}},\
g_{\mathrm{NL}},\ \tau_{\mathrm{NL}}$ from galaxy power spectra obtained in
the context of the SPHEREx all-sky survey.

In Section 2.1, we show that the squeezed limits of arbitrary higher-point
correlation functions of the primordial curvature perturbations all impart the
same characteristic scale-dependence to the bias of matter density tracers –
which is an important signal targeted by large-scale galaxy surveys that
intend to constrain LPnG using galaxy power spectra. SPHEREx galaxy power
spectra could therefore constrain or detect LPnG due to arbitrarily high-order
correlation functions. In the context of the two-field inflation model assumed
in this paper, this means that we can obtain joint MCMC forecasts for the pair
of parameters $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ using the scale-
dependent galaxy bias. However, it also means that $f_{\mathrm{NL}}$ and
$g_{\mathrm{NL}}$ are strongly degenerate w.r.t. their signature on the galaxy
power spectrum, causing their joint MCMC forecasts to be substantially
degraded compared to MCMC forecasts for the individual parameters
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$. Figure 1 and Table 2 show our joint
MCMC constraints on $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ around fiducial
values compatible with state-of-the-art constraints from CMB datasets. As
expected, the joint forecasts on $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ are
degraded compared to their single-parameter forecasts (in Table 1). The strong
degeneracy between $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ is due to the fact
that their contributions to $\Delta b_{NG}$ have the same $k$-dependence and
differ only in their dependence on the Eulerian bias, $b_{E}$. This degeneracy
is reflected in their covariance, $\sim-0.9$, which also depends on modelling
choices for $\Delta b_{NG}$, the scale-dependent bias (Eqs. (2.24), (2.25),
and (3.2).

Unlike $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$, the parameter
$\tau_{\mathrm{NL}}$ parametrises the collapsed limit of a primordial
correlation function and leads to a stochasticity in the galaxy power spectrum
which increases at the largest scales (see Section 2.2). Figure 2 and Tables 3
show our joint MCMC forecasts for $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$
around fiducial values consistent with constraints obtained from CMB datasets.
In all cases, the covariance between $f_{\mathrm{NL}}$ and
$\tau_{\mathrm{NL}}$ is $\sim-0.6$. This is a reflection of the fact that
while $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ are quite degenerate in their
impact on the galaxy power spectrum, they are not as degenerate as
$f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$ owing to the fact the
$f_{\mathrm{NL}}$ terms in $\Delta b_{NG}$ scale as $1/k^{2}$ while the
$\tau_{\mathrm{NL}}$ terms scale as $1/k^{4}$. Our results also indicate that
the SPHEREx survey can potentially constrain a large, but CMB-consistent,
value of $\tau_{\mathrm{NL}}$ tightly though with degraded constraints on
$f_{\mathrm{NL}}$. As noted in Section 4.2, a robust measurement of
$\tau_{\mathrm{NL}}$ would also indicate the presence of an additional light
field during inflation even if the parameter $f_{\mathrm{NL}}$ could not be
distinguished from zero.

In Section 4.3, we note that the constraints on $f_{\mathrm{NL}}$,
$g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ obtainable from galaxy power
spectra depend on how one chooses to model the nuisance parameters $\beta_{f}$
and $\beta_{g}$ (see Eq. (2.24)). Modelling ambiguity in $\beta_{f}$ as well
as $\beta_{g}$ reduces the utility of the galaxy power spectrum in jointly
constraining $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$. This is because the
scale-dependent bias effectively probes a weighted (redshift-dependent)
combination of the parameters $f_{\mathrm{NL}}$ and $g_{\mathrm{NL}}$, where
the relative weights are set by $\beta_{f}$ and $\beta_{g}$. In a general LPnG
scenario (e.g. Eq. (1.1) to all orders in $\mathcal{Z}$), scale-dependent bias
measurements can conservatively only help us ascertain the presence/absence of
an additional light degree of freedom (over and above the inflaton); modelling
ambiguities in the bias expansion (Eq. (2.12) preclude determining which non-
trivial squeezed limits of higher-point functions are sourcing scale-dependent
bias. To obtain model-independent constraints on different inflationary models
that generate LPnG we need to exploit higher-order galaxy correlation
measurements, which provide information complementary to measurements of
galaxy power spectra.

As for $f_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$, we note in Section 4.3 that
our joint forecasts are altered by $\lesssim 20\%$ under much weaker modelling
choices than the ones we have assumed in our primary analysis. This means that
the conclusion that SPHEREx can potentially measure $\tau_{\mathrm{NL}}$ more
precisely using galaxy power spectra at the cost of degraded constraints on
$f_{\mathrm{NL}}$ likely holds true even under more realistic model choices.
Adding galaxy bispectrum data can only make the constraints better.

Since modelling ambiguities represent a big limitation in using galaxy power
spectrum data to constrain LPnG, it would be interesting to investigate
optimal ways of mitigating their impact. In particular, the fact that they are
late-time parameters and encode the non-linear physics of galaxy formation
leads to the observation that they should couple to primordial non-
Gaussianities as well as non-Gaussianities generated through gravitational
non-linearities in the same way. This approach might help suggest novel ways
to break their degeneracy with LPnG parameters. We plan to investigate along
these lines in the future.

## Acknowledgments

We acknowledge helpful conversations with Benjamin Wallisch and Massimiliano
Lattanzi. Results in this paper were obtained using the high performance
computing system at the Institute for Advanced Computational Science at Stony
Brook University. CS is grateful for support from the NASA ATP Award
80NSSC20K0541. ML is supported by the Department of Physics and the College of
Arts and Sciences at the University of Washington and by the Department of
Energy grant DE-SC0023183. TB was supported by ICSC – Centro Nazionale di
Ricerca in High Performance 18 Computing, Big Data and Quantum Computing,
funded by European Union – NextGenerationEU.

## References

  * [1] P.D. Meerburg et al., _Primordial Non-Gaussianity_ , _Bull. Am. Astron. Soc._ 51 (2019) 107 [1903.04409].
  * [2] A. Slosar et al., _Scratches from the Past: Inflationary Archaeology through Features in the Power Spectrum of Primordial Fluctuations_ , _Bull. Am. Astron. Soc._ 51 (2019) 98 [1903.09883].
  * [3] A.H. Guth and S.Y. Pi, _Fluctuations in the New Inflationary Universe_ , _Phys. Rev. Lett._ 49 (1982) 1110.
  * [4] A.A. Starobinsky, _Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations_ , _Phys. Lett. B_ 117 (1982) 175.
  * [5] J.M. Bardeen, P.J. Steinhardt and M.S. Turner, _Spontaneous Creation of Almost Scale - Free Density Perturbations in an Inflationary Universe_ , _Phys. Rev. D_ 28 (1983) 679.
  * [6] J.M. Maldacena, _Non-Gaussian features of primordial fluctuations in single field inflationary models_ , _JHEP_ 05 (2003) 013 [astro-ph/0210603].
  * [7] X. Chen, M.-x. Huang, S. Kachru and G. Shiu, _Observational signatures and non-Gaussianities of general single field inflation_ , _JCAP_ 01 (2007) 002 [hep-th/0605045].
  * [8] C. Cheung, P. Creminelli, A.L. Fitzpatrick, J. Kaplan and L. Senatore, _The Effective Field Theory of Inflation_ , _JHEP_ 03 (2008) 014 [0709.0293].
  * [9] W.H. Kinney, _Horizon crossing and inflation with large eta_ , _Phys. Rev. D_ 72 (2005) 023515 [gr-qc/0503017].
  * [10] J. Martin, H. Motohashi and T. Suyama, _Ultra Slow-Roll Inflation and the non-Gaussianity Consistency Relation_ , _Phys. Rev. D_ 87 (2013) 023514 [1211.0083].
  * [11] K. Enqvist and M.S. Sloth, _Adiabatic CMB perturbations in pre - big bang string cosmology_ , _Nucl. Phys. B_ 626 (2002) 395 [hep-ph/0109214].
  * [12] D.H. Lyth and D. Wands, _Generating the curvature perturbation without an inflaton_ , _Phys. Lett. B_ 524 (2002) 5 [hep-ph/0110002].
  * [13] M. Sasaki, J. Valiviita and D. Wands, _Non-Gaussianity of the primordial perturbation in the curvaton model_ , _Phys. Rev. D_ 74 (2006) 103003 [astro-ph/0607627].
  * [14] G. Dvali, A. Gruzinov and M. Zaldarriaga, _A new mechanism for generating density perturbations from inflation_ , _Phys. Rev. D_ 69 (2004) 023505 [astro-ph/0303591].
  * [15] G. Dvali, A. Gruzinov and M. Zaldarriaga, _Cosmological perturbations from inhomogeneous reheating, freezeout, and mass domination_ , _Phys. Rev. D_ 69 (2004) 083505 [astro-ph/0305548].
  * [16] P. Creminelli and M. Zaldarriaga, _Single field consistency relation for the 3-point function_ , _JCAP_ 10 (2004) 006 [astro-ph/0407059].
  * [17] P. Creminelli, G. D’Amico, M. Musso and J. Norena, _The (not so) squeezed limit of the primordial 3-point function_ , _JCAP_ 11 (2011) 038 [1106.1462].
  * [18] Planck collaboration, _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9 [1905.05697].
  * [19] S. Shandera, N. Dalal and D. Huterer, _A generalized local ansatz and its effect on halo bias_ , _JCAP_ 03 (2011) 017 [1010.3722].
  * [20] D. Tseliakhovich, C. Hirata and A. Slosar, _Non-Gaussianity and large-scale structure in a two-field inflationary model_ , _Phys. Rev. D_ 82 (2010) 043531 [1004.3302].
  * [21] T. Suyama and M. Yamaguchi, _Non-gaussianity in the modulated reheating scenario_ , _Phys. Rev. D_ 77 (2008) 023505.
  * [22] K.M. Smith, M. LoVerde and M. Zaldarriaga, _A universal bound on N-point correlations from inflation_ , _Phys. Rev. Lett._ 107 (2011) 191301 [1108.1805].
  * [23] N. Dalal, O. Dore, D. Huterer and A. Shirokov, _The imprints of primordial non-gaussianities on large-scale structure: scale dependent bias and abundance of virialized objects_ , _Phys. Rev. D_ 77 (2008) 123514 [0710.4560].
  * [24] A. Slosar, C. Hirata, U. Seljak, S. Ho and N. Padmanabhan, _Constraints on local primordial non-Gaussianity from large scale structure_ , _JCAP_ 08 (2008) 031 [0805.3580].
  * [25] N. Agarwal, S. Ho and S. Shandera, _Constraining the initial conditions of the Universe using large scale structure_ , _JCAP_ 02 (2014) 038 [1311.2606].
  * [26] T. Giannantonio, A.J. Ross, W.J. Percival, R. Crittenden, D. Bacher, M. Kilbinger et al., _Improved Primordial Non-Gaussianity Constraints from Measurements of Galaxy Clustering and the Integrated Sachs-Wolfe Effect_ , _Phys. Rev. D_ 89 (2014) 023511 [1303.1349].
  * [27] B. Leistedt, H.V. Peiris and N. Roth, _Constraints on Primordial Non-Gaussianity from 800 000 Photometric Quasars_ , _Phys. Rev. Lett._ 113 (2014) 221301 [1405.4315].
  * [28] O. Doré et al., _Cosmology with the SPHEREX All-Sky Spectral Survey_ , 1412.4872.
  * [29] S. Ferraro et al., _Inflation and Dark Energy from Spectroscopy at $z>2$_, _Bull. Am. Astron. Soc._ 51 (2019) 72 [1903.09208].
  * [30] G. Cabass, M.M. Ivanov, O.H.E. Philcox, M. Simonović and M. Zaldarriaga, _Constraints on multifield inflation from the BOSS galaxy survey_ , _Phys. Rev. D_ 106 (2022) 043506 [2204.01781].
  * [31] D.J. Schlegel et al., _The MegaMapper: A Stage-5 Spectroscopic Instrument Concept for the Study of Inflation and Dark Energy_ , 2209.04322.
  * [32] D. Baumann, S. Ferraro, D. Green and K.M. Smith, _Stochastic bias from non-gaussian initial conditions_ , _Journal of Cosmology and Astroparticle Physics_ 2013 (2013) 001.
  * [33] K.M. Smith, S. Ferraro and M. LoVerde, _Halo clustering and gnl-type primordial non-gaussianity_ , _Journal of Cosmology and Astroparticle Physics_ 2012 (2012) 032.
  * [34] K.M. Smith and M. LoVerde, _Local stochastic non-Gaussianity and N-body simulations_ , _JCAP_ 11 (2011) 009 [1010.0055].
  * [35] D. Jeong, F. Schmidt and C.M. Hirata, _Large-scale clustering of galaxies in general relativity_ , _Phys. Rev. D_ 85 (2012) 023504 [1107.5427].
  * [36] P. Creminelli, J. Gleyzes, M. Simonović and F. Vernizzi, _Single-Field Consistency Relations of Large Scale Structure. Part II: Resummation and Redshift Space_ , _JCAP_ 02 (2014) 051 [1311.0290].
  * [37] L. Senatore and M. Zaldarriaga, _A Note on the Consistency Condition of Primordial Fluctuations_ , _JCAP_ 08 (2012) 001 [1203.6884].
  * [38] J.M. Bardeen, J.R. Bond, N. Kaiser and A.S. Szalay, _The Statistics of Peaks of Gaussian Random Fields_ , _Astrophys. J._ 304 (1986) 15.
  * [39] T. Tanaka and Y. Urakawa, _Dominance of gauge artifact in the consistency relation for the primordial bispectrum_ , _JCAP_ 05 (2011) 014 [1103.1251].
  * [40] E. Pajer, F. Schmidt and M. Zaldarriaga, _The Observed Squeezed Limit of Cosmological Three-Point Functions_ , _Phys. Rev. D_ 88 (2013) 083502 [1305.0824].
  * [41] R. de Putter, O. Doré and D. Green, _Is There Scale-Dependent Bias in Single-Field Inflation?_ , _JCAP_ 10 (2015) 024 [1504.05935].
  * [42] A. Barreira, G. Cabass, F. Schmidt, A. Pillepich and D. Nelson, _Galaxy bias and primordial non-gaussianity: insights from galaxy formation simulations with illustristng_ , _Journal of Cosmology and Astroparticle Physics_ 2020 (2020) 013.
  * [43] A. Barreira, _Predictions for local png bias in the galaxy power spectrum and bispectrum and the consequences for fnl constraints_ , _Journal of Cosmology and Astroparticle Physics_ 2022 (2022) 033.
  * [44] P. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela et al., _Planck 2013 results. xxiv. constraints on primordial non-gaussianity_ , _ArXiv e-prints_ (2013) .
  * [45] C. Shiveshwarkar, T. Brinckmann, M. Loverde and M. McQuinn, _Postinflationary contamination of local primordial non-Gaussianity in galaxy power spectra_ , _Phys. Rev. D_ 108 (2023) 103538 [2306.07517].
  * [46] B. Audren, J. Lesgourgues, K. Benabed and S. Prunet, _Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code_ , _JCAP_ 02 (2013) 001 [1210.7183].
  * [47] T. Brinckmann and J. Lesgourgues, _MontePython 3: boosted MCMC sampler and other features_ , _Phys. Dark Univ._ 24 (2019) 100260 [1804.07261].
  * [48] E.L. Wright et al., _The Wide-field Infrared Survey Explorer (WISE): Mission Description and Initial On-orbit Performance_ , _Astron. J._ 140 (2010) 1868 [1008.0031].
  * [49] N. Metcalfe et al., _The Pan-STARRS1 Small Area Survey 2_ , _Mon. Not. Roy. Astron. Soc._ 435 (2013) 1825 [1310.6368].
  * [50] DES collaboration, _The Dark Energy Survey_ , astro-ph/0510346.
  * [51] D. Blas, J. Lesgourgues and T. Tram, _The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes_ , _JCAP_ 07 (2011) 034 [1104.2933].
  * [52] J. Lesgourgues, _The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview_ , 1104.2932.
  * [53] J. Lesgourgues and T. Tram, _The Cosmic Linear Anisotropy Solving System (CLASS) IV: efficient implementation of non-cold relics_ , _JCAP_ 09 (2011) 032 [1104.2935].
  * [54] N. Kaiser, _Clustering in real space and in redshift space_ , _Mon. Not. Roy. Astron. Soc._ 227 (1987) 1.
  * [55] T. Sprenger, M. Archidiacono, T. Brinckmann, S. Clesse and J. Lesgourgues, _Cosmology in the era of Euclid and the Square Kilometre Array_ , _JCAP_ 02 (2019) 047 [1801.08331].
  * [56] T. Brinckmann, D.C. Hooper, M. Archidiacono, J. Lesgourgues and T. Sprenger, _The promising future of a robust cosmological neutrino mass measurement_ , _JCAP_ 01 (2019) 059 [1808.05955].
  * [57] M. Gerbino, K. Freese, S. Vagnozzi, M. Lattanzi, O. Mena, E. Giusarma et al., _Impact of neutrino properties on the estimation of inflationary parameters from current and future observations_ , _Phys. Rev. D_ 95 (2017) 043512 [1610.08830].
  * [58] S. Ferraro and K.M. Smith, _Using large scale structure to measure $f_{NL},g_{NL}$ and $\tau_{NL}$_, _Phys. Rev. D_ 91 (2015) 043506 [1408.3126].