Searching for Cosmological Collider in the Planck CMB Data

Wuhyun Sohn1, Dong-Gang Wang2, James R. Fergusson2, and E. P. S. Shellard2

1 Korea Astronomy and Space Science Institute, Daejeon 34055, South Korea

2Centre for Theoretical Cosmology, Department of Applied Mathematics and
Theoretical Physics
University of Cambridge, Wilberforce Road, Cambridge, CB3 0WA, UK

Abstract

In this paper, we present the first comprehensive CMB data analysis of
cosmological collider physics. New heavy particles during inflation can leave
imprints in the primordial correlators which are observable in today’s
cosmological surveys. This remarkable detection channel provides an
unsurpassed opportunity to probe new physics at extremely high energies. Here
we initiate the search for these relic signals in the cosmic microwave
background (CMB) data from the Planck legacy release. On the theory side,
guided by recent progress from the cosmological bootstrap, we first propose a
family of analytic bispectrum templates that incorporate the distinctive
signatures of cosmological collider physics. Our consideration includes the
oscillatory signals in the squeezed limit, the angular dependence from
spinning fields, and several new shapes from nontrivial sound speed effects.
On the observational side, we apply the recently developed pipeline, CMB
Bispectrum Estimator (CMB-BEST), to efficiently analyze the three-point
statistics and search directly for these new templates in the Planck 2018
temperature and polarization data. We report stringent CMB constraints on
these new templates. Furthermore, we perform parameter scans to search for the
best-fit values with maximum significance. For a benchmark example of collider
templates, we find $f_{\rm NL}=-91\pm 40$ at the $68\%$ confidence level.
After accounting for the look-elsewhere effect, the biggest adjusted
significance we get is $1.8\sigma$. In general, we find no significant
evidence of cosmological collider signals in the Planck data. However, this
innovative analysis demonstrates the potential for discovering new heavy
particles during inflation in forthcoming cosmological surveys.

###### Contents

  1. 1 Introduction
  2. 2 Shapes of the Cosmological Collider
    1. 2.1 PNG: A Bootstrap Recap
    2. 2.2 Shapes from massive scalar exchange
    3. 2.3 Shapes from spinning exchange
    4. 2.4 Shapes with small sound speeds
  3. 3 CMB Bispectrum Estimator
    1. 3.1 CMB bispectrum statistics
    2. 3.2 CMB bispectrum estimator
  4. 4 The Search for Collider Signals in the Planck Data
    1. 4.1 Oscillations from massive scalars
    2. 4.2 Angular dependence from spinning particles
    3. 4.3 Shapes with small sound speeds
  5. 5 Summary and Outlook
  6. A Approximated Scalar Seed

##  1 Introduction

One remarkable aspect of inflationary cosmology is that the earliest stage of
our Universe also serves as a natural laboratory for high energy physics,
where the effects of the microscopic quantum world during inflation are
imprinted on the macroscopic classical spacetime, which can be probed later
through astronomical observations [1, 2]. One of the most promising channels
for detecting signatures of the early universe lies in the higher-order
cosmological correlators of primordial fluctuations that are formed during
inflation. As these spatial correlations capture deviations from Gaussian
statistics, they are usually known as primordial non-Gaussianity (PNG), and
for this reason they have been an important observational target in on-going
and upcoming experiments, such as those mapping the cosmic microwave
background (CMB) and for large scale structure (LSS) surveys. Given that
inflation may reach the highest energy densities that are observationally
accessible in nature ($\lesssim 10^{13}$GeV), the search for PNG provides a
unique and exciting opportunity to probe fundamental physics at energy scales
far beyond the reach of terrestrial experiments.

This confrontation between observational cosmology and high energy theory
finds a particularly lucid exposition in the “Cosmological Collider Physics”
program [3, 4, 5, 6]. The analogy is drawn with particle accelerators where
interactions with intermediate particles is imprinted in the scattering
amplitudes determined in flat spacetime. For new heavy particles during
inflation, likewise, they can mediate specific correlations of primordial
fluctuations, with their mass and spin leading to distinctive signatures in
the shape of PNG. Therefore, this provides the opportunity to do particle
spectroscopy by measuring the statistics of primordial fluctuations.
Considering the rich phenomenology available in cosmological observations and
their potentially deep implications for particle physics, theoreticians have
extensively studied various possibilities of the cosmological collider physics
in the past decade [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44].

The recent development of the cosmological bootstrap has further enhanced
theoretical efforts along this direction [45, 46, 47, 48, 49, 50, 51, 52, 53,
54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72,
73]. This new approach aims to directly determine the forms of cosmological
correlators from a set of fundamental principles, such as symmetry, unitarity
and locality. For this reason, the boostrap methodology makes it possible to
carve out the theory space in a model-independent way, and achieve a
systematic classification of the possible categories for all the predictions.
In addition, the approach also provides powerful computational tools. For the
first time, we are able to compute the shape function of cosmological
colliders for any kinematics, and identify the full analytic structure of the
corresponding correlators. These latest advances have inspired a more complete
survey revealing non-Gaussianity predictions with new features and large
signals that can be tested even using currently available data.

At the observational frontier the searches for cosmological collider signals
in real data are far from exhausted and discovery potential remains. The
leading target of PNG has been the bispectrum of primordial curvature
perturbations, with most of the data analysis to date focused on simpler
templates from multi-field and single field inflation models, such as the
local and equilateral shapes. The latest and tightest constraints on these
non-Gaussianities come from the Planck CMB experiment [74] (see also [75,
76]), which included general methods to cover a much wider array of
theoretically-motivated templates [77, 78, 79], though only few were directly
inspired by cosmological colliders (see, e.g. [80]). LSS surveys have also
placed some constraints on PNG using the BOSS data [81, 82, 83], albeit being
very much weaker than those currently available from the CMB. For the
cosmological collider signal, there are already several forecast studies
relevant to future galaxy surveys [84, 85, 86], but no comprehensive CMB data
analysis has been performed yet. One major difficulty is that the shape
functions with collider signals are very complicated in general due to
oscillations and various types of special functions. This poses a number of
technical challenges for the data analysis in order to extract the signal of
physical interest.

The CMB bispectrum data analysis is, in general, computationally challenging
because a naive implementation scales as $O(\ell_{\mathrm{max}}^{5})$, with
$\ell_{\mathrm{max}}\approx 2500$ for Planck. Three main approaches have been
utilised in the literature to mitigate this problem: the KSW estimator [87,
88, 89, 90, 91], exploiting the separability of simple approximate templates
reduces high-dimensional integrals into products of lower dimensional
integrals; the Modal estimator [77, 78, 92], which expands arbitrary
bispectrum shapes to high precision using separable basis functions; and the
related binned bispectrum estimator [79, 93], which uses multipole bins to
compress the data. However, the KSW and binned estimators are not well suited
to an extensive study of cosmological colliders, because the resulting
bispectrum shapes tend to be non-separable and highly oscillatory.

In this work, we use the high-resolution CMB bispectrum estimator with a
publicly available Python code, named CMB-BEST (CMB Bispectrum ESTimator)
[94]. The formalism combines the advantages of two conventional methods (KSW
and Modal) using a flexible set of modal basis functions in the primordial
space and it has been extensively validated against the Modal bispectrum
pipeline used for the Planck analysis [78, 74]. All computationally expensive
steps are precomputed at a given modal resolution into a data file provided
with the code, so that the users can obtain the Planck 2018 constraints on
arbitrary shape functions of interest rapidly (within a minute on a laptop).
Because of CMB-BEST’s flexible mode expansions and its speed and resolution,
it is an ideal tool for studying cosmological collider signals.

Shape | Template |  ${f_{\mathrm{NL}}}$ (68% CL) | Raw S/N | Adjusted S/N | Section
---|---|---|---|---|---
Quasi-single field [3] | (2.6) | $10\pm 26$ | 0.37 | 0.12 | 4.1
Scalar exchange I | (2.2) | $11\pm 13$ | 0.86 | 0.67 | 4.1
Scalar exchange II | (2.20) | $-91\pm 40$ | 2.3 | 1.8 | 4.1
Heavy-spin exchange | (2.24) | $-59\pm 32$ | 1.9 | 1.2 | 4.2
Massive spin-2 exchange | (2.27) | $-2.1\pm 1.1$ | 1.9 | 0.90 | 4.2
Equilateral collider [59] | (2.32) | $-178\pm 72$ | 2.5 | 0.90 | 4.3
Low-speed collider [41] | (2.33) | $-9\pm 10$ | 0.89 | 0.29 | 4.3
Multi-speed PNG [64] | (2.34) | $-3.1\pm 2.3$ | 1.3 | 0.61 | 4.3
Table 1: Summary of the CMB bispectrum constraints presented in this work.
Most templates have free parameters, which were set to their best-fit values
for the ‘constraints’ column in this Table; refer to the corresponding section
for the full constraints. The adjusted signal-to-noise incorporates the look-
elsewhere effect, also detailed in the section.

On the basis of the theoretical survey presented here, we have deployed the
CMB-BEST pipeline to perform the first comprehensive CMB search for
cosmological collider signals, yielding the most precise measurements to date
of the wide array of relevant bispectrum shapes. The originality of this work
is twofold:

  * •

We derive a set of cosmological collider bispectrum templates that are
relatively simple and sufficiently accurate across the whole observational
domain. With the help of the bootstrap approach, we exploit the systematic
classification and complete analytic understanding of the bispectrum shapes to
capture the relevant possibilities and to simplify the analytical expressions.
In addition, we also endeavour to collate all the various types of shape
ansatzes given in the cosmological collider literature. Our consideration of
the shape functions incorporates: 1) the oscillatory signals from massive
fields; 2) the angular-dependence of profiles due to spinning particles; and
3) three additional new shapes caused by the effect of varying the sound
speed. One particular advantage of these simplified bispectrum templates is
that it allows us to perform parameter scans in the present data analysis.
Their analytical expressions are shown in grey boxes in the text and they can
also be directly applied to future observational surveys.

  * •

We perform a systematic search for these collider templates in the latest
Planck temperature and polarization data. Our main results are summarized in
Table 1, which lists the measurements of the collider bispectrum templates
presented in the different referenced sections of this work. These summary
statistics are given as “constraints”, but actually represent the highest
signal-to-noise measurements obtained, that is, the best-fit parameter choices
after scans over the available parameter space for which both the theoretical
predictions and observational methods are robust.

To offer a brief overall summary at the outset, we find no significant
evidence for cosmological collider signals given the resolution and
sensitivity available with the present Planck CMB data, and within the
improved theoretical approximations made for this investigation. As the first
full analysis of cosmological collider signals, we expect the methodology we
have presented here will be directly applicable to forthcoming observations,
notably for Simons Observatory and CMB-S4, but also future large-scale
structure surveys, which forecasts suggest may become competitive with the
CMB.

The structure of the paper is organized as follows. In Section 2, we first
review the classification of the primordial bispectrum shapes from the
bootstrap perspective, and then based on the exact analytical solution, we
propose a set of simplified templates for various types of cosmological
collider signals, for the convenience of the data analysis. More technical
details about how to construct the approximated shape functions are left in
Section A. In Section 3, we introduce basic setup and methodology of the
python pipeline CMB-BEST. In Section 4, we present the CMB constraints on
various shape templates of cosmological colliders using the latest Planck
data. We conclude in Section 5.

##  2 Shapes of the Cosmological Collider

In this section, we briefly review the non-Gaussian shape functions of the
primordial scalar bispectra in cosmological collider scenarios. We exploit the
full analytic understanding of the three-point functions from the bootstrap
analysis, and then propose the approximated shape templates for various types
of new physics signatures.

###  2.1 PNG: A Bootstrap Recap

The major target of our analysis is the primordial bispectrum of the curvature
perturbation $\zeta=(H/\dot{\Phi})\phi$

$\displaystyle\langle\zeta_{{\bf k}_{1}}\zeta_{{\bf k}_{2}}\zeta_{{\bf
k}_{3}}\rangle=(2\pi)^{3}\delta({\bf k}_{1}+{\bf k}_{2}+{\bf
k}_{3})\frac{18}{5}f_{\rm
NL}\frac{S(k_{1},k_{2},k_{3})}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}P_{\zeta}^{2}\leavevmode\nobreak\
,$ (2.1)

where $P_{\zeta}$ is the power spectrum of $\zeta$ and $f_{\rm NL}$ represents
the size of the non-Gaussian signal. The information of new physics during
inflation is mainly captured by the dimensionless function
$S(k_{1},k_{2},k_{3})$.

In our current understanding of PNG, if we assume the inflation theory is
(nearly) scale-invariant and weakly coupled, there are three broad classes of
scalar bispectra, which correspond to the Feynman diagrams shown in Figure 1.
One remarkable advantage of the bootstrap approach is that, without
referencing to specific models, we are allowed to systematically classify all
the possible bispectrum shapes based on the minimal assumptions. The PNG
predictions from the three leading scenarios of inflation are summarized as
follows.

  * •

Single field inflation: In this scenario, the scalar bispectrum is generated
by the self-interaction of the inflaton. As the inflaton field is normally
expected to be protected by a shift symmetry, we have derivative interactions.
The bootstrap analysis of single field scenario is presented by the boostless
bootstrap [54, 55, 56], which starts from an ansatz of rational polynomials
and then keeps imposing constraints from scale invariance, Bose symmetry,
flat-space limit and locality. In the end the shape functions from derivative
interactions can be written into the following form

$S^{\rm eq}(k_{1},k_{2},k_{3})=\frac{{\rm
Poly}_{p+3}(k_{T},e_{2},e_{3})}{k_{1}k_{2}k_{3}k_{T}^{p}}\leavevmode\nobreak\
,$ (2.2)

where $k_{T}=k_{1}+k_{2}+k_{3}$, $e_{2}=k_{1}k_{2}+k_{2}k_{3}+k_{1}k_{3}$ and
$e_{3}=k_{1}k_{2}k_{3}$ are symmetric polynomials. The index $p$ is an integer
related to the number of derivatives in the cubic vertex, while the form of
the $(p+3)$-order polynomial in the numerator is fixed by the flat-space
amplitude and the manifestly local test. From the perspective of effective
field theories, the lowest order interactions are given by $\dot{\phi}^{3}$
and $\dot{\phi}(\partial_{i}\phi)^{2}$, which may be seen as consequences of
integrating out heavy physics with masses much larger than Hubble. As in this
case the contact field interactions are local, and enhanced around (sound)
horizon-exit, the resulting bispectra are usually peaked at the equilateral
configuration with $k_{1}=k_{2}=k_{3}$, which are known as the equilateral-
type non-Gaussianity.

Figure 1: Feynman diagrams for three leading predictions of primordial
bispectra from inflation. In these diagrams, time flows from bottom to top,
and the horizontal grey lines correspond to the end of inflation. All the
external legs (black lines) are the inflaton fluctuations $\phi$ which are
associated to the curvature perturbation $\zeta$. Each diagram describes one
scenario of cosmic inflation with distinctive non-Gaussian signatures: a) the
contact diagram for single field inflation; b) the massive exchange diagram
for cosmological colliders; c) the massless exchange diagram for multi-field
inflation. This work mainly focuses on testing non-Gaussian signals from the
second diagram.

  * •

Multi-field inflation: Non-Gaussianities from additional light scalars during
inflation have been extensively studied in literature. The traditional
approach is to take the separate universe approximation, and the $\delta N$
formalism there intuitively captures the conversion process, which leads to
the famous local shape $S^{\rm
loc}=(k_{1}^{3}+k_{2}^{3}+k_{3}^{3})/(3k_{1}k_{2}k_{3})$. This type of signal
can be understood from the massless exchange diagram shown in the last diagram
in Figure 1, where the intermediate light isocurvature field can travel a long
distance before converting into the curvature perturbation. This leads to
shape functions peaked in the squeezed limit $k_{3}\ll k_{1}\simeq k_{2}$,
which are dramatically different from the single field ones. A recent
bootstrap analysis shows that the local-type PNG are consequences of IR
divergences in de Sitter space, which in field-theoretic computations lead to
mild logarithmic deviation of the shape functions from the standard local
ansatz [64].

  * •

Cosmological collider: roughly speaking, this scenario can be seen as an
interpolation between the previous two. As we have argued, when the extra
fields are heavy $m\gg H$, they quickly decay into the inflaton which leads to
contact interactions in the effective single field inflation. The opposite
situation happens for additional light scalars with $m\ll H$ where the long-
lived internal state leads to multi-field non-Gaussianities. Their differences
can be demonstrated more clearly by taking the squeezed limit of the
primordial bispectra:

$\displaystyle\lim_{k_{l}\ll k_{s}}S^{\rm eq}$ $\displaystyle\propto$
$\displaystyle\frac{k_{l}}{k_{s}}$ (2.3) $\displaystyle\lim_{k_{l}\ll
k_{s}}S^{\rm loc}$ $\displaystyle\propto$ $\displaystyle\frac{k_{s}}{k_{l}}$
(2.4)

In addition to the above single field and multi-field non-Gaussianities, the
cosmological collider captures the intermediate regime of massive states
$m\sim H$, where the squeezed bispectrum interpolates the two scalings above.
Furthermore, rich phenomenologies are expected, from which we may be able to
extract information about the mass and spin of the unknown heavy particles
during inflation. In the rest of this section, we shall elaborate on this by
looking into various possibilities of the shape functions of cosmological
colliders.

On the observation frontier, there have been many efforts for testing the
local template and the equilateral-type shape functions (normally for
equilateral and orthogonal templates), using both the CMB data and the LSS
surveys. But few data analysis has been performed for cosmological colliders
yet, even though the phenomenology there is richer. One major challenge arises
from the theory side, as the analytical forms of predicted shape functions (if
exist) are very complicated. In the rest of this section, we shall tackle this
problem by proposing well-motivated simple templates.

###  2.2 Shapes from massive scalar exchange

One characteristic signal of cosmological colliders is the oscillations in the
squeezed bispectrum. This corresponds to the generalized scaling behaviour in
the soft limit, which can already been identified in massive scalar exchanges.
Through the Feynman diagram in Figure 1b), the super-horizon decay of the
massive particle during inflation leaves the imprints in cosmological
correlators. In general, there are two classes depending on the mass of the
intermediate particle.

  * •

When the massive scalar is in the complimentary series ($0<m<3H/2$), the field
demonstrates a power-law decay with conformal time after horizon-exit. The
signatures of massive scalars were first extensively studied in the quasi-
single field (QSF) inflation111In general, the quasi-single field models
correspond to an inflation scenario with $m\sim H$ massive fields, which also
includes the $m>3H/2$ case. We use the acronym “QSF” for the $m<3H/2$ template
as introduced in [3]., whose squeezed bispectrum takes the following form [3,
4, 8]

$\lim_{k_{l}\ll k_{s}}S^{\rm
QSF}\propto\left(\frac{k_{l}}{k_{s}}\right)^{\frac{1}{2}-\nu}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ {\rm with}\leavevmode\nobreak\
\nu=\sqrt{\frac{9}{4}-\frac{m^{2}}{H^{2}}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ (2.5)

with the index $\nu\in(0,3/2)$ given by the mass. For this type of models, an
approximated shape template has been proposed in [3]

$S^{\rm
QSF}(k_{1},k_{2},k_{3})=3\sqrt{3\kappa}\frac{N_{\nu}(8\kappa)}{N_{\nu}({8}/{27})}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ {\rm with}\leavevmode\nobreak\
\kappa=\frac{k_{1}k_{2}k_{3}}{k_{T}^{3}}$ (2.6)

where $N_{\nu}$ is the Neumann function. In the squeezed limit this template
can reproduce the correct scaling behaviour in (2.5). Meanwhile for the non-
squeezed kinematics configuration, by varying the mass parameter, the full
shape provides an interpolation between two types of templates [3]: when
$\nu\rightarrow 3/2$, it goes back to the local shape from multi-field
inflation; when $\nu\rightarrow 0$ we find it approaches the constant shape.
The observational test of this template has been analyzed by using the WMAP7
data in [80], and by the Planck13 data in [75]. We shall update the constraint
using the Planck legacy release in Section 4.

  * •

The collider signal arises when the massive scalar is in the principle series
($m>3H/2$). There on super-horizon scales the scalar particle decays and
oscillates logarithmically in conformal time. Through the exchange diagram,
this produces a distinctive signature in the the scaling of the squeezed
bispectrum [5, 6, 12]

$\lim_{k_{l}\ll k_{s}}S^{\rm
col.}\propto\left(\frac{k_{l}}{k_{s}}\right)^{\frac{1}{2}+i\mu}+c.c.\propto\leavevmode\nobreak\
\left(\frac{k_{l}}{k_{s}}\right)^{{1}/{2}}\cos\left[\mu\ln\left(\frac{k_{l}}{k_{s}}\right)+\delta\right]\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
with}\leavevmode\nobreak\
\mu=\sqrt{\frac{m^{2}}{H^{2}}-\frac{9}{4}}\leavevmode\nobreak\ .$ (2.7)

From the oscillatory signal above, it becomes possible for us to do particle
spectroscopy during inflation, as the frequency of the oscillation measures
the mass of the intermediate state. Away from the squeezed limit, in general
for massive exchanges, it is very difficult to analytically derive the full
shape functions of the scalar bispectrum. While for the data analysis, it is
important to have approximated bispectrum templates, like (2.6) for the quasi-
single field scenario. In the following, we will tackle this problem with the
help of the cosmological bootstrap.

#### Bootstrapping scalar exchange templates

One particular advantage of the bootstrap approach is to provide a powerful
computational tool that allows us to gain the full analytical understanding of
bispectrum shapes for any kinematics [45, 59, 60, 64]. Here we follow the
approach developed in the boostless cosmological collider bootstrap [59] and
focus on deriving simpler analytical shape templates that can be easily
applied in data analysis.

The starting point is an auxiliary correlator called the three-point scalar
seed $\langle\varphi^{2}\phi\rangle$, where $\varphi$ is the conformally
coupled scalar ($m^{2}=2H^{2}$), and $\phi$ is the massless scalar. This
correlator is generated by the single exchange of a massive scalar $\sigma$,
with one cubic vertex $\varphi^{2}\sigma$ and one linear mixing
$\dot{\phi}\sigma$. Then in the $s$-channel where $k_{s}=k_{3}$ is associated
with the linear mixing leg, the three-point scalar seed can be written as a
function of the following momentum ratio 222In general, we can also introduce
two sound speed parameters for $\phi$ and $\sigma$ fields in the scalar seed.
Here for simplicity we focus on the situation where both $\phi$ and $\sigma$
are canonical scalars, and leave the effects of nontrivial sound speed in
Section 2.4.

$\langle\varphi^{2}\phi\rangle\sim\hat{\mathcal{I}}(u)\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
with}\leavevmode\nobreak\ u\equiv\frac{k_{3}}{k_{1}+k_{2}}\leavevmode\nobreak\
.$ (2.8)

Without doing the nested time integrals, the bootstrap approach is to identify
that from spacetime symmetries the scalar seed function $\hat{\mathcal{I}}(u)$
satisfies a boundary differential equation as follows

$\left(\Delta_{u}+\frac{m^{2}}{H^{2}}-2\right)\hat{\mathcal{I}}(u)=\frac{u}{1+u}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
with}\leavevmode\nobreak\
\Delta_{u}\equiv(1-u^{2})u^{2}\partial_{u}^{2}-2u^{3}\partial_{u}\leavevmode\nobreak\
.$ (2.9)

By solving this equation with proper boundary conditions, we fully fix the
analytical form of the scalar seed. In general, it contains two parts: one is
the homogeneous solution of (2.9) given by hypergeometric functions
${}_{2}F_{1}$; the other is the inhomogeneous part given by a power series.
Although this solution gives us the full analytical control of three-point
scalar seed $\langle\varphi^{2}\phi\rangle$, its form is rather complicated.
For the purpose of deriving simple bispectrum templates for data analysis,
here we propose the following approximated scalar seed

$\hat{\mathcal{I}}_{a}(u)=\frac{u}{\beta}{(1+u)^{-\frac{\beta}{\beta+2}}}-\frac{1}{2}\sum_{\pm}B_{\pm}\left(\frac{u}{2}\right)^{\frac{1}{2}\pm
i\mu}+\tilde{f}(u)\leavevmode\nobreak\ ,$ (2.10)

with

$\beta=\mu^{2}+\frac{1}{4}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
B_{\pm}=\frac{\pi^{3/2}}{\cosh\pi\mu}\left(1\mp\frac{i}{\sinh\pi\mu}\right)\frac{\Gamma(\frac{1}{2}\pm
i\mu)}{\Gamma(1\pm i\mu)}\leavevmode\nobreak\ .$ (2.11)

The expression of $\tilde{f}(u)$ is given in (A.8). The first two terms in
(2.10) provide a good match with the exact result for $\mu>1$, and there we
can safely take $\tilde{f}(u)=0$. The approximated scalar seed basically
captures the two major properties of the exact solution of
$\hat{\mathcal{I}}(u)$: i) the mass dependence of the non-squeezed kinematics
configuration is approximated by the first term; ii) the second term gives us
the collider signals in the squeezed limit. For $\mu>1$, we see the first term
is more dominant, while the second term is exponentially suppressed and only
contributes to the squeezed limit oscillations. Thus the final shape function
looks more like the equilateral one, while on the top of it we have small
wiggles around the corner of the squeezed configuration, as show in Figure 2.
We also notice that this simplification becomes less accurate when the scalar
mass gets close to $3H/2$ (i.e. $0<\mu<1$), and we need to introduce the last
term $\tilde{f}(u)$ to have a good approximation. We leave the subtleties of
this small mass regime and more details of $\tilde{f}(u)$ in Appendix A.

(a) $\mu=0.5$

(b) $\mu=1$

(c) $\mu=1.5$

(d) $\mu=3$

Figure 2: Equilateral-like shapes from massive scalar exchanges with different
masses. The blue one is the standard equilateral shape generated by the single
field EFT operator $\dot{\phi}^{3}$.

Then the scalar bispectrum from massive exchange can be derived from the
scalar seed using the weight shifting procedure

$\langle\phi_{\bf k_{1}}\phi_{\bf k_{2}}\phi_{\bf
k_{3}}\rangle^{\prime}\sim\frac{1}{k_{1}^{2}k_{2}^{2}}\mathcal{W}_{12}\left[k_{3}\hat{\mathcal{I}}(u)\right]+{\rm
perms}\leavevmode\nobreak\ ,$ (2.12)

where $\mathcal{W}_{12}$ are differential operators associated with the form
of the cubic vertices. The most general form of the weight-shifting operators
for boost-breaking interactions is given in [59], using which we can derive
the scalar exchange bispectrum from vertices with any number of derivatives.
Here we just focus on the following two simplest examples

$\displaystyle\mathcal{W}^{\dot{\phi}^{2}\sigma}_{12}$ $\displaystyle=$
$\displaystyle-k_{1}k_{2}\partial_{k_{12}}^{2}\leavevmode\nobreak\ ,$ (2.13)
$\displaystyle\mathcal{W}^{(\partial_{i}\phi)^{2}\sigma}_{12}$
$\displaystyle=$
$\displaystyle-\frac{1}{2k_{1}k_{2}}(k_{3}^{2}-k_{1}^{2}-k_{2}^{2}){\left(1-k_{1}\partial_{k_{1}}\right)\left(1-k_{2}\partial_{k_{2}}\right)}\leavevmode\nobreak\
.$ (2.14)

Then substituting the approximated scalar seed (2.10) into (2.12), we find the
simplified shape templates from massive scalar exchanges. The result from the
cubic vertex $\dot{\phi}^{2}\sigma$ is given by

$\displaystyle S_{\rm col.}^{\rm scalar-I}$ $\displaystyle=$
$\displaystyle\frac{2k_{1}k_{2}k_{3}}{\beta
k_{T}^{2}(k_{1}+k_{2})}\left[1+\frac{4k_{3}}{(\beta+2)(k_{1}+k_{2})}+\frac{(\beta+4)k_{3}^{2}}{(\beta+2)^{2}(k_{1}+k_{2})^{2}}\right]\left(\frac{k_{T}}{k_{1}+k_{2}}\right)^{-\frac{\beta}{\beta+2}}+\tilde{S}_{I}$
$\displaystyle+\frac{k_{1}k_{2}}{(k_{1}+k_{2})^{2}}\sqrt{\frac{\pi^{3}\beta(\beta+2)}{\mu\sinh(2\pi\mu)}}\left(\frac{k_{3}}{k_{1}+k_{2}}\right)^{1/2}\cos\left[\mu\log\left(\frac{k_{3}}{2(k_{1}+k_{2})}\right)+\delta\right]+2\leavevmode\nobreak\
{\rm perms.}\leavevmode\nobreak\ ,$

with

$\delta=\arg\left[\Gamma\left(\frac{5}{2}+i\mu\right)\Gamma\left(-i\mu\right)(1+i\sinh\pi\mu)\right]\leavevmode\nobreak\
.$ (2.16)

Examples of this shape template are shown in Figure 2. The second line in
(2.2) represents the typical oscillatory signals in the squeezed limit, while
the first line gives us an equilateral-like shape that is deformed by the mass
of the $\sigma$ field. There the second term $\tilde{S}_{I}\equiv-
k_{1}k_{2}\partial_{k_{1}}\partial_{k_{2}}\tilde{f}$ is introduced for the
small mass regime $\mu<1$ with its expression given in (• ‣ A). In the heavy
field regime $\mu>1$, the overall shape is dominated by the non-oscillatory
part, which resembles the standard equilateral shape, and simply reproduces
$S\propto k_{1}k_{2}k_{3}/k_{T}^{3}$ for the $m\gg H$ limit. On the top of it,
the second line in (2.2) gives the cosmological collider oscillations around
the squeezed configuration. The situation becomes different when the mass gets
close to $3H/2$ (or $\mu\leq 1$). There, the oscillatory part and the non-
oscillatory part are comparable with each other, and the overall shape could
show large deviations from the equilateral one.

(a) $\mu=1$

(b) $\mu=2$

(c) $\mu=5$

(d) $\mu=10$

Figure 3: Orthogonal-like shapes from massive scalar exchanges with different
masses. The blue one is the orthogonal shape generated by a linear combination
of the single field EFT operators $\dot{\phi}^{3}$ and
$\dot{\phi}(\partial_{i}\phi)^{2}$.

Similarly, using the weight-shifting operator (2.14), we derive the
approximated template from another cubic vertex
$(\partial_{i}\phi)^{2}\sigma$.

$\displaystyle S_{\rm col.}^{\rm II}$ $\displaystyle=$
$\displaystyle\frac{k_{3}(k_{3}^{2}-k_{1}^{2}-k_{2}^{2})}{\beta(k_{1}+k_{2})^{3}}\left[6-\frac{6\beta
k_{3}}{(\beta+2)k_{T}}+\frac{2\beta(\beta+1)k_{3}^{2}}{(\beta+2)^{2}k_{T}^{2}}+\frac{k_{1}^{2}+k_{2}^{2}}{k_{1}k_{2}}\left(2-\frac{\beta
k_{3}}{(\beta+2)k_{T}}\right)\right]\left(\frac{k_{T}}{k_{1}+k_{2}}\right)^{-\frac{\beta}{\beta+2}}$
(2.17)
$\displaystyle+\frac{k_{3}^{2}-k_{1}^{2}-k_{2}^{2}}{k_{1}k_{2}}\left(\frac{k_{3}}{k_{1}+k_{2}}\right)^{1/2}\left\\{\sqrt{\frac{\pi^{3}(\beta+2)}{\mu\sinh(2\pi\mu)}}\cos\left[\mu\log\left(\frac{k_{3}}{2(k_{1}+k_{2})}\right)+\delta_{1}\right]\right.$
$\displaystyle\left.+\frac{k_{1}k_{2}}{(k_{1}+k_{2})^{2}}\sqrt{\frac{\pi^{3}\beta(\beta+2)}{\mu\sinh(2\pi\mu)}}\cos\left[\mu\log\left(\frac{k_{3}}{2(k_{1}+k_{2})}\right)+\delta_{2}\right]\right\\}+\tilde{S}_{II}+2\leavevmode\nobreak\
{\rm perms.}\leavevmode\nobreak\ ,$

with

$\displaystyle\delta_{1}$ $\displaystyle=$
$\displaystyle\arg\left[\left(\Gamma\left(\frac{1}{2}+i\mu\right)+\Gamma\left(\frac{3}{2}+i\mu\right)\right)\Gamma\left(-i\mu\right)\left(i+{\frac{1}{\sinh\pi\mu}}\right)\right]\leavevmode\nobreak\
,$ (2.18) $\displaystyle\delta_{2}$ $\displaystyle=$
$\displaystyle\arg\left[\Gamma\left(\frac{5}{2}+i\mu\right)\Gamma\left(-i\mu\right)\left(i+{\frac{1}{\sinh\pi\mu}}\right)\right]\leavevmode\nobreak\
.$ (2.19)

