††thanks: Einstein Fellow

# Constraints on Multi-Field Inflation from the BOSS Galaxy Survey

Giovanni Cabass<EMAIL_ADDRESS>School of Natural Sciences, Institute for
Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA  Mikhail M. Ivanov
<EMAIL_ADDRESS>School of Natural Sciences, Institute for Advanced Study, 1
Einstein Drive, Princeton, NJ 08540, USA  Oliver H. E. Philcox  Department of
Astrophysical Sciences, Princeton University, Princeton, NJ 08540, USA  School
of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive,
Princeton, NJ 08540, USA  Marko Simonović  Theoretical Physics Department,
CERN, 1 Esplanade des Particules, Geneva 23, CH-1211, Switzerland  Matias
Zaldarriaga  School of Natural Sciences, Institute for Advanced Study, 1
Einstein Drive, Princeton, NJ 08540, USA

###### Abstract

We use redshift-space galaxy clustering data from the BOSS survey to constrain
local primordial non-Gaussianity (LPNG). This is of particular importance due
to the consistency relations, which imply that a detection of LPNG would rule
out all single-field inflationary models. Our constraints are based on the
consistently analyzed redshift-space galaxy power spectra and bispectra,
extracted from the public BOSS data with optimal window-free estimators. We
use a complete perturbation theory model including all one-loop power spectrum
corrections generated by LPNG. Our constraint on the amplitude of the local
non-Gaussian shape is $f_{\rm NL}^{\rm local}=-33\pm 28$ at 68% CL, yielding
no evidence for primordial non-Gaussianity. The addition of the bispectrum
tightens the $f_{\rm NL}^{\rm local}$ constraints from BOSS by $20\%$, and
allows breaking of degeneracies with non-Gaussian galaxy bias. These results
set the stage for the analysis of future surveys, whose larger volumes will
yield significantly tighter constraints on LPNG.

††preprint: CERN-TH-2022-055

##  1 Introduction

Inflation provides a mechanism to seed density fluctuations that we observe in
the late Universe. The physics responsible for it, which may have operated at
energies as high as $10^{16}\,{\rm GeV}$, has left observable imprints in
these density fluctuations. The observations of cosmic microwave background
(CMB) anisotropies and the distribution of galaxies in the large-scale
structure (LSS) present particularly appealing ways to probe the inflationary
epoch, and thus the physics of this high-energy regime.

There is a special class of inflationary models, in which inflation is driven
by a medium whose quantum fluctuations are the only source of the observable
overdensities. Assuming the attractor solution and the Bunch-Davies vacuum,
these single-field (or single-clock) models generically predict purely
adiabatic fluctuations with vanishing physical coupling of long-wavelength and
short-wavelength modes. This result is formalized in the well-known
consistency relations Maldacena (2003); Creminelli and Zaldarriaga (2004).
Given the bispectrum $B_{\phi}(k_{1},k_{2},k_{3})$ of the primordial Bardeen
potential $\phi$, they have the following form

$B_{\phi}(k_{1},k_{2},k_{3})\Big{|}_{k_{3}\ll
k_{1},k_{2}}=-{\frac{5}{3}}P_{\phi}(k_{3})\left[3+k_{1}\frac{\partial}{\partial
k_{1}}\right]P_{\phi}(k_{1})\;.$ (1)

This limit, when one of the wavenumbers is much smaller than the other two, is
called the squeezed limit. Eq. (1) implies that in single-field models the
only effect of the long-wavelength modes of $\phi$ on the short-scale modes is
a simple rescaling of coordinates, which is locally unobservable. Therefore,
any detection of local primordial non-Gaussianity (LPNG), i.e. a detection of
a non-vanishing amplitude in the squeezed limit of the initial bispectrum due
to physical interactions of long and short modes, would rule out single-field
inflation Creminelli and Zaldarriaga (2004).

Whilst the simplicity of single-field inflation is very appealing (and so far
supported by observations), having more than one relevant fluctuating degree
of freedom besides the inflaton is easy to achieve. Some well-known examples
are the curvaton scenario Enqvist and Sloth (2002); Lyth and Wands (2002);
Moroi and Takahashi (2001) and modulated reheating Zaldarriaga (2004). In
contrast to single-field inflation, these multi-field models can produce large
and observable LPNG.

This important distinction between the two classes of inflationary models
makes the amplitude of the initial bispectrum containing LPNG, called $f_{\rm
NL}^{\rm local}$, the key observable that we can use to make quantitative,
model-independent statements about the primordial Universe. It is defined in
the following way

$B_{\phi}(k_{1},k_{2},k_{3})=6f_{\rm NL}^{\rm
local}\Delta^{4}_{\phi}\frac{{\cal S}_{\rm
local}(k_{1},k_{2},k_{3})}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}\,\,,$ (2)

where $\Delta^{2}_{\phi}$ is the amplitude of the primordial power spectrum
$k^{3}P_{\phi}(k)=\Delta^{2}_{\phi}(k/k_{\ast})^{n_{s}-1}$ and $n_{s}$ is the
spectral index.111 _Planck_ 2018 Aghanim  _et al._ (2018) gives precise
measurements for both of these quantities: $\Delta^{2}_{\phi}\approx 1.5\times
10^{-8}$ and $n_{s}\approx 0.96$, for the pivot scale $k_{\ast}=0.05$ Mpc-1.
The local template is given by

${\cal S}_{\rm
local}(k_{1},k_{2},k_{3})=\frac{1}{3}\frac{k_{1}^{2}}{k_{2}k_{3}}+2~{}{\rm
perms.}\,\,.$ (3)

This template is such that the squeezed limit bispectrum generated by the LPNG
takes the following form

$B_{\phi}(k_{1},k_{2},k_{3})\Big{|}_{k_{3}\ll k_{1},k_{2}}=4f_{\rm NL}^{\rm
local}P_{\phi}(k_{1})P_{\phi}(k_{3})\,\,,$ (4)

which is very different from the single-field result given by Eq. (1). Generic
values of $f_{\rm NL}^{\rm local}$ in multifield models are of order one or
higher (for some counterexamples, see Vernizzi and Wands (2006); Senatore and
Zaldarriaga (2012)), making $f_{\rm NL}^{\rm local}\approx 1$ an interesting
and well-motivated observational target.

Significant efforts are aimed at measuring $f_{\rm NL}^{\rm local}$, with the
tightest current constraints coming from CMB observations. In particular, the
_Planck_ 2018 data yields $f_{\rm NL}^{\rm local}=-0.9\pm 5.1$ Akrami  _et
al._ (2020). Measurements from galaxy clustering are currently somewhat
weaker, but are expected to improve significantly with upcoming galaxy
surveys. These surveys will eventually reach the target of $f_{\rm NL}^{\rm
local}\approx 1$ (see for example Doré  _et al._ (2014); Ferraro  _et al._
(2019); Castorina  _et al._ (2020)). Almost all LSS analyses done so far use
the fact that LPNG produces the so-called scale-dependent galaxy bias Dalal
_et al._ (2008); Matarrese and Verde (2008), and therefore can be constrained
by observations of galaxy power spectra on large scales Slosar  _et al._
(2008); Xia  _et al._ (2010); Ross  _et al._ (2013); Castorina  _et al._
(2019); Mueller  _et al._ (2021). Whilst measuring the galaxy power spectrum
and looking for scale-dependent bias has the advantage of being relatively
straightforward, this may not be an optimal way to constrain LPNG from LSS
data. Indeed, as many Fisher forecasts and full likelihood analyses indicate,
the dominant source of information on $f_{\rm NL}^{\rm local}$ for the shot-
noise limited samples is the galaxy bispectrum Scoccimarro  _et al._ (2004);
Baldauf  _et al._ (2011); Doré  _et al._ (2014); Ferraro  _et al._ (2019);
Moradinezhad Dizgah  _et al._ (2021); Barreira (2020). Developing consistent
and robust pipelines to harvest this information is one of the major
milestones on the way towards achieving the tightest possible bounds on LPNG.

Performing an optimal search for LPNG in the observed galaxy bispectrum is not
a trivial task for a number of reasons. One of the main difficulties is the
survey geometry, which mixes the Fourier modes on large scales. In order to
circumvent this problem, in this paper we use recently developed optimal
window-free power spectrum and bispectrum estimators Philcox (2021a, b). In
principle, such an approach guarantees that the results are unbiased, close-
to-optimal and that all effects related to the window function are
consistently taken into account. This is particularly important for
constraints on $f_{\rm NL}^{\rm local}$, since most of the signal comes from
the largest scales in the survey, either through the scale-dependent bias or
through the squeezed triangles. An alternative is to model the effects of
window convolution when calculating the observed bispectrum. Doing this
exactly is very challenging numerically, and novel methods to tackle this
problem were developed only very recently Pardede  _et al._ (2022). On the
other hand, if the effects of the window function are modelled using
approximate treatments available in the literature (e.g., Gil-Marín  _et al._
, 2015, 2017), this can lead to biases at low $k$, leading to such bins
needing to be dropped from the analysis. This approach was used recently in
D’Amico  _et al._ (2022) to measure LPNG, and it remains unclear to what
extent the results are impacted by the approximate treatment of the window
function.

Another non-trivial task is the modeling of the galaxy bispectrum signal. This
includes all aspects of the nonlinear evolution such as the backreaction of
short-scale physics on large-scale modes, the nonlinear evolution of the BAO
signature (IR resummation), as well as a robust control over projection and
binning effects. Many years of intense theoretical efforts Scoccimarro  _et
al._ (1999); Scoccimarro (2000); Sefusatti  _et al._ (2006); Sefusatti and
Komatsu (2007); Sefusatti (2009); Baldauf  _et al._ (2015); Angulo  _et al._
(2015); Blas  _et al._ (2016a); Nadler  _et al._ (2018); Eggemeier  _et al._
(2018, 2021); Desjacques  _et al._ (2018a, b); Ivanov and Sibiryakov (2018);
Oddo  _et al._ (2020, 2021); Ivanov  _et al._ (2021a) have recently made the
incorporation of these effects possible, so that the bispectrum data can be
routinely used in cosmological parameter analyses Philcox and Ivanov (2021).

In this work, we present a search of LPNG using galaxy power spectrum and
bispectrum from the publicly available BOSS data Alam  _et al._ (2017). This
is a natural continuation of our previous analysis of PNG in single-field
inflation Cabass  _et al._ (2022a), based on the same tree-level bispectrum
model and the data cuts that were extensively tested in Ivanov  _et al._
(2021a). Importantly, the galaxy bispectrum treatment presented in Ivanov  _et
al._ (2021a) is fully systematic, i.e. there is a way to control the precision
of various effects so it can be applied to future high-precision galaxy survey
data.

The remainder of this paper is structured as follows. In Section 2 we present
key theoretical ingredients needed to extract LPNG from the galaxy clustering
data. They include the complete calculation of the one-loop galaxy power
spectrum in the presence of LPNG. Section 3 describes the data and analysis
details. We validate our pipeline on the mock galaxy clustering data in
Section 4, and then apply it to the BOSS data in Section 5. We present limits
on non-Gaussian bias parameters from the BOSS data in Sec. 6. Section 7 draws
conclusions. Additional details of the theory model are given in the Appendix
A.

##  2 Structure Formation in the presence of LPNG

In this section we present our theoretical model, which includes all necessary
terms generated by LPNG. We work in the framework of the effective field
theory of large scale structure (hereafter EFT of LSS), as described in Perko
_et al._ (2016); Ivanov  _et al._ (2020a); D’Amico  _et al._ (2019); Ivanov
_et al._ (2020b); Nishimichi  _et al._ (2020); Cabass  _et al._ (2022b) and
references therein. Since the perturbative model for structure formation has
been discussed in detail in the works cited above, we will provide only a
brief overview in what follows, focusing on the ingredients that will be
necessary to carry out the calculation of the one-loop LPNG contributions. For
dark matter and biased tracers in real space, these contributions have already
been studied in the literature, Taruya  _et al._ (2008); Assassi  _et al._
(2015a, b); Moradinezhad Dizgah  _et al._ (2021). Here, we extend these
results to the case of galaxies in redshift-space needed for comparison to
observations.

To simplify the text and formulas, we omit the superscript “local” in our
notation of the $f_{\rm NL}$ parameter in this section. We stress, however,
that all results obtained here, strictly speaking, apply only to the case of
LPNG. The analogous theory model for non-local primordial non-Gaussianity is
presented in Ref. (Cabass  _et al._ , 2022a).

###  2.1 Gaussian part

For Gaussian initial conditions, the statistical properties of the first-order
density field, $\delta^{(1)}$, are completely determined by its power spectrum
$P_{11}$:

$\langle\delta^{(1)}({\bf k})\delta^{(1)}({\bf
k}^{\prime})\rangle=(2\pi)^{3}\delta^{(3)}_{D}({\bf k}^{\prime}+{\bf
k})P_{11}(k)\,,$ (5)

where we have suppressed the explicit time dependence for brevity. In this
work, we restrict our analysis of the galaxy power spectrum to one-loop order
in the EFT of LSS, where the usual Gaussian part reads

$P_{\text{Gauss}}=P_{\rm tree}+P_{\rm 1-loop}+P_{\rm ctr}+P_{\rm stoch}\,\,,$
(6)

where $P_{\rm tree}$ is the linear term, $P_{\rm 1-loop}$ is the one-loop
correction (with Gaussian initial conditions), $P_{\rm ctr}$ is the higher
derivative term (counterterm), and $P_{\rm stoch}$ is the term that captures
galaxy stochasticity.

