# Massive-ish Particles from Small-ish Scales: Non-Perturbative Techniques for
Cosmological Collider Physics from Large-Scale Structure Surveys

Samuel Goldstein<EMAIL_ADDRESS>Department of Physics, Columbia
University, New York, NY 10027, USA  Oliver H. E. Philcox  Department of
Physics, Columbia University, New York, NY 10027, USA  Simons Society of
Fellows, Simons Foundation, New York, NY 10010, USA  J. Colin Hill  Department
of Physics, Columbia University, New York, NY 10027, USA  Lam Hui  Department
of Physics, Columbia University, New York, NY 10027, USA

###### Abstract

Massive particles produced during inflation impact soft limits of primordial
correlators. Searching for these signatures presents an exciting opportunity
to uncover the particle spectrum in the inflationary epoch. We present non-
perturbative methods to constrain intermediate-mass scalars ($0\leq m/H<3/2$,
where $H$ is the inflationary Hubble scale) produced during inflation, which
give rise to a power-law scaling in the squeezed primordial bispectrum.
Exploiting the large-scale structure consistency relations and the separate
universe approach, we derive models for the late-time squeezed matter
bispectrum and collapsed matter trispectrum sourced by these fields. To
validate our models, we run $N$-body simulations with the “Cosmological
Collider” squeezed bispectrum for two different particle masses. Our models
yield unbiased constraints on the amplitude of non-Gaussianity, $f_{\rm
NL}^{\Delta}$, from the squeezed bispectrum and collapsed trispectrum deep
into the non-linear regime ($k_{\rm max}\approx 2~{}h/{\rm Mpc}$ at $z=0$). We
assess the information content of these summary statistics, emphasizing the
importance of sample variance cancellation in the matter sector. We also study
the scale-dependent halo bias in our simulations. For mass-selected halos, the
non-Gaussian bias estimated from our simulations agrees with predictions based
on (i) separate universe simulations and (ii) universal mass functions. With
further work, these results can be used to search for inflationary massive
particle production with upcoming galaxy surveys.

##  I Introduction

In the standard model of cosmology, structures formed from quantum
fluctuations that were stretched to macroscopic scales during a period of
exponential expansion known as inflation [1, 2, 3]. Since the Hubble scale
during inflation may have been as high as $10^{14}~{}{\rm GeV}$, understanding
the physics responsible for inflation could provide insights into physics
beyond the standard model. Remarkably, the statistical properties of the
primordial fluctuations encode substantial information about their physical
origins. For instance, while the simplest single-field models of inflation
produce essentially Gaussian fluctuations [4, 5, 6, 7]111The simplest, single-
field, slow roll inflation does produce some degree of non-Gaussianity from
purely gravitational interactions [8]. However, as far as large-scale
structure is concerned, this level of non-Gaussianity can be safely ignored. ,
the presence of extra interactions and additional particles during inflation
can introduce non-Gaussianities in the primordial perturbations [9, 10]. The
search for this primordial non-Gaussianity (PNG) is a central goal of many
ongoing and upcoming cosmological surveys [11, 12, 13, 14, 15, 16].

Currently, the tightest constraints on various models of PNG come from
analyses of the cosmic microwave background (CMB) primary anisotropies [17].
While near-term CMB experiments are expected to improve constraints on PNG by
roughly a factor of $2-3$ [16], large-scale structure (LSS) surveys can, in
principle, provide more substantial improvements due to their large effective
volume and access to three-dimensional information [18, 19]. In practice,
extracting PNG from LSS surveys is a messy business because LSS data contain
non-Gaussianities due to, _e.g._ , non-linear structure formation, baryonic
complications, and redshift-space distortions. In the quasi-linear regime,
perturbative approaches have been used to robustly constrain PNG from galaxy
clustering observations [20, 21, 22, 23, 24]. Nevertheless, since LSS surveys
also provide precision measurements in the non-linear regime, it is valuable
to develop techniques to constrain PNG from non-linear LSS observables.

On small scales, non-linearities pose a significant theoretical roadblock
towards robust cosmological inference with LSS data. One possibility to
circumvent these non-linear nuisances is to exploit the LSS consistency
relations [25, 26] (see Refs. [27, 28, 29, 30, 31, 32] for more details),
which state that, in the absence of PNG and equivalence-principle-violating
physics, the ratio of the squeezed bispectrum to the long-wavelength (soft)
power spectrum is protected from inverse-powers-of-$q$ poles [33], _i.e._ ,

$\lim_{q\ll k}\frac{B(\mathbf{q},\mathbf{k})}{P(\mathbf{q})}\,\,\,\textrm{has
no {\hbox{1/q^\alpha}} poles}.$ (1)

Here, $B$ and $P$ are the equal-time, real-space matter bispectrum and power
spectrum, respectively 222In this work, we will assume $B$ and $P$ are
computed from the real-space matter field. However, we emphasize that Eq. (1)
is expected to hold for biased tracers of the matter distribution, such as
galaxies [26, 34, 29, 30, 35, 36]. Furthermore, the LSS consistency relations
take a similar form in redshift space [29, 28, 37]., $\mathbf{q},\mathbf{k}$
are the momenta (wavevectors) with $q\ll k$ (the squeezed limit), and
$\alpha>0$ is a real number. Since Eq. (1) is a non-perturbative consequence
of the symmetries of the LSS equations of motion, it applies for non-linear
$k$, outside the perturbative regime. Furthermore, Eq. (1) can be violated if
the primordial density perturbations are not Gaussian. Thus, one can constrain
PNG by searching for poles in the _non-linear_ squeezed $B/P$ ratio. Indeed,
Refs. [38, 39, 40] have used this approach to develop estimators for local PNG
that are valid in the non-linear regime. In this work, we generalize these
estimators to constrain PNG sourced by massive scalar particle exchange during
inflation.

Particles with mass $m\approx H$ produced during inflation leave a distinctive
imprint on the squeezed primordial bispectrum, characterized by power-law and
oscillatory features, depending on the particle’s mass and spin [41, 42, 43,
44]. Searching for these signatures forms the basis of the “Cosmological
Collider Physics” program [41], which aims to uncover the particle spectrum
during inflation. In the past decade, significant theoretical progress has
been made towards understanding and generalizing the Cosmological Collider
scenario [43, 45, 46, 47, 48, 42, 49, 41, 50, 44, 51, 52, 53, 54, 55, 56, 57,
58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73]. Conversely,
the observational imprints of the Cosmological Collider have only recently
begun to be explored [74, 75, 76, 77, 78, 79, 80]. The first dedicated
searches for the collider signal in cosmological data were conducted in the
past year [80, 24, 81]. Thus, there is room for improvement on the
observational front to optimally search for evidence of the Cosmological
Collider using upcoming CMB and LSS surveys.

In this work, we present a numerical investigation of the imprint of
intermediate-mass scalars333We use the term intermediate-mass to refer to
particles in the complementary series, $0\leq\frac{m}{H}<\frac{3}{2}$, which
imprint a power-law scaling in the squeezed bispectrum (see Eq. (6)). In
title-space, we prefer “massive-ish particles.” on LSS correlation functions.
Specifically, we run $N$-body simulations with the collider squeezed
bispectrum for two different particle masses. Our simulations use the same
settings as the Quijote-PNG simulations [82], but modified initial conditions.
Using the separate universe formalism, we derive non-perturbative models for
the squeezed matter bispectrum and collapsed matter trispectrum associated
with this type of PNG, which are then validated deep into the non-linear
regime using the simulations. We also use our simulations to study the large-
scale halo bias for this type of PNG. These results help pave the way towards
leveraging simulations and analytic models to search for signatures of the
Cosmological Collider with LSS surveys.

The remainder of this paper is structured as follows. In Sec. II, we review
the squeezed bispectrum in the Cosmological Collider scenario. We also derive
theoretical models for the squeezed matter bispectrum, collapsed matter
trispectrum, and scale-dependent bias associated with this scenario. In Sec.
III, we discuss our simulations and analysis pipeline. We validate our models
and simulations in Sec. IV. We conclude in Sec. V. The appendices give useful
theoretical results. In Appendix A, we present an algorithm for generating
$N$-body initial conditions with the Cosmological Collider bispectrum,
including spin and oscillations. In Appendix B, we derive the bispectrum and
trispectrum estimators used in this work.

Figure 1: The shape of the squeezed bispectrum (top) and collapsed trispectrum
(bottom) in momentum space. The squeezed bispectrum is characterized by three
external momenta, one of which, $q$, is much smaller than the other two,
$k\approx k^{\prime}$. The trispectrum is parameterized by four external
momenta and two internal momenta. In the collapsed limit, one of the internal
momenta, $q\equiv k_{12}$, is much smaller than all external momenta. Since
the soft mode is internal, the collapsed trispectrum is much less susceptible
to cosmic variance.

##  II Theoretical background

###  II.1 Conventions

We work in natural units, $c=1$. We assume a flat $\Lambda$CDM fiducial
cosmology based on the fiducial Quijote simulations [83]: $\Omega_{m}=0.3175$,
$\Omega_{b,0}=0.049$, $h=0.6711$, $n_{s}=0.9624$, and $\sigma_{8}=0.834.$ We
assume that all quantities are evaluated at a fixed redshift, $z$, which we
will typically omit in the expressions below.

We use $\delta_{m}({\bm{x}})\equiv\rho_{m}({\bm{x}})/\mathrm{bar}{\rho}_{m}-1$
to denote the matter density fluctuation and $\delta_{h}({\bm{x}})\equiv
n_{h}({\bm{x}})/\mathrm{bar}{n}_{h}-1$ to denote the halo number density
fluctuation. We typically work with the Fourier-space fields, which are
related to the configuration-space fields by

$\delta({{\bm{x}}})=\int\frac{d^{3}k}{(2\pi)^{3}}\,\delta(\mathbf{k})e^{-i\mathbf{k}\cdot{\bm{x}}}\equiv\int_{\mathbf{k}}\delta(\mathbf{k})e^{-i\mathbf{k}\cdot{\bm{x}}}.$
(2)

For a three-dimensional field $\delta(\mathbf{k})$, the power spectrum,
bispectrum, and trispectrum are defined by

$\displaystyle\langle\delta(\mathbf{k}_{1})\delta(\mathbf{k}_{2})\rangle^{\prime}_{c}$
$\displaystyle\equiv P(k_{1})\,,$ (3)
$\displaystyle\langle\delta(\mathbf{k}_{1})\delta(\mathbf{k}_{2})\delta(\mathbf{k}_{3})\rangle^{\prime}_{c}$
$\displaystyle\equiv B(k_{1},k_{2},k_{3})\,,$ (4)
$\displaystyle\langle\delta(\mathbf{k}_{1})\delta(\mathbf{k}_{2})\delta(\mathbf{k}_{3})\delta(\mathbf{k}_{4})\rangle^{\prime}_{c}$
$\displaystyle\equiv T(k_{1},k_{2},k_{3},k_{4},k_{12},k_{14})\,,$ (5)

respectively, where $\langle\dots\rangle^{\prime}_{c}$ denotes the connected
correlator with the overall momentum-conserving delta function stripped off.
The bispectrum is parameterized by the side lengths of three momentum vectors
forming a triangle. Similarly, the trispectrum is parameterized by the lengths
of the four sides and two diagonals spanning a tetrahedron. Finally, we define
the sum $\mathbf{k}_{12}\equiv\mathbf{k}_{1}+\mathbf{k}_{2}.$

In this work, we are interested in soft limits of correlation functions,
wherein one of the momenta is significantly smaller than all others. We will
always use $q$ to denote the soft mode. For the bispectrum, we analyze the
“squeezed limit,” where one of the external momenta, $q,$ is much smaller than
the other two, $k\approx k^{\prime}$. For convenience, we will often write the
bispectrum in terms of two of the three momenta as
$B(\mathbf{k}_{1},\mathbf{k}_{2})$. For the trispectrum, we focus on the
“collapsed limit,” where one of the internal momenta, $q\equiv k_{12}$, is
much smaller than all external momenta. See Fig. 1 for a graphical
representation of these limits.

###  II.2 Cosmological Collider bispectrum

The main goal of this work is to study LSS correlators in cosmologies with the
following squeezed primordial bispectrum,

$\lim_{q\ll k}B_{\Phi}(\mathbf{q},\mathbf{k})=4f_{\rm
NL}^{\Delta}\left(\frac{q}{k}\right)^{\Delta}P_{\Phi}(q)P_{\Phi}(k),$ (6)

where $f_{\rm NL}^{\Delta}$ is the amplitude of PNG,444We normalize $f_{\rm
NL}^{\Delta}$ such that $f_{\rm NL}^{\Delta=0}=f_{\rm NL}^{\rm loc.}$. Note
that this differs by a factor of three from the convention used in Green  _et
al._ [80], where their $f_{\rm NL}^{\Delta}$ is defined such that $f_{\rm
NL}^{\Delta=2}=f_{\rm NL}^{\rm eq.}$, thus $f_{\rm NL}^{\rm loc.}=3f_{\rm
NL}^{\Delta=0}$. $\Delta$ is a (possibly complex) exponent, and $P_{\Phi}$ is
the power spectrum of the primordial Bardeen potential, $\Phi$. Here, $q$ is
the “soft” mode and $k\approx k^{\prime}$ are the “hard” modes (see Fig. 1).

The squeezed bispectrum in Eq. (6) generically arises if there are
scalars555The discussion in this section can be extended to spin-$s$
particles. In this case, the squeezed bispectrum takes the following
anisotropic form [41],  $\lim_{q\ll k}B_{\Phi}(\mathbf{q},\mathbf{k})=4f_{\rm
NL}^{\Delta}\left(\frac{q}{k}\right)^{\Delta}\mathcal{L}_{s}(\hat{\mathbf{q}}\cdot\hat{\mathbf{k}})P_{\Phi}(q)P_{\Phi}(k),$
(7) where $\mathcal{L}_{s}(x)$ is the Legendre polynomial of order $s.$ The
exponent depends on both the spin and mass as follows,
$\Delta=3/2-\sqrt{\left(s-{1}/{2}\right)^{2}-m^{2}/H^{2}}\textrm{ for }s\neq
0.$ (8) We leave a study of spinning particles to a future work. with mass
comparable to the Hubble scale during inflation [41]. In this case, the
exponent is set by the mass of the field,

$\Delta=3/2-\sqrt{9/4-m^{2}/H^{2}},$ (9)

which has three interesting regimes:

  1. 1.

Light fields: For light fields ($m\ll H$), $\Delta\approx 0$, and we recover
the _local model_ associated with multi-field inflation models parameterized
by $f_{\rm NL}^{\rm loc.}$ [84].

  2. 2.

Intermediate-mass fields: For intermediate masses ($0\ll m\leq 3H/2$;
complementary series), $0<\Delta\leq 3/2$, and the squeezed bispectrum
includes a characteristic _power-law_ scaling. This squeezed bispectrum has
been analyzed in several specific models of inflation, most notably quasi-
single-field inflation (QSFI) [45, 43, 46, 42, 48, 85]; however, it is a
generic consequence of the presence of intermediate-mass particles during
inflation [41].

  3. 3.

Massive fields: For heavier fields ($m>3H/2$; principal series), $\Delta$ is
complex and there are _oscillations_ in the squeezed bispectrum. These
oscillations produce interesting signatures in LSS correlators [86, 76, 78,
77].

In what follows, we derive theoretical predictions for the squeezed matter
bispectrum, collapsed matter trispectrum, and scale-dependent bias sourced by
the primordial bispectrum in Eq. (6) for real values of $\Delta$. We defer an
investigation of oscillatory bispectra to a future work. We work in the peak-
background split formalism and follow a procedure similar to that used
previously to derive the scale-dependent bias effect for non-local PNG [87,
88, 89, 90, 77].

###  II.3 Non-perturbative model for the squeezed matter bispectrum

First, we derive a non-perturbative model for the squeezed matter bispectrum.
Our model generalizes the model derived for $f_{\rm NL}^{\rm loc.}$ in Ref.
[38] to account for the $(q/k)^{\Delta}$ momentum dependence in Eq. (6). The
central idea is to split the (squeezed) bispectrum into a term corresponding
to the gravitational-evolved contribution sourced by PNG and a term
corresponding to the contributions from gravitational non-Gaussianity that
would be present even in the absence of $f_{\rm NL}^{\Delta},$ _i.e._ ,

$B_{m}(\mathbf{q},\mathbf{k})=B_{m}^{f_{\rm
NL}^{\Delta}}(\mathbf{q},\mathbf{k})+B_{m}^{\rm
grav.}(\mathbf{q},\mathbf{k}).$ (10)

This is possible solely because the leading-order contribution to the
bispectrum associated with the squeezed primordial bispectrum in Eq. (6) is
proportional to $P_{m}(q)/q^{2-\Delta}$, which _cannot be generated by
gravitational evolution_ because of the LSS consistency relation, Eq. (1).

####  II.3.1 Primordial non-Gaussianity contribution

In the squeezed limit, the matter bispectrum can be described by the
correlator
$\langle{\delta_{m,L}(\mathbf{q})P_{m}(k|{\bm{x}})|_{\Phi_{L}}}\rangle^{\prime}$
where $\delta_{m,L}(\bm{q})$ is the long-wavelength (linear) matter density
field and $P_{m}(k|{\bm{x}})|_{\Phi_{L}}$ is the locally measured small-scale
(non-linear) matter power spectrum at position ${\bm{x}}$ in the presence of a
long-wavelength potential $\Phi_{L}$. For $q\ll k$, we can treat the long mode
as a background and expand the locally measured small-scale power spectrum as

$P_{m}(k|{\bm{x}})|_{\Phi_{L}}={}P_{m}(k)+\int_{\mathbf{p}}\,\frac{\partial
P_{m}(k)}{\partial\,\Phi_{L}(\mathbf{p})}\,\Phi_{L}(\mathbf{p})+\dots\,$ (11)

This yields the squeezed matter bispectrum

$\displaystyle\begin{split}B_{m}(\mathbf{q},\mathbf{k})&=\int_{\mathbf{p}}\frac{\partial
P_{m}(k)}{\partial\,\Phi_{L}(\mathbf{p})}\langle{\delta_{m,L}(\mathbf{q})\Phi_{L}(\mathbf{p})}\rangle\,,\\\
&=\frac{\partial
P_{m}(\mathbf{k})}{\partial\,\Phi_{L}(\mathbf{q})}\frac{P_{m}(q)}{\mathcal{M}(q,z)}\,.\end{split}$
(12)

