(dashed line) theories. For illustration purposes, we have normalised the power spectrum to the featureless one and multiplied it by a factor of 20. We have chosen the overall amplitude $\rho_0$ so that the power spectrum presents features of order a few percents for $\mu = 2$ and is compatible with current bounds from CMB observations. This figure exhibits the presence of correlated oscillations in $2k/k_0$ for the power spectrum and $3k/k_0$ for the bispectrum, characteristic of a sharp feature that eventually decays at smaller scales. It is interesting to observe that the bispectrum envelop reaches a maximum at $k/k_0 \sim 20$. These features are compatible with previous findings in the literature, including general statements about the universal running and model-dependent envelop in [172]. With $\Delta t = 0.1H^{-1}$ and our choice of masses, it is clear that the adiabaticity conditions (8.7) are not satisfied. Nonetheless, it is interesting to see that the effective theory already provides a reasonable description, modestly reproducing the correct amplitude of the feature but replicating the oscillatory profile. For larger masses $m \gg \Delta t^{-1} = 10H$, we expect the power spectrum and the bispectrum in equilateral configurations obtained with the effective theory to be indistinguishable from the ones obtained with the full theory.

## 8.3 Probing Inflationary Landscapes with Cosmological Colliders

We now study the bispectrum shape dependence at a given overall scale. Importantly, the explicit breaking of scale invariance by the time dependence of the quadratic mixing, and therefore of the cubic one by symmetry, induces unobservable contributions that should be subtracted from any bispectrum calculation. Here, we prove that after subtracting these gauge artefacts, the observable bispectrum encodes non-standard cosmological collider signals whose frequency is partially dictated by the background frequency.

**Consistency relation and observable bispectrum.** In its squeezed limit, the observable primordial bispectrum measured in cosmological data, is found by carefully subtracting the leading-order contribution coming from an effective rescaling of the background experienced by short-wavelength fluctuations, in presence of a long-wavelength $\zeta$ mode [185–187]. This gauge artefact—since a rescaling of the background does not constitute an observable quantity—which is present in soft limits of all primordial correlation functions, is determined by so-called consistency relations. For the bispectrum, squeezed configurations are famously dominated by the power spectrum scale-dependence [17, 188, 189] (see also e.g. [73, 75, 190–196] for related works)

