# Fundamental limits on constraining primordial non-Gaussianity

Alba Kalaja  P. Daniel Meerburg  Guilherme L. Pimentel  William R. Coulton

###### Abstract

We study the cosmic variance limit on constraining primordial non-Gaussianity
for a variety of theory-motivated shapes. We consider general arguments for 2D
and 3D surveys, with a particular emphasis on the CMB. A scale-invariant
$N$-point correlator can be measured with a signal-to-noise that naively
scales with the square root of the number of observed modes. This intuition
generally fails for two reasons. First, the signal-to-noise scaling is reduced
due to the blurring of the last scattering surface at short distances. This
blurring is caused by the combination of projection and damping, but the loss
of signal is not due to exponential decay, as both signal and noise are
equally damped. Second, the behavior of the $N$-point correlator in the
squeezed and collapsed (for $N>3$) limits can enhance the scaling of the
signal-to-noise with the resolution, even with a reduced range of momenta
probing these limits. We provide analytic estimates for all $N$-point
correlators. We show that blurring affects equilateral-like shapes much more
than squeezed ones. We discuss under what conditions the optimistic scalings
in the collapsed limit can be exploited. Lastly, we confirm our analytical
estimates with numerical calculations of the signal-to-noise for local,
orthogonal and equilateral bispectra, and local trispectra. We also show that
adding polarization to intensity data enhances the scaling for equilateral-
like spectra.

##  1 Introduction

The statistics of primordial fluctuations imprinted in the Cosmic Microwave
Background (CMB) and visible in the late-time clustering of galaxies is our
most important window into the very early universe. A measurement and
characterization of the primordial fluctuations is necessary to improve our
understanding of a putative period of inflation, shedding light on its
microscopic origin. As future experiments achieve higher sensitivity and
resolution, increasing the capacity to detect statistics beyond the two-point
correlation function, i.e. non-Gaussianity [1], both CMB [2, 3, 4, 5] and
large scale structure (LSS) missions [6, 7, 8, 9, 10, 11] will face challenges
associated to removing foregrounds and taking into account non-primordial
effects (see e.g. [12, 13] and [14]). Besides these systematic challenges, we
run into limitations determined by the nature of the tracer field. For
example, the CMB is fundamentally limited by two factors:

  * •

the last scattering surface is two-dimensional, thus observed anisotropies are
a 2D projection of the 3D fluctuations;

  * •

the thickness of the last scattering surface requires that we average over all
projections along the line of sight, effectively blurring the non-Gaussianity
on small scales, as shown in Fig. 1.

In this paper we build a qualitative and quantitative understanding of how
these limitations affect the signal-to-noise. We tackle the problem in two
ways: first, using a scale invariant spectrum of non-Gaussianity as
theoretical prior, we derive the scaling relation of the signal-to-noise
ratios ($S/N$) with the number of measured fluctuations (modes), i.e. the
inverse of the resolution of a given experiment. Secondly, we perform a
numerical Fisher analysis to test the consistency of our theoretical
predictions. In both cases, we will assume cosmic variance limited
observations and neglect non-primordial sources that could hinder a detection
of primordial non-Gaussianity.

In order to contextualize our motivation, let us consider the bispectrum of
temperature anisotropies. First, neglecting the effect of diffusion damping on
small scales, a simple estimate shows that the $S/N$ scales linearly with the
maximal measured multipole $\ell_{\mathrm{max}}$. This follows from scale
invariance of the bispectrum. However, when the damping effect is taken into
account in the transfer function, the $S/N$ shows a much poorer scaling and
grows with the square root of $\ell_{\mathrm{max}}$. This result was derived
in Ref. [15] in the context of equilateral non-Gaussianity. Here, we show
that, at least in these simple estimates, the problem gets worse at higher
point correlation functions, to the extreme that the $S/N$ converges (for
five- and higher point correlation functions), even if all the CMB modes are
measured. We argue that this result is a general consequence of the blurring
of non-Gaussianity at small scales.

Figure 1: The blurriness of the last scattering surface at short scales washes
out our ability to detect primordial non-Gaussianity. In the figure, we show a
cartoon that largely exaggerates this effect for illustration purposes. When
we measure the correlations within a triangle in the sky, we must average over
all possible triangles along the same line of sight. This averaging reduces
the signal to noise ratio as a function of the survey resolution $\ell_{\rm
max}$.

Furthermore, in view of recent advances in the theory of cosmological
correlators, for our theoretical estimates we go beyond equilateral and local
bispectra analysis, and consider more general shapes. In particular, around
squeezed or collapsed limits—when there is a hierarchy between distance scales
being correlated—, non-Gaussianity exhibits features similar to those of a
particle detector in a collider experiment, potentially probing new particles
that decay into the primordial fluctuations [16, 17]. While the exact scaling
is determined by the full shape of a given primordial $N$-point correlation
function—or, equivalently, the bispectrum, trispectrum, etc.111A slightly
confusing terminology is that $N$-point correlation functions are referred to
as $(N-1)$-spectra. Thus the $3$-point correlation function is the bispectrum,
the $4$-point correlation function is the trispectrum, etc.— we find that the
leading scaling can be estimated by computing only the squeezed and collapsed
limits of correlators. This greatly simplifies calculations of $S/N$ ratios,
allowing us to adapt the effect of the blurriness and to extend the argument
to all ($N-1$)-spectra. Despite the reduced phase space in restricted
kinematics, for certain shapes, the presence of damping at small scales turns
out to have little influence on the $S/N$ scaling. Intuitively, as shown in
Fig. 2, when large scale perturbations are correlated with small ones, the
signal is well defined since the possible triangles keep the same shape along
the line of sight. On the other hand, a signal coming from the correlation of
perturbations with a similar scale (which is below the damping scale) receives
contributions from different shapes due to the thickness of the last
scattering surface. Thus the average over all possible triangles along the
line of sight produces a blurriness that reduces the $S/N$.

Figure 2: For shapes of non-Gaussianity peaked around the squeezed limit,
blurriness does not reduce the $S/N$ ratio much. This is because within the
same line of sight the possible triangles do not change shape drastically. For
equilateral-like non-Gaussianity, we have to average over triangles with
various shapes and the blurriness reduces the $S/N$ significantly.

We show that, indeed, equilateral-like shapes have a more reduced $S/N$ than
shapes with different wavelengths modes correlated. In particular, we find
that squeezed bispectra and collapsed trispectra have a large $S/N$ scaling,
in the case where it is sourced by the exchange of very light particles during
inflation. These results were already shown in Ref. [18] and Refs. [19, 20]
for the bispectrum and the trispectrum respectively. However, the enhanced
scaling of the trispectrum needs to be interpreted properly, within a specific
model generating the shapes [21]. We discuss this later in the paper.

We will confirm our heuristic estimates for the bispectrum and the trispectrum
with a full-sky Fisher estimate. Here, we will use simplified templates of
non-Gaussianity, that still capture the relevant physics of the ab initio
shapes. We also include $E$-mode polarization in our estimates, and show that
they increase the $S/N$ ratio, sometimes parametrically in
$\ell_{\mathrm{max}}$. Our interpretation of this finding is that polarization
knows about the velocity field around recombination, and effectively increases
the dimensionality of the CMB to be slightly above 2D. A similar effect is
found in Ref. [22], where Rayleigh scattering is included as a tracer of
primordial fluctuations. Adding data to the primary CMB temperature modes only
improves the $\ell_{\mathrm{max}}$ scaling of shapes for which the scaling is
not already optimal (e.g. squeezed shapes, which already reach mode-counting
scaling).

The paper is organized as follows: in Sec. 2, we briefly review the shapes of
the correlators in the squeezed and collapsed limits; in Sec. 3, we show the
theoretical estimates of the $S/N$; in Sec. 4, we outline the numerical Fisher
analysis and compare the results with the theoretical ones; in Sec. 5, we
discuss the results of the paper and future prospects. In the appendices we
include additional technical details and derivations: App. A, B and D contain
a brief review of CMB statistics and $S/N$ derivation; in App. C we explicitly
show the effect of diffusion damping, and in App. E we derive multiple
squeezed and collapsed limits.

Results The signal-to-noise ratio is a function of various quantities:
$\ell_{\rm max}$ ($\ell_{\rm min}$), the smallest (largest) angular
resolutions of the survey; $f_{\rm NL}$, the size of the non-Gaussianity222In
this case, $f_{\rm NL}$ is a general parameter, not to be confused with the
amplitude of the bispectrum.; $f_{\rm sky}$, the available fraction of the sky
being measured; and the details of the shape–local, equilateral, degree of
correlation function, etc. For weak non-Gaussianity, the dependence on $f_{\rm
NL}$ and $f_{\mathrm{sky}}$ is very simple. They appear as overall factors in
the signal-to-noise. The dependence on $\ell_{\rm max}$, $\ell_{\rm min}$ (or
$k_{\rm max}$, $k_{\rm min}$ for a $3$D survey) is more complicated. A
detailed analysis of the various shapes gives us

$\left(\frac{S}{N}\right)^{2}\sim(f_{\mathrm{sky}}f_{\rm NL})^{2}\ell_{\rm
max}^{p}~{},$ (1.1)

where we omit numerical factors and the $\ell_{\rm min}$ dependence333In
principle, $\ell_{\mathrm{min}}$ is limited by the fraction of the sky
observed. For simplicity, we assume $\ell_{\mathrm{min}}$ to have a fixed
value.. The power $p$ depends on the shape of the non-Gaussianity and whether
$\ell_{\rm max}$ is above or below the damping scale $\ell_{D}$. We
parametrize a large family of non-Gaussian shapes by their behavior around
squeezed and collapsed limits, and find the corresponding scaling $p$ above
and below the damping scale. The resulting values of $p$ are summarized in
Tab. 8. Notation and conventions The modulus of a vector is given by
$|\bm{k}|\equiv k$ and $|\bm{\ell}|\equiv\ell$. We use $\bm{k}_{ij\dots
n}=\bm{k}_{I}+\bm{k}_{j}+\dots+\bm{k}_{n}$ and $\bm{\ell}_{ij\dots
n}=\bm{\ell}_{I}+\bm{\ell}_{j}+\dots+\bm{\ell}_{n}$.

We denote the primordial curvature perturbation
$\zeta(\bm{k})=\zeta_{\bm{k}}$. The various moments of the fluctuations are
defined as

$\begin{split}&\langle\zeta_{\bm{k}_{1}}\zeta_{\bm{k}_{2}}\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}_{12})P_{\zeta}(k_{1}),\\\
&\langle\zeta_{\bm{k}_{1}}\zeta_{\bm{k}_{2}}\zeta_{\bm{k}_{3}}\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}_{123})B_{\zeta}(k_{1},k_{2},k_{3}),\\\
&\langle\zeta_{\bm{k}_{1}}\zeta_{\bm{k}_{2}}\zeta_{\bm{k}_{3}}\zeta_{\bm{k}_{4}}\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}_{1234})T_{\zeta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4}),\\\
&\langle\zeta_{\bm{k}_{1}}\zeta_{\bm{k}_{2}}\cdots\zeta_{\bm{k}_{N}}\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}_{12\dots
N})F_{\zeta}(\bm{k}_{1},\bm{k}_{2},\dots,\bm{k}_{N}).\end{split}$ (1.2)

In our numerical computations we consider a flat $\Lambda$CDM cosmology, with
cosmological parameters in accordance with the latest Planck results [23],
summarized in Tab. 1.

|  $\bm{\Lambda}$CDM parameters |
---|---|---
H${}_{0}=67.66$ | $\Omega_{\mathrm{b}}h^{2}=0.02242$ |  $\sum$m${}_{\nu}=0.06$
$\Omega_{\mathrm{k}}=0$ | $\Omega_{\mathrm{c}}h^{2}=0.11993$ | $\tau=0.0561$
$n_{s}=0.9665$ | $A_{\mathrm{s}}=2.1056\times 10^{-9}$ | $r=0$
Table 1: Best-fit Planck parameters (specifically, Tab. 2 of Ref. [23] with
$TT$, $TE$, $EE+$low$E+$lensing+BAO) used in our numerical computations.

##  2 Non-Gaussianity from inflation

Cosmic inflation provides a compelling mechanism for the generation of
primordial fluctuations—as the universe expands and flattens, quantum
mechanical vacuum fluctuations are stretched to cosmological distances. During
the past decades and more recently from the Planck mission [23], measurements
of the CMB fluctuations have shown consistency with the predictions of slow-
roll models of inflation: superhorizon primordial fluctuations are nearly
scale-invariant, adiabatic and nearly Gaussian. While the single field slow-
roll model predicts a small amount of non-Gaussianity [24, 25], coming from
gravitationally mediated self-interactions of the inflaton, other models of
inflation predict stronger signals. Different physical processes, mediated by
particles with masses of order the Hubble scale during inflation, give rise to
distinctive signatures in the non-Gaussian signal. The precision calculation
of such signatures is a rich subject, much akin to the area of scattering
amplitudes in particle physics. In certain kinematical configurations, such as
squeezed and collapsed limits, the shape of non-Gaussianity reveals
information about the mass and the spin of particles mediating the
interactions among the curvature fluctuations [26, 27, 16, 17, 1]. Generally,
the $N$-point correlation functions of the primordial curvature fluctuations
are given by

$\langle\zeta_{\bm{k}_{1}}\zeta_{\bm{k}_{2}}\cdots\zeta_{\bm{k}_{N}}\rangle=\delta^{(3)}(\bm{k}_{12\dots
N})F_{\zeta}(\bm{k}_{1},\dots,\bm{k}_{N})$ (2.1)

Assuming rotation and translation invariance, the degrees of freedom of the
scalar function $F_{\zeta}(\bm{k}_{1},\cdots,\bm{k}_{N})$ will be reduced to
$3(N-2)$. All the information about the amplitude of non-Gaussianity and its
specific shape is encoded in $F_{\zeta}$. Below, we focus on the shape of
$F_{\zeta}$ and we leave out the coefficients in front of it, which are
related to the amplitude.

Initially we ignore the angular dependence of the correlators and focus on
investigating how the volume of phase space affects the $S/N$444The angular
dependence of correlators as a function of momenta most clearly gives the
signature of a spinning particle mediating the interaction among curvature
fluctuations.. We then examine specific kinematical regimes, i.e. angular
limits, and show how these limits determine the dominant scaling of the $S/N$.

When all momenta are of similar size, then we correlate wavelengths generated
at the same time during inflation, and we expect to probe contact, self-
interactions of the inflaton. Around this kinematics, non-Gaussianities are
devoid of features, and have a scaling consistent with the symmetries of the
inflationary background555 The scaling of $F_{\zeta}$ has a factor of
$k^{3(N-1)}$ for an $N$-point function. The extra factor of $k^{-3}$ comes
from the momentum conserving delta function, rendering the correlator to be
scale invariant.,

$F_{\zeta}(k_{1}\sim k_{2}\sim\cdots\sim
k_{N})\sim\frac{1}{k_{1}^{3(N-1)}}\,.$ (2.2)

The details of the shape become clearer as we dial the various momenta,
generating a hierarchy between the times in which the various fluctuations are
sourced during inflation.

There are two interesting kinematical limits of the correlator that generate a
hierarchy between momenta:

  * •

Squeezed limits—When one of the momenta is much smaller than all the others,
we probe correlations between fluctuations generated at two different times
during inflation. If there are new particles $\sigma$, perhaps massive, that
couple to the inflaton, then their imprints become clearest in the squeezed
limit, acting as mediators between fluctuations generated at earlier and later
times (see left diagram in Fig. 3). In the squeezed limit, a single leg probes
the later time, and due to momentum conservation, the long mediator has the
same wavelength as said external fluctuation. A convenient way of
parametrizing the general behavior of the squeezed limit is666Throughout, we
will drop anything but the degrees of freedom relevant to $F_{\zeta}$.

$\lim_{k_{1}\ll
k_{2,\dots,N}}F_{\zeta}(k_{1},k_{2},\dots,k_{N})\sim\frac{1}{k_{1}^{3}k_{2}^{3(N-2)}}\left(\frac{k_{1}}{k_{2}}\right)^{\Delta}+\dots,$
(2.3)

where $\dots$ are subleading terms, suppressed by powers of the small ratio
$k_{1}/k_{2}$. The leading order exponent $\Delta$ is related to the mass of
the particle in Hubble units. In detail, when inflation is very well
approximated by de Sitter space (zero-order in slow-roll), then
$\Delta=3/2-\sqrt{9/4-m_{\sigma}^{2}/H^{2}}$. We will be agnostic about the
specific microscopic origin of $\Delta$, but it is useful to keep in mind that
it is related to the mass of a mediator particle.

Figure 3: Left: Squeezed $N-$point functions. Cartoon of physical process
sourcing squeezed non-Gaussianity. Time grows along the vertical direction and
the horizontal line marks the end of inflation. The squeezed momentum
corresponds to a fluctuation with very long wavelength, sourced at early
times. A particle $\sigma$ mediates a correlation between the long mode and
the $N-1$ short modes sourced at later times. The scaling behavior of the
correlator is related to the mass (in Hubble units) of $\sigma$. Right: A
collapsed $N-$point function configuration. We assume that all momenta on each
“bundle” have similar magnitude, but the “left” and “right” momenta can be
different. In the figure, left and right momenta are sourced at similar times,
but the correlation is mediated by $\sigma$, which is generated earlier. In
the limit where the exchanged momentum is much smaller than the momenta on
each side of the diagram, the correlation function will exhibit scaling
controlled by the mass of the exchanged particle.

  * •

Collapsed limits—Instead of squeezing a side of a polygon, we can correlate
two collections of modes generated late in inflation mediated by a pair
produced at early times, with long wavelength. In order to probe such process,
we consider “collapsing” a polygon (at least a quadrilateral, i.e. the
trispectrum) by making a diagonal very small777Note that permutation
symmetries of $N$-point functions imply correlators must be the same for two
polygons with very different shapes, obtained by permuting the order of the
sides. As a consequence of this, it is reductive to assume the collapsed
limits arise only when a diagonal becomes small. In turn, collapsed
configurations correspond also to cases where none of the geometrical
diagonals vanish and it is better to say that a subset of momenta adds up to
zero. The results we obtain in Sec. 3 therefore apply more broadly, but are
derived in the limit where the configuration is chosen such that the diagonal
vanishes.. More generally, we can bundle the momenta in two groups, as shown
in the right panel in Fig. 3. Then, denoting the “left” momenta
$\\{k_{1},\cdots k_{M}\\}\sim k_{L}$ and the “right” momenta
$\\{k_{M+1},\cdots k_{N}\\}\sim k_{R}$ (with $M>1$), the mass of the mediator
particle will be read off by the scaling $\Delta$ in the collapsed limit
$k_{I}\ll k_{L},k_{R}$

$\lim_{k_{I}\ll
k_{L},k_{R}}F_{\zeta}(k_{I},k_{L},k_{R})\propto\frac{1}{k^{3(M-1)}_{L}k^{3}_{I}k^{3(N-M-1)}_{R}}\left(\frac{k_{I}^{2}}{k_{R}k_{L}}\right)^{\Delta}+\dots,$
(2.4)