Here, we used Poisson’s equation to relate the linear density field to the
long-wavelength gravitational potential,

$\delta_{m,L}(\mathbf{q},z)=\frac{2\,q^{2}\,T(q)D_{\rm
md}(z)}{3\Omega_{m}H_{0}^{2}}\Phi_{L}(\mathbf{q})\equiv\mathcal{M}(q,z)\Phi_{L}(\mathbf{q}),$

where $T(q)$ is the transfer function normalized to unity on large scales and
$D_{\rm md}(z)$ is the linear growth factor normalized to the scale factor
during matter domination, $\Omega_{m}$ is the present-day matter density, and
$H_{0}$ is the Hubble constant.

To proceed, we need to evaluate the potential derivative in Eq. (12). Since
the derivative acts on the _non-linear_ matter power spectrum, it is
challenging to compute using conventional techniques, such as perturbation
theory. Instead, the derivative can be estimated using the separate universe
approach [91], wherein we express the derivative in terms of a change in the
background cosmology to mimic the primordial mode coupling described by Eq.
(6). Specifically, according to Eq. (6) a long-wavelength potential
fluctuation on scale $\mathbf{q}$ induces a scale-dependent modification of
the small-scale linear density field [90, 77],

$\delta_{m,L}(\bm{k}|\Phi_{L}(\bm{q}))=\left(1+2f_{\rm
NL}^{\Delta}\left(\frac{q}{k}\right)^{\Delta}\Phi_{L}(\bm{q})\right)\delta_{m,L}(\bm{k}).$
(13)

Before continuing, we provide a graphical depiction of this mode coupling for
several values of $f_{\rm NL}^{\Delta}$ and $\Delta$. In Fig. 2, we show the
response of two small-scale density fluctuations with fixed amplitude, but
different momenta, to a long-wavelength potential fluctuation for several
values of $f_{\rm NL}^{\Delta}$ and $\Delta$. For Gaussian initial conditions,
the small-scale modes are uncorrelated with the long-wavelength potential. For
local non-Gaussianity, $\Delta=0$, the small-scale response is uniform. In
this case, the modulation is equivalent to changing the amplitude of
clustering, $\sigma_{8}.$ Finally, if $\Delta\neq 0$, then the modulation of
the small-scale fluctuations depends on their momenta.

Per Eq. (13), the locally measured, small-scale linear matter power spectrum
is modulated by the long-wavelength potential, $\Phi_{L}(\mathbf{q})$, as
follows,

$P_{m}^{\rm lin.}(k|{\bm{x}})|_{\Phi_{L}}=\left[1+4f_{\rm
NL}^{\Delta}\left(\frac{q}{k}\right)^{\Delta}\Phi_{L}(\mathbf{q})e^{i\mathbf{q}\cdot{\bm{x}}}\right]P_{m}^{\rm
lin.}(k).$ (14)

Treating $\mathbf{q}$ as a background, the mode coupling in Eq. (6) is
equivalent to a scale-dependent modification of the linear power spectrum,

$P_{m}^{\rm lin.}(k|\epsilon,\Delta)\equiv\left(1+2\epsilon
k^{-\Delta}\right)P_{m}^{\rm lin.}(k)$ (15)

for small parameter $\epsilon$ encoding the long-wavelength field. Thus, the
potential derivative is

$\frac{\partial P_{m}(k)}{\partial\,\Phi_{L}(\mathbf{q})}=2f_{\rm
NL}^{\Delta}q^{\Delta}\frac{\partial
P_{m}(k|\epsilon,\Delta)}{\partial\epsilon}\bigg{|}_{\epsilon=0},$ (16)

where $\partial P_{m}(k|\epsilon,\Delta)/\partial\epsilon|_{\epsilon=0}$ is
the derivative of the non-linear matter power spectrum with respect to a
change in the initial linear power spectrum described by Eq. (15). Since PNG
corresponds to mode coupling in the initial, _i.e._ , linear density field,
Eq. (16) can be evaluated in the non-linear regime using simulations with
modified initial conditions. We shall elaborate on this point later.

Figure 2: Graphical representation of the modulation of two small-scale
density fluctuations with _fixed amplitude, but different momenta,_ to a long-
wavelength gravitational potential, $\Phi_{L}({\bm{x}})$ (black-dashed) as
described by Eq. (13). For Gaussian initial conditions (top), the small-scale
modes are uncorrelated with $\Phi_{L}({\bm{x}})$. For local non-Gaussianity
(middle), the small-scale fluctuations respond uniformly to a long-wavelength
potential. This is equivalent to changing $\sigma_{8}$. Finally, for the
intermediate-mass particle scenario (bottom), the modulation depends on the
momenta of the small-scale modes.

Using Eq. (16), a non-perturbative model for the primordial contribution to
the squeezed matter bispectrum at leading order in $f_{\rm NL}^{\Delta}$ is

$\displaystyle B^{f_{\rm NL}^{\Delta}}_{m}(\mathbf{q},\mathbf{k})$
$\displaystyle=\frac{2f_{\rm
NL}^{\Delta}q^{\Delta}}{\mathcal{M}(q,z)}P_{m}(q)\frac{\partial
P_{m}(k|\epsilon,\Delta)}{\partial\epsilon}\bigg{|}_{\epsilon=0},$ (17)
$\displaystyle=\frac{3f_{\rm NL}^{\Delta}\Omega_{m}H_{0}^{2}}{D_{\rm
md}(z)}\frac{P_{m}(q)}{q^{2-\Delta}T(q)}\frac{\partial
P_{m}(k|\epsilon,\Delta)}{\partial\epsilon}\bigg{|}_{\epsilon=0}.$

We now make a couple of comments about Eq. (17). First, on large scales, the
ratio $B^{f_{\rm NL}^{\Delta}}_{m}(\mathbf{q},\mathbf{k})/P_{m}(q)\propto
1/q^{2-\Delta}$ (with $0\leq\Delta\leq 3/2$). This term cannot be generated by
gravitational evolution alone because it violates the LSS consistency
relation, Eq. (1). As a result, we can model the squeezed bispectrum
contributions sourced by gravitational non-Gaussianity separately from those
sourced by the Cosmological Collider.666This also applies for the bispectrum
templates associated with particles mass $m>3H/2$, as well as particles with
spin. Second, by expressing the potential derivative in terms of $\epsilon,$
we have written the bispectrum in a form that is of practical use. In
particular, we can compute the derivative in Eq. (17) by running simulations
with a modified initial linear power spectrum according to Eq. (15) for small
values of $\pm\epsilon$ and computing $\partial
P_{m}(k|\epsilon,\Delta)/\partial\epsilon|_{\epsilon=0}$ numerically from the
simulation output. We describe this procedure in more detail in Sec. III.1.2.

Finally, we check that Eq. (17) produces the expected behavior in known
limits. At tree level, we can compute the potential derivative from Eq. (15),
leading to

$B_{m}^{\rm tree}(\mathbf{q},\mathbf{k})=\frac{4f_{\rm
NL}^{\Delta}}{\mathcal{M}(q,z)}\left(\frac{q}{k}\right)^{\Delta}P_{m}(q)P_{m}(k),$
(18)

in agreement with Eq. (6). For local non-Gaussianity ($\Delta=0$) the
transformation in Eq. (15) is equivalent to rescaling
$\sigma_{8}\rightarrow\sigma_{8}(1+\epsilon)$. The potential derivative is
then $\partial P_{m}(k)/\partial\epsilon|_{\epsilon=0}=2\,\partial
P_{m}(k)/\partial\log\sigma_{8}^{2}$ and we recover the non-perturbative model
for the squeezed matter bispectrum in local PNG introduced in Ref. [38].

####  II.3.2 Gravitational contribution

Even in the absence of PNG, non-linear structure formation generates a
bispectrum. On sufficiently large scales, any primordial contribution of the
form given by Eq. (17) will dominate the equal-time squeezed bispectrum by the
LSS consistency relations. Nevertheless, for the range of scales probed by
cosmological surveys (and simulations) and for the values of $f_{\rm
NL}^{\Delta}\lesssim\mathcal{O}(100)$ [81] that are allowed with current
observations, the gravitational contribution is typically non-negligible.
Thus, to obtain unbiased constraints on $f_{\rm NL}^{\Delta}$ from the
squeezed matter bispectrum, we need a model also for the signal in the absence
of PNG.

Whereas the full matter bispectrum is difficult to model analytically beyond
perturbative scales, the squeezed matter bispectrum can be modeled non-
perturbatively using the response function approach [92, 91, 93, 94, 95, 33,
96, 97], _i.e._ ,

$\displaystyle B_{m}^{\rm grav.}(\bm{q},\bm{k})=$
$\displaystyle\,\left[a_{0}(k)+{a}_{2}(k)\frac{q^{2}}{k^{2}}+\cdots\right]P(k)P(q)\,,$
(19)

where $a_{0}$ and $a_{2}$ are coefficients characterizing the response of
small-scale matter clustering to a long-wavelength density perturbation. In
principle, these coefficients can be estimated from simulations (_e.g._ , [98,
97]) or estimated from theory (_e.g._ , [27, 99]); however, we follow the
approach in Refs. [33, 38], where we treat the coefficients as free parameters
that we marginalize over.777This definition of $a_{0}$ and $a_{2}$ differs
slightly from the definition used in Ref. [38], where we absorbed any factors
of $k$ and $P(k)$ into the coefficients themselves. The two approaches are
equivalent, but the convention adopted here has dimensionless coefficients. In
practice, we measure the squeezed bispectrum averaged over modes with $k_{\rm
min}<k,k^{\prime}<k_{\rm max}$; hence the response functions reduce to scalars
$\mathrm{bar}{a}_{0}$ and $\mathrm{bar}{a}_{2}$. We describe this procedure in
more detail in Sec. III.2.

Combining the contributions due to PNG, Eq. (17), and non-linear structure
formation, Eq. (19), our model for the squeezed matter bispectrum is

$B_{m}(q,k)=\frac{3f_{\rm NL}^{\Delta}\Omega_{m}H_{0}^{2}}{D_{\rm
md}(z)}\frac{P_{m}(q)}{q^{2-\Delta}T(q)}\frac{\partial
P_{m}(k|\epsilon,\Delta)}{\partial\epsilon}\bigg{|}_{\epsilon=0}+\left[a_{0}(k)+a_{2}(k)\frac{q^{2}}{k^{2}}\right]P(q)P(k)+\mathcal{O}\left((f_{\rm
NL}^{\Delta})^{2},q^{4}/k^{4}\right).$ (20)

###  II.4 Non-perturbative model for the collapsed matter trispectrum

Having derived a non-perturbative model for the squeezed matter bispectrum, we
now turn to the collapsed matter trispectrum. In this limit, the trispectrum
is a function of the magnitude of the soft mode, $q\equiv k_{12}$, as well as
two external momenta $k_{1}\approx k_{2}$ and $k_{3}\approx k_{4}$ (see Fig.
1). In the separate universe picture, the collapsed trispectrum contains the
following correlator
$\langle{P_{m}(k_{1}|{\bm{x}})|_{\Phi_{L}}P_{m}(k_{3}|{\bm{x}})|_{\Phi_{L}}}\rangle^{\prime}$
(see, _e.g._ , Refs. [100, 101]), which leads to the relation:

$\displaystyle T_{m}(q,k_{1},k_{3})$ $\displaystyle=\frac{\partial
P_{m}(k_{1})}{\partial\Phi_{L}(\mathbf{q})}\frac{\partial
P_{m}(k_{3})}{\partial\Phi_{L}(\mathbf{q})}P_{\Phi}(q),$
$\displaystyle=\frac{B_{m}(q,k_{1})B_{m}(q,k_{3})}{P_{m}(q)},$ (21)

where the second equality follows from Eq. (12). Thus, given a non-
perturbative model for the squeezed bispectrum, _i.e._ , Eq. (20), one can
readily derive a non-perturbative model for the collapsed trispectrum (though
see the discussion below). In practice, gravitational evolution or alternative
inflationary scenarios, such as $g_{\rm NL}$ PNG, produce contributions to the
trispectrum that are independent of the internal mode, $q$. For non-zero
$f_{\rm NL}^{\Delta}$, these contributions are generically suppressed on the
largest scales, where the collapsed trispectrum is dominated by the $(f_{\rm
NL}^{\Delta})^{2}P_{m}(q)/q^{4-2\Delta}$ contribution. Nevertheless, the
contact terms are non-negligible for the range of scales analyzed in this
work. To account for these contributions, we model the collapsed trispectrum
as,888For the models considered here, this equality holds on a realization-by-
realization basis. Thus, neglecting contact terms, the collapsed $T$ is set by
the _realization_ of the squeezed $B$ and soft $P$.

$T_{m}(q,k_{1},k_{3})=\frac{B_{m}(q,k_{1})B_{m}(q,k_{3})}{P_{m}(q)}+T_{0}(k_{1},k_{3}),$
(22)

where $B_{m}(q,k)$ is the squeezed bispectrum, given by Eq. (20), and
$T_{0}(k_{1},k_{3})$ encapsulates all $q$-independent contributions to the
collapsed trispectrum. Finally, we measure the trispectrum integrated over a
wide range of hard modes; thus $T_{0}(k_{1},k_{3})$ reduces to a scalar,
$\mathrm{bar}{T}_{0}$, which we marginalize over.

In some inflation models, the collapsed trispectrum can receive contributions
from exchange diagrams beyond the squeezed bispectrum. In the local model,
these contributions are conventionally parameterized by the amplitude
$\tau_{\rm NL}^{\rm loc.}$.999 We focus on $\tau_{\rm NL}^{\rm loc.}$ for
simplicity, but these results can be easily extended to non-zero values of
$\Delta$ by defining an amplitude $\tau_{\rm NL}^{\Delta}$. In full, such
additional contributions have a rich phenomenology, tracing, for example,
scalars that do not couple quadratically to the inflaton and non-longitudinal
spin-states. We ignore such effects in this work. A non-perturbative model for
the collapsed trispectrum including $f_{\rm NL}^{\rm loc.}$ and $\tau^{\rm
loc.}_{\rm NL}$ is

$T^{\rm loc.}_{m}(q,k_{1},k_{3})=\left(\frac{25}{36}\tau^{\rm loc.}_{\rm
NL}-(f_{\rm NL}^{\rm loc.})^{2}\right)\frac{36\,\Omega_{m}H_{0}^{2}}{D_{\rm
md}(z)^{2}}\frac{P_{m}(q)}{q^{4}T^{2}(q)}\frac{\partial
P_{m}(k_{1})}{\partial\log\sigma_{8}^{2}}\frac{\partial
P_{m}(k_{3})}{\partial\log\sigma_{8}^{2}}+\frac{B_{m}(q,k_{1})B_{m}(q,k_{3})}{P_{m}(q)}+T_{0}(k_{1},k_{3}),$
(23)

where the factor of $6/5$ is conventional and $B_{m}(q,k)$ is computed from
Eq. (20) with $\Delta=0.$

For local non-Gaussianity described by a quadratic rescaling of the primordial
potential ($\Phi\to\Phi+f_{\rm NL}^{\rm loc.}\Phi^{2}$), the collapsed
trispectrum is completely determined by the squeezed bispectrum, _i.e._ ,
$\tau_{\rm NL}^{\rm loc.}=\left(\frac{6}{5}f_{\rm NL}^{\rm loc.}\right)^{2}$.
In contrast, for inflationary scenarios in which multiple fields contribute to
the curvature perturbation, the amplitude of the collapsed trispectrum will be
boosted. Regardless of the inflationary model, $\tau_{\rm NL}^{\rm loc.}$ and
$f_{\rm NL}^{\rm loc.}$ must satisfy the Suyama-Yamaguchi inequality (_e.g._ ,
[102, 103, 104, 105, 42]),

$\tau_{\rm NL}^{\rm loc.}\geq\left(\frac{6}{5}f_{\rm NL}^{\rm
loc.}\right)^{2}.$ (24)

###  II.5 Scale-dependent bias

In addition to modifying higher-order correlation functions of the matter
field, the squeezed primordial bispectrum in Eq. (6) impacts the two-point
statistics of biased tracers, such as galaxies and halos. Specifically, the
correlation between the long-wavelength primordial potential and the locally
measured small-scale power spectrum induces spatial variations in the halo
number density, $n_{h}({\bm{x}})$. These variations are correlated with the
long-wavelength potential. Thus, the large-scale clustering of biased tracers
is a powerful probe of inflationary scenarios with a large squeezed primordial
bispectrum. A prime example of this is the “scale-dependent bias” in local PNG
[106, 107, 108, 109], where the large-scale bias includes a term proportional
to $f_{\rm NL}^{\rm loc.}/q^{2}.$ More general types of PNG, including the
Cosmological Collider bispectrum, induce their own scale-dependent bias [87,
88, 110, 111, 89, 90, 75, 77, 80]. Here, we review the scale-dependent bias
associated with the squeezed bispectrum in Eq. (6). For a more detailed
treatment, we refer the reader to Refs. [87, 110, 111, 89, 90, 77].

On large scales, the halo auto-spectrum and halo-matter cross-spectrum are
given by

$P_{h}(q)=\left[b_{1}+\frac{b_{\Psi,\Delta}f_{\rm
NL}^{\Delta}q^{\Delta}}{\mathcal{M}(q)}\right]^{2}P_{m}(q)+\frac{1}{\mathrm{bar}{n}_{h}}(1+\epsilon_{N}),$
(25) $P_{hm}(q)=\left[b_{1}+\frac{b_{\Psi,\Delta}f_{\rm
NL}^{\Delta}q^{\Delta}}{\mathcal{M}(q)}\right]P_{m}(q),$ (26)

respectively, where $b_{1}$ is the linear halo bias, $b_{\Psi,\Delta}$ is a
non-Gaussian bias parameter, $\mathrm{bar}{n}_{h}$ is the mean number density
of halos, and $\epsilon_{N}$ is a small parameter that encapsulates the
leading order effect of non-Poissonian shot noise. On large scales, the bias
is proportional to $f_{\rm NL}^{\Delta}/q^{2-\Delta}$ (and physically encodes
the same mechanism as the above squeezed bispectrum). Thus, in principle, one
can constrain $f_{\rm NL}^{\Delta}$ using the large-scale halo power spectrum.
In practice, since $f_{\rm NL}^{\Delta}$ is completely degenerate with the
non-Gaussian bias parameter, $b_{\Psi,\Delta}$, an estimate of
$b_{\Psi,\Delta}$ is necessary to robustly constrain $f_{\rm NL}^{\Delta}$
using the scale-dependent bias [112, 113, 114, 115, 116, 117]. In the
remainder of this section, we discuss several ways to estimate
$b_{\Psi,\Delta}$ that we will test with $N$-body simulations in Sec. IV.3.