The expression of $\tilde{S}_{II}$ is given in (• ‣ A). One can easily check
that for $m\gg H$ the shape function goes back to the single field one from
the EFT operator $\dot{\phi}(\partial_{i}\phi)^{2}$. Meanwhile, recall that in
single field inflation, although the two EFT shapes from $\dot{\phi}^{3}$ and
$\dot{\phi}(\partial_{i}\phi)^{2}$ look similar with each other, a linear
combination of the two leads to the new orthogonal shape, which is peaked both
at equilateral and folded configurations [90]. Similarly here, for
cosmological colliders, a linear combination of (2.2) and (2.17) gives us a
new shape that resembles the orthogonal one but has squeezed-limit
oscillations

$S_{\rm col.}^{\rm scalar-II}=-S_{\rm col.}^{\rm scalar-I}-0.14S_{\rm
col.}^{\rm II}\leavevmode\nobreak\ .$ (2.20)

Examples of this template are shown in Figure 3. Like in single field
inflation, (2.2) and (2.20) provide a good basis to capture various
possibilities of cosmological colliders from massive scalar exchanges.

The above analysis shows that, for masses with $\mu>1$, the cosmological
collider signals from scalar exchange are small corrections to the standard
equilateral or orthogonal shapes, which in general would serve as background
signals of PNG. This is expected, since the oscillations in the squeezed limit
are always suppressed by a Boltzmann factor $e^{-\pi\mu}$. This is simply due
to the fact that when particles are heavier than the Hubble scale, it becomes
more difficult to produce them during inflation.

As a final remark, we notice that the bootstrap analysis requires the full
shapes from massive exchanges to have both the oscillatory and background
parts. While for the purpose of preparing simple shape templates, another
possible approach is to consider the oscillatory parts only. As we have seen,
for $\mu>1$, the non-oscillatory parts mimic the bispectrum shapes from single
field inflation. Thus, we may consider certain combinations of massive
exchange and contact interactions to cancel the non-oscillatory background and
make the collider signal more manifest. For instance, in (2.2), this treatment
would leave us with the second line there as the full template, which
coincides with the equilateral collider shape in Section 2.4. This approach
may provide simpler and unified templates of cosmological colliders for
observational tests, which we shall further examine in a follow-up work. In
this paper, we will mainly take the full analytical shapes (“background +
oscillations”) from the bootstrap computation for the data anlaysis.

###  2.3 Shapes from spinning exchange

Another typical signal of cosmological collider is the angular dependence from
the exchange of spinning fields [6, 12, 45, 59] (see also [95] for an early
relevant work). This distinctive signature in the bispectrum can be seen as a
consequence of cubic interactions between the inflaton and spinning particles.
Two specific examples are

$\dot{\phi}\partial_{i_{1}i_{2}...i_{s}}\phi\sigma_{i_{1}i_{2}...i_{s}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\nabla_{\mu}\phi\nabla^{\mu}\nabla_{\mu_{1}\mu_{2}...\mu_{s}}\phi\sigma_{\mu_{1}\mu_{2}...\mu_{s}}\leavevmode\nobreak\
,$ (2.21)

where the first one comes from the boost-breaking EFT [12, 59], and the second
one is dS-invariant [6, 45]. To have an exchange bispectrum, we also need
linear mixings between the inflaton and spinning field, such as
$\partial_{i_{1}i_{2}...i_{s}}\phi\sigma_{i_{1}i_{2}...i_{s}}$. As a
consequence, only the longitudinal mode of the spin-$s$ particle would
contribute to the bispectrum, while other components with nonzero helicities
are projected out by the linear mixing. Therefore, we can simply assign the
spin-$s$ field with a fixed tensor structure
$\epsilon^{L}_{i_{1}i_{2}...i_{s}}$, which through the cubic vertices in
(2.21) leads to the universal angular-dependent profile for the bispectrum

$k_{2}^{i_{1}}k_{2}^{i_{2}}...k_{2}^{i_{s}}\epsilon^{L}_{i_{1}i_{2}...i_{s}}({\bf
k}_{3})\sim P_{s}(\hat{\bf k}_{2}\cdot\hat{\bf k}_{3})\leavevmode\nobreak\ ,$
(2.22)

where $P_{s}$ is the Legendre polynomial as a function of the angle between
${\bf k}_{2}$ and ${\bf k}_{3}$. As a result, the squeezed bispectrum in (2.7)
is generalized for spinning exchange

$\lim_{k_{l}\ll k_{s}}S^{{\rm spin}-s}\sim P_{s}(\hat{\bf k}_{l}\cdot\hat{\bf
k}_{s})\left(\frac{k_{l}}{k_{s}}\right)^{1/2}\cos\left[\mu\ln\left(\frac{k_{l}}{k_{s}}\right)+\delta\right]\leavevmode\nobreak\
,$ (2.23)

whose angular dependence measures the spin of the massive particles during
inflation. Away from the squeezed limit, the full shape expression for any
kinematics is complicated as one can imagine. Again, in order to simplify the
data analysis, we are motivated to look for approximated templates that have
similar angular-dependent and oscillatory signals.

#### Heavy-spin exchange template

A simple approach is to consider the heavy field limit of the spin-$s$
particle $m_{s}\gg H$, where the intermediate state can be integrated out and
we get a single field EFT. In this limit, the oscillations in the squeezed
bispectrum become highly suppressed, but the angular dependence remains. We
get the simple form of the heavy-spin exchange shape function in terms of
rational polynomials [86]

$S^{{\rm spin}-s}(k_{1},k_{2},k_{3})=P_{s}(\hat{\bf k}_{1}\cdot\hat{\bf
k}_{3})\frac{k_{2}}{(k_{1}k_{3})^{1-s}k_{T}^{2s+1}}\left[(2s-1)\left((k_{1}+k_{3})k_{T}+2sk_{1}k_{3}\right)+k_{T}^{2}\right]+{\rm
5\leavevmode\nobreak\ perms.}\leavevmode\nobreak\ ,$ (2.24)

which can be seen as the equilateral shape with an angular-dependent profile.
However, we notice that these shapes may also be predictions of single field
inflation with higher-derivative interactions, such as
$\dot{\phi}\partial_{i_{1}i_{2}...i_{s}}\phi\hat{\partial}_{i_{1}i_{2}...i_{s}}\phi$,
with $\hat{\partial}$ meaning that we only take the traceless part.

(a) $\mu=2$

(b) $\mu=3$

(c) $\mu=4$

(d) $\mu=5$

Figure 4: Angular-dependent shapes from spin-2 exchange template (2.27) with
different masses.

#### Massive spin-2 exchange template

To truly identify the distinctive signature of massive spinning particle, we
still need to consider the $m\sim H$ regime and construct an approximated
template for the shape functions there. As we have done for scalar exchanges,
here we follow the bootstrap approach in Ref. [59] and derive the template of
massive spin-2 exchange as an explicit example. The starting point is a
generalized version of the approximated scalar seed
$\hat{\mathcal{I}}_{a}^{(n)}(u)$, whose definition is given in (A.12), and we
leave more details in Appendix A. Then the analytical templates of spin-$s$
exchanges can be obtained by using two types differential operators

$S^{{\rm spin}-s}_{\rm col.}(k_{1},k_{2},k_{3})=P_{s}(\hat{\bf
k}_{2}\cdot\hat{\bf k}_{3}){k_{2}^{s-1}k_{3}^{-s}}\leavevmode\nobreak\
\mathcal{W}^{s}_{12}\leavevmode\nobreak\
\mathcal{D}^{(s)}_{23}\leavevmode\nobreak\
\left[k_{3}\hat{\mathcal{I}}^{(s)}\right]+{\rm 5\leavevmode\nobreak\
perms.}\leavevmode\nobreak\ ,$ (2.25)

where $\mathcal{W}^{s}_{12}=k_{1}\partial_{k_{1}}(1-k_{2}\partial_{k_{2}})$ is
another weight-shifting operator and $\mathcal{D}^{(s)}_{23}$ is the spin-
raising operator. For spin-2 exchanges, we have

$\mathcal{D}_{23}^{(2)}\equiv
k_{3}^{3}\left[\left(\partial_{k_{3}}-\frac{2}{k_{3}}\right)^{2}\partial_{k_{2}}-\frac{1}{3c_{s}^{2}}\partial_{k_{2}}^{3}\right]\leavevmode\nobreak\
,$ (2.26)

which leads to the following bispectrum template

$\displaystyle S_{\rm col.}^{s=2}$ $\displaystyle=$ $\displaystyle
P_{2}(\hat{\bf k}_{2}\cdot\hat{\bf
k}_{3})\frac{1991k_{1}k_{2}k_{3}^{3}}{(\beta+6)^{2}k_{T}^{8}}\left[(k_{1}+k_{2})^{2}(k_{1}+6k_{2})-(k_{1}+k_{2})(7k_{1}+57k_{2})k_{3}+(-7k_{1}+8k_{2})k_{3}^{2}+k_{3}^{3}\right]$
(2.27) $\displaystyle+P_{2}(\hat{\bf k}_{2}\cdot\hat{\bf
k}_{3})\frac{k_{1}k_{2}}{8(k_{1}+k_{2})^{3}}\sqrt{\frac{\pi^{3}\left[1+\sinh^{2}(\pi\mu)\right]}{2\mu\tanh(\pi\mu)}}\frac{\beta^{2}(\beta+2)^{2}}{\cosh^{2}(\pi\mu)}\left(\frac{k_{3}}{k_{1}+k_{2}}\right)^{1/2}$
$\displaystyle\times\left\\{k_{1}\cos\left[\mu\log\left(\frac{k_{3}}{2(k_{1}+k_{2})}\right)+\delta_{1}^{s=2}\right]+\sqrt{\beta+12}k_{2}\cos\left[\mu\log\left(\frac{k_{3}}{2(k_{1}+k_{2})}\right)+\delta_{2}^{s=2}\right]\right\\}$
$\displaystyle-P_{2}(\hat{\bf k}_{2}\cdot\hat{\bf
k}_{3})\frac{k_{1}k_{2}}{6(k_{1}+k_{2})^{3}}\sqrt{\frac{\pi^{3}\beta(\beta+2)(\beta+6)(\beta+12)}{\mu\tanh(\pi\mu)}}\frac{\beta(\beta+2)}{\cosh(\pi\mu)}\left(\frac{k_{3}}{2k_{1}+2k_{2}}\right)^{5/2}$
$\displaystyle\times\left\\{k_{1}\cos\left[\mu\log\left(\frac{k_{3}}{2(k_{1}+k_{2})}\right)+\delta_{3}^{s=2}\right]+\sqrt{\beta+30}k_{2}\cos\left[\mu\log\left(\frac{k_{3}}{2(k_{1}+k_{2})}\right)+\delta_{4}^{s=2}\right]\right\\}$
$\displaystyle+5\leavevmode\nobreak\ {\rm perms.}\leavevmode\nobreak\ ,$

with

$\displaystyle\delta_{1}^{s=2}$ $\displaystyle=$
$\displaystyle\arg\left[\frac{\Gamma\left(-i\mu\right)}{\Gamma\left(\frac{1}{2}-i\mu\right)}\left(1+{i\sinh\pi\mu}\right)\right]\leavevmode\nobreak\
,$ (2.28) $\displaystyle\delta_{2}^{s=2}$ $\displaystyle=$
$\displaystyle\arg\left[\frac{\Gamma\left(-i\mu\right)}{\Gamma\left(\frac{1}{2}-i\mu\right)}\left(1+{i\sinh\pi\mu}\right)(7+2i\mu)\right]\leavevmode\nobreak\
,$ (2.29) $\displaystyle\delta_{3}^{s=2}$ $\displaystyle=$
$\displaystyle\arg\left[{\Gamma\left(-i\mu\right)}{\Gamma\left(\frac{9}{2}+i\mu\right)}\left(1+{i\sinh\pi\mu}\right)\right]\leavevmode\nobreak\
,$ (2.30) $\displaystyle\delta_{4}^{s=2}$ $\displaystyle=$
$\displaystyle\arg\left[{\Gamma\left(-i\mu\right)}{\Gamma\left(\frac{9}{2}+i\mu\right)}\left(1+{i\sinh\pi\mu}\right)\left(\frac{11}{2}+i\mu\right)\right]\leavevmode\nobreak\
,$ (2.31)

Although the expression is quite complicated, we can see that the first line
here gives us the non-oscillatory part, while the rest is the collider signals
in the squeezed limit. Several examples are shown in Figure 4. Interestingly,
in addition to the angular dependence that changes the overall profile, the
squeezed limit oscillations become more manifested in the spinning exchanges.
This is because the non-oscillatory part of the template is less dominant in
(2.27), and the amplitudes of oscillations are comparably not-very-small
before getting exponentially suppressed for $\mu\gg 1$.

###  2.4 Shapes with small sound speeds

While the proposal of cosmological collider physics was first realized in
(nearly) dS-invariant theories, the size of the bispectrum there is in general
suppressed by slow-roll parameters, which makes it extremely difficult to test
these non-Gaussianity predictions. Then as in single field inflation, one of
the most natural ways to enhance the signal is to apply the effective field
theory approach, where dS boost symmetries can be broken and then small sound
speeds of perturbations generically lead to detectably large non-Gaussianity.

More interestingly, the recent bootstrap analysis of these boost-breaking
scenarios has identified new types of bispectrum shapes due to the presence of
multiple sound speeds. In general, it is possible for the inflaton $\phi$ and
the mediator field $\sigma$ to have different sound speeds: $c_{s}$ and
$c_{\sigma}$. Then with a large hierarchy between the two, intuitively
speaking, the exchange process can probe multiple sound horizon crossings
before the $\sigma$ field decays into the inflaton. As a result, the inflaton
leg with linear mixing may have a different clock from the other inflaton
legs. This in the end modulates the non-Gaussian shapes of the primordial
bispectrum. In the following, we briefly review three types of new shape
functions from cosmological colliders with nontrivial sound speed effects.

#### Equilateral collider

In conventional scenarios, the oscillatory signals mainly appear around the
squeezed corner of the full shape. This is because only there the oscillatory
parts of the shape function dominates over the equilateral-type contributions,
as we can see by comparing the scalings in (2.7) and (2.23) with the one in
(2.3). In the full analytical results from the bootstrap analysis, the
expansion parameter we use for the squeezed limit is $u\equiv
k_{3}/(k_{1}+k_{2})\ll 1$. For scenarios with two sound speeds, this parameter
is extended to be $u\equiv c_{\sigma}k_{3}/(c_{s}k_{1}+c_{s}k_{2})$, which
basically reflects the fact that the linear mixing leg $k_{3}$ inherits the
sound speed of the intermediate $\sigma$ field. For $c_{\sigma}\ll c_{s}$, one
interesting consequence is that we do not need to go to the actual squeezed
triangle configuration to have $u\ll 1$, and therefore, the oscillatory
collider signals would be able to appear even for equilateral configurations
with $k_{1}\simeq k_{2}\simeq k_{3}$. This is dubbed the equilateral collider
shape [59], and one particular template from the cubic interaction
$\dot{\phi}^{2}\sigma$ is given by

$S^{\rm
eq.col.}(k_{1},k_{2},k_{3})=\frac{k_{1}k_{2}}{(k_{1}+k_{2})^{2}}\left(\frac{k_{3}}{k_{1}+k_{2}}\right)^{1/2}\cos\left[\mu\log\left(\frac{c_{\sigma}k_{3}}{2c_{s}(k_{1}+k_{2})}\right)+\delta^{\rm
eq.col.}\right]+{\rm 2\leavevmode\nobreak\ perms.}\leavevmode\nobreak\ .$
(2.32)

This shape function is similar with the oscillatory part of the scalar
exchange template I (the second line in (2.2)), though here we have a
nontrivial sound speed ratio $c_{\sigma}/c_{s}\ll 1$, which is important for
the realization of this new shape. Meanwhile, we also notice that there is a
degeneracy between the sound speed ratio and the phase parameter $\delta^{\rm
eq.col.}$. Thus for practical purpose, once we have shifted these oscillations
outside of the squeezed limit, we can absorb $\mu\log(c_{\sigma}/(2c_{s}))$
into the phase parameter, to further simplify the template expression.

In addition to the template in (2.32), there are certainly various
possibilities for equilateral colliders, such as the ones from other types of
scalar interactions, and also from spinning exchanges with angular-dependent
profiles. In this work, we use this simplest example in the data analysis for
demonstration, and we leave a more comprehensive investigation for future
studies.

#### Low-speed collider

Once we turn on the hierarchy between the inflaton and the massive field sound
speeds, another possibility is $c_{s}\ll c_{\sigma}$. This is motivated by the
consideration in the EFT of inflation that the boost-breaking operators can
lead to a very small sound speed for the inflaton fluctuations and naturally
enhance the non-Gaussianity signal. Without loss of generality, here we can
take $c_{\sigma}=1$. The bootstrap analysis of cosmological colliders with
this setup has been performed in detail in Ref. [60], with a new bispectrum
shape named as the low-speed collider. Here we simply recap the physical
picture and the non-Gaussian phenomenology.

For a two-field system during inflation, in addition to the sound horizon
crossing of the inflaton at $k=aH/c_{s}$, there is also a mass horizon
crossing for the massive field at $k=am$. Now we consider the massive field
travels much faster than the inflaton fluctuations. When $m<H/c_{s}$, it
becomes possible that these two types of horizon crossings for different legs
would happen at the coincident time. As a result, a resonance enhancement may
happen around the mildly squeezed configuration $k_{l}/k_{s}\simeq c_{s}m/H$,
which generates a characteristic bump in the shape function. In this
situation, the oscillatory collider signal is further shifted deep into the
squeezed limit, while the resonance peak in the mildly squeezed regime may
become more dominant than the equilateral-type background part of the
bispectrum. A simple template for the low-speed collider was proposed in Ref.
[41]

$S^{{\rm LSC}}(k_{1},k_{2},k_{3})=S^{\rm
eq}(k_{1},k_{2},k_{3})+\frac{k_{1}^{2}}{3k_{2}k_{3}}\left[1+\alpha\left(\frac{k_{1}^{2}}{3k_{2}k_{3}}\right)^{2}\right]^{-1}+{\rm
2\leavevmode\nobreak\ perms.}\leavevmode\nobreak\ ,$ (2.33)

with $\alpha=c_{s}m/H$ and $S^{\rm eq}$ being the equilateral template used by
WMAP and Planck. By varying the parameter $\alpha$, the location of the low-
speed collider bump will be shifted. We present the data analysis in Section
4.3 and the constraints in Figure 15.

#### Multi-speed non-Gaussianity

This type of new phenomenology is identified in the massless exchange diagram
in Figure 1c), and can be seen as one prediction of multi-field inflation that
differs from the local shape. For multi-field models, the linear mixing
$\dot{\phi}\sigma$ in general describes the usual isocurvature-to-curvature
conversion process, and we expect infrared divergences in the massless
exchange diagrams. This provides a field-theoretical explanation for the
generation of local-type non-Gaussianities. However, if the linear mixing is
given by higher derivative interaction, such as $\dot{\phi}\dot{\sigma}$ or
$\partial^{2}\phi\sigma$, the infrared divergences disappear, and the mixed
propagators can inherit the sound speed of the intermediate field. As a
result, for massless exchanges with multiple internal lines, each legs would
have a different time of sound horizon crossing, controlled by their own sound
speed. Similar with the low-speed collider, the resonance enhancement for the
shape function may appear for arbitrary kinematic configurations depending on
the sound speed ratios. More details about this new shape can be found in Ref.
[64]. Here we use the following template333See also [96] for an earlier work
with a similar shape function but from a different setup.