where $k_{I}\equiv k_{1}+\cdots+k_{M}$. Notice that both the squeezed limit
and the collapsed limit give access to similar processes, but the squeezed
limit requires linear mixing of $\sigma$ with the inflaton, while more general
interactions are accessible via collapsed limits. Moreover, momentum
conservation still gives certain phase space in the various configurations of
left and right momenta, while keeping the diagonal collapsed.

Finally, for large enough values of $\Delta$, the non-Gaussianity quickly
decays as any hierarchy between momenta is generated. These shapes generally
arise from contact self-interactions of the primordial fluctuations during
inflation, and are referred to as equilateral non-Gaussianities888In the
effective field theory of single-field inflation, there are two distinct self-
interactions of the inflaton that produce distinct patterns of non-Gaussianity
in the bispectrum. One of them is referred to as the equilateral shape, while
the other is referred to as the orthogonal shape. Nonetheless, both of those
shapes probe the physics of self-interactions of the curvature..

A summary of the discussion above is presented in Fig. 4. A measurement of
$\Delta$ provides spectroscopic information about what sets up the initial
conditions, correlating primordial fluctuations along different scales. In the
simplest inflationary scenarios, $0<\Delta<3/2$ signals the presence of a
light (in Hubble units) mediator particle setting long range correlations
among the inflaton particles. Values $\Delta\geq 2$ probe self-interactions of
the inflaton and include equilateral non-Gaussianity. In many models, as for
example in slow-roll, quasi-single field inflation, the value of $\Delta$ in
the squeezed bispectrum is often the same as that in the collapsed
trispectrum. A notable exception to this rule is graviton exchange. In the
collapsed limit, the masslessness of the graviton is evident, and the scaling
of the trispectrum is $\Delta=0$. In the squeezed bispectrum, we only exchange
the static Newtonian potential between the fluctuations, which is effectively
a contact interaction, and the scaling is that of an equilateral shape, with
$\Delta=2$.

Figure 4: The value of $\Delta$ encodes interesting new physics of particles
generating long range forces during inflation. The plot above shows the
relation between the mass of a new particle and the scaling of the squeezed
limit of the bispectrum. For higher masses, the values of $\Delta$ become
complex (see the pink arrows), and the non-Gaussian signal develops an
oscillating pattern. Notice the “mass gap” for $3/2<\Delta<2$ – the simplest
inflationary models predict that the leading contribution to the bispectrum
should not lie on that range.

Our discussion of shapes of non-Gaussianity has not been exhaustive, but it
will suffice for the analysis below. We will assume general values of
$\Delta$, while not necessarily the same in the squeezed vs. collapsed limit.
We should mention that other shapes leave very interesting features, related
to spins of the particles, bursts of particle production, different
inflationary mechanisms, alternative scenarios for sourcing the initial
conditions etc. We leave the extension of our analysis to other shapes for
future work.

##  3 Theoretical signal-to-noise ratio estimation

In this section, we give analytical estimations of the $S/N$ of non-
Gaussianity from general $N$-point correlation functions. Our estimates assume
that our only limitation in directly probing non-Gaussianity is the finite
number of modes in the Universe, i.e. cosmic variance.

In principle, we expect the $S/N$ of the amplitude of non-Gaussianities to
scale with the square root of the total number of modes in a survey. For a 3D
survey, as in the case of an LSS mission, the total number of modes is
proportional to $N_{3\mathrm{D}}\sim 1/r_{\rm min}^{3}$, with $r_{\rm min}$
the smallest observed scale within the survey, or equivalently in Fourier
space $N_{3\mathrm{D}}\sim k_{\mathrm{max}}^{3}$.

For the CMB, the survey is confined to a 2D surface, therefore we expect the
total number of modes to be proportional to
$N_{2\mathrm{D}}\sim\ell_{\mathrm{max}}^{2}$, with $\ell_{\mathrm{max}}$ the
largest accessible multipole. However, some effects that are intrinsic to our
probes, such as the finite thickness of the last scattering surface for the
CMB, drastically change the amount of information that we can access,
underestimating it with respect to mode counting.

In Sec. 3.1, we give a general overview of the scaling of $S/N$ in terms of
$k_{\mathrm{max}}$ ($\ell_{\mathrm{max}}$) for a 3D (2D) survey. Then, we show
how damping affects the $S/N$ scaling for a CMB survey for general
$(N-1)$-spectra. Here, it is interesting to point out an important distinction
between the bispectrum and the trispectrum. In order to specify a triangle in
2D vs. 3D, three pieces of data (the sides of the triangle) are needed. In
that sense, there is no loss of freedom in going from a tomogram to a
projection. This is not true for a quadrilateral, which would require six
pieces of data (four sides, two diagonals) in 3D, while requiring five pieces
of data (four sides, one diagonal) in 2D999Strictly speaking, given four sides
and one diagonal, there are two possible quadrilaterals in two dimensions, one
concave and one convex. However, it remains true that there is no freedom to
dial the size of the second diagonal arbitrarily. (see e.g. Fig. 5). This
restriction does not seem to play an important role in the asymptotic
estimates below, but it might be important in analyzing features or more
detailed signatures of new physics in the trispectrum, where having the
freedom to change the various shapes of the quadrilateral becomes important.

Figure 5: The last scattering surface is a 2D surface where fluctuations are
projected.

Finally, in Sec. 3.2 and 3.3 we study the scaling properties of the $S/N$ in
the squeezed and collapsed limits respectively, highlighting the dependence on
the specific shape considered.

For the sake of clarity, we use the flat-sky approximation (see e.g. App. B)
to describe the CMB statistics and we limit our analysis to temperature
perturbations. These estimations are far from being rigorous, but they convey
the message without lengthy calculations.

###  3.1 General estimation: the damping effect on the $S/N$

For a 3D survey, the $(S/N)^{2}$ of an $N$-point correlation functions is
given by

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(N)}\sim\prod_{i=1}^{N}\left(\int\mathrm{d}^{3}k_{i}\right)\,\delta^{(3)}(\bm{k}_{12\dots
N})\frac{F^{2}_{\zeta}(k_{1},k_{2},\dots,k_{N})}{P(k_{1})P(k_{2})\cdots
P(k_{N})}.\end{split}$ (3.1)

Here we neglect all coefficients in front of the integral since we only focus
on the scaling in terms of $k_{\mathrm{max}}$. In general, for a 3D survey the
scale-invariant $(N-1)$-spectrum scales as $F\sim k^{-3(N-1)}$, therefore by
replacing it in Eq. (3.1), we obtain

$\left(\frac{S}{N}\right)^{2}_{(N)}\sim\int^{k_{\mathrm{max}}}\mathrm{d}^{3(N-1)}k\,\frac{k^{-6(N-1)}}{k^{-3N}}\sim\int^{k_{\mathrm{max}}}\mathrm{d}kk^{2}\sim
k_{\mathrm{max}}^{3},$ (3.2)

where we used the Dirac-delta to integrate one of the variables out.

The same argument applies to a 2D survey, where the $(S/N)^{2}$ reads as

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(N)}\sim\prod_{i=1}^{N}\left(\int\mathrm{d}^{2}\ell_{i}\right)\,\delta^{(2)}(\bm{\ell}_{12\dots
N})\frac{F^{2}(\ell_{1},\ell_{2},\dots,\ell_{N})}{C(\ell_{1})C(\ell_{2})\cdots
C(\ell_{N})},\end{split}$ (3.3)

Here we distinguish between two regimes: one where diffusion damping is
ineffective, namely where $\ell\ll\ell_{\mathrm{D}}$ with
$\ell_{\mathrm{D}}\simeq 1300$ the damping scale, the other one where
diffusion damping becomes dominant, i.e. $\ell\gg\ell_{\mathrm{D}}$. In the
former, the flat-sky angular $(N-1)$-spectrum scales as $F\sim\ell^{-2(N-1)}$
(for more details see App. B). Consequently, the $(S/N)^{2}$ scales as

$\left(\frac{S}{N}\right)^{2}_{(N)}\sim\int^{\ell_{\mathrm{max}}}\mathrm{d}^{2(N-1)}\ell\frac{\ell^{-4(N-1)}}{\ell^{-2N}}\sim\int^{\ell_{\mathrm{max}}}\mathrm{d}\ell\,\ell\sim\ell_{\mathrm{max}}^{2}.$
(3.4)

As one would intuitively expect, both 3D and 2D estimates are proportional to
the number of modes per survey, and the results are general for any
$(N-1)$-spectra.

The situation changes drastically in the $\ell\gg\ell_{\mathrm{D}}$ regime,
namely when the effect of the radiation transfer function, which we previously
neglected, is incorporated in the estimation. As we show in detail in App. C,
the effect of the transfer function can be included by replacing
$\ell^{2}\rightarrow\ell^{2}(\ell/\ell_{D})\mathrm{e}^{-(\ell/\ell_{\mathrm{D}})^{1.2}}$.
The damping in the exponential factor in the flat-sky transfer function (C.1)
cancel out in the $S/N$ (it appears equally in the numerator and the
denominator), however the line-of-sight damping, which is sensitive to the
thickness of the last scattering surface, modifies the scaling behavior of the
$(N-1)$-spectra, and ultimately of the $S/N$. Indeed, if we add this effect in
Eq. (3.4), we obtain

$\left(\frac{S}{N}\right)^{2}_{(N)}\sim\ell_{\mathrm{D}}^{N-2}\int
d^{2(N-1)}\ell\,\frac{\ell^{-6(N-1)}}{\ell^{-3N}}\sim\ell_{\rm max}^{4-N}.$
(3.5)

The $N$ dependence implies that there is a intrinsic limitation in probing the
shape of a general scale-invariant $(N-1)$-spectrum. While for $N=4$ the
scaling is logarithmic,

$\left(\frac{S}{N}\right)^{2}_{(4)}\sim\ell_{\mathrm{D}}^{2}\int^{\ell_{\mathrm{max}}}\mathrm{d}^{6}\ell\frac{\ell^{-18}}{\ell^{-12}}\sim\log\ell_{\mathrm{max}},$
(3.6)

the $(S/N)^{2}$ actually converges for $N>4$, where we expect

$\left(\frac{S}{N}\right)^{2}_{(N>4)}\sim
A+\frac{B}{\ell_{\mathrm{max}}^{p}}\to A,~{}~{}{\rm
as}~{}\ell_{\mathrm{max}}\to\infty.$ (3.7)

As a consequence, diffusion damping along the line of sight puts a fundamental
stop gap at how accurate we can distinguish signal from noise with $N>4$
correlation functions. In this sense, looking for kinematical corners, where
the scaling can improve, becomes crucial. In the following, we study the
squeezed and collapsed limits of correlation functions. We will refer to the
tracers with a reduced $S/N$ as damped tracers while keeping in mind that the
loss of information is due to a combination of projection and damping.

###  3.2 The squeezed limit

For the estimates made above, we used only the scaling properties of the
$(N-1)$-spectrum, without making any assumptions about potential divergences
around specific configurations or specific kinematic regimes. In this
subsection, we estimate the scaling of the $(S/N)^{2}$ in the squeezed regime,
stressing its dependence on the mass of mediators $\Delta$, hence on the shape
of non-Gaussianity. We show that for certain values of $\Delta$, the scaling
estimates in Sec. 3.1 turn out to be an underestimation of the squeezed
scaling. This result is trivial in the case of 3D surveys, however for high
resolution CMB maps, well beyond the damping scale $\ell_{\mathrm{D}}$, the
squeezed limit of certain shapes produces a mode-counting scaling.

We first describe in detail how to estimate the $(S/N)^{2}$ for the 3D and 2D
$N$-point correlation function, the latter in the $\ell\ll\ell_{\mathrm{D}}$
region. We then investigate the CMB bispectrum and trispectrum $(S/N)^{2}$
scaling in the damping region, $\ell\gg\ell_{\mathrm{D}}$. These results are
then generalized to the CMB $(N-1)$-spectrum.

Given the marked hierarchy between the momenta $k_{1}\ll k$, with $k=k_{2}\sim
k_{3}\sim\dots\sim k_{N}$, we can integrate $k_{1}$ from a minimum $k$-mode,
$k_{\mathrm{min}}$, up to $k_{1}=ck_{\mathrm{max}}$, with $c\ll 1$. Then the
$(S/N)^{2}$ becomes

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(N)}&\sim\int_{k_{\mathrm{min}}}^{ck_{\mathrm{max}}}\mathrm{d}^{3}k_{1}\int_{ck_{\mathrm{max}}}^{k_{\mathrm{max}}}\mathrm{d}^{3(N-2)}k\,\frac{k_{1}^{-6+2\Delta}k^{-6(N-2)-2\Delta}}{k_{1}^{-3}k^{-3(N-1)}}\\\
&\sim\int_{k_{\mathrm{min}}}^{ck_{\mathrm{max}}}\mathrm{d}k_{1}k_{1}^{2\Delta-1}\int_{ck_{\mathrm{max}}}^{k_{\mathrm{max}}}\mathrm{d}k\,k^{2-2\Delta},\end{split}$
(3.8)

where we used the squeezed $(N-1)$-spectrum defined in Eq. (2.3). The solution
of the integral depends on the value of $\Delta$, but not on the correlation
function considered, meaning that, given that shape, the results is universal
for any bispectrum, trispectrum, etc.

It is straightforward to show that, for general $\Delta>0$, the scaling is
given by mode-counting $\sim k_{\mathrm{max}}^{3}$. On the other hand, if
$\Delta=0$, which corresponds to the local type of non-Gaussianity, the
$(S/N)^{2}$ shows a more favourable scaling with the number of modes,
proportional to $\sim
k_{\mathrm{max}}^{3}\log(k_{\mathrm{max}}/k_{\mathrm{min}})$. In the extreme
case where the $N$-point correlation function diverges faster than the local
shape in the squeezed limit, namely for configurations with
$\Delta<0$101010$\Delta<0$ corresponds to mediator particles with imaginary
mass, see e.g. [28, 29] for inflationary scenarios where tachyonic
instabilities are present for a few e-folds. Furthermore, the super-squeezed
scaling might also imply a non-local theory [30]., we expect an enhancement of
$k_{\mathrm{max}}^{3}$ by a factor of
$(k_{\mathrm{max}}/k_{\mathrm{min}})^{-2\Delta}$ .

The argument applies straight away to the CMB, with a simple
$k^{3}\rightarrow\ell^{2}$ adjustment. In the region below the damping scale,
the squeezed angular $(N-1)$-spectrum is given by

$\lim_{\ell_{1}\ll\ell\ll\ell_{\mathrm{D}}}F_{\zeta}(\ell_{1},\ell_{2},\dots,\ell_{N})\propto\frac{1}{\ell_{1}^{2}\ell^{2(N-2)}}\left(\frac{\ell_{1}}{\ell}\right)^{2\Delta/3},$
(3.9)

with $\ell=\ell_{2}\sim\ell_{3}\sim\dots\sim\ell_{N}$. Since the integration
is equivalent to the one outlined for the 3D $S/N$, we report only the
results. We find

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(N)}\sim\begin{cases}\ell_{\mathrm{max}}^{2},&\Delta>0\\\
\ell_{\mathrm{max}}^{2}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{min}}),&\Delta=0.\end{cases}\end{split}$
(3.10)

As in the 3D case, for $\Delta<0$, we expect an enhancement of
$\ell_{\mathrm{max}}^{2}$ by a factor of
$(\ell_{\mathrm{max}}/\ell_{\mathrm{min}})^{-4\Delta/3}$.

| Squeezed limit
---|---
| 3D | 2D ($\ell\ll\ell_{\mathrm{D}}$)
$\Delta<0$ | $k_{\mathrm{max}}^{3}(k_{\mathrm{max}}/k_{\mathrm{min}})^{-2\Delta}$ | $\ell_{\mathrm{max}}^{2}(\ell_{\mathrm{max}}/\ell_{\mathrm{min}})^{-4\Delta/3}$
$\Delta=0$ | $k_{\mathrm{max}}^{3}\log(k_{\mathrm{max}}/k_{\mathrm{min}})$ | $\ell_{\mathrm{max}}^{2}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{min}})$
$\Delta>0$ | $k_{\mathrm{max}}^{3}$ | $\ell_{\mathrm{max}}^{2}$
Table 2: General scaling behavior in the squeezed limit for the $(S/N)^{2}$ in
3D and 2D for $N$-point correlation functions.

The logarithmic enhancement of the scaling of the squeezed local non-
Gaussianity is a well-known result in the literature, however the effect of
the squeezed regime can be appreciated more at large
$\ell\gg\ell_{\mathrm{D}}$.

As already seen in Sec. 3.1, as we go to higher multipole values, the effect
of diffusion damping along the line of sight affects the $S/N$ scaling (Eq.
(3.5)). This result has immediate consequences for the statistics of the
perturbations, limiting our capacity of extracting information on non-
Gaussianity from CMB maps.

Given the great interest in such observables, we take the squeezed limit of
the damped angular bispectrum and trispectrum, before generalizing to the
$(N-1)$-spectrum. We show that the signal greatly benefits from the squeezed
limit, especially in the case of the bispectrum.

Using Eq. (3.5) for $N=3$, we see that the general scaling estimation is given
by

$\left(\frac{S}{N}\right)^{2}_{(3)}\sim\ell_{\mathrm{D}}\int^{\ell_{\mathrm{max}}}d^{4}\ell\,\frac{\ell^{-12}}{\ell^{-9}}\sim\ell_{\mathrm{max}},$
(3.11)

thus lower than mode-counting. Let us now take the squeezed bispectrum with
the addition of the transfer function111111Here, we neglect the exponential
factor, given that it cancels out in the $S/N$ computation.

$\lim_{\ell_{1}\ll\ell_{2}\sim\ell_{3}}B(\ell_{1},\ell_{2},\ell_{3})\propto\frac{\ell_{\mathrm{D}}^{2}}{\ell_{1}^{3}\ell_{2}^{3}}\times\left(\frac{\ell_{1}}{\ell_{2}}\right)^{\Delta}.$
(3.12)

Then, the $(S/N)^{2}$ scaling becomes

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(3)}\sim\ell_{\mathrm{D}}\int_{\ell_{\mathrm{D}}}^{c\ell_{\mathrm{max}}}\mathrm{d}\ell_{1}\ell_{1}^{2\Delta-2}\int^{\ell_{\mathrm{max}}}_{c\ell_{\mathrm{max}}}\mathrm{d}\ell_{2}\ell_{2}^{1-2\Delta},\end{split}$
(3.13)

where we set the lower limit of integration to the damping scale
$\ell_{\mathrm{D}}$. We notice that there exist an interval of values
$0<\Delta<1/2$, where the scaling with the number of modes improves with
respect to $\ell_{\mathrm{max}}$ as

$\left(\frac{S}{N}\right)^{2}_{(3)}\sim\ell_{\mathrm{max}}\left(\frac{\ell_{\mathrm{max}}}{\ell_{\mathrm{D}}}\right)^{1-2\Delta}.$
(3.14)

In particular, in the case of local non-Gaussianity, we recover mode-counting
scaling $\sim\ell_{\mathrm{max}}^{2}$ [18]. On the contrary, for $\Delta>1/2$,
the squeezed region becomes subdominant with respect to the damping scale and
the scaling reduces to Eq. (3.11).