$$ \lim_{k_3 \ll k_1} \frac{\langle \zeta_{k_1} \zeta_{k_2} \zeta_{k_3} \rangle'}{P_{\zeta}(k_3)} = -P_{\zeta}(k_s) \frac{\partial \log[k_s^3 P_{\zeta}(k_s)]}{\partial \log k_s} + \mathcal{O}((k_{\ell}/k_s)^2), \quad (8.11) $$

where $P_{\zeta}(k) = \langle \zeta_k \zeta_{-k} \rangle'$, $k_{\ell} = k_3$ and $k_s = k_1 + k_{\ell}/2 = (k_1 - k_2)/2$. When computing the bispectrum in the pure scale-invariant limit, i.e. $\partial \log[k^3 P_{\zeta}(k)]/\partial \log k = n_s - 1 = 0$, the gauge artefact contribution to the bispectrum vanishes in the squeezed limit, and observable effects, which normally show up at next-to-leading order, become automatically the dominant contribution. However, in the presence of features, $n_s(k) \neq 1$ and the consistency relation (8.11) should be subtracted to obtain theoretical predictions that can be readily compared to actual cosmological data. As we will see with a concrete example supported by both an exact numerical calculation and analytical insights, it is crucial to properly subtract the consistency relations in cosmological
collider signals when features are present. Indeed, one could wrongly conclude that the observ-
able signal in the squeezed limit is hidden below a dominant and universal contribution, which
is nothing but the aforementioned gauge artefact.

When the quadratic mixing strength is weak, some analytical insight can be gained by evaluating the bispectrum shape using perturbative techniques. Using the standard in-in formalism [112] (see also [19, 197] for reviews), the bispectrum shape is composed of two contributions $S(k_1, k_2, k_3) = S_1(k_1, k_2, k_3) + S_2(k_1, k_2, k_3)$, where $S_1$ corresponds to a quadratic mixing insertion prior to the cubic interaction and $S_2$ the opposite. Specifying to isosceles configurations $k_1 = k_2 \simeq k_s$ and $k_3 = k_\ell$ where $k_s$ is the short (i.e. hard) momentum and $k_\ell$ is the long (i.e. soft) momentum, each of the two shape contributions have two fundamental permutations, $S_i = S_i^a + S_i^b$ for $i = \{1, 2\}$, corresponding to having either $k_s$ (denoted in the following permutation “a”) or $k_\ell$ (denoted in the following permutation “b”) carried by a $\zeta$ mode attached to the quadratic mixing. Note that each contribution “a” includes two additional internal permutations accounting for the exchange of $k_\ell$ and $k_s$ for both $\zeta$ modes attached to the cubic interaction. Diagrammatically and using the above notations, the total shape is given by the four following contributions$^{49}$

<figure><img src="image_13.png" /></figure>

The contributions $S_1^b$ and $S_2^b$ contain the cosmological collider signals as these channels effectively
probe the super-horizon decay of the massive field for a large hierarchy of scales $k_\ell \ll k_s$. As
we will see in what follows and for the quadratic mixing given in Eq. (8.10), the contribution
$S_1^b$ contains the dominant cosmological collider signal. On the other hand, the contributions $S_1^a$
and $S_2^a$, scaling as $\sim (k_\ell/k_s)^{-1+2n}$ (see insert below), become large and dominant in the squeezed
limit. Note that for an ever-present feature $n=0$, the latter shape contributions are of the
local type. One might expect that the cosmological collider signal is therefore completely out
of reach. However, these large contaminating signals are not observable. Indeed, translating
the bispectrum consistency relation to the dimensionless shape function for isosceles kinematic
configurations, we obtain the remarkable identity

<footnote n="49">⁴⁹Time runs from bottom to top until the end of inflation denoted by the grey line, red lines represent propagators of ζ that are independent—in this perturbative scheme—of those of σ in blue, black dots are quadratic and cubic vertices.</footnote>
as depicted in (8.12), and the integrands read

$$
\begin{align}
\mathcal{J}_1^\text{a}(x, y; \kappa, \lambda) &= \frac{\partial_t \tilde{\rho}}{H}(x) \tilde{\rho}(y) (xy)^{-\frac{1}{2}} \operatorname{Im} \left[ (1 + i\kappa x + \kappa^2(1 + ix)) e^{-i(1+\kappa)x} \right] \operatorname{Im} \left[ H_{i\mu}^{(1)}(x) H_{i\mu}^{(2)}(y) e^{-iy} \right], \\
\mathcal{J}_2^\text{a}(x, y; \kappa, \lambda) &= \tilde{\rho}(x) \frac{\partial_t \tilde{\rho}}{H}(y) (xy)^{-\frac{1}{2}} \sin(x) \operatorname{Im} \left[ H_{i\mu}^{(1)}(x) H_{i\mu}^{(2)}(y) (1 + i\kappa y + \kappa^2(1 + iy)) e^{-i(1+\kappa)y} \right]. \tag{8.17}
\end{align}
$$

Let us now estimate the scalings of these contributions in the squeezed limit $\kappa \ll 1$. The exponential $e^{-iy}$ in the integrand $\mathcal{J}_1^\text{a}$ effectively renders all contributions from the $y \ge 1$ region negligible, so that the integral is dominated by the saddle point around $y \sim 1$. As a result, one can expand in $\kappa y \ll 1$ inside the integral for $\kappa \ll 1$. Similarly inside the second nested integral for $x \le y$, one can take the $\kappa x \ll 1$ limit, which could also be seen from the factor $e^{-i(1+\kappa)x}$ killing all contributions from $(1+\kappa)x \approx x \ge 1$. At leading order in $\kappa$, the integrand $\mathcal{J}_1^\text{a}$ simplifies to

$$
\mathcal{J}_1^\text{a}(x, y; \kappa, \lambda) = - \frac{\partial_t \tilde{\rho}}{H}(x) \tilde{\rho}(y) (xy)^{-\frac{1}{2}} \sin(x) \operatorname{Im} \left[ H_{i\mu}^{(1)}(x) H_{i\mu}^{(2)}(y) e^{-iy} \right]. \quad (8.18)
$$

For the quadratic mixing (8.10), and taking into account the dependence of $\tilde{\rho}$ on $x_0 = \kappa/\lambda$, the overall scaling in the squeezed limit is found to be

$$
S_1^a(\kappa, \lambda) \sim \kappa^{-1+2n} \lambda^{2n} \operatorname{Im} \left[ \alpha_{1,+}^a \kappa^{2i\mu_c} + \alpha_{1,-}^a \kappa^{-2i\mu_c} + \alpha_{1,0}^a \right], \quad (8.19)
$$

where the coefficients $\alpha_{\pm}$ and $\alpha_0$ depend on $\mu, n$ and $\lambda$. Note that in the presence of a slowly-decaying feature in the mixing such that $n \le 1/2$, this shape contribution dominates in the squeezed limit as it grows as function of $\kappa$, even for a process involving the exchange of a heavy field. Notably, the shape oscillates in $\log \kappa$ at a frequency solely set by the background frequency $\mu_c$. Proceeding in the same way, the contribution $S_2^a$ has the same scaling in the squeezed limit. Note that these scalings and the oscillations with frequency $2\mu_c$ can be recovered using the effective single-field Lagrangian (8.4) in the corresponding regime of validity.

Additionally, from

$$
\frac{\partial \tilde{\rho}(x)}{\partial k} = \frac{x}{k^2} \left. \frac{\partial \tilde{\rho}(\tau)}{\partial \tau} \right|_{\tau = -\frac{x}{k}} = \frac{1}{k} \frac{\partial_t \tilde{\rho}(x)}{H}, \qquad (8.20)
$$

we obtain the remarkable identity

$$
\frac{\partial \mathcal{I}}{\partial k} = \frac{1}{k} (\mathcal{J}_1^\text{a} + \mathcal{J}_2^\text{a}) (1 + \mathcal{O}(\kappa^2)). \qquad (8.21)
$$

Lastly, one needs to pull the derivative out of the time integrals. Defining $I(x_0) = \int_0^\infty dx f(x, x_0) \int_x^\infty dy g(y, x_0)$ where $\mathcal{I}(x, y; x_0) = f(x, x_0)g(y, x_0)$, and differentiating with respect to $x_0$, one straightforwardly obtains

$$
\frac{\partial I}{\partial x_0} = \int_0^\infty dx \frac{\partial f}{\partial x_0}(x, x_0) \int_x^\infty dy g(y, x_0) + \int_0^\infty dx f(x, x_0) \int_x^\infty dy \frac{\partial g}{\partial x_0}(y, x_0). \quad (8.22)
$$

Altogether, we obtain the identity

$$
\sum_{i=1,2} S_i^a(\kappa, \lambda) = -\frac{1}{4\kappa} \left. \frac{\partial \log(k^3 \Delta_{\zeta}^2)}{\partial \log k} \right|_{k=k_s} + \mathcal{O}(\kappa), \quad (8.23)
$$

which is precisely the consistency relation at the level of the shape function.
<figure><img src="image_14.png" /><figcaption>Figure 12: Observable shape functions (solid lines) dictated by the mixing Lagrangian (8.3) with the time-dependent mixing $\rho(t)$ given in Eq. (8.10), in the isosceles kinematic configuration $k_1 = k_2$ at a fixed overall scale $\log(k_3/k_0) = 2$ where $k_0 = (aH)(t_0)$ and $t_0$ is the characteristic time of the feature, for a light field $\nu = 1$ (left panel) and a heavy field $\mu = 5$ (right panel), varying the background frequency $\mu_c = 5, 7$, and setting $n = 0$. The observable signal is reconstructed from the computed one (dotted line) upon numerically subtracting the consistency relation (8.13). The shaded region corresponds to the squeezed limit where the templates (8.24) are valid. We have chosen $\rho_0/H = 0.1$ to respect the current bounds on the power spectrum [32].</figcaption></figure>

**Cosmological collider signals.** We focus on the bispectrum arising from the cubic interaction $\dot{\rho}(t)\pi_c\dot{\pi}_c\sigma$ because it is the dominant one for rapid oscillations $\mu_c > 1$, with the quadratic mixing time dependence of Eq. (8.10). We show in Figure 12 the corresponding shape functions for a light and heavy field $\sigma$, for $n = 0$ and different background frequencies $\mu_c$. The dashed line corresponds to the cosmological flow results—hence containing all contributions, even the unobservable ones corresponding to $S_{1,2}^a$ in the perturbative analysis—and the solid lines correspond to the shapes that are directly measurable in cosmological data as they are found after subtracting the consistency relation given solely in terms of the power spectrum. Note that in the cosmological flow framework, both the two- and three-point correlators are computed simultaneously as the flow equations form a closed system. It therefore requires no additional computations to directly obtain the observable signals out of a single computation. Note also that our numerical approach is precise enough (both at the level of the power spectrum and of the bispectrum) to allow for a clean subtraction of dominant unobservable contributions in order to recover a subdominant scaling behaviour. This is readily visible, for example, in the light-field case for $\mu_c = 5$ where the raw signal (dashed line) presents a distorted oscillatory pattern due to the interference of various signals with different frequencies, as well as a quickly growing envelope (as $k_s/k_\ell$) in the squeezed limit. The observable signal after subtracting gauge artefacts (solid line), on the other hand, features a single frequency and is growing at a slower pace (as $(k_\ell/k_s)^{-1/2}$).

A notable feature of these signals is that their frequencies are partially set by the background frequency $\mu_c$. This offers the possibility to probe inflationary landscapes with cosmological collider signals.<fnref n="50" /> Indeed, in the perturbative analysis and considering only the shape contributions $S_1^b$

<footnote n="50">⁵⁰Instead, for cubic interactions of constant strengths, cosmological collider signals are conventional with frequencies only set by the mass of the additional field, despite the time-dependence of the quadratic mixing [200]. This can also be seen in our proof in the insert below.</footnote>
and $S_2^b$, as the other two are removed by the consistency relation, the non-analytical signal lying in the squeezed limit ($k_3 \ll k_1$) of the shape function both for the exchange of a light ($m/H \le 3/2$) and a heavy ($m/H \ge 3/2$) field is found to be

$$S(k_1, k_2, k_3) = \begin{cases} \left(\frac{k_3}{k_1}\right)^{1/2+n-\nu} \mathcal{A} \cos\left(\mu_c \log\left(\frac{k_3}{k_1}\right) + \varphi\right), & (\text{light}) \\ \left(\frac{k_3}{k_1}\right)^{1/2+n} \sum_{\pm} \mathcal{A}_{\pm} \cos\left((\mu \pm \mu_c) \log\left(\frac{k_3}{k_1}\right) + \varphi_{\pm}\right), & (\text{heavy}) \end{cases} \quad (8.24)$$

with $\mu = \sqrt{(m/H)^2 - 9/4}$ and $\nu = i\mu$. Details on the derivation can be found in the insert below. A generic feature of these signals is that the power-law scaling acquires an additional $(k_3/k_1)^n$ suppression compared to the standard one [35]. This is a direct consequence of the dilution of the feature in $(a/a_0)^{-n}$. For the exchange of a light field, in addition to the power-law scaling that depends on the mass, the shape presents oscillations in $\log(k_3/k_1)$ whose frequency is set by the background frequency $\mu_c$. The light-field template is not valid for the massless case ($\nu = 3/2$) and an ever-present feature $n = 0$ as curvature modes never freeze on super-horizon scales. In the case of a heavy-field exchange process, the power-law scaling is saturated and the cosmological collider signals present modulated frequencies in $\mu \pm \mu_c$. These beats are characteristic of interference patterns between the usual oscillations of the heavy field on super-horizon scales and an additional oscillating signal, in this case generated by the background. Cross-correlating the analytic template (8.24) with features in the power spectrum, namely their amplitudes, and the scale-dependence of the bispectrum, can remove the degeneracy between the mass of the additional field and the oscillating background frequency. Remarkably, we find numerically that the templates (8.24) hold in the strong mixing regime. Altogether, these results provide a new way of probing the landscape of inflationary backgrounds through soft limits of cosmological correlators. In particular, for light fields with $\nu > 1/2 + n$, the oscillations that modulate a growing shape envelope in the squeezed limit are expected to generate distinctive oscillations in the scale-dependent galaxy bias. It would also be interesting to find a concrete background model that reproduces these features, as our analysis entirely relies on symmetries at the level of inflationary fluctuations. We leave these open questions for an upcoming work.

**Derivation.**—In this insert, we estimate the scaling of the observable shape contributions $S_1^b$ and $S_2^b$ in the squeezed limit and derive the analytical templates (8.24). After subtracting the consistency relation at leading order in $k_\ell/k_s \ll 1$, the physical shape is given by

$$S(\kappa, \lambda) = -\frac{\pi}{4} \left(\frac{\rho_0}{H}\right)^2 \kappa^{-1} \int_0^\infty dx \int_x^\infty dy \sum_{i=1,2} \mathcal{J}_i^b(x, y; \kappa, \lambda), \quad (8.25)$$

where $\kappa \equiv k_\ell/k_s$ and $\lambda \equiv k_\ell/k_0$. The integrands read

$$\begin{aligned} \mathcal{J}_1^b(x, y; \kappa, \lambda) &= \frac{\partial_t \tilde{\rho}}{H}(x) \tilde{\rho}(y) \kappa^2 (xy)^{-\frac{1}{2}} \operatorname{Im}[(1 + ix)e^{-2ix}] \operatorname{Im}[H_{i\mu}^{(1)}(\kappa x) H_{i\mu}^{(2)}(\kappa y) e^{-i\kappa y}], \\ \mathcal{J}_2^b(x, y; \kappa, \lambda) &= \frac{\partial_t \tilde{\rho}}{H}(x) \tilde{\rho}(y) \kappa^2 (xy)^{-\frac{1}{2}} \sin(\kappa x) \operatorname{Im}[H_{i\mu}^{(1)}(\kappa x) H_{i\mu}^{(2)}(\kappa y)(1 + iy)e^{-2iy}]. \end{aligned} \quad (8.26)$$

We now estimate the scalings and behaviours of these contributions in the squeezed limit $\kappa \ll 1$.
The exponential $e^{-2ix}$ in the integrand $\mathcal{J}_1^b$ selects the integration region $x < 1$, discarding possible resonances between this exponential and the oscillating feature. Note that we have neglected the window function, which corresponds to considering $k_\ell/k_0 \gg 1$, i.e. that both the long and the short modes exit the horizon well after the characteristic scale of the feature. We can then expand in $\kappa x \ll 1$ inside the integral and use the late-time behaviour of the Hankel function $H_{i\mu}^{(1)}(z) \approx A_\mu^+ z^{i\mu} + A_\mu^- z^{-i\mu}$ for $z \ll 1$. After changing variables $y \to \tilde{y} = \kappa y$, the lower bound of the inner integral can be extended up to $\kappa x \to 0$, effectively decoupling the two integrals. Altogether, the overall scaling in the squeezed limit for the shape contribution $S_1^b$ is found to be

$$ S_{1,\text{light}}^{b}(k_1, k_2, k_3) \sim \kappa^{\frac{1}{2}+n-\nu}\lambda^{-2n}\text{Im}[\alpha_{1,+}^{b}\kappa^{i\mu_c} + \alpha_{1,-}^{b}\kappa^{-i\mu_c}], \\ S_{1,\text{heavy}}^{b}(k_1, k_2, k_3) \sim \kappa^{\frac{1}{2}+n}\lambda^{-2n}\text{Im}[\alpha_{1,++}^{b}\kappa^{i(\mu+\mu_c)} + \alpha_{1,+-}^{b}\kappa^{i(\mu-\mu_c)} \\ +\alpha_{1,-+}^{b}\kappa^{-i(\mu-\mu_c)} + \alpha_{1,--}^{b}\kappa^{-i(\mu+\mu_c)}], \qquad (8.27) $$

for the light and the heavy case, respectively. Similarly, for the integrand $\mathcal{J}_2^b$, the exponential $e^{-2iy}$ selects the region $y < 1$ and we can take the $\kappa y \ll 1$ limit inside the integrand. As $x < y$, we can also take the limit $\kappa x \ll 1$. Using the Hankel expansion and at leading order in $\kappa$, the scalings are found to be

$$ S_{2,\text{light}}^{b}(k_1, k_2, k_3) \sim \kappa^{2+2n-2\nu}\lambda^{-2n}\text{Im}[\alpha_{2,+}^{b}\kappa^{2i\mu_c} + \alpha_{2,0}^{b} + \alpha_{2,-}^{b}\kappa^{-2i\mu_c}], \\ S_{2,\text{heavy}}^{b}(k_1, k_2, k_3) \sim \kappa^{2+2n}\lambda^{-2n}\text{Im}[\alpha_{2,++}^{b}\kappa^{2i(\mu+\mu_c)} + \alpha_{2,--}^{b}\kappa^{-2i(\mu+\mu_c)} \\ +\alpha_{2,+-}^{b}\kappa^{2i(\mu-\mu_c)} + \alpha_{2,-+}^{b}\kappa^{-2i(\mu-\mu_c)} + \alpha_{2,0+}^{b}\kappa^{2i\mu_c} \\ +\alpha_{2,0-}^{b}\kappa^{-2i\mu_c} + \alpha_{2,+0}^{b}\kappa^{2i\mu} + \alpha_{2,-0}^{b}\kappa^{-2i\mu} + \alpha_{2,00}^{b}]. \qquad (8.28) $$

The contribution $S_1^b$ is therefore the dominant contribution. Defining $Ae^{i\varphi} = (\alpha_{1,+}^b - \alpha_{1,-}^{b*})/i$, and $A_\pm e^{i\varphi_\pm} = (\alpha_{1,\pm}^b - \alpha_{1,-\mp}^{b*})/i$, we obtain the templates (8.24). Setting aside intricacies related to IR $(x, y \to 0)$ divergences of the corresponding integrals that cancel when adding all contributions, the amplitudes and the phases of these signals can be computed analytically.
# 9 Conclusions and Outlook

In this paper, following a shorter letter [62], we have presented a systematic framework—the cosmological flow—to compute and study primordial correlators. This numerical approach consists in solving differential equations *in time* satisfied by the correlators, and enables us to obtain exact results in regimes that are challenging to probe analytically. We have exemplified this method with an extensive study of correlators coming from the exchange of a massive scalar particle, both at strong mixing—adopting a non-perturbative treatment of quadratic interactions—and in non scale-invariant theories. This development provides an efficient way to study the rich phenomenology of primordial non-Gaussianities as it makes predictions straightforward. In order for the community to progress faster in these challenging directions, we have decided to publicly release the numerical implementation of our new approach as a freely available resource.<fnref n="51" />

In more details, we have reformulated the computation of tree-level in-in correlators by deriving a closed set of differential equations in time, providing direct insights into the bulk dynamics of equal-time correlators. While the structure of these flow equations is universal, the specific form of the (time-dependent) coefficients is theory dependent. These theories, formulated at the level of fluctuations, can include any number of scalar degrees of freedom with arbitrary dispersion relations and masses, coupled through any type of time-dependent interactions. Our approach is convenient as it enables one to decouple from a specific background and to focus directly on observables. We have established the already-existing flow equations for the two- and three-point correlators by introducing handy diagrammatic rules. This has the advantage of following the usual Leibniz product rule of differentiation and making the contributions coming from an exchange or a contact diagram manifest. Using the diagrammatic representation of the flow equations, we have also derived the flow equations for all tree-level n-point correlators.

In a second phase, we have applied our formalism to the study of correlators emerging from a theory where the Goldstone boson of broken time translations is coupled to an additional massive scalar field. After an extensive analysis of the quadratic theory, we have provided an exhaustive dictionary of strong coupling scales, perturbativity bounds and naturalness criteria for the theory under scrutiny. Such glossary is crucial for theoretical consistency and observability of primordial signals. We have provided exact results for the size of the non-Gaussian signal in the entire parameter space, including the strong mixing regime. Then, we have studied the shape dependence, varying the quadratic mixing strength and the speed of sound. We found that the resulting shapes have interesting new phenomenology compared to the known previous cases, namely the shapes coming from a single-exchange correlator at weak mixing. In particular, we have highlighted two features that had not been appreciated enough. First, as the quadratic mixing becomes stronger, the additional field acquires an *effective mass* that sets the frequency of the cosmological collider signal. We have also presented exact results for the amplitude of the non-analytic signal, pointing out that the maximal amplitude is achieved when the mixing strength is comparable to the mass of the exchanged particle. Second, a large quadratic mixing results in a slight boost for the heavy fluctuation upon crossing the horizon, due to backreaction from the

<footnote n="51">⁵¹In this repository, we have also gathered numerical results for power spectra and non-Gaussianities computed with our code, for different theories and parameter values. We plan to turn this into a bank of theoretical data to be nurtured by the work of the community.</footnote>
curvature fluctuation. This interesting effect is intrinsic to the strong mixing regime, and creates
a phase shift in the cosmological collider signal. Differences between the weak and the strong
mixing regimes can be appreciated by examining videos of the building of the cosmological collider
signal as time passes, a unique opportunity offered by our approach. We have also unveiled a
large cosmological collider signal, built from a non-linear interaction automatically generated by
symmetry as soon as the *additional field* has a reduced intrinsic speed of sound. Finally, we have
applied the cosmological flow to non scale-invariant scenarios, where the quadratic mixing picks
up a time dependence. We have provided exact results for the power spectrum and bispectrum—
both its scale and shape dependence—in the regime where a single-field effective description
breaks down. With minimal assumptions—a non-linearly realised symmetry and an oscillating
background—, we have shown that the observable cosmological collider signal encodes information
about the background, offering a novel approach to probe the inflationary landscape. This generic
phenomenon, beyond the standard lore, breaks the direct link between the frequency of the
cosmological collider signal and the mass of an additional field.

Heading towards a universal program to generate theoretical data, the cosmological flow has a number of benefits. First, it requires no analytical nor numerical mode functions as it focuses directly on correlators. Second, it effortlessly deals with the early-time $i\epsilon$ prescription as the correlators are evolved on the real time axis. Third, the procedure of setting the flow equations given a certain Lagrangian is straightforward. Finally, the cosmological flow systematically works for all $n$-point correlators. With this formalism at hand, we have paved the way towards a systematic investigation of the rich and fascinating subject of primordial correlators. Since observational progresses are expected to be made in the near future, the coming years will lead to a wealth of observational data and constraints on the statistical properties of primordial fluctuations. It is therefore important to provide an adequate amount of theoretical benchmarks, and establish the capability for an efficient exploration of the vast landscape of possible signals.<fnref n="52" /> In this respect, the cosmological flow can both extend the reach of primordial correlators phenomenology and complement analytical computations by providing useful checks and generating intuition.

Finally, yet importantly, our work suggests a number of directions for future exploration:

*   **Small-scale phenomenology.** Primordial density fluctuations on cosmological scales are probed by observations of the CMB and LSS. Density fluctuations on smaller scales that are forged in the “dark era” of inflation are not directly observable by themselves. Yet, they shape many properties of new observables like gravitational wave backgrounds and primordial black holes. In this respect, we highlight that the cosmological flow approach is not restricted to the largest cosmological scales but is also perfectly applicable to these questions. As the inflationary dynamics involved in this context often relies on a strong breaking of scale invariance and substantial non-linear effects which are difficult to handle analytically, a systematic approach like the cosmological flow provides an essential tool to illuminate the dark era of inflation.

<footnote n="52">⁵²More ambitiously, we could include the cosmological flow routine in the already-existing chain of late-time cosmological tools. This would automate the generation of theoretical primordial data that can be directly used for CMB or LSS observables, hence extending the cosmological flow to the late-time Universe.</footnote>
* **Spinning fields.** It would be interesting to extend the cosmological flow to include spinning particles. On the one hand, in conventional setups of weak mixing and under the lamppost of weakly broken scale-invariance, it is known that the cosmological collider signal arising from the exchange of a spinning particle displays a characteristic angular dependence. On the other hand, computing correlators with external spin (for example gravitons, hence accounting for gravitational waves) is notoriously difficult because it is challenging to track the detailed time evolution of the physics in the bulk. The cosmological flow would encapsulate both cases as the flow equations couple all correlators, that need to be solved at once. More practically, implementing particles with spin within the cosmological flow requires a careful treatment of the tensorial structures of the spinning fields, and hence necessitates additional indices that would account for polarisation states. A possible direction for future work would be to relax the assumptions of weak mixing and/or scale-invariance, and to extend the phenomenology of spinning fields to regimes that have not yet been explored.

* **Path integral formulation.** It would also be preferable to derive the flow equations using the path integral formulation instead of the canonical operator formalism. The reasons for this a priori cumbersome change of perspective are at least twofold. First, most of the theories at the level of the fluctuations are defined by a Lagrangian (or equivalently an action), whereas the current formalism required switching to the Hamiltonian, expanding the Dyson series, and doing field contractions by hand order by order. Second, from a particle physics perspective, a path integral formulation of the flow approach would allow us to include gauge theories. Indeed, when working with spinning fields, one needs to deal with problems such as gauge redundancies. In flat space, these complications are avoided by using Feynman diagrams to express the $S$-matrix elements order by order in perturbation theory. Using a path integral approach and importing standard techniques from flat-space QFT, one could expect to render the implementation of fields with spin within the cosmological flow more limpid.

* **Loop-level.** Mostly for technical reasons, the large and rich phenomenology of inflationary correlators mainly comes from explicit calculations at tree level. *Per contra*, pushing the available technology further in perturbation theory and developing new techniques to compute loops is of primary importance. Practically, it would extend the reach of the cosmological collider phenomenology by including particles within and beyond the Standard Model (it is known that spin-1/2 fermions or the Higgs field contribute to primordial non-Gaussianities at loop level). Yet, aside from isolated examples, there is no systematic treatment of loops in cosmological spacetimes. That is why it would be interesting to derive the flow equations at one-loop level. One can already notice that extending the cosmological flow to loop-level would require way more computational capability. For example, computing the two-point correlators at one-loop level requires to solve the flow equations for the three-point functions (working at loop-level precisely means that we do not neglect anymore non-linearities of the three-point correlators on the two-point correlators), which requires first the solutions of the two-point correlators flow equations at tree-level (this is equivalent to solving the quadratic theory). Therefore, one is not able to solve all equations at once—as we have done so far—but rather order by order.
* **Unitarity.** More conceptually, the cosmological flow—by evolving correlators in time—brings the question of unitarity. Indeed, a unitary time evolution ensures that the total probability is conserved over time, and therefore the theory is complete in the sense that it makes self-contained predictions for observables. In flat-space QFT, important properties follow from unitarity. For example, correlation functions must *factorize* into products of lower-order ones. Within the cosmological flow, we have seen that the time evolution of correlators is *sourced* by products of lower-order ones. This unforeseen concordance implores a better understanding of the status of the flow equations and their flat-space analogues. Finally, it would be interesting to see how non-unitary effects arising from the interference between the two branches of the in-in contour—such as decoherence and dissipation—can be incorporated in the flow equations.

**Acknowledgements.** We are grateful to Angelo Caravano, Sebastián Céspedes, Xingang Chen, Paolo Creminelli, Tanguy Grall, Sadra Jazayeri, Liam McAllister, Scott Melville, David Mulryne, Toshifumi Noumi, Enrico Pajer, Gui Pimentel, David Seery, David Stefanyszyn, Xi Tong, Dong-Gang Wang, Lukas Witkowski and Yuhang Zhu for helpful discussions. We would like to especially thank David Mulryne and David Seery for enlightening conversations when this work was initiated. DW and SRP are supported by the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement No 758792, Starting Grant project GEODESI). L.P. acknowledges funding support from the Initiative Physique des Infinis (IPI), a research training program of the Idex SUPER at Sorbonne Université. This article is distributed under the Creative Commons Attribution International Licence (CC-BY 4.0).
# A Details on the Cosmological Flow