We now present our model for the galaxy bias, referring the interested reader
to Desjacques  _et al._ (2018a) for an extensive review. For the statistics
considered herein, it is sufficient to consider the galaxy density field at
cubic order. We use the following set of bias operators

$\begin{split}\delta_{g}^{(r)}=b_{1}\delta+\frac{b_{2}}{2}\delta^{2}+b_{{\cal
G}_{2}}{\cal G}_{2}+\frac{b_{3}}{6}\delta^{3}+b_{\delta{\cal
G}_{2}}\delta{\cal G}_{2}+b_{{\cal G}_{3}}{\cal
G}_{3}+b_{\Gamma_{3}}\Gamma_{3}+R^{2}_{\ast}\partial^{2}\delta\,\,.\end{split}$
(7)

The Galileon operator ${\cal G}_{2}$ is defined as
$(\partial_{i}\partial_{j}\Phi_{g})^{2}-(\partial^{2}\Phi_{g})^{2}$, where
$\Phi_{g}$ is the gravitational potential. The $\Gamma_{3}$ operator instead
is defined as ${\cal G}_{2}[\Phi_{g}]-{\cal G}_{2}[\Phi_{v}]$, where
$\Phi_{v}$ is the velocity potential. The cubic operators $\delta^{3}$,
$\delta{\cal G}_{2}$, ${\cal G}_{3}$ do not contribute to the one-loop power
spectrum after renormalization.

The redshift-space mapping for fixed line of sight $\hat{\textbf{z}}$ at order
$(\delta^{(1)})^{3}$ is given by

$\delta_{g}^{(s)}=\delta_{g}^{(r)}-\partial_{z}\big{(}u_{z}(1+\delta^{(r)}_{g})\big{)}+\frac{1}{2}\partial^{2}_{z}\big{(}u^{2}_{z}(1+\delta^{(r)}_{g})\big{)}-\frac{1}{6}\partial^{3}_{z}\big{(}u^{3}_{z}\big{)}\,\,,$
(8)

(from expanding the usual real-space to redshift-space relation), where
$u_{z}\equiv\hat{\textbf{z}}\cdot{\bf v}/{\cal H}$, ${\bf v}$ is the peculiar
velocity field and ${\cal H}$ is the conformal Hubble parameter. These can be
written in terms of the Fourier-space kernels

$\displaystyle Z_{1}({\bf k})$ $\displaystyle=$ $\displaystyle
b_{1}+f\mu^{2}\,,$ (9) $\displaystyle Z_{2}({\bf k}_{1},{\bf k}_{2})$
$\displaystyle=$
$\displaystyle\frac{b_{2}}{2}+b_{\mathcal{G}_{2}}\left(\frac{({\bf
k}_{1}\cdot{\bf k}_{2})^{2}}{k_{1}^{2}k_{2}^{2}}-1\right)+b_{1}F_{2}({\bf
k}_{1},{\bf k}_{2})+f\mu^{2}G_{2}({\bf k}_{1},{\bf k}_{2})+\frac{f\mu
k_{12}}{2}\left(\frac{\mu_{1}}{k_{1}}(b_{1}+f\mu^{2}_{2})+\frac{\mu_{2}}{k_{2}}(b_{1}+f\mu^{2}_{1})\right)\,,$

and $Z_{3}$, whose expression can be found in Section 2 of Chudaykin  _et al._
(2020), where $f$ is the logarithmic growth rate and $\mu=\hat{\bf
k}\cdot\hat{\textbf{z}}$. In this notation, the deterministic part of the
redshift-space galaxy density field can be written as

$\delta_{g}({\bf k})=Z_{1}({\bf k})\delta^{(1)}({\bf k})+\,\int_{{\bf
p}_{12}={\bf k}}Z_{2}({\bf p}_{1},{\bf p}_{2})\delta^{(1)}({\bf
p}_{1})\delta^{(1)}({\bf p}_{2})+\int_{{\bf p}_{123}={\bf k}}Z_{3}({\bf
p}_{1},{\bf p}_{2},{\bf p}_{3})\delta^{(1)}({\bf p}_{1})\delta^{(1)}({\bf
p}_{2})\delta^{(1)}({\bf p}_{3})\,\,,$ (10)

where we have introduced the following notation $\int_{{\bf p}_{1\ldots
n}={\bf k}}\equiv\int\frac{d^{3}{\bf p}_{1}}{(2\pi)^{3}}\cdots\frac{d^{3}{\bf
p}_{n}}{(2\pi)^{3}}(2\pi)^{3}\delta^{D}({\bf k}-{\bf p}_{1\ldots n})$ and
${\bf p}_{1\ldots n}\equiv{\bf p}_{1}+\cdots+{\bf p}_{n}$. We also supplement
these kernels with the appropriate redshift-space counterterms that are
omitted in (8) for clarity; these are discussed in Ref. (Ivanov  _et al._ ,
2021a).

###  2.2 LPNG-related non-linearity

LPNG affects the statistics of the galaxy overdensity in two ways. First, we
have a non-zero bispectrum for the linear matter overdensity $\delta^{(1)}$.
This generates a connected bispectrum contribution $B_{111}$,

$\langle\delta^{(1)}({\bf k}_{1})\delta^{(1)}({\bf k}_{2})\delta^{(1)}({\bf
k}_{3})\rangle\equiv(2\pi)^{3}\delta_{D}^{(3)}({\bf
k}_{123})B_{111}(k_{1},k_{2},k_{3})=(2\pi)^{3}\delta_{D}^{(3)}({\bf
k}_{123})\prod_{a=1}^{3}\mathcal{M}(k_{a})B_{\phi}(k_{1},k_{2},k_{3})\,,$ (11)

where we have introduced the transfer functions

$\delta^{(1)}({\bf k})={\cal M}(k)\phi({\bf k})\quad\Rightarrow\quad{\cal
M}(k)=\sqrt{\frac{P_{11}(k)}{P_{\phi}(k)}}\,\,.$ (12)

The initial bispectrum (11) also generates an additional loop correction to
the matter power spectrum dubbed $P_{12}$. We will discuss this term shortly.

Second, LPNG modulates the correlation between the long and short modes, which
ultimately alters the probability of galaxy formation (inducing scale-
dependent bias Dalal  _et al._ (2008); Matarrese and Verde (2008)). In order
to reproduce this effect in the perturbative galaxy bias expansion, one needs
to include new operators with the appropriate bias coefficients analogous to
the Gaussian case (7). At linear order in $f_{\rm NL}\Delta_{\phi}$ and cubic
order in $\delta^{(1)}$ these operators are given by

$\begin{split}&\delta_{g}^{\rm LPNG}(\textbf{x})=b_{\phi}f_{\rm NL}\phi({\bf
q})+b_{\phi\delta}f_{\rm NL}\phi({\bf
q})\delta(\textbf{x})+b_{\phi\delta^{2}}f_{\rm NL}\phi({\bf
q})\delta^{2}(\textbf{x})+b_{\phi{\cal G}_{2}}f_{\rm NL}\phi({\bf q}){\cal
G}_{2}(\textbf{x})\,\,.\end{split}$ (13)

Note that this expansion is valid only for LPNG. For non-local primordial non-
Gaussianity the squeezed bispectrum is typically proportional to derivatives
of $\phi$, and hence $\phi$ in the above expansion must be replaced by
appropriate higher derivative operators like $\partial^{2}\phi$ Schmidt  _et
al._ (2019). These operators appear to be higher order and hence their effect
can be neglected at the one-loop order in the EFT of LSS Cabass  _et al._
(2022a).

In contrast to Eq. (7), here we have made the argument of all relevant fields
explicit. More precisely, the Bardeen potential appearing on the right-hand
side is evaluated at the Lagrangian position ${\bf q}$ corresponding to the
Eulerian position x Assassi  _et al._ (2015a, b); Moradinezhad Dizgah  _et
al._ (2021). In order to evaluate all fields at the Eulerian coordinates we
need to Taylor expand the primordial gravitational potential. If we want to
keep all terms up to cubic order we can write

$\phi({\bf q})=\phi(\textbf{x}-{\boldsymbol{\psi}}({\bf
q}))=\phi\big{(}\textbf{x}-{\boldsymbol{\psi}}(\textbf{x}-{\boldsymbol{\psi}}(\textbf{x}))\big{)}\,\,.$
(14)

Expanding perturbatively in the displacement field ${\boldsymbol{\psi}}$ we
get

$\phi({\bf
q})=\phi-{\bf\psi}^{i}\partial_{i}\phi+\psi^{k}(\partial_{k}\psi^{i})\partial_{i}\phi+\frac{1}{2}\psi^{i}\psi^{j}\partial_{i}\partial_{j}\phi\,\,,$
(15)

where the fields on the right-hand side are all evaluated at the Eulerian
position x, and we emphasize that the displacement ${\bf\psi}$ contains both
the linear and the second-order contribution. We keep terms up to cubic order
in the expansion (15) since they are needed for the consistent calculation of
the one-loop power spectrum. Before we move on, let us comment on the omission
of higher derivative terms of the form $\partial^{2}_{{\bf q}}\phi({\bf q})$
in (13). These corrections can be straightforwardly included, see e.g. Assassi
_et al._ (2015a), but for realistic values of $f_{\rm NL}$ they are always
suppressed compared to the two-loop Gaussian contributions that we neglect
here. Therefore, we neglect the higher derivative LPNG terms in what follows.

Let us now shift our attention to redshift-space. In this case LPNG generates
additional counterterms in $\delta_{g}$ involving the matter velocity field
${\bf v}$. However, as we have just discussed, these terms can be neglected in
our analysis because they have the same order of magnitude as the higher
derivative LPNG operators. Hence, it is enough to use (8) to map the rest-
frame galaxy overdensity in presence of LPNG to redshift-space.

All in all, the Taylor expansion of $\delta_{\rm NL}$ though $\delta^{(1)}$ in
the presence of LPNG will take a form identical to Eq. (10), but with the new
kernels $Z^{\rm tot}_{n}=Z_{n}+Z^{\rm NG}_{n}$ ($n=1,2,3$), where $Z^{\rm
NG}_{n}$ are the additional PNG kernel contributions. The linear kernel is
given by

$Z^{\rm NG}_{1}(k)=b_{\phi}f_{\rm NL}\,,$ (16)

with the second kernel taking the form

$\begin{split}Z^{\rm NG}_{2}({\bf p}_{1},{\bf p}_{2})&=b_{\phi}f_{\rm
NL}\,\frac{{\bf p}_{1}\cdot{\bf
p}_{2}}{2p_{1}p_{2}}\bigg{(}\frac{p_{2}}{p_{1}}\frac{1}{{\cal
M}(p_{2})}+\frac{p_{1}}{p_{2}}\frac{1}{{\cal M}(p_{1})}\bigg{)}+b_{\phi}f_{\rm
NL}\,\frac{f\mu k}{2}\bigg{(}\frac{\mu_{1}}{p_{1}}\frac{1}{{\cal
M}(p_{2})}+\frac{\mu_{2}}{p_{2}}\frac{1}{{\cal M}(p_{1})}\bigg{)}\\\
&\;\;\;\;+b_{\phi\delta}f_{\rm NL}\,\frac{1}{2}\bigg{(}\frac{1}{{\cal
M}(p_{1})}+\frac{1}{{\cal M}(p_{2})}\bigg{)}\,\,,\end{split}$ (17)

where we have introduced

$\mu_{i}=\hat{\textbf{z}}\cdot\hat{\bf
p}_{i}\,\,,\quad\mu_{ij}=\hat{\textbf{z}}\cdot({\bf p}_{i}+{\bf
p}_{j})/\lvert{\bf p}_{i}+{\bf p}_{j}\rvert\,\,.$ (18)

For the cubic fields we find

$\begin{split}&Z^{\rm NG}_{3}({\bf p}_{1},{\bf p}_{2},{\bf
p}_{3})=b_{\phi}f_{\rm NL}\,\bigg{(}{-\frac{1}{14}}{\cal G}_{2}({\bf
p}_{1},{\bf p}_{2})\frac{({\bf p}_{1}+{\bf p}_{2})\cdot{\bf p}_{3}}{\lvert{\bf
p}_{1}+{\bf p}_{2}\rvert^{2}}\frac{1}{{\cal M}(p_{3})}+\text{$2$
perms.}\bigg{)}\\\ &\;\;\;\;+b_{\phi}f_{\rm NL}\,\bigg{(}\frac{1}{6}\frac{{\bf
p}_{1}\cdot{\bf p}_{2}}{p^{2}_{1}p^{2}_{2}}\frac{{\bf p}_{2}\cdot{\bf
p}_{3}}{{\cal M}(p_{3})}+\text{$5$ perms.}\bigg{)}+b_{\phi}f_{\rm
NL}\,\bigg{(}\frac{1}{6}\frac{{\bf p}_{1}\cdot{\bf
p}_{3}}{p^{2}_{1}p^{2}_{2}}\frac{{\bf p}_{2}\cdot{\bf p}_{3}}{{\cal
M}(p_{3})}+\text{$2$ perms.}\bigg{)}\\\ &\;\;\;\;+b_{\phi}f_{\rm NL}\,f\mu
p_{123}\bigg{(}\frac{1}{3}G_{2}({\bf p}_{1},{\bf
p}_{2})\frac{\mu_{12}}{\lvert{\bf p}_{1}+{\bf p}_{2}\rvert}\frac{1}{{\cal
M}(p_{3})}+\text{$2$ perms.}\bigg{)}+b_{\phi}f_{\rm NL}\,(f\mu
p_{123})^{2}\bigg{(}\frac{1}{6}\frac{\mu_{1}\mu_{2}}{p_{1}p_{2}}\frac{1}{{\cal
M}(p_{3})}+\text{$2$ perms.}\bigg{)}\,\,.\end{split}$ (19)