####  II.5.1 Separate universe bias parameters

In the peak-background split formalism, the non-Gaussian bias parameter
$b_{\Psi,\Delta}$ describes the response of galaxy formation to a long-
wavelength potential perturbation. Specifically, for the squeezed bispectrum
in Eq. (6), the non-Gaussian bias is (_e.g._ , [90])

$b_{\Psi,\Delta}=\frac{2\,\partial\log\mathrm{bar}{n}_{h}}{\partial\epsilon}\bigg{|}_{\epsilon=0},$
(27)

where $\mathrm{bar}{n}_{h}$ is the mean number density of halos and the
derivative is taken with respect to a change in the linear power spectrum
according to Eq. (15). Physically, Eq. (27) can be understood as the analog of
the separate universe potential derivative (Eq. 16) that we derived for the
squeezed matter bispectrum model. Indeed, $b_{\Psi,\Delta}$ can be estimated
using the same set of separate universe simulations used to compute Eq. (16).
For local PNG, $\Delta=0$, Eq. (27) becomes
$b_{\Psi,\Delta}=4\,\partial\log\mathrm{bar}{n}_{h}/\partial\log\sigma_{8}^{2}$,
which is consistent with the separate universe simulations used in previous
studies of local PNG, _e.g._ , [113, 118].

####  II.5.2 Universal mass functions

Assuming a universal halo mass function (UMF),101010In this context, universal
means that the halo mass function depends only on the mean density of the
Universe, the variance of the density field smoothed on the scale of the
Lagrangian radius of the halo, $R_{*}(M)$, and the Jacobian of the
transformation relating the smoothed variance to the halo mass,
$J\equiv|\partial\log\sigma_{R_{*}}/\partial M|.$ See Sec. 7.3 of Ref. [90]
for more details. the non-Gaussian bias is (_e.g._ , [90])

$b_{\Psi,\Delta}^{\rm UMF}=\left[b_{\phi}^{\rm
UMF}+4\left(\frac{d\log\sigma_{R_{*},-\Delta/2}^{2}}{d\log\sigma_{R_{*}}^{2}}-1\right)\right]\frac{\sigma_{R_{*},-\Delta/2}^{2}}{\sigma_{R_{*}}^{2}},$
(28)

where $b_{\phi}^{\rm UMF}=2\delta_{c}(b_{1}-1)$ is the universality prediction
for the local PNG bias parameter and $\delta_{c}=1.686$ is the threshold for
spherical collapse linearly extrapolated to $z=0$. In Eq. (28), $R_{*}\equiv
R_{*}(M)$ denotes the Lagrangian radius of the halo, and $\sigma_{R,p}$ is the
momentum-weighted variance of the linear density field smoothed on some scale
$R$. Specifically,

$\sigma_{R,p}^{2}\equiv\int\frac{d^{3}k}{(2\pi)^{3}}k^{2p}W_{R}^{2}(k)P_{m}^{\rm
lin.}(k),$ (29)

where $W_{R}(k)=3j_{1}(kR)/(kR)$ is the Fourier transform of a spherical top-
hat filter with radius $R.$

Eq. (28) is practically useful as it provides a means to compute
$b_{\Psi,\Delta}$ directly from $b_{1}$ and the halo mass (which sets
$R_{*}$)111111The logarithmic derivative in Eq. (28) can be computed
numerically using finite differences with Eq. (29) at $R_{*}\rightarrow
R_{*}\pm\epsilon$.. However, recent studies in the context of local PNG have
shown that Eq. (28) can be a poor approximation for realistic halo and galaxy
samples [113, 114, 115, 116, 117]. Thus, further work is required to
understand when, if it all, Eq. (28) can be applied to realistic galaxy
surveys.

Finally, from Eq. (25), $b_{\Psi,\Delta}$ has units of $({\rm
Mpc}/h)^{\Delta}$. Thus, by dimensional analysis, the non-Gaussian bias
parameter can be approximated as

$b_{\Psi,\Delta}\approx 2\delta_{c}(b_{1}-1)R_{*}^{\Delta}.$ (30)

This approximation for the non-Gaussian bias has been used in several recent
works [75, 80].

##  III Methodology

###  III.1 Simulations

This section describes the details of the simulations used in this work. We
run two types of simulations to analyze non-linear LSS in the Cosmological
Collider scenario. First, we run simulations with PNG in the initial
conditions described by the squeezed bispectrum in Eq. (6) for several values
of $\Delta$ and $f_{\rm NL}^{\Delta}$. Second, for each value of $\Delta$, we
run separate universe simulations where we modify the input linear power
spectrum according to Eq. (15) for small values $\pm\epsilon$ to compute the
non-linear responses in the bispectrum, trispectrum, and scale-dependent bias
models. Table 1 summarizes the simulation parameters used in this work.

Aside from modifying the initial conditions, which we describe in the
following subsections, we generate all simulations using the same settings as
the Quijote suite [83, 82].121212https://quijote-
simulations.readthedocs.io/en/latest/ We refer the reader to Ref. [83] for a
detailed discussion of the Quijote simulations, but we review the main
information here. All simulations are run using the Gadget-3 treePM code [119]
and contain $512^{3}$ particles in a $(1~{}{\rm Gpc}/h)^{3}$ comoving volume.
We assume a fiducial cosmology with $\Omega_{m}=0.3175$, $\Omega_{b}=0.049$,
$h=0.6711$, $n_{s}=0.9624$, and $\sigma_{8}=0.834$ and compute all transfer
functions using CAMB [120]. We produce halo catalogs using the Friends-of-
Friends (FoF) algorithm [121] with a linking length of $b=0.2.$

PNG Simulations
---
$~{}~{}f_{\rm NL}^{\Delta}~{}~{}$ | $~{}~{}\Delta~{}~{}$ | $N_{\rm sim}$
100 | 0 | 50
300 | 0.5 | 50
1000 | 1 | 50

Separate Universe
---
$~{}~{}\pm\epsilon~{}[{\rm Mpc}/h]^{\Delta}~{}~{}$ | $~{}~{}\Delta~{}~{}$ | $N_{\rm sim}$
0.012 | 0 | 5
0.01 | 0.5 | 5
0.001 | 1 | 5

Table 1: Parameter values for the simulations used in this work. _Left_ :
settings for simulations with the squeezed primordial bispectrum in Eq. (6).
We run 50 realizations for each value of $\Delta$. For $\Delta=0$, we use the
first 50 realizations of the LC_p simulations from the Quijote-PNG suite [82].
_Right_ : settings for the separate universe simulations used to compute
response functions. For each value of $\Delta,$ we run five realizations with
the linear power spectrum in Eq. (15) with $\pm\epsilon$ reported in the
table. For the $\Delta=0$ case, we use the first five realizations of the s8_m
and s8_p Quijote simulations [83].

####  III.1.1 Initial conditions with the collider bispectrum

We generate initial conditions for the PNG simulations based on the method
developed in Ref. [88]. Specifically, we express the primordial potential as

$\Phi(\mathbf{k})=\phi(\mathbf{k})+f_{\rm
NL}^{\Delta}\left[\Psi_{\Delta}(\mathbf{k})-\langle\Psi_{\Delta}(\mathbf{k})\rangle\right],$
(31)

where $\phi(\mathbf{k})$ is a Gaussian random field, $f_{\rm NL}^{\Delta}$ is
the amplitude of non-Gaussianity, and $\Psi_{\Delta}(\mathbf{k})$ is an
auxiliary field that is quadratic in $\phi(\mathbf{k})$. Imposing homogeneity
and isotropy, the most general form for $\Psi_{\Delta}(\mathbf{k})$ is

$\Psi_{\Delta}(\mathbf{k})=\int\limits_{\mathbf{k}_{1},\mathbf{k}_{2}}(2\pi)^{3}\delta_{D}(\mathbf{k}-\mathbf{k}_{12})\,K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})\,\phi(\mathbf{k}_{1})\,\phi(\mathbf{k}_{2}),$

where $K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})$ is a coupling kernel whose
form is determined by demanding that the three-point function of
$\Phi(\mathbf{k})$ gives the desired bispectrum at leading order in $f_{\rm
NL}^{\Delta}.$ For a given bispectrum, the coupling kernel is not uniquely
determined [87, 88]. Nevertheless, it is important to choose
$K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})$ such that loop do not alter the
$k^{n_{s}-4}$ shape of the large-scale power spectrum. In this work, we use
the following kernel:131313In practice, we account for departures from scale
invariance by replacing all exponents
$\Delta\rightarrow\frac{\Delta}{3}(4-n_{s})$ [82].

$K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})=\left[\left(\frac{k_{1}}{k_{12}}\right)^{\Delta}+\left(\frac{k_{2}}{k_{12}}\right)^{\Delta}-1\right].$
(32)

The squeezed limit ($q\equiv k_{1}\ll k_{2}\approx k_{12}\equiv k$) of the
bispectrum associated with this kernel is, up to $\mathcal{O}\left(f_{\rm
NL}^{\Delta}\right)$,

$\frac{B_{\Phi}(\mathbf{q},\mathbf{k})}{P_{\Phi}(q)P_{\Phi}(k)}=4f_{\rm
NL}^{\Delta}\bigg{[}{\left(\frac{q}{k}\right)^{\Delta}+\left(\frac{k}{q}\right)^{\Delta}\frac{P_{\Phi}(k)}{P_{\Phi}(q)}}\bigg{]}.$
(33)

If the first term above dominates over the second term, then we recover the
squeezed bispectrum in Eq. (6).141414We emphasize that Eq. (32) generates the
target _squeezed bispectrum_ , which is all that is needed for our work. More
sophisticated techniques, (_e.g._ , [122]) are necessary to simulate non-
squeezed bispectrum configurations in the collider scenario. This also
requires a template for the full shape dependence of the collider bispectrum,
which is both model-dependent and very difficult. For a scale-invariant
primordial power spectrum, the second term is suppressed relative to the first
by $((q/k)^{3-2\Delta})$. In this work, we simulate $\Delta=0.5$ and 1; thus
we can safely ignore the second contribution.151515For simulations with larger
values of $\Delta$, both contributions in Eq. (33) would need to be included
in the theoretical model. Both solutions appear in the theoretical model for
the Cosmological Collider squeezed bispectrum, where
$\Delta^{\prime}=3-\Delta$ is a second solution to the equations of motion.
Indeed, this solution must be included for fields with $m>3H/2$ to ensure the
bispectrum is real. However, the amplitude, $f_{\rm NL}^{\Delta}$, is
generally different for the two solutions. Thus, Eq. (12) applies only when
the second solution is sub-dominant. The exception is for fields with $m=3H/2$
($\Delta=3/2$), for which Eq. (32) generates the target bispectrum with an
additional factor of two. Alternatively, one could use a high-pass filtered
kernel $K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})\rightarrow
K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})\Theta(k_{12}-k_{L})$, where $\Theta$
is the step function and $k_{L}\approx 0.1~{}h/{\rm Mpc}$. This filter would
suppress the second contribution to Eq. (33).

Assuming a scale-invariant primordial power spectrum, then the one-loop
contribution to $P_{\Phi}$ is

$P_{\Phi}^{1-{\rm
loop}}(k)\propto\int_{\mathbf{p}}\left(K_{\Delta}(\mathbf{p},\mathbf{k}-\mathbf{p})\right)^{2}p^{-3}(|\mathbf{k}-\mathbf{p}|)^{-3},$
(34)

which is proportional to $k^{-2\Delta}$ at low $k$. For $\Delta<3/2$, the one-
loop corrections are suppressed relative to the linear power spectrum at low
$k$. For $\Delta=3/2$, both the $P^{\rm lin.}_{\Phi}$ and $P_{\Phi}^{1-{\rm
loop}}$ scale like $k^{-3};$ however, since the one-loop contribution is
proportional to $(f_{\rm NL}^{\Delta})^{2}$, it is subdominant. Thus, Eq. (32)
does not significantly alter $P_{\Phi}.$

Finally, we can write the auxiliary field $\Psi_{\Delta}$ in a separable form
by replacing the Dirac delta with an exponential,

$\displaystyle{{\Psi_{\Delta}(\mathbf{k})}}=\int{d^{3}{{\bm{x}}}}\,e^{i\mathbf{k}\cdot{\bm{x}}}\bigg{[}\frac{2\,\phi({\bm{x}})}{k^{\Delta}}\int_{\mathbf{p}}{p}^{\Delta}\phi(\mathbf{p})e^{-i\mathbf{p}\cdot{\bm{x}}}-\phi^{2}({\bm{x}})\bigg{]},$

which can be efficiently evaluated using FFTs. In Appendix A, we show how to
generalize this kernel to account for fields with spin and $m>3H/2.$

We modify the code 2LPTPNG161616https://cosmo.nyu.edu/roman/2LPT/ to generate
initial conditions for the collider bispectrum. We run 50 simulations with
$f_{\rm NL}^{\Delta=0.5}=300$ and $f_{\rm NL}^{\Delta=1}=1000$. To isolate the
impact of PNG, we use the same seeds as the first 50 realizations of the
Quijote and Quijote-PNG simulations. For comparison, we also analyze the same
50 realizations of the fiducial Quijote simulations, which have Gaussian
initial conditions, as well as the LC_p Quijote-PNG simulations, which have
$f_{\rm NL}^{\rm loc.}=100.$

####  III.1.2 Separate Universe Simulations

We also run separate universe simulations to compute the potential derivative
(Eq. 16) and the non-Gaussian bias (Eq. 27). As described in Sec. II.3.1,
these derivatives can be estimated numerically from the output of simulations
with initial linear power spectra rescaled by $P(k)\rightarrow(1+2\epsilon
k^{-\Delta})P(k)$ for small values of $\pm\epsilon$.171717Some initial
condition codes recompute the normalization of the input power spectrum to
match the input value of $\sigma_{8}$. In this case, one should fix the
normalization to the value computed from the fiducial, _i.e._ , $\epsilon=0$
power spectrum. We run five realizations for each value of $\pm\epsilon$,
which we find sufficient to reduce cosmic variance for the small-scale modes
analyzed here. We list our values of $\epsilon$ in Table 1. For the
simulations with $f_{\rm NL}^{\rm loc.}$, _i.e._ , $\Delta=0$, Eq. (15) is
equivalent to a change in the amplitude of the linear power spectrum. Hence,
we use the s8_m and s8_p simulations from the Quijote suite, which are run
with $\sigma_{8}=0.819$ and $0.834$, respectively.

###  III.2 Measurements and likelihood

In this section, we describe the measurements and likelihood used to study the
squeezed matter bispectrum, collapsed matter trispectrum, and scale-dependent
halo bias for the squeezed bispectrum in Eq. (6) with $N$-body simulations.
Our method is based on the bispectrum analysis presented in Goldstein  _et
al._ [38]. We refer the reader there for a more detailed description of our
modeling and measurement choices. We use the PNGolin
package181818https://github.com/samgolds/PNGolin to compute the power
spectrum, bispectrum, and trispectrum from the simulation output. See Appendix
B for a detailed discussion of our estimators.

####  III.2.1 Measurements

We aim to study imprints of the Cosmological Collider on higher-order
statistics of the matter field, as well as the two-point statistics of the
halo field. To estimate these summary statistics, we first assign the redshift
$z=0$ matter particles and halos to a $1024^{3}$ mesh using Cloud-in-Cell
interpolation with the Pylians package
[123].191919https://pylians3.readthedocs.io/en/master/ We split the FoF halo
catalog into four logarithmically spaced mass bins with $13.5<\log_{10}(M_{\rm
FoF})<14.5~{}M_{\odot}/h$, as well as one bin with $M_{\rm
FoF}>10^{14.5}~{}M_{\odot}/h$. From these fields, we compute the matter power
spectrum, halo power spectrum, and halo-matter cross-spectrum in bins of width
$2k_{f}$, where $k_{f}\equiv 2\pi/L_{\rm box}\approx 0.006~{}h/{\rm Mpc}$ is
the fundamental mode. We always exclude the fundamental mode from our analyses
to avoid numerical issues.

We analyze the squeezed bispectrum integrated over a range of hard modes.
Specifically, we compute the squeezed bispectrum for soft modes $q>k_{f}$ in
bins of width $2k_{f}$ up to $q_{\rm max}=0.05~{}h/{\rm Mpc}$ and averaged
over hard modes in some bin $k_{\rm min}<k,k^{\prime}<k_{\rm max}$.202020This
definition differs slightly from the bispectrum estimator used in Ref. [38],
which required $k_{\rm min}<k<k_{\rm max}$ and let $k^{\prime}$ take any value
allowed by the triangle inequality. Not imposing the condition $k_{\rm
min}<k^{\prime}<k_{\rm max}$ leads to a somewhat enhanced $q$ dependence in
the squeezed bispectrum. This corresponds to $N_{q}=7$ soft mode bins for each
choice of $k_{\rm min}$ and $k_{\rm max}$. We use $B_{m}(q|\,\mathcal{K})$ to
denote the integrated squeezed bispectrum as a function of the soft mode, $q$,
averaged over _all_ hard modes $k,k^{\prime}\in\mathcal{K}$. In Sec. IV.2, we
compare the information content of the squeezed bispectrum integrated over a
single bin $\mathcal{K}$ with the squeezed bispectrum split into multiple
disjoint sub-bins, _e.g._ , $B_{m}(q|\,\mathcal{K}_{a})$ and
$B_{m}(q|\,\mathcal{K}_{b})$ where
$\mathcal{K}_{a}\cup\mathcal{K}_{b}=\mathcal{K}.$

Similar to the squeezed bispectrum measurements, we compute the collapsed
trispectrum for soft modes $q\equiv k_{12}>k_{f}$ in bins of width $2k_{f}$ up
to $q_{\rm max}=0.05~{}h/{\rm Mpc}$ and averaged over all hard modes in some
bin $k_{\rm min}<k_{1},k_{2},k_{3},k_{4}<k_{\rm max}.$ We use the notation
$T_{m}(q|\,\mathcal{K})$ to denote the integrated trispectrum as a function of
the soft mode, $q$, averaged over all hard modes
$k_{1},k_{2},k_{3},k_{4}\in\mathcal{K}$. We also analyze the collapsed
trispectrum divided into two hard mode bins, $\mathcal{K}_{a}$ and
$\mathcal{K}_{b}$. In this case, there are three configurations
$T_{m}(q|\,\mathcal{K}_{a}),$ $T_{m}(q|\,\mathcal{K}_{b})$, and
$T_{m}(q|\,\mathcal{K}_{a},\mathcal{K}_{b})$, where the final configuration
corresponds to the asymmetric configuration with
$k_{1},k_{2}\in\mathcal{K}_{a}$ and $k_{3},k_{4}\in\mathcal{K}_{b}$ (see Fig.
1).