${S}^{{\rm
multi}-c_{s}}(k_{1},k_{2},k_{3})=\frac{k_{1}k_{2}k_{3}}{(c_{1}k_{1}+c_{2}k_{2}+c_{3}k_{3})^{3}}+5\leavevmode\nobreak\
{\rm perms}\leavevmode\nobreak\ ,$ (2.34)

where $c_{1,2,3}$ are different sound speed parameters. By tuning their
ratios, we are changing the location of the peak in the bispectrum shape. The
data analysis on this template is presented in Section 4.3, with results shown
in Figure 16.

##  3 CMB Bispectrum Estimator

In this section, we briefly review the current status of the CMB bispectrum
analysis, and introduce the CMB-BEST python pipeline [94]. This new high-
resolution estimator uses a flexible set of modal bases for the primordial
bispectra, and optimizes the algorithm to significantly improve the
computational efficiency. Particularly, it provides us with an ideal tool for
testing PNG signals with complicated shape functions in the Planck temperature
and polarization data.

###  3.1 CMB bispectrum statistics

If the primordial fluctuations are non-Gaussian, the two-point function (power
spectrum) no longer captures the full statistical information, and higher-
order correlation functions can be significant. As the CMB anisotropies relate
mostly linearly to the primordial perturbations, PNG is captured in higher-
order statistics of the CMB.

#### CMB bispectrum

Given the observed CMB temperature anisotropies $a_{\ell m}$, the bispectrum
statistic is constructed as the three-point function of the s as

$\displaystyle\hat{B}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}\equiv
a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}a_{\ell_{3}m_{3}}-\left[\left<a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}\right>_{\mathrm{MC}}a_{\ell_{3}m_{3}}+(2\
\mathrm{perms})\right],$ (3.1)

where the terms in brackets account for the anisotropic sky coverage and
masking effects, evaluated through ensemble averages (Monte-Carlo) from
realistic simulations.

The expected value of $\hat{B}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}$
is zero if there is no PNG. Otherwise, given a primordial bispectrum shape
$S(k_{1},k_{2},k_{3};\bm{\theta})$ with some free parameters $\bm{\theta}$,

$\displaystyle\langle\hat{B}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}\rangle$
$\displaystyle={f_{\mathrm{NL}}}\
\mathcal{G}^{{\ell_{1}\ell_{2}\ell_{3}}}_{m_{1}m_{2}m_{3}}\frac{144P_{\zeta}^{2}}{5\pi^{3}}\int
r^{2}drdk_{1}dk_{2}dk_{3}\;S(k_{1},k_{2},k_{3};\bm{\theta})\prod_{j=1}^{3}\left[j_{\ell_{j}}(k_{j}r)T_{\ell_{j}}(k_{j})\right]$
(3.2) $\displaystyle={f_{\mathrm{NL}}}\
\mathcal{G}^{{\ell_{1}\ell_{2}\ell_{3}}}_{m_{1}m_{2}m_{3}}b_{{\ell_{1}\ell_{2}\ell_{3}}}(\bm{\theta}).$
(3.3)

The Gaunt integral $\mathcal{G}$ is a geometric
factor444$\mathcal{G}^{{\ell_{1}\ell_{2}\ell_{3}}}_{m_{1}m_{2}m_{3}}\equiv\int
d^{2}\hat{\mathbf{n}}\
Y_{\ell_{1}m_{1}}(\hat{\mathbf{n}})Y_{\ell_{2}m_{2}}(\hat{\mathbf{n}})Y_{\ell_{3}m_{3}}(\hat{\mathbf{n}})$,
where $Y_{\ell m}$s denote spherical harmonics. that enforces the angular
momentum conservation. The CMB transfer function $T_{\ell}(k)$ and spherical
Bessel function $j_{\ell}(x)$ capture the evolution and projection of the
primordial perturbations to the CMB anisotropies we observe today. The reduced
bispectrum $b_{{\ell_{1}\ell_{2}\ell_{3}}}(\bm{\theta})$ contains all
observable signatures of the theory in the CMB bispectrum, up to the overall
amplitude ${f_{\mathrm{NL}}}$.

#### Likelihood

For CMB bispectrum analysis, we assume the weakly non-Gaussian regime where
$\hat{B}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}$s are multivariate
Gaussian with variance proportional to $C_{\ell_{1}}C_{\ell_{2}}C_{\ell_{3}}$
given by Wick’s theorem. The CMB bispectrum likelihood is then given by

$\displaystyle\mathcal{L}({f_{\mathrm{NL}}},\bm{\theta})$
$\displaystyle\propto\mathrm{exp}\left[-\frac{1}{2}\sum_{\ell_{j},m_{j}}\frac{1}{6\;C_{\ell_{1}}C_{\ell_{2}}C_{\ell_{3}}}\left|\hat{B}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}-{f_{\mathrm{NL}}}\
\mathcal{G}^{{\ell_{1}\ell_{2}\ell_{3}}}_{m_{1}m_{2}m_{3}}b_{{\ell_{1}\ell_{2}\ell_{3}}}(\bm{\theta})\right|^{2}\right]$
(3.4)
$\displaystyle\equiv\mathrm{exp}\left[(\mathrm{const.})+{f_{\mathrm{NL}}}P(\bm{\theta})-\frac{1}{2}f^{2}_{\mathrm{NL}}F(\bm{\theta})\right].$
(3.5)

The function $P(\bm{\theta})$ and the Fisher information $F(\bm{\theta})$ are
linear and quadratic in $b_{\ell_{1}\ell_{2}\ell_{3}}(\bm{\theta})$,
respectively.

If the CMB bispectrum had been measured with high signal-to-noise, one could
potentially perform a detailed analysis by varying both ${f_{\mathrm{NL}}}$
and $\bm{\theta}$ to constrain the model space, analogously to the CMB power
spectrum analysis where the $\Lambda$CDM model parameters are measured.
However, we are in a regime with very low signal-to-noise so far,
unfortunately. The CMB bispectrum analyses focus on fitting a single parameter
${f_{\mathrm{NL}}}$.

For a given theoretical template $S(\bm{\theta})$, the best-fit
${f_{\mathrm{NL}}}$ is given by the maximum likelihood estimator (MLE):

$\displaystyle\hat{f}_{\mathrm{NL}}|_{\bm{\theta}}=P(\bm{\theta})/F(\bm{\theta}).$
(3.6)

This estimator is unbiased:
$\langle\hat{f}_{\mathrm{NL}}\rangle={f_{\mathrm{NL}}}$ if the underlying
model (3.3) is true. Furthermore, it is optimal, having the smallest expected
variance amongst all unbiased estimators of ${f_{\mathrm{NL}}}$ given
$\bm{\theta}$:
$\mathrm{Var}[\hat{f}_{\mathrm{NL}}|_{\bm{\theta}}]=F(\bm{\theta})^{-1}$,
equal to the Cramer-Rao bound.

#### Normalization

The parametric form of (3.3) has a scaling redundancy: the theoretical
bispectrum remains unchanged under ${f_{\mathrm{NL}}}\rightarrow
c{f_{\mathrm{NL}}}$ and $S(k_{1},k_{2},k_{3})\rightarrow
S(k_{1},k_{2},k_{3})/c$ for some constant $c$. This directly influences the
${f_{\mathrm{NL}}}$ constraints quoted. Larger shape function amplitude gives
smaller ${f_{\mathrm{NL}}}$ estimates and seemingly tighter constraints. A
shape function template’s normalization therefore needs to be specified for
clarity.

In this work, the templates have been normalized at the equilateral limit
unless specified:

$\displaystyle S(k_{1},k_{2},k_{3})=1\quad\mathrm{at}\quad
k_{1}=k_{2}=k_{3}=k_{*}.$ (3.7)

The pivot value is fixed to $k_{*}\equiv 0.05\ \mathrm{Mpc}^{-1}$, although
this would not affect the cosmological collider templates which are scale-
invariant. This choice is equivalent to the Planck normalization convention
written in terms of the bispectrum of the gravitational potential $\Phi$:

$\displaystyle B_{\Phi}^{{f_{\mathrm{NL}}}=1}(k_{*},k_{*},k_{*})=6A^{2},$
(3.8)

where $A\equiv P_{\Phi}(k_{*})$ is the primordial $\Phi$ power spectrum
amplitude.

While this normalization convention provides a clear reference point to
compare the detectability of various bispectrum shapes, it is not always the
most fair one. In this convention, shapes with a relatively small equilateral
limit and larger squeezed/flattened limits often have larger amplitudes, which
lead to much tighter ${f_{\mathrm{NL}}}$ constraints. In such cases, we direct
the readers to focus on the significance
${f_{\mathrm{NL}}}/\sigma({f_{\mathrm{NL}}})$—the ratio of the estimated
${f_{\mathrm{NL}}}$ and its uncertainty—as it is insensitive to the
normalization.

#### Shape correlation

The shape and $\bm{\theta}$ are fixed in the ${f_{\mathrm{NL}}}$ estimation
shown in (3.6). Independent studies of different shapes will hence provide
different estimates of ${f_{\mathrm{NL}}}$ from the observed CMB maps.
However, the individual estimates are correlated. Suppose that two shapes
$S^{(1)}$ and $S^{(2)}$ give two estimates $\hat{f}_{\mathrm{NL}}^{(1)}$ and
$\hat{f}_{\mathrm{NL}}^{(2)}$. Their expected covariance is given by

$\displaystyle\mathrm{Cov}\left[\hat{f}_{\mathrm{NL}}^{(1)},\hat{f}_{\mathrm{NL}}^{(2)}\right]=\frac{\langle
S^{(1)},S^{(2)}\rangle_{\mathrm{CMB}}}{\langle
S^{(1)},S^{(1)}\rangle_{\mathrm{CMB}}\ \langle
S^{(2)},S^{(2)}\rangle_{\mathrm{CMB}}},$ (3.9)

where we defined the inner product $\langle\cdot\rangle_{\mathrm{CMB}}$ as a
weighted dot product in the CMB bispectrum space:

$\displaystyle\langle
S^{(1)},S^{(2)}\rangle_{\mathrm{CMB}}\equiv\sum_{\ell_{j}}\frac{h_{\ell_{1}\ell_{2}\ell_{3}}^{2}}{6\;C_{\ell_{1}}C_{\ell_{2}}C_{\ell_{3}}}b^{(1)}_{{\ell_{1}\ell_{2}\ell_{3}}}b^{(2)}_{{\ell_{1}\ell_{2}\ell_{3}}}.$
(3.10)

The reduced bispectrum $b^{(i)}_{\ell_{1}\ell_{2}\ell_{3}}$ relate to the
shape $S^{(i)}$ through (3.3). The geometric factor
$h^{2}_{\ell_{1}\ell_{2}\ell_{3}}\equiv\sum_{m_{j}}\left|\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}\right|^{2}$.
Using this notation, the Fisher information can be rewritten as
$F=\left<S,S\right>_{\mathrm{CMB}}$.

We define the correlation (or ‘cosine’) between two shapes as the Pearson
correlation coefficient of their respective $\hat{f}_{\mathrm{NL}}$s:

$\displaystyle\rho(S^{(1)},S^{(2)})\equiv\mathrm{corr}\left[\hat{f}_{\mathrm{NL}}^{(1)},\hat{f}_{\mathrm{NL}}^{(2)}\right]=\frac{\langle
S^{(1)},S^{(2)}\rangle_{\mathrm{CMB}}}{\sqrt{\langle
S^{(1)},S^{(1)}\rangle_{\mathrm{CMB}}\ \langle
S^{(2)},S^{(2)}\rangle_{\mathrm{CMB}}}}.$ (3.11)

This value lies in $[-1,1]$ by construction. A strong correlation (+1) or
anti-correlation (-1) between the two shapes implies that they probe similar
forms of the CMB bispectrum, and their ${f_{\mathrm{NL}}}$ estimates are
highly correlated. Such shapes are difficult to distinguish from the CMB
bispectrum analysis only. We also note that it is possible to have two very
distinct shape functions with a strong correlation under this metric, since
the projection from the primordial perturbations to the CMB anisotropies can
wash out the differences.

Figure 5 shows the correlation between cosmological colliders templates and
the standard local, equilateral and orthogonal templates. Note that the
templates were evaluated at their best-fit parameters; refer to Section 4 for
detailed analysis. We also take this chance to introduce the acronyms of the
new templates which will be used later. Many of the templates studied in this
work are not strongly correlated with any of the standard templates.

Figure 5: Correlation between the cosmological collider templates studied in
this work, together with the standard local, equilateral and orthogonal
templates. The correlation (defined as (3.11)) measures the similarity of
their observable CMB bispectrum predictions. Free parameters in each template
have been fixed to their best-fit values for this plot. Here we also summarize
our acronyms for templates.

###  3.2 CMB bispectrum estimator

#### Mode expansion

Direct numerical computation of (3.6) is prohibitively expensive due to the
number of terms appearing in the sum: $O(\ell_{\mathrm{max}}^{5})$ for
$\ell_{\mathrm{max}}\lesssim 2500$ for Planck. Existing formalisms utilise the
separability of the high-dimensional integral [87, 88, 89, 90, 91, 77, 78, 92]
or binning to reduce the computational complexity [79, 93]. In CMB-BEST, the
(primordial) shape function is expanded in terms of a separable basis, so that