Note that $b_{\phi\delta^{2}}$ and $b_{\phi{\cal G}_{2}}$ do not appear in
$Z^{\rm NG}_{3}$: the reason for this is that they are removed after
renormalization of $b_{1}$ and $b_{\phi}$. Moreover, the contributions from
$b_{\phi\delta}$ where either $\phi$ or $\delta$ are expanded at second order
in perturbations are also absorbed by renormalization of these two parameters,
was first proved in Ref. Moradinezhad Dizgah  _et al._ (2021) in the context
of the real space perturbation theory. Finally, let us note that compared to
the analysis of D’Amico  _et al._ (2022), we include the cubic non-Gaussian
kernel in the model, which is needed to calculate corrections to the one-loop
galaxy power spectrum induced by LPNG.

###  2.3 Stochasticity

So far we have focused on the deterministic part of the galaxy overdensity.
PNG leads to additional contributions to the stochastic part of $\delta_{g}$
as we discuss below.

Refs. Desjacques  _et al._ (2018a); Ivanov  _et al._ (2020a, 2021a) contain a
detailed description of stochastic terms in the case of Gaussian initial
conditions. As far as the tree-level bispectrum and one-loop power spectrum in
the presence of LPNG are concerned, the the full stochastic contribution is
given by Perko  _et al._ (2016); Ivanov  _et al._ (2021a)

$\delta_{g}^{\rm stoch.}({\bf
k})=\epsilon+\frac{d_{2}}{2}b_{1}[\delta\epsilon]_{\bf
k}-f\left[\partial_{z}(\epsilon u_{z})\right]_{\bf k}+f_{\rm
NL}\frac{b_{\epsilon\phi}}{2}[\epsilon\phi]_{\bf
k}+a_{0}R_{*}^{2}k^{2}\epsilon+a_{2}k_{z}\hat{z}^{i}(\epsilon^{i}+k^{i}\epsilon)+a_{4}k_{z}^{2}\hat{z}_{i}\hat{z}_{j}\epsilon^{ij}\,,$
(20)

where $\epsilon,~{}\epsilon^{i},~{}\epsilon^{ij}$ are the stochastic density,
velocity and tidal fields. The final three terms in (20) are higher derivative
stochastic contributions that are important only for the Gaussian part. The
only new LPNG contribution here is $\epsilon\phi$, where we emphasize again
that $\phi$ is evaluated at the Lagrangian position ${\bf q}$.

###  2.4 Summary of the power spectrum and bispectrum models

Once the new kernels in the presence of LPNG are obtained, it is
straightforward to compute the one-loop power spectrum and the tree-level
bispectrum. Modulo the counterterms, the deterministic part is given by

$\begin{split}P_{\text{ $1$-loop}}^{\rm tot}&=2\int_{\bf p}[Z^{\rm
tot}_{2}({\bf p},{\bf k}-{\bf p})]^{2}P_{11}(p)P_{11}(|{\bf k}-{\bf
p}|)+6Z^{\rm tot}_{1}({\bf k})P_{11}(k)\int_{\bf p}Z^{\rm tot}_{3}({\bf
k},-{\bf p},{\bf p})P_{11}(p)\\\ &\;\;\;\;+2Z^{\rm tot}_{1}({\bf k})\int_{\bf
p}Z^{\rm tot}_{2}({\bf p},{\bf k}-{\bf p})B_{111}(k,p,|{\bf k}-{\bf
p}|)\,\,,\\\ B^{\rm tot}_{\rm tree}&=Z^{\rm tot}_{1}({\bf k}_{1})Z^{\rm
tot}_{1}({\bf k}_{2})Z^{\rm tot}_{1}({\bf
k}_{3})B_{111}(k_{1},k_{2},k_{3})+2Z_{2}^{\rm tot}({\bf k}_{1},{\bf
k}_{2})Z^{\rm tot}_{1}({\bf k}_{1})Z^{\rm tot}_{2}({\bf
k}_{2})P_{11}(k_{1})P_{11}(k_{2})+\text{perms.}\,\,.\end{split}$ (21)

It is instructive to simplify this expression and break it down into separate
pieces. In this section we give explicit expressions for different terms in
Eq. (21). We focus on contributions that are linear in $f_{\rm NL}^{\rm
local}\Delta_{\phi}$. The contributions $\mathcal{O}((f_{\rm NL}^{\rm
local}\Delta_{\phi})^{2})$ can be straightforwardly computed, but turn out to
be irrelevant for our analysis (as shown in Section 2.5) except for the linear
$f_{\rm NL}^{2}$ scale-dependent bias term. We briefly discuss other
$\mathcal{O}((f_{\rm NL}^{\rm local}\Delta_{\phi})^{2})$ corrections in
Appendix A.

####  2.4.1 Power spectrum

The power spectrum has three kinds of additional contributions proportional to
$f_{\rm NL}$. At tree level we have

$P^{\rm LPNG}_{\text{tree-level}}=P_{11}^{f_{\rm NL}}+P_{11}^{f_{\rm
NL}^{2}}\,\,,$ (22)

while at one-loop order the total contribution linear in $f_{\rm NL}^{\rm
local}$ is given by

$P^{\rm LPNG}_{\text{$1$-loop}}=P_{22}^{f_{\rm NL}}+P_{13}^{f_{\rm
NL}}+P_{12}\,\,.$ (23)

The first contributions are the scale-dependent bias

$\begin{split}P_{11}^{f_{\rm NL}}(k,\mu)=2(b_{1}+f\mu^{2})b_{\phi}f_{\rm
NL}\frac{P_{11}(k)}{{\cal M}(k)}\,\,,\quad P_{11}^{f_{\rm
NL}^{2}}(k,\mu)=b^{2}_{\phi}f^{2}_{\rm NL}\frac{P_{11}(k)}{{\cal
M}^{2}(k)}\,\,.\end{split}$ (24)

Defining $P_{1\phi}(k)\equiv P_{11}(k)/{\cal M}(k)$, the $P_{22}^{f_{\rm NL}}$
contribution can be written as

$\begin{split}P_{22}^{f_{\rm NL}}({\bf k})=4f_{\rm NL}\int_{{\bf
p}}\tilde{Z}^{\rm NG}_{2}({\bf p},{\bf k}-{\bf p})Z_{2}({\bf p},{\bf k}-{\bf
p})P_{11}(p)P_{1\phi}(\lvert{\bf k}-{\bf p}\rvert)\,\,,\end{split}$ (25)

where we have introduced the new kernel

$\begin{split}\tilde{Z}^{\rm NG}_{2}({\bf p}_{1},{\bf
p}_{2})&=b_{\phi}\,\frac{{\bf p}_{1}\cdot{\bf
p}_{2}}{p_{1}^{2}}+b_{\phi}\,{f\mu
k}\frac{\mu_{1}}{p_{1}}+b_{\phi\delta}\,\,,\end{split}$ (26)

which is just a simplified version of $Z_{2}^{\rm NG}$. The use of $P_{1\phi}$
and $\tilde{Z}^{\rm NG}_{2}$ is particularly convenient for the FFTLog
evaluation of $P_{22}^{f_{\rm NL}}$, which we perform in this work following
the approach of Simonović  _et al._ (2018). $P_{13}^{f_{\rm NL}}$ is given by
the sum of three contributions:

$P_{13}^{f_{\rm NL}}=P_{13}^{f_{\rm NL}}{}^{(1)}+P_{13}^{f_{\rm
NL}}{}^{(2)}+P_{13}^{f_{\rm NL}}{}^{(3)}\,\,.$ (27)

The first is simply

$P_{13}^{f_{\rm NL}}{}^{(1)}=\frac{b_{\phi}f_{\rm NL}}{{\cal
M}(k)}\,\int_{{\bf p}}6Z_{3}({\bf k},{\bf p},{-{\bf p}})P_{11}(p)\,\,,$ (28)

with the second being

$P_{13}^{f_{\rm NL}}{}^{(2)}={-Z_{1}({\bf k})}b_{\phi}f_{\rm
NL}k^{2}(1+f^{2}\mu^{2})\sigma^{2}_{v}{P_{1\phi}(k)}\,\,,$ (29)

where $\sigma^{2}_{v}\equiv\int_{0}^{\infty}\mathrm{d}q\,P_{11}(q)/6\pi^{2}$.
This comes from the terms $\psi^{i}\psi^{j}\partial_{i}\partial_{j}\phi/2$ in
(15) and $\partial^{2}_{z}(u^{2}_{z}\phi)/2$ in the redshift-space mapping of
$b_{\phi}f_{\rm NL}\phi$, i.e. the third and fifth terms in (19). The
contribution $P_{13}^{f_{\rm NL}}{}^{(2)}$ exactly cancels with the IR limit
of the $P_{22}^{f_{\rm NL}}$ integral just like the IR limits of the $P_{13}$
and $P_{22}$ contributions in the Gaussian case. Finally, there is a term of
the form

$P_{13}^{f_{\rm NL}}{}^{(3)}=Z_{1}({\bf k})P_{11}(k)\int_{{\bf
p}}6\tilde{Z}^{\rm NG}_{3}({\bf k},{\bf p},{-{\bf p}})P_{11}(p)\,\,,$ (30)

where

$\begin{split}\tilde{Z}^{\rm NG}_{3}({\bf p}_{1},{\bf p}_{2},{\bf
p}_{3})&=b_{\phi}f_{\rm NL}\,\bigg{(}{-\frac{1}{14}}{\cal
G}_{2}(\bm{p}_{1},\bm{p}_{2})\frac{(\bm{p}_{1}+\bm{p}_{2})\cdot\bm{p}_{3}}{\lvert\bm{p}_{1}+\bm{p}_{2}\rvert^{2}}\frac{1}{{\cal
M}(p_{3})}+\text{$2$ perms.}\bigg{)}+b_{\phi}f_{\rm
NL}\,\bigg{(}\frac{1}{6}\frac{\bm{p}_{1}\cdot\bm{p}_{2}}{p^{2}_{1}p^{2}_{2}}\frac{\bm{p}_{2}\cdot\bm{p}_{3}}{{\cal
M}(p_{3})}+\text{$5$ perms.}\bigg{)}\\\ &\;\;\;\;+b_{\phi}f_{\rm NL}\,f\mu
p_{123}\bigg{(}\frac{1}{3}G_{2}({\bf p}_{1},{\bf
p}_{2})\frac{\mu_{12}}{\lvert{\bf p}_{1}+{\bf p}_{2}\rvert}\frac{1}{{\cal
M}(p_{3})}+\text{$2$ perms.}\bigg{)}\,\,.\end{split}$ (31)

Notice that the only permutations surviving are those for which the transfer
function remains inside the loop integral, i.e. it is $P_{1\phi}$ that is
integrated in $P_{13}^{f_{\rm NL}}{}^{(3)}$. The last term $P_{12}$ is given
by

$\begin{split}&P_{12}({\bf k})=12f_{\rm NL}Z_{1}({\bf k})\Delta_{\phi}\,{\cal
T}(k)\int_{{\bf p}}\big{[}{\cal S}(k,p,\lvert{\bf k}-{\bf p}\rvert)Z_{2}({\bf
p},{\bf k}-{\bf p})\big{]}{\cal T}(p){\cal T}(\lvert{\bf k}-{\bf
p}\rvert)\,\,,\end{split}$ (32)

where we have defined ${\cal T}(k)\equiv\Delta_{\phi}{\cal M}(k)/k^{2}$.

Finally, we note that the contribution from the stochastic term $\epsilon\phi$
to the one-loop power spectrum is degenerate with the stochastic shot noise
contributions we have in the zero-$f_{\rm NL}$ case. For this reason we do not
include it in the model.

Figure 1:  Left panel – “Gaussian” one-loop contributions to the power
spectrum monopole at $z=0.61$ compared with linear theory. We take $b_{1}=1$,
and the different curves have the corresponding bias parameters set to unity.
Right panel – PNG contributions to the power spectrum monopole at $z=0.61$
compared with linear theory. We take $f_{\rm NL}=100$ and $b_{1}=1$. The grey
curve shows the scale-dependent bias contribution for $b_{\phi}=1$. The
remaining curves show the different contributions ($P_{12}$ and
$P_{22}^{f_{\rm NL}}+P_{13}^{f_{\rm NL}}$) to $P^{\rm NG}_{\text{$1$-loop}}$
for unit values of the corresponding bias parameters.

The Gaussian and non-Gaussian one-loop corrections to the galaxy power
spectrum monopole are shown in the left and right panels of Fig. 1,
respectively, for $f_{\rm NL}=100$. We show all separate shapes without
multiplying them by the nuisance parameters to clearly illustrate the size of
these terms regardless of particular galaxy samples. The terms labeled “no
biases” correspond to pure matter contributions (i.e. they have $b_{1}=1$ and
all other biases set to zero, which corresponds a contribution from pure
matter). We see that some of the LPNG loops are actually as large as the
matter loops for $k\lesssim 0.1~{}h\text{Mpc}^{-1}$. Thus, these terms must be
included for consistency.

####  2.4.2 Bispectrum

Working at tree level in perturbations and at the linear order in $f_{\rm
NL}\Delta_{\phi}$, the PNG contributions to the redshift-space bispectrum are

$B^{\rm NG}_{\text{tree-level}}=f_{\rm NL}B_{211}^{f_{\rm
NL}}+B^{(s)}_{111}\,\,.$ (33)