To avoid binning effects due to the finite number of modes at low $q$ in our
simulations, we explicitly bin our theory model by computing $\langle
q^{-2+\Delta}\delta_{m}(\mathbf{q})\delta_{m}(-\mathbf{q})/T(q)\rangle$ and
$\langle q^{2}\delta_{m}(\mathbf{q})\delta_{m}(-\mathbf{q})\rangle$ using the
discrete density fields. This is described in detail in Goldstein  _et al._
[38].

####  III.2.2 Likelihood

We now describe how we fit our simulation measurements to constrain the
amplitude of PNG, $f_{\rm NL}^{\Delta}$. We focus on the likelihood for the
squeezed matter bispectrum in a single hard mode bin, which includes all of
the main ingredients without the notational complexity. We explain how to
generalize this to the other summary statistics considered in this work at the
end of this section.

Let ${\bm{\theta}}$ denote the vector of parameters that we vary in a fit,
_i.e._ , $f_{\rm NL}^{\Delta}$, $\mathrm{bar}{a}_{0}$, and
$\mathrm{bar}{a}_{2}$ for the bispectrum in a single bin. Since the main goal
of this work is to demonstrate that our non-perturbative models for the
squeezed matter bispectrum and collapsed trispectrum are unbiased, we often
use a joint likelihood in the soft matter power spectrum, $P_{m}(q)$, which
reduces cosmic variance. As described in Goldstein  _et al._ [38], the joint
likelihood for $B_{m}(q|\,\mathcal{K})$ and $P_{m}(q)$ can be written as a
likelihood in $B_{m}(q|\,\mathcal{K})$ with a parameter-dependent covariance
that depends on $P_{m}(q)$.212121In Sec. IV.2, we show constraints from a
likelihood in just $B_{m}(q|\,\mathcal{K})$. We determine these constraints
from the same procedure described in this section, _except_ we compute
$P_{m}(q)$ from linear theory and we do not use a parameter dependent
covariance.  Specifically,

$\displaystyle\mathcal{L}\propto\frac{1}{\sqrt{\det{\mathcal{C}}({\bm{\theta}})}}\bigg{[}1+\frac{\delta
B({\bm{\theta}})\cdot{\mathcal{C}}^{-1}({\bm{\theta}})\cdot\delta
B({\bm{\theta}})}{N_{\rm sim}-1}\bigg{]}^{-\frac{N_{\rm sim}}{2}}\,,$ (35)

where $\delta B({\bm{\theta}})\equiv{B}_{m}(q|\,\mathcal{K})-B_{\rm
thr.}(q,\mathcal{K}|{\bm{\theta}})$ is the difference between the measured
bispectrum and our theory model (Eq. 12),
${\mathcal{C}}_{ij}({\bm{\theta}})\equiv\langle\delta
B_{i}({\bm{\theta}})\delta B_{j}({\bm{\theta}})\rangle$ is the _parameter-
dependent_ covariance of $\delta B$, and $N_{\rm sim}$ is the number of
realizations used to estimate it. We estimate the covariance using $N_{\rm
sim}=50$ simulation realizations. To mitigate numerical instabilities, we
assume that the covariance is diagonal in $q$, which is a good approximation
over the scales we consider [39]. Finally, we use a multivariate
$t$-distribution to account for the fact that our covariance is estimated from
a finite number of simulations [124].

Figure 3: Impact of intermediate-mass scalars during inflation on the soft
limits of large-scale structure correlation functions. The solid points show
redshift $z=0$ measurements of the squeezed matter bispectrum (left),
collapsed matter trispectrum (middle), and scale-dependent halo bias (right)
from $N$-body simulations for different values of the amplitude ($f_{\rm
NL}^{\Delta}$) and particle mass ($\Delta$). To highlight the different
scalings, we divide all measurements by the power spectrum, $P_{m}(q)$, and
fix the ratios to one at $q=0.028~{}h/{\rm Mpc}$ for the bispectrum and
trispectrum and $q=0.021~{}h/{\rm Mpc}$ for the scale-dependent bias. The
squeezed primordial bispectrum in Eq. (6) leads to poles whose amplitude and
shape depend on $f_{\rm NL}^{\Delta}$ and $\Delta.$ The dashed lines show our
best-fit theory models given by Eq. (20), Eq. (II.4), and Eq. (26), for the
bispectrum, trispectrum, and the scale-dependent bias, respectively. The
shaded grey region indicates modes excluded from the halo analysis due to non-
linear biasing. See the main text for precise definitions of
$B_{m}(q|\,\mathcal{K})$ and $T_{m}(q|\,\mathcal{K}),$ as well as the scale
cuts and halo masses used in the analysis shown here.

Since sampling from a parameter-dependent covariance is costly, we evaluate
the covariance with the parameters fixed to their _maximum a posteriori_
values using the iterative maximization procedure described in Ref. [38]. In
Ref. [38], we found that this procedure has negligible impact on our results.
With the covariance determined, we sample from the parameter posterior using
the emcee package [125] assuming uniform wide priors with $-10^{5}<f_{\rm
NL}^{\Delta},\mathrm{bar}{a}_{i}<10^{5}$. We assess convergence using the
Gelman-Rubin statistic [126] with a tolerance of $|R-1|<0.01$.

To generalize the above discussion to a likelihood in the squeezed bispectrum
measured in $N_{k}$ hard mode bins, $\mathcal{K}_{1},\dots,\mathcal{K}_{n}$,
we replace the $N_{q}$-dimensional bispectrum data-vector with an $N_{q}\times
N_{k}$-dimensional data-vector of bispectrum measurements,
$B(q|\mathcal{K}_{i})$, in each hard mode bin. We estimate the covariance from
simulations assuming it is diagonal in $q$ but including off-diagonal
correlations between different $k$ bins. For each $k$ bin, we introduce a new
set of free coefficients $\mathrm{bar}{a}_{0}^{(i)}$ and
$\mathrm{bar}{a}_{2}^{(i)}$; thus we have $2N_{k}+1$ free parameters for the
bispectrum model. A similar approach was recently used in Ref. [39]. Here, we
analyze the squeezed bispectrum for $N_{k}=1,$ $2$, and $4$.

The collapsed trispectrum likelihood in a single hard mode bin takes the same
form as the single-bin bispectrum likelihood with $\delta
B({\bm{\theta}})\rightarrow\delta
T_{m}({\bm{\theta}})=T_{m}(q|\mathcal{K})-T_{\rm
thr.}(q,\mathcal{K}|{\bm{\theta}}).$ In addition to the free parameters in the
bispectrum model, the trispectrum theory model includes $\mathrm{bar}{T}_{0}$.
We vary $\mathrm{bar}{T}_{0}$ in our trispectrum analyses assuming a uniform
prior $5<\log_{10}\left(\mathrm{bar}{T}_{0}~{}({\rm Mpc}/h)^{9}\right)<15$.
For the simulation volume and scale cuts used in this work, we find that it is
sufficient to fix $\mathrm{bar}{a}_{2}=0$, while varying $f_{\rm
NL}^{\Delta},\mathrm{bar}{a}_{0}$, $\mathrm{bar}{T}_{0}$.222222If we vary
$\mathrm{bar}{a}_{2}$, then our constraints on $f_{\rm NL}^{\Delta}$ are still
unbiased. However, the error on $f_{\rm NL}^{\Delta}$ increases significantly
because of the degeneracy between $\mathrm{bar}{a}_{2}$ and
$\mathrm{bar}{T}_{0}$. Thus, we fix $\mathrm{bar}{a}_{2}=0$ for all
trispectrum analyses in this work, except for the joint analyses of the
squeezed bispectrum and collapsed trispectrum, as including the squeezed
bispectrum breaks the $\mathrm{bar}{T}_{0}-\mathrm{bar}{a}_{2}$ degeneracy.

We can generalize the collapsed trispectrum likelihood to include measurements
in $N_{k}$ hard mode bins. This procedure is identical to the generalization
to the multi-bin bispectrum likelihood, except the data vector also includes
the $\binom{N_{k}}{2}$ asymmetric trispectrum configurations,
$T(q|\mathcal{K}_{i},\mathcal{K}_{j}).$ Importantly,
$\log_{10}(\mathrm{bar}{T}_{i,j})$ is the only free parameter that needs to be
added to analyze the asymmetric configurations. Here, we analyze the collapsed
trispectrum with $N_{k}=1$ and $N_{k}=2$. Finally, we can jointly analyze the
squeezed matter bispectrum and collapsed trispectrum in multiple hard mode
bins by combining the multi-bin bispectrum and trispectrum likelihoods,
including the covariance between the squeezed bispectrum and collapsed
trispectrum at a fixed soft mode, $q.$ We fit the squeezed matter bispectrum
and collapsed trispectrum up to $q_{\rm max}=0.05~{}h/{\rm Mpc}.$

We study the scale-dependent halo bias using a joint likelihood in $P_{h}(q)$,
$P_{hm}(q)$, and $P_{m}(q)$ to reduce cosmic variance. We fix $f_{\rm
NL}^{\Delta}$ to its true value for a given simulation while varying the
linear bias ($b_{1}$), the non-Gaussian bias ($b_{\Psi,\Delta}$), and the
stochastic contribution ($\epsilon_{N}$). We assume the following uniform
priors: $0<b_{1}<5$, $-20<b_{\Psi,\Delta}<20$, and $-1<\epsilon_{N}<1$. We fit
the halo statistics up to $q_{\rm max}=0.025~{}h/{\rm Mpc}$, beyond which
higher-order bias terms are required given the high-mass sample and large
$f_{\rm NL}^{\rm loc.}$ value.

##  IV Results

Fig. 3 is one of the main results of this work: comparison of our squeezed
limit $N$-point function models to measurements from numerical simulations.
Using simulations with and without PNG, we plot measurements of the ratio of
the squeezed matter bispectrum (left), collapsed matter trispectrum (center),
and halo-matter cross-spectrum (right) relative to the soft matter power
spectrum $P_{m}(q)$ (all at redshift zero). For simulations without PNG, all
ratios approach a constant on large scales. On the other hand, for simulations
with PNG described by the squeezed bispectrum in Eq. (6), these ratios have
poles in the soft limits. The amplitude and shape of these poles are set by
$f_{\rm NL}^{\Delta}$ and $\Delta$, respectively.

The dashed lines show the best-fit theory models given by Eq. (20), Eq.
(II.4), and Eq. (26), for the squeezed matter bispectrum, collapsed matter
trispectrum, and the scale-dependent bias, respectively. For the matter
bispectrum and trispectrum, we determine the best-fit model parameters from a
joint analysis of the soft matter power spectrum, squeezed matter bispectrum,
and collapsed matter trispectrum with soft modes between
$0.006<q<0.05~{}h/{\rm Mpc}$ and hard modes averaged over a single bin with
$0.4<k<1.1~{}h/{\rm Mpc}$. For the scale-dependent halo bias, we determine the
best-fit model parameters from a joint analysis of the halo power spectrum,
halo-matter cross-spectrum, and matter power spectrum for halos with mass
between ${14.25}<\log_{10}(M_{\rm FoF})<14.5~{}M_{\odot}/h$. The scale-
dependent bias analysis excludes modes with $q>0.025~{}h/{\rm Mpc}$ because
the linear biasing model breaks down beyond these scales. The best-fit theory
predictions provide an excellent fit to the data for all summary statistics
and simulations considered, thus validating the models discussed above.

Figure 4: Non-linear contribution to the potential derivative (Eq. 16)
estimated from separate universe simulations for the three values of $\Delta$
analyzed in this work. The derivative is computed at redshift $z=0$. On large
scales, the potential derivative agrees with linear theory. On small scales,
the derivative is amplified by non-linear structure formation, particularly
for large values of $\Delta.$

###  IV.1 Validation of matter bispectrum and trispectrum models

To obtain constraints on $f_{\rm NL}^{\Delta}$ from the squeezed matter
bispectrum and collapsed matter trispectrum, we first need to estimate the
non-linear potential derivative in Eq. (16) from the separate universe
simulations. Fig. 4 shows the non-linear enhancement of the potential
derivative as a function of the wave-number for the three values of $\Delta$
analyzed in this work. The derivatives are computed numerically using the
method described in Sec. III.1.2. We average the potential derivative over
five realizations to minimize cosmic variance. On large scales, the ratio
$\partial\log P_{m}(k)/\partial\epsilon|_{\epsilon=0}/(2/k^{\Delta})$
approaches unity, as expected from linear theory. The small-scale potential
derivative is significantly enhanced due to non-linear structure formation,
particularly for large values of $\Delta.$ This suggests that the mode
coupling induced by the primordial bispectrum in Eq. (6) could significantly
impact small-scale structures for large values of $\Delta.$ However, a
detailed investigation of the shape of the potential derivative is beyond the
scope of this work.

Equipped with the separate universe potential derivatives and measurements of
the squeezed matter bispectrum and collapsed matter trispectrum from $N$-body
simulations, we turn to the main task: constraining $f_{\rm NL}^{\Delta}$.
Fig. 5 shows the mean value of $f_{\rm NL}^{\Delta}$ and its 68% confidence
interval as a function of the minimum hard momentum, $k_{\rm min}$, inferred
from the squeezed bispectrum (red) and collapsed trispectrum (blue) for the
simulations considered in this work. The horizontal dashed lines indicate the
true value of $f_{\rm NL}^{\Delta}$. For each value of $k_{\rm min},$ we
integrate over all hard modes in a narrow bin with $k_{\rm max}=k_{\rm
min}+0.06~{}h/{\rm Mpc}$. This allows us to robustly test the scale dependence
of the response functions estimated from the separate universe simulations.
Notably, any incorrect prediction of this scale dependence would significantly
bias the constraint on $f_{\rm NL}^{\Delta}$ (see Fig. 4). To assess potential
biases in our theoretical models, we use a joint likelihood with the soft
power spectrum $P_{m}(q)$, which minimizes the cosmic variance contribution
from the long-wavelength modes. Moreover, we fit the mean data-vector
estimated from $N_{\rm sim}=50$ simulations, rescaling the covariance by
$1/N_{\rm sim}.$

The top panels of Fig. 5 shows the constraints on $f_{\rm NL}^{\Delta}$
derived from the simulations with Gaussian initial conditions for different
values of $\Delta$. In all cases, we recover $f_{\rm NL}^{\Delta}=0$. The
constraints are slightly biased for $k_{\rm min}\approx 0.25~{}h/{\rm Mpc}$,
suggesting that these triangle configurations may not be sufficiently squeezed
for our bispectrum and trispectrum models to apply. Notice that the error on
$f_{\rm NL}^{\Delta}$ grows dramatically as $\Delta$ increases. This is
expected because $B_{m}(q|\mathcal{K})\propto f_{\rm
NL}^{\Delta}/q^{2-\Delta}$, meaning that for a fixed $f_{\rm NL}^{\Delta}$,
the signal becomes harder to detect at higher values of $\Delta$. Furthermore,
with larger $\Delta$, $f_{\rm NL}^{\Delta}$ becomes increasingly degenerate
with the nuisance parameter $\mathrm{bar}{a}_{0}$, resulting in an additional
penalty when marginalizing over $\mathrm{bar}{a}_{0}$. We discuss this in more
detail in Sec. IV.2.

Figure 5: Mean and 68% confidence interval on the inferred value of $f_{\rm
NL}^{\Delta}$ from a joint analysis of the squeezed matter bispectrum and
power spectrum (red), as well as the collapsed matter trispectrum and power
spectrum (blue) as a function of the _minimum_ hard mode, $k_{\rm min}.$ For
Gaussian initial conditions (top), we correctly infer $f_{\rm NL}^{\Delta}=0$
from fits to the bispectrum and trispectrum measured over a wide range of non-
linear modes. The error on $f_{\rm NL}^{\Delta}$ increases significantly for
larger values of $\Delta$. We also find unbiased constraints when analyzing
simulations with PNG described by the squeezed primordial bispectrum in Eq.
(6) (bottom). This demonstrates that the non-perturbative models for the
squeezed matter bispectrum and collapsed matter trispectrum derived in this
work can be used to constrain the Cosmological Collider deep into the non-
linear regime. Due to the non-Gaussian covariance, the constraints are highly
correlated across different $k$ bins.

The bottom panels of Fig. 5 presents the constraints on $f_{\rm NL}^{\Delta}$
derived from the simulations with PNG. In all cases, we correctly infer the
amplitude $f_{\rm NL}^{\Delta}$ up to scales of $2~{}h/{\rm Mpc}$.232323In
principle, we can go beyond $k=2~{}h/{\rm Mpc}$. However, given the numerical
resolution of our simulations, we choose to analyze only modes with
$k<2~{}h/{\rm Mpc}$. Using similar methods, Giri  _et al._ [39] constrained
$f_{\rm NL}^{\rm loc.}$ with the squeezed bispectrum and collapsed trispectrum
from the Quijote simulations up to $k=3~{}h/{\rm Mpc}.$ It is worth noting
that the constraints from the trispectrum for the $\Delta=1$ simulations are
slightly low at $k_{\rm min}\gtrsim 1.2~{}h/{\rm Mpc}$. This could be a
statistical fluctuation, or, alternatively, could arise since we fix
$\mathrm{bar}{a}_{2}=0$ in the trispectrum model thus neglecting higher-order
($q^{2}/k^{2}$) contributions. Finally, for the $\Delta=1$ simulations, the
loop contributions to the matter power spectrum are $\mathcal{O}(15\%)$, which
can lead to slight biases in the constraints.242424By eye, one may worry that
all of the constraints in Fig. 5 look suspiciously good. However, we stress
that due to the non-Gaussian squeezed bispectrum and collapsed trispectrum
covariances, the inferred value of $f_{\rm NL}^{\Delta}$ is highly correlated
as a function of $k_{\rm min}$. Furthermore, since all simulations are run
using the same seeds, the residuals will be correlated _across_ simulations.
Finally, rescaling the sample covariance by $1/N_{\rm sim}$ is only an
approximate measure of the covariance of the mean.