In this appendix, we provide additional materials on the cosmological flow that have been used in this work. In Section A.1, we give details about the Fourier summation. We then derive the initial conditions for the flow equations of the two- and three-point correlators in Section A.2 and A.3 respectively, generalising the results of [121] to arbitrary cubic interactions.

## A.1 Extended Fourier Summation

We deal in full generality with $N$ scalar degrees of freedom. Latin letters $(a, b, ...)$ run over phase-space coordinates from 1 to $2N$. Greek letters $(\alpha, \beta, ...)$ run over all fields or their conjugate momenta from 1 to $N$.

**Fourier indices.** We have seen in the main text that it is helpful to condense the notations and the Fourier integrals in the summation convention. We indicate that this extended interpretation is used by typesetting the labels to which it applies in a sans serif face. An index contraction such as $A_a B^a$ reads

$$A_a B^a = \sum_a \int \frac{d^3 k_a}{(2\pi)^3} A_a(\mathbf{k}_a) B^a(\mathbf{k}_a). \qquad (A.1)$$

Note that we position indices to respect the normal rules for covariant expressions. There is an extra complexity with $\mathbf{k}$-labels included because Fourier-space expressions sometimes produce the $\delta$-function $\delta_{ab} = (2\pi)^3 \delta_{ab} \delta^{(3)}(\mathbf{k}_a + \mathbf{k}_b)$. Within a Fourier integral, this reverses the sign of a $\mathbf{k}$-label, which we indicate by decorating the label with a bar as in $B^{\bar{a}}$. Hence,

$$A_a B^{\bar{a}} = \sum_a \int \frac{d^3 k_a}{(2\pi)^3} A_a(\mathbf{k}_a) B^a(-\mathbf{k}_a). \qquad (A.2)$$

A contraction with $\delta^a_b$ will bar an index $A^{\bar{a}} = \delta^{ab} A_b$. When dealing with field Greek indices $(\alpha, \beta)$ running from 1 to $N$, we also indicate the existence of an extended Fourier summation by typesetting the labels in sans serif face $(\alpha, \beta)$.

**Commutator in Fourier space.** With $X^a = (\varphi^\alpha, p^\beta)$ being the phase-space vector in the Heisenberg picture, the canonical commutation algebra has a compact expression

$$[X^a, X^b] = i\epsilon^{ab}, \qquad (A.3)$$

where $\epsilon^{ab} \equiv (2\pi)^3 \delta^{(3)}(\mathbf{k}_a + \mathbf{k}_b) \epsilon^{ab}$ and the $2N \times 2N$ matrix can be written in block form

$$\epsilon^{ab} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}. \qquad (A.4)$$

We also assume that phase-space indices are organised so that a block of field labels are followed by a block of momentum labels, in the same order.
## A.2 Two-point Correlator Initial Conditions

The derivation of initial conditions for the two-point correlators is done using the path integral formalism [42] instead of the canonical in-in formalism, phrased in operator language.<fnref n="53" /> In line with the point of view of the “cosmological flow”, this approach circumvents the use of mode functions.

**Path integral formalism.** The quadratic theory that we consider (see Section 3.3) is

$$S^{(2)} = \frac{1}{2} \int dt a^3 \left( \bar{\Delta}_{\alpha} (\dot{\varphi}^{\alpha})^2 - c_{\alpha}^2 \frac{k^2}{a^2} (\varphi^{\alpha})^2 \right). \qquad (A.5)$$

After integrating by parts while assuming the boundary terms vanish at infinity and introducing conformal time $d\tau = dt/a$, the above action in real space reads