These results can be generalized to higher $(N-1)$-spectrum. In Eq. (3.6), we
showed that for the trispectrum, the $(S/N)^{2}$ scaling becomes worse with
respect to the bispectrum result, being proportional to
$\sim\log{\ell_{\mathrm{max}}}$. Performing a similar integration as done for
the bispectrum, we have that the squeezed limit of the trispectrum is still
less sensitive to non-Gaussianity than the bispectrum, as shown in Tab. 3.

As for the $(N-1)$-spectrum, denoting with $\ell=\ell_{2},\dots,\ell_{N}$, we
obtain

$\left(\frac{S}{N}\right)^{2}_{(N)}=\int
d^{2(N-2)}\ell\,d^{2}\ell_{1}\,\frac{\ell_{1}^{-6+2\Delta}\ell^{-6(N-2)-2\Delta}}{\ell_{1}^{-3}\ell^{-3(N-1)}}\sim\ell_{D}^{-1+2\Delta}\ell_{\rm
max}^{5-N-2\Delta}\,,\quad\text{for }0<\Delta<1/2.$ (3.15)

Our analysis has a potential caveat — we assumed that the scaling of the
$N$-point function in the squeezed limit is “equilateral”-like for the
remaining hard momenta. The resulting $S/N$ can be an underestimate, if there
is large signal when taking multiple squeezed limits. An example of interest,
which we analyze numerically later on, is of the local trispectrum of
$g^{\mathrm{local}}_{\mathrm{NL}}$ form,

$T_{\zeta}^{g_{\mathrm{NL}}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})=\frac{54}{25}g^{\mathrm{local}}_{\mathrm{NL}}[P_{\zeta}(k_{2})P_{\zeta}(k_{3})P_{\zeta}(k_{4})+P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(k_{4})+\text{2
perms.}]$ (3.16)

In that case, if we take a squeezed limit $k_{1}\ll k_{2,3,4}$ then the
resulting function has permutations that look identical to the local
bispectrum. In particular, we can consider the double squeezed limit $k_{1}\ll
k_{2}\ll k_{3,4}$, in which case our general estimates add up and we obtain

$\left(\frac{S}{N}\right)^{2}_{g_{\rm NL}}\sim\int
d^{2}\ell_{1}\ell_{1}^{-3}\int d^{2}\ell_{2}\ell_{2}^{-3}\int
d^{2}\ell_{3}\sim\ell_{\rm max}^{2}\,.$ (3.17)

In App. E, we generalize the idea of multiple squeezed limits to the damped
$(N-1)$-spectrum and show that the scaling never exceeds mode-counting. We
conclude that the squeezed bispectrum likely remains the best observable for
non-Gaussianity, especially to test light particles in the early Universe.

| Squeezed limit $\ell\gg\ell_{\mathrm{D}}$
---|---
| Bispectrum | Trispectrum |  $(N-1)$-spectrum
$\Delta<1/2$ | $\ell_{\mathrm{max}}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})^{1-2\Delta}$ | $\ell_{\mathrm{D}}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})^{1-2\Delta}$ | $\ell_{\mathrm{D}}^{N-2}\ell_{\mathrm{max}}^{4-N}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}}^{N-2})^{1-2\Delta}$
$\Delta=1/2$ | $\ell_{\mathrm{D}}^{2}\ell_{\mathrm{max}}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})$ | $\ell_{\mathrm{D}}^{2}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})$ | $\ell_{\mathrm{D}}^{2}\ell_{\mathrm{max}}^{4-N}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})$
$\Delta>1/2$ | $\ell_{\mathrm{max}}$ | $\ell_{\mathrm{D}}^{2}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})^{1-2\Delta}$ | $\ell_{\mathrm{D}}^{N-2}\ell_{\mathrm{max}}^{4-N}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}}^{N-2})^{1-2\Delta}$
Table 3: Scaling of the $(S/N)^{2}$ in the squeezed limit for $N$-point
correlators, when the damping effect is dominant.

###  3.3 The collapsed limit

Following the same steps of Sec. 3.2, in this section we derive the
theoretical estimations for the collapsed limit of $N>3$ correlators. We show
that, for certain values of $\Delta$, the collapsed limit provides a better
scaling than mode-counting, for both 3D and 2D estimates.

Using the same notation as Sec. 2, we denote the “left” momenta
$\\{k_{1},\cdots k_{M}\\}\sim k_{L}$ and the “right” momenta
$\\{k_{M+1},\cdots k_{N}\\}\sim k_{R}$, with $M>1$. In this limit, the
diagonal $k_{I}\equiv k_{1}+\cdots+k_{M}$ is much smaller than $k_{L},k_{R}$.
Then, replacing the $(N-1)$-spectrum defined in Eq. (2.4) into the volume
survey $(S/N)^{2}$, we obtain

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(N)}\sim\int_{k_{\mathrm{min}}}^{ck_{\mathrm{max}}}\mathrm{d}^{3}k_{I}\int_{ck_{\mathrm{max}}}^{k_{\mathrm{max}}}\mathrm{d}^{3A}k_{R}\mathrm{d}^{3B}k_{\mathrm{L}}\,\frac{k_{R}^{-6A-2\Delta}k_{L}^{-6B-2\Delta}k_{I}^{-6+4\Delta}}{k_{R}^{-3(A+1)}k_{L}^{-3(B+1)}}\sim
k_{\mathrm{min}}^{4\Delta-3}k_{\mathrm{max}}^{6-4\Delta},\end{split}$ (3.18)

where $A=N-M-1$ and $B=M-1$. We find that for values $\Delta<3/4$ the scaling
increases with respect to mode-counting as

$\left(\frac{S}{N}\right)^{2}_{(N)}\sim
k_{\mathrm{max}}^{3}\left(\frac{k_{\mathrm{max}}}{k_{\mathrm{min}}}\right)^{3-4\Delta},\quad\Delta<3/4.$
(3.19)

Along this lines, given the collapsed angular $(N-1)$-spectrum

$\lim_{\ell_{I}\ll\ell_{L},\ell_{R}}F_{\zeta}(\ell_{I},\ell_{L},\ell_{R})\propto\frac{1}{\ell_{R}^{2A}\ell^{2B}_{L}\ell^{2}_{I}}\left(\frac{\ell_{I}^{2}}{\ell_{R}\ell_{L}}\right)^{2\Delta/3},$
(3.20)

we obtain
$(S/N)_{(N)}^{2}\sim\ell_{\mathrm{min}}^{8\Delta/3-2}\ell_{\mathrm{max}}^{4-8\Delta/3}$.
As for the 3D survey, $\Delta=3/4$ represents a threshold values below which
the scaling improves as

$\left(\frac{S}{N}\right)^{2}_{(N)}\sim\ell_{\mathrm{max}}^{2}\left(\frac{\ell_{\mathrm{max}}}{\ell_{\mathrm{min}}}\right)^{2-8\Delta/3},\quad\Delta<3/4.$
(3.21)

As a result, even if our estimates are highly idealized and optimistic,
collapsed $(N-1)$-spectra might be very sensitive to non-Gaussianity arisen
from a light particle driven inflationary process. In particular, local type
of non-Gaussianity produces a $\ell_{\mathrm{max}}^{4}$ scaling
($k_{\mathrm{max}}^{6}$ for a 3D survey), as was found in Refs. [19, 20]. We
summarize the results in Tab. 4.

| Collapsed limit $\ell\ll\ell_{\mathrm{D}}$
---|---
| 3D | 2D ($\ell\ll\ell_{\mathrm{D}}$)
$\Delta<3/4$ | $k_{\mathrm{max}}^{3}(k_{\mathrm{max}}/k_{\mathrm{min}})^{3-4\Delta}$ | $\ell_{\mathrm{max}}^{2}(\ell_{\mathrm{max}}/\ell_{\mathrm{min}})^{2-8\Delta/3}$
$\Delta=3/4$ | $k_{\mathrm{max}}^{3}\log(k_{\mathrm{max}}/k_{\mathrm{min}})$ | $\ell_{\mathrm{max}}^{2}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{min}})$
$\Delta>3/4$ | $k_{\mathrm{max}}^{3}$ | $\ell_{\mathrm{max}}^{2}$
Table 4: General scaling behavior in the collapsed limit for the $(S/N)^{2}$
in 3D and 2D for $N$-point correlation functions.

Given the enhanced scaling in the collapsed limit for some values of $\Delta$,
and in analogy with the squeezed bispectra, we might expect these shapes to be
only mildly affected by damping in the CMB. As an example, let us consider the
angular damped trispectrum in the collapsed limit

$T(\ell_{1},\ell_{2},\ell_{3},\ell_{4})\sim\frac{\ell_{\mathrm{D}}^{3}}{(\ell_{1}\ell_{3}\ell_{I})^{3}}\left(\frac{\ell_{I}^{2}}{\ell_{1}\ell_{3}}\right)^{\Delta}.$
(3.22)

In this case, we find
$(S/N)_{(4)}^{2}\sim\ell_{\mathrm{D}}^{2}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})^{4(1-\Delta)}$,
therefore trispectra shapes in the interval $0<\Delta<1$ have a favorable
scaling despite the effect of damping. The fact that for $\Delta=0$ yields a
scaling of $\ell_{\mathrm{max}}^{4}$, as in the $\ell\ll\ell_{\mathrm{D}}$
scenario, indicates that damping has little influence on such highly collapsed
trispectra.

These results can easily be extended to general collapsed $N$-point
correlation function

$\lim_{\ell_{I}\ll\ell_{L},\ell_{R}}F_{\zeta}(\ell_{I},\ell_{L},\ell_{R})\propto\frac{\ell_{\mathrm{D}}^{N-1}}{\ell_{R}^{3A}\ell^{3B}_{L}\ell^{3}_{I}}\left(\frac{\ell_{I}^{2}}{\ell_{R}\ell_{L}}\right)^{\Delta}.$
(3.23)

Tab. 5 summarizes the results. Similar to having nested squeezed limits, there
is the possibility of nested collapsed limits, which we investigate in App. E.
In the single collapsed limit scenario, for a massless particle ($\Delta=0$)
there is still potential to gain signal-to-noise when moving towards higher
resolution (higher $\ell_{\mathrm{max}}$), as long as $N<8$, whereas in the
general case this enhancement is dependent on the specific number of nested
collapsed limits. In the section below, we discuss in more depth the extent of
these optimistic scalings in the specific example of the trispectrum.

| Collapsed limit $\ell\gg\ell_{\mathrm{D}}$
---|---
| Trispectrum |  $(N-1)$-spectrum
$\Delta<1$ | $\ell_{\mathrm{D}}^{2}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}}^{2})^{4(1-\Delta)}$ | $\ell_{\mathrm{D}}^{N-2}\ell_{\mathrm{max}}^{4-N}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}}^{N-2})^{4(1-\Delta)}$
$\Delta=1$ | $\ell_{\mathrm{D}}^{2}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})$ | $\ell_{\mathrm{D}}^{2}\ell_{\mathrm{max}}^{4-N}\log(\ell_{\mathrm{max}}/\ell_{\mathrm{D}})$
$\Delta>1$ | $\ell_{\mathrm{D}}^{2}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}}^{2})^{4(1-\Delta)}$ | $\ell_{\mathrm{D}}^{N-2}\ell_{\mathrm{max}}^{4-N}(\ell_{\mathrm{max}}/\ell_{\mathrm{D}}^{N-2})^{4(1-\Delta)}$
Table 5: The table shows the scaling estimations for 2D in the collapsed limit
when the damping effect is dominant.

####  3.3.1 A comment on the enhanced scaling

The improved scaling of the $(N-1)$-spectra for $N>3$ in the (nested)
collapsed limit immediately raises the question whether trispectrum (and
beyond) measurements can provide more precise constraints on the primordial
interaction couplings. Here we will limit the discussion to the CMB
trispectrum, compared to the bispectrum, but the arguments can be applied more
generally to $(N-1)$-spectra, including nested collapsed limits, and to LSS
surveys in 3D. Earlier work [19] had noticed the enhanced trispectrum scaling
for local type non-Gaussianity and the question of whether the trispectrum can
better constrain $f_{\mathrm{NL}}$ within that model has been extensively
examined in Ref. [21]. In that work the authors consider the local model, so
named as the associated non-Gaussianity can be expressed using the following
local modifications to the linear perturbation

$\displaystyle\Phi(\mathbf{x})=\Phi_{L}(\mathbf{x})+f_{\mathrm{NL}}(\Phi_{L}^{2}(\mathbf{x})-\langle\Phi_{L}^{2}(\mathbf{x})\rangle)+g_{\mathrm{NL}}\Phi_{L}^{3}(\mathbf{x}),$
(3.24)

where $\Phi_{L}(\mathbf{x})$ is the linear, Gaussian perturbation. They show
that when the $S/N$ of the trispectrum starts to exceed that of the
bispectrum, additional non-Gaussian terms become important, such that the
bispectrum signal-to-noise is always larger than that of the trispectrum.
Heuristically, this can be understood knowing that in the local model the
trispectrum amplitude is $\mathcal{O}(f_{\mathrm{NL}}^{2})$ and the non-
Gaussian contributions to the trispectrum variance are of the same order.
Thus, when the trispectrum signal exceeds the Gaussian variance, so do the
non-Gaussian contributions to the variance, completely canceling out enhanced
scaling.

Eq. (3.24) has two leading contributions when computing the 4-point
correlation function. The first contribution is coming from the cubic term
$\propto g_{\rm NL}$. By definition, the associated
$(S/N)^{2}\propto\ell_{\mathrm{max}}^{2}$. Dimensionally, we can assume that
$g_{\mathrm{NL}}=f_{\mathrm{NL}}^{2}\alpha$. Therefore the only way for this
contribution to be observed with $S/N$ exceeding that of the bispectrum is to
have sufficiently large values of $\alpha$ ($\alpha>1/(f_{\rm
NL}\sqrt{A_{s}})$). The second contribution in the model has favorable scaling
in the collapsed limit, i.e. $(S/N)^{2}\propto\ell_{\mathrm{max}}^{4}$. We
generally refer to the amplitude of this contribution with
$\tau_{\mathrm{NL}}$. However, in this model $\tau_{\mathrm{NL}}\propto
f_{\mathrm{NL}}^{2}$, with no possibility to fine-tune the overall scaling
factor ($\alpha=6$). As a result, even though the $S/N$ could potentially
exceed that of the bispectrum, it requires a very large value of $\alpha$ and
the scaling of the signal-to-noise of the trispectrum does not improve over
simple mode counting. The fact that $\tau_{\mathrm{NL}}\propto
f_{\mathrm{NL}}^{2}$ also guarantees that a trispectrum will not have a
favorable scaling with respect to the signal $f_{\rm NL}$ compared to the
bispectrum. In order to exploit the improved scaling the contribution to
$\tau_{\mathrm{NL}}$ should become independent from the value of
$f_{\mathrm{NL}}$. A straightforward way to achieve this is by introducing an
additional field that couples to the inflaton.

In Ref. [31], this idea is generalized. The authors point out that even when
the bispectrum estimator is an optimal estimator (one which satisfies the
Cramér-Rao bound and is therefore the minimum possible unbiased estimator)
trispectrum measurements can add statistically independent information! This
is possible as trispectrum estimators that measure $f_{\mathrm{NL}}^{2}$ are
biased estimators for $f_{\rm NL}$, and bispectrum estimators are suboptimal
estimators for $f_{\mathrm{NL}}^{2}$. Further, they show that trispectra,
whose $S/N$ is dominated by a small number of shapes that can be expressed as
two triangles (see Fig. 3 in Ref. [32]), are highly correlated with the
bispectrum and thus contain minimal new information.

From these papers we draw two conclusions. First, there is no general theorem
that prevents trispectra (and higher point correlation functions) from
providing more constraining power than bispectra on primordial couplings.
Second, the $S/N$ computations above, which neglect non-Gaussian
contributions, can be too naive, particularly for collapsed $(N-1)$-spectra.
For strongly collapsed trispectra, the non-Gaussian contributions can dominate
the variance, reducing (or even completely negating) the information in the
trispectrum. This conclusion is however model dependent and when the relation
between the trispectrum and bispectrum amplitude deviates from the relation
demonstrated for the local model, as in Ref. [20], the trispectrum can be
highly informative and constraining.

##  4 The CMB bi- and tri-spectrum

In this section we present a numerical analysis of the $S/N$ for the
bispectrum and the trispectrum of the temperature and polarization
anisotropies. We assume measurements are cosmic variance limited and assume
unlensed spectra (see e.g. [12] for the effect of lensing and how to mitigate
it). The analysis in this section will use the full radiation transfer
functions, which we calculate numerically using the Boltzmann solver CAMB
[33]. We aim to confirm the heuristic results of Sec. 3 in the damped regime
$\ell\gg\ell_{\mathrm{D}}$. Tab. 1 summarizes the flat $\Lambda$CDM cosmology
parameters used.

###  4.1 Bispectrum

In this subsection, we aim to compute

$\left(\frac{S}{N}\right)^{2}_{(3)}=\sum_{X_{i},X^{\prime}_{i}}\sum_{\ell_{1}\leq\ell_{2}\leq\ell_{3}\leq\ell_{\mathrm{max}}}\frac{h^{\ell_{1}\ell_{2}\ell_{3}}}{f(\ell_{1},\ell_{2},\ell_{3})}b^{X_{1}X_{2}X_{3}}_{\ell_{1}\ell_{2}\ell_{3}}\left(C_{\ell_{1}}^{X_{1}X^{\prime}_{1}}\right)^{-1}\left(C_{\ell_{2}}^{X_{2}X^{\prime}_{2}}\right)^{-1}\left(C_{\ell_{3}}^{X_{3}X^{\prime}_{3}}\right)^{-1}b^{X^{\prime}_{1}X^{\prime}_{2}X^{\prime}_{3}}_{\ell_{1}\ell_{2}\ell_{3}},$
(4.1)

up to $\ell_{\mathrm{max}}=5000$. Here $X_{i}=\\{T,E\\}$ where $T,\,E$ refer
respectively to temperature and $E$-mode polarization.

The general functional form of the primordial bispectrum can render the
computation of the CMB bispectrum challenging. To mitigate this issue,
separable templates, which approximate the theoretical predictions and are
suitable for data analysis, have been proposed. In this subsection, we will
investigate the local [34], equilateral [35] and orthogonal templates [36].

The local template has a corresponding value $\Delta=0$ and is given by

$B^{\mathrm{local}}_{\zeta}(k_{1},k_{2},k_{3})=\frac{6}{5}f_{\mathrm{NL}}^{\mathrm{local}}\bigg{[}P_{\zeta}(k_{1})P_{\zeta}(k_{2})+\text{2
perms.}\bigg{]}.$ (4.2)