In conclusion, given the significant dependence of the potential derivative,
$\partial P_{m}(k|\epsilon,\Delta)/\partial\epsilon$, on both $\Delta$ and $k$
(see Fig. 4), the results of this section further validate our separate
universe models for the squeezed matter bispectrum and collapsed matter
trispectrum.

###  IV.2 Information content of the matter bispectrum and trispectrum

Having demonstrated that our non-perturbative models for the squeezed matter
bispectrum and collapsed matter trispectrum yield unbiased constraints on
$f_{\rm NL}^{\Delta}$ deep into the non-linear regime, we now assess the
information content of these summary statistics. All results in this section
use the mean bispectrum and trispectrum measurements averaged over 50
realizations and the covariance estimated for a single realization with a
$(1~{}{\rm Gpc}/h)^{3}$ volume. We first focus on the squeezed bispectrum.

Figure 6: Dependence of the 68% error on $f_{\rm NL}^{\Delta}$ on the maximum
hard mode, $k_{\rm max}$, and the number of $k$ bins for different values of
$\Delta$. For the single-bin analysis (black), increasing $k_{\rm max}$ yields
minimal improvement. In contrast, multi-bin analyses (red and blue) show
significant improvement at higher $k_{\rm max}$ due to sample variance
cancellation (see the main text for details). For larger $\Delta$,
$\sigma_{f_{\rm NL}^{\Delta}}$ is less sensitive to $k_{\rm max}$ because of
the additional $k^{-\Delta}$ scaling of the primordial squeezed bispectrum
(Eq. 6). All errors are estimated from a joint analysis of the redshift $z=0$
squeezed matter bispectrum and large-scale matter power spectrum, assuming
$k_{\rm min}=0.38~{}h/{\rm Mpc}$.

In Fig. 6, we plot the $1\sigma$ error on $f_{\rm NL}^{\Delta}$ as a function
of the maximum hard mode, $k_{\rm max}$, from a joint analysis of the squeezed
matter bispectrum and large-scale matter power spectrum for three different
values of $\Delta$. We analyze the squeezed bispectrum with soft modes between
$0.006<q<0.05~{}h/{\rm Mpc}$ and hard modes between $0.38<k<k_{\rm
max}~{}h/{\rm Mpc}$. We consider three values of $k_{\rm max}$ between 1 and
$2~{}h/{\rm Mpc}.$ For each value of $k_{\rm max}$, we show the error for an
analysis with one $k$ bin (black), two $k$ bins (red), and four $k$ bins
(blue). To emphasize how the constraints depend on the number of $k$ bins and
choice of $k_{\rm max}$, we divide all errors by the error from the single bin
analyses at $k_{\rm max}=1~{}h/{\rm Mpc}.$ The different panels correspond to
different values of $\Delta$. All errors are derived by fitting the
simulations with Gaussian initial conditions.

For the single bin analysis, the constraints improve only marginally when
increasing $k_{\rm max}$. Indeed, the most significant improvement is seen in
the local PNG analysis ($\Delta=0$), where the constraint on $f_{\rm
NL}^{\Delta}$ improves by 15% after increasing $k_{\rm max}$ from $1$ to
$2~{}h/{\rm Mpc}$. This is consistent with Refs. [38, 39], which found that
the constraints on $f_{\rm NL}^{\rm loc.}$ from the squeezed matter bispectrum
averaged over a single bin saturate quickly as a function of $k_{\rm max}$ due
to the non-Gaussian covariance [127, 128].

Figure 7: Marginalized posterior on $f_{\rm NL}^{\Delta}$ for a joint analysis
of the squeezed matter bispectrum and soft matter power spectrum with
different assumptions about the nuisance parameters, $a_{0}^{(i)}$ and
$a_{2}^{(i)}$, describing late-time physics. The error on $f_{\rm
NL}^{\Delta}$ is significantly larger for high values of $\Delta$. This is
primarily due to the degeneracy between $f_{\rm NL}^{\Delta}$ and
$\mathrm{bar}{a}_{0}^{(i)}$. For an _idealistic_ analysis with fixed
$a_{0}^{(i)}$, the error on $f_{\rm NL}^{\Delta}$ improves by a factor of 2,
3, and 5 for $\Delta=0,~{}0.5,$ and $1$, respectively. Fixing $a_{2}^{(i)}$
leads to further improvements. Table 2 lists the 68% confidence error bar for
these analyses. Marginalized error on $\bm{f_{\rm NL}^{\Delta}}$
---
$~{}~{}\Delta~{}~{}$ |  Fiducial |  Fixed $a_{0}^{(i)}$ |  Fixed $a_{0}^{(i)}$ and $a_{2}^{(i)}$
0 | 36 | 19 | 16
0.5 | 217 | 78 | 54
1 | 1300 | 243 | 128
Table 2: 68% confidence error bar on $f_{\rm NL}^{\Delta}$ from a joint
analysis of the squeezed matter bispectrum and matter power spectrum from a
$(1~{}{\rm Gpc}/h)^{3}$ volume with different assumptions about the nuisance
parameters, $a_{0}^{(i)}$ and $a_{2}^{(i)}$. The squeezed bispectrum is
measured with $0.006<q<0.05~{}h/{\rm Mpc}$ and averaged over two equal-width
hard mode bins with $0.38<k<2.01~{}h/{\rm Mpc}$. For large values of $\Delta$,
marginalizing over $a_{0}^{(i)}$ and $a_{2}^{(i)}$ significantly degrades the
constraint on $f_{\rm NL}^{\Delta}.$

When dividing the bispectrum into multiple $k$ bins (red and blue), the error
on $f_{\rm NL}^{\Delta}$ is notably more sensitive to $k_{\rm max}$. For
$\Delta=0$, we find an improvement of nearly 40% going from $k_{\rm max}=1$ to
$2~{}h/{\rm Mpc}$, consistent with previous works [39]. This improvement
arises because, as discussed extensively in Ref. [129], analyzing the
bispectrum in multiple $k$ bins cancels the cosmic variance associated with
the long mode in the bispectrum, $\delta_{m}(\mathbf{q})$. Physically, this
cosmic variance cancellation is analogous to the cancellation that occurs when
constraining $f_{\rm NL}^{\rm loc.}$ using the scale-dependent bias in a
multi-tracer survey [130, 131, 132, 133]. In our case, the locally measured
small-scale matter power spectrum in multiple different $k$-bins act as
multiple biased tracers of the long-wavelength mode. Note that we find
negligible further improvement if we use four bins instead of two.

From Fig. 6, it is evident that for larger values of $\Delta,$ the constraints
are less sensitive to the choice of $k_{\rm max}$. In linear theory, this is
expected since the bispectrum is suppressed by the $k^{-\Delta}$ scaling.
However, this argument does not necessarily hold in the non-linear regime,
where the $k$-dependence of the bispectrum is described by the potential
derivative, $\partial P_{m}(k|\epsilon,\Delta)/\partial\epsilon$. Indeed, Fig.
4 shows that while the non-linear contribution to the logarithmic derivative
decreases as a function of $k$ over $1<k<2~{}h/{\rm Mpc}$ for the $\Delta=0$
simulations, it _increases_ over these scales for the $\Delta=1$ simulations.
Nevertheless, this increase is insufficient to offset the steep $k^{-\Delta}$
contribution. Thus, the $\Delta=0.5$ and $\Delta=1$ squeezed bispectra are
still suppressed relative to the $\Delta=0$ simulations over the range of hard
modes that we analyze here.

Having analyzed the dependence of our constraints on $k_{\rm max}$ and the
number of hard mode bins, we now study how $\sigma_{f_{\rm NL}^{\Delta}}$
depends on the value of $\Delta.$ Fig. 7 shows the marginalized posterior on
$f_{\rm NL}^{\Delta}$ from a joint analysis of the squeezed matter bispectrum
and large-scale matter power spectrum for the three different values of
$\Delta$ considered in this work. We use the bispectrum estimated in two
equal-width bins with $k_{\rm min}=60k_{f}\approx 0.38~{}h/{\rm Mpc}$ and
$k_{\rm max}=320k_{f}\approx 2.01~{}h/{\rm Mpc}.$ We show constraints from our
fiducial analysis where we vary the coefficients, $\mathrm{bar}{a}_{0}^{(i)}$
and $\mathrm{bar}{a}_{2}^{(i)}$ (black), as well as two _idealistic_ analyses
where we fix $\mathrm{bar}{a}_{0}^{(i)}$ (blue) and
$\mathrm{bar}{a}_{0}^{(i)}$ and $\mathrm{bar}{a}_{2}^{(i)}$ (red), practically
assuming perfect knowledge of structure-formation physics. Table 2 lists the
68% confidence error bar for these analyses.

Figure 8: Constraint on $f_{\rm NL}^{\Delta}$ from the squeezed bispectrum
(blue) vs. a joint analysis of the squeezed bispectrum and collapsed
trispectrum (red). If the particular realization of the matter power spectrum
is known, the collapsed trispectrum does not improve the constraint on $f_{\rm
NL}^{\Delta}$ (solid lines). Conversely, the collapsed trispectrum
significantly improves constraints if the particular realization of $P_{m}(q)$
is unknown because it helps cancel cosmic variance.

In the fiducial analysis, the error on $f_{\rm NL}^{\Delta}$ is over 36 times
larger for $\Delta=1$ compared to $\Delta=0$. This significant difference in
constraining power arises from two main factors: (i) the milder scaling of the
squeezed-bispectrum pole for large values of $\Delta$, and (ii) the increased
degeneracy between $f_{\rm NL}^{\Delta}$ and our model’s nuisance parameters
for large values of $\Delta$. To assess how much information is lost upon
marginalizing over $\mathrm{bar}{a}_{0}^{(i)}$ and
$\mathrm{bar}{a}_{2}^{(i)}$, we repeat the fiducial analysis with
$\mathrm{bar}{a}_{0}^{(i)}$ and $\mathrm{bar}{a}_{2}^{(i)}$ fixed to their
_maximum a posteriori_ values, thus varying only $f_{\rm NL}^{\Delta}$. While
the error bar for the $\Delta=0$ analysis shrinks by a factor of two, the
error bar for the $\Delta=1$ analysis decreases by over a factor of ten. We
still find significant improvements in the constraints if we fix
$a_{0}^{(i)}$, but vary $a_{2}^{(i)}$. This demonstrates that marginalizing
over nuisance parameters, particularly $a_{0}^{(i)}$, significantly limits how
well we can constrain $f_{\rm NL}^{\Delta}$ using the method presented here.
The same situation arises when trying to constrain these models using the
scale-dependent bias, where marginalization over galaxy bias parameters
significantly degrades the constraints [80]. Therefore, to optimally constrain
$f_{\rm NL}^{\Delta}$ for large values of $\Delta$, it would be extremely
valuable to have informative priors for $\mathrm{bar}{a}_{0}^{(i)}$ and
$\mathrm{bar}{a}_{2}^{(i)}$. While this may be possible for matter using
theory [27, 99] or $N$-body simulations _e.g._ , [98, 97], it is significantly
more challenging for galaxies due to uncertainties in galaxy formation
(however, see Ref. [23] for recent progress on this front).

Figure 9: Constraints on $\tau_{\rm NL}^{\rm loc.}$ and $f_{\rm NL}^{\rm
loc.}$ from a joint analysis of the squeezed matter bispectrum and collapsed
matter trispectrum measured from the Quijote-PNG simulations, with $f_{\rm
NL}^{\rm loc.}=100$. The black line shows the relationship between $\tau_{\rm
NL}^{\rm loc.}$ and $f_{\rm NL}^{\rm loc.}$ for single-source local non-
Gaussianity. Using the collapsed trispectrum, we can jointly infer $f_{\rm
NL}^{\rm loc.}$ and $\tau_{\rm NL}^{\rm loc.}$. These results use the $z=0$
measurements with the covariance estimated for a $(1~{}{\rm Gpc}/h)^{3}$ box.

We now discuss the information from the collapsed trispectrum. In Fig. 8, we
show the error on $f_{\rm NL}^{\Delta}$ from a joint analysis of the squeezed
matter bispectrum and the collapsed matter trispectrum. We analyze the
bispectrum and trispectrum measurements binned in soft modes with
$0.005<q<0.05~{}h/{\rm Mpc}$, and split in two equal-width hard mode bins with
$k_{\rm min}=60k_{f}\approx 0.38~{}h/{\rm Mpc}$ and $k_{\rm
max}=320k_{f}\approx 2.01~{}h/{\rm Mpc}.$ We only show results for the
$\Delta=0.5$ simulations, but we find similar results for all simulations
considered in this work. The solid lines show the results from our fiducial
analysis where we assume a joint likelihood in the measured soft matter power
spectrum to cancel cosmic variance. For such an analysis, the constraint on
$f_{\rm NL}^{\Delta}$ improves negligibly after including the collapsed
trispectrum. As discussed in Ref. [129], analyzing the bispectrum in multiple
$k$-bins jointly with the measured matter power spectrum leads to near-perfect
cosmic variance cancellation. Thus, the collapsed matter trispectrum provides
little information on $f_{\rm NL}^{\Delta}$ beyond the squeezed matter
bispectrum and soft matter power spectrum. To further illustrate the
importance of sample variance cancellation in the matter sector, we repeat
this analysis assuming no knowledge of the particular realization of
$P_{m}(q).$ In this case, including the collapsed trispectrum (red dashed)
leads to considerably tighter constraints on $f_{\rm NL}^{\Delta}$ than an
analysis of the squeezed bispectrum alone (black dashed).

Finally, even if the collapsed trispectrum does not significantly improve the
constraint on $f_{\rm NL}^{\Delta}$ (if knowledge of the power spectrum is
available), it can be used to constrain more general models of PNG where
multiple fields contribute to the curvature perturbation. As an example of
this, Fig. 9 shows the marginalized constraints on $f_{\rm NL}^{\rm loc.}$ and
$\tau_{\rm NL}^{\rm loc.}$ derived from a joint analysis of measurements of
the squeezed matter bispectrum and collapsed trispectrum. This fit uses the
squeezed matter bispectrum and collapsed matter trispectrum measured in two
equal-width $k$-bins with $0.37<k<1.88~{}h/{\rm Mpc}$. We model the collapsed
trispectrum including the $\tau_{\rm NL}^{\rm loc.}$ contribution using Eq.
(23). We use the mean measurements from 50 realizations of the Quijote
simulations with $f_{\rm NL}^{\rm loc.}=100$ and the covariance estimated for
a $(1~{}{\rm Gpc}/h)^{3}$ volume. The solid black line shows the relationship
between $f_{\rm NL}^{\rm loc.}$ and $\tau_{\rm NL}^{\rm loc.}$ in quadratic
local non-Gaussianity. We infer the true value of $f_{\rm NL}^{\rm loc.}$ and
$\tau_{\rm NL}^{\rm loc.}.$ Thus, the collapsed trispectrum provides a
powerful probe of higher-order PNG.

###  IV.3 Scale-dependent halo bias

We now analyze the scale-dependent halo bias of our simulations. Fig. 10 shows
the $b_{1}-b_{\Psi,\Delta}$ relationship for mass-selected halos from the
$\Delta=0$ (top) and $\Delta=0.5$ (bottom) simulations at the redshift
$z=0$.252525We do not analyze the scale-dependent bias in the $\Delta=1$
simulations because of the large error on the values of $b_{\Psi,\Delta}$ fit
to these simulations. The black points denote the mean values of $b_{1}$ and
$b_{\Psi,\Delta}$, alongside their 68% confidence intervals, inferred from a
joint analysis of the large-scale halo auto-power spectrum ($P_{h}$(q)), halo-
matter cross-power spectrum $(P_{hm}(q))$, and matter power spectrum
$P_{m}(q)$ averaged over 50 realizations. For this analysis, we fix $f_{\rm
NL}^{\Delta}$ and $\Delta$ to their true values for a given simulation.
Finally, since we model the halo statistics using linear theory, we only
include modes with $0.006<q<0.025~{}h/{\rm Mpc}$.

Figure 10: Constraints on the $b_{1}-b_{\Psi,\Delta}$ relationship from PNG
simulations with $\Delta=0$ (top) and $\Delta=0.5$ (bottom). The black data
points are derived by fitting $P_{h}$, $P_{hm}$, and $P_{m}$ jointly in five
different halo mass bins with $f_{\rm NL}^{\Delta}$ and $\Delta$ fixed at
their true values. The solid red line, which is computed from separate
universe simulations using Eq. (27), is in excellent agreement with the
measured bias. The measured bias is also consistent with the universality
prediction (Eq. 28), although this prediction is likely inaccurate for more
realistic halo and galaxy samples [113, 114, 115, 116, 117]. At large values
of $b_{1}$, $b_{\Psi}^{\Delta=0.5}$ is significantly larger than
$b_{\Psi}^{\Delta=0}$, illustrating the impact of the Cosmological Collider
bispectrum on galaxy formation.

The top panel of Fig. 10 illustrates the $b_{1}-b_{\Psi,\Delta}$ relationship
for local PNG ($\Delta=0$). We find excellent agreement between the measured
non-Gaussian bias parameters and the predictions from the separate universe
simulations (solid red line) and universality (dotted black line). Note that
for $\Delta=0$, the universality prediction (Eq. 28) is equivalent to
$b_{\Psi,\Delta}=2\delta_{c}(b_{1}-1)$. The bottom panel shows the same
results for the $\Delta=0.5$ simulations. As in the $\Delta=0$ simulations,
the measured non-Gaussian bias is consistent with the theoretical predictions
from the separate universe simulations (solid red) and universality (solid
blue). Notably, for non-zero values of $\Delta$, the universality prediction
for $b_{\Psi,\Delta}$ depends on the Lagrangian radius of the halo,
$R_{*}(M)$. In this work, we estimate $R_{*}$ from the mean FoF mass of each
halo sample. Although we find that $b_{\Psi,\Delta}$ agrees with the
universality prediction for the halos considered here, we emphasize that the
main takeaway from Fig. 10 is the agreement between the separate universe
predictions and measured bias parameters. For more realistic halo and galaxy
samples, the universality prediction is likely a poor approximation of
$b_{\Psi,\Delta}$ [113, 114, 115, 116, 117].