$B_{211}^{f_{\rm NL}}$ arises from scale-dependent bias and is given by

$\begin{split}B_{211}^{f_{\rm NL}}({\bf k}_{1},{\bf k}_{2},{\bf
k}_{3})&=Z_{1}({\bf k}_{1})Z_{1}({\bf
k}_{2})b_{\delta\phi}\Big{[}P_{11}(k_{1})P_{1\phi}(k_{2})+P_{11}(k_{2})P_{1\phi}(k_{1})\Big{]}\\\
&\;\;\;\;+Z_{1}({\bf k}_{1})Z_{1}({\bf k}_{2})b_{\phi}\frac{{\bf
k}_{1}\cdot{\bf k}_{2}}{k_{1}k_{2}}\bigg{(}\frac{k_{2}}{k_{1}}\frac{1}{{\cal
M}(k_{2})}+\frac{k_{1}}{k_{2}}\frac{1}{{\cal
M}(k_{1})}\bigg{)}P_{11}(k_{1})P_{11}(k_{2})\\\ &\;\;\;\;+Z_{1}({\bf
k}_{1})Z_{1}({\bf
k}_{2})b_{\phi}\,f\mu_{12}k_{12}\bigg{(}\frac{\mu_{1}}{k_{1}}\frac{1}{{\cal
M}(k_{2})}+\frac{\mu_{2}}{k_{2}}\frac{1}{{\cal
M}(k_{1})}\bigg{)}P_{11}(k_{1})P_{11}(k_{2})\\\ &\;\;\;\;+2b_{\phi}Z_{2}({\bf
k}_{1},{\bf k}_{2})\Big{[}Z_{1}({\bf
k}_{1})P_{11}(k_{1})P_{1\phi}(k_{2})+Z_{1}({\bf
k}_{2})P_{11}(k_{2})P_{1\phi}(k_{1})\Big{]}+\text{2
permutations}\,\,.\end{split}$ (34)

$B^{(s)}_{111}$ is the standard tree-level redshift-space PNG contribution,

$B^{(s)}_{111}({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})=\prod_{a=1}^{3}Z_{1}({\bf
k}_{a})\mathcal{M}(k_{a})B_{\phi}(k_{1},k_{2},k_{3})=Z_{1}({\bf
k}_{1})Z_{1}({\bf k}_{2})Z_{1}({\bf k}_{3})\,f_{\rm
NL}\Delta_{\phi}~{}6~{}{\cal S}(k_{1},k_{2},k_{3})\,{\cal T}(k_{1}){\cal
T}(k_{1}){\cal T}(k_{1})\,\,.$ (35)

So far we have discussed the deterministic contributions. In contrast to the
power spectrum case, there is an additional stochastic contribution that is
not degenerate with the ones present also in the zero-$f_{\rm NL}$ case. This
contribution comes from (20), and takes the form

$\begin{split}&B^{f_{\rm NL}}_{\rm stoch}=f_{\rm
NL}b_{\phi}\mathcal{M}^{-1}(k_{1})Z_{1}({\bf
k}_{1})P_{11}(k_{1})P_{\epsilon}(k_{2})+f_{\rm
NL}b_{\epsilon\phi}\mathcal{M}^{-1}(k_{1})Z_{1}({\bf
k}_{1})P_{11}(k_{1})P_{\epsilon}(k_{2})+\text{$5$ perms.}\,\,,\end{split}$
(36)

where $P_{\epsilon}(k)$ is the power spectrum of $\epsilon$. At leading order
it is proportional to the constant shot noise value $\bar{n}^{-1}$. As we will
argue shortly, these terms turn out to be irrelevant for our analysis.

Before closing this section, we also note that we have implemented IR
resummation for all the LPNG terms entering the power spectra and bispectra
models, following the formalism of time-sliced perturbation theory Blas  _et
al._ (2016b, a); Ivanov and Sibiryakov (2018); Vasudevan  _et al._ (2019).
After implementing both IR resummation and the Alcock-Pazcynski projection
effects Alcock and Paczynski (1979) in our models for the tree-level
bispectrum and the one-loop power spectra, we numerically compute the Legendre
multipoles of the power spectrum and the bispectrum monopole, allowing for
robust comparison to data.

###  2.5 Behavior in a scaling universe

Let us estimate the relative importance of the different $f_{\rm NL}$
contributions. This can be done using the scaling universe approach Pajer and
Zaldarriaga (2013); Assassi  _et al._ (2015a). It is based on the fact that
the linear power spectrum in our Universe can be well approximated by a power
law: $P_{11}\propto(k/k_{\rm NL})^{n}k_{\rm NL}^{-3}$ with $n\approx-1.5$ for
quasi-linear wavenumbers $k\simeq 0.1\,h{\rm Mpc}^{-1}$. We also introduced
the nonlinear scale $k_{\rm NL}=0.5~{}h\text{Mpc}^{-1}$ at $z=0.5$.

We choose to focus on this particular range for the following reason. Given
that the leading LPNG contribution is a linear scale-dependent bias enhanced
on large scales, and the LPNG loop corrections dominate the usual Gaussian
loops at low-$k$, large scales should be crucial for our analysis. The
relative contributions of these terms diminish compared to the Gaussian loops
at small scales, but the errorbars also get smaller. This suggests that the
relative importance of the LPNG corrections should be maximal at some
intermediate wavenumber scale, which we choose we to be $k_{\rm
ref}=0.1~{}h\text{Mpc}^{-1}$, roughly in the center of the wavenumber range
that we use in the data analysis. In what follows, all estimates will be
presented for $k=k_{\rm ref}$.

Assuming that there is a single non-linear scale in the problem, the estimates
for the total dimensionless galaxy power spectrum $\Delta^{2}(k)\equiv
k^{3}P(k)$ for purely Gaussian initial conditions give

$\begin{split}\Delta^{2}(k)&={\underbrace{\left(\frac{k}{k_{\rm
NL}}\right)^{1.5}}_{P_{\text{tree}}}}+\underbrace{\left(\frac{k}{k_{\rm
NL}}\right)^{3}}_{P_{\text{$1$-loop}}}+\underbrace{\left(\frac{k}{k_{\rm
NL}}\right)^{3.5}}_{\rm ctr}+\underbrace{\left(\frac{k}{k_{\rm
NL}}\right)^{3}}_{\rm stoch}\,\,.\end{split}$ (37)

Recalling that the Bardeen potential has a nearly scale-invariant spectrum, we
get the following expressions for the LPNG terms:

$\begin{split}\Delta^{2}_{\rm NG}(k)&=\underbrace{f_{\rm
NL}\Delta_{\phi}\left(\frac{k}{k_{\rm NL}}\right)^{0.75}}_{P_{\text{tree-
level}}^{\text{NG,}~{}f_{\rm NL}}}+\underbrace{(f_{\rm
NL}\Delta_{\phi})^{2}}_{P_{\text{tree-level}}^{\text{NG,}~{}f_{\rm
NL}^{2}}}+\underbrace{f_{\rm NL}\Delta_{\phi}\left(\frac{k}{k_{\rm
NL}}\right)^{2.25}}_{P_{\text{$1$-loop}}^{\rm NG}}\,\,.\end{split}$ (38)

Evaluating these corrections at the reference scale $k_{\rm
ref}=0.1~{}h\text{Mpc}^{-1}$, we get

$\begin{split}&\Delta^{2}_{P_{\text{tree}}}\simeq 0.089\,\,,\\\
&\Delta^{2}_{P_{\text{tree-level}}^{\text{NG,}~{}f_{\rm NL}}}\simeq 1.1\cdot
10^{-2}\times\frac{f_{\rm NL}}{300}\,\,,\\\
\end{split}\qquad\begin{split}&\Delta^{2}_{P_{\text{1-loop}}}=\Delta^{2}_{P_{\text{stoch}}}\simeq
8\cdot 10^{-3}\,\,,\\\ &\Delta^{2}_{P_{\text{tree-level}}^{\text{NG,}~{}f_{\rm
NL}^{2}}}\simeq 1.3\cdot 10^{-3}\times\frac{f_{\rm NL}}{300}\,\,,\\\
\end{split}\qquad\begin{split}&\Delta^{2}_{P_{\text{ctr}}}\simeq 3.6\cdot
10^{-3}\,\,,\\\ &\Delta^{2}_{P^{\rm NG}_{\text{1-loop}}}\simeq 9.6\cdot
10^{-4}\times\frac{f_{\rm NL}}{300}\,\,.\end{split}$ (39)

As expected, we see that the scale-dependent bias contribution $P_{\text{tree-
level}}^{\rm LPNG}$ always dominates over $P_{\text{$1$-loop}}^{\rm LPNG}$ ,
and it is the main source of constraining power in the power spectrum data.
For $f_{\rm NL}\lesssim 300$ typical for our analyses we also see that the
one-loop PNG contributions are a small fraction of the “Gaussian”
$P_{\text{$1$-loop}}$.

The leading correction to the above result is given by the Gaussian two-loop
contribution, which can be estimates as

$\Delta^{2}_{\rm 2-loop}=\left(\frac{k_{\rm ref}}{k_{\rm
NL}}\right)^{4.5}\simeq 7.2\cdot 10^{-4}\,\,.$ (40)

This can be contrasted with the terms that we have dropped. Higher derivative
$f_{\rm NL}$ contributions stemming down from terms like $\partial^{2}\phi$
would be suppressed compared to other 1-loop LPNG terms that we retain in the
theory model,

$\Delta^{2}_{\langle\partial^{2}\phi\delta\rangle}=f_{\rm
NL}\Delta_{\phi}\left(\frac{k_{\rm ref}}{k_{\rm NL}}\right)^{2.75}\simeq
4.3\cdot 10^{-4}\times\left(\frac{f_{\rm NL}}{300}\right)\,\,.$ (41)

This justifies our choice of dropping these terms in Section 2.2.

We can make a similar argument for the loop terms $\mathcal{O}(f_{\rm
NL}^{2})$, which are also suppressed,

$\Delta^{2}_{\langle\phi\delta^{2}\rangle}=(f_{\rm
NL}\Delta_{\phi})^{2}\left(\frac{k_{\rm ref}}{k_{\rm NL}}\right)^{1.5}\simeq
1.1\cdot 10^{-4}\times\left(\frac{f_{\rm NL}}{300}\right)^{2}\,\,.$ (42)

All in all, our scaling universe estimates confirm that the one-loop LPNG
corrections can be important in the data analysis Moradinezhad Dizgah  _et
al._ (2021). In addition, we also need to retain the leading $f_{\rm NL}^{2}$
tree-level power spectrum contribution.

For the squeezed-limit tree-level bispectrum, where the shape of LPNG plays
the most important role, it is straightforward to see that $B_{211}^{f_{\rm
NL}}$ and $B_{111}$ scale in the same way. Their relative importance with
respect to the “Gaussian” $B^{(s)}_{211}$ is the same as that of the scale-
dependent bias piece versus $P_{11}$ in the power spectrum. The contribution
from $B^{f_{\rm NL}}_{\rm stoch}$ is suppressed in the squeezed limit, and
thus we do not include it in the analysis.

###  2.6 LPNG bias parameters

In the context of the EFT of LSS, $b_{\phi}$ and $b_{\phi\delta}$ (as well as
usual bias parameters like $b_{1}$, $b_{2}$, etc.) should be treated as free
nuisance parameters and marginalized over in data analysis. However, there are
certain phenomenological models of dark matter halo formation, which predict
$b_{\phi}$ and $b_{\phi\delta}$ parameters as a function of the linear bias
$b_{1}$. These are known as “universality relations” Desjacques  _et al._
(2018a). For the relevant LPNG bias coefficients they predict

$b_{\phi}=2\delta_{c}(b_{1}-1)\,,\quad
b_{\phi\delta}=b_{\phi}-(b_{1}-1)+\delta_{c}\left[b_{2}-\frac{8}{21}(b_{1}-1)\right]\,,$
(43)

where $\delta_{c}=1.686$. A typical approach then is to assume that the same
relationship holds true even for galaxies. The universality relation is
routinely used in most of $f_{\rm NL}$ constraints from galaxy surveys.
However, the relationships (43) fail for galaxies Barreira  _et al._ (2020)
from realistic hydrodynamical simulations. The most accurate analysis to date
gives the following fits based on the state-of-the art galaxy formation
simulations Barreira (2020, 2022):

$b_{\phi}=2\delta_{c}(b_{1}-0.55)\,,\quad
b_{\delta\phi}=3.85-9.49b_{1}+3.44b_{1}^{2}\,\,.$ (44)

We adopt relationship (44) in our baseline analysis. As a cross check, we also
repeat our analysis for the vanilla universality relations (43). In Section 6
we go beyond any assumptions on the LPNG parameters and fit
$b_{\delta\phi}f_{\rm NL}^{\rm local}$ and $b_{\phi}f_{\rm NL}^{\rm local}$
directly from the data for the first time.

##  3 Data and Analysis details

Our analysis is based on the twelfth data release (DR12) Alam  _et al._ (2017)
of the Baryon Oscillation Spectroscopic Survey (BOSS). The galaxy clustering
data covers two redshift bins with effective centers $z=0.38,\,0.61$, for the
Northern and Southern galactic caps, resulting in four independent slices. The
BOSS DR12 release contains a total of $\sim 1.2\times 10^{6}$ galaxies in a
total volume of $6\,(h^{-1}\text{Gpc})^{3}$. From each data chunk, we extract
the redshift-space power spectrum multipole moments $P_{\ell}$ ($\ell=0,2,4$),
the real space power spectrum proxy $Q_{0}$ Ivanov  _et al._ (2021b), the
redshift-space bispectrum monopoles for triangle configurations within the
range of $k_{i}\in[0.01,0.08)\,h\text{Mpc}^{-1}$ (a total of 62 bispectrum
data points per data chunk), and the BAO parameters
$\alpha_{\parallel},\alpha_{\perp}$ measured from the post-reconstructed power
spectrum data using the method of Philcox  _et al._ (2020). Both power spectra
and bispectra are measured using window-free estimators (Philcox, 2021a, b),
thus we do not need to include the survey window function in our theoretical
model.