Here
$P_{\zeta}(k)=2\pi^{2}A_{\mathrm{s}}(k/k_{\star})^{n_{\mathrm{s}}-1}k^{-3}$,
where $A_{\mathrm{s}}$ the amplitude of initial fluctuations, $n_{\mathrm{s}}$
is the scalar spectral index (see e.g. Tab. 1) and
$k_{\star}=0.05\;\mathrm{Mpc^{-1}}$ is the pivot scale. This shape is largest
in the squeezed limit, i.e. $k_{1}\ll k_{2}\sim k_{3}$, therefore according to
our analytical estimates in Sec. 3.2, we expect the scaling of the $(S/N)^{2}$
to be proportional to $\ell_{\mathrm{max}}^{2}$.

Equilateral non-Gaussianity can be captured using the following template

$\displaystyle B^{\mathrm{equil}}_{\zeta}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{equil}}\left[-P_{\zeta}(k_{1})P_{\zeta}(k_{2})-\text{2
perms.}-2P^{2/3}_{\zeta}(k_{1})P^{2/3}_{\zeta}(k_{2})P^{2/3}_{\zeta}(k_{3})+\right.$
(4.3)
$\displaystyle\left.P^{1/3}_{\zeta}(k_{1})P^{2/3}_{\zeta}(k_{2})P_{\zeta}(k_{3})+\text{5
perms.}\right],$

and peaks when $k_{1}=k_{2}=k_{3}$. It is straightforward to check that this
template has $\Delta=2$, and we expect $(S/N)^{2}$ to be proportional to
$\ell_{\mathrm{max}}$.

A third shape, orthogonal to the equilateral template, was introduced in Ref.
[36] and it is parametrized by

$\displaystyle B^{\mathrm{ortho}}_{\zeta}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{ortho}}\left[-3P_{\zeta}(k_{1})P_{\zeta}(k_{2})-\text{2
perms.}-8P^{2/3}_{\zeta}(k_{1})P^{2/3}_{\zeta}(k_{2})P^{2/3}_{\zeta}(k_{3})+\right.$
(4.4)
$\displaystyle\left.3P^{1/3}_{\zeta}(k_{1})P^{2/3}_{\zeta}(k_{2})P_{\zeta}(k_{3})+\text{5
perms.}\right].$

This shape peaks in both equilateral and flattened configurations, i.e.
$k_{1}=k_{2}+k_{3}$. The orthogonal template above has $\Delta=1$, however as
shown in Ref. [36], Eq. (4.4) is an approximation of a more numerically
challenging shape and does not have the correct scaling in the squeezed limit
($\Delta=2$). For our purposes however it suffices since it nicely completes a
set of templates with $\Delta=0$, 1, and 2. Our analytical estimate suggests
that $(S/N)^{2}$ should be proportional $\ell_{\rm max}$ because $\Delta>1/2$.
Our numerical analysis however will show that the orthogonal template produces
an improved scaling, at least over the range we explored
($\ell_{\mathrm{max}}=5000$). We will comment on this discrepancy in some
depth below.

In Fig. 6, we show the $(S/N)^{2}$ for only temperature, only $E$-mode
polarization, and both temperature and mode polarization anisotropies.

Figure 6: Numerical solution of the bispectrum $(S/N)^{2}$ given in Eq. (4.1)
in the range $\ell_{D}\leq\ell_{\mathrm{max}}\leq 5000$, where $\ell_{D}\simeq
1300$. Upper left panel: $(S/N)^{2}$ for the local shape for temperature,
$E$-mode polarization and in combination. Upper right and lower panel:
$(S/N)^{2}$ for the equilateral and orthogonal shapes respectively. Curves
shown are for $f_{\rm NL}=1$.

Guided by our analytical results in Sec. 3, we expect the $(S/N)^{2}$ to be
proportional to $\ell_{\mathrm{max}}^{p}$. Numerically, we can extract the
exponent by taking $\log$-derivative of the $(S/N)^{2}$, i.e.,

$p=\frac{\mathrm{d}}{\mathrm{d}\ell_{\mathrm{max}}}\log{[(S/N)^{2}]}\equiv\frac{\ell_{\mathrm{max}}}{(S/N)^{2}}\frac{\mathrm{d}(S/N)^{2}}{\mathrm{d}\ell_{\mathrm{max}}}.$
(4.5)

We are mostly interested in $p$ in the damped regime, so it suffices to fit
the numerical exponent with a constant based on the theoretical scaling in the
range $\ell_{\mathrm{max}}>\ell_{\mathrm{D}}$, where $\ell_{\mathrm{D}}\simeq
1300$. In Tab. 6 we report our results together with the one standard
deviation error.

| Numerical Scaling Bispectrum
---|---
| Local | Equilateral | Orthogonal
$T$ | $2.0169\pm 0.0003$ | $0.912\pm 0.003$ | $1.544\pm 0.002$
$E$ | $2.005\pm 0.002$ | $0.936\pm 0.004$ | $1.326\pm 0.004$
$T+E$ | $2.0650\pm 0.0004$ | $1.170\pm 0.003$ | $1.658\pm 0.001$
Table 6: Values of the numerical scaling $p$ inferred from Eq. (4.5) using a
least square fit and the standard deviation error. The numerical derivatives
are derived in the range $\ell_{D}\leq\ell_{\mathrm{max}}\leq 5000$. These
values are to be compared to the theoretical ones in Tab. 3.

The numerical exponents are in good agreement with theory in Tab. 3, except
for the orthogonal shape which shows a better scaling than theoretical
expectations in this range. While we do not have an analytical argument to
explain this behavior, we can derive some intuition from computing the $S/N$
in the flat-sky numerically using Eq. (D.20) and damping as in Eq. (C.1). We
then find (not shown) that while our analytical results suggests that the
effect of damping should appear for $\ell_{\mathrm{max}}>\ell_{D}$ for shapes
$\Delta>\Delta_{c}$ (we define $\Delta_{c}$ as the critical value for which
the scaling changes, i.e. $\Delta_{c}=1/2$ for the bispectrum), depending on
the value of $\Delta$ we find that equilateral and orthogonal shapes do not
converge equally towards the limiting scaling for the same
$\ell_{\mathrm{max}}$, where the orthogonal shape (but also the flat shape,
see e.g. [37] Eq. (21), which also has $\Delta=1$) converges more slowly. This
could suggest that the scaling predicted by our analytical analysis is not yet
captured when limiting our numerical analysis to $\ell_{\mathrm{max}}\leq
5000$. Another explanation could be that our heuristic analytical derivation
does not capture the details of the full radiative transfer function,
resulting in intermediate scalings for $\Delta$ close to the critical value
$\Delta_{c}$. The main point however, as confirmed by our numerical analysis,
is that for $\Delta>\Delta_{c}$ (equilateral-like shapes) the scaling is
affected much more severely by damping then for shapes where
$\Delta<\Delta_{c}$ (local-like shapes).

Our analytical calculations were strictly derived as an approximation of the
temperature transfer functions (see App. C), but polarization fluctuations are
similarly affected by the thickness of last scattering, evident from our
numerical calculations. Combining the temperature and polarization signatures
can lead to improvement of the scaling as can be read off from Tab. 6 for the
orthogonal and equilateral shapes. While we do not have a full qualitative
derivation why this happens, we can make the following heuristic argument. The
reduced scaling for shapes with $\Delta>\Delta_{c}$ is caused by the blurring
of the modes on small scales. At the level of the radiative transfer function,
the temperature and polarization modes oscillate (around zero) and are out of
phase. When combining polarization and temperature measurements, where the
zero points of each individual transfer function yielded zero response and a
loss of primordial signal, the product of these transfer functions can now
limit some of these losses around the zero points of the transfer function
(the zero points are recovered only in the auto correlation of the fields).
The net effect is a ‘de-blurring’ of the last scattering surface. While this
leads to immediate improvements on the scaling of spectra with
$\Delta>\Delta_{c}$, since blurring does very little to spectra with
$\Delta<\Delta_{c}$, no change in the scaling is anticipated. Hence the
scaling of the $S/N$ for e.g. the local shape does not benefit from adding
polarization.

###  4.2 Trispectrum

In this subsection, we will perform a numerical Fisher analysis for the
trispectrum. We will focus on the separable form of the local primordial
trispectrum [19, 38, 32]121212For a semi-analytic method to compute
cosmological angular trispectra related to the discussion on cosmological
correlators in Sec. 2, see e.g. [39].

$T^{\mathrm{local}}_{\zeta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4},\bm{K})=T^{\tau_{\mathrm{NL}}}_{\zeta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4},\bm{K})+T^{g_{\mathrm{NL}}}_{\zeta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4},\bm{K}),$
(4.6)

where

$\displaystyle
T^{\tau_{\mathrm{NL}}}_{\zeta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4},\bm{K})$
$\displaystyle=$
$\displaystyle\tau_{\mathrm{NL}}[P_{\zeta}(K)P_{\zeta}(k_{1})P_{\zeta}(k_{3})+\text{11
perms.}],$ (4.7) $\displaystyle
T^{g_{\mathrm{NL}}}_{\zeta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4},\bm{K})$
$\displaystyle=$
$\displaystyle\frac{54}{25}g^{\mathrm{local}}_{\mathrm{NL}}[P_{\zeta}(k_{2})P_{\zeta}(k_{3})P_{\zeta}(k_{4})+P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(k_{4})+\text{2
perms.}],$

with $\tau_{\mathrm{NL}}\geq(6/5f^{\mathrm{local}}_{\mathrm{NL}})^{2}$ [40].
Both the $\tau_{\rm NL}$ and the $g_{\rm NL}$ templates arise from local type
of non-Gaussianity. However, here we consider them separately because our
analytical analysis suggest these trispectra should have a different scaling
with $\ell_{\mathrm{max}}$. The $\tau_{\mathrm{NL}}$-trispectrum is largest in
the collapsed limit with $\Delta=0$, and we expect
$(S/N)^{2}\propto\ell_{\mathrm{max}}^{4}$. On the other hand, as shown in Sec.
3.2, the $g_{\mathrm{NL}}$-trispectrum is largest in the double squeezed
limit, producing a scaling of $(S/N)^{2}\propto\ell_{\mathrm{max}}^{2}$.

We compute the $S/N$

$\displaystyle\left(\frac{S}{N}\right)^{2}_{(4)}$ $\displaystyle=$
$\displaystyle\sum_{X_{i},X^{\prime}_{i}}\sum\limits_{\ell_{i}}\sum\limits_{L}\frac{1}{(2L+1)}T^{X_{1}X_{2}X_{3}X_{4}}_{c,\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)\left(\mathbb{C}^{-1}\right)_{\ell_{1}}^{X_{1}X^{\prime}_{1}}\left(\mathbb{C}^{-1}\right)_{\ell_{2}}^{X_{2}X^{\prime}_{2}}\times$
(4.9)
$\displaystyle\left(\mathbb{C}^{-1}\right)_{\ell_{3}}^{X_{3}X^{\prime}_{3}}\left(\mathbb{C}^{-1}\right)_{\ell_{4}}^{X_{4}X^{\prime}_{4}}T^{X^{\prime}_{1}X^{\prime}_{2}X^{\prime}_{3}X^{\prime}_{4}}_{c,\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L),$

for the angular trispectra in Eq. (4.7) and Eq. (4.2). The computational
requirements are challenging as in general the above equation scales as
$\mathcal{O}(r^{2}\ell_{\rm max}^{4})$. With the reformulation and
approximations described in App. A we show that we can evaluate
$(S/N)^{2}_{\tau_{\mathrm{NL}}}$ with $\mathcal{O}(r\ell_{\rm max}^{4})$ and
$(S/N)^{2}_{g_{\mathrm{NL}}}$ with $\mathcal{O}(r^{2}\ell_{\rm max}^{3})$
computations. For this reason we limit ourselves to
$\ell_{\mathrm{max}}=4000$. When combining temperature and polarization, the
number of computations increases by a factor of 16 for $\tau_{\rm NL}$ and 2
for $g_{\rm NL}$. We therefore limit our analysis of the combined temperature
and polarization modes to $\ell_{\mathrm{max}}=2000$ and
$\ell_{\mathrm{max}}=4000$ for $\tau_{\mathrm{NL}}$ and $g_{\rm NL}$
respectively. The $g_{\mathrm{NL}}$-angular trispectrum does not explicitly
dependent on the diagonal mode $L$, while the $\tau_{\mathrm{NL}}$-angular
trispectrum does and we compute $S/N$ for the $\tau_{\mathrm{NL}}$-angular
trispectrum up to $L_{\rm max}=10$. We varied $L_{\rm max}$ and found it did
not change our results significantly, which is in accordance with the findings
in Ref. [19]. Further details of the computation can be found in App. A. The
results are shown in Fig. 7.

Figure 7: Numerical solution of the trispectrum $(S/N)^{2}$, given in Eq.
(D.17) in the range $\ell_{D}\leq\ell_{\mathrm{max}}\leq 4000$, where
$\ell_{D}\simeq 1300$. Left panel: $(S/N)^{2}$ for the
$\tau_{\mathrm{NL}}$-trispectrum shape for temperature, $E$-mode polarization
and in combination. The $T+E$ solution is shown in the range
$\ell_{D}\leq\ell_{\mathrm{max}}\leq 2000$ and the dashed curve refers to the
expected solution at higher multipoles. Notice that the $T$ and $E$ show
identical contribution to the $S/N$. This is expected because for the
trispectrum in the collapsed limit the two triangles in the limit of small $L$
will be squeezed and their contribution to the trispectrum is expected to be
proportional to $\prod_{i}^{4}C^{X}_{\ell_{i}}$. This cancels exactly with the
cosmic variance noise if considering only $T$ or $E$. Right panel: $(S/N)^{2}$
$g_{\mathrm{NL}}$-trispectrum shape. Curves shown are for $\tau_{\rm NL}=1$
and $g_{\rm NL}=1$.

We repeat the analysis done for the bispectrum and extract the scaling of the
$(S/N)^{2}$ using Eq. (4.5) in the damped regime
$\ell_{\mathrm{max}}>\ell_{\mathrm{D}}$. The results are summarized in Tab. 6.
We find a scaling for the $(S/N)^{2}_{\tau_{\rm
NL}}\propto\ell_{\mathrm{max}}^{4}$ in accordance with our analytical
forecast. When combining temperature and polarization, the scaling gets a
little worse, but the change is caused by the limited range in
$\ell_{\mathrm{max}}$ over which the scaling is determined. If we limit
ourselves in the same range in temperature, we obtain a consistent scaling.
Our analysis shows that the scaling of $(S/N)^{2}_{\tau_{\rm NL}}$ is
saturated for the damped CMB, and the scaling can not be improved when
combining temperature and polarization. This confirms the results we found for
the bispectrum, where the local bispectrum does not appear to benefit when
combining the data, because damping has a little effect on the scaling. Tab. 4
suggests that there is no difference between the damped an undamped tracers in
the collapsed limit when considering a shape with $\Delta=0$ for the
trispectrum.

| Numerical Scaling Trispectrum
---|---
| Local $\tau_{\mathrm{NL}}$ | Local $g_{\mathrm{NL}}$
$T$ | $4.105\pm 0.001$ | $2.16\pm 0.07$
$E$ | $4.136\pm 0.002$ | $2.10\pm 0.05$
$T+E$ | $4.070\pm 0.001$ | $2.23\pm 0.05$
Table 7: Values of the numerical scaling $p$ inferred as in Eq. (4.5) using a
least square fit, together with the standard deviation error. The numerical
derivatives are derived in the range $\ell_{D}\leq\ell_{\mathrm{max}}\leq
4000$, with the exception of combined temperature and polarization data for
the $\tau_{\mathrm{NL}}$-trispectrum which is computed in
$\ell_{D}\leq\ell_{\mathrm{max}}\leq 2000$. The $T+E$ scaling for the
$\tau_{\mathrm{NL}}$-trispectrum is lower that the temperature one because the
range considered is smaller. Indeed, when considering
$\ell_{D}\leq\ell_{\mathrm{max}}\leq 2000$, we find the temperature
$(S/N)^{2}$ has $p=4.058$. The error on the $g_{\mathrm{NL}}$-trispectrum is
large due to the few data points used to infer the scaling.

For $g_{\rm NL}$ we obtain a scaling close to $(S/N)^{2}_{g_{\rm
NL}}\propto\ell_{\mathrm{max}}^{2}$. The numerical scaling is a little higher
than the scaling predicted by our analytical analysis, but has a fairly large
std, because of the low number of samples we were able to compute. The scaling
is not expected to benefit much from adding temperature and polarization
measurements as we do not expect the scaling of the $S/N$ of the squeezed
bispectrum to exceed the scaling of the squeezed trispectrum. Hence, at best
for an undamped CMB we anticipate $(S/N)^{2}_{g_{\rm
NL}}\propto\ell_{\mathrm{max}}^{2}\log\ell_{\mathrm{max}}/\ell_{\mathrm{min}}$.

##  5 Discussion and conclusions

In this paper we built a quantitative and intuitive understanding of scaling
(in inverse units of resolution) of the signal-to-noise of primordial
$(N-1)$-spectra. Statistical fluctuations in the field that sources the
observed density and radiation fluctuations in the universe have the potential
to provide evidence for the dynamics of inflation and allow us to probe new
physics through relic signatures of new particles that modify these
$(N-1)$-spectra. From a theoretical point of view, future constraints on these
spectra and/or the detection of any of these signatures would be extremely
valuable, on par with the discovery of the constituent particles of the
standard model, like the Higgs. It is therefore of immediate relevance how we
expect future constraints to improve as a function of resolution. In this work
we considered the ideal limits in the noiseless case, while assuming that
there are no other sources of non-Gaussianity that could hinder a detection
(either via signal confusion or by introducing extra variance). Hence all
results presented in this paper should be considered in this context, i.e.
they present an ideal scenario. With more realistic assumptions, the scalings
derived are likely to get worse and constraints and/or detections will become
more challenging. At the same time, we also limit our analytical estimates to
just a single tracer field, e.g. in the case of the CMB we only consider the
CMB temperature field. As shown in our numerical forecasts, it turns out that
adding polarization improves scalings in the CMB, which would mean that our
analytical findings are pessimistic compared to the case when multiple tracers
are combined.

Naively, the scaling of the signal-to-noise for all spectra grows
proportionally to the highest observed multipole. Our simplified setup gave us
insight into how such scalings come about, and how the naive intuition fails.
Specifically, we have shown that, due to damping and projection effects,
scaling is much poorer, even converging for higher-point correlation
functions. Nonetheless, for certain scalings around squeezed and collapsed
limits, there is a smaller loss of information when modes are damped. The most
important results are summarized in Tab. 8. We have compared our analytic
estimates to a full-sky numerical analysis, including the additional signal of
polarization.

Let us summarize the main points of our analysis:

  * •

The signal-to-noise scaling is reduced due to blurring of the last scattering
surface at short distances. This blurring is caused by the combination of
projection and damping, but the loss of signal is not due to exponential
decay, as both signal and noise are equally damped.

  * •

The behavior of the $N$-point correlator in the squeezed limit or the
collapsed limit for $N>3$ can significantly impact the scaling of the signal-
to-noise with the number of modes, even with a reduced range of momenta
probing this limit.

  * •

Blurring affects equilateral-like $N$-point correlators much more than
squeezed ones.

  * •

We confirm that our analytical estimates capture the scaling in the limit
$\ell\gg\ell_{D}$ after comparing these to CMB forecasts using the full
radiation transfer equations for the bispectrum and trispectrum.

  * •