In the bottom panel of Fig. 10, the lower dashed black line shows the
universality prediction for local PNG ($\Delta=0$). For large values of
$b_{1}$, the non-Gaussian bias for the $\Delta=0.5$ cosmology is significantly
enhanced relative to that of the $\Delta=0$ cosmology, akin to the non-linear
enhancement of the potential derivatives seen in Fig. 4. The upper dashed
black line in the same panel is a simple prediction for the non-Gaussian bias
based on dimensional analysis, where
$b_{\Psi,\Delta}=2\delta_{c}(b_{1}-1)R_{*}^{\Delta}$. This prediction
overestimates the non-Gaussian bias for the halos considered here.
Nevertheless, it captures the approximate shape of the $b_{1}-b_{\Psi,\Delta}$
relationship for this halo sample.

Stepping back, these findings demonstrate the utility of leveraging
simulations with analytic calculations to study the influence of intermediate-
mass particles during inflation on the large-scale halo bias. Using separate
universe simulations, we can accurately predict the non-Gaussian bias
parameters for specific halo samples. These results have several potential
applications, particularly in the context of generalizing studies of
$b_{\phi}$ in local PNG, to studies of $b_{\Psi,\Delta}$ for the models
considered here. Given the significant variation in the
$b_{\Psi,\Delta}(b_{1})$ relationship for the $\Delta=0.5$ simulations, one
could carefully choose galaxy samples to optimally constrain $f_{\rm
NL}^{\Delta}.$ Additionally, exploring higher-order statistics of the halo
field by generalizing the methods we developed for the matter field could be
interesting. We choose not to study the information content of the halo field
in this work because previous work has found that, for the number densities
used in the Quijote simulations, the squeezed halo bispectrum provides little
information beyond the scale-dependent halo bias [82, 134]. Nevertheless,
comparing the constraining power of the scale-dependent halo bias with other
summary statistics for the models considered here would be an interesting
follow-up.

##  V Conclusion

Fields with mass $m\approx H$ present during inflation produce distinct
signatures in the squeezed bispectrum of the primordial curvature perturbation
[41, 42, 43, 44]. Searching for these shapes of PNG presents a unique
opportunity to uncover the spectrum of particles present during inflation. In
this work, we have presented non-perturbative techniques to search for light
and intermediate-mass fields ($0\leq m<3H/2$) present during inflation using
LSS data. Our main results are summarized as follows:

  * •

We run $N$-body simulations with PNG sourced by intermediate-mass scalars
during inflation, using the same settings as Quijote-PNG [83, 82], thus
complementing existing $N$-body simulations with PNG. These simulations could
be a valuable tool for future studies of the Cosmological Collider.

  * •

Using the separate universe approach, we derive non-perturbative models for
the squeezed matter bispectrum (Eq. 17) and collapsed matter trispectrum (Eq.
22) associated with this type of PNG. We use these models to estimate the
amplitude of PNG, $f_{\rm NL}^{\Delta}$, from $N$-body simulations deep into
the non-linear regime ($k_{\rm max}=2~{}h/{\rm Mpc}$ at $z=0$).

  * •

For light fields ($\Delta=0$), the constraint on $f_{\rm NL}^{\Delta}$
improves by $\approx 40\%$ after increasing $k_{\rm max}$ from $1$ to
$2~{}h/{\rm Mpc}$. The improvement is smaller for heavier fields $(\Delta>0$)
because the bispectrum is suppressed by $\approx k^{-\Delta}$. Additionally,
including the collapsed trispectrum provides a negligible improvement in our
constraint on $f_{\rm NL}^{\Delta}$ (due to sample-variance cancellation in
the bispectrum), but allows us to constrain higher-order PNG, such as
$\tau_{\rm NL}^{\rm loc.}$.

  * •

The error on $f_{\rm NL}^{\Delta}$ grows significantly as $\Delta$ increases.
This is due to both the milder scaling of the pole, as well as the increased
degeneracy between $f_{\rm NL}^{\Delta}$ and the nuisance parameters in our
model for $\Delta>0.$

  * •

We study the scale-dependent halo-bias sourced by intermediate-mass particles
during inflation. The measured non-Gaussian bias parameter agrees with
theoretical predictions from (i) separate universe simulations and (ii)
assuming universal mass functions, although the second prediction will likely
fail for realistic galaxy samples.

There are several interesting follow-ups to this analysis. On the simulation
side, it would be interesting to extend our formalism to spinning fields,
which lead to an orientation-dependent bispectrum, as well as fields with
$m>3H/2$, which induce oscillations in the squeezed bispectrum. These squeezed
bispectrum configurations could be simulated by generalizing the coupling
kernel in Eq. (32). We discuss this in more detail in Appendix A. Previous
works have presented methods to search for spinning fields using galaxy shape
statistics [74, 79]. It would be interesting to quantify how competitive these
are with constraints from the squeezed bispectrum and collapsed trispectrum.
Similarly, the methods presented here could be used to search for tachyons
during inflation, which produce an even stronger scaling in the squeezed
bispectrum than local PNG [72].

On the observational front, we could apply these methods to the (CMB or
galaxy) weak lensing bispectrum, although our previous work has shown that
this does not provide competitive constraints on $f_{\rm NL}^{\rm loc.}$ for
realistic upcoming surveys [40]. The situation is expected be even worse for
large values of $\Delta$ due to the severe penalty from marginalizing over
nuisance parameters. Thus, it would be useful to have informative priors on
the nuisance parameters in this model. For the non-linear matter statistics
considered here, the main issue arises from marginalizing over
$\mathrm{bar}{a}_{0}^{(i)}$. In principle, this can be computed from separate
universe simulations [98, 97], or approximated from theory [92], although more
work is necessary to assess the robustness of these methods. It is also
possible that one can break these degeneracies by including non-squeezed
bispectrum configurations, although the methods presented here only work in
the squeezed limit, and the non-squeezed signatures are more easily mimicked
by single-field effects.

Given that this work focused primarily on constraining $f_{\rm NL}^{\Delta}$
using the squeezed matter bispectrum and collapsed matter trispectrum, it
would be interesting to extend our models for the matter field statistics to
describe the real-space and redshift-space squeezed galaxy bispectrum and
collapsed galaxy trispectrum. In particular, it would be useful to quantify
the information content of each of these statistics, as it remains to be seen
whether higher-order statistics of the halo field provide any significant
improvement over constraints on $f_{\rm NL}^{\Delta}$ from the scale-dependent
bias [134]. Alternatively, given the challenges posed by large-angle
systematics in galaxy surveys, a potentially interesting follow-up would be to
use our models to derive non-perturbative quadratic-estimator approaches for
constraining PNG, which can reconstruct the long-wavelength modes from small-
scale correlations [135, 136].

There are also several interesting follow-ups investigating the impact of the
Cosmological Collider on small-scale structure. For example, it could be
useful to have a more physical understanding of the behavior of the non-linear
response functions and non-Gaussian bias parameters estimated from our
separate universe simulations.262626It could be interesting to study the
impact of the Cosmological Collider bispectrum on the halo mass function and
halo density profiles. This would require higher-resolution simulations than
those used here, as well as a more careful procedure to generate the initial
conditions to account for non-squeezed bispectrum configurations. Such a study
would complement the recent findings of Ref. [137], which used simulations to
analyze the impact of scale-dependent PNG on non-linear structure. Although we
note that the squeezed limit of the bispectrum used in Ref. [137] is enhanced
by the hard mode $k$ (and independent of the soft mode $q$), whereas the
template considered here is suppressed by the ratio $(q/k)^{\Delta}$. A good
understanding of $b_{\Psi,\Delta}$ from, _e.g._ , hydrodynamical simulations,
could significantly improve our ability to constrain these models, especially
in multi-tracer scenarios [138]. On a similar note, the simulations developed
here could be used to test more formal aspects of the non-Gaussian bias
parameters, such as their running [139].

Finally, the ultimate dream of Cosmological Collider Physics is not to
constrain the amplitude, $f_{\rm NL}^{\Delta}$, but rather to determine the
shape of the pole, $\Delta$, and consequently the particle mass. Ideally, this
should be done as a joint analysis with other shapes of $f_{\rm NL}$, which
parameterize the background of inflationary particles. Unfortunately, as there
is currently no detection of PNG, there is no background to marginalize over.
Thus, just like in true collider physics, we are left anxiously waiting.

###### Acknowledgements.

We are grateful to Fabian Schmidt for insightful discussions and to Paco
Villaescusa-Navarro for help running the Quijote simulations. OHEP is a Junior
Fellow of the Simons Society of Fellows and acknowledges Aperol for providing
creative inspiration. JCH acknowledges support from NSF grant AST-2108536,
NASA grant 80NSSC22K0721, NASA grant 80NSSC23K0463, DOE grant DE-SC0011941,
the Sloan Foundation, and the Simons Foundation. LH acknowledges support from
the DOE DE-SC011941 and a Simons Fellowship in Theoretical Physics. The
authors acknowledge the Texas Advanced Computing Center
(TACC)272727http://www.tacc.utexas.edu at The University of Texas at Austin
for providing computational resources that have contributed to the research
results reported within this paper. We acknowledge computing resources from
Columbia University’s Shared Research Computing Facility project, which is
supported by NIH Research Facility Improvement Grant 1G20RR030893-01, and
associated funds from the New York State Empire State Development, Division of
Science Technology and Innovation (NYSTAR) Contract C090171, both awarded
April 15, 2010.

## Data Availability

The simulations used in this work are publicly available. More details can be
found at https://quijote-simulations.readthedocs.io/en/latest/collider.html.