We use the data cuts222Note that we use $k_{\rm max}$ for $P_{\ell}$ that is
slightly larger than that adopted in Ref. Philcox  _et al._ (2021a). This is
because the particular choice of the data cut in that paper $k_{\rm
max}=0.2~{}h\text{Mpc}^{-1}$ was based on detecting biases in the cosmological
parameter posteriors. However, in contrast to Philcox  _et al._ (2021a), here
we fix all cosmological parameters, in which case the fit to $f_{\rm NL}^{\rm
local}$ is unbiased up to somewhat larger $k_{\rm max}$. We stress that this
choice is not essential for the purposes of our work, as the $f_{\rm NL}^{\rm
local}$ constraints are dominated by the linear LPNG bias and hence are
saturated at large scales.  $k_{\rm max}=0.25~{}h\text{Mpc}^{-1}$, $k_{\rm
min}=0.01~{}h\text{Mpc}^{-1}$ for $P_{\ell}$ and $k_{\rm
min}=0.25~{}h\text{Mpc}^{-1}$, $k_{\rm max}=0.27~{}h\text{Mpc}^{-1}$ for
$Q_{0}$, so that the two statistics are largely independent. We use lower
$k_{\rm max}$ for $Q_{0}$ because the two-loop corrections can be non-
negligible compared to the PNG contributions for $k_{\rm
max}>0.3~{}h\text{Mpc}^{-1}$. Note that our choice $k_{\rm
min}=0.01~{}h\text{Mpc}^{-1}$ for both the power spectra and bispectra is
conservative. We remove the first bin in order to limit systematic effects
related to stellar contamination and residual radial and atmospheric
systematics, as well as integral constraints.333 As shown in Fig. 19 of (Kalus
_et al._ , 2019), weight-based approaches to removing large-scale systematics
(such as those applied in (Ross  _et al._ , 2017) and herein) produce
comparable results to more sophisticated methods on comparatively large
scales. For $k\lesssim 0.01~{}h\text{Mpc}^{-1}$, the differences between
different approaches become significant (and the systematics become larger
than the statistical errors, if uncorrected), thus these modes are excised
from our analysis.

The power spectra and bispectra used in this work are extracted using the
window-free estimators Philcox (2021a, b). The covariances for our total
datavector
$\\{P_{0},P_{2},P_{4},Q_{0},B_{0},\alpha_{\parallel},\alpha_{\perp}\\}$ for
each data chunk are extracted from a suite of 2048 MultiDark-Patchy mocks
Kitaura  _et al._ (2016), using the standard empirical covariance matrix
estimator.444See Wadekar and Scoccimarro (2019); Wadekar  _et al._ (2020);
Philcox  _et al._ (2021b); Byun  _et al._ (2021); Biagetti  _et al._ (2021)
for alternative covariance matrix estimation techniques.

Figure 2:  Marginalized constraint on $f_{\rm NL}^{\rm local}$ from the mean
of $84$ Nseries simulations, with a total volume approximately $40\times$
larger than that of BOSS. Note that we do not rescale the covariance to the
BOSS volume, but use that appropriate for the entire Nseries volume, allowing
a robust probe of theoretical systematics. Here, we find $f_{\rm NL}^{\rm
local}=-4.9\pm 5.0$ at $68\%$ CL.

Figure 3:  Marginalized constraints on local-type primordial non-Gaussianity
from the BOSS power spectrum (blue) and power spectrum plus bispectrum (red).
We find $f_{\rm NL}^{\rm local}=9^{+33}_{-35}$ and $-33\pm 28$ in the two
cases respectively at $68\%$ CL, with the bispectrum tightening the
constraints by $\approx 20\%$. These are the main results of this work.

Our full-shape analysis matches the ones of Ivanov  _et al._ (2020a); Philcox
_et al._ (2020); Philcox and Ivanov (2021); Cabass  _et al._ (2022a). Unlike
these works, we explicitly fix all cosmological parameters to the Planck 2018
priors Aghanim  _et al._ (2018). This is done because in this work we are
interested only in the constraints on $f_{\rm NL}^{\rm local}$ from the BOSS
survey. Formally, this corresponds to a combination of the CMB power spectra
and BOSS power+bispectrum data. Thus, in our MCMC analysis we only fit $f_{\rm
NL}^{\rm local}$ with an infinitely large flat prior, plus the Gaussian EFT
nuisance parameters (encompassing biases, stochasticity, and counterterms).

As discussed in Section 2.6, our fiducial analysis fixes the PNG bias
coefficients to values predicted by the dark matter halo relations as
functions of the corresponding linear bias $b_{1}$ for each data chunk. This
choice is optional. In principle, we can fit both $b_{\phi}$ and
$b_{\phi\delta}$ directly from the data, but the current limits on these
parameters are not very constraining, as shown in Section 6. Therefore, they
are fully consistent with the fits from simulations (44), making it reasonable
to fix them for the primary purposes of this paper.

Our analysis is based on the publicly available class-pt code Chudaykin  _et
al._ (2020). Since we do not vary cosmology in this study, we compute the full
one-loop power spectrum corrections including the LPNG terms only once,
utilizing the Planck cosmology, and only vary the bias parameters and $f_{\rm
NL}^{\rm local}$ in the likelihood. For the Nseries mocks we recompute the
relevant templates to match the Nseries fiducial cosmology. We plan to
implement the full cosmology-dependent LPNG calculation in a future update of
class-pt, using which we will systematically study the sensitivity of the
$f_{\rm NL}^{\rm local}$ constraints to uncertainties in cosmological
parameters. Our Markov Chain Monte Carlo (MCMC) analysis is run with the
Montepython code (Brinckmann and Lesgourgues, 2019) and is based on the
previously-used public likelihoods.555Available at
github.com/oliverphilcox/full_shape_likelihoods.

##  4 Validation on Mocks

As a validation test, we apply our pipeline to Nseries mock catalogs. These
catalogs were used by the BOSS collaboration for internal validation tests
Alam  _et al._ (2017). The suite consists of 84 semi-independent simulation
boxes. The Nseries mocks are designed to reproduce the clustering signal of
the high-z NGC BOSS sample. Each box has a similar effective volume and mean
effective redshift $z=0.56$. We fit the mean of 84 Nseries boxes with the
covariance of one box, divided by 84. Effectively, this is equivalent to
fitting a dataset which is $\approx 40$ times larger than the BOSS survey.
Just like in the actual BOSS data analysis, we fix all the cosmological
parameters (to the true values used in the simulations), and vary only $f_{\rm
NL}^{\rm local}$ and nuisance parameters in the fit.

The Nseries mocks were produced for Gaussian initial conditions, which we can
recover with our pipeline. Indeed, we find $f_{\rm NL}^{\rm local}$ consistent
with zero,

$f_{\rm NL}^{\rm local}=-4.9\pm 5.0\quad\text{at 68\%\,CL}\,\,,$ (45)

with the 1d marginalized posterior shown in Fig. 3. Note that the mean is
expected to differ from zero to be different from zero by $(1-2)~{}\sigma$ due
to random fluctuations. This also gives us an estimate of the theory
systematic error, $\Delta f_{\rm NL}^{\rm local}|_{\rm syst}\lesssim 5$, which
is less than $0.2\sigma$ of the actual BOSS 1d marginalized statistical error
(rescaling by the square-root of the volume ratio) .

We stress that the main goal of our validation test is to estimate the bias
due to the theory systematic error. An alternative approach for validation is
to fit the mean of the mocks with the covariance that matches the overall
volume of the BOSS survey. This test, however, does not allow one to assess
the theory bias because the posterior distribution in that case is affected by
prior volume effects (arising from the priors necessarily imposed on nuisance
parameters), which can be as large as the actual theory bias. This obscures
the estimation of theory bias and can lead to wrong conclusions on the
validity of the fitting pipeline. For example, the prior volume effects
exactly cancel the theory bias on $\sigma_{8}$ for BOSS-like mocks Ivanov  _et
al._ (2020a); Chudaykin  _et al._ (2021a); Philcox and Ivanov (2021). Thus, if
one fits the mean of the mocks with the covariance matching the BOSS survey
volume, one can erroneously conclude that the theory model is valid even for
$k_{\rm max}=0.3~{}h\text{Mpc}^{-1}$, whereas fitting the same data with the
actual covariance of the simulation suggests that the theory systematic bias
on $\sigma_{8}$ becomes sizable in the analysis of multipoles $P_{\ell}$ for
$k_{\rm max}>0.20~{}h\text{Mpc}^{-1}$.

##  5 Results for BOSS

Figure 4:  Left panel: residual variations of the galaxy power spectrum
monopole $P_{0}$ w.r.t. variations of $f_{\rm NL}^{\rm local}$. These
variations have distinctive shape dependence and therefore can be constrained
by the data. Right panel: residual variations of the galaxy bispectrum
monopole w.r.t. variations of $f_{\rm NL}^{\rm local}$. Black, orange and
green dots denote the squeezed ($k_{3}=0.015~{}h\text{Mpc}^{-1}$,
$k_{1},k_{2}>k_{3}$), equilateral ($k_{1}=k_{2}=k_{3}$), and flattened
($k_{2}=k_{3}$, $2k_{2}=k_{1}+0.015~{}h\text{Mpc}^{-1}$) triangle
configurations, respectively.

Here, we present the main results from the combined analysis of four BOSS data
chunks. We start with the power spectra + BAO data alone, i.e. without the
bispectrum. In this case we obtain

$\begin{split}\text{$P_{\ell}$+$Q_{0}$+BAO:}\quad f_{\rm NL}^{\rm
local}=9_{-35}^{+33}\,\,,\quad\text{(68\%~{}CL)}\,\,,\quad-57<f_{\rm NL}^{\rm
local}<78\,\,,\quad\text{(95\%~{}CL)}\,\,.\end{split}$ (46)

The marginalized posterior distribution is shown in Fig. 3. These results are
similar, but somewhat stronger than those obtained previously from BOSS DR9,
$-45<f_{\rm NL}^{\rm local}<195$ (95% CL) Ross  _et al._ (2013). This suggests
that the constraints in the power spectrum are dominated by the linear scale-
dependent bias from large scales. Indeed, repeating our analysis for $k_{\rm
min}=0.05~{}h\text{Mpc}^{-1}$ we find constraints that are worse by a factor
of four,

$\text{ $k_{\rm min}=0.05~{}h\text{Mpc}^{-1}$~{}~{}:}\quad f_{\rm NL}^{\rm
local}=-120_{-140}^{+100}\,\,,\quad\text{(68\%~{}CL)}\,\,,\quad-353<f_{\rm
NL}^{\rm local}<140\,\,,\quad\text{(95\%~{}CL)}\,\,.$ (47)

These results are consistent with the expectation that LPNG constraints are
dominated by the scale-dependent bias, which is sensitive to the lowest
available bin in the survey (see also Moradinezhad Dizgah  _et al._ (2021)).

Note that there are important differences between our analysis and that of
Ross  _et al._ (2013). That work was based on the monopole power spectrum
moment of the BOSS DR9 CMASS sample and had a lower scale cut. In contrast to
that, we use all three power spectrum multipole moments plus the real space
proxy of the complete BOSS DR12 data sample, but impose a conservative scale
cut $k_{\rm min}=0.01~{}h\text{Mpc}^{-1}$, significantly reducing any
systematics caused by observational effects, such as galactic foregrounds.
Nevertheless, the results of the two analyses are pleasingly consistent.

The addition of the bispectrum monopole shrinks the errorbar on $f_{\rm
NL}^{\rm local}$ by $\simeq 20\%$,

$\text{$P_{\ell}$+$Q_{0}$+BAO+$B_{0}$:}\quad f_{\rm NL}^{\rm local}=-33\pm
28\,\,,\quad\text{(68\%~{}CL)}\,,\quad-88<f_{\rm NL}^{\rm
local}<23\,\,,\quad\text{(95\%~{}CL)},.$ (48)

We do not find any evidence for LPNG: the 95% CL limits are consistent with
zero. The final posterior distribution is presented in Fig. 3.

To estimate the dependence of our results on the LPNG priors, we have repeated
our analysis assuming the universality relations (43) instead of the more
accurate simulation-calibrated fits (44). We found somewhat weaker bounds,

$\text{Prior~{}\eqref{eq:bphi0}:}\quad f_{\rm NL}^{\rm local}=-50\pm
40\,\,,\quad\text{(68\%~{}CL)}\,,\quad-130<f_{\rm NL}^{\rm
local}<30\,\quad\text{(95\%CL)}\,.$ (49)

This weakening of the constraints is expected, since the universality
relations underpredict the actual values of $b_{\phi}$ by $\sim 30\%$ compared
to (44) for $b_{1}\approx 2$. Note that when the universality bias relations
(43) are used, the relative impact of the bispectrum is somewhat stronger: it
tightens the constraints by $\simeq 30\%$. The result without the bispectrum
in this case is $f_{\rm NL}^{\rm local}=64.7_{-60}^{+52}$.