This comparison also shows that adding polarization data will generally
improve the scaling. The improvement is not observed for spectra that are not
already close to mode-counting, such as squeezed spectra.

  * •

Comparison between analytical and numerical results also shows that while our
analytical results seem to suggest a critical value for $\Delta$ above which
scaling changes, this cutoff is less clear when applying the full radiative
transfer equations. If we numerically compute the signal-to-noise in the flat
sky case using Eq. (D.20) and damping as in Eq. (C.1) we do find a hierarchy
between different shapes at low multipoles, depending on the value of
$\Delta$, but scaling does converge to the critical scaling at high $\ell_{\rm
max}$. We conclude that either we have to go to higher $\ell_{\rm max}$ in our
full radiative transfer numerical calculations to reach critical scaling,
given that our analytical estimates consider scales much larger than the
damping scale, or that the full sky radiative transfer functions introduce
complexity to the scaling that we can not capture with our simplified
estimates for shapes that are near to the critical value $\Delta_{c}$.

  * •

We recover a result that was found earlier that the collapsed limit of the
trispectrum can have a signal-to-noise ratio that scales better with
resolution than the bispectrum. Our analytical analysis has allowed us to
extend this to $(N-1)$-spectra for $N>4$. The enhanced scaling is a function
of $\Delta$ (in general) and $N$ (for damped tracers ). We have elaborated on
this finding in some depth in Sec. 3.3.1. The enhancement saturates above
$N=4$ even for undamped tracers for a single collapsed limit. However for
multiple collapsed limit the enhancement can be larger (see App. E).

  * •

Our analysis assumed that the $S/N$ can be estimated by taking specific
kinematic limits, but ignores the remaining angular dependence of the
$(N-1)$-spectra. While we empirically find this to be valid for the local and
equilateral bispectrum and the local trispectrum, it is possible that our
treatment could miss scalings due to phase space limits not captured by the
equilateral, collapsed or squeezed limits. For example, the folded limit in
the bispectrum $k_{i}+k_{j}\sim k_{k}$ $(i\neq j\neq k)$ could be responsible
for the some of the improved scaling we find for the orthogonal template.
Keeping explicit angular dependence will complicate the analytical
computations and we will leave this for future investigations.

| 2D | 3D | Notes |
---|---|---|---|---
| Undamped | Damped |
Scaling | $\ell_{\rm max}^{2}$ | $\ell_{\rm max}^{4-N}$ | $k_{\rm max}^{3}$ | For $N=4$ and $\ell\gg\ell_{\mathrm{D}}$, $(S/N)^{2}\sim\log\ell_{\rm max}$ |
Squeezed | $\ell_{\rm max}^{2}$ | $\ell_{\rm max}^{5-N-2\Delta_{S}}$ | $k_{\rm max}^{3}$ |  | Undamped $\log$-enhancement for $\Delta_{S}=0$.
---
Damped enhancement for $\Delta_{S}<1/2$.
Collapsed | $\ell_{\rm max}^{4-8\Delta_{C}/3}$ | $\ell_{\rm max}^{8-N-4\Delta_{C}}$ | $k_{\rm max}^{6-4\Delta_{C}}$ |  | Undamped enhancement for $\Delta_{C}<3/4$.
---
Damped enhancement for $\Delta_{C}<1$.
Table 8: Summary of the theoretical estimations of the $(S/N)^{2}$ scaling
provided in Sec. 3. The analysis does not include nested squeezed or collapsed
limits, where the scaling might change, as showed for the
$g_{\mathrm{NL}}$-trispectrum.

An important overarching conclusion is that equilateral-like $(N-1)$-spectra
are harder to constrain, because there is very little signal in the squeezed
and collapsed limits. For an undamped equilateral-like tracer in 3-dimensions,
mode-counting is applicable, but no enhancement is possible in the collapsed
limit. If information is damped, the difference between squeezed and
equilateral shapes becomes even more pronounced. While this loss of
information is guaranteed in the CMB, it is imaginable that something similar
could happen in large scale structure measurements. For example, a promising
avenue to improve constraints on primordial non-Gaussianity is intensity
mapping [41, 42]. Clearly, intensity mapping has a lot of similarities with
the CMB, as one can imagine that each tomographic bin forms a last scattering
surface (with thickness $\Delta z$) of the photons that make up the intensity
signal. For that reason, if the bin is too wide, it is possible that
information will be lost due to artificial blurring, and this loss would be
most evident in equilateral $(N-1)$-spectra. There is in principle no
fundamental limitation on the width of a tomographic bin. But it does suggest
that one should carefully consider observational and data analysis strategies
to make sure you are not throwing away information when looking for these
signatures in the (processed) data.

The analysis in this paper was limited to only the temperature and $E$-mode
polarization signal. However, it is expected that the scaling derived from
kinematic limits will not change when $B$-mode polarization is considered in
the cosmic variance limit and assuming no lensing. If however you aim to
constrain a cosmological correlator which contains at least one tensor mode,
and the tracer for those tensor modes is a $T$-mode or $E$-mode (see e.g. [43,
44, 45, 46]), the scalar modes in those tracers would quickly overtake the
tensor contribution and the scaling should saturate on small angular scales.

Results obtained in this paper are relevant only for weak non-Gaussianity,
where we can expand the non-Gaussian field such as in e.g. Eq. (3.24). In case
the perturbations are strongly non-Gaussian or in cases where the non-
Gaussianities become manifest in the tails of the distribution, such as
proposed in Ref. [47] and explored in data in Ref. [48], the derived scalings
will not apply.

Currently the best constraints on primordial non-Gaussianity come from the CMB
bispectrum. In the near future, several CMB experiments will improve on these
constraints, predominantly131313In addition the sensitivity to polarization
will improve which will double the number of modes and, as our analysis shows,
when combined with temperature measurements can result in non-trivial
improvements on shapes with $\Delta_{c}>1/2$. by increasing the spatial
resolution, reaching to higher $\ell_{\rm max}$ [4, 2, 49]. An often quoted
threshold of any type of non-Gaussianity is to reach $\sigma_{f_{\rm NL}}\sim
1$ [50]. Given the current bounds on orthogonal and equilateral non-
Gaussianity $\sigma_{f_{\rm NL}}\sim\mathcal{O}(50)$, reaching that threshold
will be challenging if not downright impossible141414An estimate suggests that
this would require $\ell_{\rm max}>25000$ for orthogonal and $\ell_{\rm
max}>100000$ for equilateral, requiring a CMB dish of $\mathcal{O}(100)$
meters in diameter with a focal plane loaded with detectors. This assumes we
would have unconstrained access to primary modes, which we know is not the
case already for $\ell_{\rm max}\sim$ few thousand. We should note that this
extrapolation relies on the orthogonal template of Eq. (4.4). The correct
orthogonal shape should have $\Delta=2$ (see App. B of Ref. [36]) and the
expected scaling would be similar to the equilateral shape.  with CMB
measurements alone. While this is a somewhat pessimistic reading of our
analysis, the fact that the bispectrum is not the only measure of non-
Gaussianity, and that spectra with $N>3$ can exhibit very favorable scaling
with resolution, implies that the CMB can certainly contribute to the search
for primordial non-Gaussianity in general. For example, graviton exchange
trispectra are qualitatively similar to $\tau_{\rm NL}$-like non-Gaussianity,
while the corresponding bispectra are equilateral-like. Of course, this
example poses a challenging detection regardless, due to the overall
coefficient being Planck suppressed. It would be interesting to consider other
examples displaying a similar mismatch of scaling behavior of the collapsed
trispectrum vs. squeezed bispectrum.

For large scale structure, the situation is generally more optimistic if we
consider cosmic variance limits. For all $(N-1)$-spectra at least mode
counting applies and enhancement is possible in the collapsed limit. At the
same time, we are always hindered by loss of information on small scales due
to non-linearities which introduce a non-linear scale $k_{\rm NL}$ beyond
which primordial information will be hard to extract. Spectra that explicitly
couple small to large scale modes, such as the local bispectrum, introduce
observational effects in the bias of the large scale structure. These effects
will likely allow us to reach the threshold of $\sigma_{f_{\rm NL}}\sim 1$ for
these type of non-Gaussianities. In the absence of such a coupling, our
constraints will rely on measurement of $(N-1)$-spectra (with $N>2$). Ideally
a measurement of these $(N-1)$-spectra aims to avoid both limitations of
blurring and the non-linear scale, for example by mapping the density field
out to very high redshifts, such as possible with 21cm measurements [51, 52].
While many challenges lie ahead in measuring those fluctuations, it is perhaps
the only path towards measurement of $\sigma_{f_{\rm NL}}\sim
1$151515Obviously, for $f_{\rm NL}>1$, a detection will become more likely.
for shapes with $\Delta>\Delta_{c}$.

Acknowledgments We are grateful to Daniel Baumann, Anthony Challinor, Thomas
Flöss, Daniel Green, Eiichiro Komatsu, Joel Meyers, Giorgio Orlando, Enrico
Pajer, Antonio Riotto, Leonardo Senatore and Eva Silverstein for useful
discussions. We would like to thank Enrico Pajer for detailed comments on the
manuscript. A.K. and P.D.M. acknowledge support from the Netherlands
organization for scientific research (NWO) VIDI grant (dossier 639.042.730).
The work of G.P. is funded by NWO, and is part of the Delta-ITP consortium.
W.R.C. acknowledges support from the UK Science and Technology Facilities
Council (grant number ST/N000927/1).

##  Appendix A Cosmic microwave background: full-sky

In this appendix we review and summarize existing results in the literature
for the full-sky bispectrum and the trispectrum.

The CMB temperature and polarization anisotropies are expressed in terms of
the $a^{X}_{\ell m}$ spherical harmonics coefficients

$a^{X}_{\ell
m}=4\pi(-i)^{\ell}\int\frac{\mathrm{\mathrm{d}^{3}}k}{(2\pi)^{3}}\Delta_{\mathrm{\ell}}^{X}(k)\zeta(\bm{k})Y_{\ell
m}(\hat{\bm{k}}),$ (A.1)

where $\hat{\bm{n}}$ is a direction in the sky and $X=T,E$ represent the CMB
temperature and $E$-mode polarization respectively. The transfer function
$\Delta_{\mathrm{\ell}}^{X}$ encodes the linear evolution which relates
temperature and polarization anisotropies to the primordial curvature
perturbations. For our analysis, we obtain the full radiative transfer
functions using the publicly available code CAMB [33].

From the $a^{X}_{\ell m}$ coefficients, we can define the rotational-invariant
angular power spectrum

$\langle
a^{X_{1}}_{\ell_{1}m_{1}}a^{X_{2}}_{\ell_{2}m_{2}}\rangle=\delta_{\ell_{1}\ell_{2}}\delta_{m_{1}m_{2}}C^{X_{1}X_{2}}_{\ell_{1}},$
(A.2)

with

$C_{\ell}^{X_{1}X_{2}}=\frac{2}{\pi}\int\mathrm{d}k\,k^{2}P_{\zeta}(k)[\Delta_{\mathrm{\ell}}^{X}(k)]^{2}.$
(A.3)

Likewise, the three-point correlation function can be conveniently expressed
in terms of the so-called “reduced” bispectrum, which contains the physical
information about non-Gaussian sources, and a geometrical factor [34]

$\langle
a^{X_{1}}_{\ell_{1}m_{1}}a^{X_{2}}_{\ell_{2}m_{2}}a^{X_{3}}_{\ell_{3}m_{3}}\rangle=\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}b^{X_{1}X_{2}X_{3}}_{\ell_{1}\ell_{2}\ell_{3}},$
(A.4)

where the Gaunt integral
$\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}$ has a known
solution given by

$\begin{split}\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}&=\int
d\Omega_{r}Y_{\ell_{1}m_{1}}(\hat{r})Y_{\ell_{2}m_{2}}(\hat{r})Y_{\ell_{3}m_{3}}(\hat{r})\\\
&=\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell_{3}+1)}{4\pi}}\left(\begin{matrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{matrix}\right)\left(\begin{matrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
m_{1}&m_{2}&m_{3}\end{matrix}\right).\end{split}$ (A.5)

The two matrices are Wigner 3-$j$ symbols. The “reduced” bispectrum
$b^{X_{1}X_{2}X_{3}}_{\ell_{1}\ell_{2}\ell_{3}}$ is connected to the
primordial bispectrum $B_{\zeta}(k_{1},k_{2},k_{3})$ (see Eq. (1.2)) via

$\begin{split}&b^{X_{1}X_{2}X_{3}}_{\ell_{1}\ell_{2}\ell_{3}}=\left(\frac{2}{\pi}\right)^{3}\int\mathrm{d}rr^{2}\prod_{i=1}^{3}\left(\int\mathrm{d}k_{i}\,k_{i}^{2}\Delta^{X_{i}}_{\ell_{i}}(k_{i})j_{\ell_{i}}(k_{i}r)\right)B_{\zeta}(k_{1},k_{2},k_{3}),\end{split}$
(A.6)

with $j_{\ell}(kr)$ the spherical Bessel function which appear through the
Rayleigh expansion formula. In principle, given any primordial bispectrum, we
can compute the reduced bispectrum using Eq. (A.6). However, it would involve
all possible combinations of $\ell$ and $X$, which is computationally
expensive (scaling as $\ell^{5}$).

If $B_{\zeta}$ is separable, that is, it can be expressed in terms of a
product of functions that depend on a single momentum $k_{i}$, a fast and
efficient way to estimate the reduced bispectrum was introduced by the authors
of Ref. [53], the so-called Komatsu-Spergel-Wandelt (KSW) estimators. If we
define the radial functions as

$\begin{split}&\alpha_{\ell}^{X}(r)\equiv\frac{2}{\pi}\int\mathrm{d}k\,k^{2}\Delta^{X}_{\ell}(k)j_{\ell}(kr),\quad\beta_{\ell}^{X}(r)\equiv\frac{2}{\pi}\int\mathrm{d}k\,k^{2}\Delta^{X}_{\ell}(k)j_{\ell}(kr)P_{\zeta}(k),\quad\\\
&\gamma_{\ell}^{X}(r)\equiv\frac{2}{\pi}\int\mathrm{d}k\,k^{2}\Delta^{X}_{\ell}(k)j_{\ell}(kr)P^{1/3}_{\zeta}(k),\quad\delta_{\ell}^{X}(r)\equiv\frac{2}{\pi}\int\mathrm{d}k\,k^{2}\Delta^{X}_{\ell}(k)j_{\ell}(kr)P^{2/3}_{\zeta}(k),\end{split}$
(A.7)

then, using the templates introduced in Eqs. (4.2), (4.3) and (4.4), we obtain
the reduced bispectra

$\begin{split}&b^{X_{1}X_{2}X_{3},\mathrm{local}}_{\ell_{1}\ell_{2}\ell_{3}}=\frac{6}{5}f^{\mathrm{local}}_{\mathrm{NL}}\int\mathrm{d}r\,r^{2}\left[\alpha^{X_{1}}_{\ell_{1}}(r)\beta^{X_{2}}_{\ell_{2}}(r)\beta^{X_{3}}_{\ell_{3}}(r)+{\rm
2\;perms}\right],\end{split}$ (A.8)
$\begin{split}b^{X_{1}X_{2}X_{3},\mathrm{equil}}_{\ell_{1}\ell_{2}\ell_{3}}=\frac{18}{5}f^{\mathrm{equil}}_{\mathrm{NL}}&\int\mathrm{d}r\,r^{2}\bigg{[}-\alpha^{X_{1}}_{\ell_{1}}(r)\beta^{X_{2}}_{\ell_{2}}(r)\beta^{X_{3}}_{\ell_{3}}(r)-{\rm
2\;perms}\\\
&-2\delta^{X_{1}}_{\ell_{1}}(r)\delta^{X_{2}}_{\ell_{2}}(r)\delta^{X_{3}}_{\ell_{3}}(r)+\gamma^{X_{1}}_{\ell_{1}}(r)\delta^{X_{2}}_{\ell_{2}}(r)\beta^{X_{3}}_{\ell_{3}}(r)+{\rm
5\;perms}\bigg{]},\end{split}$ (A.9)
$\begin{split}b^{X_{1}X_{2}X_{3},\mathrm{ortho}}_{\ell_{1}\ell_{2}\ell_{3}}=\frac{18}{5}f^{\mathrm{ortho}}_{\mathrm{NL}}&\int\mathrm{d}r\,r^{2}\bigg{[}-3\alpha^{X_{1}}_{\ell_{1}}(r)\beta^{X_{2}}_{\ell_{2}}(r)\beta^{X_{3}}_{\ell_{3}}(r)-{\rm
2\;perms}\\\
&-8\delta^{X_{1}}_{\ell_{1}}(r)\delta^{X_{2}}_{\ell_{2}}(r)\delta^{X_{3}}_{\ell_{3}}(r)+3\gamma^{X_{1}}_{\ell_{1}}(r)\delta^{X_{2}}_{\ell_{2}}(r)\beta^{X_{3}}_{\ell_{3}}(r)+{\rm
5\;perms}\bigg{]}.\end{split}$ (A.10)

Similarly, the trispectrum is defined as the connected part of the four point
correlation function of temperature and polarization anisotropies in Eq. (A.1)

$\begin{split}\langle
a^{X_{1}}_{\ell_{1}m_{1}}a^{X_{2}}_{\ell_{2}m_{2}}a^{X_{3}}_{\ell_{3}m_{3}}a^{X_{4}}_{\ell_{4}m_{4}}\rangle_{c}&=\sum_{LM}(-1)^{M}\left(\begin{matrix}\ell_{1}&\ell_{2}&L\\\
m_{1}&m_{2}&-M\end{matrix}\right)\left(\begin{matrix}\ell_{3}&\ell_{4}&L\\\
m_{3}&m_{4}&M\end{matrix}\right)T^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L),\end{split}$
(A.11)

where we denote with $\bm{X}=X_{1}X_{2}X_{3}X_{4}$ to lighten the notation.
The trispectrum generically consists of the connected part, $T_{c}$, which
contains the non-Gaussian signatures, and the unconnected part, $T_{G}$, which
contains only the angular power spectrum and is ignored here. Using
permutation symmetry, we may write the connected part as

$\begin{split}{T_{c}}^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L)=P^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L)+(2L+1)\sum\limits_{L^{\prime}}\bigg{(}&(-1)^{\ell_{2}+\ell_{3}}\begin{Bmatrix}\ell_{1}&\ell_{2}&L\\\
\ell_{4}&\ell_{3}&L^{\prime}\end{Bmatrix}P^{\ell_{1}\ell_{3},\bm{X}}_{\ell_{2}\ell_{4}}(L^{\prime})+\\\
&+(-1)^{L+L^{\prime}}\begin{Bmatrix}\ell_{1}&\ell_{2}&L\\\
\ell_{3}&\ell_{4}&L^{\prime}\end{Bmatrix}P^{\ell_{1}\ell_{4},\bm{X}}_{\ell_{3}\ell_{2}}(L^{\prime})\bigg{)},\end{split}$
(A.12)

where each unique pairing of the multipoles implies 4 permutations