$\displaystyle
S(\bm{\theta})(k_{1},k_{2},k_{3})=\sum_{n\leftrightarrow(p_{1},p_{2},p_{3})}\alpha_{n}(\bm{\theta})\;q_{p_{1}}(k_{1})q_{p_{2}}(k_{2})q_{p_{3}}(k_{3}),$
(3.12)

where the one-dimensional basis functions $q_{p}(k)$ are Legendre polynomials
of order up to $p_{\mathrm{max}}=30$, forming a three-dimensional basis of
size $\sim 5000$ up to symmetry. The basis function choices are flexible in
principle, but we have found that our Legendre basis accurately covers the
range of templates studied in this work.

Thanks to the separability of the terms appearing in (3.12), the three-
dimensional integrals $\int dk_{1}dk_{2}dk_{3}$ appearing in (3.3) split into
three one-dimensional integrals, dramatically reducing the computational
complexity. Under the CMB-BEST formalism, only the expansion (3.12) in the
primordial space is required, and the CMB bispectrum likelihood is given by

$\displaystyle\log\mathcal{L}({f_{\mathrm{NL}}},\bm{\theta})=\log\mathcal{L}_{0}+{f_{\mathrm{NL}}}\sum_{n}\beta_{n}\alpha_{n}(\bm{\theta})-\frac{1}{2}f^{2}_{\mathrm{NL}}\sum_{n,n^{\prime}}\gamma_{nn^{\prime}}\alpha_{n}(\bm{\theta})\alpha_{n^{\prime}}(\bm{\theta}).$
(3.13)

Here, $\mathcal{L}_{0}$ is a constant corresponding to the likelihood of
${f_{\mathrm{NL}}}=0$, or vanishing bispectrum. The quantities $\beta_{n}$ and
$\gamma_{nn^{\prime}}$ only need to be computed once per basis and CMB map.
The ones corresponding to the Legendre basis and the Planck 2018 CMB map (both
temperature-only and with polarization) are provided together with the public
CMB-BEST code. The best-fit ${f_{\mathrm{NL}}}$ can be written as

$\displaystyle\hat{f}_{\mathrm{NL}}|_{\bm{\theta}}=\frac{\sum_{n}\beta_{n}\alpha_{n}(\bm{\theta})}{\sum_{n,n^{\prime}}\gamma_{nn^{\prime}}\alpha_{n}(\bm{\theta})\alpha_{n^{\prime}}(\bm{\theta})}.$
(3.14)

We use 160 Planck full focal plane simulations with Gaussian initial
conditions for two purposes: 1) to estimate the ensemble average in the linear
term of (3.1) and 2) to compute the best-fit ${f_{\mathrm{NL}}}$s on each map
to evaluate the estimation uncertainty: $\sigma(\hat{f}_{\mathrm{NL}})$.

#### Significance

We define the significance, or the signal-to-noise ratio, as the ratio of the
estimated ${f_{\mathrm{NL}}}$ and its estimation uncertainty:

$\displaystyle{\sigma_{\mathrm{SNR}}}|_{\bm{\theta}}\equiv\left|\frac{\hat{f}_{\mathrm{NL}}|_{\bm{\theta}}}{\sigma(\hat{f}_{\mathrm{NL}}|_{\bm{\theta}}).}\right|.$
(3.15)

The larger the ${\sigma_{\mathrm{SNR}}}$, the more significant the evidence
towards a non-zero ${f_{\mathrm{NL}}}$ for the given template. Formally put,
under the null hypothesis of vanishing PNG ($H_{0}:\ {f_{\mathrm{NL}}}=0$),
the alternative hypothesis of a non-zero bispectrum signal ($H_{1}:\
{f_{\mathrm{NL}}}\neq 0$) can be tested. For a fixed value of $\theta$, the
corresponding $p$-value is given by
$p=1-\mathrm{erf}({\sigma_{\mathrm{SNR}}}/\sqrt{2})$ under the likelihood
(3.13). A test with a significance level of $\alpha=0.05$, for example,
rejects $H_{0}$ whenever ${\sigma_{\mathrm{SNR}}}>2$.

It is common to test multiple templates or a template with varying parameters
$\bm{\theta}$ simultaneously. In such cases, one obtains a set of
${\sigma_{\mathrm{SNR}}}|_{\bm{\theta}}$, which provides a higher probability
of finding at least one large number by chance. One may make a false detection
by simply ‘looking elsewhere’ until a large significance is found. Such look-
elsewhere effect must be accounted for in any extensive parametric searches
for PNG.

We follow the ideas presented in [97] and consider the statistical
distribution of the raw maximum significance
$\sigma_{\mathrm{max}}\equiv\mathrm{max}_{\bm{\theta}}\left[{\sigma_{\mathrm{SNR}}}|_{\bm{\theta}}\right]$.
This distribution is then used to compute the look-elsewhere-adjusted
significance $\tilde{\sigma}_{\mathrm{SNR}}$. More precisely, given the raw
maximum significance $\sigma_{\mathrm{max}}$, the $p$-value of the joint
hypothesis test is given by

$\displaystyle
p(\sigma_{\mathrm{max}})=1-\mathrm{Prob}\left[{\sigma_{\mathrm{SNR}}}|_{\bm{\theta}}\leq\sigma_{\mathrm{max}}\
\mathrm{for\ all}\ \bm{\theta}\right],$ (3.16)

which gives the adjusted significance level

$\displaystyle\tilde{\sigma}_{\mathrm{SNR}}\equiv\sqrt{2}\,{\mathrm{erf}}^{-1}\left(1-p(\sigma_{\mathrm{max}})\right),$
(3.17)

where $\mathrm{erf}^{-1}$ is the inverse error function.

In practice, we compute $\tilde{\sigma}_{\mathrm{SNR}}$ for a template scan
with $\bm{\theta}=\bm{\theta}_{i}$ for $i=1,\cdots,N$ as follows. Under the
null hypothesis of vanishing bispectrum, the theoretical correlation between
significances is given by (3.11):
$C_{ij}\equiv\rho\left(S(\bm{\theta}_{i}),S(\bm{\theta}_{j})\right)=F_{ij}/\sqrt{F_{ii}F_{jj}}$
(no summation assumed), where the Fisher matrix
$F_{ij}=\sum_{n,n^{\prime}}\gamma_{nn^{\prime}}\alpha_{n}(\bm{\theta}_{i})\alpha_{n^{\prime}}(\bm{\theta}_{j})$
follows from the likelihood (3.13). We create $10^{5}$ realisations of the
multivariate normal distribution $\mathcal{N}(\mathbf{0},C)$, and compute the
maximum significance for each realisation. Next, we find the fraction of
realisations that have their maximum significance greater than the
observational value $\sigma_{\mathrm{max}}$. This fraction accurately
estimates $p(\sigma_{\mathrm{max}})$, which can then be used to compute
$\tilde{\sigma}_{\mathrm{SNR}}$ through (3.17).

##  4 The Search for Collider Signals in the Planck Data

In this section, we apply the CMB-BEST for data analysis and present
constraints on cosmological colliders based on various templates introduced in
Section 2. All constraints are based on Planck 2018 temperature and
polarisation data.

The main objective of the data analysis shown here is twofold. First, we
provide explicit bounds on the ${f_{\mathrm{NL}}}$ for various templates. This
serves as a reference for theoretical model-building and can rule out models
that predict larger amplitudes. Second, we search for signatures of
cosmological colliders by studying the significance of non-zero
${f_{\mathrm{NL}}}$ for various templates. Statistically significant detection
would have important implications for early universe physics.

The shape templates are normalized so that $S(k_{*},k_{*},k_{*})=1$ unless
explicitly specified, as discussed in Section 3. The uncertainties on the
${f_{\mathrm{NL}}}$ constraints are computed from 160 Planck simulated maps
with Gaussian initial conditions (${f_{\mathrm{NL}}}=0$). All
${f_{\mathrm{NL}}}$ constraints quoted explicitly are at $68\%$ confidence
level.

###  4.1 Oscillations from massive scalars

We start with the analysis of three bispectrum templates from massive scalar
exchanges proposed in Section 2.2.

Figure 6: Left: Constraints on the quasi-single field inflation template (2.6)
for allowed values of the index ${0<\nu<3/2}$. The black solid line indicates
the best-fit ${f_{\mathrm{NL}}}$ value for each value of $\nu$ considered
independently. The 68% ($1\sigma$) and 95% ($2\sigma$) confidence levels are
shown. Right: The raw significance ${\sigma_{\mathrm{SNR}}}$ is shown as a
function of $\nu$, maxmised at $\nu=0.14$ with maximum value $0.37$.

#### Quasi-single field inflation

The non-oscillatory template (2.6) is generated from massive scalars in
complementary series, and has one free mass parameter $\nu\in[0,3/2]$ in
addition to the size of the bispectrum $f_{\rm NL}^{\rm QSF}$. The constraints
on this template at each value of $\nu$ are shown in Figure 6. The best-fit
values of ${f_{\mathrm{NL}}}$ are shown together with $1\sigma$ and $2\sigma$
confidence levels.

Note that this template has been constrained previously in the Planck 2013
analysis [75], and our constraints are consistent with these results. However,
we do not perform the same Monte-Carlo Markov Chain (MCMC) analysis as shown
in [75] and refrain from setting a prior distribution on $\nu$. Instead, the
independent constraints for each $\nu$ are summarised in Figure 6.

The constraints are given by $9\pm 26$ at $\nu=0$ and $-2\pm 12$ at $\nu=1.5$
at $68\%$ confidence level (CL hereafter). The error bars are tighter for
larger $\nu$ because our normalization convention fixes the equilateral limit;
the shapes with larger $\nu$ have larger squeezed-limit contributions to the
bispectrum, and larger bispectrum templates get tighter constraints on their
amplitudes ${f_{\mathrm{NL}}}$.

Overall, having $f_{\mathrm{NL}}^{\mathrm{QSF}}=0$ is entirely consistent with
the data. The righthand plot of Figure 6 shows the significance (signal-to-
noise) of ${f_{\mathrm{NL}}}^{\mathrm{QSF}}$ as a function of $\nu$, which
takes the maximum value of $0.37$ at $\nu=0.14$. This value is rather small.
To put this into context, we studied the theoretical distribution of the
maximum significance ($0\leq\nu\leq 1.5$) statistic. The maximum significance
is greater than or equal to 0.37 about 90% of the time (or less than 0.37
about 10% of the time). Our results are fully consistent with having Gaussian
initial conditions.

#### Massive scalar exchange

Two templates (2.2) and (2.20) correspond to the equilateral and orthogonal
types of massive scalar bispectra, respectively. Both shapes contain a single
free mass parameter $\mu$. For $\mu\gg 1$, they reproduce the standard
equilateral and orthogonal shapes because the oscillations are suppressed by
the Boltzmann factor $e^{-\pi\mu}$. In the bispectrum data analysis, we scan
through $0<\mu<6$.

Figure 7: Constraints on the scalar exchange templates (2.2) (Left) and (2.20)
(Right) for the index $\mu$ varied from $0$ to $6$. The black solid line
indicates the best-fit ${f_{\mathrm{NL}}}$ value for each value of $\mu$
considered independently, and the 68% and 95% confidence levels are shown in
shaded colours. The amplitude ${f_{\mathrm{NL}}}$ is constrained at the
$\sigma({f_{\mathrm{NL}}})\sim 10-50$ level.

Figure 7 shows the constraints for the two scalar exchange templates. On the
left-hand side is the first template (2.2). As before, the best-fit value of
${f_{\mathrm{NL}}}$ and the $1\sigma$ and $2\sigma$ error bounds are shown for
each $\mu$. The constraints are given by ${f_{\mathrm{NL}}}=11\pm 13$ at
$\mu=0$, widening to ${f_{\mathrm{NL}}}=13\pm 18$ at $\mu=1$, and converging
around ${f_{\mathrm{NL}}}=16\pm 33$ for $\mu>5$, which simply agrees with the
constraint on the equilateral shape. Vanishing PNG (${f_{\mathrm{NL}}}=0$)
lies within $1\sigma$ of the contours everywhere; the data is fully consistent
with zero bispectrum amplitude for this template.

Shown on the right-hand side of Figure 7 is the scalar exchange template II
(2.20) constraints. The shape function here approaches the orthogonal template
when we increase the mass. The constraints are ${f_{\mathrm{NL}}}=-24\pm 45$
at $\mu=0$, shifted to ${f_{\mathrm{NL}}}=-15.7\pm 49$ at $\mu=1$. The lowest
best-fit value is reached for ${f_{\mathrm{NL}}}=-111\pm 53$ at $\mu=1.89$,
and slowly changes to ${f_{\mathrm{NL}}}=-59\pm 29$ at $\mu=6$. Unlike the
equilateral-type template earlier, ${f_{\mathrm{NL}}}=0$ is outside the
$2\sigma$ contour for some parts of the parameter space around $\mu=2.1$

Figure 8 shows the significance of signals for the two scalar exchange
templates. For the first template, the significance is less than 1 for all
values of $\mu$ and is statistically consistent with random fluctuations
around $0$. We obtain the maximum significance of $0.86$ about $50\%$ of the
time in a Universe with Gaussian initial conditions, corresponding to the
look-elsewhere-effect-adjusted significance of $0.67$.

Meanwhile, the orthogonal-type template is more favoured by the Planck CMB
data. The maximum significance of $2.3$ is obtained at $\mu=2.13$, which is
adjusted to give $\tilde{\sigma}_{\mathrm{SNR}}=1.8$. This is not sufficiently
strong statistically to reject ${f_{\mathrm{NL}}}=0$, but we note that this is
the largest adjusted significance found for the cosmological collider shapes
so far.

Figure 8: Significance of the best-fit ${f_{\mathrm{NL}}}$ values for massive
scalar exchange templates (2.2) and (2.20). Left: The raw significance
${\sigma_{\mathrm{SNR}}}$ is shown as a function of $\mu$ for the two
templates. Right: The scalar-II template’s maximum significance of 2.3 at
$\mu=2.1$ is compared with the theoretical distribution of the statistic. We
find the significance as big as 2.3 $7.8\%$ of the time and the adjusted
significance of 1.8.

###  4.2 Angular dependence from spinning particles

Next, we turn to analyze the two types of bispectrum templates with angular
dependence proposed in Section 2.3.

Figure 9: Correlation between the standard local, equilateral and orthogonal
templates and the heavy-spin exchange templates (2.24). The correlation is
defined in (3.11) and measures the similarity of their observable CMB
bispectrum predictions, with $1$ and $-1$ for perfect correlation and anti-
correlation, respectively.

#### Heavy-spin exchange

The template (2.24) probes the characteristic angular dependence from spinning
particles with $m_{s}\gg H$, which can be seen as the predictions of single
field interactions after integrating out a heavy spin-$s$ field. The spin $s$
enters the template as the order of the Legendre polynomial, and higher spins
lead to more complex angular behaviours. Figure 9 shows the correlation
between the heavy-spin exchange templates and the standard local, equilateral
and orthogonal templates. Spin $0$ and $1$ cases are indistinguishable from
the local and equilateral templates with tiny differences due to small
spectral tilt, as expected from the formula (2.24). Spins $2$-$4$ show strong
to moderate correlations with the orthogonal template and each other. Spin 6
probes a relatively new shape different from the standard templates.

Spin |  ${f_{\mathrm{NL}}}$ constraint | Significance
---|---|---
0 | $-1.5\pm 5.8$ | 0.26
1 | $-8\pm 45$ | 0.18
2 | $-18\pm 10$ | 1.8
3 | $-59\pm 32$ | 1.8
4 | $-17\pm 16$ | 1.1
6 | $-1.9\pm 6.7$ | 0.28
Table 2: Constraints on the heavy-spin exchange template (2.24), for different
values of spin. Errors are quoted at $68\%$ CL. The significance
${\sigma_{\mathrm{SNR}}}$ for each template is shown.

The constraints on ${f_{\mathrm{NL}}}$ are shown in Table 2. We limit
ourselves to smaller and even spins based on Occam’s razor. We note that this
choice is somewhat arbitrary and one can in principle study larger values of
$s$. The ${f_{\mathrm{NL}}}$ errors for spins $0$ and $1$ are close to the
ones for the local and equilateral templates, but the errors tend to decrease
for higher spins. This is because we normalize the templates at the
equilateral limit, while the Legendre polynomials evaluated at
$\hat{\mathbf{k}}_{1}\cdot\hat{\mathbf{k}}_{3}=1/2$ tend to decrease with the
order $s$. The higher spin templates tend to have larger amplitudes, leading
to stronger bounds on ${f_{\mathrm{NL}}}$.