Figure 5:  1d and 2d marginalized posteriors for $f_{\rm NL}^{\rm local}$ and
the normalized LPNG bias parameters $b_{\phi}^{(i)}f_{\rm NL}^{\rm local}$,
$b_{\phi\delta}^{(i)}f_{\rm NL}^{\rm local}$, extracted from the BOSS galaxy
power spectra and bispectra data. $i=1,2,3,4$ corresponds to the NGCz3, SGCz3,
NGCz1, and SGCz1 BOSS data slices, respectively.

It is instructive to study where our constraints originate from. To this end
we illustrate the effect of the variation of $f_{\rm NL}^{\rm local}$ on the
galaxy power spectrum and bispectrum monopole, showing the corresponding
residuals in Fig. 4. We focus on the NGC z3 data chunk, and show the residuals
w.r.t. the best-fit model with $f_{\rm NL}^{\rm local}=0$. The bispectrum
residuals are shown as a function of the triangle index which is defined by
the bin center $(k_{1},k_{2},k_{3})$. The bin centers here satisfy
$k_{i}\in[0.015,0.075]~{}h\text{Mpc}^{-1}$, encompassing all wavenumbers with
$k\in(0.01,0.08)~{}h\text{Mpc}^{-1}$. We additionally mark the squeezed,
equilateral, and flattened triangles with black, orange, and green dots on the
x-axis (note that our notation differs somewhat from Cabass  _et al._
(2022a)). The shaded region corresponds to data errors.

Let us first focus on the power spectrum. We observe that the largest
deviations take place both on large scales (due to the linear scale-dependent
bias), and on short scales due to the one-loop LPNG corrections. This is
another indication that the one-loop LPNG corrections should be included in
data analysis. Of course, their effect is washed out, to some extent, by
marginalization over the standard EFT nuisance parameters, which are important
at short scales. Nevertheless, the LPNG nonlinear corrections must to be
included for the sake of consistency, and additionally, the degeneracies are
greatly reduced in practice by the inclusion of higher-order statistics.

Let us now move on to the bispectrum. We see that the PNG contribution has a
very significant configuration-dependence. The LPNG terms peak at the squeezed
triangles. This behavior is very different from the typical effect of the
standard non-linear galaxy bias Cabass  _et al._ (2022a), which illustrates
that this particular type of non-Gaussianity does not strongly suffer from
degeneracy with the unknown galaxy formation details.

Let us compare our results with the $f_{\rm NL}^{\rm local}$ measurements
available in the literature. Our limit is a factor of ten worse than the
Planck 2018 constraint, $f_{\rm NL}^{\rm local}=-0.9\pm 5.1$ (68% CL) Akrami
_et al._ (2020). As we stressed before, our limit is better than the one
obtained from BOSS DR9 Ross  _et al._ (2013), which is equivalent to
$\sigma_{f_{\rm NL}^{\rm local}}\approx 60$. The main reasons for this
improvement are new data, the complete theory model for the power spectrum,
more accurate priors for LPNG bias parameters, and the large-scale galaxy
bispectrum, which is quite sensitive to the scale-dependent bias signal. Note
that our measurement has a precision somewhat worse but comparable with the
eBOSS quasars $f_{\rm NL}^{\rm local}=-12\pm 21$ (68% CL) Castorina  _et al._
(2019); Mueller  _et al._ (2021) (which boast a much longer redshift baseline,
and thus a substantially lower $k_{\rm min}$), and with WMAP, $f_{\rm NL}^{\rm
local}=37\pm 20$ (68% CL) Bennett  _et al._ (2013). We find comparable results
to those from an independent analysis of the BOSS power spectrum and
bispectrum (using a partial one-loop theory model for the latter statistic)
(D’Amico  _et al._ , 2022): $f_{\rm NL}^{\rm local}=-30\pm 29$, though, as
noted above, our analysis differs due to the use of a fully consistent theory
model and complete treatment of the survey window, allowing larger-scale
information to be robustly included. Finally, our measurements are somewhat
better than the ones coming from the UV luminosity function $f_{\rm NL}^{\rm
local}=71^{+426}_{-237}$ (95% CL) Sabti  _et al._ (2021a, b, 2022), although
they include information from scales with $k>0.3~{}h\text{Mpc}^{-1}$, which we
do not consider in our study.

##  6 LPNG bias parameters

The galaxy bispectrum allows us to measure $f_{\rm NL}^{\rm local}$ separately
from $b_{\phi}f_{\rm NL}^{\rm local}$. At the power spectrum level this is,
essentially impossible, because the constraints are dominated by scale-
dependent bias controlled by combination $b_{\phi}f_{\rm NL}^{\rm local}$
Moradinezhad Dizgah  _et al._ (2021). However, the bispectrum allows us to
extract $f_{\rm NL}^{\rm local}$ directly from the $B_{111}^{(s)}$ shape that
is generated by the matter clustering and not the LPNG bias. In this section
we present the constraints on $b_{\phi}f_{\rm NL}^{\rm local}$ and $f_{\rm
NL}^{\rm local}$ independently. Moreover, we also present constraints on the
quadratic scale-dependent bias term $b_{\phi\delta}f_{\rm NL}^{\rm local}$,
which shows up in the tree-level galaxy bispectrum with LPNG. We found that
$f_{\rm NL}^{\rm local}$ can be quite large in our chains if $b_{\phi}$ and
$b_{\phi\delta}$ are not fixed, and therefore we have included the $(f_{\rm
NL}^{\rm local})^{2}$ corrections to the tree-level galaxy bispectrum. For
these pieces, we keep the additional $b_{\phi^{2}}$ bias fixed to the
prediction of the universality relation, as discussed in Appendix A.

We fit $f_{\rm NL}^{\rm local}$ plus parameters $b_{\phi}f_{\rm NL}^{\rm
local}$, $b_{\phi\delta}f_{\rm NL}^{\rm local}$ for each independent data
chunk. Our results are displayed in Fig. 5 and Table 1. We observe that the
BOSS data can constrain the LPNG bias parameter only at the level
$\sigma_{b_{\phi}f_{\rm NL}^{\rm local}}\sim 3\times 10^{2}$ and
$\sigma_{b_{\phi\delta}f_{\rm NL}^{\rm local}}\sim 3\times 10^{3}$. We see
that most of the posteriors are compatible with zero values of corresponding
parameters within $95\%$ CL. However, the parameters $f_{\rm NL}^{\rm
local},b_{\phi}^{(1)}f_{\rm NL}^{\rm local},b_{\phi\delta}^{(3)}f_{\rm
NL}^{\rm local}$, overlap with zero only within $99\%$ CL of the marginalized
posterior. Inspecting the 2d marginalized contours (shown in Fig. 5) suggests
that this is a result of degeneracies between $f_{\rm NL}^{\rm local}$ and the
LPNG bias parameter combinations. Note that the resulting posteriors are also
significantly non-Gaussian, which implies that having non-zero $f_{\rm
NL}^{\rm local}$ at $95\%$ CL does not actually imply a detection at a
significance level equivalent to that of a Gaussian-distributed parameter at
$2\sigma$, as evidenced by the lack of detection of LPNG in the fiducial
analysis.

Param | best-fit | mean$\pm\sigma$ | 95% lower | 95% upper
---|---|---|---|---
$f_{\rm NL}^{\rm local}$ | $-720$ | $-676_{-250}^{+150}$ | $-1080$ | $-210$
$b_{\phi}^{(1)}f_{\rm NL}^{\rm local}$ | $-878$ | $-740_{-250}^{+280}$ | $-1280$ | $-210$
$b_{\phi\delta}^{(1)}f_{\rm NL}^{\rm local}$ | $1900$ | $4300_{-2700}^{+2700}$ | $-950$ | $9500$
$b_{\phi}^{(2)}f_{\rm NL}^{\rm local}$ | $170$ | $160_{-360}^{+390}$ | $-600$ | $900$
$b_{\phi\delta}^{(2)}f_{\rm NL}^{\rm local}$ | $5400$ | $5400_{-3400}^{+3200}$ | $-1100$ | $12000$
$b_{\phi}^{(3)}f_{\rm NL}^{\rm local}$ | $180$ | $350_{-250}^{+270}$ | $-170$ | $860$
$b_{\phi\delta}^{(3)}f_{\rm NL}^{\rm local}$ | $4400$ | $5000_{-2200}^{+2100}$ | $800$ | $9300$
$b_{\phi}^{(4)}f_{\rm NL}^{\rm local}$ | $26$ | $120_{-390}^{+420}$ | $-700$ | $900$
$b_{\phi\delta}^{(4)}f_{\rm NL}^{\rm local}$ | $-22$ | $290_{-2700}^{+2600}$ | $-5000$ | $5700$
Table 1: 1d marginalized limits for $f_{\rm NL}^{\rm local}$ and the
normalized LPNG bias parameters $b_{\phi}^{(i)}f_{\rm NL}^{\rm local}$,
$b_{\phi\delta}^{(i)}f_{\rm NL}^{\rm local}$, extracted from the BOSS galaxy
power spectra and bispectra data. $i=1,2,3,4$ corresponds to the NGCz3, SGCz3,
NGCz1, and SGCz1 BOSS data slices, respectively.

##  7 Conclusions

We have presented constraints on local primordial non-Gaussianity from the
BOSS full-shape galaxy clustering data. The two main novelties of our analysis
are (a) we use the full one-loop power spectrum model that includes all
necessary non-linear one-loop corrections generated by LPNG, and (b) we
include the consistently analyzed galaxy bispectrum, incorporating a full
treatment of all relevant theoretical and observational effects. We have found
that the latter improves $f_{\rm NL}^{\rm local}$ constraints by $20\%$
compared to the power spectrum analysis. Thus, our paper extends and
complements previous works on LPNG from the galaxy clustering data.

There are many ways in which our analysis can be improved. First, the
$k$-range can be expanded, including additional information from both small
and large scales. The maximum wavenumber used in the analysis $k_{\rm max}$
can be significantly enhanced by the addition of the redshift-space galaxy
two-loop power spectrum, one-loop bispectrum, as well as the tree-level
trispectrum. Partial calculations of these observables already exist in the
literature, e.g. Bertolini  _et al._ (2016); Eggemeier  _et al._ (2018);
Konstandin  _et al._ (2019), and we plan to incorporate them in our future
analyses. Going beyond perturbative analysis, additional constraints on
$f_{\rm NL}^{\rm local}$ can be obtained from the nonlinear regime using
consistency relations for LSS Creminelli  _et al._ (2013). They guarantee that
the local shape of the bispectrum in the squeezed limit is protected by the
equivalence principle and it can be though of as a feature which is very
distinct from anything that can be produced by the astrophysical processes.
Such feature can be extracted even when the short modes are in the nonlinear
regime, marginalizing over the standard nonlinear physics using theoretical
error Baldauf  _et al._ (2016); Chudaykin  _et al._ (2021b), similarly to what
was done in extracting the BAO feature from the broadband in Philcox  _et al._
(2020). In addition, we plan to increase the $k$-range also on the lower end,
by including modes with $k<0.01~{}h\text{Mpc}^{-1}$ that are omitted in the
present analysis. This is particularly important in order to enhance
constraints from the scale-dependent bias. A simple Fisher forecast indicates
that the errors on $f_{\rm NL}^{\rm local}$ from BOSS can improve by a factor
of two by including all low-$k$ modes. Whilst straightforward from a
theoretical point of view, this will require a detailed study of large-scale
systematics, such as the integral constraint (both global and radial),
foreground stars, atmospheric effects, seeing, and galactic extinction, which
can produce large-scale radial and angular distortions, e.g. (Ross  _et al._ ,
2017; Kalus  _et al._ , 2019). Such work will be of particular importance as
the survey volume increases, and the range of $f_{\rm NL}$ parameters allowed
by data tightens. Another important point to keep in mind is that the Gaussian
power spectra and bispectra likelihood, which is used in our analysis, is not
valid for the lowest $k$ bins and it can skew the constraints on $f_{\rm
NL}^{\rm local}$. For some recent reflections on how to deal with this
problem, see Wang  _et al._ (2019).

Second, it would be interesting to study the dependence of the result on the
priors on EFT nuisance parameters. Previous works Wadekar  _et al._ (2020);
Cabass  _et al._ (2022a) have found that marginalization over Gaussian
nuisance parameters leads to a very significant degradation of parameter
errorbars. It is important to understand to what extent this can be avoided
with priors on nuisance parameters extracted from high fidelity simulations.
We have also shown that our current constraints significantly rely on using
universality-like relationships for the LPNG bias parameters
$b_{\phi},b_{\phi\delta}$. We will study if using relationships is accurate
enough for simulated data with injected LPNG. A similar analysis was done in
Ref. Moradinezhad Dizgah  _et al._ (2021) for the case of real space halo
clustering.