## References

  * Guth [1981] A. H. Guth, Phys. Rev. D 23, 347 (1981).
  * Linde [1982] A. D. Linde, Phys. Lett. B 108, 389 (1982).
  * Albrecht and Steinhardt [1982] A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48, 1220 (1982).
  * 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 .
  * Creminelli  _et al._ [2011] P. Creminelli, G. D’Amico, M. Musso, and J. Norena, JCAP 11, 038 (2011), arXiv:1106.1462 [astro-ph.CO] .
  * Pajer  _et al._ [2013] E. Pajer, F. Schmidt, and M. Zaldarriaga, Phys. Rev. D 88, 083502 (2013), arXiv:1305.0824 [astro-ph.CO] .
  * Cabass  _et al._ [2017] G. Cabass, E. Pajer, and F. Schmidt, JCAP 01, 003 (2017), arXiv:1612.00033 [hep-th] .
  * Meerburg  _et al._ [2019] P. D. Meerburg _et al._ , (2019), arXiv:1903.04409 [astro-ph.CO] .
  * Achúcarro  _et al._ [2022] A. Achúcarro _et al._ , (2022), arXiv:2203.08128 [astro-ph.CO] .
  * Levi  _et al._ [2013] M. Levi _et al._ (DESI), (2013), arXiv:1308.0847 [astro-ph.CO] .
  * Laureijs  _et al._ [2011] R. Laureijs _et al._ (EUCLID), (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Mandelbaum  _et al._ [2018] R. Mandelbaum _et al._ (LSST Dark Energy Science), (2018), arXiv:1809.01669 [astro-ph.CO] .
  * Doré  _et al._ [2014] O. Doré _et al._ (SPHEREx), (2014), arXiv:1412.4872 [astro-ph.CO] .
  * Ade  _et al._ [2019] P. Ade _et al._ (Simons Observatory), JCAP 02, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Abazajian  _et al._ [2016] K. N. Abazajian _et al._ (CMB-S4), (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Akrami  _et al._ [2020] Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A9 (2020), arXiv:1905.05697 [astro-ph.CO] .
  * Sailer  _et al._ [2021] N. Sailer, E. Castorina, S. Ferraro, and M. White, JCAP 12, 049 (2021), arXiv:2106.09713 [astro-ph.CO] .
  * Cabass  _et al._ [2023] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonovic, and M. Zaldarriaga, Phys. Lett. B 841, 137912 (2023), arXiv:2211.14899 [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] .
  * Cabass  _et al._ [2022a] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, Phys. Rev. D 106, 043506 (2022a), arXiv:2204.01781 [astro-ph.CO] .
  * Cabass  _et al._ [2022b] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, Phys. Rev. Lett. 129, 021301 (2022b), arXiv:2201.07238 [astro-ph.CO] .
  * Ivanov  _et al._ [2024] M. M. Ivanov, C. Cuesta-Lazaro, S. Mishra-Sharma, A. Obuljen, and M. W. Toomey, (2024), arXiv:2402.13310 [astro-ph.CO] .
  * Cabass  _et al._ [2024] G. Cabass, O. H. E. Philcox, M. M. Ivanov, K. Akitsu, S.-F. Chen, M. Simonović, and M. Zaldarriaga, (2024), arXiv:2404.01894 [astro-ph.CO] .
  * Kehagias and Riotto [2013] A. Kehagias and A. Riotto, Nucl. Phys. B 873, 514 (2013), arXiv:1302.0130 [astro-ph.CO] .
  * Peloso and Pietroni [2013] M. Peloso and M. Pietroni, JCAP 05, 031 (2013), arXiv:1302.0223 [astro-ph.CO] .
  * Valageas [2014] P. Valageas, Phys. Rev. D 89, 123522 (2014), arXiv:1311.4286 [astro-ph.CO] .
  * Creminelli  _et al._ [2014a] P. Creminelli, J. Gleyzes, M. Simonović, and F. Vernizzi, JCAP 02, 051 (2014a), arXiv:1311.0290 [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] .
  * Horn  _et al._ [2014] B. Horn, L. Hui, and X. Xiao, JCAP 09, 044 (2014), arXiv:1406.0842 [hep-th] .
  * Horn  _et al._ [2015] B. Horn, L. Hui, and X. Xiao, JCAP 09, 068 (2015), arXiv:1502.06980 [hep-th] .
  * Marinucci  _et al._ [2020] M. Marinucci, T. Nishimichi, and M. Pietroni, JCAP 07, 054 (2020), arXiv:2005.09574 [astro-ph.CO] .
  * Esposito  _et al._ [2019] A. Esposito, L. Hui, and R. Scoccimarro, Phys. Rev. D 100, 043536 (2019), arXiv:1905.11423 [astro-ph.CO] .
  * Kehagias  _et al._ [2014] A. Kehagias, J. Noreña, H. Perrier, and A. Riotto, Nucl. Phys. B 883, 83 (2014), arXiv:1311.0786 [astro-ph.CO] .
  * Simonović [2014] M. Simonović,  _Cosmological Consistency Relations_ , Ph.D. thesis, SISSA, Trieste (2014).
  * Kehagias  _et al._ [2015] A. Kehagias, A. Moradinezhad Dizgah, J. Noreña, H. Perrier, and A. Riotto, JCAP 08, 018 (2015), arXiv:1503.04467 [astro-ph.CO] .
  * Creminelli  _et al._ [2014b] P. Creminelli, J. Gleyzes, L. Hui, M. Simonović, and F. Vernizzi, JCAP 06, 009 (2014b), arXiv:1312.6074 [astro-ph.CO] .
  * Goldstein  _et al._ [2022] S. Goldstein, A. Esposito, O. H. E. Philcox, L. Hui, J. C. Hill, R. Scoccimarro, and M. H. Abitbol, Phys. Rev. D 106, 123525 (2022), arXiv:2209.06228 [astro-ph.CO] .
  * Giri  _et al._ [2023] U. Giri, M. Münchmeyer, and K. M. Smith, (2023), arXiv:2305.03070 [astro-ph.CO] .
  * Goldstein  _et al._ [2024] S. Goldstein, O. H. E. Philcox, J. C. Hill, A. Esposito, and L. Hui, Phys. Rev. D 109, 043515 (2024), arXiv:2310.12959 [astro-ph.CO] .
  * Arkani-Hamed and Maldacena [2015] N. Arkani-Hamed and J. Maldacena, (2015), arXiv:1503.08043 [hep-th] .
  * Assassi  _et al._ [2012] V. Assassi, D. Baumann, and D. Green, JCAP 11, 047 (2012), arXiv:1204.4207 [hep-th] .
  * Chen and Wang [2010a] X. Chen and Y. Wang, JCAP 04, 027 (2010a), arXiv:0911.3380 [hep-th] .
  * Lee  _et al._ [2016] H. Lee, D. Baumann, and G. L. Pimentel, JHEP 12, 040 (2016), arXiv:1607.03735 [hep-th] .
  * Chen and Wang [2010b] X. Chen and Y. Wang, Phys. Rev. D 81, 063511 (2010b), arXiv:0909.0496 [astro-ph.CO] .
  * Baumann and Green [2012] D. Baumann and D. Green, Phys. Rev. D 85, 103520 (2012), arXiv:1109.0292 [hep-th] .
  * Chen and Wang [2012] X. Chen and Y. Wang, JCAP 09, 021 (2012), arXiv:1205.0160 [hep-th] .
  * Noumi  _et al._ [2013] T. Noumi, M. Yamaguchi, and D. Yokoyama, JHEP 06, 051 (2013), arXiv:1211.1624 [hep-th] .
  * Ghosh  _et al._ [2014] A. Ghosh, N. Kundu, S. Raju, and S. P. Trivedi, JHEP 07, 011 (2014), arXiv:1401.1426 [hep-th] .
  * Dimastrogiovanni  _et al._ [2016] E. Dimastrogiovanni, M. Fasiello, and M. Kamionkowski, JCAP 02, 017 (2016), arXiv:1504.05993 [astro-ph.CO] .
  * Kehagias and Riotto [2017] A. Kehagias and A. Riotto, JCAP 07, 046 (2017), arXiv:1705.05834 [hep-th] .
  * An  _et al._ [2018] H. An, M. McAneny, A. K. Ridgway, and M. B. Wise, JHEP 06, 105 (2018), arXiv:1706.09971 [hep-ph] .
  * Kumar and Sundrum [2018] S. Kumar and R. Sundrum, JHEP 05, 011 (2018), arXiv:1711.03988 [hep-ph] .
  * Baumann  _et al._ [2018] D. Baumann, G. Goon, H. Lee, and G. L. Pimentel, JHEP 04, 140 (2018), arXiv:1712.06624 [hep-th] .
  * Bordin  _et al._ [2018] L. Bordin, P. Creminelli, A. Khmelnitsky, and L. Senatore, JCAP 10, 013 (2018), arXiv:1806.10587 [hep-th] .
  * Kumar and Sundrum [2019] S. Kumar and R. Sundrum, JHEP 04, 120 (2019), arXiv:1811.11200 [hep-ph] .
  * Goon  _et al._ [2019] G. Goon, K. Hinterbichler, A. Joyce, and M. Trodden, JHEP 10, 182 (2019), arXiv:1812.07571 [hep-th] .
  * Hook  _et al._ [2020] A. Hook, J. Huang, and D. Racco, JHEP 01, 105 (2020), arXiv:1907.10624 [hep-ph] .
  * Kumar and Sundrum [2020] S. Kumar and R. Sundrum, JHEP 04, 077 (2020), arXiv:1908.11378 [hep-ph] .
  * Liu  _et al._ [2020] T. Liu, X. Tong, Y. Wang, and Z.-Z. Xianyu, JHEP 04, 189 (2020), arXiv:1909.01819 [hep-ph] .
  * Wang and Xianyu [2020] L.-T. Wang and Z.-Z. Xianyu, JHEP 02, 044 (2020), arXiv:1910.12876 [hep-ph] .
  * Bodas  _et al._ [2021] A. Bodas, S. Kumar, and R. Sundrum, JHEP 02, 079 (2021), arXiv:2010.04727 [hep-ph] .
  * Lu  _et al._ [2021] Q. Lu, M. Reece, and Z.-Z. Xianyu, JHEP 12, 098 (2021), arXiv:2108.11385 [hep-ph] .
  * Pinol  _et al._ [2023] L. Pinol, S. Aoki, S. Renaux-Petel, and M. Yamaguchi, Phys. Rev. D 107, L021301 (2023), arXiv:2112.05710 [hep-th] .
  * Cui and Xianyu [2022] Y. Cui and Z.-Z. Xianyu, Phys. Rev. Lett. 129, 111301 (2022), arXiv:2112.10793 [hep-ph] .
  * Reece  _et al._ [2023] M. Reece, L.-T. Wang, and Z.-Z. Xianyu, Phys. Rev. D 107, L101304 (2023), arXiv:2204.11869 [hep-ph] .
  * Chen  _et al._ [2022] X. Chen, R. Ebadi, and S. Kumar, JCAP 08, 083 (2022), arXiv:2205.01107 [hep-ph] .
  * Chen  _et al._ [2023] X. Chen, J. Fan, and L. Li, JHEP 12, 197 (2023), arXiv:2303.03406 [hep-ph] .
  * Jazayeri  _et al._ [2023] S. Jazayeri, S. Renaux-Petel, and D. Werth, JCAP 12, 035 (2023), arXiv:2307.01751 [hep-th] .
  * Chakraborty and Stout [2024a] P. Chakraborty and J. Stout, JHEP 02, 021 (2024a), arXiv:2310.01494 [hep-th] .
  * Chakraborty and Stout [2024b] P. Chakraborty and J. Stout, JHEP 03, 149 (2024b), arXiv:2311.09219 [hep-th] .
  * McCulloch  _et al._ [2024] C. McCulloch, E. Pajer, and X. Tong, JHEP 05, 262 (2024), arXiv:2401.11009 [hep-th] .
  * Quintin  _et al._ [2024] J. Quintin, X. Chen, and R. Ebadi, (2024), arXiv:2405.11016 [astro-ph.CO] .
  * Schmidt  _et al._ [2015] F. Schmidt, N. E. Chisari, and C. Dvorkin, JCAP 10, 032 (2015), arXiv:1506.02671 [astro-ph.CO] .
  * Gleyzes  _et al._ [2017] J. Gleyzes, R. de Putter, D. Green, and O. Doré, JCAP 04, 002 (2017), arXiv:1612.06366 [astro-ph.CO] .
  * Moradinezhad Dizgah and Dvorkin [2018] A. Moradinezhad Dizgah and C. Dvorkin, JCAP 01, 010 (2018), arXiv:1708.06473 [astro-ph.CO] .
  * Cabass  _et al._ [2018] G. Cabass, E. Pajer, and F. Schmidt, JCAP 09, 003 (2018), arXiv:1804.07295 [astro-ph.CO] .
  * Moradinezhad Dizgah  _et al._ [2018] A. Moradinezhad Dizgah, H. Lee, J. B. Muñoz, and C. Dvorkin, JCAP 05, 013 (2018), arXiv:1801.07265 [astro-ph.CO] .
  * Akitsu  _et al._ [2021] K. Akitsu, T. Kurita, T. Nishimichi, M. Takada, and S. Tanaka, Phys. Rev. D 103, 083508 (2021), arXiv:2007.03670 [astro-ph.CO] .
  * Green  _et al._ [2023] D. Green, Y. Guo, J. Han, and B. Wallisch, (2023), arXiv:2311.04882 [astro-ph.CO] .
  * Sohn  _et al._ [2024] W. Sohn, D.-G. Wang, J. R. Fergusson, and E. P. S. Shellard, (2024), arXiv:2404.07203 [astro-ph.CO] .
  * Coulton  _et al._ [2022] W. R. Coulton, F. Villaescusa-Navarro, D. Jamieson, M. Baldi, G. Jung, D. Karagiannis, M. Liguori, L. Verde, and B. D. Wandelt, (2022), arXiv:2206.15450 [astro-ph.CO] .
  * Villaescusa-Navarro  _et al._ [2020] F. Villaescusa-Navarro _et al._ , Astrophys. J. Suppl. 250, 2 (2020), arXiv:1909.05273 [astro-ph.CO] .
  * Komatsu and Spergel [2001] E. Komatsu and D. N. Spergel, Phys. Rev. D 63, 063002 (2001), arXiv:astro-ph/0005036 .
  * Sefusatti  _et al._ [2012] E. Sefusatti, J. R. Fergusson, X. Chen, and E. P. S. Shellard, JCAP 08, 033 (2012), arXiv:1204.6318 [astro-ph.CO] .
  * Cyr-Racine and Schmidt [2011] F.-Y. Cyr-Racine and F. Schmidt, Phys. Rev. D 84, 083505 (2011), arXiv:1106.2806 [astro-ph.CO] .
  * Schmidt and Kamionkowski [2010] F. Schmidt and M. Kamionkowski, Phys. Rev. D 82, 103002 (2010), arXiv:1008.0638 [astro-ph.CO] .
  * Scoccimarro  _et al._ [2012] R. Scoccimarro, L. Hui, M. Manera, and K. C. Chan, Phys. Rev. D 85, 083002 (2012), arXiv:1108.5512 [astro-ph.CO] .
  * Schmidt  _et al._ [2013] F. Schmidt, D. Jeong, and V. Desjacques, Phys. Rev. D 88, 023515 (2013), arXiv:1212.0868 [astro-ph.CO] .
  * Desjacques  _et al._ [2018] V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rept. 733, 1 (2018), arXiv:1611.09787 [astro-ph.CO] .
  * Wagner  _et al._ [2015a] C. Wagner, F. Schmidt, C.-T. Chiang, and E. Komatsu, Mon. Not. Roy. Astron. Soc. 448, L11 (2015a), arXiv:1409.6294 [astro-ph.CO] .
  * Valageas  _et al._ [2017] P. Valageas, A. Taruya, and T. Nishimichi, Phys. Rev. D 95, 023504 (2017), arXiv:1610.00993 [astro-ph.CO] .
  * Chiang  _et al._ [2014] C.-T. Chiang, C. Wagner, F. Schmidt, and E. Komatsu, JCAP 05, 048 (2014), arXiv:1403.3411 [astro-ph.CO] .
  * Wagner  _et al._ [2015b] C. Wagner, F. Schmidt, C.-T. Chiang, and E. Komatsu, JCAP 08, 042 (2015b), arXiv:1503.03487 [astro-ph.CO] .
  * Chiang  _et al._ [2017] C.-T. Chiang, A. M. Cieplak, F. Schmidt, and A. Slosar, JCAP 06, 022 (2017), arXiv:1701.03375 [astro-ph.CO] .
  * Barreira and Schmidt [2017] A. Barreira and F. Schmidt, JCAP 06, 053 (2017), arXiv:1703.09212 [astro-ph.CO] .
  * Biagetti  _et al._ [2022a] M. Biagetti, J. Calles, L. Castiblanco, K. González, and J. Noreña, (2022a), arXiv:2212.11940 [astro-ph.CO] .
  * Chiang [2017] C.-T. Chiang, Phys. Rev. D 95, 123517 (2017), arXiv:1701.03374 [astro-ph.CO] .
  * Nishimichi and Valageas [2014] T. Nishimichi and P. Valageas, Phys. Rev. D 90, 023546 (2014), arXiv:1402.3293 [astro-ph.CO] .
  * Lewis [2011] A. Lewis, JCAP 10, 026 (2011), arXiv:1107.5431 [astro-ph.CO] .
  * Kenton and Mulryne [2016] Z. Kenton and D. J. Mulryne, JCAP 10, 035 (2016), arXiv:1605.03435 [astro-ph.CO] .
  * Suyama and Yamaguchi [2008] T. Suyama and M. Yamaguchi, Phys. Rev. D 77, 023505 (2008), arXiv:0709.2545 [astro-ph] .
  * Sugiyama  _et al._ [2011] N. S. Sugiyama, E. Komatsu, and T. Futamase, Phys. Rev. Lett. 106, 251301 (2011), arXiv:1101.3636 [gr-qc] .
  * Smith  _et al._ [2011] K. M. Smith, M. Loverde, and M. Zaldarriaga, Phys. Rev. Lett. 107, 191301 (2011), arXiv:1108.1805 [astro-ph.CO] .
  * Baumann  _et al._ [2013] D. Baumann, S. Ferraro, D. Green, and K. M. Smith, JCAP 05, 001 (2013), arXiv:1209.2173 [astro-ph.CO] .
  * 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] .
  * Desjacques  _et al._ [2009] V. Desjacques, U. Seljak, and I. Iliev, Mon. Not. Roy. Astron. Soc. 396, 85 (2009), arXiv:0811.2748 [astro-ph] .
  * Desjacques  _et al._ [2011a] V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rev. D 84, 061301 (2011a).
  * Desjacques  _et al._ [2011b] V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rev. D 84, 063512 (2011b).
  * Reid  _et al._ [2010] B. Reid, L. Verde, K. Dolag, S. Matarrese, and L. Moscardini, JCAP 2010, 013 (2010), arXiv:1004.1637 [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 [2022a] A. Barreira, JCAP 01, 033 (2022a), arXiv:2107.06887 [astro-ph.CO] .
  * Lazeyras  _et al._ [2023] T. Lazeyras, A. Barreira, F. Schmidt, and V. Desjacques, JCAP 01, 023 (2023), arXiv:2209.07251 [astro-ph.CO] .
  * Barreira [2022b] A. Barreira, (2022b), arXiv:2205.05673 [astro-ph.CO] .
  * Barreira and Krause [2023] A. Barreira and E. Krause, (2023), arXiv:2302.09066 [astro-ph.CO] .
  * Hadzhiyska  _et al._ [2024] B. Hadzhiyska, L. H. Garrison, D. J. Eisenstein, and S. Ferraro, Phys. Rev. D 109, 103530 (2024), arXiv:2402.10881 [astro-ph.CO] .
  * Springel [2005] V. Springel, Mon. Not. Roy. Astron. Soc. 364, 1105 (2005), arXiv:astro-ph/0505010 .
  * Lewis  _et al._ [2000] A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), arXiv:astro-ph/9911177 .
  * Davis  _et al._ [1985] M. Davis, G. Efstathiou, C. S. Frenk, and S. D. M. White, Astrophys. J. 292, 371 (1985).
  * Regan  _et al._ [2012] D. M. Regan, M. M. Schmittfull, E. P. S. Shellard, and J. R. Fergusson, Phys. Rev. D 86, 123524 (2012), arXiv:1108.3813 [astro-ph.CO] .
  * Villaescusa-Navarro [2018] F. Villaescusa-Navarro, “Pylians: Python libraries for the analysis of numerical simulations,” Astrophysics Source Code Library, record ascl:1811.008 (2018), ascl:1811.008 .
  * Sellentin and Heavens [2016] E. Sellentin and A. F. Heavens, Mon. Not. Roy. Astron. Soc. 456, L132 (2016), arXiv:1511.05969 [astro-ph.CO] .
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. Pac. 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Gelman and Rubin [1992] A. Gelman and D. B. Rubin, Statist. Sci. 7, 457 (1992).
  * Biagetti  _et al._ [2022b] M. Biagetti, L. Castiblanco, J. Noreña, and E. Sefusatti, JCAP 09, 009 (2022b), arXiv:2111.05887 [astro-ph.CO] .
  * Flöss  _et al._ [2022] T. Flöss, M. Biagetti, and P. D. Meerburg, (2022), arXiv:2206.10458 [astro-ph.CO] .
  * de Putter [2018] R. de Putter, (2018), arXiv:1802.06762 [astro-ph.CO] .
  * Seljak [2009] U. Seljak, Phys. Rev. Lett. 102, 021302 (2009), arXiv:0807.1770 [astro-ph] .
  * McDonald and Seljak [2009] P. McDonald and U. Seljak, JCAP 10, 007 (2009), arXiv:0810.0323 [astro-ph] .
  * Hamaus  _et al._ [2012] N. Hamaus, U. Seljak, and V. Desjacques, Phys. Rev. D 86, 103513 (2012), arXiv:1207.1102 [astro-ph.CO] .
  * Ferraro and Smith [2015] S. Ferraro and K. M. Smith, Phys. Rev. D 91, 043506 (2015), arXiv:1408.3126 [astro-ph.CO] .
  * Giri  _et al._ [2022] U. Giri, M. Münchmeyer, and K. M. Smith, (2022), arXiv:2205.12964 [astro-ph.CO] .
  * Darwish  _et al._ [2021] O. Darwish, S. Foreman, M. M. Abidi, T. Baldauf, B. D. Sherwin, and P. D. Meerburg, Phys. Rev. D 104, 123520 (2021), arXiv:2007.08472 [astro-ph.CO] .
  * Wang and Jeong [2024] Z. Wang and D. Jeong, JCAP 07, 020 (2024), arXiv:2312.17321 [astro-ph.CO] .
  * Baldi  _et al._ [2024] M. Baldi _et al._ , (2024), arXiv:2407.06641 [astro-ph.CO] .
  * Sullivan  _et al._ [2023] J. M. Sullivan, T. Prijon, and U. Seljak, JCAP 08, 004 (2023), arXiv:2303.08901 [astro-ph.CO] .
  * Nikolis  _et al._ [2024] C. Nikolis, H. Rubira, and F. Schmidt, (2024), arXiv:2405.21002 [astro-ph.CO] .
  * Higuchi [1987] A. Higuchi, Nucl. Phys. B 282, 397 (1987).
  * Coulton  _et al._ [2023] W. R. Coulton, O. H. E. Philcox, and F. Villaescusa-Navarro, (2023), arXiv:2306.11782 [astro-ph.CO] .
  * Gualdi  _et al._ [2021] D. Gualdi, S. Novell-Masot, H. Gil-Marin, and L. Verde, JCAP 01, 015 (2021), arXiv:2009.02290 [astro-ph.CO] .
  * Smith  _et al._ [2015] K. M. Smith, L. Senatore, and M. Zaldarriaga, (2015), arXiv:1502.00635 [astro-ph.CO] .
  * Shen  _et al._ [2024] D. Shen, E. Schaan, and S. Ferraro, (2024), arXiv:2402.04309 [astro-ph.CO] .

##  Appendix A Cosmological Collider initial conditions including spin and
massive-er particles

In this appendix, we describe how to generate $N$-body initial conditions with
the Cosmological Collider squeezed bispectrum template for fields with spin,
as well as masses $m>3/2H$. Importantly, the transformations derived here can
be efficiently computed using Fast Fourier Transformations (FFTs).

###  A.1 Intermediate-mass spinning particles

For intermediate-mass fields with even spin, $s$,282828The squeezed bispectrum
suppressed by an additional power of $(q/k)$ for odd values of $s$ [76, 78]
and thus ignored herein. the squeezed bispectrum depends on the angle between
the momenta $\mathbf{q}$ and $\mathbf{k}$,

$\lim_{q\ll k}B_{\Phi}(\mathbf{q},\mathbf{k})=4f_{\rm
NL}^{\Delta}\left(\frac{q}{k}\right)^{\Delta}\mathcal{L}_{s}(\hat{\mathbf{q}}\cdot\hat{\mathbf{k}})P_{\Phi}(q)P_{\Phi}(k),$
(36)

where $\mathcal{L}_{s}(x)$ is the Legendre polynomial of order $s$ and the
(real) exponent is,292929Unless $\Delta=3/2$, the trispectrum in such models
takes a more complex form than that of Eq. (II.4) [41], due to the massive
field polarization states. Treatment of this is beyond the scope of this
paper.

$\Delta=\begin{cases}3/2-\sqrt{9/4-m^{2}/H^{2}}&\textrm{ for }s=0,\\\
3/2-\sqrt{\left(s-{1}/{2}\right)^{2}-m^{2}/H^{2}}&\textrm{ for }s\neq
0.\end{cases}$ (37)

Following Sec. III.1.1, we can simulate this squeezed bispectrum using the
kernel

$\displaystyle\begin{split}K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})&=\left[\left(\frac{k_{1}}{k_{12}}\right)^{\Delta}+\left(\frac{k_{2}}{k_{12}}\right)^{\Delta}-1\right]\mathcal{L}_{s}(\hat{\mathbf{k}_{1}}\cdot\hat{\mathbf{k}_{2}})\,,\\\
&=\left(\frac{4\pi}{s+1}\right)\sum_{m=-s}^{s}\left[\left(\frac{k_{1}}{k_{12}}\right)^{\Delta}+\left(\frac{k_{2}}{k_{12}}\right)^{\Delta}-1\right]Y_{sm}(\hat{\mathbf{k}_{1}})\,Y_{sm}^{*}(\hat{\mathbf{k}_{2}})\,,\end{split}$
(38)

where we separated the kernel using the spherical harmonic addition
theorem.303030For a fixed value of $s$, one could also derive a coupling
kernel by writing out the Legendre polynomial of order $s$ explicitly instead
of using spherical harmonics. Ref. [79] does this for massless spin-2 fields.
The auxiliary field in Eq. (31) can be efficiently evaluated using FFTs

$\Psi_{\Delta}(\mathbf{k})=\frac{4\pi}{2s+1}\int
d^{3}x\,e^{i\mathbf{k}\cdot\mathbf{x}}\sum\limits_{m=-s}^{s}(-1)^{m}\left[\frac{2}{k^{\Delta}}\mathcal{F}_{s,m}^{\,\Delta}({\bm{x}})\mathcal{F}_{s,-m}^{\,0}({\bm{x}})-\mathcal{F}_{s,m}^{\,0}({\bm{x}})\mathcal{F}_{s,-m}^{\,0}({\bm{x}})\right],$
(39)

where

$\mathcal{F}_{s,m}^{\,\gamma}({\bm{x}})\equiv\int_{\mathbf{p}}p^{\gamma}\,\phi(\mathbf{p})Y_{sm}(\hat{\mathbf{p}})\,e^{-i\mathbf{p}\cdot\mathbf{x}}.$
(40)