The largest significance level of $1.8\sigma$ is achieved at spin 3, closely
tracing that of the orthogonal template. After accounting for the look-
elsewhere effect from studying six templates ($s=0,1,2,3,4,6$), the adjusted
significance is around $1.2\sigma$.

#### Massive spin-2 exchange

Figure 10: Constraints on the massive spin-2 exchange template (2.27) for the
mass parameter in range ${1<\mu<8}$. Each template at $\mu$ has been
normalized based on its Fisher information so that the uncertainties are close
to $1$. Note the error bars are still estimated from Gaussian simulations and
do not exactly equal the optimal value of $1$.

Figure 11: Significance of the best-fit ${f_{\mathrm{NL}}}$ values for the
massive spin-2 exchange template. Left: The raw significance
${\sigma_{\mathrm{SNR}}}$ is shown as a function of $\mu$. Right: The maximum
significance value of 1.99 at $\mu=4.83$ is compared with the theoretical
distribution of the maximum statistic. One finds the significance as big as
1.99 roughly $29\%$ of the time.

The signature of a massive spinning particle with $m\sim H$ is captured in our
template (2.27), where we focus on the case with $s=2$. This template has both
angular-dependent profiles and squeezed-limit oscillations, which can not be
generated in single-field inflation models. We constrain this shape for a
range of the mass parameter $\mu$.

We note that the default normalization convention fails for this highly
oscillatory template, because the equilateral limit vanishes for some values
of $\mu$. While it is possible to normalize using different limits (e.g. set
$S(k_{*},k_{*},2k_{*})=1$), we found that the template’s overall amplitude
varies drastically in most conventions and creates widely varying
uncertainties. To avoid unnecessary confusion and to make a fairer comparison
of all values of $\mu$ involved, we rescale the template using the Fisher
information at each $\mu$. Under this convention, the estimation uncertainty
should equal $1$ when the ${f_{\mathrm{NL}}}$ estimator is optimal. The
rescaled constraints are shown in Figure 10. The maximum significance is
achieved at $\mu=4.83$ for which the rescaled ${f_{\mathrm{NL}}}=-2.1\pm 1.1$,
as shown in Figure 11. After accounting for the look-elsewhere effect, the
adjusted significance is 1.1.

###  4.3 Shapes with small sound speeds

#### Equilateral collider shape

The template in (2.32) depends on the mass parameter $\mu$, the sound speed
ratio $c_{s}/c_{\sigma}$, and the phase parameter $\delta^{\mathrm{eq.col.}}$.
However, as noted in section 2.4, phenomenologically there are only two free
parameters $\mu$ and $\phi_{\rm
EC}\equiv\mu\log(c_{\sigma}/c_{s})+\delta^{\mathrm{eq.col.}}$. Furthermore,
note that $S^{\mathrm{eq.col.}}(\phi_{\rm EC})=\cos\phi_{\rm EC}\
S^{\mathrm{eq.col.}}(\phi_{\rm EC}=0)+\sin\phi_{\rm EC}\
S^{\mathrm{eq.col.}}(\phi_{\rm EC}=\pi/2)$ by trigonometric identities. This
allows the CMB bispectrum constraints to be placed quickly for all $\phi_{\rm
EC}$ with only two evaluations required at $\phi_{\rm EC}=0$ and $\pi/2$.

Figure 12: Left: Constraints on the equilateral collider template (2.32) for
different mass parameters, when the phase parameter $\phi$ is fixed to the
respective best-fit values. The black solid line indicates the best-fit
${f_{\mathrm{NL}}}$ value at each $\mu$, while the contours represent the 68%
($1\sigma$) and 95% ($2\sigma$) confidence levels. Right: The maximum
significance of the best-fit ${f_{\mathrm{NL}}}$ value for the equilateral
collider template, compared with the theoretical distribution. One finds the
significance as big as 2.5 roughly $37\%$ of the time, so this is fully
consistent with random fluctuation.  Figure 13: Significance of the best-fit
${f_{\mathrm{NL}}}$ values for the equilateral collider template (2.32), as a
function of the mass parameter $\mu$ and phase $\phi_{\rm
EC}=\mu\log(c_{\sigma}/c_{s})+\delta^{\mathrm{eq.col.}}$. Top: The
significance ${\sigma_{\mathrm{SNR}}}$ for each $\mu$ and $\phi$ is shown,
where darker colours correspond to higher significance. Bottom: The
significance for which $\phi_{\rm EC}$ fixed to the best-fit value at each
$\mu$. Overall, the maximum significance of 2.5 is attained at $\mu=3.94$ and
$\phi_{\rm EC}=2.60$.

The normalization of this template is tricky as the shape function at the
equilateral limit becomes $\propto\cos(\mu\log(1/4)+\phi_{\rm EC})$, which can
vanish for some parameter choices. We therefore normalize the shape as if the
oscillatory part equals $1$, which is equivalent to multiplying (2.32) by a
factor of $2^{5/2}/3$, or setting
$S^{\mathrm{eq.col.}}(k_{*},k_{*},k_{*})=\cos(\mu\log(1/4)+\phi_{\rm EC})$.

We first show the constraints on ${f_{\mathrm{NL}}}$ when $\phi$ is fixed to
the best-fit values at each $\mu$ in the left panel of Figure 12. Starting
with ${f_{\mathrm{NL}}}=-680\pm 520$ at $\mu=1$, the maximum significance is
reached at $\mu=3.94$ where ${f_{\mathrm{NL}}}=-178\pm 72$, eventually
reducing to ${f_{\mathrm{NL}}}=-29\pm 21$ around $\mu=8$.

The significance is shown as a function of both $\mu$ and $\phi_{\rm EC}$ in
Figure 13. At a given value of $\mu$, the significance is an oscillatory
function of $\phi$ which is always maximised (darkest colour) once in
$[0,\pi]$. The maximum significance of $2.5$ is obtained at $\mu=3.94$ and
$\phi_{\rm EC}=2.60$, which is a quite significant signal by itself. However,
after accounting for the vast parameter space scanned within the 2-parameter
search of equilateral collider signals, the adjusted significance reduces to
$0.90$, as shown in the right panel of Figure 12.

#### Low-speed collider

Figure 14: Correlation between the low-speed collider template (LSC, (2.33)),
multi-speed non-Gaussianity template (MS, (2.34)), and the three standard
templates. The correlation is based on the late-time inner product and
quantifies the similarity between two templates. As expected from the
construction, LSC (MS) templates at $\alpha=10^{-3}$ ($c=10^{-3}$) and
$\alpha=1$ ($c=1$) strongly correlate with the local and equilateral
templates, respectively. Intermediate values of $\alpha$ ($c$) interpolates
between the two.

Figure 15: Constraints on the low-speed collider template (2.33) for the
parameter $\alpha$ varied in $[10^{-3},1]$. Left: The best-fit
${f_{\mathrm{NL}}}$ estimate (black) and the 68% and 95% confidence levels
(green) are computed independently for each $\alpha$. The amplitude
${f_{\mathrm{NL}}}$ is consistent with $0$ for all $\alpha$, while being
constrained at the $\sigma({f_{\mathrm{NL}}})\sim 10-40$ level. Right: The
significance ${\sigma_{\mathrm{SNR}}}$ as a function of $\alpha$. The maximum
of $0.89$ is obtained at $\alpha=0.011$.

The template (2.33) has one free parameter $\alpha\equiv c_{s}m/H$. For
$\alpha=0$ and $1$, it is strongly correlated with the local and equilateral
template respectively, as shown in Figure 14. At intermediate values of
$\alpha$, the template probes a shape that is in between the two with a
moderate correlation with the orthogonal template. Most notable changes in the
template appear at small $\alpha$ values, so we present the constraints in log
scale in $\alpha$. The template at $\alpha=0$ is not shown here but is very
close to $\alpha=10^{-3}$.

The constraints are shown in Figure 15. As we can see, here we have
${f_{\mathrm{NL}}}=-0.4\pm 8.3$ at $\alpha=10^{-3}$ and
${f_{\mathrm{NL}}}=-2\pm 45$ at $\alpha=1$. The maximum significance is
obtained at $\alpha=0.011$, where ${f_{\mathrm{NL}}}=-9\pm 10$, which remains
consistent with ${f_{\mathrm{NL}}}=0$ within $1\sigma$ level.

Figure 16:  Constraints on the multi-speed non-Gaussianity template (2.34).
Top Left: the best-fit $f_{\rm NL}$ with $1\sigma$ and $2\sigma$ contours for
$c_{1}=c_{2}\in[10^{-3},1]$, while $c_{3}=1$ is fixed by the overall
normalization. Top Right: The raw significance ${\sigma_{\mathrm{SNR}}}$ is
shown as a function of the two sound speeds $c_{1}$ and $c_{2}$. Bottom Left:
For $c_{1}=c_{2}$, the maximum significance of 1.3 is reached at
$(c_{1},c_{2},c_{3})=(0.001,0.021,1)$. Bottom Right: The maximum significance
of 1.3 is compared with its theoretical distribution. One finds the
significance as large as 1.3 roughly $54\%$ of the time, so the result is
consistent with statistical fluctuation.

#### Multi-speed non-Gaussianities

The template (2.34) depends on the three sound speed parameters $c_{1}$,
$c_{2}$ and $c_{3}$, while the overall normalization fixes one of them.
Without loss of generality, here we assume that $c_{1},c_{2}\leq c_{3}=1$ and
normalize it so that $S(k_{*},k_{*},k_{*})=1$. As shown in Figure 14, similar
to the low-speed collider shape, the multi-speed template has a big overlap
with the local ansatz for $c_{1},c_{2}\rightarrow 0$, and becomes correlated
with the equilateral one for $c_{1},c_{2}\simeq 1$.

The top left panel of Figure 16 shows the constraints for
$c_{1}=c_{2}\in[10^{-3},1]$. At $c_{1}=c_{2}=10^{-3}$,
${f_{\mathrm{NL}}}=(0.9\pm 1.6)\times 10^{-2}$, while $c_{1}=c_{2}=1$ gives
${f_{\mathrm{NL}}}=16\pm 34$. Note that the error bars approach zero in the
limit of small sound speed $c_{1}=c_{2}\ll c_{3}$. This is an artefact of our
normalization convention of setting $S=1$ at the equilateral limit. In the
squeezed limit, the multi-speed template’s amplitude is strongly boosted
compared to that of the standard local template, resulting in a larger overall
bispectrum and hence tighter ${f_{\mathrm{NL}}}$ bounds.

The top right panel of Figure 16 shows the significance where $c_{1}$ and
$c_{2}$ are allowed to differ. In fact, the maximum significance is attained
at $(c_{1},c_{2},c_{3})=(0.001,0.021,1)$ which has $c_{1}\neq c_{2}$, as the
colour plot shows. At this point, the constraint is given by
${f_{\mathrm{NL}}}=-3.1\pm 2.3$. Note that this significance of $1.3$ is
reduced to $0.61$ after adjusting for the look-elsewhere effect, as shown in
the bottom right panel.

##  5 Summary and Outlook

The cosmological collider physics programe provides a unique window on the
particle content of the extremely high energy environment of the primordial
Universe. It has demonstrated rich phenomenology in PNG, some of which is
plausibly detectable in cosmological observations today. In this paper, we use
the legacy release of the Planck CMB observation to constrain the various
possible PNG shapes from cosmological colliders. Our analysis combines new
theoretical understanding from the cosmological bootstrap, with the recently
developed PNG-testing pipeline CMB-BEST. A brief summary of the outcomes is
listed below.

Based on the bootstrap analysis, we first proposed a set of analytic
bispectrum templates for various types of cosmological collider signals. Our
investigation incorporates the standard predictions with squeezed-limit
oscillations and angular-dependent profiles from massive exchange processes
with and without spins. Inspired by the full analytical solutions in the
bootstrap computation, we propose the approximated scalar seeds and applied
weight-shifting operators to derive the simplified templates for scalar
bispectra. These approximate templates successfully reflect the main features
of cosmological colliders and significantly enrich the shapes of available
PNG. In addition, we also collect the recently proposed new types of shapes
that arise from nontrivial sound speeds. These analytical templates (shown in
grey boxes in Section 2) provide a set of bases for the primordial bispectra
that can also be applied in future observational tests for PNG.

Next, we ran the data analysis of these shape templates of cosmological
colliders using the CMB-BEST pipeline. We use the temperature and polarization
data from the legacy release of Planck 2018\. For each template, we scanned
the parameter space to search for best-fit values of $f_{\rm NL}$ and the
corresponding 1$\sigma$ and $2\sigma$ constraints. We also computed the
signal-to-noise ratios to identify the shapes with maximum significance. The
most favoured shape from the Planck data is given by the massive scalar
exchange II template (2.20) with the mass parameter $\mu=2.13$, for which we
find $f_{\rm NL}=-91\pm 40$ at the $68\%$ confidence level. After accounting
the look-elsewhere effect, the largest adjusted significance is $1.8\sigma$.
For all the shapes, $f_{\rm NL}=0$ is still consistent with the data. The main
constraints are summarized in Table 1. We expect the analysis here to set the
stage for cosmological collider probes in future surveys, which are expected
to substantially improve sensitivity and thus have the potential for exciting
new discoveries.

The hunt for PNG is a challenging task that needs joint and collaborative
efforts between theoretical studies and observational tests. In the current
work, we have endeavoured to overcome the barriers separating the two.
Although no evidence of cosmological collider signals is found in the Planck
2018 data, at this first attempt, our analysis opens up interesting
opportunities for future exploration. Here we list several examples:

  * •

In the theory part, the starting point of our bootstrap analysis is the
requirement for nearly scale-invariance and weak coupling, which at the
leading order leads to the Feynman diagrams in Figure 1 and the corresponding
bispectrum shapes of cosmological colliders. However, other possibilities of
PNG may appear if we violate these two simple assumptions. For instance, in
some models the dominant collider signals arise from loop-level processes; and
it is possible to break the scale-invariance during inflation for feature
models. There are many examples in the literature, and they may lead to new
types of signatures not captured by our templates here. Meanwhile, as we have
shown, CMB-BEST using the general Modal methodology is a convenient tool for
testing the non-standard shapes of PNG in the Planck data. It is encouraging
to perform the data analysis for other PNG predictions.

  * •

Second, it will be important to further simplify the bispectrum templates of
cosmological colliders. The current choices can represent various typical
signatures of massive particles during inflation, however, several of them
have rather complicated analytic expressions, and there are still correlations
between these. To search for the primordial signals more efficiently in
current and future observations, we are motivated to propose a simple and
complete basis of bispectrum templates for various predictions of PNG.

  * •

Third, in CMB-BEST, our current choice of basis provides an accurate
construction for the primordial templates of cosmological colliders. An
alternative approach is to decompose the CMB bispectrum for the numerical
computations. It would be interesting to continute to cross-validate the
CMB=BEST results of this analysis with the Modal estimator used in the main
Planck analysis. Meanwhile, upcoming CMB experiments such as the Simons
Observatory and CMB-S4 will provide pristine measurements of the CMB
polarization. Forecasts on the improvements in the cosmological collider
constraints from these experiments would be illuminating, especially given the
notable (but not yet significant) signals found in our analysis. We leave
these for future investigation.

#### Note added:

During the preparation of this paper, a relevant work [98] appeared on arXiv,
which performed the observational test of cosmological colliders using the
BOSS data of galaxy surveys. While we look at the same topic, there are two
main differences in the analysis. First, Ref. [98] translated the Planck
constraints on the equilateral and orthogonal templates into bounds on
collider signals, and mainly focused on the LSS analysis; while we perform the
full CMB analysis of the PNG templates of collider signals directly using the
Planck data. Second, Ref. [98] studied the PNG signals from massive scalar
exchanges. In addition to that, we scanned the parameter space and also
incorporated collider signals from spinning exchanges. For the scalar
exchanges, we have found tighter constraints as shown in Section 4.1.

#### Acknowledgements

We would like to thank Giovanni Cabass, Xingang Chen, Guilherme Pimentel,
Petar Suman for inspiring discussions. DGW is supported by a Rubicon
Postdoctoral Fellowship awarded by the Netherlands Organisation for Scientific
Research (NWO), and partially by the STFC Consolidated Grants ST/X000664/1 and
ST/P000673/1. Part of this work was undertaken on the Cambridge CSD3 part of
the STFC DiRAC HPC Facility (www.dirac.ac.uk) funded by BEIS capital funding
via STFC Capital Grants ST/P002307/1 and ST/R002452/1 and STFC Operations
Grant ST/R00689X/1.