$\begin{split}P^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L)=&t^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L)+(-1)^{2L+\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}}t^{\ell_{2}\ell_{1},\bm{X}}_{\ell_{3}\ell_{4}}(L)+(-1)^{L+\ell_{3}+\ell_{4}}t^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{4}\ell_{3}}(L)+\\\
&+(-1)^{L+\ell_{1}+\ell_{2}}t^{\ell_{2}\ell_{1},\bm{X}}_{\ell_{4}\ell_{3}}(L).\end{split}$
(A.13)

Here, the matrix is the Wigner 6-$j$ symbol, which is a summation over the
product of four Wigner 3-$j$ symbols, and
$t^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L)$ is the reduced
trispectrum, defined as

$\begin{split}t^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L)=&\left(\frac{2}{\pi}\right)^{5}h_{\ell_{1}\ell_{2}L}h_{\ell_{3}\ell_{4}L}\int\mathrm{d}r_{1}\mathrm{d}r_{2}r_{1}^{2}r_{2}^{2}\int
K\,K^{2}j_{L}(Kr_{1})j_{L}(Kr_{2})\times\\\
&\times\prod_{i=1}^{4}\left(\int\mathrm{d}k_{i}\,k_{i}j_{\ell_{i}}(k_{i}r_{i})\Delta^{X_{i}}_{\ell_{i}}(k_{i})\right)T_{\zeta}(k_{1},k_{2},k_{3},k_{4}),\end{split}$
(A.14)

with

$h_{\ell_{i}\ell_{j}L}=\sqrt{\frac{(2\ell_{i}+1)(2\ell_{j}+1)(2L+1)}{4\pi}}\left(\begin{matrix}\ell_{i}&\ell_{j}&L\\\
0&0&0\end{matrix}\right).$ (A.15)

The Wigner 3-$j$ symbol guarantees that two sides of the quadrilateral and the
diagonal form a triangle. Given the definition in Eq. (A.12), we denote a
particular ordering of the multipoles as

$\displaystyle\tau^{\ell_{1}\ell_{2}\ell_{3}\ell_{4}\bm{X}}_{m_{1}m_{2}m_{3}m_{4}}=\sum_{LM}(-1)^{M}\begin{pmatrix}\ell_{1}&\ell_{2}&L\\\
m_{1}&m_{2}&-M\end{pmatrix}\begin{pmatrix}\ell_{3}&\ell_{4}&L\\\
m_{3}&m_{4}&M\end{pmatrix}t^{\ell_{1}\ell_{2},\bm{X}}_{\ell_{3}\ell_{4}}(L),$
(A.16)

and observe that the remaining combinations are given by

$\displaystyle
T^{\ell_{1}\ell_{2}\ell_{3}\ell_{4},\bm{X}}_{m_{1}m_{2}m_{3}m_{4}}=$
$\displaystyle\tau^{\ell_{1}\ell_{2}\ell_{3}\ell_{4},\bm{X}}_{m_{1}m_{2}m_{3}m_{4}}+\tau^{\ell_{2}\ell_{1}\ell_{3}\ell_{4},\bm{X}}_{m_{2}m_{1}m_{3}m_{4}}+\tau^{\ell_{1}\ell_{2}\ell_{4}\ell_{3},\bm{X}}_{m_{1}m_{2}m_{4}m_{3}}+\tau^{\ell_{2}\ell_{1}\ell_{4}\ell_{3},\bm{X}}_{m_{2}m_{1}m_{4}m_{3}}+$
$\displaystyle\tau^{\ell_{1}\ell_{3}\ell_{2}\ell_{4},\bm{X}}_{m_{1}m_{3}m_{2}m_{4}}+\tau^{\ell_{1}\ell_{3}\ell_{4}\ell_{2},\bm{X}}_{m_{1}m_{3}m_{4}m_{2}}+\tau^{\ell_{3}\ell_{1}\ell_{2}\ell_{4},\bm{X}}_{m_{3}m_{1}m_{2}m_{4}}+\tau^{\ell_{3}\ell_{1}\ell_{4}\ell_{2},\bm{X}}_{m_{3}m_{1}m_{4}m_{2}}+$
$\displaystyle\tau^{\ell_{1}\ell_{4}\ell_{3}\ell_{2},\bm{X}}_{m_{1}m_{4}m_{3}m_{2}}+\tau^{\ell_{4}\ell_{1}\ell_{3}\ell_{2},\bm{X}}_{m_{4}m_{1}m_{3}m_{2}}+\tau^{\ell_{1}\ell_{4}\ell_{2}\ell_{3},\bm{X}}_{m_{1}m_{4}m_{2}m_{3}}+\tau^{\ell_{4}\ell_{1}\ell_{2}\ell_{3},\bm{X}}_{m_{4}m_{1}m_{2}m_{3}}.$
(A.17)

We may now write the reduced trispectrum as a combination of the radial
functions defined in Eq. (A.7). Let us consider the $\tau_{\mathrm{NL}}$ type
of trispectrum of Eq. (4.7). Then

$\begin{split}t^{\tau_{\mathrm{NL}},\bm{X}}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)=\tau_{\mathrm{NL}}h_{\ell_{1}\ell_{2}L}h_{\ell_{3}\ell_{4}L}\int\mathrm{d}r_{1}\mathrm{d}r_{2}r_{1}^{2}r_{2}^{2}\,[&b_{L}(r_{1},r_{2})\beta^{X_{1}}_{\ell_{1}}(r_{1})\alpha^{X_{2}}_{\ell_{2}}(r_{1})\beta^{X_{3}}_{\ell_{3}}(r_{2})\alpha^{X_{4}}_{\ell_{4}}(r_{2})+\\\
&+\text{11 perms.}],\end{split}$ (A.18)

where we defined

$b_{L}(r_{1},r_{2})=\frac{2}{\pi}\int\mathrm{d}K\,K^{2}j_{L}(Kr_{1})j_{L}(Kr_{2})P_{\zeta}(K).$
(A.19)

While Eq. (A.19) has no exact solution for a general power spectrum, for a
scale invariant power spectrum $P_{\zeta}(K)=2\pi^{2}A_{\mathrm{s}}K^{-3}$, we
find

$\displaystyle b_{L}(r_{1},r_{2})$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}A_{\mathrm{s}}}{2}\left(\begin{array}[]{ccc}\frac{2\left(\frac{r_{2}}{r_{1}}\right)^{L}\Gamma(L)\,_{2}\tilde{F}_{1}\left(-\frac{1}{2},L;L+\frac{3}{2};\frac{r_{2}^{2}}{r_{1}^{2}}\right)}{\sqrt{\pi}}&r_{2}<r_{1}\\\
\frac{2\left(\frac{r_{1}}{r_{2}}\right)^{L}\Gamma(L)\,_{2}\tilde{F}_{1}\left(-\frac{1}{2},L;L+\frac{3}{2};\frac{r_{1}^{2}}{r_{2}^{2}}\right)}{\sqrt{\pi}}&r_{2}>r_{1}\\\
\frac{4}{\pi}\frac{1}{L(1+L)}&r_{2}=r_{1}\\\ \end{array}\right).$ (A.23)

Here $\tilde{F}_{1}$ is a hypergeometric function. The solution above is
numerically unstable, because for large arguments its value relies on
cancellation of large numbers. However, further investigation shows that for
$L\leq 10$, $b_{L}(r_{1},r_{2})$ is approximately constant over the last
scattering surface. Since this would considerably reduce the computational
cost of $t^{\tau_{\mathrm{NL}},\bm{X}}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)$,
we take this approximation when computing the $S/N$ for this trispectrum.
Therefore, the $\tau_{\mathrm{NL}}$ angular reduced trispectrum is given by

$\begin{split}t^{\tau_{\mathrm{NL}},\bm{X}}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)=\frac{2\pi
A_{\mathrm{s}}\tau_{\mathrm{NL}}}{L(1+L)}h_{\ell_{1}\ell_{2}L}h_{\ell_{3}\ell_{4}L}\int\mathrm{d}r_{1}\mathrm{d}r_{2}r_{1}^{2}r_{2}^{2}\,[&\beta^{X_{1}}_{\ell_{1}}(r_{1})\alpha^{X_{2}}_{\ell_{2}}(r_{1})\beta^{X_{3}}_{\ell_{3}}(r_{2})\alpha^{X_{4}}_{\ell_{4}}(r_{2})+\\\
&\text{11 perms.}].\end{split}$ (A.24)

Along these lines, the $g^{\mathrm{local}}_{\mathrm{NL}}$ reduced trispectrum
is given by

$\begin{split}t^{g_{\mathrm{NL}},\bm{X}}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)=\frac{54}{25}g^{\mathrm{local}}_{\mathrm{NL}}h_{\ell_{1}\ell_{2}L}h_{\ell_{3}\ell_{4}L}\int&\mathrm{d}r_{1}\mathrm{d}r_{2}r_{1}^{2}r_{2}^{2}\,[a_{L}(r_{1},r_{2})\alpha^{X_{1}}_{\ell_{1}}(r_{1})\beta^{X_{2}}_{\ell_{2}}(r_{1})\beta^{X_{3}}_{\ell_{3}}(r_{2})\beta^{X_{4}}_{\ell_{4}}(r_{2})+\\\
&a_{L}(r_{1},r_{2})\beta^{X_{1}}_{\ell_{1}}(r_{1})\beta^{X_{2}}_{\ell_{2}}(r_{1})\beta^{X_{3}}_{\ell_{3}}(r_{2})\alpha^{X_{4}}_{\ell_{4}}(r_{2})+\text{2
perms.}],\end{split}$ (A.25)

with

$a_{L}(r_{1},r_{2})=\frac{2}{\pi}\int\mathrm{d}K\,K^{2}j_{L}(Kr_{1})j_{L}(Kr_{2})=\frac{1}{r_{1}^{2}}\delta(r_{1}-r_{2}),$
(A.26)

Replacing Eq. (A.26) in
$t^{g_{\mathrm{NL}},\bm{X}}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)$, then we
are left with

$\begin{split}t^{g_{\mathrm{NL}},\bm{X}}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)=\frac{54}{25}g^{\mathrm{local}}_{\mathrm{NL}}h_{\ell_{1}\ell_{2}L}h_{\ell_{3}\ell_{4}L}\int&\mathrm{d}r_{1}r_{1}^{2}\,[\alpha^{X_{1}}_{\ell_{1}}(r_{1})\beta^{X_{2}}_{\ell_{2}}(r_{1})\beta^{X_{3}}_{\ell_{3}}(r_{2})\beta^{X_{4}}_{\ell_{4}}(r_{2})+\\\
&\beta^{X_{1}}_{\ell_{1}}(r_{1})\beta^{X_{2}}_{\ell_{2}}(r_{1})\beta^{X_{3}}_{\ell_{3}}(r_{2})\alpha^{X_{4}}_{\ell_{4}}(r_{2})+\text{2
perms.}].\end{split}$ (A.27)

##  Appendix B Cosmic microwave background: flat-sky

In the flat-sky approximation we ignore the curvature of the sky, which is
equivalent to approximate the sphere in the neighborhood of a point by the
tangent plane at that point. As a consequence, the spherical harmonic
expansion of a perturbation is reduced to a simple Fourier transform [54, 55].
In this way, it will be easier to appreciate the effect of diffusion damping
on the CMB statistics. For simplicity, we consider only temperature
anisotropies $\Delta T/T(\hat{\bm{n}})$.

In the line-of-sight approach [56], these are given by

$\begin{split}\frac{\Delta
T}{T}(\hat{\bm{n}})&=\int\frac{\mathrm{\mathrm{d}^{3}}k}{(2\pi)^{3}}\,\zeta(\bm{k})\int_{0}^{\tau_{0}}\mathrm{d}\tau\mathrm{e}^{i\bm{k}\cdot\hat{\bm{n}}(\tau-\tau_{0})}S(k,\tau),\end{split}$
(B.1)

where $S(k,\tau)$ is the CMB source function, which encodes all the
information about metric perturbations and photon fluctuations, and $\tau$ is
the conformal time and $\tau_{0}$ refers to the conformal time today. The
projection on a plane perpendicular to the line of sight is given by [57]

$\begin{split}\frac{\Delta
T}{T}(\hat{\bm{n}}^{\perp})&=\int\mathrm{d}r\,F(r)\frac{\Delta
T}{T}(\hat{\bm{n}}^{\perp},r)\\\
&=\int\mathrm{d}r\,F(r)\int\frac{\mathrm{\mathrm{d}^{3}}k}{(2\pi)^{3}}\,\zeta(\bm{k})\int_{0}^{\tau_{0}}\mathrm{d}\tau\mathrm{e}^{i\bm{k}\cdot\hat{\bm{n}}(\tau-\tau_{0})}S(k,\tau),\end{split}$
(B.2)

with $F(r)$ the projection kernel. Here we consider
$F(r)=\delta(r-r_{\mathrm{rec}})$, where $r_{\mathrm{rec}}$ is the distance
from the last scattering surface. Accordingly, we define the flat-sky
temperature anisotropies as

$\begin{split}a(\bm{\ell})&=\int\mathrm{\mathrm{d}^{2}}x^{\perp}\frac{\Delta
T}{T}(\hat{\bm{n}}^{\perp})\mathrm{e}^{-i\bm{\ell}\cdot\hat{\bm{n}}^{\perp}}\\\
&=\int\mathrm{\mathrm{d}^{2}}\hat{n}^{\perp}\int\frac{\mathrm{\mathrm{d}^{3}}k}{(2\pi)^{3}}\,\zeta(\bm{k})\int_{0}^{\tau_{0}}\mathrm{d}\tau\,\mathrm{e}^{i[\bm{k}^{\perp}(\tau-\tau_{0})-\bm{\ell}]\cdot\hat{\bm{n}}^{\perp}}\mathrm{e}^{ik^{\parallel}(\tau-\tau_{0})}S(k,\tau),\end{split}$
(B.3)

where we used the decomposition $\bm{k}=(\bm{k}^{\perp},k^{\parallel})$ such
that $k=\sqrt{(k^{\parallel})^{2}+\bm{k}^{\perp}\cdot\bm{k}^{\perp}}$, and
$\bm{\ell}$ is a $2$D wavevector. Finally, since

$\int\mathrm{\mathrm{d}^{2}}\hat{n}^{\perp}\,\mathrm{e}^{i(\bm{k}^{\perp}(\tau_{1}-\tau_{0})-\bm{\ell})\cdot\hat{\bm{n}}^{\perp}}=(2\pi)^{2}\delta^{(2)}(\bm{k}^{\perp}(\tau-\tau_{0})-\bm{\ell}),$
(B.4)

the temperature anisotropies can be rewritten as

$\begin{split}a(\bm{\ell})=\int\frac{\mathrm{\mathrm{d}^{3}}k}{(2\pi)^{3}}\,\zeta(\bm{k})\Delta_{\ell}^{T}(k^{\parallel}),\end{split}$
(B.5)

with the flat-sky transfer function $\Delta_{\ell}^{T}(k^{\parallel})$ given
by

$\Delta_{\ell}^{T}(k^{\parallel})=(2\pi)^{2}\int_{0}^{\tau_{0}}\mathrm{d}\tau\,\delta^{(2)}(\bm{k}^{\perp}(\tau-\tau_{0})-\bm{\ell})\mathrm{e}^{ik^{\parallel}(\tau-\tau_{0})}S(k,\tau).$
(B.6)

At this point, we can introduce the flat-sky power spectrum

$\langle
a(\bm{\ell}_{1})a(\bm{\ell}_{2})\rangle=(2\pi)^{2}\delta^{(2)}(\bm{\ell}_{1}+\bm{\ell}_{2})C(\ell_{1}).$
(B.7)

In the simplest scenario, where the radiation transfer effect is negligible,
the temperature anisotropies read as

$\begin{split}a(\bm{\ell})=(2\pi)^{2}\int\frac{\mathrm{\mathrm{d}^{3}}k}{(2\pi)^{3}}\,\zeta(\bm{k})\delta^{(2)}(\bm{k}^{\perp}r_{\mathrm{rec}}-\bm{\ell})\mathrm{e}^{ik^{\parallel}r_{\mathrm{rec}}}.\end{split}$
(B.8)

With a few calculations, we obtain $C(\ell)\propto\ell^{-2}$. Given that the
scope of this work focuses on the scaling relations, we omitted all the
coefficients. In a similar fashion, we derive the flat-sky bispectrum

$\langle
a(\bm{\ell}_{1})a(\bm{\ell}_{2})a(\bm{\ell}_{3})\rangle=(2\pi)^{2}\delta^{(2)}(\bm{\ell}_{1}+\bm{\ell}_{2}+\bm{\ell}_{3})B(\ell_{1},\ell_{2},\ell_{3}),$
(B.9)

where

$B(\ell_{1},\ell_{2},\ell_{3})=\prod_{i=1}^{3}\left(\int\frac{\mathrm{\mathrm{d}^{3}}k_{i}}{(2\pi)^{3}}\,\Delta^{T}_{\ell_{i}}(k^{\parallel}_{i})\right)\,\delta^{(3)}(\bm{k}_{123})B_{\zeta}(k_{1},k_{2},k_{3}).$
(B.10)

As an example, we consider here the local bispectrum of Eq. (4.2). Replacing
Eq. (B.8), we obtain

$\begin{split}B(\ell_{1},\ell_{2},\ell_{3})\approx
A_{\mathrm{s}}^{2}f_{\mathrm{NL}}\left[\frac{1}{\ell_{1}^{2}\ell_{2}^{2}}+\text{2
perms.}\right].\end{split}$ (B.11)

Generally, we may define an angular flat-sky $(N-1)$-spectrum as

$\langle a(\bm{\ell}_{1})\cdots
a(\bm{\ell}_{N})\rangle=(2\pi)^{2}\delta^{(2)}(\bm{\ell}_{1\dots
N})F(\ell_{1},\dots,\ell_{N}).$ (B.12)

The scale invariance of the $(N-1)$-spectrum implies that
$F\sim\ell^{-2(N-1)}$.

##  Appendix C The effect of damping on CMB statistics

Let us now include the effect of diffusion damping at small scales by
introducing [18, 15, 58]

$\Delta_{\ell}^{T}(k^{\parallel})\approx\delta^{(2)}(\bm{k}^{\perp}(\tau_{\mathrm{rec}}-\tau_{0})-\bm{\ell})\mathrm{e}^{-1/2(\ell/\ell_{\mathrm{D}})^{1.2}}\mathrm{e}^{-1/2(k^{\parallel}/k_{\mathrm{D}})^{1.2}},$
(C.1)

where $\ell_{\mathrm{D}}\simeq 1300$ is the damping scale, namely the scale
above which diffusion damping dominates over gravitational pull, and the
exponential cuts the integral off at the corresponding damping mode
$k_{\mathrm{D}}$. The exponent $1.2$ comes from the study of Ref. [59] and is
approximate. Its precise value is irrelevant for the work presented in this
paper. This transfer function highlights the fact that the damping effect is
3-dimensional: it has a $k^{\perp}\propto\ell$ and a line of sight
$k^{\parallel}$ component. The former integrates out, effectively becoming a
multiplicative transfer function term. However, the line of sight damping,
which probes the thickness of the last scattering surface, remains a
convolution, coupling $k^{\parallel}$ modes. The net result is a change in the
scaling behavior of the correlators.