$$S = -\frac{1}{2} \int d\tau \, d^3x \, a^2 \varphi^\alpha \mathcal{O}_\alpha \varphi^\alpha, \quad \text{with} \quad \mathcal{O}_\alpha = \bar{\Delta}_\alpha \left( \partial_\tau^2 + 2 \frac{a'}{a} \partial_\tau \right) - c_\alpha^2 \partial_i^2, \qquad (A.6)$$

where a prime denotes derivative with respect to conformal time. It can be shown that the partition function can be recast into two copies of path integrals, one going forward in time for the time-ordered factor in Eq. (2.7), and the other going backward in time for the anti-time ordered factor. As a consequence, two sets of field configurations $\varphi_+^\alpha$ and $\varphi_-^\alpha$ are introduced<fnref n="54" />, and the in-in generating functional can be written in the following way

$$Z = \int [\mathcal{D}\varphi_+^\alpha \mathcal{D}\varphi_-^\alpha] \exp \left\{ -\frac{i}{2} \int_{\tau_0}^{\tau} d\tau d^3x a^2 \bar{\varphi}^\alpha \begin{pmatrix} \mathcal{O} & -\mathcal{O} \\ -\mathcal{O} & \mathcal{O} \end{pmatrix}_\alpha \varphi^\alpha \right\} \\ \times \prod_{\alpha, x} \delta \left( \varphi_+^\alpha(\tau, x) - \varphi_-^\alpha(\tau, x) \right), \qquad (A.7)$$

where, in this formula, $\varphi^\alpha = (\varphi_+^\alpha, \varphi_-^\alpha)$ and $\bar{\varphi}^\alpha$ denotes the transpose vector. The time $\tau_0$ should be set in the infinite past and $\tau$ is the time at which we want to initialise the two-point correlation functions. The $\delta$-function enforces that the fields coincide at the upper bound $\varphi_+^\alpha(\tau, x) = \varphi_-^\alpha(\tau, x)$. We define the time-ordered two-point function of fields as

$$\Delta_{++}^{\alpha\beta}(\tau_1, x; \tau_2, y) = \langle T \varphi_+^\alpha(\tau_1, x) \varphi_+^\beta(\tau_2, y) \rangle, \qquad (A.8)$$

with similar definitions for $\Delta_{+-}$, $\Delta_{-+}$ and $\Delta_{--}$. The theory being quadratic and diagonal, it is an easy task to compute the $\Delta_{\pm\pm}$. They are obtained by formally inverting the quadratic operator $\mathcal{O}$,

$$\begin{pmatrix} \mathcal{O} \\ -\mathcal{O} \end{pmatrix} \begin{pmatrix} \Delta_{++} & \Delta_{+-} \\ \Delta_{-+} & \Delta_{--} \end{pmatrix} = -\frac{i}{a^2} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \delta(\tau_1 - \tau_2) \delta^{(3)}(x - y). \qquad (A.9)$$

<footnote n="53">⁵³We will closely follow the derivation in [145].</footnote>

<footnote n="54">⁵⁴One also needs to introduce two sets of conjugate momenta $p_+^\alpha$ and $p_-^\alpha$ but the resulting generating functional is quadratic in the momenta so that the path integral is Gaussian and can be evaluated exactly. In this way one can trade a Hamiltonian formulation for a Lagrangian one.</footnote>
These inhomogeneous equations can be solved exactly in Fourier space where the spatial dependence becomes diagonal. With the proper Bunch-Davies vacuum boundary conditions and imposing that $\Delta_{++}$ (in Fourier space) oscillates with positive frequency in the far past, the solution is given by

$$ \Delta_{++}^{\alpha\beta} = (2\pi)^3 \delta^{\alpha\beta} \delta^{(3)}(\mathbf{k}_1 + \mathbf{k}_2) \bar{\Delta}_{\alpha} \frac{H^2}{2c_{\alpha}^3 k^3} (1 + ic_{\alpha}k\tau_1)(1 - ic_{\alpha}k\tau_2)e^{ic_{\alpha}k(\tau_2 - \tau_1)}, \quad (A.10) $$

assuming $\tau_2 < \tau_1$. Here, $k = |\mathbf{k}_1| = |\mathbf{k}_2|$, and there is no sum on the repeated $\alpha$ index. We also recall that the index position is meaningless because the field space for fluctuations is flat i.e. field-space indices are lowered and raised with the flat Euclidean metric $\delta_{\alpha\beta}$. The remaining time-ordered two-point correlators can then be found. The $\delta$-function in (A.7) enforces that $\Delta_{-+}(\tau, \tau_2) = \Delta_{++}(\tau, \tau_2)$ for all $\tau_2$ when $\tau_1 = \tau$ is the external time. From the properties of the path integral, $\Delta_{--}$ and $\Delta_{+-}$ are Hermitian conjugates of $\Delta_{++}$ and $\Delta_{-+}$, respectively.

**Equal-time correlators.** From the above solutions, we can get the equal-time two-point correlators of fields by setting $\tau_1 = \tau_2 = \tau$. Trading the time-dependence for a scale factor dependence using $\tau = -1/aH$, we obtain

$$ \langle \varphi_{\mathbf{k}}^{\alpha}(\tau) \varphi_{-\mathbf{k}}^{\beta}(\tau) \rangle' = \delta^{\alpha\beta} \bar{\Delta}_{\alpha} \frac{H^2}{2c_{\alpha}^3 k^3} \left( 1 + \frac{c_{\alpha}^2 k^2}{a^2 H^2} \right), \quad (A.11) $$

where the prime denotes that we have extracted the momentum conservation $\delta$-function. Similarly one can derive the equal-time (rescaled) momentum-field and momentum-momentum two-point correlators. Eventually, one finds the two-point correlator initial conditions in Eq. (3.31), expressed in terms of the diagonal matrix $\Delta_{\alpha\beta} = \delta^{\alpha\beta} \bar{\Delta}_{\alpha}$ (no sum).

### A.3 Three-point Correlator Initial Conditions

Let us now use the results of the previous section to derive the initial conditions for the three-point correlators. We are interested in computing the three-point correlators at some initial conformal time $\tau_0$ in the deep past with the Hamiltonian (3.6). Expanding the exponentials to leading-order in Eq. (2.7) leads to

$$ \langle \mathbf{X}^a(\tau_0, \mathbf{k}_1) \mathbf{X}^b(\tau_0, \mathbf{k}_2) \mathbf{X}^c(\tau_0, \mathbf{k}_3) \rangle = i \left\langle \int_{-\infty}^{\tau_0} d\tau \left[ H_I(\tau), \mathbf{X}^a(\tau_0, \mathbf{k}_1) \mathbf{X}^b(\tau_0, \mathbf{k}_2) \mathbf{X}^c(\tau_0, \mathbf{k}_3) \right] \right\rangle, \quad (A.12) $$

where the expectation value is evaluated with respect to the usual vacuum of the free theory and $H_I$ is the full cubic Hamiltonian that can be read of Eq. (3.6). To condensate notations, we will now on write the momentum dependency of the operators as subscripts. Assuming it is clear at which time the operators are evaluated, we will also omit the time dependency of the operators. By Fourier transforming the operators in $H_I$, the first term in the commutator becomes

$$ \langle \mathbf{X}_{\mathbf{k}_1}^a \mathbf{X}_{\mathbf{k}_2}^b \mathbf{X}_{\mathbf{k}_3}^c \rangle = i \int_{-\infty}^{\tau_0} d\tau H_{\text{def}} \int \frac{\Pi_i}{(2\pi)^9} \frac{d^3 q_i}{(2\pi)^3} \delta^{(3)}(\Sigma_i q_i) \langle X_{\mathbf{q}_1}^{\text{d}} X_{\mathbf{q}_2}^{\text{e}} X_{\mathbf{q}_3}^{\text{f}} X_{\mathbf{k}_1}^a X_{\mathbf{k}_2}^b X_{\mathbf{k}_3}^c \rangle. \quad (A.13) $$
We have written the integrals over Fourier modes to be explicit—this time not using the Fourier summation convention—to be clear about momentum dependencies. We can now perform Wick contractions between the different operators to rewrite the expectation value as a combination of two-point functions. Discarding disconnected contributions, we obtain

$$ \langle \mathbf{X}_{k_1}^a \mathbf{X}_{k_2}^b \mathbf{X}_{k_3}^c \rangle = i \int_{-\infty}^{\tau_0} d\tau H_{\text{def}} \int \frac{\Pi_i \mathrm{d}^3 q_i}{(2\pi)^9} (2\pi)^3 \delta^{(3)}(\Sigma_i \mathbf{q}_i) \left\{ \langle X_{q_1}^d X_{k_1}^a \rangle \langle X_{q_2}^e X_{k_2}^b \rangle \langle X_{q_3}^f X_{k_3}^c \rangle + \text{5 perms} \right\}. \quad (A.14) $$

We now need to specify the form of $H_{\text{def}}$. For pedagogical reasons before treating the general case, we start with a simple example to illustrate how the machinery works.

**Preliminaries.** Let us treat the concrete case of a cubic field-field-field interaction with no derivative for simplicity. Using the notations in Eq. (3.6), the cubic Hamiltonian reads

$$ H_{\text{def}} \mathbf{X}^d \mathbf{X}^e \mathbf{X}^f = A_{\alpha\beta\gamma} \varphi^\alpha \varphi^\beta \varphi^\gamma, \quad (A.15) $$

where we recall that the Latin indices run over the phase-space variables and the Greek indices run over field-space variables. Taking into account the possible non-canonical normalisation of the fields in the $\Delta_{\alpha\beta}$ tensor and using the time-ordered two-point function (A.10) derived in the previous section, the three-point initial correlator reads

$$ \langle \varphi_{\mathbf{k}_1}^{\alpha} \varphi_{\mathbf{k}_2}^{\beta} \varphi_{\mathbf{k}_3}^{\gamma} \rangle = 6i (2\pi)^3 \delta^{(3)}(\Sigma_i \mathbf{k}_i) \frac{H^6}{8e_3^3} \Delta_{\alpha\alpha} \Delta_{\beta\beta} \Delta_{\gamma\gamma} \times (1 + ik_1\tau_0)(1 + ik_2\tau_0)(1 + ik_3\tau_0)e^{-ik_t\tau_0} \\ \times \left\{ \int_{-\infty}^{\tau_0} d\tau \Delta_{\alpha\alpha}(\tau) \Delta_{\beta\beta}(\tau) \Delta_{\gamma\gamma}(\tau) A_{\alpha\beta\gamma}(\tau) (1 - ik_1\tau)(1 - ik_2\tau)(1 - ik_3\tau)e^{ik_t\tau} \right\} + \text{c.c.}, \quad (A.16) $$

where $k_t = k_1 + k_2 + k_3$. For convenience in later calculations, we also introduce the other symmetric polynomials

$$ e_2 \equiv k_1 k_2 + k_1 k_3 + k_2 k_3, \quad e_3 \equiv k_1 k_2 k_3. \quad (A.17) $$

Note that one needs to carry the $\Delta$ dependency both for external (at time $\tau_0$) and internal (at time $\tau$) operators. There is no summation for the repeated indices $\alpha, \beta$ and $\gamma$. If we assume that the background parameters are slowly varying near $\tau = \tau_0$, then we can Taylor expand the third-order vertex function to first order $[\Delta_{\alpha\alpha}\Delta_{\beta\beta}\Delta_{\gamma\gamma}A_{\alpha\beta\gamma}](\tau) \approx [\Delta_{\alpha\alpha}\Delta_{\beta\beta}\Delta_{\gamma\gamma}A_{\alpha\beta\gamma}](\tau_0)$ and remove it from the time integral. In the end, the three-point correlator is given by

$$ \langle \varphi_{\mathbf{k}_1}^{\alpha} \varphi_{\mathbf{k}_2}^{\beta} \varphi_{\mathbf{k}_3}^{\gamma} \rangle = 6 \Delta_{\alpha\alpha} \Delta_{\beta\beta} \Delta_{\gamma\gamma} i (2\pi)^3 \delta^{(3)}(\Sigma_i \mathbf{k}_i) \frac{H^6}{8e_3^3} \Delta_{\alpha\alpha} \Delta_{\beta\beta} \Delta_{\gamma\gamma} A_{\alpha\beta\gamma} \\ \times (1 + ik_1\tau_0)(1 + ik_2\tau_0)(1 + ik_3\tau_0)e^{-ik_t\tau_0} \\ \times \int_{-\infty}^{\tau_0} d\tau \left\{ (1 - ik_1\tau)(1 - ik_2\tau)(1 - ik_3\tau)e^{ik_t\tau} \right\} + \text{c.c.} \quad (A.18) $$
In the first line, we have collected all “constant” terms. These terms include the normalisation of the external and internal operators encoded in the $\Delta_{\alpha\beta}$ tensor and the third-order vertex kernel from the considered cubic interaction. The remaining factor appears in all three-point correlators. The second line contains the external polynomial and the third line the internal polynomial. This simple example can now be generalised to the computation of all three-point correlators including all cubic interactions.

**External polynomials.** The external polynomial depends on the correlator considered. There are four different types of such polynomials. Formally setting $\tau_1 = \tau_0$ and $\tau_2 = 0$ in the time-ordered two-point function (A.10), we see that the operators $\varphi^\alpha$ and $p^\alpha$ contribute the following polynomials

$$ \varphi^{\alpha} \rightarrow (1 + ik\tau_0)e^{-ik\tau_0}, \quad p^{\alpha} \rightarrow k^2\tau_0 e^{-ik\tau_0}. \qquad (A.19) $$

Using these expressions, it is an easy task to derive the external polynomials

$$ \begin{aligned} \langle \varphi_{k_1}^{\alpha} \varphi_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle & \rightarrow (1 + ik_t \tau_0 - e_2 \tau_0^2 - ie_3 \tau_0^3) e^{-ik_t \tau_0}, \\ \langle p_{k_1}^{\alpha} \varphi_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle & \rightarrow (k_1^2 \tau_0 + ik_1^2 (k_2 + k_3) \tau_0^2 - e_3 k_1 \tau_0^3) e^{-ik_t \tau_0}, \\ \langle p_{k_1}^{\alpha} p_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle & \rightarrow (k_1^2 k_2^2 \tau_0 + ie_3 k_1 k_2 \tau_0^2) e^{-ik_t \tau_0}, \\ \langle p_{k_1}^{\alpha} p_{k_2}^{\beta} p_{k_3}^{\gamma} \rangle & \rightarrow e_3^2 \tau_0^3 e^{-ik_t \tau_0}. \end{aligned} \qquad (A.20) $$

**Internal polynomials.** The internal polynomials depend on the cubic interaction we consider. From Eq. (3.6), we see that we have four different cubic vertex integrals to perform, each of them depending on the third-order kernels

$$ \begin{aligned} A_{\alpha\beta\gamma} &= (2\pi)^3 \delta^{(3)}(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) A_{\alpha\beta\gamma}, \\ B_{\alpha\beta\gamma} &= (2\pi)^3 \delta^{(3)}(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) B_{\alpha\beta\gamma}, \\ C_{\alpha\beta\gamma} &= (2\pi)^3 \delta^{(3)}(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) C_{\alpha\beta\gamma}, \\ D_{\alpha\beta\gamma} &= (2\pi)^3 \delta^{(3)}(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) D_{\alpha\beta\gamma}. \end{aligned} \qquad (A.21) $$

When choosing an initial time well inside the horizon, we want to keep only the highest-order terms in $\tau$. Although the computation of the internal polynomials is straightforward, we still need to deal with higher-order derivative interactions. Thus in Fourier space, the third-order kernels in Eq. (A.21) above may contain explicit $(k^2/a^2)^n$ dependencies for some integer $n$ that needs to be tracked when evaluating the integrals. In order to take into account the scale factor dependency in such terms, we split the kernels according to the number of spatial derivatives they contain in the following way

$$ \begin{aligned} A_{\alpha\beta\gamma} &= A_{\alpha\beta\gamma}^{(0)} + \frac{1}{a^2} A_{\alpha\beta\gamma}^{(2)} + \frac{1}{a^4} A_{\alpha\beta\gamma}^{(4)} + \dots, & B_{\alpha\beta\gamma} &= B_{\alpha\beta\gamma}^{(0)} + \frac{1}{a^2} B_{\alpha\beta\gamma}^{(2)} + \frac{1}{a^4} B_{\alpha\beta\gamma}^{(4)} + \dots, \\ C_{\alpha\beta\gamma} &= C_{\alpha\beta\gamma}^{(0)} + \frac{1}{a^2} C_{\alpha\beta\gamma}^{(2)} + \frac{1}{a^4} C_{\alpha\beta\gamma}^{(4)} + \dots, & D_{\alpha\beta\gamma} &= D_{\alpha\beta\gamma}^{(0)} + \frac{1}{a^2} D_{\alpha\beta\gamma}^{(2)} + \frac{1}{a^4} D_{\alpha\beta\gamma}^{(4)} + \dots, \end{aligned} \qquad (A.22) $$

where we have absorbed all momentum dependency in the newly defined kernels. We indicate with the superscript $2n$ ($n \ge 0$) the number of spatial derivatives that appears in the cubic
interaction. Let us start with the internal polynomial involving the $A_{\alpha\beta\gamma}$ tensor elements. For $n \ge 0$, it reads

$$ 
\begin{aligned}
\underline{A}_{\alpha\beta\gamma}^{(2n)} &= \frac{1}{2} \int_{-\infty}^{\tau_0} d\tau a^4 \frac{A_{\alpha\beta\gamma}^{(2n)}}{a^{2n}} (1 - ik_1\tau)(1 - ik_2\tau)(1 - ik_3\tau) e^{ik_t\tau} + \text{5 perms} \\
&= \frac{-A_{\alpha\beta\gamma}^{(2n)}}{2i^{3-2n}H^{4-2n}} \left\{ \frac{e_3}{k_t^{2n}} \Gamma(2n, -ik_t\tau_0) + \frac{e_2}{k_t^{2n-1}} \Gamma(2n-1, -ik_t\tau_0) \right. \\
&\qquad \left. + \frac{1}{k_t^{2n-2}} \Gamma(2n-2, -ik_t\tau_0) + \frac{1}{k_t^{2n-3}} \Gamma(2n-3, -ik_t\tau_0) \right\} + \text{5 perms},
\end{aligned}
\quad (A.23) $$

where $\Gamma(m, z)$ for some integer $m$ and a complex entry $z$ is the upper incomplete $\Gamma$-function. Details about this function can be found in the following insert.

Incomplete $\Gamma$-function.—The incomplete $\Gamma$-function is defined as

$$ \Gamma(s, z) = \int_{z}^{+\infty} t^{s-1} e^{-t} dt, \quad (A.24) $$

with $s$ being a complex parameter with $\mathcal{R}(s) > 0$. This function is defined without restrictions on the integration paths in the complex plane. In our case, we are interested in the case where $s$ is an integer (positive or negative) and $z$ is purely imaginary. For $s = n \ge 1$, the incomplete $\Gamma$-function can be expressed as a finite sum

$$ \Gamma(n, z) = (n - 1)! \sum_{j=0}^{n-1} \frac{z^j}{j!} e^{-z}, \quad (A.25) $$

with the special value $\Gamma(1, z) = e^{-z}$. For a fixed $n \le 0$, $\Gamma(n, z)$ has an asymptotic expansion – with infinite terms – at large $z$ of the form

$$ \Gamma(n, z) = z^{n-1} e^{-z} \left( \sum_{j=0}^{m-1} \frac{u_j}{z^j} + R_m(n, z) \right), \quad m = 1, 2, \dots, \quad (A.26) $$

where the series coefficients $u_j$ are defined as

$$ u_j = (-1)^j (1-n)_j = (n-1)(n-2)\dots(n-j), \quad (A.27) $$

and $R_m(n, z)$ is the remainder that decays as $\mathcal{O}(1/z^m)$.

From inspecting Eq. (A.23), we see that $\underline{A}^{(0)}$ and $\underline{A}^{(2)}$ have an infinite number of terms when expanding at large $|\tau_0|$. Because one cannot implement an infinite number of terms numerically, we choose to keep only the two dominant powers of $\tau_0$. Bearing in mind that this internal polynomial will eventually be multiplied with the external polynomial, this series truncation is an approximation. Indeed, the corrections to the leading behaviour of higher-order derivative interactions will not be of the same order as the leading behaviour of lower-order in derivative
interactions. Keeping only the two leading-order terms in $\tau$, the $\underline{A}$ internal polynomials read

$$ \underline{A}_{\alpha\beta\gamma}^{(2n)} = \frac{A_{\alpha\beta\gamma}^{(2n)}}{2H^{4-2n}k_t} \left\{ e_3\tau_0^{2n-1} + i\tau_0^{2n-2} \left[ e_2 + \frac{(2n-1)e_3}{k_t} \right] + \mathcal{O}(\tau_0^{2n-3}) \right\} e^{ik_t\tau_0} \quad (A.28) $$
+ 5 perms.

Note that it is important to keep track of factors of $i$ because the three-point correlator is real. The integral for the internal polynomial involving the $B$ tensor for $n \ge 0$ reads

$$ \begin{aligned} \underline{B}_{\alpha\beta\gamma}^{(2n)} &= \frac{1}{2} \int_{-\infty}^{\tau_0} d\tau a^4 \frac{B_{\alpha\beta\gamma}^{(2n)}}{a^{2n}} \frac{1}{a} (1 - ik_1\tau)(1 - ik_2\tau)k_3^2\tau e^{ik_t\tau} + 5 \text{ perms} \\ &= \frac{-B_{\alpha\beta\gamma}^{(2n)}}{2i^{3-2n}H^{3-2n}} \left\{ \frac{e_3k_3}{k_t^{2n+1}} \Gamma(2n+1, -ik_t\tau_0) + \frac{(k_1+k_2)k_3^2}{k_t^{2n}} \Gamma(2n, -ik_t\tau_0) \right. \\ & \qquad \left. + \frac{k_3^2}{k_t^{2n-1}} \Gamma(2n-1, -ik_t\tau_0) \right\} + 5 \text{ perms.} \end{aligned} \quad (A.29) $$

When expanding the $\Gamma$-function, only $\underline{B}^{(0)}$ has an infinite number of terms. Keeping only the two leading-order terms, one obtains

$$ \underline{B}_{\alpha\beta\gamma}^{(2n)} = \frac{-B_{\alpha\beta\gamma}^{(2n)}}{2H^{3-2n}k_t} \left\{ i e_3 \tau_0^{2n} - \frac{k_3}{k_t} [2n e_3 + k_t(k_1+k_2)k_3] \tau_0^{2n-1} + \mathcal{O}(\tau_0^{2n-2}) \right\} e^{ik_t \tau_0} \quad (A.30) $$
+ 5 perms.

Similarly, the integral for the internal polynomial involving the $C$ tensor for $n \ge 0$ is

$$ \begin{aligned} \underline{C}_{\alpha\beta\gamma}^{(2n)} &= \frac{1}{2} \int_{-\infty}^{\tau_0} d\tau a^4 \frac{C_{\alpha\beta\gamma}^{(2n)}}{a^{2n}} \frac{1}{a^2} k_1^2 \tau k_2^2 \tau (1 - ik_3\tau) e^{ik_t\tau} + 5 \text{ perms} \\ &= \frac{C_{\alpha\beta\gamma}^{(2n)} i^{2n+1}}{2H^{2-2n}} \left\{ -\frac{e_3 k_1 k_2}{k_t^{2n+2}} \Gamma(2n+2, -ik_t\tau_0) - \frac{k_1^2 k_2^2}{k_t^{2n+1}} \Gamma(2n+1, -ik_t\tau) \right\} \\ & \qquad + 5 \text{ perms.} \end{aligned} \quad (A.31) $$

We note that $\underline{C}$ has a finite number of terms for all $n \ge 0$. However for consistency we still keep only the first two leading-order terms, giving

$$ \underline{C}_{\alpha\beta\gamma}^{(2n)} = \frac{-C_{\alpha\beta\gamma}^{(2n)}}{2H^{2-2n}k_t} \left\{ e_3 k_1 k_2 \tau_0^{2n+1} + i k_1^2 k_2^2 \left( 1 + \frac{(2n+1)k_3}{k_t} \right) \tau_0^{2n} + \mathcal{O}(\tau_0^{2n-1}) \right\} e^{ik_t \tau_0} \quad (A.32) $$
+ 5 perms,

where the first line is an exact result. Finally for the internal polynomial involving the $D$ tensor
for $n \ge 0$, we obtain

$$ 
\begin{aligned}
\underline{D}_{\alpha\beta\gamma}^{(2n)} &= \frac{1}{2} \int_{-\infty}^{\tau_0} d\tau a^4 \frac{D_{\alpha\beta\gamma}^{(2n)}}{a^{2n}} \frac{1}{a^3} e_3^2 \tau^3 e^{ik_t \tau} + 5 \text{ perms} \quad (A.33) \\
&= \frac{D_{\alpha\beta\gamma}^{(2n)}}{2H^{1-2n}} \frac{e_3^2 i^{2n+3}}{k_t^{2n+3}} \Gamma(2n+3, -ik_t\tau_0) + 5 \text{ perms},
\end{aligned}
$$

from which we get

$$ \underline{D}_{\alpha\beta\gamma}^{(2n)} = \frac{D_{\alpha\beta\gamma}^{(2n)}}{2H^{1-2n}k_t} e_3^2 \left\{ i\tau_0^{2n+2} - \frac{(2n+2)\tau_0^{2n+1}}{k_t} + \mathcal{O}(\tau_0^{2n}) \right\} e^{ik_t\tau_0} + 5 \text{ perms.} \quad (A.34) $$

Like the previous case, $\underline{D}$ has a finite number of terms for all $n \ge 0$. We have now all the necessary ingredients to obtain the initial conditions for the three-point correlators. Combining the external polynomials in Eq. (A.20), the truncated internal polynomials in Eqs. (A.28)-(A.30)-(A.32) and (A.34), and taking into account the common overall constant, the final results are the ones given in (3.34), (3.35), (3.36) and (3.37).

**Tensorial structure.** We have omitted the tensorial structure coming from the field normalisations in (3.34), (3.35), (3.36) and (3.37) to cleanse the expressions. Combinations of $\Delta$ tensors can be included a posteriori by performing the following replacement on the correlators (no summation on repeated indices)

$$ 
\begin{aligned}
\langle \varphi_{k_1}^{\alpha} \varphi_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle_0 &\rightarrow \Delta_{\alpha\alpha} \Delta_{\beta\beta} \Delta_{\gamma\gamma} \langle \varphi_{k_1}^{\alpha} \varphi_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle_0, \\
\langle \mathbf{p}_{k_1}^{\alpha} \varphi_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle_0 &\rightarrow \Delta_{\alpha\alpha}^{-1} \Delta_{\beta\beta} \Delta_{\gamma\gamma} \langle \mathbf{p}_{k_1}^{\alpha} \varphi_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle_0, \\
\langle \mathbf{p}_{k_1}^{\alpha} \mathbf{p}_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle_0 &\rightarrow \Delta_{\alpha\alpha}^{-1} \Delta_{\beta\beta} \Delta_{\gamma\gamma}^{-1} \langle \mathbf{p}_{k_1}^{\alpha} \mathbf{p}_{k_2}^{\beta} \varphi_{k_3}^{\gamma} \rangle_0, \\
\langle \mathbf{p}_{k_1}^{\alpha} \mathbf{p}_{k_2}^{\beta} \mathbf{p}_{k_3}^{\gamma} \rangle_0 &\rightarrow \Delta_{\alpha\alpha}^{-1} \Delta_{\beta\beta}^{-1} \Delta_{\gamma\gamma}^{-1} \langle \mathbf{p}_{k_1}^{\alpha} \mathbf{p}_{k_2}^{\beta} \mathbf{p}_{k_3}^{\gamma} \rangle_0,
\end{aligned}
\quad (A.35) $$

and by redefining the third-order vertex functions as follows

$$ 
\begin{aligned}
A_{\alpha\beta\gamma}^{(2n)} &\rightarrow \Delta_{\alpha\alpha} \Delta_{\beta\beta} \Delta_{\gamma\gamma} A_{\alpha\beta\gamma}^{(2n)}, \\
B_{\alpha\beta\gamma}^{(2n)} &\rightarrow \Delta_{\alpha\alpha} \Delta_{\beta\beta} \Delta_{\gamma\gamma}^{-1} B_{\alpha\beta\gamma}^{(2n)}, \\
C_{\alpha\beta\gamma}^{(2n)} &\rightarrow \Delta_{\alpha\alpha}^{-1} \Delta_{\beta\beta}^{-1} \Delta_{\gamma\gamma} C_{\alpha\beta\gamma}^{(2n)}, \\
D_{\alpha\beta\gamma}^{(2n)} &\rightarrow \Delta_{\alpha\alpha}^{-1} \Delta_{\beta\beta}^{-1} \Delta_{\gamma\gamma}^{-1} D_{\alpha\beta\gamma}^{(2n)}.
\end{aligned}
\quad (A.36) $$

for all $n \ge 0$. These manipulations are necessary if one or several fields are not canonically normalised.
## B Light Limit of the Dynamics

In this appendix, we first derive a simple formula for the power spectrum in the light field and weak mixing regime of the phase diagram. This reproduces the observed enhancement in Figure 4. We then comment on the effect of the non-linear conversion from $\pi$ to $\zeta$ on the bispectrum.

**Linear Dynamics.** In Section 4.2, we obtained exact results for the Goldstone boson power spectrum in the entire phase space $(m/H, \rho/H)$. Our result shows a characteristic enhancement of $\Delta_{\zeta}^2$ in the light limit $m/H \ll 1$ and at weak mixing $\rho/H \ll 1$. Here, we explicitly show that this enhancement is a consequence of the light $\sigma$ field not decaying fast enough on super-horizon scales, therefore leaving an imprint in $\Delta_{\zeta}^2$ which strongly depends on the number of super-horizon e-folds. A very similar analysis has been applied to the spinning case in [201] so that we will closely follow their derivation. We set $c_s = 1$ for simplicity throughout this appendix.

Our starting point is the quadratic Lagrangian in (4.6) that leads to the following coupled equations of motion

$$ \begin{aligned} \ddot{\pi}_c + 3H\dot{\pi}_c + \frac{k^2}{a^2}\pi_c &= -\rho[3H\sigma + \dot{\sigma}], \\ \ddot{\sigma} + 3H\dot{\sigma} + \left(\frac{k^2}{a^2} + m^2\right)\sigma &= \rho\dot{\pi}_c. \end{aligned} \qquad (B.1) $$

We are interested in obtaining the power spectrum at late times so we want to find the late-time behaviour of the mode functions. Hence, we look for a solution of the form $\pi_c = \pi_0\tau^{\Delta}$ and $\sigma = \sigma_0\tau^{\Delta}$. Plugging this ansatz in the system and neglecting the gradient term leads to

$$ \begin{pmatrix} \Delta(\Delta - 3) & \frac{\rho}{H}(3 - \Delta) \\ \frac{\rho}{H}\Delta & \Delta(\Delta - 3) + \frac{m^2}{H^2} \end{pmatrix} \begin{pmatrix} \pi_0 \\ \sigma_0 \end{pmatrix} = 0. \qquad (B.2) $$

The requirement of a vanishing determinant gives us the following scaling dimensions

$$ \Delta_{-}^{\pi_c} = 0, \quad \Delta_{+}^{\pi_c} = 3, \quad \Delta_{\pm}^{\sigma} = \frac{3}{2} \pm \sqrt{\frac{9}{4} - \frac{m_{\text{eff}}^2}{H^2}}, \qquad (B.3) $$

with $m_{\text{eff}}^2 = m^2 + \rho^2$. Plugging back in the matrix equation the found solutions and trying to solve for $\pi_0$ and $\sigma_0$ gives us constraints satisfied by the mode amplitudes:

* $\Delta = 0$ enforces $\sigma_0 = 0$ but leaves $\pi_0$ underdetermined. In fact, the precise amplitude of this massless mode is fixed by the quantisation condition in the far past.

* $\Delta = \frac{3}{2} \pm \sqrt{9/4 - m_{\text{eff}}^2/H^2}$ gives $\frac{\pi_0}{\sigma_0} = \frac{\rho}{H} \frac{\Delta_{-}^{\sigma}}{\Delta_{-}^{\sigma} - \Delta_{+}^{\sigma}}$ and $\frac{\sigma_0}{\pi_0} = \frac{H}{\rho} \Delta_{-}^{\sigma}$. Both solutions are equivalent. This solution is the manifestation of the exponential dilution of $p_{\pi_c}$ with time as the Universe inflates. Indeed at late times neglecting gradient terms, one can show that the conjugate momentum to the field $\pi_c$ satisfies

$$ \dot{p}_{\pi_c} + 3H p_{\pi_c} = 0, \qquad (B.4) $$

so that it decays as $p_{\pi_c} \propto 1/a^3$. Close to the infinite future boundary, we then have
$\dot{\pi}_c \approx -\rho\sigma$. If the leading behaviour of the massive field is $\sigma \sim \tau^{\Delta_{\pm}^{\sigma}}$, then after integration one finds $\pi_c \sim \text{cst} + \frac{\rho}{H}\frac{\sigma}{\Delta_{\pm}^{\sigma}}$ which gives us the found constraint on the mode amplitudes. The constant is found by imposing the quantisation condition.

* $\Delta = 3$ enforces $\frac{\rho}{H}\pi_0 + \frac{m^2}{H^2}\sigma_0 = 0$. In the $M = \frac{m^2}{\rho} \gg H$ limit, we obtain $\sigma_0 = 0$ if the field is massive. This is expected because a massive field in de-Sitter does not contain any $\Delta = 3$ mode. In the $M = \frac{m^2}{\rho} \ll H$, we obtain $\pi_0 = 0$. This is not surprising as the $\sigma$ field can be integrated out, leading to a modified dispersion relation for $\pi_c$. One can explicitly check that the exact form of the $\pi_c$ mode function$^{55}$ in this regime does not contain a $\Delta = 3$ mode.

Imposing the previous constraints on the mode amplitudes, the general solution for the system at late times is

$$ \begin{pmatrix} \pi_c \\ \sigma \end{pmatrix} = \pi_0^- \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \pi_0^+ (-k\tau)^3 \begin{pmatrix} \frac{m^2}{H^2} \\ -\frac{\rho}{H} \end{pmatrix} + \sigma_0^- (-k\tau)^{\Delta_-} \begin{pmatrix} \frac{\rho}{H} \Delta_+^\sigma \\ \frac{m_{\text{eff}}^2}{H^2} \end{pmatrix} + \sigma_0^+ (-k\tau)^{\Delta_+^\sigma} \begin{pmatrix} \frac{\rho}{H} \\ \Delta_+^\sigma \end{pmatrix}, \quad (B.5) $$

where $\pi_0^-, \pi_0^+, \sigma_0^-$ and $\sigma_0^+$ are integration constants. When $\sigma$ is effectively light $m_{\text{eff}} \ll H$ but not strictly speaking massless, its $\sigma_0^-$ mode decays very slowly on super-horizon scales with scaling dimension $\Delta_-^\sigma \approx \frac{m_{\text{eff}}^2}{3H^2}$. Only considering the two lightest modes, the solution for the system is

$$ \begin{pmatrix} \pi_c \\ \sigma \end{pmatrix} \approx \pi_0^- \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \sigma_0^- (-k\tau)^{\frac{m_{\text{eff}}^2}{3H^2}} \begin{pmatrix} \frac{3\rho}{H} \\ \frac{m_{\text{eff}}^2}{H^2} \end{pmatrix}. \quad (B.6) $$

The $\sigma_0^-$ mode has a component in the $\sigma$ direction that is parametrically smaller than its component in the $\pi_c$ direction in the regime of interest. We now impose that at horizon crossing, the $\sigma_0^-$ mode is dominated by the component in the $\sigma$ direction

$$ \pi_c|_{\star} = 0, \quad \sigma|_{\star} = A_0 (-k\tau)^{\frac{m_{\text{eff}}^2}{3H^2}}, \quad (B.7) $$

where $\star$ means that the quantity is evaluated at horizon crossing $k = aH$. This condition enforces that we chose $A_0 = \frac{m_{\text{eff}}^2}{H^2}\sigma_0^-$ and $\pi_0^- = -\frac{3\rho}{H}\sigma_0^-$ in order to have an exact cancellation between the large $\pi_c$ modes at horizon crossing, i.e. this amounts to only study the sourcing by $\sigma$ on $\pi_c$ on super-horizon scales. The solution for $\pi_c$ is then determined up to an overall constant

$$ \pi_c = -A_0 \frac{3\rho H}{m_{\text{eff}}^2} \left(1 - e^{-\frac{m_{\text{eff}}^2}{3H^2} \Delta N}\right), \quad (B.8) $$

where $\Delta N = N - N_{\star}$ is the number of e-folds elapsed since horizon crossing. The power spectrum correction is found to be

$$ \frac{\Delta_{\zeta}^2}{\Delta_{\zeta,0}^2} - 1 \propto \left[ \frac{3\rho H}{m_{\text{eff}}^2} \left( 1 - e^{-\frac{m_{\text{eff}}^2}{3H^2} \Delta N} \right) \right]^2. \quad (B.9) $$

<footnote>$^{55}$The leading term of $z^{5/4}H_{5/4}^{(1)}(z^2)$ as $z \to 0$ is $\propto z^4$.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 13: Correction to the power spectrum (normalised to its single-field value for $\rho = 0$) in the effectively light regime $m_{\text{eff}} \ll H$ for $\Delta N = 60$, $\Delta N = 40$, and $\Delta N = 20$, where $\Delta N$ is the number of super-horizon e-folds. The solid and dashed lines correspond to the exact numerical results and the analytical approximation in (B.9), respectively. The amplitude of the approximation (B.9) was determined to match the exact result at $\rho/H = 10^{-2}$. The mismatch between the solid and dashed lines for $\rho/H \gtrsim 0.5$ indicates that the field $\sigma$ is effectively massive enough so that the approximation breaks down.</figcaption></figure>

At weak coupling $\rho \ll H$ when the light mode does not have enough time to decay, this formula reproduces well the leading correction found in perturbative calculations [202]

$$ \frac{\Delta_{\zeta}^2}{\Delta_{\zeta,0}^2} - 1 \propto \frac{\rho^2}{H^2} \Delta N^2, \qquad (B.10) $$

where the precise coefficient has been computed in [154, 203]. However when the effective mass becomes of order $m_{\text{eff}} \approx H$, this analytical approximation breaks down and does not capture the modified dispersion relation behaviour. This is illustrated in Figure 13 where we show a slice of the phase diagram at fixed bare mass $m/H = 0.1$ with the exact and approximated power spectrum for different values of $\Delta N$. A notable feature is that the amplitude and the location of the peak is well reproduced by the analytical approximation (B.9). This qualitatively explains the boost of the power spectrum in the case of a light $\sigma$ field. Note that the finite duration of inflation induces an intrinsic violation of scale-invariance, and that a complete prediction even just for the power spectrum would require to describe the reheating epoch.

**Non-linear conversion and the bispectrum.** The non-linear conversion from $\pi$ to $\zeta$ en- genders additional terms in the bispectrum that should be added to the “intrinsic” bispectrum generated by cubic interactions. These terms are displayed in Eq. (4.12) for any theory in the super-horizon limit. In the decoupling limit, only the second term in the first line contributes. Specifying to equilateral configurations $k_1 = k_2 = k_3$ and in the theory (4.6), the size of non- Gaussianities needs to be corrected by

$$ f_{\text{NL}} \rightarrow f_{\text{NL}} - \frac{5}{3} \frac{\rho}{H} \frac{\sum_{\mathcal{R}}^{\pi_c \sigma}}{\sum_{\pi_c \pi_c}}. \qquad (B.11) $$
<figure><img src="image_16.png" /><figcaption>Figure 14: Correction to the “intrinsic” $f_{\text{NL}}$ parameter coming from the non-linear conversion from $\pi$ to $\zeta$, see Eq. (4.12), in the limit $\dot{H} \to 0$. We have set $N = 60$ e-folds of super-horizon evolution. The dashed line corresponds to $M = m^2/\rho = H$, see Section 4.2.2 for more details.</figcaption></figure>

The exact correction is displayed in Figure 14 in the phase space $(m/H, \rho/H)$. As anticipated, this contribution only affects the light-field regime, where $\sigma$ decays slowly on super-horizon scales. In the parameter space under display, its amplitude is comparable to the expected gravitational floor ($\sim 10^{-2}$), and therefore negligible for current and near-future cosmological surveys. In the effectively light field limit, one expects the bispectrum to resemble a local shape in its squeezed limit. However, predicting the observable part of it would require subtracting the artefact coming from a rescaling of the background by a long wavelength perturbation, including the usual consistency relation that is proportional to the tilt of the power spectrum. All this depends on the arbitrary choice of the duration of inflation, and is also degenerate with a description of reheating. A full study of these effects is beyond the scope of this paper.
C Perturbativity Bounds

In this appendix, for the reader's convenience, we collect all the perturbativity bounds on the coupling constants of the theory (4.6) found in Section 5, both in the weak and strong mixing regimes.<fnref n="56" /> We recall that these bounds follow from imposing that the strong coupling scales of irrelevant operators be higher than the Hubble scale, or imposing standard perturbativity criteria on relevant/marginal operators. We have coloured in red the most stringent bounds.

Weak Mixing Regime

<table><thead><tr><th>Operators<br>(unitary<br>gauge)</th><th>δg<sup>00</sup>σ</th><th>δg<sup>00</sup>σ<sup>2</sup></th><th>(δg<sup>00</sup>)<sup>2</sup> σ</th><th>σ<sup>3</sup></th></tr></thead><tbody><tr><td>Interactions</td><td>π̇<sub>c</sub>σ, (∂<sub>i</sub>π<sub>c</sub>)<sup>2</sup>σ</td><td>π̇<sub>c</sub>σ<sup>2</sup>, π̇<sub>c</sub><sup>2</sup>σ<sup>2</sup>,<br>(∂<sub>i</sub>π<sub>c</sub>)<sup>2</sup>σ<sup>2</sup></td><td>π̇<sub>c</sub><sup>2</sup>σ, π̇<sub>c</sub>(∂<sub>μ</sub>π<sub>c</sub>)<sup>2</sup>σ, π̇<sub>c</sub><sup>4</sup>σ,<br>π̇<sub>c</sub><sup>2</sup>(∂<sub>i</sub>π<sub>c</sub>)<sup>2</sup>σ, (∂<sub>i</sub>π<sub>c</sub>)<sup>4</sup>σ</td><td>σ<sup>3</sup></td></tr><tr><td>Λ<sub>*</sub></td><td>c<sub>s</sub><sup>3/2</sup> f<sub>π</sub><sup>2</sup> / ρ</td><td>c<sub>s</sub><sup>3/4</sup> f<sub>π</sub> / √α</td><td>c<sub>s</sub><sup>-1/2</sup> f<sub>π</sub><sup>2</sup> / ρ</td><td>×</td></tr><tr><td>H ≲ Λ<sub>*</sub></td><td>ρ/H ≲ c<sub>s</sub><sup>3/2</sup> / (2πΔ<sub>ζ</sub>)</td><td>α̃ ≲ c<sub>s</sub><sup>-1/8</sup> / (4πΔ<sub>ζ</sub>)</td><td>ρ̃/H ≲ c<sub>s</sub><sup>-1/2</sup> / (2πΔ<sub>ζ</sub>)</td><td>×</td></tr><tr><td>Perturbativity<br>of relevant or<br>marginal<br>operators</td><td>ρ/H ≲ c<sub>s</sub><sup>-1/2</sup></td><td>α̃ ≲ c<sub>s</sub><sup>-1/2</sup> / (2πΔ<sub>ζ</sub><sup>1/2</sup>)</td><td>×</td><td>μ/H ≲ 1</td></tr></tbody></table>

Strong Mixing Regime

<table><thead><tr><th>Operators<br>(unitary<br>gauge)</th><th>δg<sup>00</sup>σ</th><th>δg<sup>00</sup>σ<sup>2</sup></th><th>(δg<sup>00</sup>)<sup>2</sup> σ</th><th>σ<sup>3</sup></th></tr></thead><tbody><tr><td>Interactions</td><td>π̇<sub>c</sub>σ, (θ̃<sub>i</sub>π̃<sub>c</sub>)<sup>2</sup>σ</td><td>π̇<sub>c</sub>σ<sup>2</sup>, π̇<sub>c</sub><sup>2</sup>σ<sup>2</sup>, (θ̃<sub>i</sub>π̃<sub>c</sub>)<sup>2</sup>σ<sup>2</sup></td><td>π̇<sub>c</sub><sup>2</sup>σ, π̇<sub>c</sub>(∂<sub>μ</sub>π̃<sub>c</sub>)<sup>2</sup>σ, π̇<sub>c</sub><sup>4</sup>σ,<br>π̇<sub>c</sub><sup>2</sup>(θ̃<sub>i</sub>π̃<sub>c</sub>)<sup>2</sup>σ, (θ̃<sub>i</sub>π̃<sub>c</sub>)<sup>4</sup>σ</td><td>σ<sup>3</sup></td></tr><tr><td>Λ<sub>*</sub></td><td>c<sub>s</sub> (f<sub>π</sub><sup>8</sup> / ρ<sup>5</sup>)<sup>1/3</sup></td><td>(ρ/f<sub>π</sub><sup>4</sup> / α<sup>2</sup>)<sup>1/5</sup>, (c<sub>s</sub><sup>2</sup>f<sub>π</sub><sup>4</sup> / α<sup>2</sup>ρ)<sup>1/3</sup></td><td>(ρ<sup>3</sup>f<sub>π</sub><sup>8</sup> / c<sub>s</sub>ρ<sup>4</sup>)<sup>1/7</sup></td><td>×</td></tr><tr><td>H ≲ Λ<sub>*</sub></td><td>ρ/H ≲ c<sub>s</sub><sup>3/4</sup> / (κ<sup>1/2</sup>Δ<sub>ζ</sub>)</td><td>α̃ ≲ κ<sup>3/4</sup> / (2c<sub>s</sub><sup>3/8</sup>Δ<sub>ζ</sub><sup>3/2</sup>),<br>α̃ ≲ κ<sup>1/2</sup> / (2c<sub>s</sub><sup>1/4</sup>Δ<sub>ζ</sub>)</td><td>ρ̃/H ≲ c<sub>s</sub><sup>1/2</sup> / (κΔ<sub>ζ</sub><sup>2</sup>)</td><td>×</td></tr><tr><td>Perturbativity<br>of relevant or<br>marginal<br>operators</td><td>×</td><td>α̃ ≲ κ<sup>1/2</sup> / (2c<sub>s</sub><sup>1/4</sup>Δ<sub>ζ</sub>)</td><td>×</td><td>μ/H ≲ c<sub>s</sub><sup>-3/4</sup> (ρ/H)<sup>3/4</sup></td></tr></tbody></table>

<footnote n="56">⁵⁶We define $κ = 2Γ(5/4)^2/π^3 \approx 0.053$, see e.g. (4.29).</footnote>
References

[1] A. A. Starobinsky, “Spectrum of relict gravitational radiation and the early state of the universe,” *JETP Lett.* **30** (1979) 682–685.

[2] A. A. Starobinsky, “A New Type of Isotropic Cosmological Models Without Singularity,” *Phys. Lett. B* **91** (1980) 99–102.

[3] A. H. Guth, “Inflationary universe: A possible solution to the horizon and flatness problems,” *Phys. Rev. D* **23** (Jan, 1981) 347–356. https://link.aps.org/doi/10.1103/PhysRevD.23.347.

[4] A. H. Guth and S.-Y. Pi, “Fluctuations in the new inflationary universe,” *Phys. Rev. Lett.* **49** (Oct, 1982) 1110–1113. https://link.aps.org/doi/10.1103/PhysRevLett.49.1110.

[5] A. D. Linde, “A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems,” *Phys. Lett. B* **108** (1982) 389–393.

[6] A. D. Linde, “Chaotic Inflation,” *Phys. Lett. B* **129** (1983) 177–181.

[7] A. D. Linde, “ETERNAL CHAOTIC INFLATION,” *Mod. Phys. Lett. A1* (1986) 81.

[8] A. D. Linde, “Eternally Existing Selfreproducing Chaotic Inflationary Universe,” *Phys. Lett. B* **175** (1986) 395–400.

[9] V. F. Mukhanov and G. V. Chibisov, “Quantum Fluctuations and a Nonsingular Universe,” *JETP Lett.* **33** (1981) 532–535.

[10] S. W. Hawking, “The Development of Irregularities in a Single Bubble Inflationary Universe,” *Phys. Lett. B* **115** (1982) 295.

[11] A. A. Starobinsky, “Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations,” *Phys. Lett. B* **117** (1982) 175–178.

[12] A. H. Guth and S. Y. Pi, “Fluctuations in the New Inflationary Universe,” *Phys. Rev. Lett.* **49** (1982) 1110–1113.

[13] A. D. Linde, “Scalar Field Fluctuations in Expanding Universe and the New Inflationary Universe Scenario,” *Phys. Lett. B* **116** (1982) 335–339.

[14] A. Vilenkin and L. H. Ford, “Gravitational Effects upon Cosmological Phase Transitions,” *Phys. Rev. D* **26** (1982) 1231.

[15] V. F. Mukhanov, “Gravitational Instability of the Universe Filled with a Scalar Field,” *JETP Lett.* **41** (1985) 493–496.

[16] M. Sasaki, “Large Scale Quantum Fluctuations in the Inflationary Universe,” *Prog. Theor. Phys.* **76** (1986) 1036.

[17] J. M. Maldacena, “Non-Gaussian features of primordial fluctuations in single field inflationary models,” *JHEP* **05** (2003) 013, arXiv:astro-ph/0210603.

[18] P. Creminelli, “On non-Gaussianities in single-field inflation,” *JCAP* **10** (2003) 003, arXiv:astro-ph/0306122.

[19] X. Chen, “Primordial Non-Gaussianities from Inflation Models,” *Adv. Astron.* **2010** (2010) 638979, arXiv:1002.1416 [astro-ph.CO].

[20] S. Renaux-Petel, “Primordial non-Gaussianities after Planck 2015: an introductory review,” *Comptes Rendus Physique* **16** (2015) 969–985, arXiv:1508.06740 [astro-ph.CO].
[21] P. D. Meerburg et al., "Primordial Non-Gaussianity," *Bull. Am. Astron. Soc.* **51** no. 3, (2019) 107, arXiv:1903.04409 [astro-ph.CO].

[22] K. Abazajian et al., "CMB-S4 Science Case, Reference Design, and Project Plan," arXiv:1907.04473 [astro-ph.IM].

[23] M. Alvarez et al., "Testing Inflation with Large Scale Structure: Connecting Hopes with Reality," arXiv:1412.4671 [astro-ph.CO].

[24] A. Achúcarro et al., "Inflation: Theory and Observations," arXiv:2203.08128 [astro-ph.CO].

[25] **Snowmass Cosmic Frontier 5 Topical Group** Collaboration, K. S. Karkare, A. M. Dizgah, G. K. Keating, P. Breysse, and D. T. Chung, "Snowmass 2021 Cosmic Frontier White Paper: Cosmology with Millimeter-Wave Line Intensity Mapping," in *2022 Snowmass Summer Study*, 3, 2022. arXiv:2203.07258 [astro-ph.CO].

[26] S. Ferraro, N. Sailer, A. Slosar, and M. White, "Snowmass2021 Cosmic Frontier White Paper: Cosmology and Fundamental Physics from the three-dimensional Large Scale Structure," arXiv:2203.07506 [astro-ph.CO].

[27] D. Baumann and D. Green, "Equilateral Non-Gaussianity and New Physics on the Horizon," *JCAP* **09** (2011) 014, arXiv:1102.5343 [hep-th].

[28] **Cosmic Visions 21 cm** Collaboration, R. Ansari et al., "Inflation and Early Dark Energy with a Stage II Hydrogen Intensity Mapping experiment," arXiv:1810.09572 [astro-ph.CO].

[29] J. B. Muñoz, Y. Ali-Haïmoud, and M. Kamionkowski, "Primordial non-gaussianity from the bispectrum of 21-cm fluctuations in the dark ages," *Phys. Rev. D* **92** no. 8, (2015) 083508, arXiv:1506.04152 [astro-ph.CO].

[30] X. Chen, P. D. Meerburg, and M. Münchmeyer, "The Future of Primordial Features with 21 cm Tomography," *JCAP* **09** (2016) 023, arXiv:1605.09364 [astro-ph.CO].

[31] **Snowmass 2021 Cosmic Frontier 5 Topical Group** Collaboration, A. Liu, L. Newburgh, B. Saliwanchik, and A. Slosar, "Snowmass2021 Cosmic Frontier White Paper: 21cm Radiation as a Probe of Physics Across Cosmic Ages," arXiv:2203.07864 [astro-ph.CO].

[32] **Planck** Collaboration, Y. Akrami et al., "Planck 2018 results. X. Constraints on inflation," *Astron. Astrophys.* **641** (2020) A10, arXiv:1807.06211 [astro-ph.CO].

[33] **BICEP, Keck** Collaboration, P. A. R. Ade et al., "Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season," *Phys. Rev. Lett.* **127** no. 15, (2021) 151301, arXiv:2110.00483 [astro-ph.CO].

[34] X. Chen and Y. Wang, "Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation," *Phys. Rev. D* **81** (2010) 063511, arXiv:0909.0496 [astro-ph.CO].

[35] X. Chen and Y. Wang, "Quasi-Single Field Inflation and Non-Gaussianities," *JCAP* **04** (2010) 027, arXiv:0911.3380 [hep-th].

[36] T. Noumi, M. Yamaguchi, and D. Yokoyama, "Effective field theory approach to quasi-single field inflation and effects of heavy fields," *JHEP* **06** (2013) 051, arXiv:1211.1624 [hep-th].

[37] N. Arkani-Hamed and J. Maldacena, "Cosmological Collider Physics," arXiv:1503.08043 [hep-th].
[38] X. Chen, M. H. Namjoo, and Y. Wang, “Quantum Primordial Standard Clocks,” JCAP 02 (2016) 013, arXiv:1509.03930 [astro-ph.CO].

[39] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Loop Corrections to Standard Model Fields in Inflation,” JHEP 08 (2016) 051, arXiv:1604.07841 [hep-th].

[40] H. Lee, D. Baumann, and G. L. Pimentel, “Non-Gaussianity as a Particle Detector,” JHEP 12 (2016) 040, arXiv:1607.03735 [hep-th].

[41] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Standard Model Background of the Cosmological Collider,” Phys. Rev. Lett. 118 no. 26, (2017) 261302, arXiv:1610.06597 [hep-th].

[42] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Schwinger-Keldysh Diagrammatics for Primordial Perturbations,” JCAP 12 (2017) 006, arXiv:1703.10166 [hep-th].

[43] H. An, M. McAneny, A. K. Ridgway, and M. B. Wise, “Quasi Single Field Inflation in the non-perturbative regime,” JHEP 06 (2018) 105, arXiv:1706.09971 [hep-ph].

[44] A. V. Iyer, S. Pi, Y. Wang, Z. Wang, and S. Zhou, “Strongly Coupled Quasi-Single Field Inflation,” JCAP 01 (2018) 041, arXiv:1710.03054 [hep-th].

[45] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Neutrino Signatures in Primordial Non-Gaussianities,” JHEP 09 (2018) 022, arXiv:1805.02656 [hep-ph].

[46] S. Lu, Y. Wang, and Z.-Z. Xianyu, “A Cosmological Higgs Collider,” JHEP 02 (2020) 011, arXiv:1907.07390 [hep-th].

[47] L.-T. Wang and Z.-Z. Xianyu, “In Search of Large Signals at the Cosmological Collider,” JHEP 02 (2020) 044, arXiv:1910.12876 [hep-ph].

[48] S. Alexander, S. J. Gates, L. Jenks, K. Koutrolikos, and E. McDonough, “Higher Spin Supersymmetry at the Cosmological Collider: Sculpting SUSY Rilles in the CMB,” JHEP 10 (2019) 156, arXiv:1907.05829 [hep-th].

[49] A. Hook, J. Huang, and D. Racco, “Minimal signatures of the Standard Model in non-Gaussianities,” Phys. Rev. D 101 no. 2, (2020) 023519, arXiv:1908.00019 [hep-ph].

[50] T. Liu, X. Tong, Y. Wang, and Z.-Z. Xianyu, “Probing P and CP Violations on the Cosmological Collider,” JHEP 04 (2020) 189, arXiv:1909.01819 [hep-ph].

[51] L.-T. Wang and Z.-Z. Xianyu, “Gauge Boson Signals at the Cosmological Collider,” JHEP 11 (2020) 082, arXiv:2004.02887 [hep-ph].

[52] L.-T. Wang, Z.-Z. Xianyu, and Y.-M. Zhong, “Precision calculation of inflation correlators at one loop,” JHEP 02 (2022) 085, arXiv:2109.14635 [hep-ph].

[53] N. Maru and A. Okawa, “Non-Gaussianity from X, Y gauge bosons in Cosmological Collider Physics,” arXiv:2101.10634 [hep-ph].

[54] Y. Cui and Z.-Z. Xianyu, “Probing Leptogenesis with the Cosmological Collider,” Phys. Rev. Lett. 129 no. 11, (2022) 111301, arXiv:2112.10793 [hep-ph].

[55] L. Pinol, S. Aoki, S. Renaux-Petel, and M. Yamaguchi, “Inflationary flavor oscillations and the cosmic spectroscopy,” Phys. Rev. D 107 no. 2, (2023) L021301, arXiv:2112.05710 [hep-th].

[56] Q. Lu, M. Reece, and Z.-Z. Xianyu, “Missing scalars at the cosmological collider,” JHEP 12 (2021) 098, arXiv:2108.11385 [hep-ph].
[57] X. Tong, Y. Wang, and Y. Zhu, “Cutting rule for cosmological collider signals: a bulk evolution perspective,” JHEP 03 (2022) 181, arXiv:2112.03448 [hep-th].

[58] X. Tong and Z.-Z. Xianyu, “Large spin-2 signals at the cosmological collider,” JHEP 10 (2022) 194, arXiv:2203.06349 [hep-ph].

[59] Z. Qin and Z.-Z. Xianyu, “Phase information in cosmological collider signals,” JHEP 10 (2022) 192, arXiv:2205.01692 [hep-th].

[60] S. Jazayeri and S. Renaux-Petel, “Cosmological bootstrap in slow motion,” JHEP 12 (2022) 137, arXiv:2205.10340 [hep-th].

[61] G. L. Pimentel and D.-G. Wang, “Boostless cosmological collider bootstrap,” JHEP 10 (2022) 177, arXiv:2205.00013 [hep-th].

[62] D. Werth, L. Pinol, and S. Renaux-Petel, “Cosmological Flow of Primordial Correlators,” arXiv:2302.00655 [hep-th].

[63] X. Tong, Y. Wang, C. Zhang, and Y. Zhu, “BCS in the Sky: Signatures of Inflationary Fermion Condensation,” arXiv:2304.09428 [hep-th].

[64] S. Jazayeri, S. Renaux-Petel, and D. Werth, “Shapes of the Cosmological Low-Speed Collider,” arXiv:2307.01751 [hep-th].

[65] S. Jazayeri, S. Renaux-Petel, X. Tong, D. Werth, and Y. Zhu, “Parity Violation from Emergent Non-Locality During Inflation,” arXiv:2308.11315 [hep-th].

[66] E. Sefusatti, J. R. Fergusson, X. Chen, and E. P. S. Shellard, “Effects and Detectability of Quasi-Single Field Inflation in the Large-Scale Structure and Cosmic Microwave Background,” JCAP 08 (2012) 033, arXiv:1204.6318 [astro-ph.CO].

[67] J. Norena, L. Verde, G. Barenboim, and C. Bosch, “Prospects for constraining the shape of non-Gaussianity with the scale-dependent bias,” JCAP 08 (2012) 019, arXiv:1204.6324 [astro-ph.CO].

[68] P. D. Meerburg, M. Münchmeyer, J. B. Muñoz, and X. Chen, “Prospects for Cosmological Collider Physics,” JCAP 03 050, arXiv:1610.06559 [astro-ph.CO].

[69] A. Moradinezhad Dizgah and C. Dvorkin, “Scale-Dependent Galaxy Bias from Massive Particles with Spin during Inflation,” JCAP 01 (2018) 010, arXiv:1708.06473 [astro-ph.CO].

[70] N. Bartolo, A. Kehagias, M. Liguori, A. Riotto, M. Shiraishi, and V. Tansella, “Detecting higher spin fields through statistical anisotropy in the CMB and galaxy power spectra,” Phys. Rev. D 97 no. 2, (2018) 023503, arXiv:1709.05695 [astro-ph.CO].

[71] A. Moradinezhad Dizgah, H. Lee, J. B. Muñoz, and C. Dvorkin, “Galaxy Bispectrum from Massive Spinning Particles,” JCAP 05 (2018) 013, arXiv:1801.07265 [astro-ph.CO].

[72] K. Kogai, K. Akitsu, F. Schmidt, and Y. Urakawa, “Galaxy imaging surveys as spin-sensitive detector for cosmological colliders,” JCAP 03 (2021) 060, arXiv:2009.05517 [astro-ph.CO].

[73] N. Kundu, A. Shukla, and S. P. Trivedi, “Constraints from Conformal Symmetry on the Three Point Scalar Correlator in Inflation,” JHEP 04 (2015) 061, arXiv:1410.2606 [hep-th].

[74] D. Anninos, T. Anous, D. Z. Freedman, and G. Konstantinidis, “Late-time Structure of the Bunch-Davies De Sitter Wavefunction,” JCAP 11 (2015) 048, arXiv:1406.5490 [hep-th].

[75] N. Kundu, A. Shukla, and S. P. Trivedi, “Ward Identities for Scale and Special Conformal Transformations in Inflation,” JHEP 01 (2016) 046, arXiv:1507.06017 [hep-th].
[76] A. Shukla, S. P. Trivedi, and V. Vishal, “Symmetry constraints in inflation, α-vacua, and the three point function,” JHEP 12 (2016) 102, arXiv:1607.08636 [hep-th].

[77] P. Benincasa, “Wavefunctionals/S-matrix techniques in de Sitter,” in 21st Hellenic School and Workshops on Elementary Particle Physics and Gravity. 3, 2022. arXiv:2203.16378 [hep-th].

[78] P. Benincasa, “Amplitudes meet Cosmology: A (Scalar) Primer,” arXiv:2203.15330 [hep-th].

[79] P. Benincasa and W. J. T. Bobadilla, “Physical representations for scattering amplitudes and the wavefunction of the universe,” SciPost Phys. 12 no. 6, (2022) 192, arXiv:2112.09028 [hep-th].

[80] P. Benincasa, A. J. McLeod, and C. Vergu, “Steinmann Relations and the Wavefunction of the Universe,” Phys. Rev. D 102 (2020) 125004, arXiv:2009.03047 [hep-th].

[81] P. Benincasa, “From the flat-space S-matrix to the Wavefunction of the Universe,” arXiv:1811.02515 [hep-th].

[82] N. Arkani-Hamed, P. Benincasa, and A. Postnikov, “Cosmological Polytopes and the Wavefunction of the Universe,” arXiv:1709.02813 [hep-th].

[83] E. Pajer, “Building a Boostless Bootstrap for the Bispectrum,” JCAP 01 (2021) 023, arXiv:2010.12818 [hep-th].

[84] J. Bonifacio, E. Pajer, and D.-G. Wang, “From amplitudes to contact cosmological correlators,” JHEP 10 (2021) 001, arXiv:2106.15468 [hep-th].

[85] E. Pajer, D. Stefanyszyn, and J. Supeł, “The Boostless Bootstrap: Amplitudes without Lorentz boosts,” JHEP 12 (2020) 198, arXiv:2007.00027 [hep-th]. [Erratum: JHEP 04, 023 (2022)].

[86] J. M. Maldacena and G. L. Pimentel, “On graviton non-Gaussianities during inflation,” JHEP 09 (2011) 045, arXiv:1104.2846 [hep-th].

[87] H. Goodhew, S. Jazayeri, and E. Pajer, “The Cosmological Optical Theorem,” JCAP 04 (2021) 021, arXiv:2009.02898 [hep-th].

[88] S. Melville and E. Pajer, “Cosmological Cutting Rules,” JHEP 05 (2021) 249, arXiv:2103.09832 [hep-th].

[89] H. Goodhew, S. Jazayeri, M. H. Gordon Lee, and E. Pajer, “Cutting cosmological correlators,” JCAP 08 (2021) 003, arXiv:2104.06587 [hep-th].

[90] S. Jazayeri, E. Pajer, and D. Stefanyszyn, “From locality and unitarity to cosmological correlators,” JHEP 10 (2021) 065, arXiv:2103.08649 [hep-th].

[91] L. Di Pietro, V. Gorbenko, and S. Komatsu, “Analyticity and unitarity for cosmological correlators,” JHEP 03 (2022) 023, arXiv:2108.01695 [hep-th].

[92] D. Meltzer, “The inflationary wavefunction from analyticity and factorization,” JCAP 12 no. 12, (2021) 018, arXiv:2107.10266 [hep-th].

[93] M. Hogervorst, J. a. Penedones, and K. S. Vaziri, “Towards the non-perturbative cosmological bootstrap,” JHEP 02 (2023) 162, arXiv:2107.13871 [hep-th].

[94] D. Baumann, W.-M. Chen, C. Duaso Pueyo, A. Joyce, H. Lee, and G. L. Pimentel, “Linking the singularities of cosmological correlators,” JHEP 09 (2022) 010, arXiv:2106.05294 [hep-th].

[95] T. Heckelbacher, I. Sachs, E. Skvortsov, and P. Vanhove, “Analytical evaluation of cosmological correlation functions,” JHEP 08 (2022) 139, arXiv:2204.07217 [hep-th].
[96] J. Bonifacio, H. Goodhew, A. Joyce, E. Pajer, and D. Stefanyszyn, “The graviton four-point function in de Sitter space,” JHEP 06 (2023) 212, arXiv:2212.07370 [hep-th].

[97] S. A. Salcedo, M. H. G. Lee, S. Melville, and E. Pajer, “The Analytic Wavefunction,” JHEP 06 (2023) 020, arXiv:2212.08009 [hep-th].

[98] S. Agui-Salcedo and S. Melville, “The Cosmological Tree Theorem,” arXiv:2308.00680 [hep-th].

[99] C. Duaso Pueyo and E. Pajer, “A Cosmological Bootstrap for Resonant Non-Gaussianity,” arXiv:2311.01395 [hep-th].

[100] N. Arkani-Hamed, D. Baumann, H. Lee, and G. L. Pimentel, “The Cosmological Bootstrap: Inflationary Correlators from Symmetries and Singularities,” JHEP 04 (2020) 105, arXiv:1811.00024 [hep-th].

[101] D. Baumann, D. Green, A. Joyce, E. Pajer, G. L. Pimentel, C. Sleight, and M. Taronna, “Snowmass White Paper: The Cosmological Bootstrap,” in 2022 Snowmass Summer Study. 3, 2022. arXiv:2203.08121 [hep-th].

[102] C. Sleight, “A Mellin Space Approach to Cosmological Correlators,” JHEP 01 (2020) 090, arXiv:1906.12302 [hep-th].

[103] C. Sleight and M. Taronna, “Bootstrapping Inflationary Correlators in Mellin Space,” JHEP 02 (2020) 098, arXiv:1907.01143 [hep-th].

[104] C. Sleight and M. Taronna, “From AdS to dS exchanges: Spectral representation, Mellin amplitudes, and crossing,” Phys. Rev. D 104 no. 8, (2021) L081902, arXiv:2007.09993 [hep-th].

[105] Z. Qin and Z.-Z. Xianyu, “Helical inflation correlators: partial Mellin-Barnes and bootstrap equations,” JHEP 04 (2023) 059, arXiv:2208.13790 [hep-th].

[106] Z.-Z. Xianyu and H. Zhang, “Bootstrapping one-loop inflation correlators with the spectral decomposition,” JHEP 04 (2023) 103, arXiv:2211.03810 [hep-th].

[107] Z. Qin and Z.-Z. Xianyu, “Closed-form formulae for inflation correlators,” JHEP 07 (2023) 001, arXiv:2301.07047 [hep-th].

[108] Z. Qin and Z.-Z. Xianyu, “Inflation correlators at the one-loop order: nonanalyticity, factorization, cutting rule, and OPE,” JHEP 09 (2023) 116, arXiv:2304.13295 [hep-th].

[109] Z. Qin and Z.-Z. Xianyu, “Nonanalyticity and On-Shell Factorization of Inflation Correlators at All Loop Orders,” arXiv:2308.14802 [hep-th].

[110] D. Stefanyszyn, X. Tong, and Y. Zhu, “Cosmological Correlators Through the Looking Glass: Reality, Parity, and Factorisation,” arXiv:2309.07769 [hep-th].

[111] Z.-Z. Xianyu and J. Zang, “Inflation Correlators with Multiple Massive Exchanges,” arXiv:2309.10849 [hep-th].

[112] S. Weinberg, “Quantum contributions to cosmological correlations,” Phys. Rev. D 72 (2005) 043514, arXiv:hep-th/0506236.

[113] D. J. Mulryne, D. Seery, and D. Wesley, “Moment transport equations for non-Gaussianity,” JCAP 01 (2010) 024, arXiv:0909.2256 [astro-ph.CO].

[114] D. J. Mulryne, D. Seery, and D. Wesley, “Moment transport equations for the primordial curvature perturbation,” JCAP 04 (2011) 030, arXiv:1008.3159 [astro-ph.CO].
[115] M. Dias and D. Seery, “Transport equations for the inflationary spectral index,” Phys. Rev. D 85 (2012) 043519, arXiv:1111.6544 [astro-ph.CO].

[116] D. Seery, D. J. Mulryne, J. Frazer, and R. H. Ribeiro, “Inflationary perturbation theory is geometrical optics in phase space,” JCAP 09 (2012) 010, arXiv:1203.2635 [astro-ph.CO].

[117] G. J. Anderson, D. J. Mulryne, and D. Seery, “Transport equations for the inflationary trispectrum,” JCAP 10 (2012) 019, arXiv:1205.0024 [astro-ph.CO].

[118] J. Elliston, D. Seery, and R. Tavakol, “The inflationary bispectrum with curved field-space,” JCAP 11 (2012) 060, arXiv:1208.6011 [astro-ph.CO].

[119] D. J. Mulryne, “Transporting non-Gaussianity from sub to super-horizon scales,” JCAP 09 (2013) 010, arXiv:1302.3842 [astro-ph.CO].

[120] M. Dias, J. Frazer, and D. Seery, “Computing observables in curved multifield models of inflation—A guide (with code) to the transport method,” JCAP 12 (2015) 030, arXiv:1502.03125 [astro-ph.CO].

[121] M. Dias, J. Frazer, D. J. Mulryne, and D. Seery, “Numerical evaluation of the bispectrum in multiple field inflation—the transport approach with code,” JCAP 12 (2016) 033, arXiv:1609.00379 [astro-ph.CO].

[122] D. J. Mulryne and J. W. Ronayne, “PyTransport: A Python package for the calculation of inflationary correlation functions,” J. Open Source Softw. 3 no. 23, (2018) 494, arXiv:1609.00381 [astro-ph.CO].

[123] D. Seery, “CppTransport: a platform to automate calculation of inflationary correlation functions,” arXiv:1609.00380 [astro-ph.CO].

[124] J. W. Ronayne and D. J. Mulryne, “Numerically evaluating the bispectrum in curved field-space— with PyTransport 2.0,” JCAP 01 (2018) 023, arXiv:1708.07130 [astro-ph.CO].

[125] X. Chen, R. Easther, and E. A. Lim, “Large Non-Gaussianities in Single Field Inflation,” JCAP 06 (2007) 023, arXiv:astro-ph/0611645.

[126] X. Chen, R. Easther, and E. A. Lim, “Generation and Characterization of Large Non-Gaussianities in Single Field Inflation,” JCAP 04 (2008) 010, arXiv:0801.3295 [astro-ph].

[127] F. Arroja, A. E. Romano, and M. Sasaki, “Large and strong scale dependent bispectrum in single field inflation from a sharp feature in the mass,” Phys. Rev. D 84 (2011) 123503, arXiv:1106.5384 [astro-ph.CO].

[128] D. K. Hazra, L. Sriramkumar, and J. Martin, “BINGO: A code for the efficient computation of the scalar bi-spectrum,” JCAP 05 (2013) 026, arXiv:1201.0926 [astro-ph.CO].

[129] H. Funakoshi and S. Renaux-Petel, “A Modal Approach to the Numerical Calculation of Primordial non-Gaussianities,” JCAP 02 (2013) 002, arXiv:1211.3086 [astro-ph.CO].

[130] P. Adshead, W. Hu, and V. Miranda, “Bispectrum in Single-Field Inflation Beyond Slow-Roll,” Phys. Rev. D 88 no. 2, (2013) 023507, arXiv:1303.7004 [astro-ph.CO].

[131] V. Assassi, D. Baumann, D. Green, and L. McAllister, “Planck-Suppressed Operators,” JCAP 01 (2014) 033, arXiv:1304.5226 [hep-th].

[132] V. Sreenath, D. K. Hazra, and L. Sriramkumar, “On the scalar consistency relation away from slow roll,” JCAP 02 (2015) 029, arXiv:1410.0252 [astro-ph.CO].
[133] P. Clarke and E. P. S. Shellard, “Probing inflation with precision bispectra,” JCAP 08 (2021) 002, arXiv:2012.08546 [astro-ph.CO].

[134] D. H. Tran, Y. Wang, J. Yang, and Y. Zhu, “Numerical prescriptions of early-time divergences of the in-in formalism,” JCAP 07 no. 07, (2022) 014, arXiv:2202.06350 [hep-th].

[135] D. Seery and J. E. Lidsey, “Primordial non-Gaussianities in single field inflation,” JCAP 06 (2005) 003, arXiv:astro-ph/0503692.

[136] X. Chen, M.-x. Huang, S. Kachru, and G. Shiu, “Observational signatures and non-Gaussianities of general single field inflation,” JCAP 01 (2007) 002, arXiv:hep-th/0605045.

[137] D. Langlois, S. Renaux-Petel, D. A. Steer, and T. Tanaka, “Primordial perturbations and non-Gaussianities in DBI and general multi-field inflation,” Phys. Rev. D 78 (2008) 063523, arXiv:0806.0336 [hep-th].

[138] P. Creminelli, M. A. Luty, A. Nicolis, and L. Senatore, “Starting the Universe: Stable Violation of the Null Energy Condition and Non-standard Cosmologies,” JHEP 12 (2006) 080, arXiv:hep-th/0606090.

[139] C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, “The Effective Field Theory of Inflation,” JHEP 03 (2008) 014, arXiv:0709.0293 [hep-th].

[140] L. Senatore and M. Zaldarriaga, “The Effective Field Theory of Multifield Inflation,” JHEP 04 (2012) 024, arXiv:1009.2093 [hep-th].

[141] R. Arnowitt, S. Deser, and C. W. Misner, “Dynamical structure and definition of energy in general relativity,” Phys. Rev. 116 (Dec, 1959) 1322-1330. https://link.aps.org/doi/10.1103/PhysRev.116.1322.

[142] S. Garcia-Saenz, L. Pinol, and S. Renaux-Petel, “Revisiting non-Gaussianity in multifield inflation with curved field space,” JHEP 01 (2020) 073, arXiv:1907.10403 [hep-th].

[143] L. Pinol, “Multifield inflation beyond $N_{field} = 2$: non-Gaussianities and single-field effective theory,” JCAP 04 (2021) 002, arXiv:2011.05930 [astro-ph.CO].

[144] N. Arkani-Hamed, P. Creminelli, S. Mukohyama, and M. Zaldarriaga, “Ghost inflation,” JCAP 04 (2004) 001, arXiv:hep-th/0312100.

[145] S. Butchers and D. Seery, “Numerical evaluation of inflationary 3-point functions on curved field space—with the transport method \& CppTransport,” JCAP 07 (2018) 031, arXiv:1803.10563 [astro-ph.CO].

[146] F. Piazza and F. Vernizzi, “Effective field theory of cosmological perturbations,” Classical and Quantum Gravity 30 no. 21, (Oct, 2013) 214007. https://doi.org/10.1088/2F0264-9381%2F30%2F21%2F214007.

[147] D. Baumann and D. Green, “Signatures of Supersymmetry from the Early Universe,” Phys. Rev. D 85 (2012) 103520, arXiv:1109.0292 [hep-th].

[148] S. R. Behbahani, A. Dymarsky, M. Mirbabayi, and L. Senatore, “(Small) Resonant non-Gaussianities: Signatures of a Discrete Shift Symmetry in the Effective Field Theory of Inflation,” JCAP 12 (2012) 036, arXiv:1111.3373 [hep-th].

[149] D.-G. Wang, G. L. Pimentel, and A. Achúcarro, “Bootstrapping multi-field inflation: non-Gaussianities from light scalars revisited,” JCAP 05 (2023) 043, arXiv:2212.14035 [astro-ph.CO].
[150] S. Cremonini, Z. Lalak, and K. Turzyński, “Strongly coupled perturbations in two-field inflationary models,” *Journal of Cosmology and Astroparticle Physics* **2011** no. 03, (Mar, 2011) 016–016.
http://dx.doi.org/10.1088/1475-7516/2011/03/016.

[151] R. Bravo, G. A. Palma, and S. Riquelme, “A Tip for Landscape Riders: Multi-Field Inflation Can Fulfill the Swampland Distance Conjecture,” *JCAP* **02** (2020) 004, arXiv:1906.05772 [hep-th].

[152] A. J. Tolley and M. Wyman, “The Gelaton Scenario: Equilateral non-Gaussianity from multi-field dynamics,” *Phys. Rev. D* **81** (2010) 043502, arXiv:0910.1853 [hep-th].

[153] A. Achucarro, V. Atal, S. Cespedes, J.-O. Gong, G. A. Palma, and S. P. Patil, “Heavy fields, reduced speeds of sound and decoupling during inflation,” *Phys. Rev. D* **86** (2012) 121301, arXiv:1205.0710 [hep-th].

[154] S. Pi and M. Sasaki, “Curvature Perturbation Spectrum in Two-field Inflation with a Turning Trajectory,” *JCAP* **10** (2012) 051, arXiv:1205.0161 [hep-th].

[155] T. Grall and S. Melville, “Inflation in motion: unitarity constraints in effective field theories with (spontaneously) broken Lorentz symmetry,” *JCAP* **09** (2020) 017, arXiv:2005.02366 [gr-qc].

[156] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. IX. Constraints on primordial non-Gaussianity,” *Astron. Astrophys.* **641** (2020) A9, arXiv:1905.05697 [astro-ph.CO].

[157] E. Castillo, B. Koch, and G. Palma, “On the integration of fields and quanta in time dependent backgrounds,” *JHEP* **05** (2014) 111, arXiv:1312.3338 [hep-th].

[158] R. Gwyn, G. A. Palma, M. Sakellariadou, and S. Sypsas, “Effective field theory of weakly coupled inflationary models,” *JCAP* **04** (2013) 004, arXiv:1210.3020 [hep-th].

[159] J.-O. Gong, S. Pi, and M. Sasaki, “Equilateral non-Gaussianity from heavy fields,” *JCAP* **11** (2013) 043, arXiv:1306.3691 [hep-th].

[160] M. Junaid and D. Pogosyan, “Geometrical measures of non-Gaussianity generated from single field Inflationary models,” *Phys. Rev. D* **92** no. 4, (2015) 043505, arXiv:1506.01212 [hep-th].

[161] N. Bartolo, D. Cannone, and S. Matarrese, “The Effective Field Theory of Inflation Models with Sharp Features,” *JCAP* **10** (2013) 038, arXiv:1307.3483 [astro-ph.CO].

[162] P. Adshead and W. Hu, “Bounds on nonadiabatic evolution in single-field inflation,” *Phys. Rev. D* **89** no. 8, (2014) 083531, arXiv:1402.1677 [astro-ph.CO].

[163] D. Cannone, N. Bartolo, and S. Matarrese, “Perturbative Unitarity of Inflationary Models with Features,” *Phys. Rev. D* **89** no. 12, (2014) 127301, arXiv:1402.2258 [astro-ph.CO].

[164] S. Cespedes, V. Atal, and G. A. Palma, “On the importance of heavy fields during inflation,” *JCAP* **05** (2012) 008, arXiv:1201.4848 [hep-th].

[165] A. Achucarro, J.-O. Gong, S. Hardeman, G. A. Palma, and S. P. Patil, “Effective theories of single field inflation when heavy fields matter,” *JHEP* **05** (2012) 066, arXiv:1201.6342 [hep-th].

[166] A. Achucarro, V. Atal, B. Hu, P. Ortiz, and J. Torrado, “Inflation with moderately sharp features in the speed of sound: Generalized slow roll and in-in formalism for power spectrum and bispectrum,” *Physical Review D* **90** no. 2, (Jul, 2014).
https://doi.org/10.1103/2Fphysrevd.90.023511.

[167] A. Achúcarro, J.-O. Gong, G. A. Palma, and S. P. Patil, “Correlating features in the primordial spectra,” *Phys. Rev. D* **87** no. 12, (2013) 121301, arXiv:1211.5619 [astro-ph.CO].
[168] J.-O. Gong, K. Schalm, and G. Shiu, “Correlating correlation functions of primordial perturbations,” Phys. Rev. D 89 no. 6, (2014) 063540, arXiv:1401.4402 [astro-ph.CO].

[169] A. Achucarro, V. Atal, B. Hu, P. Ortiz, and J. Torrado, “Inflation with moderately sharp features in the speed of sound: Generalized slow roll and in-in formalism for power spectrum and bispectrum,” Phys. Rev. D 90 no. 2, (2014) 023511, arXiv:1404.7522 [astro-ph.CO].

[170] G. A. Palma, “Untangling features in the primordial spectra,” JCAP 04 (2015) 035, arXiv:1412.5615 [hep-th].

[171] A. Achucarro, J.-O. Gong, S. Hardeman, G. A. Palma, and S. P. Patil, “Features of heavy physics in the CMB power spectrum,” JCAP 01 (2011) 030, arXiv:1010.3693 [hep-ph].

[172] X. Chen, “Primordial Features as Evidence for Inflation,” JCAP 01 (2012) 038, arXiv:1104.1323 [hep-th].

[173] X. Chen, “Fingerprints of Primordial Universe Paradigms as Features in Density Perturbations,” Phys. Lett. B 706 (2011) 111–115, arXiv:1106.1635 [astro-ph.CO].

[174] T. Noumi and M. Yamaguchi, “Primordial spectra from sudden turning trajectory,” JCAP 12 (2013) 038, arXiv:1307.7110 [hep-th].

[175] R. Saito, M. Nakashima, Y.-i. Takamizu, and J. Yokoyama, “Resonant Signatures of Heavy Scalar Fields in the Cosmic Microwave Background,” JCAP 11 (2012) 036, arXiv:1206.2164 [astro-ph.CO].

[176] X. Chen and M. H. Namjoo, “Standard Clock in Primordial Density Perturbations and Cosmic Microwave Background,” Phys. Lett. B 739 (2014) 285–292, arXiv:1404.1536 [astro-ph.CO].

[177] X. Chen, M. H. Namjoo, and Y. Wang, “Models of the Primordial Standard Clock,” JCAP 02 (2015) 027, arXiv:1411.2349 [astro-ph.CO].

[178] M. Braglia, X. Chen, and D. K. Hazra, “Primordial standard clock models and CMB residual anomalies,” Phys. Rev. D 105 no. 10, (2022) 103523, arXiv:2108.10110 [astro-ph.CO].

[179] M. Braglia, X. Chen, D. K. Hazra, and L. Pinol, “Back to the features: assessing the discriminating power of future CMB missions on inflationary models,” JCAP 03 (2023) 014, arXiv:2210.07028 [astro-ph.CO].

[180] J. Chluba, J. Hamann, and S. P. Patil, “Features and New Physical Scales in Primordial Observables: Theory and Observation,” Int. J. Mod. Phys. D 24 no. 10, (2015) 1530023, arXiv:1505.01834 [astro-ph.CO].

[181] A. Slosar et al., “Scratches from the Past: Inflationary Archaeology through Features in the Power Spectrum of Primordial Fluctuations,” Bull. Am. Astron. Soc. 51 no. 3, (2019) 98, arXiv:1903.09883 [astro-ph.CO].

[182] C. Gordon, D. Wands, B. A. Bassett, and R. Maartens, “Adiabatic and entropy perturbations from inflation,” Phys. Rev. D 63 (2000) 023506, arXiv:astro-ph/0009131.

[183] S. Groot Nibbelink and B. J. W. van Tent, “Scalar perturbations during multiple field slow-roll inflation,” Class. Quant. Grav. 19 (2002) 613–640, arXiv:hep-ph/0107272.

[184] S. Renaux-Petel and K. Turzyński, “Geometrical Destabilization of Inflation,” Phys. Rev. Lett. 117 no. 14, (2016) 141301, arXiv:1510.01281 [astro-ph.CO].

[185] T. Tanaka and Y. Urakawa, “Dominance of gauge artifact in the consistency relation for the primordial bispectrum,” JCAP 05 (2011) 014, arXiv:1103.1251 [astro-ph.CO].
[186] T. Baldauf, U. Seljak, L. Senatore, and M. Zaldarriaga, “Galaxy Bias and non-Linear Structure Formation in General Relativity,” JCAP 10 (2011) 031, arXiv:1106.5507 [astro-ph.CO].

[187] E. Pajer, F. Schmidt, and M. Zaldarriaga, “The Observed Squeezed Limit of Cosmological Three-Point Functions,” Phys. Rev. D 88 no. 8, (2013) 083502, arXiv:1305.0824 [astro-ph.CO].

[188] P. Creminelli and M. Zaldarriaga, “Single field consistency relation for the 3-point function,” JCAP 10 (2004) 006, arXiv:astro-ph/0407059.

[189] C. Cheung, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, “On the consistency relation of the 3-point function in single field inflation,” JCAP 02 (2008) 021, arXiv:0709.0295 [hep-th].

[190] P. Creminelli, G. D'Amico, M. Musso, and J. Norena, “The (not so) squeezed limit of the primordial 3-point function,” JCAP 11 (2011) 038, arXiv:1106.1462 [astro-ph.CO].

[191] P. Creminelli, J. Noreña, and M. Simonović, “Conformal consistency relations for single-field inflation,” JCAP 07 (2012) 052, arXiv:1203.4595 [hep-th].

[192] L. Senatore and M. Zaldarriaga, “A Note on the Consistency Condition of Primordial Fluctuations,” JCAP 08 (2012) 001, arXiv:1203.6884 [astro-ph.CO].

[193] V. Assassini, D. Baumann, and D. Green, “On Soft Limits of Inflationary Correlation Functions,” JCAP 11 (2012) 047, arXiv:1204.4207 [hep-th].

[194] K. Hinterbichler, L. Hui, and J. Khoury, “An Infinite Set of Ward Identities for Adiabatic Modes in Cosmology,” JCAP 01 (2014) 039, arXiv:1304.5527 [hep-th].

[195] G. Cabass, E. Pajer, and F. Schmidt, “How Gaussian can our Universe be?,” JCAP 01 (2017) 003, arXiv:1612.00033 [hep-th].

[196] T. Suyama, Y. Tada, and M. Yamaguchi, “Local observer effect on the cosmological soft theorem,” PTEP 2020 no. 11, (2020) 113E01, arXiv:2008.13364 [astro-ph.CO].

[197] Y. Wang, “Inflation, Cosmic Perturbations and Non-Gaussianities,” Commun. Theor. Phys. 62 (2014) 109–166, arXiv:1303.1523 [hep-th].

[198] P. Creminelli, J. Noreña, and M. Simonović, “Conformal consistency relations for single-field inflation,” Journal of Cosmology and Astroparticle Physics 2012 no. 07, (Jul, 2012) 052–052. https://doi.org/10.1088/2F1475-7516%2F2012%2F07%2F052.

[199] L. Berezhiani and J. Khoury, “Slavnov-Taylor Identities for Primordial Perturbations,” JCAP 02 (2014) 003, arXiv:1309.4461 [hep-th].

[200] X. Chen, R. Ebadi, and S. Kumar, “Classical Cosmological Collider Physics and Primordial Features,” arXiv:2205.01107 [hep-ph].

[201] L. Bordin, P. Creminelli, A. Khmelnitsky, and L. Senatore, “Light Particles with Spin in Inflation,” JCAP 10 (2018) 013, arXiv:1806.10587 [hep-th].

[202] A. Achúcarro, V. Atal, C. Germani, and G. A. Palma, “Cumulative effects in inflation with ultra-light entropy modes,” JCAP 02 (2017) 013, arXiv:1607.08609 [astro-ph.CO].

[203] X. Chen and Y. Wang, “Quasi-Single Field Inflation with Large Mass,” JCAP 09 (2012) 021, arXiv:1205.0160 [hep-th].