##  Appendix A Approximated Scalar Seed

In this appendix, we show the details about the approximated scalar seed
proposed in (2.10) and justify it as a sufficiently accurate approximation of
the exact solution $\hat{\mathcal{I}}(u)$.

Let us first review the exact solution of the boundary differential equation
for the three-point scalar seed (2.9). Details of its derivation can be found
in Ref. [59]. The solution has two parts

$\mathcal{\hat{I}}(u)=\mathcal{\hat{S}}(u)+\mathcal{\hat{H}}(u)\leavevmode\nobreak\
,$ (A.1)

where the particular piece is given by

$\mathcal{\hat{S}}=\sum_{n=0}^{\infty}c_{n}u^{n+1}\leavevmode\nobreak\
\leavevmode\nobreak\ {\rm with}\leavevmode\nobreak\ \leavevmode\nobreak\
c_{n}=\sum_{m=0}^{\lfloor
n/2\rfloor}\frac{(-1)^{n}n!/(n-2m)!}{\left[\left(n+\frac{1}{2}\right)^{2}+\mu^{2}\right]\left[\left(n-\frac{3}{2}\right)^{2}+\mu^{2}\right]...\left[\left(n+\frac{1}{2}-2m\right)^{2}+\mu^{2}\right]},$
(A.2)

and the homogeneous solution is given by

$\mathcal{\hat{H}}=-\frac{1}{2}\sum_{\pm}B_{\pm}\left(\frac{u}{2}\right)^{\frac{1}{2}\pm
i\mu}{}_{2}F_{1}\left[\frac{1}{4}\pm\frac{i\mu}{2},\frac{3}{4}\pm\frac{i\mu}{2};1\pm{i\mu};u^{2}\right]\leavevmode\nobreak\
,$ (A.3)

with

$B_{\pm}=\frac{\pi^{3/2}}{\cosh\pi\mu}\left(1\mp\frac{i}{\sinh\pi\mu}\right)\frac{\Gamma(\frac{1}{2}\pm
i\mu)}{\Gamma(1\pm i\mu)}\leavevmode\nobreak\ .$ (A.4)

One main advantage of the exact solution is that we can fully understand its
singularity structure, which can help us identify where the new physics
effects would arise. For instance, there are two nontrivial limits

$\lim_{u\rightarrow-1}\mathcal{\hat{I}}=-\frac{1}{4}\log^{2}(1+u)\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \lim_{u\rightarrow
0}\mathcal{\hat{I}}=-\frac{1}{2}\sum_{\pm}B_{\pm}\left(\frac{u}{2}\right)^{\frac{1}{2}\pm
i\mu}\leavevmode\nobreak\ ,$ (A.5)

where the first one gives us the total-energy pole, and the second one
corresponds to the collider signal as squeezed-limit oscillations. However,
for a practical purpose, this exact solution is less helpful, as it is rather
complicated and both two parts contain unphysical singularities at
$u\rightarrow 1$

$\lim_{u\rightarrow
1}\mathcal{\hat{S}}=-\frac{\pi}{2\cosh\pi\mu}\log(1-u)\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\lim_{u\rightarrow
1}\mathcal{\hat{H}}=\frac{\pi}{2\cosh\pi\mu}\log(1-u)\leavevmode\nobreak\ .$
(A.6)

In the final solution these two logarithmic poles would cancel each other and
leave the scalar seed function regular in the folded limit
$k_{3}=k_{1}+k_{2}$. But in practice, it is difficult to deal with power
series with infinite sums and we would also like to get rid of complicated
hypergeometric functions in a shape template. For example, if we want to plot
the exact solution for $0<u\leq 1$, in principle one needs to include all the
terms in the power series to obtain a converging result. One way out is to use
the exact solution (A.1) in the regime $u\in(0,0.5)$, and then for
$u\in(0.5,1)$ we solve the boundary differential equation (2.9) around $u=1$
with another converging series expansion, and glue the two solutions at
$u=0.5$. These glued solutions are shown by dashed lines in Figure 17 , which
we use as exact results for comparison.

Figure 17: The approximated scalar seed $\hat{\mathcal{I}}$ and
$\hat{\mathcal{I}}/u$ with different masses. The dashed lines are the exact
solutions from gluing. We adopt different normalizations in these two plots
for demonstration.

Next, let us try to find an approximation for the exact solution of
$\hat{\mathcal{I}}$, which should be regular at $u=1$ and can match the exact
result with good accuracy. We first notice that apart from the singularity at
$u=1$, the two parts of the exact solution basically capture two behaviours of
the scalar seed: the particular solution provides a $\mu$-dependent
equilateral-like shape for the non-squeezed momentum configurations; the
homogeneous solution gives us the oscillations in the squeezed limit. This
leads our analysis to two different cases.

  * •

For $\mu>1$, the homogeneous solution has an overall Boltzmann suppression
which leaves the particular solution more dominant for non-squeezed
configurations. Based on this observation, we propose the approximated scalar
seed in (2.10), which we copy paste here

$\hat{\mathcal{I}}_{a}(u)=\frac{u}{\beta}{(1+u)^{-\frac{\beta}{\beta+2}}}-\frac{1}{2}\sum_{\pm}B_{\pm}\left(\frac{u}{2}\right)^{\frac{1}{2}\pm
i\mu}\leavevmode\nobreak\ .$ (A.7)

Roughly speaking, the first term there presents an overall shape, while on the
top of it the second term gives the small oscillations in the squeezed-limit.
In Figure 17 we compare this template with the exact result for different
masses. One can also check that for $m\gg H$, the second term in (2.10) is
exponentially suppressed by the Boltzmann factor $e^{-\pi\mu}$, and the first
term leads to $\hat{\mathcal{I}}_{a}\sim u/(1+u)$ which corresponds to the
result of integrating out the $\sigma$ field in the heavy field limit.

  * •

For $0<\mu\leq 1$, the homogeneous and particular solutions become comparable
for most kinematic configurations, and one needs precise cancellation between
these two to generate the correct result. In this case, the two terms in (A.7)
do not agree well with the exact solution. We need to add an extra piece as

$\hat{f}(u)=-\frac{u}{6\cosh(\pi\mu)}\left[(1-2\mu^{4})\log^{2}(1+u)+\mu^{2}u^{1+\frac{8\mu^{2}}{1+8\mu^{2}}}\log(1+u)\right]\leavevmode\nobreak\
.$ (A.8)

Figure 17 also shows the comparison for $\mu=0.2,1$ between
$\hat{\mathcal{I}}_{a}(u)$ and the exact solutions from gluing. The extra term
(A.8) is important for these two results to match in the small mass regime
$m\simeq 3H/2$, and leads to nontrivial contributions in the final bispectrum.
Here we also present the extra terms in the scalar exchange templates (2.2)
and (2.17):

$\displaystyle\tilde{S}_{I}$ $\displaystyle\equiv$ $\displaystyle-
k_{1}k_{2}\partial_{k_{1}}\partial_{k_{2}}\tilde{f}$ $\displaystyle=$
$\displaystyle\frac{k_{1}k_{2}k_{3}}{6\cosh(\pi\mu)(k_{1}+k_{2})^{3}}\left[-{2(2\mu^{4}-1)}\left(\frac{k_{3}^{2}}{k_{T}^{2}}+\frac{k_{3}(4k_{T}-k_{3})}{k_{T}^{2}}\log\left(\frac{k_{T}}{k_{1}+k_{2}}\right)+\log^{2}\left(\frac{k_{T}}{k_{1}+k_{2}}\right)\right)\right.$
$\displaystyle\left.+\mu^{2}\left(\frac{k_{3}}{k_{1}+k_{2}}\right)^{\frac{1+16\mu^{2}}{1+8\mu^{2}}}\left(\frac{k_{3}(6k_{T}-k_{3}+8\mu^{2}(8k_{T}-k_{3}))}{(1+8\mu^{2})k_{T}^{2}}+2\frac{3+68\mu^{2}+384\mu^{4}}{(1+8\mu^{2})^{2}}\log\left(\frac{k_{T}}{k_{1}+k_{2}}\right)\right)\right]$
$\displaystyle\tilde{S}_{II}$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{2k_{1}k_{2}}(k_{3}^{2}-k_{1}^{2}-k_{2}^{2}){\left(1-k_{1}\partial_{k_{1}}\right)\left(1-k_{2}\partial_{k_{2}}\right)}\tilde{f}$
$\displaystyle=$
$\displaystyle\frac{k_{3}\left(k_{1}^{2}+k_{2}^{2}-k_{3}^{2}\right)}{12\cosh(\pi\mu)k_{1}k_{2}(k_{1}+k_{2})^{4}}\left[2\left(2\mu^{4}-1\right)(k_{1}+k_{2})\left((k_{1}^{2}+k_{2}^{2}+3k_{1}k_{2})\log^{2}\left(\frac{k_{T}}{k_{1}+k_{2}}\right)\right.\right.$
$\displaystyle\left.+\frac{k_{3}}{k_{T}^{2}}\left(k_{1}^{3}+k_{2}^{3}+(k_{1}^{2}+k_{2}^{2})k_{3}+7k_{1}k_{2}(k_{1}+k_{2})+5k_{1}k_{2}k_{3}\right)\log\left(\frac{k_{T}}{k_{1}+k_{2}}\right)+\frac{k_{1}k_{2}k_{3}^{2}}{k_{T}^{2}}\right)$
$\displaystyle+\mu^{2}k_{3}\left(\frac{k_{3}}{k_{1}+k_{2}}\right)^{\frac{8\mu^{2}}{8\mu^{2}+1}}\left(\frac{32\mu^{2}+3}{8\mu^{2}+1}\left(\frac{k_{1}k_{2}}{8\mu^{2}+1}-k_{1}^{2}-k_{2}^{2}-5k_{1}k_{2}\right)\log\left(\frac{k_{T}}{k_{1}+k_{2}}\right)\right.$
$\displaystyle\left.\left.+\frac{2}{8\mu^{2}+1}\frac{k_{1}k_{2}k_{3}}{k_{T}}-\frac{k_{3}}{k_{T}^{2}}\Big{(}k_{1}^{3}+k_{2}^{3}+(k_{1}^{2}+k_{2}^{2})k_{3}+11k_{1}k_{2}(k_{1}+k_{2})+9k_{1}k_{2}k_{3}\Big{)}\right)\right]\leavevmode\nobreak\
.$

#### Extension to higher order

A generalized version of the approximated scalar seed is needed for deriving
the massive spinning exchange templates (see Section 2.3 and Equation (2.25)
there). These seed functions $\hat{\mathcal{I}}^{(n)}$ satisfy boundary
differential equations similar to (2.9) but with higher-order source terms
[59]

$\left(\Delta_{u}+\mu^{2}+\frac{1}{4}\right)\hat{\mathcal{I}}^{(n)}(u)=\left(\frac{u}{1+u}\right)^{n+1}\leavevmode\nobreak\
,$ (A.11)

with $n\leq 0$ being integers. For $n=0$, we return to the scalar seed
function in (2.8) that is used for the massive scalar exchange bispectrum. For
$n>0$, this index is associated to the spin of the mediated particle as shown
in (2.25). The exact solutions of $\hat{\mathcal{I}}^{(n)}$ and details about
their analytic structures can be found in Ref. [59]. Here we propose the
following simplified approximation

$\hat{\mathcal{I}}^{(n)}_{a}=\frac{u^{n+1}}{\beta}{(1+u)^{-\frac{(n+1)\beta}{\beta+2(n+1)}}}-\frac{1}{2}\sum_{\pm}B^{(n)}_{\pm}\left(\frac{u}{2}\right)^{\frac{1}{2}\pm
i\mu}\leavevmode\nobreak\ ,$ (A.12)

with

$\beta^{(n)}=\mu^{2}+\left(\frac{1}{2}+n\right)^{2}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
B^{(n)}_{\pm}=\frac{\pi^{1/2}}{2^{n}}\frac{\Gamma\left(\frac{1}{2}+n+i\mu\right)\Gamma\left(\frac{1}{2}+n-i\mu\right)}{\Gamma\left(n+1\right)}\left(1\mp\frac{i}{\sinh\pi\mu}\right)\frac{\Gamma(\frac{1}{2}\pm
i\mu)}{\Gamma(1\pm i\mu)}\leavevmode\nobreak\ .$ (A.13)

Like (A.7), these approximated scalar seeds provide accurate templates for the
exact solutions with $\mu>1$. In the main text, we substitute (A.12) into
(2.25) to derive the massive spin-2 exchange template in (2.27).