The squeezed bispectrum associated with Eq. (38) is

$\displaystyle\lim_{q\ll k}B_{\Phi}(\mathbf{q},\mathbf{k})$
$\displaystyle=4f_{\rm
NL}^{\Delta}\left[\left(\frac{q}{k}\right)^{\Delta}{\mathcal{L}_{s}(\hat{\mathbf{q}}\cdot\hat{\mathbf{k}})}P_{\Phi}(q)P_{\Phi}(k)+\left(\frac{k}{q}\right)^{\Delta}P_{\Phi}(k)P_{\Phi}(k)\right],$
(41) $\displaystyle\approx 4f_{\rm
NL}^{\Delta}\left[\left(\frac{q}{k}\right)^{\Delta}{\mathcal{L}_{s}(\hat{\mathbf{q}}\cdot\hat{\mathbf{k}})}+\left(\frac{q}{k}\right)^{3-\Delta}\right]P_{\Phi}(q)P_{\Phi}(k),$
(42)

where the second line assumes $P_{\Phi}(k)$ is scale-invariant. Thus, the
kernel generates the target squeezed bispectrum if $\Delta\ll 3/2.$ On the
other hand, particles with $s\neq 0$ must satisfy the Higuchi bound,
$m^{2}/H^{2}\geq s(s-1)$ [140], which implies $\Delta\geq 1$. Therefore, the
second term in the squeezed bispectrum non-negligible for the majority of the
physical parameter space for spin-$s$ particles. One can suppress this
unwanted contribution by using a high-pass filtered kernel
$K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})\rightarrow
K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})\Theta(k_{12}-k_{L})$, where $\Theta$
is the step function and $k_{L}$ is above all $q$-values of interest (e.g.,
$k_{L}\approx 0.1~{}h/{\rm Mpc}$ in this work). Alternatively, if
$\Delta=3/2$, one can get rid of the unwanted contribution by subtracting off
the kernel for intermediate-mass scalars with $\Delta=3/2$, _i.e._ ,
subtracting Eq. (32) with an additional factor of $1/2$.

###  A.2 Oscillatory bispectra from massive scalars

For particles with mass $m>3H/2$, the squeezed bispectrum oscillates
logarithmically as

$\lim_{q\ll k}B_{\Phi}(\mathbf{q},\mathbf{k})=4f_{\rm
NL}^{\Delta}\left(\frac{q}{k}\right)^{3/2}\cos\left(\nu\log\left(\frac{q}{k}\right)+\delta\right)\mathcal{L}_{s}(\hat{\mathbf{q}}\cdot\hat{\mathbf{k}})P_{\Phi}(q)P_{\Phi}(k),$
(43)

where $\nu=\sqrt{m^{2}/H^{2}-9/4}$ for scalars and
$\nu=\sqrt{m^{2}/H^{2}-(s-1/2)^{2}}$ for $s\neq 0.$ The phase ($\delta$) is a
function of the mass and spin, and depends on the precise inflationary
coupling (see [44, 81]). For simplicity, we focus on scalar fields, but note
that these results can be generalized to spinning fields using the procedure
developed for intermediate-mass spinning particles. For scalar fields, we
consider the following kernel:

$\displaystyle\begin{split}K_{\Delta}(\mathbf{k}_{1},\mathbf{k}_{2})&=\left\\{\left(\frac{k_{1}}{k_{12}}\right)^{3/2}\cos\left(\nu\log\left(\frac{k_{1}}{k_{12}}\right)+\delta\right)+\left(\frac{k_{2}}{k_{12}}\right)^{3/2}\cos\left(\nu\log\left(\frac{k_{2}}{k_{12}}\right)+\delta\right)-\cos(\delta)\right\\},\\\
&=\frac{1}{2}\left\\{e^{i\delta}\left[\left(\frac{k_{1}}{k_{12}}\right)^{3/2+i\nu}+\left(\frac{k_{2}}{k_{12}}\right)^{3/2+i\nu}-1\right]+e^{-i\delta}\left[\left(\frac{k_{1}}{k_{12}}\right)^{3/2-i\nu}+\left(\frac{k_{2}}{k_{12}}\right)^{3/2-i\nu}-1\right]\right\\},\end{split}$
(44)

where the second line is written in a separable form that can be easily
evaluated using FFTs. For scale-invariant $P_{\Phi}(k)$, the squeezed
bispectrum associated with Eq. (44) is, up to $\mathcal{O}\left(f_{\rm
NL}^{\Delta}\right)$,

$\displaystyle\begin{split}\lim_{q\ll
k}B_{\Phi}(\mathbf{q},\mathbf{k})=&4f_{\rm
NL}^{\Delta}\left(\frac{q}{k}\right)^{3/2}\left[\cos\left(\nu\log\left(\frac{q}{k}\right)+\delta\right)+\cos\left(\nu\log\left(\frac{q}{k}\right)-\delta\right)\right]P_{\Phi}(q)P_{\Phi}(k).\end{split}$
(45)

Similar to the intermediate-mass case, we can apply a high-pass filter to the
kernel to suppress the second term in Eq. (45). Alternatively, one could
always fix $\delta=0$ and multiply the kernel by $1/2$. Although this
$\delta=0$ template lacks physical motivation, simulations with $\delta=0$
would likely suffice to study the impact of oscillatory bispectra on large-
scale structure in the non-linear regime.

##  Appendix B Bispectrum and trispectrum estimators

In this appendix, we derive estimators for the 3D real-space bispectrum and
trispectrum from a simulation volume. These estimators can be efficiently
computed by Fast Fourier Transformations (FFTs). We implement these estimators
in the publicly available code
PNGolin.313131https://github.com/samgolds/PNGolin

###  B.1 Bispectrum estimator

The bispectrum is parameterized by three external momenta, $k_{1}$, $k_{2}$,
and $k_{3}$. To estimate the bispectrum from a density field,
$\delta_{\mathbf{k}}\equiv\delta(\mathbf{k})$, we introduce a window
$W_{k_{i}}(p)$ which is equal to $1$ if $p$ is in $k_{i}$ and zero otherwise.
A bispectrum estimator is then

$\displaystyle\hat{B}(k_{1},k_{2},k_{3})\propto\int\limits_{\mathbf{p}_{1},\dots,\mathbf{p}_{3}}\bigg{[}$
$\displaystyle
W_{k_{1}}(p_{1})W_{k_{2}}(p_{2})W_{k_{3}}(p)\,\delta_{\mathbf{p}_{1}}\delta_{\mathbf{p}_{2}}\delta_{\mathbf{p}_{3}}(2\pi)^{3}\delta_{D}^{(3)}(\mathbf{p}_{123})\bigg{]},$
(46)

where we have ignored the overall normalization, which can be computed by re-
evaluating the integral with all density fields replaced by the identity. To
proceed, we write the Dirac delta as an exponential integral,

$\delta_{D}(\mathbf{k})=\frac{1}{(2\pi)^{3}}\int
d^{3}xe^{-i\mathbf{k}\cdot{\bm{x}}}.$ (47)

The bispectrum estimator is then

$\hat{B}(k_{1},k_{2},k_{3})\propto\int
d^{3}x\left[\left(\int_{\mathbf{p}_{1}}e^{-i\mathbf{p}_{1}\cdot\mathbf{x}}\,W_{k_{1}}(p_{1})\,\delta_{\mathbf{p}_{1}}\right)\left(\int_{\mathbf{p}_{2}}e^{-i\mathbf{p}_{2}\cdot\mathbf{x}}\,W_{k_{2}}(p_{2})\,\delta_{\mathbf{p}_{2}}\right)\left(\int_{\mathbf{p}_{3}}e^{-i\mathbf{p}_{3}\cdot\mathbf{x}}\,W_{k_{3}}(p_{3})\,\delta_{\mathbf{p}_{3}}\right)\right],$
(48)

which can be efficiently evaluated using FFTs. In this work, we analyze the
squeezed bispectrum with $k_{2}=k_{3}$, making the last two inverse Fourier
transforms in Eq. (48) equal.

###  B.2 Trispectrum estimator

Our trispectrum estimator is based on the parity-even trispectrum estimator
described in Appendix A of Ref. [141]. The trispectrum can be parameterized by
four external momenta, $k_{1}$, $k_{2}$, $k_{3}$, and $k_{4}$, as well as the
two internal momenta, $K\equiv k_{12}$ and $K^{\prime}=k_{14}$ (see Fig. 1).
The full trispectrum estimator is not separable, thus we integrate over
$k_{14}$ and estimate the trispectrum with external momenta in bins
$k_{1},\dots k_{4}$ and internal momentum $k_{12}\in K$.323232See Ref. [142]
for an estimator that is also integrated over $k_{12}.$ An estimator for the
total four-point function, _i.e._ , including disconnected contributions, is

$\displaystyle\hat{T}_{\rm
tot.}(k_{1},k_{2},k_{3},k_{4},K)\propto\int\limits_{{\bm{P}}}W_{K}(P)\int\limits_{\mathbf{p}_{1},\dots,\mathbf{p}_{4}}\bigg{[}$
$\displaystyle
W_{k_{1}}(p_{1})W_{k_{2}}(p_{2})W_{k_{3}}(p_{3})W_{k_{4}}(p_{4})\,\delta_{\mathbf{p}_{1}}\delta_{\mathbf{p}_{2}}\delta_{\mathbf{p}_{3}}\delta_{\mathbf{p}_{4}}$
(49)
$\displaystyle(2\pi)^{3}\delta_{D}^{(3)}(\mathbf{p}_{12}-{\bm{P}})(2\pi)^{3}\delta_{D}^{(3)}(\mathbf{p}_{34}+{\bm{P}})\bigg{]},$

where we have again ignored the overall normalization, which can be computed
by re-evaluating the integral with all density fields replaced by the
identity. Eq. (49) can be efficiently evaluated using FFTs by replacing the
Dirac delta functions with exponential integrals, yielding

$\hat{T}_{\rm
tot.}(k_{1},k_{2},k_{3},k_{4},K)\propto\int_{{\bm{P}}}W_{K}(P)\left(\int
d^{3}x\,e^{-i{\bm{P}}\cdot{\bm{x}}}\delta_{W_{k_{1}}}({\bm{x}})\delta_{W_{k_{2}}}({\bm{x}})\right)\left(\int
d^{3}y\,e^{i{\bm{P}}\cdot{\bm{y}}}\delta_{W_{k_{3}}}({\bm{y}})\delta_{W_{k_{4}}}({\bm{y}})\right),$
(50)

where
$\delta_{W_{k_{i}}}({\bm{x}})\equiv\int_{\mathbf{p}}W_{k_{i}}(p)\delta_{\mathbf{p}}e^{i\mathbf{p}\cdot{\bm{x}}}$.

Figure 11: Validation of trispectrum estimator on 200 simulated Gaussian
random fields with a scale-invariant power spectrum. The left panel shows the
connected (red) and disconnected (blue) contributions to the four-point
function estimator (Eq. (49)). The disconnected contributions are large and
must be subtracted to estimate the trispectrum. The trispectrum (red) is
consistent with zero. The right panel shows the ratio of the trispectrum
($T_{\rm conn.}$) to the total four-point function ($T_{\rm tot.})$ for
several estimates of the disconnected term. The optimal estimator (black) has
lower bias and lower variance than the naive estimators (blue and red). Note
that since we generate the Gaussian random fields with a _known_ power
spectrum, the bias is due to finite grid size effects. We discuss this in more
detail in the main text.

Eq. (49) includes disconnected contributions that must be subtracted to
estimate the trispectrum. In the limit of mild non-Gaussianity, the optimal
estimator for the trispectrum is [143, 144],

$\hat{T}_{\rm
opt}=\delta^{4}-6\,\delta^{2}\langle\delta^{2}\rangle+3\langle\delta^{2}\rangle^{2}.$
(51)

Therefore, we can derive an optimal estimator for the trispectrum by replacing
$\delta^{4}\rightarrow\delta^{4}-6\,\delta^{2}\langle\delta^{2}\rangle+3\langle\delta^{2}\rangle^{2}$
in Eq. (49). The $3\langle\delta^{2}\rangle^{2}$ term corresponds to a
_realization-independent_ disconnected contribution, which is given by

$\hat{T}_{\rm
disc.}^{\,3\langle\delta^{2}\rangle^{2}}\propto\int\limits_{{\bm{P}}}W_{K}(P)\left(\prod\limits_{i=1}^{4}\int\limits_{\mathbf{p}_{i}}W_{k_{i}}(p_{i})\right)\left[P(p_{1})P(p_{3})\,\delta_{D}^{(3)}(\mathbf{p}_{12})\,\delta_{D}^{(3)}(\mathbf{p}_{34})+2~{}{\rm
perm.}\right]\delta_{D}^{(3)}(\mathbf{p}_{12}-{\bm{P}})\delta_{D}^{(3)}(\mathbf{p}_{34}+{\bm{P}}),$
(52)

where $P(p_{i})$ is calculated from a fiducial power spectrum and we have
omitted the momentum dependence on the left-hand side for clarity. Assuming
${\bm{K}}\neq 0$, the first term in Eq. (52) is zero. The remaining two
permutations are

$\displaystyle\hat{T}_{\rm disc.}^{\,3\langle\delta^{2}\rangle^{2}}$
$\displaystyle\propto\int\limits_{{\bm{P}}}W_{K}(P)\left(\prod\limits_{i=1}^{4}\int\limits_{\mathbf{p}_{i}}W_{k_{i}}(p_{i})\right)\left[P(p_{1})P(p_{2})\,\delta_{D}^{(3)}(\mathbf{p}_{13})\,\delta_{D}^{(3)}(\mathbf{p}_{24})\right]\delta_{D}^{(3)}(\mathbf{p}_{12}-{\bm{P}})\delta_{D}^{(3)}(\mathbf{p}_{34}+{\bm{P}})+(3\leftrightarrow
4),$
$\displaystyle\propto\int_{{\bm{P}}}W_{K}(P)\left(\int_{\mathbf{p}_{1}}W_{k_{1}}(p_{1})W_{k_{3}}(p_{1})P(p_{1})\right)\left(\int_{\mathbf{p}_{2}}W_{k_{2}}(p_{2})W_{k_{4}}(p_{2})P(p_{2})\right)\left[(2\pi)^{3}\delta_{D}^{(3)}(\mathbf{p}_{12}-{\bm{P}})\right]^{2}+(3\leftrightarrow
4).$

To proceed, we approximate the squared Delta function using
$((2\pi)^{3}\delta^{(3)}_{D}(\mathbf{k}))^{2}=V\times(2\pi)^{3}\delta^{(3)}_{D}(\mathbf{k})$
(see, _e.g._ , [142]) and rewrite the final Delta function with an exponential
integral, yielding

$\hat{T}_{\rm
disc.}^{\,3\langle\delta^{2}\rangle^{2}}=\frac{V}{N_{k_{1},k_{2},k_{3},k_{4},K}}\int_{{\bm{P}}}W_{K}(P)\int
d^{3}x\,e^{-i{\bm{P}}\cdot{\bm{x}}}\left[F^{P}_{13}({\bm{x}})F^{P}_{24}({\bm{x}})+F^{P}_{14}({\bm{x}})F^{P}_{23}({\bm{x}})\right],$
(53)

where

$F^{P}_{ij}({\bm{x}})\equiv\int_{\mathbf{p}}W_{k_{i}}(p)W_{k_{j}}(p)P(p)e^{i\mathbf{p}\cdot\bm{x}},$
(54)

and $N_{k_{1},k_{2},k_{3},k_{4},K}$ is the number of tetrahedra in a bin.

Following the same reasoning, the $6\delta^{2}\langle\delta^{2}\rangle$ term
corresponds to a _realization-dependent_ disconnected contribution,

$\displaystyle\hat{T}_{\rm
disc.}^{6\,\delta^{2}\langle\delta^{2}\rangle}=\frac{V}{N_{k_{1},k_{2},k_{3},k_{4},K}}\int_{{\bm{Q}}}W_{E}(Q)\int
d^{3}x\,e^{-i{\bm{Q}}\cdot{\bm{x}}}[$ $\displaystyle
F^{P}_{13}({\bm{x}})F^{\delta}_{24}({\bm{x}})+F^{\delta}_{13}({\bm{x}})F^{P}_{24}({\bm{x}})$
$\displaystyle+F^{P}_{14}({\bm{x}})F^{\delta}_{23}({\bm{x}})+F^{\delta}_{14}({\bm{x}})F^{P}_{23}({\bm{x}})],$
(55)

where

$F^{\delta}_{ij}({\bm{x}})\equiv\int_{\mathbf{k}}W_{k_{i}}(p)W_{k_{j}}(p)|\delta_{\mathbf{p}}|^{2}e^{i\mathbf{k}\cdot\bm{x}}.$
(56)

To summarize, an optimal estimator for the binned trispectrum is $\hat{T}_{\rm
opt}=\hat{T}_{\rm tot.}-\hat{T}_{\rm disc.}^{6\,\delta^{2}}-\hat{T}_{\rm
disc.}^{\,3\langle\delta^{2}\rangle^{2}}$, where each term can be written in
terms of FFTs using Eqs. (49), (53), and (B.2).

In Fig. 11, we validate our trispectrum estimator using simulations of 200
Gaussian random fields with a scale-invariant power spectrum. We measure the
trispectrum for all bins with $4k_{f}<k_{1},k_{2},k_{3},k_{4},k_{12}<10k_{f}$
and width $\Delta k_{i}=1k_{f},$ where $k_{f}$ is the fundamental mode. The
left panel shows the connected (red) and disconnected (blue) contributions to
the four-point function. The trispectrum, _i.e._ , the connected contribution,
is consistent with zero as expected for a Gaussian random field. On the other
hand, the disconnected terms are quite large and must be carefully subtracted
to obtain an unbiased estimate of the trispectrum.

The right panel of Fig. 11 shows the ratio of the estimated trispectrum to the
total four-point function using three different estimates of the disconnected
contribution. The optimal estimator (black) has lower bias and variance than
the naive estimators (blue and red). Since we generate the Gaussian random
fields with a known fiducial power spectrum, the observed bias is solely due
to numerical effects. Specifically, we simulate the Gaussian random fields
with a relatively low resolution ($N_{\rm grid}=128$), thus there is a 2-3%
discrepancy between the measured power spectrum and fiducial power spectrum
for the bins analyzed here. This small bias can lead to more significant