Third, one may include redshift-space multipoles of the galaxy bispectrum
beyond the monopole moment, which we considered in this work Scoccimarro  _et
al._ (1999). Fourth, one should perform a systematic sensitivity forecast for
future surveys like DESI Aghamousa  _et al._ (2016), Euclid Amendola  _et al._
(2018), and MegaMapper Schlegel  _et al._ (2019). In particular, just based on
the ratio of BOSS and DESI volumes, one may expect improvements by a factor of
three, i.e. reaching $\sigma_{f_{\rm NL}^{\rm local}}\approx 10$. Some
forecasts have already been performed e.g. Karagiannis  _et al._ (2018);
Sailer  _et al._ (2021), but most of them have been based on simplistic
assumptions about the theoretical modeling of the power spectrum and the
bispectrum, and it is rare for forecasts to include both the power spectrum
and bispectrum in combination. It will be interesting to see if the inclusion
of all necessary non-linear corrections can impact the conclusions of these
works, i.e. to perform a fully consistent forecast similar to Chudaykin and
Ivanov (2019); Sailer  _et al._ (2021). Fifth, another important ingredient is
a systematic study of the properties of galaxy samples that will be targeted
by future surveys, e.g. emission line galaxies admit higher $k_{\rm max}$ and
therefore better LPNG measurements can be obtained from this sample Ivanov
(2021). Finally, it would be interesting to extend our analysis to the case of
projected statistics, which is motivated by photometric surveys like SPHEREx
Doré  _et al._ (2014) and the Vera Rubin observatory Ivezić  _et al._ (2019).
This analysis will naturally require including relativistic and full-sky
corrections, which can impact constraints on LPNG Camera  _et al._ (2015);
Alonso and Ferreira (2015); Di Dio  _et al._ (2017); Castorina and di Dio
(2022), and exploring to what extent the analysis based on correlation
functions is optimal and how does it compare to recent results obtained using
forward modeling Andrews  _et al._ (2022).

## Acknowledgments

We are grateful to Kazuyuki Akitsu and Azadeh Moradinezhad Dizgah for valuable
discussions. We are especially grateful to Alexandre Barreira for useful
discussions about priors on non-Gaussian biases and careful reading of the
manuscript. GC acknowledges support from the Institute for Advanced Study. The
work of MMI has been supported by NASA through the NASA Hubble Fellowship
grant #HST-HF2-51483.001-A awarded by the Space Telescope Science Institute,
which is operated by the Association of Universities for Research in
Astronomy, Incorporated, under NASA contract NAS5-26555. OHEP thanks Will
Coulton and Kendrick Smith for useful discussions and acknowledges support
from the Simons Foundation. Parameter estimates presented in this paper have
been obtained with the class-pt Boltzmann code Chudaykin  _et al._ (2020) (see
also Blas  _et al._ (2011)) interfaced with the Montepython MCMC sampler
Audren  _et al._ (2013); Brinckmann and Lesgourgues (2019). The triangle plots
are generated with the getdist package666
https://getdist.readthedocs.io/en/latest/ Lewis (2019).

##  Appendix A Higher order PNG corrections to the galaxy bispectrum

In this section we present, for completeness, the bispectrum corrections at
$\mathcal{O}(f_{\rm NL}^{2})$. To obtain them, we need first include a new
bias operator

$\delta_{g}^{\text{real space}}\supset\frac{1}{2}f_{\rm
NL}^{2}b_{\phi^{2}}\phi^{2}\,.$ (50)

The full $\mathcal{O}(f_{\rm NL}^{2})$ kernel is then given by

$\begin{split}&B_{NG}^{f_{\rm NL}^{2}}=Z_{1}({\bf k}_{1})Z_{1}({\bf
k}_{2})Z_{1}({\bf k}_{3})b_{\phi}\left(\frac{1}{Z_{1}({\bf
k}_{1})\mathcal{M}(k_{1})}+\frac{1}{Z_{1}({\bf
k}_{2})\mathcal{M}(k_{2})}+\frac{1}{Z_{1}({\bf
k}_{3})\mathcal{M}(k_{3})}\right)B_{111}(k_{1},k_{2},k_{3})\\\
&+\Bigg{[}b_{\phi}^{2}\left(\frac{Z_{1}({\bf
k}_{1})}{\mathcal{M}(k_{2})}+\frac{Z_{1}({\bf
k}_{2})}{\mathcal{M}(k_{1})}\right)\frac{({\bf k}_{1}\cdot{\bf
k}_{2})}{k_{1}k_{2}}\left(\frac{k_{1}}{k_{2}\mathcal{M}(k_{1})}+\frac{k_{2}}{k_{1}\mathcal{M}(k_{2})}\right)+fb^{2}_{\phi}\mu
k\left(\frac{\mu_{1}}{k_{1}\mathcal{M}(k_{2})}+\frac{\mu_{2}}{k_{2}\mathcal{M}(k_{1})}\right)\left(\frac{Z_{1}({\bf
k}_{1})}{\mathcal{M}(k_{2})}+\frac{Z_{1}({\bf
k}_{2})}{\mathcal{M}(k_{1})}\right)\\\ &+2b_{\phi}^{2}\frac{Z_{2}({\bf
k}_{1},{\bf
k}_{2})}{\mathcal{M}(k_{1})\mathcal{M}(k_{2})}+b_{\phi\delta}b_{\phi}\left(\frac{Z_{1}({\bf
k}_{1})}{M(k_{2})}+\frac{Z_{1}({\bf
k}_{2})}{\mathcal{M}(k_{1})}\right)\left(\frac{1}{\mathcal{M}(k_{1})}+\frac{1}{\mathcal{M}(k_{2})}\right)+b_{\phi^{2}}\frac{Z_{1}({\bf
k}_{1})Z_{1}({\bf
k}_{2})}{2\mathcal{M}(k_{1})\mathcal{M}(k_{2})}\Bigg{]}P_{11}(k_{1})P_{11}(k_{2})\,\,.\end{split}$
(51)

Finally, the universality relations for the quadratic operators dictate

$\begin{split}&b_{\phi^{2}}=4\delta_{c}(b_{2}^{L}\delta_{c}-2b_{1}^{L})\,\,,\quad
b_{\phi\delta}=b_{\phi}+(-b_{1}^{L}+\delta_{c}b_{2}^{L})=b_{\phi}-(b_{1}-1)+\delta_{c}\left[b_{2}-\frac{8}{21}(b_{1}-1)\right]\,\,,\end{split}$
(52)

where $b_{1}^{L}\equiv b_{1}-1$ and $b_{2}^{L}$ are Lagrangian bias
coefficients.