## References

  * [1] P. D. Meerburg et al., “Primordial Non-Gaussianity,” arXiv:1903.04409 [astro-ph.CO].
  * [2] A. Achúcarro et al., “Inflation: Theory and Observations,” arXiv:2203.08128 [astro-ph.CO].
  * [3] X. Chen and Y. Wang, “Quasi-Single Field Inflation and Non-Gaussianities,” JCAP 04 (2010) 027, arXiv:0911.3380 [hep-th].
  * [4] D. Baumann and D. Green, “Signatures of Supersymmetry from the Early Universe,” Phys. Rev. D 85 (2012) 103520, arXiv:1109.0292 [hep-th].
  * [5] T. Noumi, M. Yamaguchi, and D. Yokoyama, “Effective field theory approach to quasi-single field inflation and effects of heavy fields,” JHEP 06 (2013) 051, arXiv:1211.1624 [hep-th].
  * [6] N. Arkani-Hamed and J. Maldacena, “Cosmological Collider Physics,” arXiv:1503.08043 [hep-th].
  * [7] X. Chen and Y. Wang, “Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation,” Phys. Rev. D 81 (2010) 063511, arXiv:0909.0496 [astro-ph.CO].
  * [8] V. Assassi, D. Baumann, and D. Green, “On Soft Limits of Inflationary Correlation Functions,” JCAP 11 (2012) 047, arXiv:1204.4207 [hep-th].
  * [9] X. Chen and Y. Wang, “Quasi-Single Field Inflation with Large Mass,” JCAP 09 (2012) 021, arXiv:1205.0160 [hep-th].
  * [10] S. Pi and M. Sasaki, “Curvature Perturbation Spectrum in Two-field Inflation with a Turning Trajectory,” JCAP 10 (2012) 051, arXiv:1205.0161 [hep-th].
  * [11] X. Chen, M. H. Namjoo, and Y. Wang, “Quantum Primordial Standard Clocks,” JCAP 02 (2016) 013, arXiv:1509.03930 [astro-ph.CO].
  * [12] H. Lee, D. Baumann, and G. L. Pimentel, “Non-Gaussianity as a Particle Detector,” JHEP 12 (2016) 040, arXiv:1607.03735 [hep-th].
  * [13] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Standard Model Background of the Cosmological Collider,” Phys. Rev. Lett. 118 (2017) no. 26, 261302, arXiv:1610.06597 [hep-th].
  * [14] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Standard Model Mass Spectrum in Inflationary Universe,” JHEP 04 (2017) 058, arXiv:1612.08122 [hep-th].
  * [15] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Schwinger-Keldysh Diagrammatics for Primordial Perturbations,” JCAP 12 (2017) 006, arXiv:1703.10166 [hep-th].
  * [16] A. Kehagias and A. Riotto, “On the Inflationary Perturbations of Massive Higher-Spin Fields,” JCAP 07 (2017) 046, arXiv:1705.05834 [hep-th].
  * [17] S. Kumar and R. Sundrum, “Heavy-Lifting of Gauge Theories By Cosmic Inflation,” JHEP 05 (2018) 011, arXiv:1711.03988 [hep-ph].
  * [18] H. An, M. McAneny, A. K. Ridgway, and M. B. Wise, “Quasi Single Field Inflation in the non-perturbative regime,” JHEP 06 (2018) 105, arXiv:1706.09971 [hep-ph].
  * [19] H. An, M. McAneny, A. K. Ridgway, and M. B. Wise, “Non-Gaussian Enhancements of Galactic Halo Correlations in Quasi-Single Field Inflation,” Phys. Rev. D 97 (2018) no. 12, 123528, arXiv:1711.02667 [hep-ph].
  * [20] D. Baumann, G. Goon, H. Lee, and G. L. Pimentel, “Partially Massless Fields During Inflation,” JHEP 04 (2018) 140, arXiv:1712.06624 [hep-th].
  * [21] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Neutrino Signatures in Primordial Non-Gaussianities,” JHEP 09 (2018) 022, arXiv:1805.02656 [hep-ph].
  * [22] S. Kumar and R. Sundrum, “Seeing Higher-Dimensional Grand Unification In Primordial Non-Gaussianities,” JHEP 04 (2019) 120, arXiv:1811.11200 [hep-ph].
  * [23] L. Bordin, P. Creminelli, A. Khmelnitsky, and L. Senatore, “Light Particles with Spin in Inflation,” JCAP 10 (2018) 013, arXiv:1806.10587 [hep-th].
  * [24] S. Kim, T. Noumi, K. Takeuchi, and S. Zhou, “Heavy Spinning Particles from Signs of Primordial Non-Gaussianities: Beyond the Positivity Bounds,” JHEP 12 (2019) 107, arXiv:1906.11840 [hep-th].
  * [25] S. Alexander, S. J. Gates, L. Jenks, K. Koutrolikos, and E. McDonough, “Higher Spin Supersymmetry at the Cosmological Collider: Sculpting SUSY Rilles in the CMB,” JHEP 10 (2019) 156, arXiv:1907.05829 [hep-th].
  * [26] L.-T. Wang and Z.-Z. Xianyu, “In Search of Large Signals at the Cosmological Collider,” JHEP 02 (2020) 044, arXiv:1910.12876 [hep-ph].
  * [27] D.-G. Wang, “On the inflationary massive field with a curved field manifold,” JCAP 01 (2020) 046, arXiv:1911.04459 [astro-ph.CO].
  * [28] S. Aoki and M. Yamaguchi, “Disentangling mass spectra of multiple fields in cosmological collider,” JHEP 04 (2021) 127, arXiv:2012.13667 [hep-th].
  * [29] Q. Lu, M. Reece, and Z.-Z. Xianyu, “Missing scalars at the cosmological collider,” JHEP 12 (2021) 098, arXiv:2108.11385 [hep-ph].
  * [30] L.-T. Wang, Z.-Z. Xianyu, and Y.-M. Zhong, “Precision calculation of inflation correlators at one loop,” JHEP 02 (2022) 085, arXiv:2109.14635 [hep-ph].
  * [31] X. Tong, Y. Wang, and Y. Zhu, “Cutting rule for cosmological collider signals: a bulk evolution perspective,” JHEP 03 (2022) 181, arXiv:2112.03448 [hep-th].
  * [32] Y. Cui and Z.-Z. Xianyu, “Probing Leptogenesis with the Cosmological Collider,” arXiv:2112.10793 [hep-ph].
  * [33] X. Tong and Z.-Z. Xianyu, “Large Spin-2 Signals at the Cosmological Collider,” arXiv:2203.06349 [hep-ph].
  * [34] M. Reece, L.-T. Wang, and Z.-Z. Xianyu, “Large-Field Inflation and the Cosmological Collider,” arXiv:2204.11869 [hep-ph].
  * [35] X. Chen, R. Ebadi, and S. Kumar, “Classical cosmological collider physics and primordial features,” JCAP 08 (2022) 083, arXiv:2205.01107 [hep-ph].
  * [36] Z. Qin and Z.-Z. Xianyu, “Phase information in cosmological collider signals,” JHEP 10 (2022) 192, arXiv:2205.01692 [hep-th].
  * [37] X. Niu, M. H. Rahat, K. Srinivasan, and W. Xue, “Gravitational Wave Probes of Massive Gauge Bosons at the Cosmological Collider,” arXiv:2211.14331 [hep-ph].
  * [38] X. Chen, J. Fan, and L. Li, “New inflationary probes of axion dark matter,” JHEP 12 (2023) 197, arXiv:2303.03406 [hep-ph].
  * [39] Z.-Z. Xianyu and J. Zang, “Inflation correlators with multiple massive exchanges,” JHEP 03 (2024) 070, arXiv:2309.10849 [hep-th].
  * [40] P. Chakraborty and J. Stout, “Light scalars at the cosmological collider,” JHEP 02 (2024) 021, arXiv:2310.01494 [hep-th].
  * [41] S. Jazayeri, S. Renaux-Petel, and D. Werth, “Shapes of the cosmological low-speed collider,” JCAP 12 (2023) 035, arXiv:2307.01751 [hep-th].
  * [42] S. Aoki, T. Noumi, F. Sano, and M. Yamaguchi, “Analytic formulae for inflationary correlators with dynamical mass,” JHEP 24 (2020) 073, arXiv:2312.09642 [hep-th].
  * [43] C. McCulloch, E. Pajer, and X. Tong, “A Cosmological Tachyon Collider: Enhancing the Long-Short Scale Coupling,” arXiv:2401.11009 [hep-th].
  * [44] Y.-P. Wu, “The cosmological collider in $R^{2}$ inflation,” arXiv:2404.05031 [astro-ph.CO].
  * [45] N. Arkani-Hamed, D. Baumann, H. Lee, and G. L. Pimentel, “The Cosmological Bootstrap: Inflationary Correlators from Symmetries and Singularities,” JHEP 04 (2020) 105, arXiv:1811.00024 [hep-th].
  * [46] D. Baumann, C. Duaso Pueyo, A. Joyce, H. Lee, and G. L. Pimentel, “The cosmological bootstrap: weight-shifting operators and scalar seeds,” JHEP 12 (2020) 204, arXiv:1910.14051 [hep-th].
  * [47] D. Baumann, C. Duaso Pueyo, A. Joyce, H. Lee, and G. L. Pimentel, “The Cosmological Bootstrap: Spinning Correlators from Symmetries and Factorization,” SciPost Phys. 11 (2021) 071, arXiv:2005.04234 [hep-th].
  * [48] N. Arkani-Hamed, P. Benincasa, and A. Postnikov, “Cosmological Polytopes and the Wavefunction of the Universe,” arXiv:1709.02813 [hep-th].
  * [49] P. Benincasa, “From the flat-space S-matrix to the Wavefunction of the Universe,” arXiv:1811.02515 [hep-th].
  * [50] C. Sleight, “A Mellin Space Approach to Cosmological Correlators,” JHEP 01 (2020) 090, arXiv:1906.12302 [hep-th].
  * [51] C. Sleight and M. Taronna, “Bootstrapping Inflationary Correlators in Mellin Space,” JHEP 02 (2020) 098, arXiv:1907.01143 [hep-th].
  * [52] H. Goodhew, S. Jazayeri, and E. Pajer, “The Cosmological Optical Theorem,” JCAP 04 (2021) 021, arXiv:2009.02898 [hep-th].
  * [53] S. Céspedes, A.-C. Davis, and S. Melville, “On the time evolution of cosmological correlators,” JHEP 02 (2021) 012, arXiv:2009.07874 [hep-th].
  * [54] E. Pajer, “Building a Boostless Bootstrap for the Bispectrum,” JCAP 01 (2021) 023, arXiv:2010.12818 [hep-th].
  * [55] S. Jazayeri, E. Pajer, and D. Stefanyszyn, “From locality and unitarity to cosmological correlators,” JHEP 10 (2021) 065, arXiv:2103.08649 [hep-th].
  * [56] J. Bonifacio, E. Pajer, and D.-G. Wang, “From amplitudes to contact cosmological correlators,” JHEP 10 (2021) 001, arXiv:2106.15468 [hep-th].
  * [57] S. Melville and E. Pajer, “Cosmological Cutting Rules,” JHEP 05 (2021) 249, arXiv:2103.09832 [hep-th].
  * [58] H. Goodhew, S. Jazayeri, M. H. Gordon Lee, and E. Pajer, “Cutting cosmological correlators,” JCAP 08 (2021) 003, arXiv:2104.06587 [hep-th].
  * [59] G. L. Pimentel and D.-G. Wang, “Boostless cosmological collider bootstrap,” JHEP 10 (2022) 177, arXiv:2205.00013 [hep-th].
  * [60] S. Jazayeri and S. Renaux-Petel, “Cosmological Bootstrap in Slow Motion,” arXiv:2205.10340 [hep-th].
  * [61] D. Baumann, D. Green, A. Joyce, E. Pajer, G. L. Pimentel, C. Sleight, and M. Taronna, “Snowmass White Paper: The Cosmological Bootstrap,” in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.08121 [hep-th].
  * [62] Z. Qin and Z.-Z. Xianyu, “Helical Inflation Correlators: Partial Mellin-Barnes and Bootstrap Equations,” arXiv:2208.13790 [hep-th].
  * [63] S. A. Salcedo, M. H. G. Lee, S. Melville, and E. Pajer, “The Analytic Wavefunction,” JHEP 06 (2023) 020, arXiv:2212.08009 [hep-th].
  * [64] D.-G. Wang, G. L. Pimentel, and A. Achúcarro, “Bootstrapping multi-field inflation: non-Gaussianities from light scalars revisited,” JCAP 05 (2023) 043, arXiv:2212.14035 [astro-ph.CO].
  * [65] Z. Qin and Z.-Z. Xianyu, “Closed-form formulae for inflation correlators,” JHEP 07 (2023) 001, arXiv:2301.07047 [hep-th].
  * [66] D. Stefanyszyn, X. Tong, and Y. Zhu, “Cosmological Correlators Through the Looking Glass: Reality, Parity, and Factorisation,” arXiv:2309.07769 [hep-th].
  * [67] C. Duaso Pueyo and E. Pajer, “A Cosmological Bootstrap for Resonant Non-Gaussianity,” arXiv:2311.01395 [hep-th].
  * [68] S. Céspedes, A.-C. Davis, and D.-G. Wang, “On the IR Divergences in de Sitter Space: loops, resummation and the semi-classical wavefunction,” arXiv:2311.17990 [hep-th].
  * [69] N. Arkani-Hamed, D. Baumann, A. Hillman, A. Joyce, H. Lee, and G. L. Pimentel, “Kinematic Flow and the Emergence of Time,” arXiv:2312.05300 [hep-th].
  * [70] N. Arkani-Hamed, D. Baumann, A. Hillman, A. Joyce, H. Lee, and G. L. Pimentel, “Differential Equations for Cosmological Correlators,” arXiv:2312.05303 [hep-th].
  * [71] A. Bzowski, P. McFadden, and K. Skenderis, “Renormalisation of IR divergences and holography in de Sitter,” arXiv:2312.17316 [hep-th].
  * [72] J. Chen and B. Feng, “Towards Systematic Evaluation of de Sitter Correlators via Generalized Integration-By-Parts Relations,” arXiv:2401.00129 [hep-th].
  * [73] S. Melville and G. L. Pimentel, “A de Sitter S-matrix from amputated cosmological correlators,” arXiv:2404.05712 [hep-th].
  * [74] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. IX. Constraints on primordial non-Gaussianity,” arXiv:1905.05697 [astro-ph.CO].
  * [75] Planck Collaboration, P. A. R. Ade et al., “Planck 2013 Results. XXIV. Constraints on primordial non-Gaussianity,” Astron. Astrophys. 571 (2014) A24, arXiv:1303.5084 [astro-ph.CO].
  * [76] Planck Collaboration, P. A. R. Ade et al., “Planck 2015 results. XVII. Constraints on primordial non-Gaussianity,” Astron. Astrophys. 594 (2016) A17, arXiv:1502.01592 [astro-ph.CO].
  * [77] J. R. Fergusson, M. Liguori, and E. P. S. Shellard, “General CMB and Primordial Bispectrum Estimation I: Mode Expansion, Map-Making and Measures of f_NL,” Phys. Rev. D 82 (2010) 023502, arXiv:0912.5516 [astro-ph.CO].
  * [78] J. R. Fergusson, M. Liguori, and E. P. S. Shellard, “The CMB Bispectrum,” JCAP 12 (2012) 032, arXiv:1006.1642 [astro-ph.CO].
  * [79] M. Bucher, B. Van Tent, and C. S. Carvalho, “Detecting bispectral acoustic oscillations from inflation using a new flexible estimator: Detecting bispectral acoustic oscillations from inflation using a new flexible estimator,” MNRAS 407 (2010) no. 4, 2193–2206.
  * [80] E. Sefusatti, J. R. Fergusson, X. Chen, and E. P. S. Shellard, “Effects and Detectability of Quasi-Single Field Inflation in the Large-Scale Structure and Cosmic Microwave Background,” JCAP 08 (2012) 033, arXiv:1204.6318 [astro-ph.CO].
  * [81] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, “Constraints on Single-Field Inflation from the BOSS Galaxy Survey,” Phys. Rev. Lett. 129 (2022) no. 2, 021301, arXiv:2201.07238 [astro-ph.CO].
  * [82] G. D’Amico, M. Lewandowski, L. Senatore, and P. Zhang, “Limits on primordial non-Gaussianities from BOSS galaxy-clustering data,” arXiv:2201.11518 [astro-ph.CO].
  * [83] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, “Constraints on multifield inflation from the BOSS galaxy survey,” Phys. Rev. D 106 (2022) no. 4, 043506, arXiv:2204.01781 [astro-ph.CO].
  * [84] P. D. Meerburg, M. Münchmeyer, J. B. Muñoz, and X. Chen, “Prospects for Cosmological Collider Physics,” JCAP 03 (2017) 050, arXiv:1610.06559 [astro-ph.CO].
  * [85] A. Moradinezhad Dizgah and C. Dvorkin, “Scale-Dependent Galaxy Bias from Massive Particles with Spin during Inflation,” JCAP 01 (2018) 010, arXiv:1708.06473 [astro-ph.CO].
  * [86] A. Moradinezhad Dizgah, H. Lee, J. B. Muñoz, and C. Dvorkin, “Galaxy Bispectrum from Massive Spinning Particles,” JCAP 05 (2018) 013, arXiv:1801.07265 [astro-ph.CO].
  * [87] E. Komatsu, D. N. Spergel, and B. D. Wandelt, “Measuring primordial non-Gaussianity in the cosmic microwave background,” Astrophys. J. 634 (2005) 14–19, arXiv:astro-ph/0305189.
  * [88] P. Creminelli, A. Nicolis, L. Senatore, M. Tegmark, and M. Zaldarriaga, “Limits on non-Gaussianities from WMAP data,” Journal of Cosmology and Astroparticle Physics 2006 (2006) no. 05, 004.
  * [89] A. P. S. Yadav, E. Komatsu, and B. D. Wandelt, “Fast estimator of primordial non-gaussianity from temperature and polarization anisotropies in the cosmic microwave background,” The Astrophysical Journal 664 (2007) no. 2, 680.
  * [90] L. Senatore, K. M. Smith, and M. Zaldarriaga, “Non-Gaussianities in Single Field Inflation and their Optimal Limits from the WMAP 5-year Data,” JCAP 01 (2010) 028, arXiv:0905.3746 [astro-ph.CO].
  * [91] K. M. Smith and M. Zaldarriaga, “Algorithms for bispectra: Forecasting, optimal analysis, and simulation,” Mon. Not. Roy. Astron. Soc. 417 (2011) 2–19, arXiv:astro-ph/0612571.
  * [92] J. R. Fergusson, “Efficient optimal non-Gaussian CMB estimators with polarisation,” Phys. Rev. D 90 (2014) no. 4, 043533, arXiv:1403.7949 [astro-ph.CO].
  * [93] M. Bucher, B. Racine, and B. van Tent, “The binned bispectrum estimator: template-based and non-parametric CMB non-Gaussianity searches,” JCAP 05 (2016) 055, arXiv:1509.08107 [astro-ph.CO].
  * [94] W. Sohn, J. R. Fergusson, and E. P. S. Shellard, “High-resolution CMB bispectrum estimator with flexible modal bases,” Phys. Rev. D 108 (2023) no. 6, 063504, arXiv:2305.14646 [astro-ph.CO].
  * [95] M. Shiraishi, E. Komatsu, M. Peloso, and N. Barnaby, “Signatures of anisotropic sources in the squeezed-limit bispectrum of the cosmic microwave background,” JCAP 05 (2013) 002, arXiv:1302.3056 [astro-ph.CO].
  * [96] S. Renaux-Petel, S. Mizuno, and K. Koyama, “Primordial fluctuations and non-Gaussianities from multifield DBI Galileon inflation,” JCAP 11 (2011) 042, arXiv:1108.0305 [astro-ph.CO].
  * [97] J. R. Fergusson, H. F. Gruetjen, E. P. S. Shellard, and M. Liguori, “Combining power spectrum and bispectrum measurements to detect oscillatory features,” Phys. Rev. D 91 (2015) no. 2, 023502.
  * [98] G. Cabass, O. H. E. Philcox, M. M. Ivanov, K. Akitsu, S.-F. Chen, M. Simonović, and M. Zaldarriaga, “BOSS Constraints on Massive Particles during Inflation: The Cosmological Collider in Action,” arXiv:2404.01894 [astro-ph.CO].