Indeed, using Eq. (C.1) in the power spectrum expression, we obtain

$\begin{split}C(\ell)\approx\frac{A_{\mathrm{s}}}{\ell^{2}}\frac{\mathrm{e}^{-(\ell/\ell_{\mathrm{D}})^{1.2}}}{\sqrt{1+(\ell/\ell_{\mathrm{D}})^{2}}}\approx\frac{A_{\mathrm{s}}}{\ell^{3}}\mathrm{e}^{-(\ell/\ell_{\mathrm{D}})^{1.2}}\end{split}$
(C.2)

where the last passage holds for $\ell\gg\ell_{\mathrm{D}}$. As expected, the
angular power spectrum changes scaling, from $\ell^{-2}$ to $\ell^{-3}$.
Likewise, for the local bispectrum it is straightforward to show that

$\begin{split}B(\ell_{1},\ell_{2},\ell_{3})\approx
A_{\mathrm{s}}^{2}f_{\mathrm{NL}}\mathrm{e}^{-(\ell_{1}^{1.2}+\ell_{2}^{1.2}+\ell_{3}^{1.2})/\ell_{\mathrm{D}}^{1.2}}\left[\frac{1}{\ell_{1}^{2}\sqrt{1+(\ell_{1}/\ell_{\mathrm{D}})^{2}}}\frac{1}{\ell_{2}^{2}\sqrt{1+(\ell_{2}/\ell_{\mathrm{D}})^{2}}}+\text{2
perms.}\right]\end{split}$ (C.3)

and for $\ell\gg\ell_{\mathrm{D}}$ we obtain

$\begin{split}B(\ell_{1},\ell_{2},\ell_{3})\approx
A_{\mathrm{s}}^{2}f_{\mathrm{NL}}\mathrm{e}^{-(\ell_{1}^{1.2}+\ell_{2}^{1.2}+\ell_{3}^{1.2})/\ell_{\mathrm{D}}^{1.2}}\left[\frac{1}{\ell_{1}^{3}\ell_{2}^{3}}+\text{2
perms.}\right]\end{split}$ (C.4)

In general, the expected scaling can be obtained by a simple replacement
$\ell^{2}\rightarrow\ell^{2}(\ell/\ell_{D})\mathrm{e}^{-(\ell/\ell_{\mathrm{D}})^{1.2}}$,
but the exponential factor appears equally in the numerator and the
denominator in the $S/N$ computation, therefore it cancels out.

##  Appendix D Optimal estimator for non-Gaussianity and signal-to-noise

###  D.1 Large scale structure

Optimal estimators for the amplitude of non-Gaussianity can be constructed
from the bispectrum and trispectrum of the primordial perturbations, and
extended to the general $(N-1)$-spectrum (1.2).

In the weak non-Gaussian limit and assuming statistical isotropy and
homogeneity, it has been shown that the optimal estimator coming from the
bispectrum is given by [60, 61]

$\begin{split}\widehat{\mathcal{E}}^{(3)}=\frac{1}{\mathcal{N}_{(3)}}\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{3}}{(2\pi)^{3}}\,\frac{\langle\zeta^{\mathrm{th}}_{\mathbf{k}_{1}}\zeta^{\mathrm{th}}_{\mathbf{k}_{2}}\zeta^{\mathrm{th}}_{\mathbf{k}_{3}}\rangle}{P(k_{1})P(k_{2})P(k_{3})}[\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\zeta_{\mathbf{k}_{3}}-3\langle\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\rangle\zeta_{\mathbf{k}_{3}}].\end{split}$
(D.1)

Here the primordial curvature perturbation $\zeta_{\mathbf{k}}$ comes from
observations, while $\zeta^{\mathrm{th}}_{\mathbf{k}}$ is the theoretically
motivated one. The linear term
$\langle\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\rangle\zeta_{\mathbf{k}_{3}}$
accounts for any anisotropic effect due to systematics in the dataset. This
terms is irrelevant for the simple Fisher forecasts we present here, and we
will neglect this term from hereon. For the estimator to be unbiased, the
normalization factor $\mathcal{N}$ must be given by the inverse Fisher
information

$\begin{split}\mathcal{N}=F^{-1}\equiv
V(2\pi)^{3}\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{3}}{(2\pi)^{3}}\,\delta^{(3)}(\bm{k}_{123})\frac{B^{2}_{\zeta}(k_{1},k_{2},k_{3})}{P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(k_{3})},\end{split}$
(D.2)

where $P(k)$ and $B_{\zeta}(k_{1},k_{2},k_{3})$ are the power spectrum and the
bispectrum respectively (Eq. (1.2)), and $V$ is the volume of the survey.

Since an optimal estimator also saturates the Cramèr-Rao inequality, its
minimum error is given by the Fisher information itself, thus the signal-to-
noise ratio

$\sigma(\widehat{\mathcal{E}})=F^{-1/2}=(S/N)^{-1}.$ (D.3)

From now on we will use only the $S/N$,

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(3)}=V(2\pi)^{3}\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{3}}{(2\pi)^{3}}\,\delta^{(3)}(\bm{k}_{123})\frac{B^{2}_{\zeta}(k_{1},k_{2},k_{3})}{P(k_{1})P(k_{2})P(k_{3})},\end{split}$
(D.4)

with its connection to the variance of the estimator explicit in Eq. (D.3).

Likewise, we define the trispectrum optimal estimator

$\begin{split}\widehat{\mathcal{E}}^{(4)}=\frac{1}{\mathcal{N}_{(4)}}\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{3}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{4}}{(2\pi)^{3}}\,&\frac{\langle\zeta^{\mathrm{th}}_{\mathbf{k}_{1}}\zeta^{\mathrm{th}}_{\mathbf{k}_{2}}\zeta^{\mathrm{th}}_{\mathbf{k}_{3}}\zeta^{\mathrm{th}}_{\mathbf{k}_{4}}\rangle_{c}}{P(k_{1})P(k_{2})P(k_{3})P(k_{4})}[\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\zeta_{\mathbf{k}_{3}}\zeta_{\mathbf{k}_{4}})-\\\
&-6\langle\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\rangle\zeta_{\mathbf{k}_{3}}\zeta_{\mathbf{k}_{4}}+3\langle\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\rangle\langle\zeta_{\mathbf{k}_{3}}\zeta_{\mathbf{k}_{4}})],\end{split}$
(D.5)

and the $(N-1)$-spectrum estimator

$\begin{split}\widehat{\mathcal{E}}^{(N)}=\frac{1}{\mathcal{N}_{(N)}}\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}&\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\cdots\frac{\mathrm{d}^{3}k_{N}}{(2\pi)^{3}}\,\frac{\langle\zeta^{\mathrm{th}}_{\mathbf{k}_{1}}\zeta^{\mathrm{th}}_{\mathbf{k}_{2}}\cdots\zeta^{\mathrm{th}}_{\mathbf{k}_{N}}\rangle_{c}}{P(k_{1})P(k_{2})\cdots
P(k_{N})}[\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\cdots\zeta_{\mathbf{k}_{N}}-\\\
&-A\langle\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\rangle\zeta_{\mathbf{k}_{3}}\cdots\zeta_{\mathbf{k}_{N}}+B\langle\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\rangle\langle\zeta_{\mathbf{k}_{3}}\zeta_{\mathbf{k}_{4}}\rangle\zeta_{\mathbf{k}_{5}}\cdots\zeta_{\mathbf{k}_{N}}+\dots],\end{split}$
(D.6)

where $A$ and $B$ are coefficient that take into account permutations of
terms, and the dots refer to products of two-point correlation functions that
arise as $N$ increases. These terms are necessary to generalize to the case of
incomplete sample coverage, inhomogeneous noise and removing any Gaussian
noise bias . However for the purpose of this a simple Fisher forecast, we can
neglect them.

Using Eq. (D.3), the $(S/N)^{2}$ of the trispectrum is given by

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(4)}=V^{(4)}(2\pi)^{3}\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{3}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{4}}{(2\pi)^{3}}\,\delta^{(3)}(\bm{k}_{1234})\frac{T^{2}_{\zeta}(k_{1},k_{2},k_{3},k_{4})}{P(k_{1})P(k_{2})P(k_{3})P(k_{4})},\end{split}$
(D.7)

and in the case of $(N-1)$-spectrum we obtain

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(N)}=V^{(N)}(2\pi)^{3}\int\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}}\cdots\frac{\mathrm{d}^{3}k_{N}}{(2\pi)^{3}}\,\delta^{(3)}(\bm{k}_{12\dots
N})\frac{F^{2}_{\zeta}(k_{1},k_{2},\dots,k_{N})}{P(k_{1})P(k_{2})\cdots
P(k_{N})},\end{split}$ (D.8)

###  D.2 Cosmic microwave background: full-sky analysis

Similarly to the LSS case, we can define an optimal estimator to extract non-
Gaussian information from CMB data. Let us consider the bispectrum optimal
estimator, given by [35]

$\begin{split}\widehat{\mathcal{E}}^{(3)}=\frac{1}{F^{-1}}&\sum_{X_{i},X^{\prime}_{i}}\sum_{\ell_{i},m_{i}}\sum_{\ell^{\prime}_{i},m^{\prime}_{i}}\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}b^{X_{1}X_{2}X_{3},\mathrm{th}}_{\ell_{1}\ell_{2}\ell_{3}}\times\\\
\times&\bigg{\\{}\left[\left(\mathbb{C}^{-1}_{\ell_{1}m_{1},\ell_{1}^{\prime}m_{1}^{\prime}}\right)^{X_{1}X^{\prime}_{1}}a^{X^{\prime}_{1}}_{\ell_{1}^{\prime}m_{1}^{\prime}}\left(\mathbb{C}^{-1}_{\ell_{2}m_{2},\ell_{2}^{\prime}m_{2}^{\prime}}\right)^{X_{2}X^{\prime}_{2}}a^{X^{\prime}_{2}}_{\ell_{2}^{\prime}m_{2}^{\prime}}\left(\mathbb{C}^{-1}_{\ell_{3}m_{3},\ell_{3}^{\prime}m_{3}^{\prime}}\right)^{X_{3}X^{\prime}_{3}}a^{X^{\prime}_{3}}_{\ell_{3}^{\prime}m_{3}^{\prime}}\right]-\\\
&-\left[\left(\mathbb{C}^{-1}_{\ell_{1}m_{1},\ell_{2}m_{2}}\right)^{X_{1}X_{2}}\left(\mathbb{C}^{-1}_{\ell_{3}m_{3},\ell_{3}^{\prime}m_{3}^{\prime}}\right)^{X_{3}X^{\prime}_{3}}a^{X^{\prime}_{3}}_{\ell_{3}^{\prime}m_{3}^{\prime}}+\mathrm{cyclic}\right]\bigg{\\}},\end{split}$
(D.9)

where $b_{\ell_{1}\ell_{2}\ell_{3}}^{\rm th}$ is the theoretical reduced
bispectrum and $a^{X_{1}}_{\ell_{1}m_{1}}$ are observed/simulated multipoles
and $F^{-1}$ the inverse Fisher information. $\mathbb{C}^{-1}$ is the inverse
of the covariance matrix, which is given by a block matrix

$\mathbb{C}=\left(\begin{array}[]{cc}\mathbb{C}^{TT}&\mathbb{C}^{TE}\\\
\mathbb{C}^{ET}&\mathbb{C}^{EE}\end{array}\right),$ (D.10)

where the blocks represent the full $TT$, $TE$, and $EE$ covariance matrices,
with $\mathbb{C}^{ET}$ being the transpose of $\mathbb{C}^{TE}$. We neglect
the linear term in the estimator (D.9) since it is proportional to the
monopole when rotational invariance is valid. We also work in the “diagonal
covariance” approximation [62], which reduces the computation of the inverse
covariance matrix to a $2\times 2$ matrix

$\mathbb{C}_{\ell}=\begin{pmatrix}C^{TT}_{\ell}&&C^{TE}_{\ell}\\\
C^{ET}_{\ell}&&C^{EE}_{\ell}\end{pmatrix},$ (D.11)

and $C_{\ell}$ is the CMB power spectrum. Thus, the estimator (D.9) reads as

$\begin{split}\widehat{\mathcal{E}}^{(3)}=&\frac{1}{F^{-1}}\sum_{X_{i},X^{\prime}_{i}}\sum_{\ell_{i},m_{i}}\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}b^{X_{1}X_{2}X_{3},\mathrm{th}}_{\ell_{1}\ell_{2}\ell_{3}}(\mathbb{C}^{-1})_{\ell_{1}}^{X_{1}X^{\prime}_{1}}(\mathbb{C}^{-1})_{\ell_{2}}^{X_{2}X^{\prime}_{2}}(\mathbb{C}^{-1})_{\ell_{3}}^{X_{3}X^{\prime}_{3}}a^{X^{\prime}_{1}}_{\ell_{1}m_{1}}a^{X^{\prime}_{2}}_{\ell_{2}m_{2}}a^{X^{\prime}_{3}}_{\ell_{3}m_{3}}.\end{split}$
(D.12)

Since the estimator (D.9) is optimal, the variance on it is given by
$\sigma(\widehat{\mathcal{E}}^{(3)})=(F^{-1})^{1/2}=(S/N)^{-1}$, with the
$(S/N)^{2}$ given by

$\left(\frac{S}{N}\right)^{2}_{(3)}=\sum_{X_{i},X^{\prime}_{i}}\sum_{\ell_{i}}h^{\ell_{1}\ell_{2}\ell_{3}}b^{X_{1}X_{2}X_{3}}_{\ell_{1}\ell_{2}\ell_{3}}(\mathbb{C}^{-1})_{\ell_{1}}^{X_{1}X^{\prime}_{1}}(\mathbb{C}^{-1})_{\ell_{2}}^{X_{2}X^{\prime}_{2}}(\mathbb{C}^{-1})_{\ell_{3}}^{X_{3}X^{\prime}_{3}}b^{X^{\prime}_{1}X^{\prime}_{2}X^{\prime}_{3}}_{\ell_{1}\ell_{2}\ell_{3}},$
(D.13)

where

$h^{\ell_{1}\ell_{2}\ell_{3}}=\sum_{m_{i}}(\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}})^{2}=\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell_{3}+1)}{4\pi}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
0&0&0\end{pmatrix}^{2}.$ (D.14)

Here, the Wigner-$3j$ matrix in $h^{\ell_{1}\ell_{2}\ell_{3}}$ enforces the
triangle condition between the triplet $(\ell_{1},\ell_{2},\ell_{3})$. In that
sense, the domain of calculation is a tetrahedron. However, since the problem
is symmetric under permutation of indices, we can assume
$\ell_{1}\leq\ell_{2}\leq\ell_{3}$ without any loss in generality and with a
significant gain in computational speed. Therefore, the $(S/N)^{2}$ becomes

$\left(\frac{S}{N}\right)^{2}_{(3)}=\sum_{X_{i},X^{\prime}_{i}}\sum_{\ell_{1}\leq\ell_{2}\leq\ell_{3}}\frac{h^{\ell_{1}\ell_{2}\ell_{3}}}{f(\ell_{1},\ell_{2},\ell_{3})}b^{X_{1}X_{2}X_{3}}_{\ell_{1}\ell_{2}\ell_{3}}(C_{\ell_{1}}^{X_{1}X^{\prime}_{1}})^{-1}(C_{\ell_{2}}^{X_{2}X^{\prime}_{2}})^{-1}(C_{\ell_{3}}^{X_{3}X^{\prime}_{3}})^{-1}b^{X^{\prime}_{1}X^{\prime}_{2}X^{\prime}_{3}}_{\ell_{1}\ell_{2}\ell_{3}}.$
(D.15)

The function $f(\ell_{1},\ell_{2},\ell_{3})$ takes values 1, 2, and 6 when all
$\ell$’s are different, two of them are same, and all are same, respectively.

Likewise, we define an optimal estimator for the trispectrum as [38]

$\displaystyle\widehat{\mathcal{E}}^{(4)}$ $\displaystyle=$
$\displaystyle\frac{1}{F^{-1}}\sum_{X_{i},X^{\prime}_{i}}\sum_{\ell_{i},m_{i}}\langle
a_{\ell_{1}m_{1}}^{X_{1}}a_{\ell_{2}m_{2}}^{X_{2}}a_{\ell_{3}m_{3}}^{X_{3}}a_{\ell_{4}m_{4}}^{X_{4}}\rangle_{c}\bigg{[}\left(\mathbb{C}^{-1}\right)_{\ell_{1}}^{X_{1}X^{\prime}_{1}}\left(\mathbb{C}^{-1}\right)_{\ell_{2}}^{X_{2}X^{\prime}_{2}}\left(\mathbb{C}^{-1}\right)_{\ell_{3}}^{X_{3}X^{\prime}_{3}}\left(\mathbb{C}^{-1}\right)_{\ell_{4}}^{X_{4}X^{\prime}_{4}}\times$
(D.16) $\displaystyle\times
a_{\ell_{1}m_{1}}^{X^{\prime}_{1}}a_{\ell_{2}m_{2}}^{X^{\prime}_{2}}a_{\ell_{3}m_{3}}^{X^{\prime}_{3}}a_{\ell_{4}m_{4}}^{X^{\prime}_{4}}\bigg{]},$

where again we drop the linear terms.

As for the bispectrum estimator, in order to have an unbiased estimator, the
normalization must be the inverse Fisher matrix. Then the estimator variance,
the $(S/N)^{2}$, reads

$\displaystyle\left(\frac{S}{N}\right)^{2}_{(4)}$ $\displaystyle=$
$\displaystyle\sum_{X_{i},X^{\prime}_{i}}\sum\limits_{\ell_{i}}\sum\limits_{L}\frac{1}{(2L+1)}T^{X_{1}X_{2}X_{3}X_{4}}_{c,\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L)\left(\mathbb{C}^{-1}\right)_{\ell_{1}}^{X_{1}X^{\prime}_{1}}\left(\mathbb{C}^{-1}\right)_{\ell_{2}}^{X_{2}X^{\prime}_{2}}\times$
(D.17)
$\displaystyle\left(\mathbb{C}^{-1}\right)_{\ell_{3}}^{X_{3}X^{\prime}_{3}}\left(\mathbb{C}^{-1}\right)_{\ell_{4}}^{X_{4}X^{\prime}_{4}}T^{X^{\prime}_{1}X^{\prime}_{2}X^{\prime}_{3}X^{\prime}_{4}}_{c,\ell_{1}\ell_{2}\ell_{3}\ell_{4}}(L).$

The $L$ term at the denominator comes from the summation of Wigner-$3j$
matrices

$\begin{split}&\sum_{m_{1}m_{2}}\begin{pmatrix}\ell_{1}&\ell_{2}&L\\\
m_{1}&m_{2}&-M\end{pmatrix}\begin{pmatrix}\ell_{1}&\ell_{2}&L^{\prime}\\\
m_{1}&m_{2}&-M^{\prime}\end{pmatrix}=\frac{\delta_{LL^{\prime}}\delta_{MM^{\prime}}}{2L+1},\\\
&\sum_{M}\sum_{m_{3}m_{4}}\begin{pmatrix}\ell_{3}&\ell_{4}&L\\\
m_{3}&m_{4}&M\end{pmatrix}\begin{pmatrix}\ell_{3}&\ell_{4}&L\\\
m_{3}&m_{4}&M\end{pmatrix}=1.\end{split}$ (D.18)