## References

  * Maldacena (2003) J. M. Maldacena, JHEP 05, 013 (2003), arXiv:astro-ph/0210603 .
  * Creminelli and Zaldarriaga (2004) P. Creminelli and M. Zaldarriaga, JCAP 10, 006 (2004), arXiv:astro-ph/0407059 .
  * Enqvist and Sloth (2002) K. Enqvist and M. S. Sloth, Nucl. Phys. B 626, 395 (2002), arXiv:hep-ph/0109214 .
  * Lyth and Wands (2002) D. H. Lyth and D. Wands, Phys. Lett. B 524, 5 (2002), arXiv:hep-ph/0110002 .
  * Moroi and Takahashi (2001) T. Moroi and T. Takahashi, Phys. Lett. B 522, 215 (2001), [Erratum: Phys.Lett.B 539, 303–303 (2002)], arXiv:hep-ph/0110096 .
  * Zaldarriaga (2004) M. Zaldarriaga, Phys. Rev. D 69, 043508 (2004), arXiv:astro-ph/0306006 .
  * Aghanim  _et al._ (2018) N. Aghanim _et al._ (Planck), (2018), arXiv:1807.06209 [astro-ph.CO] .
  * Vernizzi and Wands (2006) F. Vernizzi and D. Wands, JCAP 05, 019 (2006), arXiv:astro-ph/0603799 .
  * Senatore and Zaldarriaga (2012) L. Senatore and M. Zaldarriaga, JHEP 04, 024 (2012), arXiv:1009.2093 [hep-th] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A9 (2020), arXiv:1905.05697 [astro-ph.CO] .
  * Doré  _et al._ (2014) O. Doré _et al._ , (2014), arXiv:1412.4872 [astro-ph.CO] .
  * Ferraro  _et al._ (2019) S. Ferraro _et al._ , (2019), arXiv:1903.09208 [astro-ph.CO] .
  * Castorina  _et al._ (2020) E. Castorina _et al._ , (2020), arXiv:2002.05072 [astro-ph.IM] .
  * Dalal  _et al._ (2008) N. Dalal, O. Dore, D. Huterer, and A. Shirokov, Phys. Rev. D 77, 123514 (2008), arXiv:0710.4560 [astro-ph] .
  * Matarrese and Verde (2008) S. Matarrese and L. Verde, Astrophys. J. Lett. 677, L77 (2008), arXiv:0801.4826 [astro-ph] .
  * Slosar  _et al._ (2008) A. Slosar, C. Hirata, U. Seljak, S. Ho, and N. Padmanabhan, JCAP 08, 031 (2008), arXiv:0805.3580 [astro-ph] .
  * Xia  _et al._ (2010) J.-Q. Xia, M. Viel, C. Baccigalupi, G. De Zotti, S. Matarrese, and L. Verde, Astrophys. J. Lett. 717, L17 (2010), arXiv:1003.3451 [astro-ph.CO] .
  * Ross  _et al._ (2013) A. J. Ross, W. J. Percival, A. Carnero, G.-b. Zhao, M. Manera, A. Raccanelli, E. Aubourg, D. Bizyaev, H. Brewington, J. Brinkmann, J. R. Brownstein, A. J. Cuesta, L. A. N. da Costa, D. J. Eisenstein, G. Ebelke, H. Guo, J.-C. Hamilton, M. V. Magaña, E. Malanushenko, V. Malanushenko, C. Maraston, F. Montesano, R. C. Nichol, D. Oravetz, K. Pan, F. Prada, A. G. Sánchez, L. Samushia, D. J. Schlegel, D. P. Schneider, H.-J. Seo, A. Sheldon, A. Simmons, S. Snedden, M. E. C. Swanson, D. Thomas, J. L. Tinker, R. Tojeiro, and I. Zehavi, MNRAS 428, 1116 (2013), arXiv:1208.1491 [astro-ph.CO] .
  * Castorina  _et al._ (2019) E. Castorina _et al._ , JCAP 09, 010 (2019), arXiv:1904.08859 [astro-ph.CO] .
  * Mueller  _et al._ (2021) E.-M. Mueller _et al._ , (2021), arXiv:2106.13725 [astro-ph.CO] .
  * Scoccimarro  _et al._ (2004) R. Scoccimarro, E. Sefusatti, and M. Zaldarriaga, Phys. Rev. D 69, 103513 (2004), arXiv:astro-ph/0312286 .
  * Baldauf  _et al._ (2011) T. Baldauf, U. Seljak, and L. Senatore, JCAP 04, 006 (2011), arXiv:1011.1513 [astro-ph.CO] .
  * Moradinezhad Dizgah  _et al._ (2021) A. Moradinezhad Dizgah, M. Biagetti, E. Sefusatti, V. Desjacques, and J. Noreña, JCAP 05, 015 (2021), arXiv:2010.14523 [astro-ph.CO] .
  * Barreira (2020) A. Barreira, JCAP 12, 031 (2020), arXiv:2009.06622 [astro-ph.CO] .
  * Philcox (2021a) O. H. E. Philcox, Phys. Rev. D 103, 103504 (2021a), arXiv:2012.09389 [astro-ph.CO] .
  * Philcox (2021b) O. H. E. Philcox, Phys. Rev. D 104, 123529 (2021b), arXiv:2107.06287 [astro-ph.CO] .
  * Pardede  _et al._ (2022) K. Pardede, F. Rizzo, M. Biagetti, E. Castorina, E. Sefusatti, and P. Monaco, (2022), arXiv:2203.04174 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2015) H. Gil-Marín, J. Noreña, L. Verde, W. J. Percival, C. Wagner, M. Manera, and D. P. Schneider, Mon. Not. Roy. Astron. Soc. 451, 539 (2015), arXiv:1407.5668 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2017) H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura, S. A. Rodríguez-Torres, and M. D. Olmstead, Mon. Not. Roy. Astron. Soc. 465, 1757 (2017), arXiv:1606.00439 [astro-ph.CO] .
  * D’Amico  _et al._ (2022) G. D’Amico, M. Lewandowski, L. Senatore, and P. Zhang, (2022), arXiv:2201.11518 [astro-ph.CO] .
  * Scoccimarro  _et al._ (1999) R. Scoccimarro, H. M. P. Couchman, and J. A. Frieman, Astrophys. J. 517, 531 (1999), arXiv:astro-ph/9808305 [astro-ph] .
  * Scoccimarro (2000) R. Scoccimarro, Astrophys. J. 544, 597 (2000), arXiv:astro-ph/0004086 .
  * Sefusatti  _et al._ (2006) E. Sefusatti, M. Crocce, S. Pueblas, and R. Scoccimarro, Phys. Rev. D74, 023522 (2006), arXiv:astro-ph/0604505 [astro-ph] .
  * Sefusatti and Komatsu (2007) E. Sefusatti and E. Komatsu, Phys. Rev. D 76, 083004 (2007), arXiv:0705.0343 [astro-ph] .
  * Sefusatti (2009) E. Sefusatti, Phys. Rev. D 80, 123002 (2009), arXiv:0905.0717 [astro-ph.CO] .
  * Baldauf  _et al._ (2015) T. Baldauf, L. Mercolli, M. Mirbabayi, and E. Pajer, JCAP 1505, 007 (2015), arXiv:1406.4135 [astro-ph.CO] .
  * Angulo  _et al._ (2015) R. E. Angulo, S. Foreman, M. Schmittfull, and L. Senatore, JCAP 1510, 039 (2015), arXiv:1406.4143 [astro-ph.CO] .
  * Blas  _et al._ (2016a) D. Blas, M. Garny, M. M. Ivanov, and S. Sibiryakov, JCAP 1607, 028 (2016a), arXiv:1605.02149 [astro-ph.CO] .
  * Nadler  _et al._ (2018) E. O. Nadler, A. Perko, and L. Senatore, JCAP 02, 058 (2018), arXiv:1710.10308 [astro-ph.CO] .
  * Eggemeier  _et al._ (2018) A. Eggemeier, R. Scoccimarro, and R. E. Smith, (2018), arXiv:1812.03208 [astro-ph.CO] .
  * Eggemeier  _et al._ (2021) A. Eggemeier, R. Scoccimarro, R. E. Smith, M. Crocce, A. Pezzotta, and A. G. Sánchez, (2021), arXiv:2102.06902 [astro-ph.CO] .
  * Desjacques  _et al._ (2018a) V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rept. 733, 1 (2018a), arXiv:1611.09787 [astro-ph.CO] .
  * Desjacques  _et al._ (2018b) V. Desjacques, D. Jeong, and F. Schmidt, JCAP 1812, 035 (2018b), arXiv:1806.04015 [astro-ph.CO] .
  * Ivanov and Sibiryakov (2018) M. M. Ivanov and S. Sibiryakov, JCAP 1807, 053 (2018), arXiv:1804.05080 [astro-ph.CO] .
  * Oddo  _et al._ (2020) A. Oddo, E. Sefusatti, C. Porciani, P. Monaco, and A. G. Sánchez, JCAP 03, 056 (2020), arXiv:1908.01774 [astro-ph.CO] .
  * Oddo  _et al._ (2021) A. Oddo, F. Rizzo, E. Sefusatti, C. Porciani, and P. Monaco, JCAP 11, 038 (2021), arXiv:2108.03204 [astro-ph.CO] .
  * Ivanov  _et al._ (2021a) M. M. Ivanov, O. H. E. Philcox, T. Nishimichi, M. Simonović, M. Takada, and M. Zaldarriaga, (2021a), arXiv:2110.10161 [astro-ph.CO] .
  * Philcox and Ivanov (2021) O. H. E. Philcox and M. M. Ivanov, (2021), arXiv:2112.04515 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Cabass  _et al._ (2022a) G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, (2022a), arXiv:2201.07238 [astro-ph.CO] .
  * Perko  _et al._ (2016) A. Perko, L. Senatore, E. Jennings, and R. H. Wechsler, (2016), arXiv:1610.09321 [astro-ph.CO] .
  * Ivanov  _et al._ (2020a) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, JCAP 05, 042 (2020a), arXiv:1909.05277 [astro-ph.CO] .
  * D’Amico  _et al._ (2019) G. D’Amico, J. Gleyzes, N. Kokron, D. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, (2019), arXiv:1909.05271 [astro-ph.CO] .
  * Ivanov  _et al._ (2020b) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, Phys. Rev. D 101, 083504 (2020b), arXiv:1912.08208 [astro-ph.CO] .
  * Nishimichi  _et al._ (2020) T. Nishimichi, G. D’Amico, M. M. Ivanov, L. Senatore, M. Simonović, M. Takada, M. Zaldarriaga, and P. Zhang, Phys. Rev. D 102, 123541 (2020), arXiv:2003.08277 [astro-ph.CO] .
  * Cabass  _et al._ (2022b) G. Cabass, M. M. Ivanov, M. Lewandowski, M. Mirbabayi, and M. Simonović, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.08232 [astro-ph.CO] .
  * Taruya  _et al._ (2008) A. Taruya, K. Koyama, and T. Matsubara, Phys. Rev. D 78, 123534 (2008), arXiv:0808.4085 [astro-ph] .
  * Assassi  _et al._ (2015a) V. Assassi, D. Baumann, E. Pajer, Y. Welling, and D. van der Woude, JCAP 11, 024 (2015a), arXiv:1505.06668 [astro-ph.CO] .
  * Assassi  _et al._ (2015b) V. Assassi, D. Baumann, and F. Schmidt, JCAP 12, 043 (2015b), arXiv:1510.03723 [astro-ph.CO] .
  * Chudaykin  _et al._ (2020) A. Chudaykin, M. M. Ivanov, O. H. E. Philcox, and M. Simonović, Phys. Rev. D 102, 063533 (2020), arXiv:2004.10607 [astro-ph.CO] .
  * Schmidt  _et al._ (2019) F. Schmidt, F. Elsner, J. Jasche, N. M. Nguyen, and G. Lavaux, JCAP 01, 042 (2019), arXiv:1808.02002 [astro-ph.CO] .
  * Simonović  _et al._ (2018) M. Simonović, T. Baldauf, M. Zaldarriaga, J. J. Carrasco, and J. A. Kollmeier, JCAP 1804, 030 (2018), arXiv:1708.08130 [astro-ph.CO] .
  * Blas  _et al._ (2016b) D. Blas, M. Garny, M. M. Ivanov, and S. Sibiryakov, JCAP 1607, 052 (2016b), arXiv:1512.05807 [astro-ph.CO] .
  * Vasudevan  _et al._ (2019) A. Vasudevan, M. M. Ivanov, S. Sibiryakov, and J. Lesgourgues, JCAP 09, 037 (2019), arXiv:1906.08697 [astro-ph.CO] .
  * Alcock and Paczynski (1979) C. Alcock and B. Paczynski, Nature 281, 358 (1979).
  * Pajer and Zaldarriaga (2013) E. Pajer and M. Zaldarriaga, JCAP 08, 037 (2013), arXiv:1301.7182 [astro-ph.CO] .
  * Barreira  _et al._ (2020) A. Barreira, G. Cabass, F. Schmidt, A. Pillepich, and D. Nelson, JCAP 12, 013 (2020), arXiv:2006.09368 [astro-ph.CO] .
  * Barreira (2022) A. Barreira, JCAP 01, 033 (2022), arXiv:2107.06887 [astro-ph.CO] .
  * Ivanov  _et al._ (2021b) M. M. Ivanov, O. H. E. Philcox, M. Simonović, M. Zaldarriaga, T. Nishimichi, and M. Takada, (2021b), arXiv:2110.00006 [astro-ph.CO] .
  * Philcox  _et al._ (2020) O. H. E. Philcox, M. M. Ivanov, M. Simonović, and M. Zaldarriaga, JCAP 05, 032 (2020), arXiv:2002.04035 [astro-ph.CO] .
  * Philcox  _et al._ (2021a) O. H. E. Philcox, Z. Slepian, J. Hou, C. Warner, R. N. Cahn, and D. J. Eisenstein, (2021a), arXiv:2105.08722 [astro-ph.IM] .
  * Kalus  _et al._ (2019) B. Kalus, W. J. Percival, D. J. Bacon, E. M. Mueller, L. Samushia, L. Verde, A. J. Ross, and J. L. Bernal, Mon. Not. Roy. Astron. Soc. 482, 453 (2019), arXiv:1806.02789 [astro-ph.CO] .
  * Ross  _et al._ (2017) A. J. Ross _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 464, 1168 (2017), arXiv:1607.03145 [astro-ph.CO] .
  * Kitaura  _et al._ (2016) F.-S. Kitaura _et al._ , Mon. Not. Roy. Astron. Soc. 456, 4156 (2016), arXiv:1509.06400 [astro-ph.CO] .
  * Wadekar and Scoccimarro (2019) D. Wadekar and R. Scoccimarro, (2019), arXiv:1910.02914 [astro-ph.CO] .
  * Wadekar  _et al._ (2020) D. Wadekar, M. M. Ivanov, and R. Scoccimarro, Phys. Rev. D 102, 123521 (2020), arXiv:2009.00622 [astro-ph.CO] .
  * Philcox  _et al._ (2021b) O. H. E. Philcox, M. M. Ivanov, M. Zaldarriaga, M. Simonovic, and M. Schmittfull, Phys. Rev. D 103, 043508 (2021b), arXiv:2009.03311 [astro-ph.CO] .
  * Byun  _et al._ (2021) J. Byun, A. Oddo, C. Porciani, and E. Sefusatti, JCAP 03, 105 (2021), arXiv:2010.09579 [astro-ph.CO] .
  * Biagetti  _et al._ (2021) M. Biagetti, L. Castiblanco, J. Noreña, and E. Sefusatti, (2021), arXiv:2111.05887 [astro-ph.CO] .
  * Brinckmann and Lesgourgues (2019) T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. 24, 100260 (2019), arXiv:1804.07261 [astro-ph.CO] .
  * Chudaykin  _et al._ (2021a) A. Chudaykin, K. Dolgikh, and M. M. Ivanov, Phys. Rev. D 103, 023507 (2021a), arXiv:2009.10106 [astro-ph.CO] .
  * Bennett  _et al._ (2013) C. L. Bennett, D. Larson, J. L. Weiland, N. Jarosik, G. Hinshaw, N. Odegard, K. M. Smith, R. S. Hill, B. Gold, M. Halpern, E. Komatsu, M. R. Nolta, L. Page, D. N. Spergel, E. Wollack, J. Dunkley, A. Kogut, M. Limon, S. S. Meyer, G. S. Tucker, and E. L. Wright, ”Astrophys. J.” 208, 20 (2013), arXiv:1212.5225 [astro-ph.CO] .
  * Sabti  _et al._ (2021a) N. Sabti, J. B. Muñoz, and D. Blas, JCAP 01, 010 (2021a), arXiv:2009.01245 [astro-ph.CO] .
  * Sabti  _et al._ (2021b) N. Sabti, J. B. Muñoz, and D. Blas, (2021b), arXiv:2110.13161 [astro-ph.CO] .
  * Sabti  _et al._ (2022) N. Sabti, J. B. Muñoz, and D. Blas, Phys. Rev. D 105, 043518 (2022), arXiv:2110.13168 [astro-ph.CO] .
  * Bertolini  _et al._ (2016) D. Bertolini, K. Schutz, M. P. Solon, and K. M. Zurek, JCAP 06, 052 (2016), arXiv:1604.01770 [astro-ph.CO] .
  * Konstandin  _et al._ (2019) T. Konstandin, R. A. Porto, and H. Rubira, JCAP 11, 027 (2019), arXiv:1906.00997 [astro-ph.CO] .
  * Creminelli  _et al._ (2013) P. Creminelli, J. Noreña, M. Simonović, and F. Vernizzi, JCAP 12, 025 (2013), arXiv:1309.3557 [astro-ph.CO] .
  * Baldauf  _et al._ (2016) T. Baldauf, M. Mirbabayi, M. Simonović, and M. Zaldarriaga, (2016), arXiv:1602.00674 [astro-ph.CO] .
  * Chudaykin  _et al._ (2021b) A. Chudaykin, M. M. Ivanov, and M. Simonović, Phys. Rev. D 103, 043525 (2021b), arXiv:2009.10724 [astro-ph.CO] .
  * Wang  _et al._ (2019) M. S. Wang, W. J. Percival, S. Avila, R. Crittenden, and D. Bianchi, Mon. Not. Roy. Astron. Soc. 486, 951 (2019), arXiv:1811.08155 [astro-ph.CO] .
  * Aghamousa  _et al._ (2016) A. Aghamousa _et al._ (DESI), (2016), arXiv:1611.00036 [astro-ph.IM] .
  * Amendola  _et al._ (2018) L. Amendola _et al._ , Living Rev. Rel. 21, 2 (2018), arXiv:1606.00180 [astro-ph.CO] .
  * Schlegel  _et al._ (2019) D. J. Schlegel _et al._ , (2019), arXiv:1907.11171 [astro-ph.IM] .
  * Karagiannis  _et al._ (2018) D. Karagiannis, A. Lazanu, M. Liguori, A. Raccanelli, N. Bartolo, and L. Verde, Mon. Not. Roy. Astron. Soc. 478, 1341 (2018), arXiv:1801.09280 [astro-ph.CO] .
  * Sailer  _et al._ (2021) N. Sailer, E. Castorina, S. Ferraro, and M. White, (2021), arXiv:2106.09713 [astro-ph.CO] .
  * Chudaykin and Ivanov (2019) A. Chudaykin and M. M. Ivanov, JCAP 11, 034 (2019), arXiv:1907.06666 [astro-ph.CO] .
  * Ivanov (2021) M. M. Ivanov, Phys. Rev. D 104, 103514 (2021), arXiv:2106.12580 [astro-ph.CO] .
  * Ivezić  _et al._ (2019) v. Ivezić _et al._ (LSST), Astrophys. J. 873, 111 (2019), arXiv:0805.2366 [astro-ph] .
  * Camera  _et al._ (2015) S. Camera, M. G. Santos, and R. Maartens, Mon. Not. Roy. Astron. Soc. 448, 1035 (2015), [Erratum: Mon.Not.Roy.Astron.Soc. 467, 1505–1506 (2017)], arXiv:1409.8286 [astro-ph.CO] .
  * Alonso and Ferreira (2015) D. Alonso and P. G. Ferreira, Phys. Rev. D 92, 063525 (2015), arXiv:1507.03550 [astro-ph.CO] .
  * Di Dio  _et al._ (2017) E. Di Dio, H. Perrier, R. Durrer, G. Marozzi, A. Moradinezhad Dizgah, J. Noreña, and A. Riotto, JCAP 03, 006 (2017), arXiv:1611.03720 [astro-ph.CO] .
  * Castorina and di Dio (2022) E. Castorina and E. di Dio, JCAP 01, 061 (2022), arXiv:2106.08857 [astro-ph.CO] .
  * Andrews  _et al._ (2022) A. Andrews, J. Jasche, G. Lavaux, and F. Schmidt, (2022), arXiv:2203.08838 [astro-ph.CO] .
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 1107, 034 (2011), arXiv:1104.2933 [astro-ph.CO] .
  * Audren  _et al._ (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302, 001 (2013), arXiv:1210.7183 [astro-ph.CO] .
  * Lewis (2019) A. Lewis, (2019), arXiv:1910.13970 [astro-ph.IM] .