###  D.3 Cosmic microwave background: flat-sky analysis

The $S/N$ in the flat-sky approximation can be easily derived from the full-
sky definition (D.15) (for temperature only) by replacing

$\sum_{\ell_{i}}\rightarrow\int\mathrm{d^{2}}\ell_{i},\quad
h^{\ell_{1}\ell_{2}\ell_{3}}b_{\ell_{1}\ell_{2}\ell_{3}}\rightarrow\delta^{(2)}(\bm{\ell}_{123})B(\ell_{1},\ell_{2},\ell_{3}),\quad
C^{-1}_{\ell_{i}}\rightarrow C^{-1}(\ell_{i}).$ (D.19)

Therefore, for the bispectrum we have

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(3)}=\frac{f_{\mathrm{sky}}}{6\pi}\frac{1}{(2\pi)^{2}}\int\mathrm{d}^{2}\ell_{1}\mathrm{d}^{2}\ell_{2}\mathrm{d}^{2}\ell_{3}\delta^{(2)}(\bm{\ell}_{123})\frac{B^{2}(\ell_{1},\ell_{2},\ell_{3})}{C(\ell_{1})C(\ell_{2})C(\ell_{3})},\end{split}$
(D.20)

where $f_{\mathrm{sky}}$ refers to the portion of the observed sky. Along
these lines we can define the trispectrum flat-sky $S/N$ as

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(4)}=\frac{f_{\mathrm{sky}}}{6\pi}\frac{1}{(2\pi)^{2}}\int\mathrm{d}^{2}\ell_{1}\mathrm{d}^{2}\ell_{2}\mathrm{d}^{2}\ell_{3}\mathrm{d}^{2}\ell_{4}\delta^{(2)}(\bm{\ell}_{1234})\frac{T^{2}(\ell_{1},\ell_{2},\ell_{3},\ell_{4})}{C(\ell_{1})C(\ell_{2})C(\ell_{3})C(\ell_{4})},\end{split}$
(D.21)

and the $S/N$ from the $N$-point correlation function

$\begin{split}\left(\frac{S}{N}\right)^{2}_{(N)}=\frac{f_{\mathrm{sky}}}{6\pi}\frac{1}{(2\pi)^{2}}\int\mathrm{d}^{2}\ell_{1}\mathrm{d}^{2}\ell_{2}\cdots\mathrm{d}^{2}\ell_{N}\delta^{(2)}(\bm{\ell}_{12\dots
N})\frac{F^{2}(\ell_{1},\ell_{2},\dots,\ell_{N})}{C(\ell_{1})C(\ell_{2})\cdots
C(\ell_{N})}.\end{split}$ (D.22)

##  Appendix E Multiple squeezed and collapsed limits

In this appendix, we derive the effect of multiple squeezed or collapsed
limits on the $S/N$. In our numerical analysis we encountered the
$g_{\mathrm{NL}}$-trispectrum, where the scaling of the $S/N$ is determined by
nested squeezed limits which results in a $S/N$ that scales as mode-counting.
We will show that nested squeezed limits never exceed mode counting (besides
the $\log$-enhancement), whereas, as already seen in Sec. 3.3 for a single
collapsed limit, multiple collapsed limits can go beyond mode counting. Our
calculations are highly idealistic and simplifying, therefore the enhanced
scalings should be carefully interpreted within a specific model.

First, we consider the double squeezed limit of the damped
trispectrum161616The higher number of squeezed limits does not affect the
scaling of $3$D and undamped $2$D spectra.

$\lim_{\ell_{1}\ll\ell_{2}\ll\ell_{3,4}}T(\ell_{1},\ell_{2},\ell_{3},\ell_{4})\propto\frac{\ell_{\mathrm{D}}^{3}}{(\ell_{1}\ell_{2}\ell_{3})^{3}}\left(\frac{\ell_{1}\ell_{2}}{\ell_{3}^{2}}\right)^{\Delta},$
(E.1)

then the leading scaling in $\ell_{\mathrm{max}}$ is given by

$\begin{split}\left(\frac{S}{N}\right)_{(4)}^{2}\sim\ell_{\mathrm{D}}^{2}\int_{\ell_{\mathrm{D}}}^{c\ell_{\mathrm{max}}}d^{2}\ell_{1}d^{2}\ell_{2}\int_{c\ell_{\mathrm{max}}}^{\ell_{\mathrm{max}}}d^{2}\ell_{3}\frac{\ell_{1}^{-6+2\Delta}\ell_{2}^{-6+2\Delta}\ell_{3}^{-6-4\Delta}}{\ell_{1}^{-3}\ell_{2}^{-3}\ell_{3}^{-6}}\sim\ell_{\mathrm{max}}^{2-4\Delta},\quad\Delta>0.\end{split}$
(E.2)

As shown in Sec. 3.2, for local shapes, i.e. $\Delta=0$, the signal-to-noise
scales as $\ell_{\mathrm{max}}^{2}$. In general, there is an improvement with
respect to the expected scaling $\log\ell_{\mathrm{max}}$ for values
$0\leq\Delta<1/2$.

It is possible to extend the argument to $S$ squeezed limits,
$\ell_{1}\ll\ell_{2}\ll\dots\ll\ell_{S}\ll\ell_{S+1}\sim\dots\sim\ell_{N}$.
Notice that the maximum number of squeezed limits is $S_{\mathrm{max}}=N-2$.
In this limit, the $(N-1)$-spectrum reads

$F_{\zeta}(\ell_{1},\dots,\ell_{N})\sim\frac{\ell_{\mathrm{D}}^{N-1}}{\ell_{1}^{3S}\ell^{3(N-S-1)}}\left(\frac{\ell_{1}}{\ell}\right)^{S\Delta}\sim\ell_{\mathrm{D}}^{N-1}\ell_{1}^{-3S+S\Delta}\ell_{3}^{-3(N-S-1)-S\Delta},$
(E.3)

therefore the $(S/N)^{2}$ scales as

$\begin{split}\left(\frac{S}{N}\right)^{2}\sim\int_{\ell_{\mathrm{D}}}^{c\ell_{\mathrm{max}}}d^{2S}\ell_{1}\int_{c\ell_{\mathrm{max}}}^{\ell_{\mathrm{max}}}d^{2(N-1-S)}\ell\,\frac{\ell_{1}^{-6S+2S\Delta}\ell^{-6(N-S-1)-2S\Delta}}{\ell_{1}^{-3S}\ell^{-3(N-S)}}\sim\ell_{\mathrm{max}}^{4-N+S(1-2\Delta)}.\end{split}$
(E.4)

When $S=1$, we recover the single-squeezed limit scaling
$\ell_{\mathrm{max}}^{5-N-2\Delta}$. Moreover, we notice that in the interval
of values $0<\Delta\leq 1/2$, the scaling gets enhanced. Indeed, assuming that
$(S/N)^{2}\propto\ell_{\mathrm{max}}^{p}$, the multiple squeezed limit allows
values of $p>0$ for $N<S(1-2\Delta)+4$. Knowing that the possible values of
$S$ are between $0\leq S\leq N-2$, we may conclude that the maximum scaling
attainable is $p_{\mathrm{max}}=2-2\Delta(N-2)$, therefore it never exceeds
mode counting.

The computation of multiple collapsed limits follows trivially from the result
for a singularly collapsed computation in Eq. (3.18). In this expression, the
integrals from the left and right momenta have the exact form of a product of
the $(S/N)^{2}$ for a $A+1$ spectrum and a $B+1$ spectrum. Thus the scaling
for multiply collapsed (or mixed squeezed-collapsed) can be obtained by
multiplying the limits for the $A+1$ spectrum by the $B+1$ spectrum. Further,
we note that this can be repeated in a nested fashion. E.g., consider a double
collapsed $(N-1)$-spectrum in $3$D, then

$\displaystyle\left(\frac{S}{N}\right)^{2}_{(N)}\sim$
$\displaystyle\int_{k_{\mathrm{min}}}^{ck_{\mathrm{max}}}\mathrm{d}^{3}k_{I}k_{I}^{-6+4\Delta}\int_{k_{\mathrm{min}}}^{ck_{\mathrm{max}}}\mathrm{d}^{3}k_{I^{\prime}}k_{I^{\prime}}^{-6+4\Delta}\int_{ck_{\mathrm{max}}}^{k_{\mathrm{max}}}\mathrm{d}^{3A}k_{R}\frac{k_{R}^{-6A-4\Delta}}{k_{R}^{-3(A+1)}}$
$\displaystyle\times\int_{ck_{\mathrm{max}}}^{k_{\mathrm{max}}}\mathrm{d}^{3B_{1}}k_{R^{\prime}}\,\frac{k_{R^{\prime}}^{-6B_{1}-2\Delta}}{k_{R^{\prime}}^{-3(B_{1}+1)}}\times\int_{ck_{\mathrm{max}}}^{k_{\mathrm{max}}}\mathrm{d}^{3B_{2}}k_{L^{\prime}}\,\frac{k_{L^{\prime}}^{-6B_{2}-2\Delta}}{k_{L^{\prime}}^{-3(B_{2}+1)}}\sim
k_{\mathrm{max}}^{9-8\Delta}.$ (E.5)

where we split $B=B_{1}+B_{2}$. Here, we expect an enhanced scaling for
$\Delta<3/4$. As for the multiple squeezed limits, we can extend the argument
to general $C$ nested collapsed limits where we obtain $(S/N)^{2}\sim
k_{\mathrm{max}}^{3+C(3-4\Delta)}$. In an analogous manner, for the undamped
CMB $(N-1)$-spectrum we find
$(S/N)^{2}\sim\ell_{\mathrm{max}}^{2+2C(1-4\Delta/3)}$, while the damped
regime produces a scaling given by

$\displaystyle\left(\frac{S}{N}\right)^{2}_{(N)}\sim\ell_{\mathrm{max}}^{4-N+4C(1-\Delta)}.$
(E.6)

We recover the scaling for a single collapsed limit,
$(S/N)^{2}\sim\ell_{\mathrm{max}}^{8-N-4\Delta}$. Unlike the case of multiple
squeezed limits, the maximum number of nested collapsed limits $C_{\rm
max}(N)$ depends on the specific $N$-point correlation function considered,
e.g. $C_{\mathrm{max}}(6)=C_{\mathrm{max}}(7)=2$, while the $8$-point
correlation function has $C_{\mathrm{max}}(8)=3$.

Note that while these scalings might look impressive (in particular for an
undamped tracer), as for the single collapsed limit, they must be critically
assessed within a model. It will be challenging to verify these enhanced
scalings in a full Fisher analysis, as the naive computational cost to
evaluate $S/N$ of the $N$-point correlator scales as $\ell_{\rm max}^{2N-1}$.

## References

  * [1] P. D. Meerburg et al., “Primordial Non-Gaussianity”, arXiv:1903.04409.
  * [2] P. Ade et al., “The Simons Observatory: science goals and forecasts”, Journal of Cosmology and Astroparticle Physics 2019 no. 02, (Feb, 2019) 056–056, arXiv:1808.07445.
  * [3] T. Matsumura et al., “Mission Design of LiteBIRD”, Journal of Low Temperature Physics 176 no. 5-6, (Sep, 2014) 733–740, arXiv:1311.2847.
  * [4] K. N. Abazajian et al., “CMB-S4 Science Book, First Edition”, arXiv:1610.02743.
  * [5] S. Hanany et al., “PICO: Probe of Inflation and Cosmic Origins”, arXiv:1902.10541.
  * [6] L. Amendola, S. Appleby, A. Avgoustidis, D. Bacon, T. Baker, M. Baldi, N. Bartolo, A. Blanchard, C. Bonvin, and et al., “Cosmology and fundamental physics with the Euclid satellite”, Living Reviews in Relativity 21 no. 1, (Apr, 2018) , arXiv:1606.00180.
  * [7] P. Marshall et al., “LSST Science Collaborations Observing Strategy White Paper: ”Science-driven Optimization of the LSST Observing Strategy””, Aug, 2017.
  * [8] D. Collaboration, A. Aghamousa, et al., “The DESI Experiment Part I: Science,Targeting, and Survey Design”, 2016.
  * [9] J. Green et al., “Wide-Field InfraRed Survey Telescope (WFIRST) Final Report”, arXiv:1208.4012.
  * [10] D. J. Bacon et al., “Cosmology with Phase 1 of the Square Kilometre Array Red Book 2018: Technical specifications and performance forecasts”, Publications of the Astronomical Society of Australia 37 (2020) .
  * [11] O. Doré et al., “Cosmology with the SPHEREX All-Sky Spectral Survey”, 2015\.
  * [12] W. R. Coulton, P. D. Meerburg, D. G. Baker, S. Hotinli, A. J. Duivenvoorden, and A. van Engelen, “Minimizing gravitational lensing contributions to the primordial bispectrum covariance”, Phys. Rev. D 101 (Jun, 2020) 123504, arXiv:1912.07619.
  * [13] J. C. Hill, “Foreground biases on primordial non-Gaussianity measurements from the CMB temperature bispectrum: Implications for Planck and beyond”, Physical Review D 98 no. 8, (Oct, 2018) .
  * [14] O. Darwish, S. Foreman, M. M. Abidi, T. Baldauf, B. D. Sherwin, and P. D. Meerburg, “Density reconstruction from biased tracers and its application to primordial non-Gaussianity”, arXiv:2007.08472 [astro-ph.CO].
  * [15] N. Bartolo and A. Riotto, “On the non-Gaussianity from recombination”, Journal of Cosmology and Astroparticle Physics 2009 no. 03, (Mar, 2009) 017–017.
  * [16] N. Arkani-Hamed and J. Maldacena, “Cosmological Collider Physics”, arXiv:1503.08043.
  * [17] H. Lee, D. Baumann, and G. L. Pimentel, “Non-Gaussianity as a particle detector”, Journal of High Energy Physics 2016 no. 12, (Dec, 2016) .
  * [18] D. Babich and M. Zaldarriaga, “Primordial bispectrum information from CMB polarization”, Phys. Rev. D 70 (Oct, 2004) 083005.
  * [19] N. Kogo and E. Komatsu, “Angular trispectrum of CMB temperature anisotropy from primordial non-Gaussianity with the full radiation transfer function”, Physical Review D 73 no. 8, (Apr, 2006) , arXiv:astro-ph/0602099.
  * [20] L. Bordin and G. Cabass, “Probing higher-spin fields from inflation with higher-order statistics of the CMB”, Journal of Cosmology and Astroparticle Physics 2019 no. 06, (Jun, 2019) 050–050, arXiv:1902.09519.
  * [21] P. Creminelli, L. Senatore, and M. Zaldarriaga, “Estimators for local non-Gaussianities”, Journal of Cosmology and Astroparticle Physics 2007 no. 03, (Mar, 2007) 019–019, arXiv:astro-ph/0606001.
  * [22] W. R. Coulton, B. Beringue, and P. D. Meerburg, “The primordial information content of Rayleigh Anisotropies”, arXiv:2010.10481.
  * [23] The Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters”, arXiv:1807.06209.
  * [24] J. Maldacena, “Non-gaussian features of primordial fluctuations in single field inflationary models”, Journal of High Energy Physics 2003 no. 05, (May, 2003) 013–013, arXiv:astro-ph/0210603.
  * [25] V. Acquaviva, N. Bartolo, S. Matarrese, and A. Riotto, “Gauge-invariant second-order perturbations and non-Gaussianity from inflation”, Nuclear Physics B 667 no. 1, (2003) 119 – 148, arXiv:astro-ph/0209156.
  * [26] X. Chen, B. Hu, M. xin Huang, G. Shiu, and Y. Wang, “Large primordial trispectra in general single field inflation”, Journal of Cosmology and Astroparticle Physics 2009 no. 08, (Aug, 2009) 008–008.
  * [27] D. Baumann and D. Green, “Signature of supersymmetry from the early universe”, Physical Review D 85 no. 10, (May, 2012) , arXiv:1109.0292.
  * [28] A. Linde, “Hybrid inflation”, Physical Review D 49 no. 2, (Jan, 1994) 748–754, arXiv:astro-ph/9307002.
  * [29] S. Renaux-Petel and K. Turzyński, “Geometrical Destabilization of Inflation”, Physical Review Letters 117 no. 14, (Sep, 2016) , arXiv:1510.01281.
  * [30] E. Pajer, “Building a Boostless Bootstrap for the Bispectrum”, 2020.
  * [31] M. Kamionkowski, T. L. Smith, and A. Heavens, “CMB bispectrum, trispectrum, non-Gaussianity, and the Cramer-Rao bound”, Physical Review D 83 no. 2, (Jan, 2011) , arXiv:1010.0251.
  * [32] L. Smith, Kendrick M. ad Senatore and M. Zaldarriaga, “Optimal analysis of the CMB trispectrum”, arXiv:1502.00635.
  * [33] A. Lewis, A. Challinor, and A. Lasenby, “Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models”, The Astrophysical Journal 538 no. 2, (Aug, 2000) 473–476, arXiv:astro-ph/9911177.
  * [34] E. Komatsu and D. N. Spergel, “Acoustic signatures in the primary microwave background bispectrum”, Phys. Rev. D 63 (Feb, 2001) 063002, arXiv:astro-ph/0005036.
  * [35] 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 no. 05, (May, 2006) 004–004, arXiv:astro-ph/0509029.
  * [36] L. Senatore, K. M. Smith, and M. Zaldarriaga, “Non-Gaussianities in single field inflation and their optimal limits from the WMAP 5-year data”, Journal of Cosmology and Astroparticle Physics 2010 no. 01, (Jan, 2010) 028–028, arXiv:0905.3746.
  * [37] J. R. Fergusson and E. P. S. Shellard, “Shape of primordial non-Gaussianity and the CMB bispectrum”, Physical Review D 80 no. 4, (Aug, 2009) .
  * [38] D. M. Regan, E. P. S. Shellard, and J. R. Fergusson, “General CMB and primordial trispectrum estimation”, Phys. Rev. D 82 (Jul, 2010) 023520, arXiv:1004.2915.
  * [39] H. Lee and C. Dvorkin, “Cosmological angular trispectra and non-Gaussian covariance”, Journal of Cosmology and Astroparticle Physics 2020 no. 05, (May, 2020) 044–044.
  * [40] T. Suyama and M. Yamaguchi, “Non-Gaussianity in the modulated reheating scenario”, Physical Review D 77 no. 2, (Jan, 2008) , arXiv:0709.2545.
  * [41] D. Karagiannis, A. Slosar, and M. Liguori, “Forecasts on Primordial non-Gaussianity from 21 cm Intensity Mapping experiments”, arXiv:1911.03964.
  * [42] The PUMA, A. Slosar et al., “Packed Ultra-wideband Mapping Array (PUMA): A Radio Telescope for Cosmology and Transients”, arXiv:1907.12559 [astro-ph.IM].