# The Cosmological Flow: A Systematic Approach to Primordial Correlators

Lucas Pinol, Sébastien Renaux-Petel, Denis Werth

► To cite this version:

Lucas Pinol, Sébastien Renaux-Petel, Denis Werth. The Cosmological Flow: A Systematic Approach to Primordial Correlators. Journal of Cosmology and Astroparticle Physics, 2025, 2025 (02), pp.019. 10.1088/1475-7516/2025/02/019. [hal-04356198](https://hal.science/hal-04356198v1)

HAL Id: hal-04356198

https://hal.science/hal-04356198v1

Submitted on 16 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# The Cosmological Flow:

## A Systematic Approach to Primordial Correlators

Lucas Pinol,¹,² Sébastien Renaux-Petel,² and Denis Werth²

¹ Laboratoire de Physique de l'École Normale Supérieure, ENS, CNRS, Université PSL, Sorbonne Université, Université Paris Cité, F-75005, Paris, France

² Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France

### Abstract

The time evolution of primordial fluctuations conceals a wealth of insights into the high-energy physics at play during the earliest moments of our Universe, which is ultimately encoded in late-time spatial correlation functions. However, the conventional procedure to compute them is technically challenging, and a complete dictionary mapping the landscape of inflationary theories and the corresponding observable signatures is not yet available. In this paper, we develop a framework to compute tree-level cosmological correlators based on following their time evolution from their origin as quantum zero-point fluctuations to the end of inflation. From first principles, the structure of the bulk time evolution imposes a set of universal differential equations in time satisfied by equal-time correlators. We automatise the process of systematically solving these equations. This allows us to accurately capture all physical effects and obtain exact results in theories formulated at the level of inflationary fluctuations that include any number of degrees of freedom with arbitrary dispersion relations and masses, coupled through any time-dependent interactions. We then illustrate the power of this formalism by exploring the phenomenology of cosmological correlators emerging from the interaction with a massive scalar field. After an extensive analysis of the quadratic theory and classifying perturbativity bounds, we study both the size and the shape dependence of non-Gaussianities in the entire parameter space, including the strong mixing regime. We present novel characteristics of cosmological collider signals in (would be) single-, double-, and triple-exchange three-point correlators. In the presence of primordial features, after subtracting gauge artefacts unavoidably generated by a breaking of scale-invariance, we show that soft limits of cosmological correlators offer a new possibility to probe the inflationary landscape. Finally, we provide templates to search for in future cosmological surveys.
# Contents

<table><tr><td>1 Introduction</td><td>4</td></tr><tr><td>2 Following the Time Evolution of Correlators</td><td>8</td></tr><tr><td>2.1 Equal-time Primordial Correlators</td><td>8</td></tr><tr><td>2.2 Two-point Correlators</td><td>11</td></tr><tr><td>2.3 Three-point Correlators</td><td>12</td></tr><tr><td>2.4 Higher-order Correlators</td><td>14</td></tr><tr><td>2.5 Resumming Quadratic Mixings</td><td>18</td></tr><tr><td>3 Cosmological Flow</td><td>19</td></tr><tr><td>3.1 Landscape of Theories</td><td>19</td></tr><tr><td>3.2 Flow Equations</td><td>21</td></tr><tr><td>3.3 Initial Conditions</td><td>26</td></tr><tr><td>Applications</td><td>32</td></tr><tr><td>4 Inflationary Fluctuations</td><td>32</td></tr><tr><td>4.1 Goldstone Description</td><td>32</td></tr><tr><td>4.2 Quadratic Theory</td><td>36</td></tr><tr><td>4.2.1 Effective mass and resonance at strong mixing</td><td>36</td></tr><tr><td>4.2.2 Phase diagram</td><td>39</td></tr><tr><td>5 Strong Coupling Scales and Naturalness</td><td>43</td></tr><tr><td>5.1 Weak Mixing</td><td>44</td></tr><tr><td>5.2 Strong Mixing</td><td>46</td></tr><tr><td>5.3 Radiative Corrections</td><td>48</td></tr><tr><td>6 Size of non-Gaussianities beyond Weak Mixing</td><td>49</td></tr><tr><td>6.1 Non-local EFT</td><td>50</td></tr><tr><td>6.2 Phenomenology</td><td>50</td></tr><tr><td>7 Cosmological Collider Phenomenology</td><td>56</td></tr><tr><td>7.1 Strongly Mixed Colliders</td><td>57</td></tr><tr><td>7.2 Sound-Speed Colliders</td><td>59</td></tr><tr><td>7.3 Cosmological Collider Flow</td><td>62</td></tr><tr><td>8 Primordial Features</td><td>63</td></tr><tr><td>8.1 Effective Description of Features</td><td>63</td></tr><tr><td>8.2 Correlated Sharp Features</td><td>65</td></tr><tr><td>8.3 Probing Inflationary Landscapes with Cosmological Colliders</td><td>67</td></tr></table>
<table><tr><td>9 Conclusions and Outlook</td><td>74</td></tr><tr><td>A Details on the Cosmological Flow</td><td>78</td></tr><tr><td>A.1 Extended Fourier Summation</td><td>78</td></tr><tr><td>A.2 Two-point Correlator Initial Conditions</td><td>79</td></tr><tr><td>A.3 Three-point Correlator Initial Conditions</td><td>80</td></tr><tr><td>B Light Limit of the Dynamics</td><td>86</td></tr><tr><td>C Perturbativity Bounds</td><td>90</td></tr><tr><td>References</td><td>91</td></tr></table>
# 1 Introduction

Time plays a crucial role in cosmology. It orchestrates the gravitational collapse of cosmological structures that has shaped our Universe from the commencement to its present-day stage. While we only have access to late-time observables, such as statistical properties of cosmological structures, our challenge is to understand the cosmic history and explain the correlations we see. In this sense, cosmology is about tracing the evolution of the Universe back in time or, instead going forward in time, tracking down the “cosmological flow” of the observed correlations.

According to the current paradigm, the early Universe underwent an inflationary phase [1–8]. Primordial quantum fluctuations were stretched to cosmological scales and were ultimately imprinted as spatial correlations that reside on the future boundary of inflationary spacetime, marking the beginning of the subsequent hot Big Bang evolution [9–16]. These equal-time correlators of primordial density fluctuations and tensor modes encode valuable information about inflation, and—for connected $n$-point correlators with $n \ge 3$—are conventionally called non-Gaussianities [17–21]. Recent years have shown increasing interests in computing these correlators as they ultimately hold the key to high-energy physics.

On observational grounds, it is expected that substantial progress is to be made in the coming years. Compared to future CMB experiments that would lead to interesting yet modest improvements on bounds for non-Gaussianities [22], the qualitative leap ahead of us lies in the constraining power of three-dimensional large-scale structure surveys [23–26]. Indeed, a host of galaxy surveys ambition to improve current constraints by one order of magnitude in the near future and reach an important threshold marking the border between a weakly or strongly coupled description of the inflationary background dynamics ($f_{\text{NL}} \sim 1$) [18, 27]. In the longer term, 21cm surveys will have the raw statistical power to improve limits by several orders of magnitude [28–31], reaching the gravitational floor ($f_{\text{NL}} \sim 10^{-2}$) [17].

The phenomenological interest in primordial correlators lies in the fact that the energy scale of inflation—the Hubble scale—can be as high as $10^{14}$ GeV [32, 33]. This energy makes the study of inflationary correlators a unique opportunity to learn about potential new physics at the highest reachable energies, much higher than those achievable in ground-based experiments. This one-of-a-kind window to fundamental physics has opened a cosmological collider physics program [34–65]. Essentially, even very heavy particles have been spontaneously produced during inflation, whose subsequent decays lead to imprints in inflationary correlators. These specific signatures—whose prospects for detection are bright [66–72] although the experimental challenge to observe them is enormous—could inform us about the inflationary field content, including the mass spectra of particles, their spins and sound speeds.

On the theory side, the study of inflationary correlators has triggered recent developments both in deriving analytical closed-form solutions for correlators, and revealing formal aspects such as their analytic structure. These advancements mirror our understanding of flat-space scattering amplitudes or that of boundary correlators in anti-de Sitter (AdS) space, see e.g. [73–99]. The motivations for such progress lie in the fact that, on cosmological scales probed in the CMB and LSS, inflation is reliably described by weakly-coupled fluctuations on top of a fixed quasi-de Sitter (dS) space. New techniques, such as solving boundary equations in external momenta [100] (see
also the review [101] and the references therein) and AdS-inspired techniques [102–104] have led
to analytical results for the tree-level dS-invariant four-point exchange correlator mediated by a
massive field with integer spin. These closed-form expressions have also been reproduced using
partial Mellin-Barnes representations [59, 105] and extended to boost-breaking interactions and
reduced sound speeds [60, 61]. Such methods opened up new analytical results for more complex
diagrams, e.g. at loop-level or including a chemical potential [106–111].

A priori, we do have a well-defined method for evaluating equal-time correlators *without approximation*: the acclaimed in-in formalism [112]. From first principles, calculations can be carried to arbitrary orders in perturbation theory, including loop contributions. Having a complete formalism at hand, the question is then the following: how to efficiently and systematically compute inflationary correlators given *any* inflationary theory?

The standard procedure comes with its own set of challenges. The background time depen-
dence leads to complicated (or even unknown) mode functions that are rapidly oscillating in the
infinite past, and perturbative calculations require us to compute difficult time integrals over the
entire inflationary evolution. As a consequence, primarily due to technical considerations, the
vast landscape of inflationary theories currently remains beyond the scope of existing predictions,
potentially leading to a biased interpretation of forthcoming data. A more pragmatic attitude,
evading the intricacy of Feynman-like diagram computations, is to directly focus on correlators
as they are the observables today. Assuming that primordial correlators emerge from a universal
initial vacuum, they only depend on the specific theory active in the bulk of spacetime, or equiv-
alently, only depend on how they “flow” towards the end of inflation. Our aim is to *automate*
this picture.

In this paper, which follows the shorter letter [62], we present a systematic approach to
compute primordial correlators, tracing their cosmological flow, see Figure 1. This method is
based on computing inflationary equal-time correlators by solving differential equations in *time*
governing their time evolution throughout the entirety of spacetime during inflation, from their
origin as quantum fluctuations in the deep past to the end of inflation. It takes into account all
physical effects at tree-level without approximation, in any theory conveniently formulated at the

<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic diagram of the cosmological flow. Inflationary correlators are tracked from initial quantum fluctuations in the asymptotic past to the end of inflation where they become spatial correlations on the reheating surface. The systematic resolution of their bulk time evolution is achieved by solving differential equations in time satisfied by the correlators.</figcaption></figure>
level of inflationary fluctuations, which can include an arbitrary number of degrees of freedom with any dispersion relations and masses, coupled through any time-dependent interactions. Systematically solving the cosmological flow of inflationary correlators does not require far-past UV regulators as correlators are evolved on the real-time axis, eliminates the need for mode functions, and bypasses the computation of intricate bulk time integrals.

Being able to have direct access to the bulk time evolution of correlators is particularly useful in conjunction with analytical works and to extend the present reach of inflationary correlator phenomenology. It enables one to easily shed light on characteristic time scales at play and give new insights on identifying various physical regimes in complex scenarios that cannot be reached analytically. As an application, we show that the frequency of the cosmological collider signal may not be solely dictated by the mass of an additional field. Specifically, at strong mixing or in the presence of features—emerging from time-dependent couplings—, we derive new cosmological collider signals that break the direct link relating frequency and mass. For the numerous possibilities that our approach provides, we have made our numerical code CosmoFlow<fnref n="1" /> available for the community. With this tool as a first building-block, we pave the way for a far-reaching program of exploring the rich physics of inflation. Although we focus on tree-level three-point correlators of scalars in this work, the cosmological flow allows for straight extensions, for example to compute higher-order correlators or to include spinning fields.

Finally, similar to what automated tools brought to the particle physics community, we believe that the cosmological flow will offer new possibilities for the study of inflation. This shift in perspective enables us to move beyond the technical achievement of computing inflationary correlators within a limited subset of the landscape of theories and instead embark on a journey to explore the physics of inflation in full generality.

**Outline.** In Section 2, we first review how primordial correlators are computed from a bulk perspective using the standard in-in formalism. Then, we derive already-known differential equations in time satisfied by the two- and three-point correlators. To do so, we introduce convenient diagrammatic representations of these equations that enable us to derive these flow equations for any *n*-point correlators. In Section 3, we concretely apply this approach to compute the two- and three-point correlators in generic theories of inflationary fluctuations. Specifically, we show how the theory dependence is readily encoded in the flow equations and derive initial conditions for the correlators. The rest of the paper is dedicated to applications. In Section 4, we construct a general theory of inflationary fluctuations where the Goldstone boson of broken time translations is coupled to a massive scalar field and provide a comprehensive study of the quadratic theory in the full parameter space. In Section 5, we derive bounds on couplings for this theory to remain weakly coupled and discuss naturalness. In Sections 6 and 7, for the first time, we provide exact results for the size and shape of non-Gaussianities and a complete understanding of the physics at play in the entire parameter space, including the strong mixing regime. In Section 8, we showcase that our method can systematically account for time-dependent couplings. We study both features that cannot be described by an effective single-field theory, and the corresponding cosmological collider signals. We show, with minimal assumptions, that the link between the

<footnote n="1">¹github.com/deniswerth/CosmoFlow</footnote>
frequency of the cosmological collider signal and the mass of an additional particle can be explicitly broken, leading to striking new observable signatures. Our conclusions are presented in Section 9. For the enthusiastic readers, appendices contain additional details. The two parts of the paper can be read independently. In particular, readers who want to skip the details of the cosmological flow derivation and only want to see it in action can skip Sections 2 and 3 and go straight to the applications.

**Relation to previous works.** The possibility of using differential equations in time to follow primordial correlators on super-horizon scales has been first recognised in [113–118]. This method has been extended to sub-horizon scales in [119] where its equivalence with the in-in formalism has been demonstrated. The modern transport approach has been implemented and publicly released as freely-available software packages in [120–124]. As such, the cosmological flow is built upon the pre-established transport approach. Notably, out of courtesy to the authors, most of the cosmological flow implementation and notations have been borrowed from [121]. Other numerical schemes, see e.g. [125–134] rely on directly computing the in-in integrals with numerical mode functions. All the existing approaches are either tied to non-linear sigma models, governing both the inflationary background and the fluctuations, or not automatic as the form of the bulk time integrals is strongly theory dependent. For instance, even the single-field effective field theory of inflation at leading order in derivatives is not encapsulated in existing automated tools, as some operators simply do not exist in non-linear sigma models. Conversely, as it is implemented directly at the level of fluctuations, the cosmological flow is applicable to any theories in a systematic manner, including of course the subset mentioned above, but offering a much richer and direct access to the physics of inflation.

**Notation and conventions.** We will use natural units in which the speed of light $c$ and Planck's constant $\hbar$ are set to unity $c = \hbar = 1$, with reduced Planck mass $M_{\text{pl}} = 1/\sqrt{8\pi G} = 2.4 \times 10^{18}$ GeV, where $G$ denotes Newton's constant. Our metric signature is $(-+++)$.

Spatial three-dimensional vectors are written in boldface $\mathbf{k}$. We will use Greek letters $(\mu, \nu, \dots)$ for spacetime indices, and Latin letters $(i, j, \dots)$ for spatial indices as usual. Overdots and primes will denote derivatives with respect to cosmic (physical) time $t$ and conformal time $\tau$ defined by $d\tau = dt/a$, respectively. Our Fourier convention is

$$f_{\mathbf{k}} = \int d^3x f(x) e^{i\mathbf{k} \cdot \mathbf{x}}, \quad f(x) = \int \frac{d^3k}{(2\pi)^3} f_{\mathbf{k}} e^{-i\mathbf{k} \cdot \mathbf{x}}. \qquad (1.1)$$

A prime on a correlator is defined to mean that we drop the momentum conserving delta function

$$\langle \mathcal{O}_{\mathbf{k}_1} \dots \mathcal{O}_{\mathbf{k}_4} \rangle = (2\pi)^3 \delta^{(3)}(\mathbf{k}_1 + \dots + \mathbf{k}_4) \langle \mathcal{O}_{\mathbf{k}_1} \dots \mathcal{O}_{\mathbf{k}_4} \rangle'. \qquad (1.2)$$

The extended Fourier summation—used in the formalism presented in this paper—will be introduced throughout the main text when needed. Details are given in Appendix A.1.
## 2 Following the Time Evolution of Correlators

We will begin our study with a general discussion of inflationary correlators in light of the cosmological flow. Our focus will be to derive differential equations in time satisfied by cosmological correlators. Specifically, starting from the in-in formalism in 2.1, we will derive the so-called flow equations for the two-point correlators in 2.2 and three-point correlators in 2.3. We then show how the flow equations for the n-point correlators can be found using a diagrammatic representation in 2.4.

### 2.1 Equal-time Primordial Correlators

We will derive flow equations for cosmological correlators starting from the in-in formalism, which we briefly review here, setting up relevant notations and definitions.

**Heisenberg-picture operators.** Our aim is to compute the vacuum expectation value of operators $O(t)$ made up of products of various fields at a given time $t$, $\langle\Omega|O(t)|\Omega\rangle$, where $|\Omega\rangle$ is the vacuum of the full interacting theory. We will consider a generic theory of $N$ degrees of freedom $\varphi^\alpha$ (with Greek letters running from 1 to $N$) and remain voluntarily agnostic throughout this section on the precise form of the theory. It may contain arbitrary interactions and time-dependent coupling constants. The theory is described by a Hamiltonian $H(\varphi^\alpha, p^\beta)$ which is a functional of the phase-space coordinates i.e. the field-space coordinates $\varphi^\alpha$ and their conjugate momenta $p^\beta$.

**Interaction-picture operators.** Following the standard procedure, we go from the Heisenberg to the interaction picture, therefore splitting the Hamiltonian into a free and an interacting part (we will describe later the implications of various choices of splitting)

$$H(\varphi^{\alpha}, \mathbf{p}^{\alpha}) = H_0(\varphi^{\alpha}, \mathbf{p}^{\alpha}) + H_I(\varphi^{\alpha}, \mathbf{p}^{\alpha}). \qquad (2.1)$$

Let us also introduce the interaction-picture operators $\varphi^\alpha$ and $p^\beta$ defined in terms of the Heisenberg-picture operators $\varphi^\alpha$ and $p^\beta$ by

$$\varphi^{\alpha} \equiv \mathcal{U}^{\dagger}\varphi^{\alpha}\mathcal{U}, \quad p^{\alpha} \equiv \mathcal{U}^{\dagger}\mathbf{p}^{\alpha}\mathcal{U}, \qquad (2.2)$$

where $\mathcal{U}$ is a unitary operator that will be explicitly given below. Inserting the identity operator inside the expectation value of the Heisenberg-picture operator of interest $O(\varphi^\alpha, \mathbf{p}^\beta)$, one obtains

$$\langle\Omega|O(\varphi^{\alpha}, \mathbf{p}^{\beta})|\Omega\rangle = \langle\Omega|\mathcal{U}O(\varphi^{\alpha}, \mathbf{p}^{\beta})\mathcal{U}^{\dagger}|\Omega\rangle. \qquad (2.3)$$

For practical purposes, we gather all fields $\varphi^\alpha$ and conjugate momenta $p^\beta$ in a phase-space vector$^2$ $X^a \equiv (\varphi^\alpha, \mathbf{p}^\beta)$, and we do the same for the interaction-picture operators $X^a \equiv (\varphi^\alpha, p^\beta)$. Now,

<footnote>$^2$The Latin index $a$ runs over all phase-space coordinates from 1 to $2N$. We will also assume that such indices are organised so that a block of field labels is followed by a block of momentum labels, in the same order.</footnote>
by choosing the unitary operator $\mathcal{U}$ to satisfy the evolution equation

$$ \frac{d\mathcal{U}}{dt} = iH_I(\mathbf{X}^a)\mathcal{U} = i\mathcal{U}H_I(\mathbf{X}^a), \qquad (2.4) $$

one can verify that the interacting-picture fields evolve with the free Hamiltonian

$$ \frac{dX^a}{dt} = i[H_0(X^b), X^a]. \qquad (2.5) $$

The solution to Eq. (2.4) approaching identity in the asymptotic past is well known and is formally given by Dyson's formula

$$ \mathcal{U} = \bar{T} \exp \left( i \int_{-\infty^{+}}^{t} H_I(t') dt' \right), \qquad (2.6) $$

where we have denoted $H_I(t') = H_I[X^a(t')]$. The interacting Hamiltonians evaluated at $t'$ in the power series expansion of the exponential need to be written in time increasing order, with $\bar{T}$ being the anti-time ordering operator. For the Hermitian conjugate operator $\mathcal{U}^\dagger$, the time ordered operator $\bar{T}^\dagger = T$ rewrites the interacting Hamiltonians in decreasing time order. The lower limit $-\infty^+$ accounts for deforming the integration contour above the real axis at early times, with the fields appearing in the integral defined by analytical continuation $t \to t(1 + i\epsilon)$. Inside correlators, this $i\epsilon$ prescription enables us to evaluate the expectation value in the interacting vacuum by adiabatically switching off the interactions in the distant past. This procedure makes it possible to match the free theory vacuum $|0\rangle$ onto the full theory vacuum $|\Omega\rangle$ in the infinite past. Collecting Eqs. (2.3) and (2.6), we end up with the well-known in-in formula [112]

$$ \langle\Omega|\mathcal{O}(\mathbf{X}^a)|\Omega\rangle = \langle 0|\left[\bar{T} e^{i \int_{-\infty^{+}}^{t} H_I(t') dt'}\right] \mathcal{O}(X^a)\left[\bar{T} e^{-i \int_{-\infty^{-}}^{t} H_I(t') dt'}\right]|0\rangle. \qquad (2.7) $$

**Equations of motion.** In order to go further, we need to specify a certain form for the Hamiltonian. Since we are going to essentially work in Fourier space, we will adopt the extended Fourier summation convention that consists in using sans serif indices to denote a sum including integrals over Fourier modes, making the following expressions compact and manageable. Hence, in what follows, an index contraction 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 (2.8) $$

This summation convention has been used and detailed in [121]. For completeness, we give more details in Appendix A.1. Using this notation, we take the Hamiltonian to be

$$ H = \frac{1}{2!} H_{ab} \mathbf{X}^a \mathbf{X}^b + \frac{1}{3!} H_{abc} \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c + \frac{1}{4!} H_{abcd} \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \mathbf{X}^d + \dots. \qquad (2.9) $$

Without loss of generality, we take $H_{ab}, H_{abc}, \dots$ to be symmetric under the exchange of any indices.<fnref n="3" /> These tensors encode the various interactions in the theory, and conventionally en-

<footnote n="3">³In a quantum theory, operator ordering of the Hamiltonian a priori matters. Canonical quantisation consisting in promoting the phase-space variables of the classical Hamiltonian to quantum operators is ambiguous, as an</footnote>
compass functions of time and of the various momenta for derivative interactions. It is worth emphasising that this form of the Hamiltonian is completely general and captures all known effective field theories for scalar degrees of freedom at the level of inflationary fluctuations.<fnref n="4" /> The fully non-linear equations of motion then read

$$ 
\begin{aligned}
\frac{d\mathbf{X}^a}{dt} &= i[H, \mathbf{X}^a] \\
&= \epsilon^{ac} H_{cb} \mathbf{X}^b + \frac{1}{2!} \epsilon^{ad} H_{dbc} \mathbf{X}^b \mathbf{X}^c + \frac{1}{3!} \epsilon^{ae} H_{ebcd} \mathbf{X}^b \mathbf{X}^c \mathbf{X}^d + \dots \quad (2.10) \\
&= u^a{}_b \mathbf{X}^b + \frac{1}{2!} u^a{}_{bc} \mathbf{X}^b \mathbf{X}^c + \frac{1}{3!} u^a{}_{bcd} \mathbf{X}^b \mathbf{X}^c \mathbf{X}^d + \dots,
\end{aligned}
$$

where in the second line we have introduced the tensor $\epsilon^{ab}$ that is defined so that the commutation relation is written in the following compact form<fnref n="5" />

$$ [\mathbf{X}^a, \mathbf{X}^b] = i\epsilon^{ab}, \quad (2.11) $$

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}. \quad (2.12) $$

The third line of (2.10) should be seen as a definition for the tensors $u^a{}_b, u^a{}_{bc}, u^a{}_{bcd}, \dots$. We also have used the fact that the Hamiltonian has been taken to be fully symmetric<fnref n="6" /> and $\epsilon^{ab}$ is antisymmetric. Written in this form, it is clear that the equations of motion (2.10) encode both the full evolution of $\mathbf{X}^a$ and the commutation relations. In practice though, they are hard—if not impossible—to solve. One must therefore choose a simpler “free” Hamiltonian $H_0$ to evolve the interaction-picture operators, thus resorting to a perturbative description of the interactions encoded in $H_I = H - H_0$, as described by the in-in formula in Eq. (2.7).

**Tree-level cosmological correlators.** We choose the free Hamiltonian to be the full quadratic Hamiltonian, $H_0 = \frac{1}{2!} H_{ab} \mathbf{X}^a \mathbf{X}^b$, with the interacting part of the Hamiltonian $H_I$ being given by the cubic and higher orders. Doing so, the equations of motion verified by the interaction-picture fields and momenta are linear. Still, all quadratic mixings between the various fields and momenta are taken into account in a non-perturbative manner, as we will show more clearly in

<footnote>infinite number of quantum theories that reduce to the classical one in the limit $\hbar \to 0$ may be defined. In the context of quantum field theory though, ambiguities in the operator ordering eventually show up as different types of UV divergences to be renormalised. Since we will not perform renormalisation explicitly—which amounts to assume we directly work with the renormalised theory—we will overlook the operator ordering ambiguity in this work. We can therefore safely and without loss of generality symmetrise the classical Hamiltonian and then promote the fields and momenta to quantum operators. We wish to especially thank David Mulryne and David Seery for helpful discussions on this matter.</footnote>

<footnote n="4">⁴Higher-order time derivative interactions can lead to equations of motion with more than two time derivatives. As this situation signals the presence of an additional dynamical degree of freedom in the theory, one can introduce it in Eq. (2.9) so that the Hamiltonian always leads to equations of motion with no more than two time derivatives.</footnote>

<footnote n="5">⁵More details on this commutator in Fourier space is presented in Appendix A.1.</footnote>

<footnote n="6">⁶Note that this implies that the tensors $u^a{}_{bc}, u^a{}_{bcd}, \dots$ are symmetric under the exchange of any lower indices.</footnote>
2.5. With this splitting choice, the equations of motion for the interaction-picture operators read

$$ \frac{dX^a}{dt} = u^a{}_b X^b, \qquad (2.13) $$

which makes it explicit that $X^a$ evolves with the full quadratic Hamiltonian. By now, we have derived all the necessary fundamental elements to write all possible n-point correlation functions. Up to three-point correlators, expanding the exponentials in the in-in formula (2.7) and working at tree-level yields

$$ \begin{aligned} \langle \Omega | \mathbf{X}^a \mathbf{X}^b | \Omega \rangle &= \langle 0 | X^a X^b | 0 \rangle , \\ \langle \Omega | \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c | \Omega \rangle &= \langle 0 | \frac{i}{3!} \int_{-\infty}^{t} dt' H_{\text{def}} [X^d X^e X^f, X^a X^b X^c] | 0 \rangle . \end{aligned} \qquad (2.14) $$

From now on, we will not write $|\Omega\rangle$ nor $|0\rangle$ to make the notations less cluttered, assuming it is clear from the context whether the operators are in the Heisenberg ($X^a$) or interaction picture ($X^a$). We retain the commutator form of the correlators because it will be more convenient in the following development. Note that no mention of mode functions is needed in this language.

## 2.2 Two-point Correlators

We start by deriving the flow equations for the two-point functions $\langle X^a X^b \rangle$. This object is the equal-time tree-level two-point correlator, taking into account the full quadratic Hamiltonian. In what follows, it will be more convenient to adopt a diagrammatic representation of such objects. We use a white dot to denote an external field insertion at the time $t$

$$ \mathcal{O}^{a} = X^{a}(t), \qquad (2.15) $$

and a black dot to denote a single vertex insertion in the bulk at a time $t' < t$ which needs to be integrated over. The two-point correlator can be represented as

$$ \langle X^a X^b \rangle = \mathcal{O}^{a} \mathcal{O}^{b} . \qquad (2.16) $$

We will discuss various choices of interaction schemes and the resummation of quadratic mixings in 2.5, which will make the double-line notation clear. Since the states are time independent, one can formally differentiate (2.16) with respect to time by acting on the interaction-picture operators $X^a$ inside the correlator in (2.14). We then use the equation of motion (2.13) to find the following closed system

$$ \begin{aligned} \frac{d}{dt} \mathcal{O}^{a} \mathcal{O}^{b} &= \mathcal{O}^{a} \mathcal{O}^{b} + \mathcal{O}^{a} \mathcal{O}^{b} \\ &= u^{a}{}_{c} \mathcal{O}^{c} \mathcal{O}^{b} + u^{b}{}_{c} \mathcal{O}^{a} \mathcal{O}^{c}, \end{aligned} \qquad (2.17) $$

where we have colored in blue the parts of the correlator that have been differentiated with respect to time. Note that the introduced diagrammatic representation follows the usual Leibniz
product rule of differentiation. The sum over repeated indices makes it explicit that—as usual in quantum physics—we sum over all possible diagrams. The flow equations for the two-point correlators are then

$$ \frac{d}{dt} \langle \mathbf{X}^a \mathbf{X}^b \rangle = u^a_c \langle \mathbf{X}^c \mathbf{X}^b \rangle + u^b_c \langle \mathbf{X}^a \mathbf{X}^c \rangle. \qquad (2.18) $$

Of course, these equations are completely equivalent—in the sense that they encode the same physics—to the linear equations of motion supplemented by the quantisation condition. Consequently, they correctly capture all physical effects arising from quadratic operators in the theory. Note also that (2.18) couples all two-point correlators through the tensor $u^a_b$, including mixed propagators and correlators which contain conjugate momenta.

## 2.3 Three-point Correlators

Let us now derive the flow equations for the three-point correlators. From (2.14), we see that such correlators carry two time dependencies: (i) the external operators, and (ii) the upper limit of the integral over a bulk vertex, in the following denoted by a black dot. Diagrammatically, differentiating the three-point correlators with respect to time is then represented by

$$ \frac{d}{dt} = + + + + . \qquad (2.19) $$

Similar to the case of two-point correlators in the previous section, using the equation of motion (2.13), differentiating an external operator leads to

$$ = u^a_d . \qquad (2.20) $$

The other term—coming from deriving the cubic bulk vertex—requires a few manipulations. We give a formal derivation of this contribution in the insert below. In terms of diagrams, one needs to differentiate all combinations of a bulk vertex together with an external operator insertion, resulting in cutting the diagram in two-point correlators. The resulting two-point correlators
Consequently, one can use Wick's theorem to contract the various operators. Retaining only connected three-point functions, one obtains

$$ \frac{i}{3!} \langle H_{\text{def}} [X^d X^e X^f, X^a X^b X^c] \rangle = u^a_{\text{de}} \langle \mathbf{X}^b \mathbf{X}^d \rangle \langle \mathbf{X}^c \mathbf{X}^e \rangle + (2 \text{ perms}), \quad (2.25) $$

which, in the end, is given by the diagrammatic representation in (2.21).

**Linearity.** These equations are linear. This specific structure allows the flow of each kinematic configuration to be tracked independently. However, the source terms of these equations are non-linear, in the sense that the three-point correlators are sourced by a product of two two-point correlators. As a result, Eq. (2.22) does not form a closed system and one needs to solve for both the two- and three-point correlators. Intuitively, solving the two-point correlators is equivalent to solving the mode functions that are needed to compute higher-point correlators. The non-linear structure of the source is directly related to the fact that fields are quantum operators. Indeed, going from a three-point correlator—that is composed of six interaction-picture operators—to a product of two-point correlators—each composed of two operators—, we have used the commutation relation to reduce the number of operator insertions. We also remind the reader that working at tree-level means that non-linearities of the three-point correlators on the two-point correlators are neglected.

## 2.4 Higher-order Correlators

We now turn to higher-order correlators. We will see that the diagrammatic rules we introduced enable us to easily find the flow equations for any tree-level n-point correlator.

**Four-point correlators.** We first illustrate the machinery with four-point correlators. In what follows, we will only consider the connected parts of higher-point correlators, hence defining them by subtracting the disconnected contribution. In full generality, a four-point correlator can be contact-like—with quartic interactions encoded in $H_{abcd}$—or exchange-like—with two cubic interactions from $H_{abc}$. Both contributions are formally written

$$ \langle \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \mathbf{X}^d \rangle = \frac{i}{4!} \langle \int_{-\infty}^{t} dt' H_{\text{efgh}}(t') [X^e X^f X^g X^h(t'), X^a X^b X^c X^d] \rangle \\ - \frac{1}{(3!)^2} \langle \int_{-\infty}^{t} dt' dt'' H_{\text{efg}}(t'') H_{\text{hij}}(t') [X^e X^f X^g(t''), [X^h X^i X^j(t'), X^a X^b X^c X^d] \rangle. \quad (2.26) $$

We first treat the case of contact-like four-point correlators. As for the three-point correlators in 2.3, such correlators carry a time dependence in external operator insertions, and in the upper limit of the integral over the bulk vertex. Following the Leibniz product rule of differentiation,
which, in the end, is given by the diagrammatic representation in (2.29).

The case of exchange-like four-point correlators follows the same line. The diagrammatic representation of differentiating an exchange-like four point correlator is

$$ \frac{d}{dt} \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.5] \node (a) at (0,0) {$a$}; \node (b) at (0,1) {$b$}; \node (c) at (1,0) {$c$}; \node (d) at (1,1) {$d$}; \draw[fill] (a) circle (0.1) node[left] {$a$}; \draw[fill] (b) circle (0.1) node[left] {$b$}; \draw[fill] (c) circle (0.1) node[left] {$c$}; \draw[fill] (d) circle (0.1) node[left] {$d$}; \draw[fill] (a) circle (0.1) node[right] {$b$}; \draw[fill] (b) circle (0.1) node[right] {$d$}; \draw[fill] (c) circle (0.1) node[right] {$c$}; \draw[fill] (d) circle (0.1) node[right] {$d$}; \end{tikzpicture} = \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.5] \node (a) at (0,0) {$a$}; \node (b) at (0,1) {$b$}; \node (c) at (1,0) {$c$}; \node (d) at (1,1) {$d$}; \draw[fill] (a) circle (0.1) node[left] {$a$}; \draw[fill] (b) circle (0.1) node[left] {$b$}; \draw[fill] (c) circle (0.1) node[left] {$c$}; \draw[fill] (d) circle (0.1) node[left] {$d$}; \draw[fill] (a) circle (0.1) node[right] {$b$}; \draw[fill] (b) circle (0.1) node[right] {$d$}; \draw[fill] (c) circle (0.1) node[right] {$c$}; \draw[fill] (d) circle (0.1) node[right] {$d$}; \end{tikzpicture} + (3 \text{ perms}) + \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.5] \node (a) at (0,0) {$a$}; \node (b) at (0,1) {$b$}; \node (c) at (1,0) {$c$}; \node (d) at (1,1) {$d$}; \draw[fill] (a) circle (0.1) node[left] {$a$}; \draw[fill] (b) circle (0.1) node[left] {$b$}; \draw[fill] (c) circle (0.1) node[left] {$c$}; \draw[fill] (d) circle (0.1) node[left] {$d$}; \draw[fill] (a) circle (0.1) node[right] {$b$}; \draw[fill] (b) circle (0.1) node[right] {$d$}; \draw[fill] (c) circle (0.1) node[right] {$c$}; \draw[fill] (d) circle (0.1) node[right] {$d$}; \end{tikzpicture}. \quad (2.33) $$

Differentiating with respect to external operator insertions introduces the usual $u^a_b$ tensor. Differentiating the bulk vertices is given by contracting a three-point correlator with a two-point correlator using the tensor $u^a_{bc}$ in the following way

$$ \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.5] \node (a) at (0,0) {$a$}; \node (b) at (0,1) {$b$}; \node (c) at (1,0) {$c$}; \node (d) at (1,1) {$d$}; \draw[fill] (a) circle (0.1) node[left] {$a$}; \draw[fill] (b) circle (0.1) node[left] {$b$}; \draw[fill] (c) circle (0.1) node[left] {$c$}; \draw[fill] (d) circle (0.1) node[left] {$d$}; \draw[fill] (a) circle (0.1) node[right] {$b$}; \draw[fill] (b) circle (0.1) node[right] {$d$}; \draw[fill] (c) circle (0.1) node[right] {$c$}; \draw[fill] (d) circle (0.1) node[right] {$d$}; \end{tikzpicture} = \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.5] \node (a) at (0,0) {$a$}; \node (b) at (0,1) {$b$}; \node (c) at (1,0) {$c$}; \node (d) at (1,1) {$d$}; \draw[fill] (a) circle (0.1) node[left] {$a$}; \draw[fill] (b) circle (0.1) node[left] {$b$}; \draw[fill] (c) circle (0.1) node[left] {$c$}; \draw[fill] (d) circle (0.1) node[left] {$d$}; \draw[fill] (a) circle (0.1) node[right] {$b$}; \draw[fill] (b) circle (0.1) node[right] {$d$}; \draw[fill] (c) circle (0.1) node[right] {$c$}; \draw[fill] (d) circle (0.1) node[right] {$d$}; \end{tikzpicture} + (2 \text{ perms}) + (3 \text{ perms}), \quad (2.34) $$

where the 2 permutations are the number of ways to choose a pair of indices where one is summed over and the other is not summed over (exchanging the two dummy indices is already taken into account in the counting), and the 3 permutations are the usual cyclic permutations $(a \to b \to c \to d)$. The second line then contains 12 terms. This diagrammatic rule is proven in the following insert.

*Cubic exchange vertex.*—From (2.26), the term coming from deriving the cubic exchange vertex reads

$$ \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.5] \node (a) at (0,0) {$a$}; \node (b) at (0,1) {$b$}; \node (c) at (1,0) {$c$}; \node (d) at (1,1) {$d$}; \draw[fill] (a) circle (0.1) node[left] {$a$}; \draw[fill] (b) circle (0.1) node[left] {$b$}; \draw[fill] (c) circle (0.1) node[left] {$c$}; \draw[fill] (d) circle (0.1) node[left] {$d$}; \draw[fill] (a) circle (0.1) node[right] {$b$}; \draw[fill] (b) circle (0.1) node[right] {$d$}; \draw[fill] (c) circle (0.1) node[right] {$c$}; \draw[fill] (d) circle (0.1) node[right] {$d$}; \end{tikzpicture} = -\frac{1}{(3!)^2} \langle \int_{-\infty}^{t} dt' H_{\text{efg}}(t') H_{\text{hij}}[X^e X^f X^g(t'), [X^h X^i X^j, X^a X^b X^c X^d]] \rangle. \quad (2.35) $$

We first start by examining the inner commutator where the operators $X^a$ are evaluated at $t$. After rearranging it and introducing the tensor $u^a_{bc}$, we find

$$ -\frac{1}{3!} H_{\text{hij}} [X^h X^i X^j, X^a X^b X^c X^d] = \frac{i}{2!} u^a_{hi} X^h X^i X^b X^c X^d + (3 \text{ perms}). \quad (2.36) $$

Inside the expectation value, we now Wick contract a pair of external operators where one is the
operator whose index is summed over. Among the six such possibilities, half of them lead to the same contribution after relabelling the indices. This cancels the factor 2!. In the end, after identifying the three-point correlators, we end up with

$$ -\frac{1}{(3!)^2} \langle \int_{-\infty}^{t} dt' H_{\text{efg}}(t') H_{\text{hij}} [X^e X^f X^g(t'), [X^h X^i X^j, X^a X^b X^c X^d]] \rangle = \qquad (2.37) \\ u^a_{ef} \langle X^c X^e \rangle \langle X^b X^d X^f \rangle + (2 \text{ perms}) + (3 \text{ perms}), $$

where the 2 permutations account for choosing a pair among (e, b, c, d) with one index being e, and the 3 permutations are the usual cyclic permutations (a → b → c → d).

Collecting the various terms in Eqs. (2.27, 2.29, 2.33, 2.34), the flow equations for the connected four-point correlator are

$$ \frac{d}{dt} \langle X^a X^b X^c X^d \rangle = u^a_{\text{e}} \langle X^e X^b X^c X^d \rangle + u^a_{\text{efg}} \langle X^b X^e \rangle \langle X^c X^f \rangle \langle X^d X^g \rangle + (3 \text{ perms}) \qquad (2.38) \\ + u^a_{\text{ef}} \langle X^b X^e \rangle \langle X^c X^d X^f \rangle + (11 \text{ perms}). $$

Similar equations were found in [117] for super-horizon interactions. Here they are derived in a generic context for the first time. It is interesting to note that contact-like four-point diagrams are only sourced by products of two-point functions. Removing the last term on the right-hand side—that come from deriving the bulk vertices of an exchange-like diagram—, the above equations form a closed system along with the flow equations for the two-point correlators only. Contrariwise, exchange-like four-point diagrams are sourced by three-point correlators, reflecting their intrinsic complexity. One would then need to solve all flow equations for the two-, three-, and four-point correlators, which still form a closed triangular system.

**Diagrammatic rules for n-point correlators.** Having treated the simplest exchange-like diagram appearing in four-point correlators, we now state—without formal proof—general diagrammatic rules to derive the flow equations for any tree-level n-point correlators. These rules are the following:

* Draw a tree-level graph assigning an operator $X^a(t)$ to each of the n external field, denoted with a white dot. Bulk vertices, denoted with black dots, should be at least of cubic order because quadratic interactions are fully resummed. Hence, the external operators are connected to bulk vertices with doubled lines. Time dependence appears in external operators and in the upper limit of the integral over bulk vertices, considered as a whole.

* Following the Leibniz product rule of differentiation, taking the time derivative of the drawn diagram equals the sum of initial-like diagrams where time dependence has been isolated, in our notations coloured in blue.

* Differentiating an external operator $X^a$ introduces the tensor $u^a_b$ which contracts – b being
$\bar{M}_{\alpha\beta} = \bar{M}_{(\alpha\beta)}$. It contains the kinetic gradient terms and the mass matrix. The non-symmetric $\bar{I}_{\alpha\beta}$ tensor captures the remaining possible linear mixings among the fields. For the cubic part of the action, the tensors $\bar{A}_{\alpha\beta\gamma} = \bar{A}_{(\alpha\beta\gamma)}$ and $\bar{D}_{\alpha\beta\gamma} = \bar{D}_{(\alpha\beta\gamma)}$ should be symmetrised over all indices, and $\bar{B}_{\alpha\beta\gamma} = \bar{B}_{(\alpha\beta\gamma)}$ and $\bar{C}_{\alpha\beta\gamma} = \bar{C}_{(\alpha\beta\gamma)}$ over the first two indices, with corresponding exchange of momenta. No further assumptions than the ones stated above are made.<fnref n="10" />

**Conjugate momenta.** Given the action (3.1), the conjugate momenta to the fields are defined by the following functional derivative

$$ p_{\alpha}(t) = \frac{\delta S}{\delta \dot{\varphi}^{\alpha}(t)}, \quad \text{with} \quad \frac{\delta \varphi^{\alpha}(t)}{\delta \varphi^{\beta}(t')} = \delta_{\beta}^{\alpha} \delta(t-t'). \qquad (3.2) $$

At cubic order, it is enough to restrict ourselves to linear order to derive the conjugate momenta and use $\mathcal{H}^{(3)} = -\mathcal{L}^{(3)}$, where $\dot{\varphi}^{\alpha}$ should be expressed in terms of the linear momenta. Indeed, higher-order terms would give quartic terms in the Hamiltonian density, not relevant to compute three-point correlators. However, for the sake of completeness, we derive the conjugate momenta at second order (which would be necessary, e.g., for finding the quartic Hamiltonian),

$$ p_{\alpha} = a^3 \left( \bar{\Delta}_{\bar{\alpha}\beta} \dot{\varphi}^{\beta} + \bar{I}_{\bar{\alpha}\beta} \varphi^{\beta} + \frac{1}{2} \bar{B}_{\beta\gamma\bar{\alpha}} \varphi^{\beta} \varphi^{\gamma} + \bar{C}_{\bar{\alpha}\beta\gamma} \dot{\varphi}^{\beta} \varphi^{\gamma} + \frac{3}{2} \bar{D}_{\bar{\alpha}\beta\gamma} \dot{\varphi}^{\beta} \dot{\varphi}^{\gamma} \right), \qquad (3.3) $$

where bar indices indicate that the sign of the corresponding momentum $k$ has been reversed i.e. $T^{\bar{\alpha}} = T^{\alpha}(-k_{\alpha})$ where $T$ can be a tensor or a field. This follows from the fact that Fourier-space expressions can produce the delta function $\delta_{\alpha\beta} = (2\pi)^3 \delta_{\alpha\beta} \delta^{(3)}(k_{\alpha} + k_{\beta})$. Consistently, a contraction with $\delta_{\beta}^{\alpha}$ bars an index $\delta_{\beta}^{\alpha} T^{\beta} = T^{\bar{\alpha}}$. More details can be found in Appendix A.1.

**Hamiltonian.** The Hamiltonian—viewed as a functional of $(\varphi^{\alpha}, p_{\beta})$—is defined from the Lagrangian by a Legendre transform

$$ H = \int dt (p_{\alpha} \dot{\varphi}^{\alpha} - \mathcal{L}). \qquad (3.4) $$

From now on, we rescale the momenta $p_{\alpha} \rightarrow a^3 p_{\alpha}$ to avoid an extra overall factor of $a^3$ coming from $\mathcal{L}$ once the field time derivatives have been replaced by their conjugate momenta. Inverting Eq. (3.3) perturbatively gives

$$ \dot{\varphi}^{\alpha} = \bar{\nabla}^{\alpha} \delta \left[ \begin{aligned} & p_{\delta} - \bar{I}_{\delta\beta} \varphi^{\beta} - \left( \frac{1}{2} \bar{B}_{\beta\gamma\delta} - \bar{C}_{\delta\sigma\gamma} \bar{\nabla}^{\sigma\tau} \bar{I}_{\tau\beta} + \frac{3}{2} \bar{D}_{\delta\rho\mu} \bar{\nabla}^{\rho\tau} \bar{I}_{\tau\beta} \bar{\nabla}^{\mu\sigma} \bar{I}_{\sigma\gamma} \right) \varphi^{\beta} \varphi^{\gamma} \\ & - \left( \bar{C}_{\delta\sigma\gamma} \bar{\nabla}^{\sigma\beta} - \frac{3}{2} \bar{D}_{\delta\tau\rho} \bar{\nabla}^{\tau\beta} \bar{\nabla}^{\rho\sigma} \bar{I}_{\sigma\gamma} \right) p_{\beta} \varphi^{\beta} - \frac{3}{2} \bar{D}_{\delta\tau\sigma} \bar{\nabla}^{\tau\beta} \bar{D}^{\sigma\gamma} p_{\beta} p_{\gamma} + \dots \end{aligned} \right], \qquad (3.5) $$

<footnote n="10">¹⁰Note that the cubic interaction parameterised by the tensor $\bar{D}_{\alpha\beta\gamma}$ has not been considered in [121] as it does not appear in the action for the fluctuations in non-linear sigma models. It does appear in generalised $P(X^{IJ}, \phi^K)$ models though, where $X^{IJ} = -\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi^I\partial_{\nu}\phi^J$ are mixed kinetic terms [135–137], and is also present in the effective field theory of inflationary fluctuations [36,138–140]. The generic approach that we are adopting therefore encompasses all these classes of models directly at the level of fluctuations.</footnote>
where the ellipses denote higher-order terms, irrelevant in our case. We have denoted $\bar{\nabla}^{\alpha\beta}$ to be the inverse matrix satisfying $\bar{\nabla}^{\alpha\beta}\bar{\Delta}_{\beta\gamma} = \delta^{\alpha}_{\gamma}$. Rearranging the various terms, the Hamiltonian reads

$$H = \frac{1}{2} \int dt a^3 \left( \begin{aligned} & \Delta_{\alpha\beta} \boldsymbol{p}^{\alpha} \boldsymbol{p}^{\beta} - M_{\alpha\beta} \varphi^{\alpha} \varphi^{\beta} - 2I_{\alpha\beta} \varphi^{\alpha} \boldsymbol{p}^{\beta} \\ & - A_{\alpha\beta\gamma} \varphi^{\alpha} \varphi^{\beta} \varphi^{\gamma} - B_{\alpha\beta\gamma} \varphi^{\alpha} \varphi^{\beta} \boldsymbol{p}^{\gamma} - C_{\alpha\beta\gamma} \boldsymbol{p}^{\alpha} \boldsymbol{p}^{\beta} \varphi^{\gamma} - D_{\alpha\beta\gamma} \boldsymbol{p}^{\alpha} \boldsymbol{p}^{\beta} \boldsymbol{p}^{\gamma} \end{aligned} \right), \quad (3.6)$$

where the anew defined tensors are

$$\begin{aligned} \Delta_{\alpha\beta} &= \bar{\nabla}^{\bar{\alpha}\bar{\beta}}, & M_{\alpha\beta} &= \bar{M}_{\alpha\beta} - \bar{\nabla}^{\bar{\gamma}\bar{\delta}} \bar{I}_{\bar{\gamma}\alpha} \bar{I}_{\bar{\delta}\beta}, & I_{\alpha\beta} &= \bar{\nabla}^{\bar{\beta}\gamma} \bar{I}_{\bar{\alpha}\gamma}, \\ A_{\alpha\beta\gamma} &= \bar{A}_{\alpha\beta\gamma} - \bar{B}_{\alpha\beta\delta} \bar{\nabla}^{\delta\bar{\sigma}} \bar{I}_{\bar{\sigma}\gamma} + \bar{C}_{\delta\sigma\gamma} \bar{\nabla}^{\delta\bar{\tau}} \bar{I}_{\bar{\tau}\alpha} \bar{\nabla}^{\sigma\bar{\rho}} \bar{I}_{\bar{\rho}\beta} - \bar{D}_{\delta\sigma\rho} \bar{\nabla}^{\delta\mu} \bar{\nabla}^{\sigma\gamma} \bar{\nabla}^{\rho\kappa} \bar{I}_{\bar{\mu}\alpha} \bar{I}_{\bar{\nu}\beta} \bar{I}_{\bar{\kappa}\gamma}, \\ B_{\alpha\beta\gamma} &= \bar{B}_{\alpha\beta\delta} \bar{\nabla}^{\delta\bar{\gamma}} - 2\bar{C}_{\delta\sigma\alpha} \bar{\nabla}^{\delta\gamma} \bar{\nabla}^{\sigma\bar{\rho}} \bar{I}_{\bar{\rho}\beta} + 3\bar{D}_{\delta\sigma\rho} \bar{\nabla}^{\delta\bar{\gamma}} \bar{\nabla}^{\sigma\bar{\mu}} \bar{\nabla}^{\rho\bar{\nu}} \bar{I}_{\bar{\mu}\beta} \bar{I}_{\bar{\nu}\alpha}, \\ C_{\alpha\beta\gamma} &= \bar{C}_{\delta\sigma\gamma} \bar{\nabla}^{\delta\bar{\alpha}} \bar{\nabla}^{\sigma\bar{\beta}} - 3\bar{D}_{\delta\sigma\rho} \bar{\nabla}^{\alpha\delta} \bar{\nabla}^{\bar{\alpha}\bar{\sigma}} \bar{\nabla}^{\rho\bar{\tau}} \bar{I}_{\bar{\tau}\gamma}, \\ D_{\alpha\beta\gamma} &= \bar{D}_{\delta\sigma\rho} \bar{\nabla}^{\alpha\delta} \bar{\nabla}^{\beta\bar{\sigma}} \bar{\nabla}^{\gamma\bar{\rho}}. \end{aligned} \quad (3.7)$$

Notice that the internal target space of the fluctuations is flat. Consequently, the tensor indices are raised and lowered using the delta function. This form of the Hamiltonian—rather than the action (3.1)—can also be taken as a definition of the theory. Also, note that non-linear contributions to the momenta cancel in the cubic Hamiltonian. It is therefore enough to only consider linear momenta in the computation of three-point correlators.

## 3.2 Flow Equations

Given the explicit form of the Hamiltonian (3.6), it is an easy task to derive the $u^a_b = \epsilon^{ac}H_{cb}$ and $u^a_{bc} = \epsilon^{ad}H_{dbc}$ tensors entering in the flow equations (2.18) and (2.22). For the quadratic Hamiltonian, after carefully symmetrising the $H_{ab} = H_{(ab)}$ tensor, we obtain

$$H_{ab} = \begin{pmatrix} -M_{\alpha\beta} & -I_{\alpha\beta} \\ -I_{\beta\alpha} & \Delta_{\alpha\beta} \end{pmatrix}, \quad (3.8)$$

where the a, b indices are arranged in such a way that the field indices are followed by momentum indices in the same order. The desired $u^a_b$ tensor is then

$$u^a_b = \begin{pmatrix} \mathbf{0} & \mathbf{1} \\ -\mathbf{1} & \mathbf{0} \end{pmatrix} \begin{pmatrix} -M^{\bar{\alpha}} & -I^{\bar{\alpha}} \\ -I^{\bar{\alpha}} & \Delta^{\bar{\alpha}} \end{pmatrix} + \begin{pmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -3H\delta^{\bar{\alpha}} \end{pmatrix} = \begin{pmatrix} I_{\beta}^{\bar{\alpha}} & \Delta^{\bar{\alpha}} \\ M^{\bar{\alpha}} & I^{\bar{\alpha}} - 3H\delta^{\bar{\alpha}} \end{pmatrix}. \quad (3.9)$$
The extra factor $-3H\delta_{\beta}^{\bar{\alpha}}$ comes from the momentum rescaling.<fnref n="11" /> Similarly, the cubic Hamiltonian reads

$$H_{abc} = \left\{ \begin{pmatrix} -3A_{\alpha\beta\gamma} & -B_{\alpha\gamma\beta} \\ -B_{\gamma\beta\alpha} & -C_{\alpha\beta\gamma} \\ -B_{\alpha\beta\gamma} & -C_{\gamma\beta\alpha} \\ -C_{\alpha\gamma\beta} & -3D_{\alpha\beta\gamma} \end{pmatrix} \right\}, \qquad (3.11)$$

where the index c labels the first (resp. second) matrix within the braces $\{\dots\}$ if it is a field (resp. a momentum) index, and similarly a (resp. b) labels the row (resp. column) in each $2 \times 2$ block matrix. We recall that $H_{abc}$ is constructed to be symmetric under the exchange of any pair of indices hence the multiple presence of some tensors. The $u^a_{bc}$ tensor is obtained by explicitly performing the matrix multiplication

$$u^a_{bc} = \left\{ \begin{pmatrix} \mathbf{0} & \mathbf{1} \\ -\mathbf{1} & \mathbf{0} \end{pmatrix} \begin{pmatrix} -3A^{\bar{\alpha}}_{\beta\gamma} & -B^{\bar{\alpha}}_{\gamma\beta} \\ -B^{\bar{\alpha}}_{\gamma\beta} & -C^{\bar{\alpha}}_{\beta\gamma} \end{pmatrix} \right\} = \left\{ \begin{pmatrix} -B_{\gamma\beta}^{\bar{\alpha}} & -C^{\bar{\alpha}}_{\gamma\beta} \\ 3A^{\bar{\alpha}}_{\beta\gamma} & B^{\bar{\alpha}}_{\gamma\beta} \\ -C^{\bar{\alpha}}_{\gamma\beta} & 3D^{\bar{\alpha}}_{\beta\gamma} \\ B^{\bar{\alpha}}_{\beta\gamma} & C^{\bar{\alpha}}_{\gamma\beta} \end{pmatrix} \right\}. \qquad (3.12)$$

As done in [121], we further notice that the particular arrangement of barred indices—namely the upper $\bar{\alpha}$ index for each tensors—enables us to extract a delta function, hence defining $k$-dependent tensor coefficients $u^a_b$ and $u^a_{bc}$ that we label with normal indices $a, b, \dots$ in the following way

$$u^a_b = (2\pi)^3 \delta^{(3)}(\mathbf{k}_a - \mathbf{k}_b) u^a_b(\mathbf{k}_a, \mathbf{k}_b), \qquad (3.13)$$

$$u^a_{bc} = (2\pi)^3 \delta^{(3)}(\mathbf{k}_a - \mathbf{k}_b - \mathbf{k}_c) u^a_{bc}(\mathbf{k}_a, \mathbf{k}_b, \mathbf{k}_c). \qquad (3.14)$$

It is important to stress that in writing expressions with explicit momentum dependence, the index labelling the momentum is associated with the corresponding tensor index. Explicitly, $k_a$ is associated with $a$ and so forth. For example if one would want to exchange $b$ and $c$ in $u^a_{bc}$, the corresponding momenta should be exchanged too. Furthermore, due to the presence of the delta function and by isotropy, the $u^a_b$ tensor coefficients only depend on the magnitude of the momentum, $u^a_b(k)$ where $k = |\mathbf{k}_a| = |\mathbf{k}_b|$.

**Flow equations.** Turning to the two- and three-point correlation functions, statistical isotropy and homogeneity imply that one can write

$$\langle X^a X^b \rangle = (2\pi)^3 \delta^{(3)}(\mathbf{k}_a + \mathbf{k}_b) \Sigma^{ab}(k), \qquad (3.14) \\ \langle X^a X^b X^c \rangle = (2\pi)^3 \delta^{(3)}(\mathbf{k}_a + \mathbf{k}_b + \mathbf{k}_c) B^{abc}(k_a, k_b, k_c),$$

$$\frac{d\mathbf{p}_\alpha}{dt} = i[H, \mathbf{p}_\alpha] - 3H\mathbf{p}_\alpha. \qquad (3.10)$$

The additional non-canonical factor can be absorbed in the $u^a_b$ tensor as done in the main text.

<footnote n="11">¹¹Due to the momentum rescaling $p_\alpha \to a^3 p_\alpha$, the Heisenberg equation is modified</footnote>
where for the two-point correlator we have set $k = |\mathbf{k}_a| = |\mathbf{k}_b|$. We do not explicitly write the time dependency of the above correlators but it should be understood that they are evaluated at the same time $t$. The three-point correlators $B^{abc}(k_a, k_b, k_c)$ are real, as we will explain below, while $\Sigma^{ab}$ is complex and can be split into real and imaginary parts $\Sigma_R^{ab} = \Sigma_R^{ab} + i\Sigma_I^{ab}$. Because the flow equation (2.18) is linear and the $u^a_c$ have real elements, one can solve independently for $\Sigma_R^{ab}$ and $\Sigma_I^{ab}$, obeying

$$ \frac{d\Sigma_{R,I}^{ab}(k)}{dt} = u_c^a(k) \Sigma_{R,I}^{cb}(k) + u_c^b(k) \Sigma_{R,I}^{ac}(k). \qquad (3.15) $$

Actually, the commutation relation—valid at all times—fully fixes the imaginary part of the two-point correlators, independently of any theory. To see this, note that the real-space operators we consider are Hermitian. In Fourier space, this implies that

$$ \mathbf{X}^{a\dagger} = \mathbf{X}^{\bar{a}}, \qquad (3.16) $$

where $\dagger$ is the Hermitian conjugate operation, and we recall that $\bar{a}$ indicates that the sign of the corresponding momentum is reversed. From this, one can first derive symmetry properties satisfied by the two-point correlators. Indeed, $\langle \mathbf{X}^a \mathbf{X}^b \rangle^* = \langle (\mathbf{X}^a \mathbf{X}^b)^\dagger \rangle = \langle \mathbf{X}^{\bar{b}} \mathbf{X}^{\bar{a}} \rangle = \langle \mathbf{X}^b \mathbf{X}^a \rangle$, where we used (3.14) in the last step. This enforces the real part to be symmetric and the imaginary part to be anti-symmetric

$$ \Sigma_R^{ab} = \Sigma_R^{ba}, \quad \Sigma_I^{ab} = -\Sigma_I^{ba}. \qquad (3.17) $$

Second, introducing the commutator in $\langle \mathbf{X}^a \mathbf{X}^b \rangle$ and taking the imaginary part, one finds $\Sigma_I^{ab} = \frac{\epsilon^{ab}}{2a^3}$, where the tensor $\epsilon^{ab}$ is defined in (2.11). As a result, the imaginary part either decays in $1/a^3$ for correlators involving a field and a conjugate momentum of the same kind (recall that the physical momenta are rescaled by $a^3$) or is strictly zero (for all other correlators). Therefore, there is no need to solve the flow equation for $\Sigma_I^{ab}$, although this equation can all the same be solved as a numerical consistency check.

For three-point correlators, similar manipulations give

$$ \langle \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \rangle^* - \langle \mathbf{X}^{\bar{a}} \mathbf{X}^{\bar{b}} \mathbf{X}^{\bar{c}} \rangle = -i\epsilon^{ab} \langle \mathbf{X}^{\bar{c}} \rangle - i\epsilon^{ac} \langle \mathbf{X}^{\bar{b}} \rangle - i\epsilon^{bc} \langle \mathbf{X}^{\bar{a}} \rangle. \qquad (3.18) $$

The right-hand side vanishes for connected correlators, for which none of the wavevectors individually vanishes. Defining more precisely the bispectra as the connected three-point correlators, this leads to $B^{abc}(\mathbf{k}_a, \mathbf{k}_b, \mathbf{k}_c)^* - B^{abc}(-\mathbf{k}_a, -\mathbf{k}_b, -\mathbf{k}_c) = 0$. Then, statistical isotropy, implying that $B^{abc}$ in (3.14) only depends on the norms of the wavevectors, enforces that $B^{abc}$ is real. As for its time evolution, it is important to keep track of both $\Sigma_R^{ab}$ and $\Sigma_I^{ab}$ because the flow equations—being sourced by a non-linear term—requires the knowledge of the complex two-point correlators. Taking into account that the tensor coefficients $u^a_{bc}$ are real, we arrive at
$$ \frac{d B^{abc}(k_a, k_b, k_c)}{dt} = u^a_d(k_a) B^{dbc}(k_a, k_b, k_c) + u^b_d(k_b) B^{adc}(k_a, k_b, k_c) + u^c_d(k_c) B^{abd}(k_a, k_b, k_c) \\ + u^a_{de}(\mathbf{k}_a, \mathbf{k}_b, \mathbf{k}_c) \Sigma_R^{db}(k_b) \Sigma_R^{ec}(k_c) - u^a_{de}(\mathbf{k}_a, \mathbf{k}_b, \mathbf{k}_c) \Sigma_I^{db}(k_b) \Sigma_I^{ec}(k_c) \\ + u^b_{de}(\mathbf{k}_b, \mathbf{k}_a, \mathbf{k}_c) \Sigma_R^{ad}(k_a) \Sigma_R^{ec}(k_c) - u^b_{de}(\mathbf{k}_b, \mathbf{k}_a, \mathbf{k}_c) \Sigma_I^{ad}(k_a) \Sigma_I^{ec}(k_c) \tag{3.19} \\ + u^c_{de}(\mathbf{k}_c, \mathbf{k}_a, \mathbf{k}_b) \Sigma_R^{ad}(k_a) \Sigma_R^{be}(k_b) - u^c_{de}(\mathbf{k}_c, \mathbf{k}_a, \mathbf{k}_b) \Sigma_I^{ad}(k_a) \Sigma_I^{be}(k_b). $$

where we have deliberately written all the permutations to avoid ambiguity in the index ordering. This form of the flow equations for the three-point functions has been derived in [121].

We have seen that the imaginary parts of two-point correlators are fixed from first principles, and that the ones of connected three-point correlators vanish. In this respect, it is instructive to also consider four-point correlators. Systematically inserting commutators, one finds

$$ \mathrm{Im} \langle \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \mathbf{X}^d \rangle = \frac{1}{2i} [\langle \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \mathbf{X}^d \rangle - \langle \mathbf{X}^\bar{a} \mathbf{X}^\bar{b} \mathbf{X}^\bar{c} \mathbf{X}^\bar{d} \rangle] \tag{3.20} \\ + \frac{1}{2} (\epsilon^{cd} \langle \mathbf{X}^a \mathbf{X}^b \rangle + \epsilon^{ab} \langle \mathbf{X}^d \mathbf{X}^c \rangle + \epsilon^{bd} \langle \mathbf{X}^a \mathbf{X}^c \rangle + \epsilon^{ac} \langle \mathbf{X}^d \mathbf{X}^b \rangle + \epsilon^{bc} \langle \mathbf{X}^a \mathbf{X}^d \rangle + \epsilon^{ad} \langle \mathbf{X}^c \mathbf{X}^b \rangle). $$

The second line vanishes for connected four-point functions, and simply corresponds to the imaginary part of the correlator that one would deduce from Wick theorem in terms of two-point correlators. Note, however, that we did not use a perturbative scheme to derive this. Only the first line remains for connected correlators of interest. In that case, and contrary to three-point kinematic above, $(\mathbf{k}_a, \dots, \mathbf{k}_d)$ and $(-\mathbf{k}_a, \dots, -\mathbf{k}_d)$ need not be related by a rotation, so that the imaginary parts of connected four-point correlators do not vanish in general. This happens for theories in which parity is violated.

**Derivative interactions.** The entries of the cubic tensors $A_{\alpha\beta\gamma}$, $B_{\alpha\beta\gamma}$, $C_{\alpha\beta\gamma}$ and $D_{\alpha\beta\gamma}$ are in general time and momentum dependent. The momentum dependency comes from considering spatial derivative interactions, that give simple factors of momenta after Fourier transform. For instance, let us take the contact interaction $\mathcal{L}/a^3 \supset \frac{\lambda}{2}\varphi(\partial_i\varphi)^2/a^2$. In Fourier space, the corresponding cubic Hamiltonian gives the following tensor coefficient $A_{\varphi k_1 \varphi k_2 \varphi k_3} = \lambda \mathbf{k}_1 \cdot \mathbf{k}_2/a^2 + 2$ perms. In general, these tensor elements can be more complicated, and can for example scale as inverse powers of the momenta. Such momentum dependencies arise from non-local (in space) interactions, to which the cosmological flow has already been applied in [64]. In this way, the cosmological flow formalism makes it convenient to implement various types of interactions without much complications.

When computing the in-in integrals using the Schwinger-Keldysh formalism, the time derivatives appearing in a vertex should be directly applied to the attached propagators. In the cosmological flow formalism, time derivative interactions can be expressed as contact vertices between fields and momenta. For example, the interaction $\mathcal{L}/a^3 \supset -\frac{\lambda}{2}\varphi^2\dot{\varphi}$ gives the tensor coefficient $B_{\varphi k_1 \varphi k_2 p k_3} = -\lambda$, where $p$ is the (rescaled) conjugate momentum<fnref n="12" /> to the field $\varphi$. Therefore, it is an easy task to treat time derivative interactions with the formalism presented in this paper.

<footnote n="12">¹²We have assumed, for illustration purposes, that the quadratic theory is simple enough so that $p = \dot{\varphi}$.</footnote>
**Boundary terms.** When defining a theory, the action is sometimes supplemented by a boundary term. This is the case for example in the ADM formalism [141] where the Gibbons-Hawking-York boundary term needs to be added to the Einstein-Hilbert action. Another celebrated example is Maldacena’s calculation of the cubic Lagrangian in single-field inflation [17,135] and its more recent generalisation to non-linear sigma models in [142,143]. There, the cubic Lagrangian for the fluctuations is supplemented by a boundary term whose various terms come from rendering the size of interactions manifest by using the equation of motions, equivalently performing multiple integrations by parts. Although these boundary terms do not contribute to the equations of motion, they do contribute to the three-point correlators. We show here how to include these terms in correlators from first principles and in the most general case.

A generic boundary term appears in the form of a total derivative at the level of the action $S = S_{\text{bulk}} + S_{\partial}$. In Fourier space, we write the boundary term $S_{\partial}$ as

$$S_{\partial} = \int dt \mathcal{L}_{\partial}, \quad \text{with} \quad \mathcal{L}_{\partial} = \frac{d}{dt} \left( \frac{1}{3!} \Omega_{\text{abc}} \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \right), \qquad (3.21)$$

where in general the tensor $\Omega$ has time and momentum dependent elements. It is taken to be fully symmetric $\Omega_{\text{abc}} = \Omega_{(abc)}$. With the cubic Hamiltonian, one can formally write the in-in formula for the three-point correlators

$$\begin{aligned} \langle \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \rangle & \supset -i \langle 0 | \int_{-\infty}^{t} dt' \left[ \frac{d}{dt'} \left( \frac{1}{3!} \Omega_{\text{def}} X^d X^e X^f \right), X^a X^b X^c \right] | 0 \rangle \qquad (3.22) \\ & = -i \langle 0 | \frac{1}{3!} \Omega_{\text{def}} \left[ X^d X^e X^f, X^a X^b X^c \right] | 0 \rangle , \end{aligned}$$

where in the last line all the interaction-picture operators are evaluated at the time $t$. As we have done previously for the $u^a_b$ and $u^a_{bc}$ tensors, we define the tensor elements by extracting the momentum conserving delta functions which follows from that of the correlator $\Omega_{\text{abc}} = (2\pi)^3 \delta^{(3)}(\mathbf{k}_a + \mathbf{k}_b + \mathbf{k}_c) \Omega_{abc}(\mathbf{k}_a, \mathbf{k}_b, \mathbf{k}_c)$. As usual, the $i\epsilon$ prescription in the far past switches off the interactions. The equation shows that a general boundary term gives a local contribution in time for the three-point correlators. The fields being free in the interaction picture, we use Wick theorem to express the contribution into explicit two-point correlators. After using the symmetry properties of $\Omega_{\text{abc}}$ and retaining only connected contributions, one obtains

$$\langle \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \rangle \supset -i \Omega_{\text{def}} \left( \langle \mathbf{X}^d \mathbf{X}^a \rangle \langle \mathbf{X}^e \mathbf{X}^b \rangle \langle \mathbf{X}^f \mathbf{X}^c \rangle - \langle \mathbf{X}^a \mathbf{X}^d \rangle \langle \mathbf{X}^b \mathbf{X}^e \rangle \langle \mathbf{X}^c \mathbf{X}^f \rangle \right). \qquad (3.23)$$

We can further use the commutator to simplify the relation, obtaining

$$\begin{aligned} \langle \mathbf{X}^a \mathbf{X}^b \mathbf{X}^c \rangle \supset & -\Omega_{\text{def}} \left( \epsilon^{ad} \langle \mathbf{X}^e \mathbf{X}^b \rangle \langle \mathbf{X}^f \mathbf{X}^c \rangle + \epsilon^{be} \langle \mathbf{X}^a \mathbf{X}^d \rangle \langle \mathbf{X}^c \mathbf{X}^f \rangle \right. \\ & \left. \qquad + \epsilon^{cf} \langle \mathbf{X}^a \mathbf{X}^d \rangle \langle \mathbf{X}^e \mathbf{X}^b \rangle \right). \end{aligned} \qquad (3.24)$$

keeping in mind that three-point correlators are real, one finds that the boundary term only needs
to be added a posteriori:

$$B^{abc}(k_a, k_b, k_c) \rightarrow B_{\text{bulk}}^{abc}(k_a, k_b, k_c) + B_{\partial}^{abc}(k_a, k_b, k_c), \qquad (3.25)$$

with the contribution coming from the boundary term being

$$B_{\partial}^{abc}(k_a, k_b, k_c) = -\Omega_{\text{def}}(-\mathbf{k}_a, -\mathbf{k}_b, -\mathbf{k}_c) \left[ \begin{aligned} & \epsilon^{ad} \Sigma_R^{be}(k_b) \Sigma_R^{fc}(k_c) - \epsilon^{ad} \Sigma_I^{be}(k_b) \Sigma_I^{fc}(k_c) \\ & + \epsilon^{be} \Sigma_R^{ad}(k_a) \Sigma_R^{cf}(k_c) - \epsilon^{be} \Sigma_I^{ad}(k_a) \Sigma_I^{cf}(k_c) \\ & + \epsilon^{cf} \Sigma_R^{ad}(k_a) \Sigma_R^{be}(k_b) - \epsilon^{cf} \Sigma_I^{ad}(k_a) \Sigma_I^{be}(k_b) \end{aligned} \right]. \qquad (3.26)$$

Note that in general, only correlators of light degrees of freedom survive until the end of inflation. Also, terms proportional to spatial gradients are suppressed by the exponential of the time elapsed in e-fold since the corresponding mode exited the horizon, and are negligible for cosmological scales. In practice, a lot of the listed terms do not contribute significantly to the three-point correlators at the end of the time evolution. However for completeness we have derived the most general form for such boundary terms.

## 3.3 Initial Conditions

To render the system complete, one must provide initial conditions for $\Sigma_R^{ab}$, $\Sigma_I^{ab}$ and $B^{abc}$. Within the cosmological flow formalism, these initial conditions can be derived analytically provided one initialises the correlators sufficiently in the deep past. Indeed, by looking at the Hamiltonian in Eq. (3.6), we see that the quadratic theory is composed of a part governing the free evolution of the fields and one encoding non-trivial quadratic mixings. These quadratic interactions constitute one reason why analytical methods to compute late-time correlators fail in the most general cases. However, one can always choose an initial time deep enough inside the horizon for the modes of interest, so that the theory is dominated by the decoupled quadratic theory. In this regime, all the various couplings and masses can be neglected so that the theory approaches that of a set of uncoupled degrees of freedom, and initial conditions can be computed analytically. This is the very idea of asymptotically reaching the vacuum state.

### Two-point correlation functions

We begin with the second-order action in Fourier space,

$$S^{(2)} = \frac{1}{2} \int dt a^3 \left( \bar{\Delta}_{\alpha\beta} \dot{\varphi}^{\alpha} \dot{\varphi}^{\beta} + \bar{M}_{\alpha\beta} \varphi^{\alpha} \varphi^{\beta} + 2 \bar{I}_{\alpha\beta} \varphi^{\alpha} \dot{\varphi}^{\beta} \right), \qquad (3.27)$$

where $\bar{M}_{\alpha\beta}$ encompasses the kinetic gradient terms as well as other non-kinetic terms. We isolate the mass matrix $\bar{M}_{\alpha\beta}$ by writing

$$\bar{M}_{\alpha\beta} = (2\pi)^3 \delta^{(3)}(\mathbf{k}_{\alpha} + \mathbf{k}_{\beta}) (\bar{\mathcal{K}}_{\alpha\beta} + \bar{\mathcal{M}}_{\alpha\beta}), \qquad (3.28)$$
where all gradient terms are included in the matrix $\bar{\mathcal{K}}_{\alpha\beta}$. Without loss of generality, it can be written as

$$ \bar{\mathcal{K}}_{\alpha\beta} = -\frac{k^2}{a^2} c_{\alpha\beta}^{(2)} - \frac{k^4}{a^4} \frac{c_{\alpha\beta}^{(4)}}{\Lambda^2} - \frac{k^6}{a^6} \frac{c_{\alpha\beta}^{(6)}}{\Lambda^4} + \dots, \qquad (3.29) $$

with the ellipses denoting terms of higher order in $k$ that have been omitted, and $\Lambda$ is some energy scale to render the elements of the $c_{\alpha\beta}^{(2n)}$ matrix dimensionless. Specifically, $c_{\alpha\beta}^{(2)}$ encapsulates all the sound speeds (squared) for the various fields.

**Mass matrix and quadratic mixings.** In general, the mass matrix $\bar{\mathcal{M}}_{\alpha\beta}$ renders analytical treatment delicate because it models complex processes such as energy exchange between various fields at the quadratic level. It has been recently stressed that these interactions were analogous to neutrino and quark oscillations, as inflationary flavor and mass eigenstates generically mix with one other, leading to a new interesting phenomenology [55]. However, because the gradient terms in $\bar{\mathcal{K}}_{\alpha\beta}$ become increasingly important in the infinite past, they will eventually dominate the elements of $\bar{\mathcal{M}}_{\alpha\beta}$ provided that the (time-dependent) eigenvalues of the mass matrix are bounded from above. If one chooses an initial time such that contributions from the mass matrix become effectively irrelevant, each mode can be treated as massless. Therefore, we ignore the mass matrix for the remaining of this section. Similarly, one can always initialise the correlators deep enough in the past so that the gradient terms also dominate the quadratic couplings in the $\bar{I}_{\alpha\beta}$ matrix. We therefore ignore the quadratic coupling $\bar{I}_{\alpha\beta}$. In short, at sufficiently early times, the full quadratic theory will approach that of a set of massless and uncoupled degrees of freedom. Indeed, one can always diagonalise both kinetic terms with time and spatial derivatives in the same basis as they are both definite positive quadratic forms. Therefore, in the following, we will write in a non-covariant way the matrix elements $\bar{\Delta}_{\alpha\beta} = \delta^{\alpha\beta}\bar{\Delta}_{\alpha}$ and $\bar{\mathcal{K}}_{\alpha\beta} = \delta^{\alpha\beta}\bar{\mathcal{K}}_{\alpha}$.

**Linear dispersion relations.** For the remainder of this section, we assume that the kinetic terms for the fields only contain two spatial derivatives. We include the possibility of having different sound speeds for the fields by considering a generic diagonal matrix $c_{\alpha\beta}^{(2)} = \text{diag}(c_1^2, c_2^2, \dots)$. We highlight that the premise of our argument of analytically deriving initial conditions is not altered by this choice. However, other quadratic and diagonal theories such as theories with higher-order derivative terms (like ghost inflation [144]) require a different treatment. The quadratic theory we consider 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 (3.30) $$

**Initial two-point correlators.** Referring the interested reader to the Appendix A.2 for the derivation, the two-point correlator initial conditions are found to be

$$ \Sigma_{\mathcal{R}}^{ab} = \frac{1}{2a^2 c_{\alpha} k} \begin{pmatrix} \bar{\Delta}_{\alpha\beta} \left(1 + \frac{a^2 H^2}{c_{\alpha}^2 k^2}\right) & -H \delta^{\alpha\beta} \\ -H \delta^{\alpha\beta} & c_{\alpha}^2 k^2 \bar{\Delta}_{\alpha\beta}^{-1} \end{pmatrix}, \quad \Sigma_{\mathcal{I}}^{ab} = \frac{1}{2a^3} \begin{pmatrix} 0 & \delta^{\alpha\beta} \\ -\delta^{\alpha\beta} & 0 \end{pmatrix}, \qquad (3.31) $$

where all time-dependent quantities should be evaluated at the initial time. Note that in [121] the authors neglected the constant $\propto \frac{H^2}{2c_{\alpha}^3 k^3}$ mode in the $\langle \varphi^{\alpha} \varphi^{\beta} \rangle$ initial correlator. Even though
this term is negligible on subhorizon scales, we include it to improve numerical stability. We also recall that $\Sigma_{\mathcal{I}}^{ab}$ is completely fixed by the commutator, so that its initial condition in (3.31) is actually valid at all times. As an aside, notice that for small sound speeds $c_{\alpha} \ll 1$, the amplitude of $\langle \varphi^{\alpha} \varphi^{\beta} \rangle$ is enhanced whereas the amplitude of $\langle \mathbf{p}^{\alpha} \mathbf{p}^{\beta} \rangle$ is damped, intuitively reflecting the fact that the fields do not propagate in the limit $c_{\alpha} \to 0$.

## Three-point correlation functions

The derivation of initial three-point correlators is more involved than that of initial two-point correlators. Computing correlators in the deep past instead of in the far future requires some foresight. Indeed, examining the Hamiltonian (3.6), the tensors $A_{\alpha\beta\gamma}$, $B_{\alpha\beta\gamma}$, $C_{\alpha\beta\gamma}$ and $D_{\alpha\beta\gamma}$ encode cubic interactions that can be of arbitrary order in derivatives, each gradient (squared) term giving a term of order $(k/a)^2 \sim (k\tau)^2$ that grows exponentially as we approach the limit $|\tau| \to \infty$. Therefore, terms with higher-order derivative will be the dominant contributions of the various correlators. Specifically, let us extract the momentum-conserving delta function of the cubic tensors to define the associated kernels

$$A_{\alpha\beta\gamma} = (2\pi)^3 \delta^{(3)}(\mathbf{k}_{\alpha} + \mathbf{k}_{\beta} + \mathbf{k}_{\gamma}) A_{\alpha\beta\gamma}, \qquad (3.32)$$

and equivalent definitions for the tensors $B_{\alpha\beta\gamma}$, $C_{\alpha\beta\gamma}$ and $D_{\alpha\beta\gamma}$. We further decompose the above kernels by extracting the scale factor dependency that needs to be considered in the in-in integral, according to the number of spatial derivatives it contains, as follows<fnref n="13" />

$$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, \qquad (3.33)$$

and similarly for the other kernels. We have absorbed all momentum dependency in the newly defined kernels, and indicate with the superscript $2n$ ($n \ge 0$) the number of spatial derivatives that appears in the cubic interaction. For a given cubic interaction, the higher the number of spatial derivatives, the (exponentially) faster these terms grow on subhorizon scales. Consequently, it is essential to keep track of every order in powers of $\tau$ to avoid creating a large hierarchy between cubic interactions. By explicitly computing the in-in integrals, we generalise the results of [121, 145] to include all possible cubic interactions to all order in (spatial) derivatives. The full derivation is a bit lengthy, so we jump directly to the final results without loss of continuity, leaving all the details in Appendix A.3.

<footnote n="13">¹³Non-local interactions characterised by inverse gradient terms are completely negligible in this regime compared to usual derivative interactions so we do not consider this possibility.</footnote>
Three fields initial correlator:

$$
\begin{equation}
\begin{split}
\langle \varphi_{\mathbf{k}_1}^{\alpha} \varphi_{\mathbf{k}_2}^{\beta} \varphi_{\mathbf{k}_3}^{\gamma} \rangle_0 ={}& \frac{(2\pi)^3 \delta^{(3)}(\Sigma_i \mathbf{k}_i)}{4a^4 e_3 k_t} \left\{ \right. \\
& \left. \frac{A_{\alpha\beta\gamma}^{(0)}}{2} a^2 + \frac{A_{\alpha\beta\gamma}^{(2)}}{2} \left[ 1 + \frac{a^2 H^2}{e_3^2} \left( e_2 \left( e_2 + \frac{e_3}{k_t} \right) \right) - k_t e_3 \right] \right. \\
& + \frac{A_{\alpha\beta\gamma}^{(4)}}{2} \left[ \frac{1}{a^2} + \frac{H^2}{e_3^2} \left( e_2 \left( e_2 + \frac{3e_3}{k_t} \right) \right) - k_t e_3 \right] + \dots - \frac{B_{\alpha\beta\gamma}^{(0)}}{2} \frac{a^2 H}{k_1 k_2} [e_2 - k_3(k_1 + k_2)] \right. \\
& - \frac{B_{\alpha\beta\gamma}^{(2)}}{2} \frac{H}{k_1 k_2} \left[ e_2 - \frac{1}{k_t} (2e_3 + k_t(k_1 + k_2)k_3) \right] - \frac{B_{\alpha\beta\gamma}^{(4)}}{2} \frac{H}{a^2 k_1 k_2} \left[ e_2 - \frac{1}{k_t} (4e_3 + k_t(k_1 + k_2)k_3) \right] \\
& - \frac{C_{\alpha\beta\gamma}^{(0)}}{2} \left[ k_1 k_2 + a^2 H^2 \left( \frac{e_2}{k_3^2} \left( 1 + \frac{k_3}{k_t} \right) \right) - \frac{k_t}{k_3} \right] - \frac{C_{\alpha\beta\gamma}^{(2)}}{2} \left[ \frac{k_1 k_2}{a^2} + H^2 \left( \frac{e_2}{k_3^2} \left( 1 + \frac{3k_3}{k_t} \right) \right) - \frac{k_t}{k_3} \right] \\
& - \frac{C_{\alpha\beta\gamma}^{(4)}}{2} \left[ \frac{k_1 k_2}{a^4} + \frac{H^2}{a^2} \left( \frac{e_2}{k_3^2} \left( 1 + \frac{5k_3}{k_t} \right) \right) - \frac{k_t}{k_3} \right] + \dots + \frac{D_{\alpha\beta\gamma}^{(0)}}{2} H \left[ e_2 - \frac{2e_3}{k_t} \right] \\
& + \frac{D_{\alpha\beta\gamma}^{(2)}}{2} \frac{H}{a^2} \left[ e_2 - \frac{4e_3}{k_t} \right] + \frac{D_{\alpha\beta\gamma}^{(4)}}{2} \frac{H}{a^4} \left[ e_2 - \frac{6e_3}{k_t} \right] + \dots + \text{5 perms} \left. \right\}.
\end{split}
\tag{3.34}
\end{equation}
$$

Two fields and one momentum initial correlator:

$$
\begin{equation}
\begin{split}
\langle \boldsymbol{p}_{\mathbf{k}_1}^{\alpha} \boldsymbol{\varphi}_{\mathbf{k}_2}^{\beta} \boldsymbol{\varphi}_{\mathbf{k}_3}^{\gamma} \rangle_0 ={}& \frac{-(2\pi)^3 \delta^{(3)}(\Sigma_i \mathbf{k}_i) H}{4a^3 e_3^2 k_t} \left\{ \right. \\
& \left. -k_1 \left[ -\frac{A_{\alpha\beta\gamma}^{(0)}}{2} a^2 \left(e_2 - \frac{e_3}{k_t}\right) - \frac{A_{\alpha\beta\gamma}^{(2)}}{2} \left(e_2 + \frac{e_3}{k_t}\right) \right. \right. \\
& \left. - \frac{A_{\alpha\beta\gamma}^{(4)}}{2} \frac{1}{a^2} \left(e_2 + \frac{3e_3}{k_t}\right) + \dots + \frac{B_{\alpha\beta\gamma}^{(0)}}{2H} e_3 k_3 + \frac{B_{\alpha\beta\gamma}^{(2)}}{2H} \frac{e_3 k_3}{a^2} + \frac{B_{\alpha\beta\gamma}^{(4)}}{2H} \frac{e_3 k_3}{a^4} + \dots \right. \\
& \left. + \frac{C_{\alpha\beta\gamma}^{(0)}}{2} k_1^2 k_2^2 \left(1 + \frac{k_3}{k_t}\right) + \frac{C_{\alpha\beta\gamma}^{(2)}}{2} \frac{k_1^2 k_2^2}{a^2} \left(1 + \frac{3k_3}{k_t}\right) + \frac{C_{\alpha\beta\gamma}^{(4)}}{2} \frac{k_1^2 k_2^2}{a^4} \left(1 + \frac{5k_3}{k_t}\right) + \dots \right. \\
& \left. \left. - \frac{D_{\alpha\beta\gamma}^{(0)}}{2H} \frac{e_3^2}{a^2} - \frac{D_{\alpha\beta\gamma}^{(2)}}{2H} \frac{e_3^2}{a^4} - \frac{D_{\alpha\beta\gamma}^{(4)}}{2H} \frac{e_3^2}{a^6} + \dots + \text{5 perms} \right] - k_1^2(k_2 + k_3) \left[ \frac{A_{\alpha\beta\gamma}^{(0)}}{2} a^2 \right. \\
& \left. + \frac{A_{\alpha\beta\gamma}^{(2)}}{2} + \frac{A_{\alpha\beta\gamma}^{(4)}}{2} \frac{1}{a^2} + \dots + \frac{B_{\alpha\beta\gamma}^{(0)}}{2} a^2 H \frac{(k_1 + k_2)k_3^2}{e_3} + \frac{B_{\alpha\beta\gamma}^{(2)}}{2} \frac{Hk_3}{e_3 k_t} (2e_3 + k_t(k_1 + k_2)k_3) \right. \\
& \left. + \frac{B_{\alpha\beta\gamma}^{(4)}}{2} \frac{H}{e_3 k_t} \frac{k_3}{a^2} (4e_3 + k_t(k_1 + k_2)k_3) + \dots - \frac{C_{\alpha\beta\gamma}^{(0)}}{2} k_1 k_2 - \frac{C_{\alpha\beta\gamma}^{(2)}}{2} \frac{k_1 k_2}{a^2} - \frac{C_{\alpha\beta\gamma}^{(4)}}{2} \frac{k_1 k_2}{a^4} + \dots \right. \\
& \left. \left. - D_{\alpha\beta\gamma}^{(0)} \frac{He_3}{k_t} - 2D_{\alpha\beta\gamma}^{(0)} \frac{H}{a^2} \frac{e_3}{k_t} - 3D_{\alpha\beta\gamma}^{(0)} \frac{H}{a^4} \frac{e_3}{k_t} + \text{5 perms} \right] \left. \right\}.
\end{split}
\tag{3.35}
\end{equation}
$$
Two momenta and one field initial correlator:

$$
\begin{equation}
\begin{aligned}
\langle \boldsymbol{p}_{\boldsymbol{k}_1}^{\alpha} \boldsymbol{p}_{\boldsymbol{k}_2}^{\beta} \varphi_{\boldsymbol{k}_3}^{\gamma} \rangle_0 = & \frac{-(2\pi)^3 \delta^{(3)}(\Sigma_i \boldsymbol{k}_i)}{4a^4 e_3 k_t} \left\{ k_1 k_2 \right. \\
& \left. \left[ \begin{aligned}[t]
& \frac{A_{\alpha\beta\gamma}^{(0)}}{2} a^2 + \frac{A_{\alpha\beta\gamma}^{(2)}}{2} + \frac{A_{\alpha\beta\gamma}^{(4)}}{2} \frac{1}{a^2} + \dots \\
& + \frac{B_{\alpha\beta\gamma}^{(0)}}{2} a^2 H \frac{(k_1 + k_2)k_3^2}{e_3} + \frac{B_{\alpha\beta\gamma}^{(2)}}{2} H \frac{k_3}{e_3 k_t} (2e_3 + k_t(k_1 + k_2)k_3) + \dots \\
& + \frac{B_{\alpha\beta\gamma}^{(4)}}{2} \frac{H}{a^2} \frac{k_3}{e_3 k_t} (4e_3 + k_t(k_1 + k_2)k_3) + \dots \\
& - D_{\alpha\beta\gamma}^{(0)} \frac{He_3}{k_t} - 2D_{\alpha\beta\gamma}^{(2)} \frac{H}{a^2} \frac{e_3}{k_t} - 3D_{\alpha\beta\gamma}^{(4)} \frac{H}{a^4} \frac{e_3}{k_t} + \dots + \text{5 perms}
\end{aligned} \right] \\
& + k_1^2 k_2^2 \left[ \begin{aligned}[t]
& \frac{A_{\alpha\beta\gamma}^{(0)}}{2} \frac{a^4 H^2}{e_3} \left(e_2 - \frac{e_3}{k_t}\right) + \frac{A_{\alpha\beta\gamma}^{(2)}}{2} \frac{a^2 H^2}{e_3} \left(e_2 + \frac{e_3}{k_t}\right) + \frac{A_{\alpha\beta\gamma}^{(4)}}{2} \frac{H^2}{e_3} \left(e_2 + \frac{3e_3}{k_t}\right) + \dots \\
& - \frac{B_{\alpha\beta\gamma}^{(0)}}{2} a^2 H k_3 - \frac{B_{\alpha\beta\gamma}^{(2)}}{2} H k_3 - \frac{B_{\alpha\beta\gamma}^{(4)}}{2} \frac{H}{a^2 k_3} + \dots \\
& - \frac{C_{\alpha\beta\gamma}^{(0)}}{2} a^2 H^2 \frac{k_1^2 k_2^2}{e_3^2} \left(1 + \frac{k_3}{k_t}\right) + \dots \\
& - \frac{C_{\alpha\beta\gamma}^{(2)}}{2} H^2 \frac{k_1^2 k_2^2}{e_3^2} \left(1 + \frac{3k_3}{k_t}\right) - \frac{C_{\alpha\beta\gamma}^{(4)}}{2} \frac{H^2}{a^2} \frac{k_1^2 k_2^2}{e_3^2} \left(1 + \frac{5k_3}{k_t}\right) + \frac{D_{\alpha\beta\gamma}^{(0)}}{2} H + \frac{D_{\alpha\beta\gamma}^{(2)}}{2} \frac{H}{a^2} \\
& + \frac{D_{\alpha\beta\gamma}^{(4)}}{2} \frac{H^2}{a^4} + \dots + \text{5 perms}
\end{aligned} \right] \left. \right\} .
\end{aligned}
\tag{3.36}
\end{equation}
$$

Three momenta initial correlator:

$$
\begin{equation}
\begin{aligned}
\langle \boldsymbol{p}_{\boldsymbol{k}_1}^{\alpha} \boldsymbol{p}_{\boldsymbol{k}_2}^{\beta} \boldsymbol{p}_{\boldsymbol{k}_3}^{\gamma} \rangle_0 = & \frac{(2\pi)^3 \delta^{(3)}(\Sigma_i \boldsymbol{k}_i) H}{4a^3 e_3 k_t} \left\{ \right. \\
& \left. \frac{A_{\alpha\beta\gamma}^{(0)}}{2} a^2 \left(e_2 - \frac{e_3}{k_t}\right) + \frac{A_{\alpha\beta\gamma}^{(2)}}{2} \left(e_2 + \frac{e_3}{k_t}\right) \right. \\
& + \frac{A_{\alpha\beta\gamma}^{(4)}}{2} \frac{1}{a^2} \left(e_2 + \frac{3e_3}{k_t}\right) + \dots - \frac{B_{\alpha\beta\gamma}^{(0)}}{2H} e_3 k_3 - \frac{B_{\alpha\beta\gamma}^{(2)}}{2H} \frac{e_3 k_3}{a^2} - \frac{B_{\alpha\beta\gamma}^{(4)}}{2H} \frac{e_3 k_3}{a^4} + \dots \\
& - \frac{C_{\alpha\beta\gamma}^{(0)}}{2} k_1^2 k_2^2 \left(1 + \frac{k_3}{k_t}\right) - \frac{C_{\alpha\beta\gamma}^{(2)}}{2} \frac{k_1^2 k_2^2}{a^2} \left(1 + \frac{3k_3}{k_t}\right) - \frac{C_{\alpha\beta\gamma}^{(4)}}{2} \frac{k_1^2 k_2^2}{a^4} \left(1 + \frac{5k_3}{k_t}\right) + \dots \\
& \left. + \frac{D_{\alpha\beta\gamma}^{(0)}}{2H} \frac{e_3^2}{a^2} + \frac{D_{\alpha\beta\gamma}^{(2)}}{2H} \frac{e_3^2}{a^4} + \frac{D_{\alpha\beta\gamma}^{(4)}}{2H} \frac{e_3^2}{a^6} + \dots + \text{5 perms} \right\},
\end{aligned}
\tag{3.37}
\end{equation}
$$

where we have introduced the symmetric polynomials $e_1 \equiv k_t = k_1 + k_2 + k_3$, $e_2 \equiv k_1 k_2 + k_1 k_3 + k_2 k_3$, and $e_3 \equiv k_1 k_2 k_3$. We have coloured in blue the powers of the scale factor to render the size of the various terms in the subhorizon regime manifest. The momentum dependencies of the various kernels are not written explicitly but these need to be regarded as functions of $\mathbf{k}_1, \mathbf{k}_2$ and $\mathbf{k}_3$ associated with each index $\alpha, \beta$ and $\gamma$ respectively. When writing the permutations for the kernels $A_{\alpha\beta\gamma}^{(2n)}, B_{\alpha\beta\gamma}^{(2n)}, C_{\alpha\beta\gamma}^{(2n)}$ and $D_{\alpha\beta\gamma}^{(2n)}$, one should also exchange the associated momenta. Additional manipulations are needed when the fields are not canonically normalised. We leave this detail to Appendix A.3. In the end, note that these initial conditions are completely fixed
by the tensor elements, or equivalently by the Hamiltonian.

An alternative valid procedure consists in initialising the three-point correlators to zero and letting them build up as the cubic interactions are switched on adiabatically. This method essentially encapsulates the essence and *defines* the $i\epsilon$ prescription. Such prescription is simple and useful when dealing with other interactions and higher-order correlators. We have checked numerically for several examples that it is equivalent to using the initial conditions (3.34)–(3.37).
# Applications

## 4 Inflationary Fluctuations

With the cosmological flow method at hand, we begin our exploration with the study of primordial fluctuations in the line of the EFT of inflationary fluctuations [139]. We will present an exhaustive study of the quadratic theory of the Goldstone boson of broken time translations coupled to an additional massive field. We will also provide new insight on regimes that are not accessible by analytical means, which includes the strong mixing regime.

### 4.1 Goldstone Description

First, we will construct and study a multi-field theory at the level of the fluctuations built upon the EFT of (multi-field) inflation. We refer the reader to the original papers [138–140] or review [146] for more details. This framework provides a unified way of studying inflationary correlators because it is independent of the microscopic details of the theory that gives rise to the quasi-de Sitter background. The unavoidable degree of freedom describing inflationary fluctuations is the Goldstone boson associated with the spontaneous breaking of time-translation invariance.

**Goldstone boson sector.** In the unitary gauge, where the spatial part of the metric reads $g_{ij} = a^2(t)e^{2\zeta(t,x)}\delta_{ij}$ after disregarding tensor modes, the most general local action is

$$ S_{\pi} = \int d^4x \sqrt{-g} \left[ \frac{1}{2} M_{\text{pl}}^2 R + M_{\text{pl}}^2 \dot{H} g^{00} - M_{\text{pl}}^2 (3H^2 + \dot{H}) \right. \\ \left. + \frac{1}{2!} M_2^4(t) (\delta g^{00})^2 + \frac{1}{3!} M_3^4(t) (\delta g^{00})^3 + \dots \right. \\ \left. - \frac{1}{2} \bar{M}_1^3(t) \delta g^{00} \delta K_{\mu}^{\mu} - \frac{1}{2} \bar{M}_2^2(t) (\delta K_{\mu}^{\mu})^2 - \frac{1}{2} \bar{M}_3^2(t) \delta K^{\mu\nu} \delta K_{\mu\nu} + \dots \right] , \quad (4.1) $$

where $\delta g^{00} = g^{00} + 1$, $\delta K_{\mu\nu} = K_{\mu\nu} - a^2 H h_{\mu\nu}$ with $h_{\mu\nu}$ the induced spatial metric, and $M_{2,3}$ and $\bar{M}_{1,2,3}$ are time-dependent mass scales. The dots stand for higher-order terms in the fluctuations or with more spatial derivatives. This effective action encapsulates a large class of single-field models of inflation. In the following, we will be interested in the case where $\dot{H}$ and the mass scales $M_{2,3}$ and $\bar{M}_{1,2,3}$ are slowly time-varying,<fnref n="14" /> so that our treatment captures the behaviour of primordial fluctuations in the limit of approximate scale invariance.

To render the dynamics of this theory more explicit, we introduce the Goldstone boson $\pi$ associated with the spontaneous breaking of time-translation invariance, i.e. $\pi$ restores full diffeomorphism invariance by transforming as $\pi \to \pi - \xi(t, x)$ under a time reparameterisation $t \to t + \xi(t, x)$ such that $t + \pi$ is invariant. This is known as the Stückelberg trick. In this way, the field $\pi$ parametrises adiabatic perturbations corresponding to a local shift in time for the matter fields. A precise definition of the gauge in which $\pi$ appears is often not mentioned, as this

<footnote n="14">¹⁴We will explore the possibility of strongly time-varying mass scales in Section 8, leading to features.</footnote>
becomes irrelevant when considering the decoupling limit. In the following, we will also concentrate on theories and observables for which the mixing with gravity is negligible. For definiteness, however, we define $\pi$ in the spatially flat gauge, with $g_{ij} = a^2(t)\delta_{ij}$, so that the link between $\pi$ and $\zeta$ is unambiguous. The decoupling limit $M_{\text{pl}} \to \infty$ and $\dot{H} \to 0$ is taken while keeping the product $M_{\text{pl}}^2\dot{H}$ fixed. In this way, the inflationary correlators determined by the decoupled dynamics of $\pi$ are accurate up to slow-roll corrections $\epsilon = -\dot{H}/H^2$. After introducing the Goldstone boson and in the decoupling limit, $\delta g^{00}$ becomes $\delta g^{00} \to -2\dot{\pi} - \dot{\pi}^2 + (\partial_i\pi)^2/a^2$. Implementing this procedure in the action (4.1), neglecting terms that involve the extrinsic curvature<fnref n="15" /> and only keeping the leading deformation of the slow-roll action, one obtains the following single-field Lagrangian for the $\pi$-sector

$$ \mathcal{L}_{\pi}/a^3 = \frac{M_{\text{pl}}^2 |\dot{H}|}{c_s^2} \left[ \dot{\pi}^2 - c_s^2 \frac{(\partial_i \pi)^2}{a^2} + (1 - c_s^2) \left( \dot{\pi}^3 - \dot{\pi} \frac{(\partial_i \pi)^2}{a^2} \right) - \frac{4}{3} M_3^4 \frac{c_s^2}{M_{\text{pl}}^2 |\dot{H}|} \dot{\pi}^3 \right], \quad (4.2) $$

where $c_s^{-2} = 1 - 2M_2^4/M_{\text{pl}}^2|\dot{H}|$ is the intrinsic speed of sound for the propagation of $\pi$, generated by the lowest-order Lorentz-breaking operator proportional to $M_2$. This effective theory makes the non-linearly realised symmetry manifest, i.e. a small sound speed in the quadratic theory necessarily engenders a large $\dot{\pi}(\partial_i\pi)^2$ cubic self-interaction.

**Additional field sector.** We are interested in coupling the adiabatic mode $\pi$ to an additional massive scalar degree of freedom $\sigma$ which may have significant self-interactions. We consider the Lagrangian up to dimension-4 operators for the $\sigma$-sector

$$ \mathcal{L}_{\sigma}/a^3 = -\frac{1}{2}(\partial_{\mu}\sigma)^2 - \frac{1}{2}m^2\sigma^2 - \mu\sigma^3. \quad (4.3) $$

At cubic order, the $\sigma$ self-interactions are completely characterised by $\mu$.<fnref n="16" />

**Mixing sector.** Our interest is the coupling between the two sectors, both converting intrinsic non-linearities of $\sigma$ to the $\pi$ sector via a quadratic mixing, and generating non-linear interactions amongst them. In the unitary gauge, we consider the following mixing action

$$ S_{\pi-\sigma} = - \int d^4 x \sqrt{-g} \left[ \tilde{M}_1^3(t) \delta g^{00} \sigma + \tilde{M}_2^2(t) \delta g^{00} \sigma^2 + \tilde{M}_3^3(t) (\delta g^{00})^2 \sigma + \dots \right], \quad (4.4) $$

<footnote n="15">¹⁵Operators that involve the extrinsic curvature—being terms with higher derivatives—are suppressed by extra powers of the cutoff of the theory and hence are irrelevant at low energy. Even though these terms can become important in some cases—e.g. in ghost inflation [144]—we do not treat them in this paper. Nonetheless, they can be easily implemented in the cosmological flow formalism.</footnote>

<footnote n="16">¹⁶The dimension-4 operator $\dot{\sigma}\sigma^2$ is redundant with $\sigma^3$ upon integration by parts.</footnote>
where $\tilde{M}_{1,2,3}$ are some—in general time-dependent—mass scales.<fnref n="17" /> Taking the decoupling limit and introducing $\pi$, the theory up to cubic interactions reads

$$ \mathcal{L}_{\pi-\sigma}/a^3 = 2\tilde{M}_1^3 \dot{\pi}\sigma + \tilde{M}_1^3 (\partial_\mu \pi)^2 \sigma + 2\tilde{M}_2^2 \dot{\pi}\sigma^2 - 4\tilde{M}_3^3 \dot{\pi}^2 \sigma. \quad (4.5) $$

Similarly to the $\pi$-sector treated previously, we note that the Wilson coefficient of the $(\partial_i \pi)^2 \sigma$ operator is completely fixed by the $\dot{\pi}\sigma$ operator that appears in the quadratic Lagrangian, again a consequence of the fact that $\pi$ non-linearly realises time diffeomorphisms.

**Full theory.** Up to cubic order including self-interactions in the $\pi$-sector, the complete Lagrangian for the coupled $\pi$-$\sigma$ system that we consider is then

$$ \begin{aligned} \mathcal{L}/a^3 = & -\frac{1}{2} \left[ -\dot{\pi}_c^2 + c_s^2 \frac{(\partial_i \pi_c)^2}{a^2} \right] - \frac{1}{2} \left[ (\partial_\mu \sigma)^2 + m^2 \sigma^2 \right] + \rho \dot{\pi}_c \sigma \\ & - \lambda_1 \dot{\pi}_c \frac{(\partial_i \pi_c)^2}{a^2} - \lambda_2 \dot{\pi}_c^3 - \mu \sigma^3 - \frac{1}{2} \alpha \dot{\pi}_c \sigma^2 - \frac{1}{2\Lambda_1} \frac{(\partial_i \pi_c)^2}{a^2} \sigma - \frac{1}{2\Lambda_2} \dot{\pi}_c^2 \sigma, \end{aligned} \quad (4.6) $$

where we have defined the canonically normalised field $\pi_c^2 = c_s^{-3} f_\pi^4 \pi^2$ with $f_\pi^4 \equiv 2M_{\text{pl}}^2 |\dot{H}| c_s$ being the symmetry breaking scale,<fnref n="18" /> and redefined the coupling constants

$$ \begin{aligned} \rho &\equiv 2c_s^{3/2} \frac{\tilde{M}_1^3}{f_\pi^2}, & \lambda_1 &\equiv c_s^{3/2} \frac{c_s^2 - 1}{2f_\pi^2}, & \lambda_2 &\equiv -\frac{c_s^{3/2}}{2f_\pi^2} \left[ 1 - c_s^2 - \frac{8}{3} \frac{M_3^4}{f_\pi^4} c_s^3 \right], \\ \alpha &\equiv -4c_s^{3/2} \frac{H^2}{f_\pi^2} \tilde{\alpha}^2, & \Lambda_1^{-1} &\equiv c_s^{3/2} \frac{\rho}{f_\pi^2}, & \Lambda_2^{-1} &\equiv c_s^{3/2} \frac{\tilde{\rho}}{f_\pi^2}, \end{aligned} \quad (4.7) $$

with $\tilde{\alpha} = \frac{\tilde{M}_2}{H}$ and $\tilde{\rho} = 8c_s^{3/2} \frac{\tilde{M}_3^3}{f_\pi^2} - \rho$. These couplings have mass dimensions $[\Lambda_{1,2}] = [\mu] = [\rho] = 1$, $[\alpha] = 0$, and $[\lambda_{1,2}] = -2$. As we will see in Section 6, the couplings $\rho$, $\tilde{\rho}$ and $\tilde{\alpha}$ enter in the genuine typical size of each contribution to the dimensionless shape function of the bispectrum in the case $c_s = 1$. Both $\rho$ and $\Lambda_1$ are correlated since they are determined by the same parameter $\tilde{M}_1$, which is a consequence of the non-linearly realised time-diffeomorphism invariance. Similarly, $M_2$ governs both the Goldstone boson sound speed $c_s$ and $\lambda_1$. Note that there is a priori no model-building requirement on the size of the quadratic coupling $\rho/H$, and we crucially allow this mixing parameter to be order one or larger, which from now on we call the strong mixing regime.

**Cosmological flow implementation.** In order to implement the flow equations for the fields $\varphi^\alpha = (\pi_c, \sigma)$, one needs to write the theory in a Hamiltonian form. The rescaled linear conjugate

<footnote n="17">¹⁷The most general mixing operators are of the form $(\delta g^{00})^m O_n$ where $O_n$ is an operator of dimension $n$ function of the following combinations: $\sigma$, $g^{0\mu}\partial_\mu\sigma$ and $g^{\mu\nu}\partial_\mu\partial_\nu\sigma$. These operators have been fully classified in [147] up to dimension-6 operators. We consider an example of such operators in Section 7.2. Note that we do not impose additional symmetries on the $\sigma$ field to make it naturally light because we will be mostly interested in the cases where $\sigma$ is (effectively) heavier than the Hubble scale. Theories with additional symmetries imposed on the $\sigma$-sector have been treated in [140].</footnote>

<footnote n="18">¹⁸For energies above $f_\pi$, one should consider an infinite number of operators in the Lagrangian because higher-order terms—that we truncated—become relevant. At energies below $f_\pi$, a description in terms of Goldstone boson is appropriate.</footnote>
momenta associated with the fields are

$$p_{\pi_c} \equiv \frac{1}{a^3} \frac{\delta S}{\delta \dot{\pi}_c} = \dot{\pi}_c + \rho \sigma, \quad p_{\sigma} \equiv \frac{1}{a^3} \frac{\delta S}{\delta \dot{\sigma}} = \dot{\sigma}. \qquad (4.8)$$

After performing the Legendre transform, the Hamiltonian can be arranged in the form (3.6) and the identification of the various tensors yields

$$\begin{align}
\Delta_{\alpha\beta} &= \delta_{\alpha\beta}, & M_{\alpha\beta} &= \begin{pmatrix} -c_s^2 \frac{k^2}{a^2} & 0 \\ 0 & -(\frac{k^2}{a^2} + m^2 + \rho^2) \end{pmatrix}, & I_{\alpha\beta} &= \begin{pmatrix} 0 & \rho \\ 0 & 0 \end{pmatrix}, \nonumber \\
A_{\sigma\sigma\sigma} &= -2\mu + \alpha\rho - \frac{\rho^2}{\Lambda_2} + 2\lambda_2\rho^3, & A_{\pi_c\pi_c\sigma} &= \left(\frac{1}{\Lambda_1} - 2\lambda_1\rho\right) \frac{\mathbf{k}_{\pi_c} \cdot \mathbf{k}_{\pi_c}}{a^2}, \tag{4.9} \\
B_{\sigma\sigma p_{\pi_c}} &= -\alpha + 2\frac{\rho}{\Lambda_2} - 6\lambda_2\rho^2, & B_{\pi_c\pi_c p_{\pi_c}} &= 2\lambda_1 \frac{\mathbf{k}_{\pi_c} \cdot \mathbf{k}_{\pi_c}}{a^2}, \nonumber \\
C_{p_{\pi_c} p_{\pi_c} \sigma} &= -\frac{1}{\Lambda_2} + 6\lambda_2\rho, & D_{p_{\pi_c} p_{\pi_c} p_{\pi_c}} &= -2\lambda_2. \nonumber
\end{align}$$

Deriving these tensors fully fixes the flow equations and the initial conditions for the two- and three-point functions.

**From $\pi$ to $\zeta$.** One needs to relate the various correlators to primordial observables, namely correlators of the curvature perturbation $\zeta$. Actually, this link is only geometrical and can be expressed simply in terms of the Goldstone boson $\pi$, independently of the specific theory that one considers. Note that the link between $\pi$ and the canonically normalised field $\pi_c$ is theory-dependent. Our aim is only to compute the primordial power spectrum and bispectrum, i.e. the late-time two- and three-point correlators of $\zeta$ at the end of inflation, when all modes are well outside the horizon. In this regime, the relation between $\pi$ and $\zeta$ reads (see e.g. [148])

$$\zeta(t, \boldsymbol{x}) = \log \left( \frac{a(t + \xi(t, \boldsymbol{x}))}{a(t)} \right) \quad \text{with} \quad \pi(t + \xi(t, \boldsymbol{x}), \boldsymbol{x}) + \xi(t, \boldsymbol{x}) = 0, \qquad (4.10)$$

where the implicit equation defining $\xi$ is found by working out the time-diffeomorphism mapping the spatially flat gauge to the unitary one. Beyond linear order, a spatial diffeomorphism is also required for the mapping, but its effect is negligible when all modes are well outside the horizon. Using Eq. (4.10), one obtains

$$\zeta = -H\pi + H\pi\dot{\pi} + \frac{1}{2}\dot{H}\pi^2 + \mathcal{O}(\pi^3), \qquad (4.11)$$

up to quadratic order, which is sufficient to work out the power spectrum and bispectrum. Explicitly, one finds

$$\begin{align}
\langle \zeta_k \zeta_{-k} \rangle' &= H^2 \langle \pi_k \pi_{-k} \rangle', \nonumber \\
\langle \zeta_{k_1} \zeta_{k_2} \zeta_{k_3} \rangle' &= -H^3 \langle \pi_{k_1} \pi_{k_2} \pi_{k_3} \rangle' + H^3 [\Sigma^{\pi\pi}(k_1)\Sigma_{\mathcal{R}}^{\pi\pi}(k_2) + 5 \text{ perms}] \tag{4.12} \\
&\quad + H^2 \dot{H} [\Sigma^{\pi\pi}(k_1)\Sigma^{\pi\pi}(k_2) + 2 \text{ perms}], \nonumber
\end{align}$$
where all correlators should be evaluated at the end of inflation. A few comments are in order. First, the result (4.12) is valid in any theory, be it single-clock or with multiple degrees of freedom. Second, no decoupling limit or slow-roll approximation has been considered here, only that modes are well outside the horizon. Note also that even in single-clock inflation, $\pi$ is strictly speaking not massless and is time-dependent even outside the horizon. However, the time dependence of $H$, $\pi$ and $\dot{\pi}$ are such that the correlators of $\zeta$ themselves become time-independent in that case.

Eventually, in the applications considered in this work, where any time-dependence of $H$ is neglected, only the first line is relevant in the expression of the bispectrum. Notice also that the surviving correction to the bispectrum of $\pi$ can be rewritten in terms of correlators of fields and conjugate momenta, see Eq. (4.8). This shift of perspective is useful as even when multi-field interactions keep building up on super-horizon scales, the cross power spectrum involving the momentum of $\pi$ quickly decays and only remains the contribution from the mixed $\pi-\sigma$ power spectrum, see [142] for a related discussion in the context of the contribution of boundary terms in multi-field inflation. This situation occurs only when the field $\sigma$ is almost massless. However, in this situation, $\zeta$ keeps evolving until the end of inflation, which effectively breaks scale invariance and makes the theory non-predictive without a complete description of reheating. We discuss this particular case in Appendix B, see also [149] for a bootstrap perspective on this parameter space. When $\sigma$ is effectively heavy, it quickly decays on super-horizon scales, an adiabatic limit is then reached by the end of inflation, and the bispectrum of $\zeta$ is directly proportional to the one of $\pi$. The large majority of our applications lie in this regime.

In the single-field limit $\rho = 0$, the amplitude of curvature perturbations $\zeta$ is fixed by the ratio of $H$ and $f_\pi$

$$ \Delta_{\zeta,0}^2 = \frac{k^3}{2\pi^2} \langle \zeta_k \zeta_{-k} \rangle' = \frac{1}{4\pi^2} \left( \frac{H}{f_\pi} \right)^4 . \qquad (4.13) $$

The measured amplitude of the dimensionless power spectrum on CMB scales is $\Delta_{\zeta}^2 = 2.2 \times 10^{-9}$ [32].

## 4.2 Quadratic Theory

We begin the study of this theory by examining the time evolution of the various two-point functions. In this section, we start by introducing the notion of effective mass that emerges from a sizeable quadratic mixing. Then, we provide exact results—derived for the first time using the cosmological flow—for the power spectrum in all regimes of interest including the strong mixing regime. This provides a comprehensive analysis of the entire quadratic phase space.

### 4.2.1 Effective mass and resonance at strong mixing

Writing the theory in the Hamiltonian form makes it explicit that the quadratic mixing induces an effective mass for the $\sigma$ field

$$ m_{\text{eff}}^2 = m^2 + \rho^2, \qquad (4.14) $$

that can be appreciated from the form of the $M$ matrix in Eq. (4.9). This is a consequence of considering the full linear momenta, i.e. the conjugate momentum of $\pi_c$ is corrected due to the
quadratic mixing.

**Late-time leading behaviour.** To understand the origin of the effective mass at energy scales $\omega \ll H$, it is sufficient to look at the leading late-time behaviour of the $\pi_c - \sigma$ system. Looking for solutions of the form (see Appendix B)

$$ \lim_{\tau \to 0} \begin{pmatrix} \pi_c(\tau, \mathbf{x}) \\ \sigma(\tau, \mathbf{x}) \end{pmatrix} = \begin{pmatrix} \pi_c^+(\mathbf{x}) \\ \sigma^+(\mathbf{x}) \end{pmatrix} \tau^{\Delta_+} + \begin{pmatrix} \pi_c^-(\mathbf{x}) \\ \sigma^-(\mathbf{x}) \end{pmatrix} \tau^{\Delta_-}, \qquad (4.15) $$

where $\tau$ is the conformal time, one finds four modes with scaling dimensions

$$ \Delta_{\perp}^{\pi_c} = 0, \quad \Delta_{+}^{\pi_c} = 3, \quad \Delta_{\pm}^{\sigma} = \frac{3}{2} \pm i\mu_{\text{eff}}, \qquad (4.16) $$

with $\mu_{\text{eff}}^2 = m_{\text{eff}}^2/H^2 - 9/4$, which makes it manifest that the system is composed of two modes corresponding to a massless field ($\Delta_{\pm}^{\pi_c}$) and a decaying massive field ($\Delta_{\pm}^{\sigma}$) with effective mass $m_{\text{eff}}$.

**Flat-space leading behaviour.** At sufficient high energies $\omega \gg H$,<fnref n="19" /> the dynamics of the theory can be approximated by the flat-space limit $H \to 0$, in which case the system can be solved analytically up to Hubble friction corrections.<fnref n="20" /> Here, we explicitly show that, even if the quadratic mixing cannot a priori be treated perturbatively, the analytical solution of the quadratic theory follows from resumming an infinite number of quadratic mixing insertions in the propagators of the two degrees of freedom. Focusing on the field $\sigma$, its dynamics is affected by the surrounding $\pi_c$ medium that interacts with it, leading to a self-energy correction. At tree-level as we only consider quadratic interactions, it reads

$$ \Pi(\omega, \mathbf{k}) = \begin{array}{c} \bullet \quad \bullet \\ \text{---} \\ \text{---} \\ \text{---} \\ \text{---} \end{array} = \frac{-i\rho^2\omega^2}{\omega^2 - c_s^2k^2 + i\epsilon}, \qquad (4.18) $$

that is found from the following elementary Feynman rules

$$ G_{\pi_c}^0(\omega, \mathbf{k}) = \begin{array}{c} \bullet \quad \bullet \\ \text{---} \\ \text{---} \\ \text{---} \\ \text{---} \end{array} = \frac{i}{\omega^2 - c_s^2 k^2 + i\epsilon}, \qquad \bullet \quad \bullet \\ \text{---} \end{array} = -i\rho\omega, \qquad (4.19) $$

where $G_{\pi_c}^0$ (resp. $G_{\sigma}^0$) is the free $\pi_c$ (resp. $\sigma$) Green two-point function. From the bare propagator, the dressed propagator $G_{\sigma}(\omega, \mathbf{k})$ is found by computing the associated Dyson series summing over all one-particle irreducible diagrams—in our case the tree-level contribution. Summing the

$$ \langle \text{T}(\varphi_{\mathbf{k}}^T \varphi_{-\mathbf{k}}) \rangle = \frac{i}{(\omega^2 - c_s^2 k^2)(\omega^2 - k^2 - m^2) - \rho^2\omega^2 + i\epsilon} \begin{pmatrix} \omega^2 - k^2 - m^2 & -i\rho\omega \\ i\rho\omega & \omega^2 - c_s^2 k^2 \end{pmatrix}. \qquad (4.17) $$

<footnote n="19">¹⁹In the strong mixing regime $\rho \gg H$, the flat-space limit is reached for energies such that $\omega \gg \rho$.</footnote>

<footnote n="20">²⁰By looking at plane-wave solutions of the form $\pi_c \propto e^{i\omega t}$ and $\sigma \propto e^{i\omega t}$, the Feynman propagator for the fields $\varphi^\alpha = (\pi_c, \sigma)$ is found by inverting the quadratic operator</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Time evolution of the rescaled σ two-point correlation function in Hubble units for a fixed effective mass such that $\mu_{\text{eff}} = 1.5$ and $c_s = 1$, varying the bare mass $m$. The dotted black line corresponds to the analytical prediction in the decoupled regime $\rho/H = 0$. The vertical dotted line shows the energy scale at which the correlators cross the horizon $-k\tau = 1$.</figcaption></figure>

geometrical series, we have

$$ 
\begin{aligned}
G_{\sigma}(\omega, \mathbf{k}) &= \text{---} + \text{---} + \text{---} + \text{---} + \dots \\
&= G_{\sigma}^{0}(\omega, \mathbf{k}) \sum_{n=0}^{\infty} [\Pi(\omega, \mathbf{k}) G_{\sigma}^{0}(\omega, \mathbf{k})]^{n} = \frac{i}{\omega^{2} - k^{2} - m^{2} - \rho^{2} \left[ \frac{\omega^{2}}{\omega^{2} - c_{s}^{2} k^{2}} \right] + i\epsilon},
\end{aligned} 
\quad (4.20) 
$$

which agrees with (4.17). Non-trivial plane-wave solutions for the propagating modes are found at the poles of the dressed propagator

$$ \omega_{\pm}^{2} = \frac{1+c_{s}^{2}}{2}k^{2} + \frac{m_{\text{eff}}^{2}}{2} \pm \frac{1}{2}\sqrt{m_{\text{eff}}^{4} + (1-c_{s}^{2})^{2}k^{4} + 2(1+c_{s}^{2})k^{2}\left[\rho^{2} + \frac{1-c_{s}^{2}}{1+c_{s}^{2}}m^{2}\right]}. \quad (4.21) $$

We wish to highlight several features of this result by examining various limits. At very high momenta $k \gg \rho, m$, the dispersion relation boils down to $\omega_{-}^{2} \approx c_{s}^{2}k^{2}$ and $\omega_{+} \approx k^{2}$ which describes two independent massless free fields, i.e. the vacuum state in the asymptotic past limit of the inflationary spacetime. In the strong mixing regime $\rho \gg k, m$, one finds

$$ \omega_{-}^{2} \approx \frac{k^{2}}{\rho^{2}}(c_{s}^{2}k^{2} + m^{2}), \quad \omega_{+}^{2} \approx m_{\text{eff}}^{2} \approx \rho^{2}, \quad (4.22) $$

which describes a massless degree of freedom and a massive degree of freedom with mass $m_{\text{eff}}^{2} = m^{2} + \rho^{2}$. We note that in the UV limit $k \gg m$, the massless degree of freedom has a modified dispersion relation $\omega \propto k^{2}$. In the rest of this section, we set $c_{s} = 1$ to avoid clutter.

**Resonance at strong mixing.** The appearance of the effective mass both in the deep UV and in the late IR suggests that a similar resummation in the intermediate regime—around
horizon crossing—is compelling. Figure 2 shows the time evolution of the $\sigma$ two-point correlation function rescaled by $a^3$ for a fixed effective mass $\mu_{\text{eff}} = 1.5$ while varying the bare mass (hence also adjusting the quadratic mixing). At late times $\omega \ll H$, as explained by the previous late-time analysis, one can appreciate that the correlators oscillate at the same frequency set by $\mu_{\text{eff}}$. The amplitude of the oscillations and the mean value of the two-point function grow as the bare mass $m$ decreases. Furthermore, another notable feature is that the oscillations present a phase shift as the quadratic mixing constant $\rho/H$ increases. This is due to a noticeable phenomenon, namely the amplification of the massive correlator soon after Hubble crossing, which induces a time delay before the late-time asymptotic regime is reached. This effect can be explained by examining the strong mixing regime of the $\pi_c - \sigma$ dynamics:<fnref n="21" />

(i) At the effective mass-shell horizon $|k\tau| \sim \mu_{\text{eff}}$, the massive field starts to decay as $a^{-3/2}$,

(ii) At the $\rho$-horizon crossing $|k\tau| \sim \sqrt{\rho/H}$, the Goldstone boson freezes, as we will see in the next section.

At strong quadratic mixing, the event (i) occurs before (ii).<fnref n="22" /> Therefore, when the massive field starts to decay, it transfers its power to the Goldstone boson, creating the $\sqrt{\rho/H}$ amplification of the power spectrum (see next section). Because the dynamics is dominated by the quadratic mixing, the Goldstone boson backreacts on the massive field after horizon crossing $\omega = H$, creating the visible boost in Figure 2. The precise dynamics can only be accessed numerically. This interesting phenomenon is also visible in the cosmological collider signal, as we will see in Section 7.

### 4.2.2 Phase diagram

Fixing $c_s = 1$, the quadratic theory depends on two variables $m$ and $\rho$. The entire phase space is not accessible analytically. In Figure 4, we show the dimensionless power spectrum as a function of $\rho/H$ and $m/H$, computed by solving the flow equations. This result is exact and is consistent with previous studies [150,151]. To get some physical intuition, one can understand the qualitative behaviour of the quadratic theory in some asymptotic limits that we summarise in Figure 3. In the following, we give insights on these regimes.

**Modified dispersion relation.** When the mixing is sufficiently strong, the quadratic interaction dominates the dynamics, and one can drop the kinetic terms $\dot{\pi}_c^2$ and $\dot{\sigma}^2$ in (4.6). The equation of motion for the massive field in Fourier space becomes $(k^2/a^2 + m^2)\sigma = \rho\dot{\pi}_c$, which, if we also neglect the mass term $m^2 \ll k^2/a^2$ at relevant times, simplifies to

$$ \frac{k^2}{a^2} \sigma = \rho \dot{\pi}_c. \qquad (4.23) $$

<footnote n="21">²¹We recall that we have set $c_s = 1$ so that the following arguments are valid for a unity intrinsic sound speed only.</footnote>

<footnote n="22">²²The reverse case, i.e. $\sqrt{\rho/H} > \mu_{\text{eff}}$, is reached only in a very limited region of the parameter space.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: Phase diagram of the quadratic theory in the $(m/H, \rho/H)$ plane, for $c_s = 1$, that delimits three distinct regions: (i) the modified dispersion relation regime ($M \ll H$ and $\rho \gg H$ where $M = m^2/\rho$), (ii) the reduced speed of sound regime ($m \gg H/\tilde{c}_s$ with $\tilde{c}_s^{-2} = 1 + \rho^2/m^2$), and (iii) the light field regime ($m_{\text{eff}} \ll H$). One should notice that the regime with a significantly reduced speed of sound, $\tilde{c}_s \ll 1$, is rather limited as it requires $\rho \gg m$ in region (ii). This is shown by the grey dashed line $\rho = m$. Of course, the borderlines of these regions should not be taken as exact but rather as faint boundaries: the deeper in the regions, the better the approximations. More details on each regime are given below.</figcaption></figure>

In this case, the equation of motion for the Goldstone boson is found to be

$$ \ddot{\pi}_c + 5H\dot{\pi}_c + \frac{k^4}{\rho^2 a^4} \pi_c = 0. \qquad (4.24) $$

Therefore, the theory can be described by an effective single-field theory—because the field $\sigma$ is non-propagating i.e. its equation of motion is a pure constraint—with a modified (quadratic) dispersion relation $\omega = k^2/\rho$, inferred from (4.24). Let us now derive the precise conditions of validity of this regime. First, we have considered that the quadratic mixing $\rho\dot{\pi}_c\sigma$ dominates over the kinetic term $\dot{\sigma}^2$. Using the constraint (4.23) and the dispersion relation, this condition translates into the following bound

$$ \frac{\rho}{H} \gtrsim 1. \qquad (4.25) $$

Imposing that the mixing term also dominates the kinetic term $\dot{\pi}_c^2$ leads to the same bound.<fnref n="23" /> Similarly, neglecting the mass term $m^2\sigma^2$ leads to<fnref n="24" />

$$ M \equiv \frac{m^2}{\rho} \lesssim H, \qquad (4.26) $$

where we see the appearance of a new energy scale $M$. These two bounds define the modified dispersion relation regime, as depicted in Figure 3. In this regime, the field $\sigma$ plays the role of the conjugate momentum of the massless field $p_{\pi_c} \approx \rho\sigma$, and the effective single-field quadratic

<footnote n="23">²³We work here in the limit $c_s = 1$. However, in the general case, the most stringent bound is instead found to be $\rho/H \gtrsim 1/c_s$ [64].</footnote>

<footnote n="24">²⁴In the general case with $c_s < 1$, the bound is found to be $c_s m^2 / \rho \lesssim H$ [64].</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Right panel: Phase diagram of the dimensionless power spectrum $\Delta_{\zeta}^2$ (in units of $\Delta_{\zeta,0}^2$). We have set $N = 60$ e-folds of super-horizon evolution. The dashed line corresponds to the $M = m^2/\rho = H$ delimitation, isolating the modified dispersion relation regime ($M \ll H$ and $\rho \gg H$) and the reduced speed of sound regime ($M \gg H$) (see Figure 3 for more details). Left panel: Slices of the phase diagram for $\rho/H = 10$, $\rho/H = 1$, and $\rho/H = 0.1$ as a function of the bare mass $m/H$ in Hubble units. The grey dots mark the $M/H = 1$ delimitation. Note that the quadratic interaction between $\pi_c$ and $\sigma$ alters the symmetry breaking scale $f_{\pi} \approx 58H \times (\Delta_{\zeta}^2/\Delta_{\zeta,0}^2)^{1/4}$.</figcaption></figure>

Lagrangian becomes

$$ \mathcal{L}^{(2)}/a^3 \approx \rho \dot{\pi}_c \sigma - \frac{1}{2} \frac{(\partial_i \pi_c)^2}{a^2} - \frac{1}{2} \frac{(\partial_i \sigma)^2}{a^2}, \qquad (4.27) $$

which leads to the equation of motion (4.24) after using (4.23). This can be solved exactly, yielding the following mode function for the light field after imposing the Bunch-Davies vacuum [27]

$$ \pi_{c,k}(\tau) = \sqrt{\frac{\pi}{4} \frac{H}{\rho} \frac{H}{\sqrt{2k^3}}} (-k\tau)^{5/2} H_{5/4}^{(1)} \left( \frac{1}{2} \frac{H}{\rho} (k\tau)^2 \right) \xrightarrow{k\tau \to 0} e^{-i\frac{\pi}{2}} \frac{2\Gamma(5/4)}{\sqrt{\pi}} \frac{H}{k^{3/2}} \left( \frac{\rho}{H} \right)^{1/4}. \quad (4.28) $$

From the analytical expression of the mode function, or simply evaluating the non-linear dispersion relation at $\omega \sim H$, one observes the emergence of a new time scale in this regime, given by $|k\tau| \sim \sqrt{\rho/H} \gtrsim 1$, hence the $\pi_c$ mode freezes before the usual horizon crossing $|k\tau| \sim 1$. This regime leads to the amplification of the power spectrum

$$ \Delta_{\zeta}^2 = \frac{2\Gamma(5/4)^2}{\pi^3} \left(\frac{H}{f_{\pi}}\right)^4 \left(\frac{\rho}{H}\right)^{1/2}. \qquad (4.29) $$

We numerically confirm this expected asymptotic behaviour as $M \ll H$, and the fact that the power spectrum becomes independent of $m/H$, hence reaching a plateau in the left panel of Figure 4 for $\rho/H = 10$.

**Reduced speed of sound.** When the bare mass of $\sigma$ becomes large $m \gg H$ while maintaining $M \gg H$, the field $\sigma$ can be integrated out by setting it on-shell$^{25}$ ($-\Box + m^2)\sigma = \rho\dot{\pi}_c$ [152, 153].

<footnote>$^{25}$We will come back to the regime of validity of this regime and give more details in Section 8.1.</footnote>
The quadratic Lagrangian then inherits the imprint of this heavy mode in an infinite number of irrelevant operators which—at leading order—induces the well-known speed of sound<fnref n="26" /> for $\pi_c$

$$ \tilde{c}_s^{-2} = 1 + \frac{\rho^2}{m^2}. \qquad (4.30) $$

In this effective theory, the massless mode crosses the sound horizon $|\tilde{c}_s k \tau| \sim 1$ and freezes while its dispersion relation is linear $\omega = \tilde{c}_s k$. The power spectrum becomes

$$ \Delta_{\zeta}^2 = \frac{1}{4\pi^2 \tilde{c}_s} \left(\frac{H}{f_{\pi}}\right)^4. \qquad (4.31) $$

We will see in Section 8.1 that the actual regime of validity of this regime is $m \gg H/\tilde{c}_s$. However, as depicted in Figure 3, the actual region in the phase diagram where $\tilde{c}_s \ll 1$, which typically leads to enhanced equilateral non-Gaussianities, is limited.

Although one might naively think that $\rho \lesssim H$ is the correct condition on the quadratic coupling to have weak mixing, it has been shown that this condition can be extended to $\rho \lesssim m$ for $m \gtrsim H$, see [60, 64]. In this regime, $\rho < m$, that we will from now on refer to as the weak mixing, one finds that the correction to the power spectrum scales as $\rho^2/m^2$. This correction has also been computed in [35, 154] in the full two-field theory for $c_s = 1$, see also [40] for numerical and [60] for analytical computations in the case of an intrinsic reduced speed of sound $c_s < 1$. At strong mixing $\rho \gtrsim m$, the speed of sound is significantly reduced and the power spectrum scales as $\rho/m$.

The mass dependency of the correction, scaling as $(H/m)^2$, reflects the imprint of the heavy field $\sigma$. Note that it also induces higher-dimensional operators suppressed by additional powers of $(H/m)^2$. We numerically confirm this mass dependence and that this scaling extends to would-be non-perturbative values of the quadratic mixing $\rho/H > 1$. In fact, we find that the perturbative regime holds up to $\rho \sim m > H$, as can be understood by the emergence of the effective mass (4.14) when going beyond perturbative treatment of the quadratic mixing, and as also shown by the regime of validity of the expansion in (4.31).

We will show in Section 6 that these two regimes, of a modified dispersion relation and of a reduced speed of sound, can be captured by a unified effective theory where the heavy field is integrated out in a (spatially) non-local manner.

**Light additional field.** Our numerical analysis is also applicable in the light $\sigma$ field regime where $m_{\text{eff}} \lesssim H$. In this regime, the field $\sigma$ decays at a very slow pace on super-horizon scales. The massless field $\pi_c$ is therefore sourced by $\sigma$ on super-horizon scales so that the value of $\Delta_{\zeta}^2$ is sensitive to the amount of $e$-folds elapsed from the moment at which the mode crosses the horizon and the end of inflation. The corresponding enhancement of the power spectrum is well visible in the lower left part of the phase diagram in Figure 4. In Appendix B, we derive the

<footnote n="26">²⁶The Goldstone boson can in principle have an additional speed of sound $c_s$ governed by the coefficient $M_2$ in (4.1). For simplicity, this intrinsic speed of sound has been set to unity in this section. However, note that the total speed of sound can be further reduced when the heavy field is integrated out: $\tilde{c}_s^{-2} = c_s^{-2}(1 + \rho^2/m^2)$.</footnote>
following analytical approximation that reproduces this amplification

$$ \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, \qquad (4.32) $$

where $\Delta N$ is the number of e-folds elapsed since horizon crossing.

# 5 Strong Coupling Scales and Naturalness

It is important to ensure that the effective description of fluctuations that we study is under theoretical control and identify the natural hierarchy of scales of the system. Because the theory contains non-renormalisable interactions, the theory will become strongly coupled at a certain energy scale.<fnref n="27" /> In this section, we derive the strong coupling scale, or unitarity bound, in the weak and strong quadratic mixing regimes. We gather all the found bounds in Appendix C. These bounds translate into bounds on the size of non-Gaussianities, see Section 6, but only give a parametric estimation, possibly missing large dimensionless numbers.

**Goldstone boson sector.** The effective description of primordial fluctuations in (4.6) is valid up to some UV cutoff scale $\Lambda_*$, called the strong coupling scale, beyond which the perturbative description breaks down. When $c_s = 1$, analogous to the chiral Lagrangian describing the dynamics of pions, the strong coupling scale of the effective theory is $\Lambda_* = 4\pi f_\pi$.<fnref n="28" /> However, when the sound speed $c_s$ is sufficiently small, it is well-known that $\Lambda_*$ is not too far from the Hubble scale. Indeed, it is given by [27, 139]<fnref n="29" />

$$ \Lambda_*^4 = 2\pi f_\pi^4 \frac{c_s^4}{1 - c_s^2}, \qquad (5.1) $$

which, in the limit $c_s \ll 1$, gives $\Lambda_* \sim c_s f_\pi$. Imposing $\Lambda_* > H$ for the EFT to be valid—otherwise fluctuations cannot be initialised in the proper vacuum inside the horizon—and assuming the single-field limit (i.e. $\rho/H = 0$) gives a theoretical lower bound on the sound speed $c_s \ge 0.01$.<fnref n="30" /> The additional degree of freedom $\sigma$ changes the physical description, and interactions between $\pi_c$ and $\sigma$ will modify $\Lambda_*$. Yet, the quadratic mixing $\dot{\pi}_c\sigma$ makes the analysis more subtle. Therefore, we examine separately the weak mixing regime and the modified dispersion relation regime—in the following loosely denoted as strong mixing—corresponding to the blue region in the phase diagram 3.

<footnote n="27">²⁷This scale can be derived by finding the maximum energy at which the tree-level scattering of the fields is unitary, through partial wave expansion [155]. Here, we give an heuristic derivation of these strong coupling scales based on power counting which, of course, misses order-one factors.</footnote>

<footnote n="28">²⁸After reintroducing the Goldstone boson, the action (4.2) without truncating to a certain order in fields is an expansion in $\pi_c/f_\pi$ with $\pi_c \sim \omega$, and higher dimensional operators become relevant if the effective dimensionless coupling constant $\omega/f_\pi$ is of order one.</footnote>

<footnote n="29">²⁹This strong coupling scale is associated with the operator $(\partial_i \pi_c)^4$.</footnote>

<footnote n="30">³⁰Assuming pure single-field $c_s$-theory, Planck constraints on primordial non-Gaussianities give a more stringent lower bound on the speed of sound: $c_s \ge 0.021$ (95% CL) [156].</footnote>
## 5.1 Weak Mixing

In the weak mixing regime, the quadratic theory is dominated by the standard kinetic terms. The free fields can then be considered decoupled as a first approximation. Usually, the strong coupling scale $\Lambda_*$ is derived from operators constructed out of canonically-normalised fields with a Lorentz invariant kinetic term. Here, the non-relativistic dispersion relation of $\pi_c$ makes it impossible to put time and space on the same footing, so that one needs to consider energy and momentum separately.<fnref n="31" />

Let us first consider the cubic interaction coming from $\delta g^{00}\sigma$ in the unitary gauge, so that the action we examine—together with the quadratic part but neglecting the quadratic mixing—reads

$$S \supset \frac{1}{2} \int dt d^3x a^3 \left[ \dot{\pi}_c^2 - c_s^2 \frac{(\partial_i \pi_c)^2}{a^2} + \dot{\sigma}^2 - \frac{(\partial_i \sigma)^2}{a^2} - m^2 \sigma^2 - c_s^{3/2} \frac{\rho}{f_\pi^2} \frac{(\partial_i \pi_c)^2}{a^2} \sigma \right]. \quad (5.2)$$

In order to estimate how various terms in this action scale with energy and momentum, we need to choose a reference. Let us fix the momentum $k$ carried by the propagating modes, whose energies should be understood as derived concepts depending on their dispersion relations. The energy of $\pi_c$ then reads $\omega_\pi = c_s k$, whereas the energy of $\sigma$ is $\omega_\sigma = k$, from which we deduce that $\omega_\pi = c_s \omega_\sigma$. In the following, we will express all quantities in terms of the energy of the Goldstone boson $\omega \equiv \omega_\pi$ because it is the energy probed experimentally at the Hubble scale $\omega \sim H$. Therefore, the kinetic term $\dot{\pi}_c^2$ scales as $\sim \omega^2 \pi_c^2$—so does the gradient term—, and the kinetic term $\dot{\sigma}^2$ scales as $\sim \omega_\sigma^2 \sigma^2 \sim c_s^{-2} \omega^2 \sigma^2$, and likewise for the gradient term. Yet, a question needs to be answered. What is the typical size of the fluctuations? By dimensional analysis, $\pi_c$ scales as $\omega^{-1/2} k^{3/2} \sim c_s^{-3/2} \omega$. Similarly, $\sigma$ scales as $\omega_\sigma \sim c_s^{-1} \omega$. Combining all scalings, we obtain that the Goldstone boson kinetic term scales as $\sim c_s^{-3} \omega^4$, and that of $\sigma$ scales as $c_s^{-4} \omega^4$. We notice that, at small sound speed, the kinetic term of $\sigma$ is the dominant one in the quadratic action.<fnref n="32" /> No strong coupling scale is associated with the quadratic mixing operator $\dot{\pi}_c \sigma$ because it is of dimension 3. However, the theory is under control in the weak mixing regime if this interaction can be treated perturbatively. Requiring that this interaction—that scales as $c_s^{-5/2} \rho \omega^3$—is smaller than the kinetic term of the Goldstone boson gives

$$\frac{\rho}{H} \lesssim c_s^{-1/2}. \qquad (5.3)$$

Now, the cubic interaction scales as $c_s^{-9/2} \rho/f_\pi^2 \omega^5$. Requiring that this term should be smaller than the kinetic term of the Goldstone boson—giving the more stringent bound—sets the strong

<footnote n="31">³¹In single-field theories, up to a redefinition of spacetime coordinates, fake Lorentz invariance can be restored and standard dimensional analysis can then be applied, see [27]. However, this trick does not work in our situation with several degrees of freedom propagating at different sound speeds.</footnote>

<footnote n="32">³²Note also that this hierarchy is reversed when the massive fluctuations propagate slower than the Goldstone boson fluctuations, $c_\sigma < c_\pi$, a situation discussed in Section 7.2.</footnote>
coupling scale to<fnref n="33" />

$$ \Lambda_{\star} = c_s^{3/2} \frac{f_{\pi}^2}{\rho}. \qquad (5.4) $$

Imposing $H \lesssim \Lambda_{\star}$ leads to $\rho/H \lesssim c_s^{3/2}/2\pi\Delta_{\zeta}$. This bound is less stringent than the previous one, and therefore (5.3) defines the weak mixing regime. Additionally, it has been shown that $\rho \lesssim m$ is also required in order to treat the quadratic mixing perturbatively, see Section 4.2.2. Depending on the value of $c_s$, this bound can be stronger.

Next, we consider the interactions coming from $\delta g^{00}\sigma^2$. After introducing the Goldstone boson, the interactions we examine are

$$ S \supset \frac{1}{2} \int dt \, d^3x \, a^3 \left[ -\alpha \dot{\pi}_c \sigma^2 + \frac{c_s^{3/2}}{2} \frac{\alpha}{f_{\pi}^2} (\partial_{\mu} \pi_c)^2 \sigma^2 \right]. \qquad (5.5) $$

The interaction $\dot{\pi}_c \sigma^2$ is of dimension 4 so no strong coupling scale is associated with it. However, requiring that the size of this term in the action be smaller than the kinetic term of $\pi_c$ to be in the perturbative regime leads<fnref n="34" /> to

$$ \alpha \lesssim c_s^{1/2}, \quad \text{i.e.} \quad \tilde{\alpha}^2 \lesssim c_s^{-1/4} \times (2\pi\Delta_{\zeta})^{-1}. \qquad (5.6) $$

Repeating the power counting for the dimension-6 operator gives the following strong coupling scale<fnref n="35" />

$$ \Lambda_{\star} = c_s^{5/4} \frac{f_{\pi}}{\sqrt{\alpha}}. \qquad (5.7) $$

Requiring that primordial fluctuations are weakly coupled at the Hubble scale $\Lambda_{\star} > H$ gives the constraint $\alpha \lesssim \frac{c_s^{5/4}}{2\pi\Delta_{\zeta}}$, equivalently $\tilde{\alpha}^2 \lesssim c_s^{-1/4} (4\pi\Delta_{\zeta})^{-2}$, on the coupling constant, which is less stringent than the previous one.

Lastly, we consider the interactions coming from $(\delta g^{00})^2 \sigma$, which in terms of the Goldstone

<footnote n="33">³³Let us comment on the mass term $m^2\sigma^2$, that we neglected so far. This term changes the dispersion relation of $\sigma$ to $\omega_{\sigma}^2 = k^2 + m^2$, which leads to $\omega_{\sigma}^2 = \omega^2/c_s^2 + m^2$. The kinetic term of $\sigma$ including this mass term therefore scales as $\sim (\omega^2/c_s^2 + m^2)\sigma^2$. The typical size of $\sigma$ is found to be $\sigma \sim \omega/c_s [1 + (c_s m/\omega)^2]^{-1/4}$. In the end, the kinetic term of $\sigma$ scales as $\sim c_s^{-4}\omega^4 [1 + (c_s m/\omega)^2]^{1/2}$. Note that we recover the scaling presented in the main text in the limit $\omega \gg c_s m$. Fixing the energy scale to the Hubble scale, we deduce that the mass term is negligible when $m/H \ll c_s^{-1}$, corresponding to the low-speed collider regime [60, 64]. Still, we note that properly taking into account the mass term increases the size of the kinetic term of $\sigma$. For $c_s < 1$ and including the mass of $\sigma$, the smallest contribution of the quadratic theory is still the kinetic term of the Goldstone boson. Therefore, we can safely discard the mass term in our analysis because it would play no role in the derivation of the strong coupling scales.</footnote>

<footnote n="34">³⁴In the standard Lorentz invariant case, the perturbative criteria on marginal operators is that the dimensionless couplings have to be less than unity. Here, the non-relativistic dispersion relation makes things more subtle, introducing a sound speed dependence in the perturbative criterion.</footnote>

<footnote n="35">³⁵Note that, for $c_s \ll 1$, the gradient term in $(\partial_{\mu}\pi_c)^2$ dominates over the $\dot{\pi}_c^2$ term, hence giving the strong coupling scale by requiring that it should be smaller than the kinetic term of $\pi_c$. However, if $c_s > 1$, the hierarchy is reversed and the strong coupling scale is found by comparing the time derivative interaction term to the kinetic term of $\sigma$. This gives $\Lambda_{\star} = c_s^{-1/4} f_{\pi}/\sqrt{\alpha}$.</footnote>
boson gives the following interactions

$$ S \supset \frac{1}{2} \int dt \, d^3 x \, a^3 \left[ -\frac{1}{\Lambda_2} \frac{\dot{\pi}_c^2}{\sigma} + \frac{c_s^{3/2}}{f_\pi^2} \frac{\dot{\pi}_c}{\Lambda_2} \frac{(\partial_\mu \pi_c)^2}{\sigma} - \frac{c_s^3}{4f_\pi^4} \frac{1}{\Lambda_2} \frac{\dot{\pi}_c^4}{\sigma} + \frac{c_s^3}{2f_\pi^4} \frac{1}{\Lambda_2} \frac{\dot{\pi}_c^2}{a^2} \frac{(\partial_i \pi_c)^2}{\sigma} - \frac{c_s^3}{4f_\pi^4} \frac{1}{\Lambda_2} \frac{(\partial_i \pi_c)^4}{a^4} \sigma \right], \quad (5.8) $$

where we have defined $\Lambda_2^{-1} = 8c_s^3 \tilde{M}_3^3 / f_\pi^4$. Note that this definition differs from (4.7). The reason is that the interaction $\dot{\pi}_c^2 \sigma$ is also generated by $\delta g^{00} \sigma$, so that its coupling constant depends on both $\tilde{M}_1$ and $\tilde{M}_3$. Here, we isolate interactions coming from $(\delta g^{00})^2 \sigma$ only, effectively setting $\tilde{M}_1 = 0$. Estimating the size of all the interactions in (5.8) gives several strong coupling scales. The smallest one is associated with the interaction $\dot{\pi}_c^2 \sigma$ and reads

$$ \Lambda_{\star} = c_s \Lambda_2. \quad (5.9) $$

Requiring $H < \Lambda_{\star}$ gives the bound $H/\Lambda_2 \lesssim c_s$, equivalently $\tilde{\rho}/H \lesssim c_s^{-1/2} (2\pi\Delta_{\zeta})^{-1}$.

Finally, the operator $\sigma^3$ being relevant and the field $\sigma$ having a relativistic dispersion relation, the criteria for a perturbative treatment of this interaction is simply $\mu/H \lesssim 1$.

## 5.2 Strong Mixing

At large quadratic mixing, as discussed in Section 4.2, the mixing term $\dot{\pi}_c \sigma$ dominates over the kinetic terms $\dot{\pi}_c^2$ and $\dot{\sigma}^2$. We recall that, in the general case with $c_s \le 1$, the condition of validity for this regime is [64]

$$ \frac{\rho}{H} \gtrsim c_s^{-1}, \quad \text{and} \quad c_s \frac{m^2}{\rho} \lesssim H. \quad (5.10) $$

The theory being effectively described by a single degree of freedom, it is possible to rescale the spatial coordinates and the fields to canonically normalise the propagating degree of freedom. This is completely equivalent to the power counting method used previously. By performing the rescaling $x \to \tilde{x} = \rho^{1/2} c_s^{-1/2} x$, $\pi_c \to \tilde{\pi}_c = \rho^{-1/4} c_s^{5/4} \pi_c$, and $\sigma \to \tilde{\sigma} = \rho^{-1/4} c_s^{1/4} \sigma$, we obtain for the interactions coming from $\delta g^{00} \sigma$—and the quadratic part but neglecting the time derivative terms—the following action

$$ S \supset \frac{1}{2} \int dt \, d^3 x \, a^3 \left[ \dot{\tilde{\pi}}_c \tilde{\sigma} - \frac{(\tilde{\partial}_i \tilde{\pi}_c)^2}{a^2} - \frac{(\tilde{\partial}_i \tilde{\sigma})^2}{a^2} - M^2 \tilde{\sigma}^2 - c_s^{-3/4} \frac{\rho^{5/4}}{f_\pi^2} \frac{(\tilde{\partial}_i \tilde{\pi}_c)^2}{a^2} \tilde{\sigma} \right], \quad (5.11) $$

where $M \equiv c_s \frac{m^2}{\rho}$. By dimensional analysis, we have $[\tilde{\partial}] = 1/2$ and $[\tilde{\pi}_c] = [\tilde{\sigma}] = 3/4$ in mass units. This implies that the strong coupling scale is

$$ \Lambda_{\star} = c_s \left( \frac{f_{\pi}^8}{\rho^5} \right)^{1/3}, \quad (5.12) $$
which generalises the strong coupling scale found in [27] for the $c_s = 1$ case. Using (4.29) and imposing $H \lesssim \Lambda_*$, we obtain the following constraint$^{36}$

$$ \frac{\rho}{H} \lesssim c_s^{3/4} \frac{\kappa^{1/2}}{\Delta_{\zeta}}, \qquad (5.13) $$

where $\kappa = 2\Gamma(5/4)^2/\pi^3 \approx 0.053$. Importantly, this bound shows that the theory can be strongly mixed without being strongly coupled.

Let us now consider the interactions coming from $\delta g^{00}\sigma^2$. After performing the same rescaling, the action we examine is

$$ S \supset -\frac{1}{2} \int dt \, d^3\tilde{x} \, a^3 \left[ c_s^{-1/4} \alpha \rho^{-3/4} \dot{\tilde{\pi}}_c \tilde{\sigma}^2 + \frac{1}{2} \frac{\alpha}{f_{\pi}^2} \left( \rho^{-1/2} \dot{\tilde{\pi}}_c^2 \tilde{\sigma}^2 - c_s^{-1} \rho^{1/2} \frac{(\tilde{\partial}_i \tilde{\pi}_c)^2}{a^2} \tilde{\sigma}^2 \right) \right]. \quad (5.14) $$

We notice that the would-be Lorentz invariant interaction $(\partial_{\mu}\pi_c)^2\sigma^2$ is broken due to the rescaling of spatial (but not time) coordinates. Therefore, we expect that different strong coupling scales will be associated with $\dot{\tilde{\pi}}_c^2\tilde{\sigma}^2$ and $(\tilde{\partial}_i\tilde{\pi}_c)^2\tilde{\sigma}^2$. The interaction $\dot{\tilde{\pi}}_c^2\tilde{\sigma}^2$ can be treated as a perturbative correction to the quadratic action as long as

$$ \alpha \lesssim c_s^{1/4} \left(\frac{\rho}{H}\right)^{3/4}, \quad \text{i.e.} \quad \tilde{\alpha}^2 \lesssim \frac{c_s^{-1/2}}{4} \frac{\kappa}{\Delta_{\zeta}^2}. \qquad (5.15) $$

The smallest strong coupling scale—associated with the interaction $(\tilde{\partial}_i\tilde{\pi}_c)^2\tilde{\sigma}^2$—is

$$ \Lambda_{*} = \left(\frac{c_s}{\alpha}\right)^{2/3} \left(\frac{f_{\pi}^4}{\rho}\right)^{1/3}, \qquad (5.16) $$

which gives the following constraint to be in the perturbative regime at energies $\omega \sim H$

$$ \alpha \lesssim c_s \left(\frac{\rho}{H}\right)^{-1/2} \left(\frac{f_{\pi}}{H}\right)^2, \quad \text{i.e.} \quad \tilde{\alpha}^2 \lesssim \frac{c_s^{-1/2}}{4} \frac{\kappa}{\Delta_{\zeta}^2}. \qquad (5.17) $$

Interestingly, both conditions (5.15) and (5.17) give the same perturbativity bounds. Therefore, it is taken to be the criteria for a consistent perturbative treatment of all interactions coming from $\delta g^{00}\sigma^2$.

Next, we consider the interactions coming from $(\delta g^{00})^2\sigma$. Following the same strategy as we did previously, we find that the smallest strong coupling scale is the one associated with the interaction $\dot{\tilde{\pi}}_c^2\tilde{\sigma}$, which is

$$ \Lambda_{*} = (c_s^5 \rho^3 \Lambda_2^4)^{1/7}. \qquad (5.18) $$

<footnote>$^{36}$Note that imposing the validity of the decoupling limit—that is neglecting the mixing with gravity through terms like $\delta g^{00}\tilde{\pi}$—also gives a lower bound on $\rho/H$.</footnote>
Requiring $H \lesssim \Lambda_{\star}$ gives us the following condition

$$ \frac{H}{\Lambda_2} \lesssim c_s^{5/4} \left(\frac{\rho}{H}\right)^{3/4}, \quad \text{i.e.} \quad \frac{\tilde{\rho}}{H} \lesssim c_s^{1/2} \frac{\kappa}{\Delta_{\zeta}^2}. \qquad (5.19) $$

Lastly, the condition for perturbativity of the interaction $\sigma^3$ is found to be $\mu/H \lesssim c_s^{-3/4} (\rho/H)^{3/4}$.37

## 5.3 Radiative Corrections

Finally, we derive bounds on the size of the coupling constants based on the requirement that the theory is stable under radiative corrections, hence requiring the theory to be technically natural. These bounds are not strict, and rely on several assumptions: (i) we only require that the mass of the field $\sigma$ does not receive large loop corrections38 $\delta m^2 \lesssim m^2 \sim H^2$, (ii) we restrict ourselves to one-loop order, therefore neglecting self-interactions of $\pi_c$. Because the UV-divergent diagrams are dominated by the UV dynamics of the system $\omega > H$, the fields—even at strong mixing provided that the loops mostly receive contributions for energies $\omega > \rho$—are considered decoupled. Moreover, this enables us to use flat-space propagators already encountered in Section 4.2.1. Divergent diagrams are cut off (in energy) at the scale where the theory breaks down, namely the strong coupling scale associated with the Goldstone sector $\Lambda_{\star} \sim c_s f_{\pi}$.39

• The interaction $\sigma^3$ generates the following mass shift

$$ \delta m^2 = \mathbf{---\circ---} \sim \mu^2 \int \frac{d\omega d^3k}{(\omega^2 - k^2 - m^2)^2} \sim \mu^2, \qquad (5.20) $$

hence giving a finite loop correction. Naturalness of the mass of the field $\sigma$ therefore requires $\mu \lesssim H$.

• Let us now consider the mass correction arising from the interaction fixed by the non-linearly realised symmetry $(\partial_i \pi_c)^2 \sigma$. It reads

$$ \delta m^2 = \mathbf{---\circ---} \sim c_s^3 \frac{\rho^2}{f_{\pi}^4} \int d\omega d^3k \frac{k^4}{(\omega^2 - c_s^2 k^2)^2} \sim \rho^2, \qquad (5.21) $$

where we have set the cutoff in energy to $\Lambda_{\omega} = c_s f_{\pi}$ and the cutoff in momentum to $\Lambda_k = c_s^{-1} \Lambda_{\omega}$ to take into account the non-relativistic dispersion relation of $\pi_c$. Hence, the naturalness criterion we consider requires the quadratic mixing to be weak $\rho/H \lesssim 1$. We see that a large quadratic mixing and a small mass—and consequently the modified dispersion relation regime—is not technically natural, and can only be reached with fine-tuning.

<footnote n="37">37 We have checked that the requirement that no group velocity exceeds the speed of light gives no further constraints on the parameters.</footnote>

<footnote n="38">38 Similar bounds were derived in [40, 147]. Of course, all couplings receive corrections under renormalisation but the entire analysis is beyond the scope of this paper.</footnote>

<footnote n="39">39 In principle, the UV-divergent integrals should be cut off at the most stringent strong coupling scale, see Section 5.1. However, these scales are strongly parameter dependent, and we consider to set the UV-cutoff to $\Lambda_{\star} \sim c_s f_{\pi}$ for definiteness [147]. Note that in the case of a reduced sound speed $c_s \ll 1$, it is the most stringent bound. Let us reiterate that the naturalness conditions we derive here are contingent upon several assumptions, and therefore, they should not be regarded as absolute.</footnote>
• Next, we consider the interaction $\hat{\pi}_c^2\sigma$ which gives the following mass correction

$$ \delta m^2 = \mathbf{O}(\mathbf{O}) \sim \frac{1}{\Lambda_2^2} \int d\omega d^3k \frac{\omega^4}{(\omega^2 - c_s^2 k^2)^2} \sim c_s \frac{f_\pi^4}{\Lambda_2^2}. \qquad (5.22) $$

In the weak mixing regime, we use (4.13) to write the condition for radiative stability in term of the dimensionless power spectrum. It reads

$$ \frac{H}{\Lambda_2} \lesssim \frac{2\pi\Delta_{\zeta}}{c_s^{1/2}}, \quad \text{i.e.} \quad \frac{\tilde{\rho}}{H} \lesssim c_s^{-2}. \qquad (5.23) $$

• Finally, we consider the radiative correction induced by the interaction $\hat{\pi}_c\sigma^2$. This reads

$$ \delta m^2 = \mathbf{O}(\mathbf{O}) \sim \alpha^2 \int d\omega d^3k \frac{\omega^2}{(\omega^2 - c_s^2 k^2)(\omega^2 - k^2 - m^2)} \sim \alpha^2 \frac{f_\pi^2}{c_s}. \qquad (5.24) $$

In the weak mixing regime, the naturalness condition reads

$$ \alpha \lesssim c_s^{1/2} (2\pi\Delta_{\zeta})^{1/2}, \quad \text{i.e.} \quad \tilde{\alpha}^2 \lesssim c_s^{-1} (2\pi\Delta_{\zeta})^{-1/2}. \qquad (5.25) $$

# 6 Size of non-Gaussianities beyond Weak Mixing

The cosmological flow provides a unique opportunity to systematically study non-Gaussian signals in the strong mixing regime by fully treating the quadratic coupling. Most previous studies have examined the bispectrum in particular regimes. In this section, we numerically solve the flow equations for the bispectrum and inspect their phenomenological implications in the entire parameter space. In particular, we provide results in regions of the parameter space that are not easily accessible by analytical means.

Following standard conventions, we define the dimensionless shape function of the bispectrum S such that

$$ \langle \zeta_{\mathbf{k}_1} \zeta_{\mathbf{k}_2} \zeta_{\mathbf{k}_3} \rangle' \equiv (2\pi)^4 \frac{S(k_1, k_2, k_3)}{(k_1 k_2 k_3)^2} \Delta_{\zeta}^4. \qquad (6.1) $$

At weak quadratic mixing $\rho \lesssim m$, the $\pi_c - \sigma$ conversion process can be treated perturbatively. At leading order, the cubic interactions we consider lead to single-, double-, and triple-exchange diagrams, respectively for $(\partial_i \pi_c)^2\sigma$ and $\hat{\pi}_c^2\sigma$, $\hat{\pi}_c\sigma^2$, and $\sigma^3$. In the strong mixing regime $\rho \gtrsim m$, the quadratic coupling needs to be treated non-perturbatively. As a consequence, the notion of diagram becomes meaningless as one needs to sum up all contributions with any number of quadratic coupling insertions. In this sense, we present “non-perturbative” results. The shape function we compute should be understood as the resummed one. In the following, we quantify the size of non-Gaussianities by the usual parameter

$$ f_{\text{NL}} \equiv \frac{10}{9} S(k, k, k), \qquad (6.2) $$

where the shape function is evaluated in the equilateral configuration $k_1 = k_2 = k_3 = k$. The shapes generated here typically resemble the equilateral template, so that it makes sense to
compare $f_{\text{NL}}$ to the corresponding constraint $f_{\text{NL}}^{\text{eq}} = -26 \pm 47$ (68% CL, statistical) [156].

## 6.1 Non-local EFT

Before exposing exact numerical results and for later convenience, it is instructive to realise that one can derive an effective theory that provides a unified picture covering a large part of the phase diagram in Figure 3. In fact, as long as $m_{\text{eff}} \gg H$, the phase diagram can be described by the following non-local single-field theory [60, 64, 157, 158]

$$ \mathcal{L}/a^3 = \frac{1}{2}\dot{\pi}_c \left[ 1 + \frac{\rho^2}{-\partial_i^2/a^2 + m^2} \right] \dot{\pi}_c - \frac{1}{2} \frac{(\partial_i \pi_c)^2}{a^2} - \mu \rho^3 \left[ (-\partial_i^2/a^2 + m^2)^{-1} \dot{\pi}_c \right]^3 \\ - \frac{1}{2} \alpha \rho^2 \dot{\pi}_c \left[ (-\partial_i^2/a^2 + m^2)^{-1} \dot{\pi}_c \right]^2 - \frac{\rho}{2} \left[ \frac{1}{\Lambda_1} \frac{(\partial_i \pi_c)^2}{a^2} + \frac{1}{\Lambda_2} \dot{\pi}_c^2 \right] (-\partial_i^2/a^2 + m^2)^{-1} \dot{\pi}_c, \qquad (6.3) $$

where we have omitted self-interactions of $\pi_c$ parameterised by $\lambda_{1,2}$ in (4.6). We have also fixed the intrinsic speed of sound of $\pi_c$ to unity for this Section; $c_s = 1$. This Lagrangian results from integrating out the massive field $\sigma = (-\Box + m^2)^{-1} \rho \dot{\pi}_c$ and neglecting the time derivatives in the d’Alembert operator (see Eq. (B.1)). As such, Eq. (6.3) should be understood as the leading order contribution in a time derivative expansion, higher-order terms being encoded in an infinite number of operators coming from

$$ \frac{1}{-\Box + m^2} = \frac{1}{-\partial_i^2/a^2 + m^2} \sum_{n=0}^{\infty} (-1)^n \left[ (\partial_t^2 + 3H\partial_t)(-\partial_i^2/a^2 + m^2)^{-1} \right]^n . \qquad (6.4) $$

Let us now take the flat-space limit to develop our intuition. One immediately sees that the above expansion is valid whenever the low-energy modes obey a non-relativistic dispersion relation such that $\omega^2 \ll k^2 + m^2$. The quadratic part of the Lagrangian (6.3) gives the following dispersion relation

$$ \omega^2 = \frac{k^2}{1 + \frac{\rho^2}{k^2 + m^2}}. \qquad (6.5) $$

In the limit $k \ll m$, this boils down to $\omega = \tilde{c}_s k$ with $\tilde{c}_s$ defined in (4.30). At strong mixing such that $\rho \gg k \gg m$, one recovers $\omega = k^2/\rho$. Therefore, the non-local effective theory interpolates between the reduced speed of sound and the modified dispersion relation regimes. However, being intrinsically single-field, this theory misses the particle production and therefore does not describe the cosmological collider signal in the squeezed limit of the three-point correlators (see Section 7). We will see later that (6.3) is very convenient to estimate the amplitude of the bispectrum using dimensional analysis.

## 6.2 Phenomenology

*Conversion from $\pi$ to $\zeta$.* Independently of cubic interactions in the Lagrangian (4.6), the non-linear relation between $\pi$ and $\zeta$ generates a contribution to the bispectrum, see the last two terms in Eq. (4.12) and the discussion below. It can be computed solely from the linear dynamics, and it is exponentially suppressed in all parameter space of the $m-\rho$ phase diagram, except in the
<figure><img src="image_6.png" /><figcaption>Figure 5: Shape function phase diagram for the $(\partial_i \pi_c)^2 \sigma$ (top left), $\dot{\pi}_c^2 \sigma$ (top right), $\dot{\pi}_c \sigma^2$ (bottom left), and $\sigma^3$ (bottom right) interactions, in the equilateral configuration $k_1 = k_2 = k_3 = k$. The dashed line corresponds to the $M=H$ limit, isolating the modified dispersion relation regime $M/H \ll 1$ and the reduced speed of sound regime $M/H \gg 1$. The dotted thin line corresponds to $S(k, k, k) = 0$ where the shape function changes sign. We have set $N = 60$ e-folds of super-horizon evolution. Slices of the phase diagram for various interactions can be found in Figure 6.</figcaption></figure>

effectively light field regime $m_{\text{eff}} \ll H$. Explicit results, see Figure 14 in Appendix B, show that this conversion contribution to $f_{\text{NL}}$ is at most of order $10^{-2}$ in the parameter space we display, for $N = 60$ e-folds of super-horizon evolution. Such a small value is negligible for today's cosmology, but should be taken into account in general.

In the following, we concentrate on the “intrinsic” bispectrum generated by cubic interactions in the Lagrangian (4.6), and which can give sizeable contributions. Figure 5 shows the $f_{\text{NL}}$ parameter in the $(m/H, \rho/H)$ plane for the various interactions. For easy comparison and better visualisation, we also present slices of the phase diagrams in Figure 6. For each interaction, we consider the three regions of the phase diagram in turn.

($\partial_i \pi_c)^2 \sigma$ interaction. First, we consider the bispectrum generated by the interaction $(\partial_i \pi_c)^2 \sigma$. In the weak mixing regime $\rho \lesssim m$, the size of the single-exchange diagram can be estimated (see
insert below), yielding$^{40}$

$$f_{\text{NL}}^{(\partial_i \pi_c)^2 \sigma} \sim \left(\frac{\rho}{H}\right)^2. \qquad (6.6)$$

It is found that this scaling indeed extends up to values of the quadratic coupling $\rho \sim m$ as we increase the mass. Hence the leading perturbative result still gives a correct estimation for $\rho/H \gtrsim 1$ when the field $\sigma$ is heavy enough. The associated non-Gaussian signal cannot reach order unity.

As function of the bare mass $m/H$, the bispectrum signal presents a suppression whose leading contribution scales as $(H/m)^2$ for large mass.$^{41}$ This comes from the leading imprint of higher-dimension operators when the heavy field is integrated out, as can be explicitly seen in (6.3) in the limit $\partial_i^2/a^2 \ll m^2$. The transition to this EFT regime is found to occur at $m/H \gtrsim 1$, as can be observed in the upper panel of Figure 6. The hilltop visible around $m/H \sim 1$ is due to the fact that the intrinsic bispectrum decreases slower than the power spectrum squared as function of $m/H$, resulting in a boosted shape function.

In the strong mixing regime $\rho \gtrsim m$, the relevant mass scale to describe the dynamics is $M = m^2/\rho$ as explained in Section 4.2. Therefore, the transition to the EFT regime is shifted to higher values of the bare mass $m/H \gtrsim \sqrt{\rho/H}$. In the modified dispersion relation regime, one can use (6.3) and estimate the size of the non-Gaussianity, obtaining

$$f_{\text{NL}}^{(\partial_i \pi_c)^2 \sigma} \sim \frac{\rho}{H}. \qquad (6.7)$$

We show in the insert below how this scaling is found using dimensional analysis. The found scaling is confirmed numerically. For the same reason as explained previously, this regime is reached for higher values of the quadratic mixing as the heavy field mass increases. Importantly, note that the size of this non-Gaussian signal monotonically increases as a function of the quadratic mixing, and can be large.

*Estimating the size of the bispectrum.*—In this insert, we show how one can easily recover the scalings (6.6) and (6.7) using purely dimensional analysis arguments and the non-local effective Lagrangian (6.3). This theory being single-field, we can estimate the size of the bispectrum using

$$f_{\text{NL}} \sim \frac{\mathcal{L}_3}{\mathcal{L}_2} \Delta_{\zeta}^{-1}, \qquad (6.8)$$

evaluated at the energy scale that we probe during inflation i.e. $\omega \sim H$. In the weak mixing regime, using the linear dispersion relation (we recall that we have fixed $c_s = 1$ for simplicity), we obtain

$$f_{\text{NL}}^{(\partial_i \pi_c)^2 \sigma} \sim \left(\frac{\rho}{H}\right)^2, \qquad (6.9)$$

where we have evaluated the gradient term in $(-\partial_i^2/a^2 + m^2)^{-1}$ at the Goldstone boson horizon crossing $|k\tau| \sim 1$ and assumed $m \sim H$. In the strong mixing regime, using $\dot{\pi}_c \sim H\pi_c$ and the modified

<footnote>$^{40}$Note that the non-linearly realised symmetry fixes this cubic interaction coupling constant to $f_\pi/\Lambda_1 = \rho/f_\pi$.</footnote>

<footnote>$^{41}$The signal also contains other terms suppressed by the exponential Boltzmann factor $e^{-\pi\mu}$ where $\mu^2 = (m/H)^2 - 9/4$. However for large mass, the $(H/m)^2$ suppression dominates in the equilateral configuration.</footnote>
dispersion relation $\omega = k^2/\rho$, one obtains

$$f_{\text{NL}}^{(\partial_i \pi_c)^2 \sigma} \sim \frac{\rho}{H} \left(\frac{H}{f_{\pi}}\right)^2 \frac{\pi_c}{H} \Delta_{\zeta}^{-1}. \qquad (6.10)$$

We finally use $\Delta_{\zeta} \sim (H/f_{\pi})^2(\rho/H)^{1/4}$ (see e.g. (4.29)) to obtain $\pi_c/H \sim (\rho/H)^{1/4}$ which, in the end, gives the desired scaling. A similar analysis can be performed for the other interactions to find the scalings (6.13), (6.15) and (6.17). In these cases, the “propagator” $(-\partial_i^2/a^2 + m^2)^{-1}$ can be estimated by evaluating the gradient term at the Goldstone boson $\rho$-horizon $|k\tau| \sim \sqrt{\rho/H}$, leading to

$$(k^2/a^2 + m^2)/H^2 \sim \rho/H. \qquad (6.11)$$

When the field is effectively light $m/H \ll 1$ and $\rho/H < 1$, we see that the $f_{\text{NL}}$ parameter does not present an amplification due to super-horizon sourcing, visible in the phase diagram of the power spectrum in Figure 4. We recall that this enhancement depends on the number of super-horizon e-folds until the end of inflation. In fact, the bispectrum is indeed enhanced but this boost is balanced by the amplification of the power spectrum to give no visible effect in the $f_{\text{NL}}$ parameter.

**$\dot{\pi}_c^2 \sigma$ interaction.** At weak mixing $\rho \lesssim m$, the size of the non-Gaussianity arising from the interaction $\dot{\pi}_c^2 \sigma$ is$^{42}$

$$f_{\text{NL}}^{\dot{\pi}_c^2 \sigma} \sim \frac{\rho}{H} \times \frac{\tilde{\rho}}{H}, \qquad (6.12)$$

which we confirm numerically up to $\rho \sim m$, similarly to the previous interaction. At strong mixing $\rho \gtrsim m$, the EFT regime is reached for higher values of the bare mass $m/H \gtrsim \sqrt{\rho/H}$ for the same reason as the interaction $(\partial_i \pi_c)^2 \sigma$. Accounting for perturbativity bounds on the couplings, we can accommodate for large non-Gaussianities in this regime. In the modified dispersion relation regime, the size of the non-Gaussianity is

$$f_{\text{NL}}^{\dot{\pi}_c^2 \sigma} \sim \frac{\tilde{\rho}}{H} \left(\frac{\rho}{H}\right)^{-1}, \qquad (6.13)$$

which, for $\rho/H \gg 1$, gives a lower signal compared to that in the weakly coupled regime. Therefore, for this interaction, a large quadratic mixing does not necessarily imply large non-Gaussianities. In the lower panel of Figure 6, one can appreciate that this interaction is damped at large quadratic mixing, as opposed to the spatial derivative interaction. In fact, from the non-local Lagrangian (6.3) and the modified dispersion relation $\omega = k^2/\rho$, it is easy to realise that interactions with strictly negative powers of spatial gradients—including the “propagators” $(k^2/a^2 + m^2)^{-1}$—are suppressed at strong quadratic mixing.

In the effectively light field case and as can be better seen in Figure 6, the shape function reaches a plateau as the boost in the bispectrum is compensated by that of the power spectrum. The transition to this regime, either from the $\rho$ or the $m$ directions, is delimited by the $f_{\text{NL}}$

<footnote>$^{42}$Due to spontaneously broken de Sitter boosts, the covariance of the cubic interaction $(\partial_{\mu}\pi_c)^2 \sigma$ is violated. The non-linearly realised symmetry only fixes the spatial derivative interaction, see Eq. (4.5).</footnote>
parameter changing sign, as represented by the dotted line in Figure 5.

Recall that because we consider interactions that break de Sitter boosts, we have split $(\partial_{\mu}\pi_{c})^{2}\sigma$ into $(\partial_{i}\pi_{c})^{2}\sigma$ and $\dot{\pi}_{c}^{2}\sigma$, the coefficient of the first one being fixed in terms of $\rho/H$. We have shown that the time derivative interaction leads to a completely different phenomenology than the spatial derivative interaction, a fact that was not enough appreciated before.

<figure><img src="image_7.png" /><figcaption>Figure 6: Slices of the shape function phase diagrams as function of the bare mass $m/H$ for $\rho/H = 0.1$ (top left) and $\rho/H = 2$ (top right), and as function of the quadratic mixing $\rho/H$ for $m/H = 1$ (bottom left) and $m/H = 2$ (bottom right), in the equilateral configuration $k_{1} = k_{2} = k_{3} = k$. For better visualisation, the shape generated by the interaction $\sigma^{3}$ has been multiplied by $\Delta_{\zeta}$.</figcaption></figure>

### $\dot{\pi}_{c}\sigma^{2}$ interaction.
At weak mixing $\rho \lesssim m$, the non-Gaussian signal arising from the interaction $\dot{\pi}_{c}\sigma^{2}$ scales as

$$f_{\text{NL}}^{\dot{\pi}_c \sigma^2} \sim \tilde{\alpha}^2 \left(\frac{\rho}{H}\right)^2. \qquad (6.14)$$

We observe that this scaling is valid up to $\rho \sim m$. In the modified dispersion relation regime, the amplitudes scales as

$$f_{\text{NL}}^{\dot{\pi}_c \sigma^2} \sim \tilde{\alpha}^2 \left(\frac{\rho}{H}\right)^{-1}, \qquad (6.15)$$

which again gives a comparable amplitude compared to that in the weakly coupled regime. As shown in Figure 6, the cosmological flow approach enables us to effortlessly probe the intermediate
regime where analytical computations fail. A large non-Gaussian signal in the strong mixing regime is possible while remaining under perturbative control.

In the light field regime, one can appreciate once more that the shape function is almost constant, for the same reason as for the previous interaction.

*$\sigma^3$ interaction.* At leading order in the weak mixing regime, the interaction $\sigma^3$ leads to the triple-exchange diagram whose amplitude is

$$f_{\text{NL}}^{\sigma^3} \sim \frac{\mu}{H} \left(\frac{\rho}{H}\right)^3 \Delta_{\zeta}^{-1}. \qquad (6.16)$$

It is well-known that this amplitude of non-Gaussianity is boosted due to the factor $\Delta_{\zeta}^{-1} \sim 10^5$, which leads to the largest non-Gaussian signal. This is explained by the fact that, for $\mu \sim H$ and $\rho \sim H$, the entire non-Gaussian signal in the $\sigma$-sector should be transferred to the $\zeta$-sector. The actual size of non-Gaussianities in the bispectrum being $f_{\text{NL}}\Delta_{\zeta}$, one should have $f_{\text{NL}}\Delta_{\zeta} \sim (\mu/H)(\rho/H)^3$. As function of $m/H$, the intrinsic bispectrum scales as $(H/m)^6$ [159]—it drops faster than $(H/m)^2$, leading to no characteristic enhancement around $m/H \sim 1$ in the shape function, as one can observe in the upper panel of Figure 6. When the bare mass increases, the modified dispersion relation regime is found for higher values of the quadratic mixing. We indeed confirm numerically that the scaling (6.16) extends up to $\rho/H \sim 1$ when increasing the bare mass. For $M/H \ll 1$ and $\rho/H \gg 1$, the non-Gaussian signal scales as

$$f_{\text{NL}}^{\sigma^3} \sim \frac{\mu}{H} \left(\frac{\rho}{H}\right)^{-3/4} \Delta_{\zeta}^{-1}, \qquad (6.17)$$

which—although a priori suppressed by $\rho/H$—gives a larger contribution than in the weak mixing regime for typical values of the quadratic coupling. Here, as well as for other interactions, the change of sign of the shape function for $m/H \sim 2$ can be attributed to a cancellation between two contributions resembling equilateral shapes and giving rise to an orthogonal shape.

In the effectively light field regime the $f_{\text{NL}}$ parameter is decreasing as the bispectrum is not enhanced like the power spectrum.

**Summary.** Accounting for perturbativity bounds derived in Section 5 and the exact size of non-Gaussianities presented in Figure 5, one can obtain sizable non-Gaussianities—typically $f_{\text{NL}} \gg 1$—in the strong mixing regime for all interactions. At weak mixing, all interactions except $(\partial_i\pi_c)^2\sigma$ can lead to $f_{\text{NL}} \sim \mathcal{O}(1)$ signals. Typically, the signal is largest for a mass of $\sigma$ of order $H$ or less. We have checked that, for each of these four cubic interactions, both in the weak and the strong mixing regimes, saturating the perturbativity bounds universally gives $f_{\text{NL}}\Delta_{\zeta} \lesssim 1$, as expected.

**Self-regularised bispectrum.** It is worth insisting on the convenience of the cosmological flow approach when computing the bispectrum. When using analytical mode functions in calculations treating the quadratic mixing as a perturbative interaction—therefore necessarily narrowing down the range of applications to the weak mixing regime—the integrals appearing in the in-in formula are very difficult to compute for at least two reasons: (i) they involve several nested integrals,
and (ii) they usually do not converge in the IR or UV, and need to be regularised.

Traditionally, the convergence in the UV is achieved by considering the *iϵ*-prescription—slightly tilting the integration contour in the imaginary direction in the complex plane—that turns highly oscillatory mode functions into damped exponentials.<fnref n="43" /> This regularisation requires the use of the factorised form of the in-in formula. In the IR, the factorised form of the in-in formula leads to individual integrals each containing spurious divergences. The convergence is instead more easily achieved by using the commutator form of the in-in formula. As a consequence, one needs to split the in-in formula into two contributions, one converging explicitly in the UV and the other one in the IR. This mixed form requires introducing a fiducial intermediate time as a regulator (see [35] for more details).

Going beyond a perturbative treatment of the mixing by numerically solving the matrix mode functions of the full quadratic theory removes the need to perform nested integrals to compute the bispectrum. As for the issue of the UV convergence, it can then be addressed by Wick-rotating the theory, i.e. solving the mode functions in Euclidean time before evaluating the corresponding time integral defining the bispectrum [43, 131]. However, this method has the drawback that if one wants to determine the time-evolution of correlators during inflation, one needs to evaluate different integrals for each time considered, whereas this is automatically obtained from solving in one go the set of cosmological flow equations in our approach. Furthermore, going beyond the bispectrum, nested integrals describing exchange processes necessarily reappear.

The cosmological flow completely bypasses all these complications because: (i) this method does not require the use of neither analytical nor numerical mode functions, (ii) it relies on evolving the correlators on the real time axis, hence there is no need to Wick rotate the integrand, and (iii) it systematically works for all tree-level *n*-point functions.

# 7 Cosmological Collider Phenomenology

The squeezed limit of the bispectrum is believed to be a robust probe of the particle content during inflation. One major advantage of the cosmological flow is that it provides a generic framework to study the squeezed limit of the bispectrum in a systematic and exact manner. In this section, we present several cosmological collider signals, including the strong mixing regime for which it is important to take into account the notion of *effective mass* $m_{\text{eff}}^2 = m^2 + \rho^2$ instead of the bare mass. Moreover, we show that a phase shift and a non-trivial amplitude of the oscillations can be seen as a signature of primordial fluctuations undergoing strong mixing. We also give novel insight on the role of the sound speeds on cosmological collider signals. We finish by exposing the cosmological collider flow—i.e. the squeezed limit of the bispectrum as function of time—that enables us to shed light on the characteristic time scales of the dynamics.

In this section, for definiteness, we will focus on the bispectrum arising from the $(\partial_t \pi_c)^2 \sigma$ and $\dot{\pi}_c \sigma^2$ interactions (except in 7.2), and we will be interested in the case where the field $\sigma$ is heavy; $\mu_{\text{eff}} > 0$.

<footnote n="43">⁴³Other methods can be used to speed up the numerical convergence in the UV, see e.g. [125, 126, 134, 160]. Note also that such a difficulty is bypassed as well in the modal approach [129, 133].</footnote>
## 7.1 Strongly Mixed Colliders

We first look at how the bispectrum shape in the squeezed limit depends on the quadratic mixing. In this limit, the bispectrum is sensitive to the heavy field mass, as it governs the frequency of the oscillations. Figure 7 shows the bispectrum shapes as functions of the momentum ratio in the isosceles-triangle configuration, for fixed $\mu_{\text{eff}}$ but varying the quadratic mixing (hence adjusting the bare mass). We observe that the cosmological collider signals oscillate at the same frequency. Indeed as we have seen in Section 4.2, the relevant notion of mass for this system is the effective mass $m_{\text{eff}}$ and not the bare mass $m$. As a consequence, a light but strongly mixed field and a heavy but weakly mixed field cannot be distinguished based only on the signal frequency. For illustration purposes and to focus on the oscillation frequencies only, the signals have been normalised to unity in the equilateral configuration (see Section 6). A notable feature is the phase shift of the signal to smaller momentum ratios for larger quadratic mixings. This constitutes a striking example that having complete predictions for the signal—amplitude, frequency, phase, as

<figure><img src="image_8.png" /><figcaption>Figure 7: Shape functions for the $(\partial_i \pi_c)^2 \sigma$ interaction (left panel) and $\pi_c \sigma^2$ interaction (right panel) in the isosceles-triangle configuration $k_1 = k_2$, for $\mu_{\text{eff}} = 2$ (top) and $\mu_{\text{eff}} = 5$ (bottom) varying the quadratic mixing $\rho/H$. For illustration purposes, we have normalised the amplitudes of the signals to be the same in the equilateral configuration. For $\mu_{\text{eff}} = 5$, we also have set $c_s = 0.1$ to amplify the oscillations of the cosmological collider signal (signals for different values for the sound speed—including $c_s = 1$—are presented in figure 9).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Shape functions for the $(\partial_i \pi_c)^2 \sigma$ interaction (left panel), and $\pi_c \sigma^2$ interaction (right panel) in the isosceles-triangle configuration $k_1 = k_2$, for a fixed bare mass $m/H = 2$ and varying the quadratic mixing $\rho/H$.</figcaption></figure>

well as contaminations from equilateral shapes—is important to disentangle the bare mass from the effective one, and will be needed to correctly infer the model parameters from the cosmological data. As discussed in Section 4.2, this phase shift is due to the resonance of the heavy field before it starts oscillating.

Because of the arbitrary normalisation of the shape function, Figure 7 does not exhibit the size of non-Gaussianities nor the amplitude of the cosmological collider signals. To shed light on these two aspects, we show in Figure 8 the shape function in the isosceles-triangle configuration for different values of the quadratic mixing $\rho/H$ while fixing the bare mass $m/H$. In the equilateral configuration $k_3/k_1 = 1$, the parametric dependence of $f_{\text{NL}}$ on $\rho/H$ has been discussed in Section 6. Yet, a noticeable feature of the cosmological collider signals is the strong dependence of its amplitude on the quadratic coupling $\rho/H$ that is not monotonic. In the weak mixing regime $\rho \lesssim m$, the amplitude grows with $\rho$, after which it decreases in the strong mixing regime. At strong mixing for the $(\partial_i \pi_c)^2 \sigma$ interaction, this indicates that although the size of the total signal is amplified, signatures of the cosmological collider signal are harder to detect, as it is exponentially suppressed by the effective mass $e^{-\pi \mu_{\text{eff}}}$.

To sum up, the cosmological collider (denoted CC in the following equation) signal in the squeezed limit can be written as

$$ \lim_{k_3/k_1 \to 0} S^{\text{CC}}(k_1, k_1, k_3) \sim \left(\frac{k_3}{k_1}\right)^{1/2} \mathcal{A}\left(\frac{m}{H}, \frac{\rho}{H}\right) \cos\left[\mu_{\text{eff}} \log\left(\frac{k_3}{k_1}\right) + \delta\left(\frac{m}{H}, \frac{\rho}{H}\right)\right]. \quad (7.1) $$

It qualitatively describes the squeezed limit even beyond the analytical lamppost of the weakly mixed regime, with an amplitude and phase information that can be efficiently obtained from the cosmological flow approach. The $m/H$-dependence of the amplitude $\mathcal{A}$ and the phase $\delta$ has been computed in [37, 40, 55, 59] in the weakly mixed regime and only for the single-exchange diagram. It should be noted that for realistically squeezed configurations like the ones presented in the figures, the shape function does not oscillate around zero. This means that it always contains a contamination part—analytic in momenta—which can be described by the (possibly non-local)
single-field EFT.

Eventually, let us highlight that if we ever detect the cosmological collider signal, it is more likely to be attributed to a substantial mixing $\rho \sim m$ for which the non-analytic signal is the largest. This illustrates one of the main aspects motivating the development of the cosmological flow approach: interesting physics lying beyond the lamppost of analytical calculations can be missed, therefore biasing the interpretation of data.

## 7.2 Sound-Speed Colliders

Next, we present the effect of varying the speed of sound $c_s$ on the cosmological collider signal. More precisely, throughout the paper, we have used $c_s$ to denote the speed of sound of the Goldstone boson $\pi_c$. However, in this section exclusively, we allow the field $\sigma$ to also have a non-unity speed of sound $c_\sigma$, and $c_s$ will actually denote the ratio of the sound speed of $\pi_c$ and the massive field: $c_s \equiv c_\pi/c_\sigma$. Therefore, $c_s < 1$ means that the Goldstone boson propagates slower than the field $\sigma$, and the opposite for $c_s > 1$, which by no means implies superluminal propagation. Crucially, as a consequence of the non-linearly realised symmetry, allowing $c_\sigma \neq 1$ induces additional mixings between $\pi_c$ and $\sigma$ that one is forced to consider.

Indeed, a non trivial propagation speed for $\sigma$ can only be generated by the unitary-gauge operator $(g^{0\mu}\partial_\mu\sigma)^2$, invariant under time-dependent spatial diffeomorphism [140]. Therefore, as long as $c_\sigma \neq 1$, one should consider it in the theory, i.e.

$$S \supset - \int d^4 x \sqrt{-g} \tilde{e} (g^{0\mu} \partial_{\mu} \sigma)^2, \qquad (7.2)$$

where $\tilde{e}$ is a dimensionless coupling that we assume to be constant. Reintroducing the Goldstone boson in the decoupling limit via $g^{0\mu}\partial_\mu\sigma \to -(1+\tilde{\pi})\dot{\sigma} + \frac{1}{a^2}\partial_i\pi\partial_i\sigma$, and discarding quartic terms in the fields—not relevant for the bispectrum—, one obtains$^{44}$

$$\mathcal{L}/a^3 \supset -\tilde{e} \left[ \dot{\sigma}^2 + 2\tilde{\pi}\dot{\sigma}^2 - 2\dot{\sigma} \frac{(\partial_i \pi)(\partial_i \sigma)}{a^2} + \dots \right]. \qquad (7.3)$$

Because the second term is also generated by the operator $\delta g^{00}(g^{0\mu}\partial_\mu\sigma)^2$, its coupling constant is not fully fixed by $\tilde{e}$. We therefore do not consider it in the following. By contrast, the interaction $\dot{\sigma}(\partial_i\pi)(\partial_i\sigma)$ can only be generated by the operator (7.2). Hence, in the same way as a non-trivial sound speed for $\pi$ universally comes with the interaction $\tilde{\pi}(\partial_i\pi)^2$ whose size is controlled by $c_\pi$, a propagation speed for $\sigma$ different from unity universally comes with the interaction in $\dot{\sigma}(\partial_i\pi)(\partial_i\sigma)$ whose size is fixed by $c_\sigma$. The only difference is that the observable non-Gaussian signal it generates also depends on the $\sigma$ to $\pi$ transfer, i.e. on the quadratic mixing between the

<footnote>$^{44}$Under a time diffeomorphism $t \to \tilde{t} = t + \pi(t, \mathbf{x})$, the four-dimensional partial derivative (inversely) transforms as $\tilde{\partial}_\mu \to \partial_\mu = \partial_\mu \tilde{x}^\nu \tilde{\partial}_\nu = \partial_\mu (t+\pi) \tilde{\partial}_0 + \partial_\mu \tilde{x}^i \tilde{\partial}_i$. The temporal component therefore transforms as $\tilde{\partial}_0 \to (1+\tilde{\pi})\tilde{\partial}_0$ which implies $\partial_0 \to \tilde{\partial}_0 = \frac{\partial_0}{1+\tilde{\pi}}$. The spatial ones transform as $\tilde{\partial}_i \to \partial_i = \tilde{\partial}_i + (\partial_i\pi)\tilde{\partial}_0$ which implies $\partial_i \to \tilde{\partial}_i = \partial_i - \frac{\partial_i\pi}{1+\tilde{\pi}}\partial_0$. The off-diagonal metric components transform as $g^{0i}(t, \mathbf{x}) \to \tilde{g}^{0i}(\tilde{t}, \tilde{\mathbf{x}}) = \partial_\mu(t+\pi)\partial_\nu\tilde{x}^i g^{\mu\nu}(t, \mathbf{x}) = \frac{1}{a^2}\partial_i\pi$ where we have ultimately evaluated the metric on the background. Collecting everything together, and noting that $g^{00}(t, \mathbf{x}) \to \tilde{g}^{00}(\tilde{t}, \tilde{\mathbf{x}}) = -(1+\tilde{\pi})^2 + \frac{1}{a^2}(\partial_i\pi)^2$ (see Section 4.1), we obtain the following exact transformation $g^{0\mu}\partial_\mu\sigma \to -(1+\tilde{\pi})\dot{\sigma} + \frac{1}{a^2}(\partial_i\pi)(\partial_i\sigma)$, valid to all orders in the fields.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 9: Shape functions for the $(\partial_i \pi_c)^2 \sigma$ interaction (left panel) and $\dot{\sigma}(\partial_i \pi_c)(\partial_i \sigma)$ interaction (right panel) in the isosceles-triangle configuration $k_1 = k_2$, for $c_\sigma = 0.1$, $\mu_{\text{eff}} = 4$ and $\rho/H = 4$, and with $c_s = 0.1$ (top), $c_s = 1$ (middle), and $c_s = 10$ (bottom). We recall that $c_s = c_\pi/c_\sigma$ and that both interactions are fully fixed by the non-linearly realised symmetry (either by $\rho$ or $c_\sigma$).</figcaption></figure>

two fields.

Explicitly, when adding the effects of (7.2) to the action considered in Section 4.1, one obtains
the following Lagrangian

$$ \mathcal{L}/a^3 = \frac{1}{2} \left[ \dot{\pi}_c^2 - c_s^2 \frac{(\partial_i \pi_c)^2}{a^2} \right] + \frac{1}{2} \left[ \dot{\sigma}_c^2 - \frac{(\partial_i \sigma_c)^2}{a^2} - m^2 \sigma_c^2 \right] + \rho \dot{\pi}_c \sigma_c \\ - \frac{c_s^{3/2}}{2f_{\pi}^2} \rho \frac{(\partial_i \pi_c)^2}{a^2} - \frac{c_s^{3/2}}{f_{\pi}^2} \left( \frac{1}{c_{\sigma}^2} - 1 \right) \dot{\sigma}_c \frac{\partial_i \pi_c \partial_i \sigma_c}{a^2}, \qquad (7.4) $$

where we concentrated on the new effects coming from $c_{\sigma} \neq 1$, with

$$ c_{\sigma}^{-2} = 1 - 2\tilde{e}, \quad c_s = c_{\pi}/c_{\sigma}, \quad f_{\pi}^4 = 2M_{\text{pl}}^2 |\dot{H}| c_{\pi}, \quad \rho = 2c_{\pi}^{3/2} c_{\sigma} \frac{\tilde{M}_1^3}{f_{\pi}^2}, \qquad (7.5) \\ \pi_c = c_s^{-3/2} f_{\pi}^2 \pi, \quad \sigma_c = c_{\sigma}^{1/2} \sigma, \quad m \to c_{\sigma} m. $$

Here, we rescaled spatial coordinates $\tilde{x} = x/c_{\sigma}$. This does not affect the physics in our scale-invariant theories, hence we omit the tilde for simplicity. In this way, the quadratic Lagrangian in terms of canonically normalised fields takes the same form as the one studied in the rest of the paper, simply allowing for $c_s > 1$, and where the mass $m$ in Eq. (7.4) has been rescaled by a factor $c_{\sigma}$ compared to the bare mass in (4.3). We leave for future work the study of the full phase diagram and non-Gaussianities, varying $m, \rho$ and $c_s$. Nonetheless, an important aspect made transparent in the Lagrangian (7.4) is worth highlighting: a small sound speed $c_{\sigma}$ for the additional field boosts the overall amplitude of non-Gaussianities by a factor $\sim 1/c_{\sigma}^2$. Indeed, at fixed $m, \rho$ and $c_s$, the sound speed $c_{\sigma}$ only appears in the overall multiplicative factors, $1/c_{\sigma}^2$ and $1/c_{\sigma}^2 - 1$ respectively, setting the sizes of the cubic interactions in $(\partial_i \pi_c)^2 \sigma$ and $\dot{\sigma}(\partial_i \pi_c)(\partial_i \sigma)$. The corresponding rather large non-Gaussian signals are well visible in Figure 9, where we show the shape functions as a function of the squeezing parameter $k_3/k_1$ at strong mixing $\rho/H = 4$, for $c_{\sigma} = 0.1$ and for $c_s \in \{0.1, 1, 10\}$.

For the $(\partial_i \pi_c)^2 \sigma$ interaction on the left panel, the amplitude of the non-analytic part of the signal, i.e. the cosmological collider signal, increases when $c_s$ decreases. In the weak mixing regime, the scaling of the amplitude with $c_s$ is known analytically. For $c_s = 1$, it is known that the amplitude of the non-analytic part of the signal is suppressed by $e^{-\pi \mu_{\text{eff}}}$.<fnref n="45" /> However, for $c_s \ll 1$, this suppression takes the form $e^{-(\pi/2+c_s)\mu_{\text{eff}}}$ [60], hence leading to an amplification of the cosmological collider signal. We find indeed that this tendency extends to the strong mixing regime, and that for $c_s > 1$ the signal is further damped. On the contrary, the amplitude of the cosmological collider signal arising from the $\dot{\sigma}(\partial_i \pi_c)(\partial_i \sigma)$ interaction is mildly dependent on $c_s$.

For $c_s > 1$ and in the weak mixing regime, the cosmological collider signal is shifted towards larger momentum ratios, closer to the equilateral configuration [61], as shown explicitly for single-exchange diagrams (see [60] for an explanation of this phenomenon related to the various time scales involved). This phenomenon is not apparent at strong mixing for our parameter choices. However, more importantly, in this regime, the total observable signal is dominated by the (would-be) double-exchange diagram, coming from the interaction unavoidably generated by the non-

<footnote n="45">⁴⁵These scalings were derived in the weak mixing regime, and for single-exchange diagrams only. In Section 4.2, we have given evidence that these scalings extend to the strong mixing regime albeit with the replacement $\mu \to \mu_{\text{eff}}$. It should also be mentioned that this Boltzmann suppression is only the leading behaviour of a more complicated scaling of the amplitude.</footnote>
linearly realised symmetry in the presence of a non-trivial sound speed for the additional field, a
point that was not appreciated before.

7.3 Cosmological Collider Flow

Following the time evolution of the three-point correlators for a fixed momentum configuration
allows us to reconstruct the flow of the cosmological collider signal. This is shown in Figure 10
where we depict the shape function in the isosceles-triangle configuration as function of time.
The measure of time is taken to be the number of e-folds elapsed since sound horizon crossing
of the long mode: N − N<sub>ℓ</sub>, with N<sub>ℓ</sub> such that k<sub>ℓ</sub> = k<sub>3</sub> = a(N<sub>ℓ</sub>)H(N<sub>ℓ</sub>). One can notice that the
cosmological collider signal builds up progressively as function of time, starting with configura-
tions close to the equilateral one. This representation makes it clear that the squeezed limit of
the three-point correlator is a probe of the super-horizon evolution of the additional massive field
σ between N<sub>ℓ</sub> and N. This example showcases the utility of the cosmological flow approach as
it enables us to reveal the physics of various processes by analysing the characteristic time scales
at play, allowing for a direct access to the bulk time evolution of correlators. A video of the
cosmological collider flow can be found here.<sup>46</sup> It is instructing to notice the different behaviours
between the weakly and strongly mixed cases. At weak mixing, the system is underdamped and
the signal exhibits temporal oscillations about its late-time asymptotic form before converging.
At strong mixing, the two fields are so bound that the late-time signal is asymptotically reached
with overdamped oscillations.

<figure><img src="image_11.png" /><figcaption>Figure 10: Time evolution of the shape function for the $\dot{\pi}_c\sigma^2$ interaction in the isosceles-triangle configuration $k_1 = k_2$, for $\mu_{\text{eff}} = 5$, $\rho/H = 0.1$ and $c_s = 0.1$. The measure of time $N - N_\ell$ is defined to be the time in e-folds elapsed since sound-horizon crossing of the long mode $k_\ell = k_3$. We have normalised the amplitude of the signal to unity in the equilateral configuration for illustration purposes. The dotted black line corresponds to the shape function at $N - N_\ell = 20$, a time by which the signal has become stationary.</figcaption></figure>

<footnote n="46"><sup>46</sup>github.com/deniswerth/Cosmological-Collider-Flow</footnote>
# 8 Primordial Features

We have illustrated that the cosmological flow is a generic framework that enables one to study non-Gaussian phenomenology in scale-invariant cases that are not easily accessible by analytical means. However, realistic models of inflation generically lead to features in the scale dependence of primordial correlation functions. In this section, we will show that features can be accounted for systematically, both in the scale-dependence and the shape-dependence of correlators. As a proof of concept to illustrate the utility of our approach, we will first perform an exact calculation of the power spectrum and the bispectrum scale-dependence in the case of a time-varying quadratic mixing, and compare it with the predictions given by an effective single-field theory with a time-varying speed of sound. We will then show how features are encoded in the bispectrum squeezed limit, giving a new way of probing the inflationary landscape through soft limits of cosmological correlators, and therefore extending the reach of cosmological collider phenomenology.

## 8.1 Effective Description of Features

Let us first discuss how time-varying background quantities arise in the effective description of fluctuations presented in Section 4.1. We refer the reader to this section for details.

**Non-linearly realised symmetry.** The mixing action in the unitary gauge (4.4) is parameterised by a set of mass scales $\tilde{M}_i(t)$. After reintroducing the Goldstone boson $\pi$, expanding the mass scales to first order around the inflationary attractor gives

$$ \tilde{M}_i(t) \rightarrow \tilde{M}_i(t + \pi) \approx \tilde{M}_i(t) + \dot{\tilde{M}}_i(t)\pi + \dots \qquad (8.1) $$

If the time scales of variation of the $\tilde{M}_i(t)$ are too small, i.e. if the features that are considered are too sharp, the expansion (8.1) may break down. This corresponds to a violation of perturbative unitarity, see e.g [148,161–163]. Concentrating for simplicity on $c_s = 1$ and at weak mixing, the (minimum) time scale of variation $\Delta t$ of any time-dependent quantity should obey $\Delta t > 1/f_\pi$. This comes from requiring that the Taylor expansion in (8.1) is valid, i.e. $\pi < \Delta t$, with $\pi = \pi_c/f_\pi^2$ and $\pi_c \sim \omega$, at the energy $\omega \sim 1/\Delta t$ excited by the feature. Only a time-dependence in $\tilde{M}_1$ generates an additional cubic interaction compared to the ones present in (4.5), that reads

$$ \mathcal{L}_{\pi-\sigma}/a^3 \supset 6\tilde{M}_1^2(t)\dot{\tilde{M}}_1(t)\pi\dot{\pi}\sigma. \qquad (8.2) $$

For this reason, we will only consider a time dependence for $\tilde{M}_1$ throughout this section, setting $\tilde{M}_2 = \tilde{M}_3 = 0$. We also set the intrinsic sound speed of $\pi$ to unity $c_s = 1$. In terms of the canonically normalised field $\pi_c$, the mixing Lagrangian we study is

$$ \mathcal{L}_{\pi-\sigma}/a^3 = \rho(t)\dot{\pi}_c\sigma - \frac{\rho(t)}{2f_{\pi}^2}(\partial_{\mu}\pi_c)^2\sigma + \frac{\dot{\rho}(t)}{f_{\pi}^2}\pi_c\dot{\pi}_c\sigma, \qquad (8.3) $$

where we have defined the quadratic mixing $\rho(t) = 2\tilde{M}_1^3(t)/f_\pi^2$. Notice that operators weighted by higher derivatives of $\rho$ appear if $\tilde{M}_1(t + \pi)$ is expanded to higher order around $\tilde{M}_1(t)$. We discard such terms as they appear at quartic order (and higher) in the fields. Also note that we
have neglected any time variation of the Hubble parameter $H$ and the mass of the field $\sigma$ in this effective description of fluctuations. We will come back to this point in the next section.

**Effective single-field theory.** When the field $\sigma$ is sufficiently massive, one can integrate it out and obtain an effective theory for $\pi$. Indeed, solving the linear equation of motion for $\sigma$ to leading order in powers of $\Box/m^2$, that is $\sigma = \rho\dot{\pi}_c/m^2$, and plugging the solution back in the original Lagrangian leads to

$$ \mathcal{L}/a^3 = \frac{1}{2\tilde{c}_s^2(t)}\dot{\pi}_c^2 - \frac{1}{2}\frac{(\partial_i\pi_c)^2}{a^2} - \frac{1}{2f_\pi^2}\left(\frac{1}{\tilde{c}_s^2(t)} - 1\right)\dot{\pi}_c(\partial_\mu\pi_c)^2 - \frac{\dot{\tilde{c}}_s(t)}{f_\pi^2\tilde{c}_s(t)}\pi_c\dot{\pi}_c^2, \quad (8.4) $$

where the time-dependent sound speed $\tilde{c}_s(t)$ is defined in (4.30). Let us now discuss the regime of validity of this effective theory. Clearly, neglecting the gradient term compared to the mass term $\partial_i^2/a^2 \ll m^2$ means that (8.4) describes the dynamics of low-energy modes with physical momenta $k/a$ such that $k^2/a^2 \ll m^2$. In fact, it can be shown that requiring that the derivative expansion in powers of $\Box/m^2$ be under control implies a set of generalised adiabaticity conditions$^{47}$ [142,164–166]

$$ \frac{1}{\tilde{c}_s} \frac{d^n \tilde{c}_s}{dt^n} \ll m^n, \quad (8.7) $$

for $n \ge 1$. These conditions indicate that the time scale of variation $\Delta t$ of all time-varying quantities are bounded from below by the mass scale of the heavy field $\Delta t \gg m^{-1}$. In particular, this leaves room for the speed of sound to vary significantly on a time scale smaller than $H^{-1}$ but still larger than $m^{-1}$ if the field $\sigma$ is heavy enough.

In addition to the previous adiabaticity conditions, since the dispersion relation of the low-energy modes is $\omega^2 = \tilde{c}_s^2 k^2$, the derivative expansion is reliable if [27,150]

$$ \omega \ll \omega_{\text{new}} = \tilde{c}_s m, \quad (8.8) $$

where $\omega_{\text{new}}$ marks the energy scale at which higher-order derivative terms should be taken into account. Now, in order for the effective theory to have predictable power, one must require that $\omega_{\text{new}}$ be well above the Hubble scale. Indeed, if $\omega_{\text{new}} \ll H$, the effective theory would describe modes that are outside the sound horizon $k/a \ll m \ll H/\tilde{c}_s$. Such modes could not be initialised in the proper vacuum without the knowledge of the underlying UV completion. Therefore, for

$$ \Box \dot{\pi}_c = (1 - \tilde{c}_s^2) \frac{\partial_i^2 \dot{\pi}_c}{a^2} - 4 \tilde{c}_s^2 \frac{\dot{\tilde{c}}_s}{\tilde{c}_s} \frac{\partial_i^2 \pi_c}{a^2} + H^2 \left[ 6 \frac{\dot{\tilde{c}}_s}{\tilde{c}_s} - 2 \left( \frac{\dot{\tilde{c}}_s}{H\tilde{c}_s} \right)^2 + 3 \frac{\dot{H}}{H^2} - 2 \frac{\ddot{\tilde{c}}_s}{H^2 \tilde{c}_s} \right] \dot{\pi}_c. \quad (8.5) $$

Requiring $\Box\dot{\pi}_c \ll m^2\dot{\pi}_c$, so that the derivative expansion is under perturbative control, imposes some restrictions. Looking at the third term gives the following conditions

$$ \frac{\dot{\tilde{c}}_s}{\tilde{c}_s} \ll m, \quad \epsilon \ll \left(\frac{m}{H}\right)^2, \quad \frac{\ddot{\tilde{c}}_s}{\tilde{c}_s} \ll m^2, \quad (8.6) $$

whereas the first and the second terms are negligible as long as $k^2/a^2 \ll m^2$ (see e.g. [142] for more details).

<footnote>$^{47}$The first two conditions for $n=1,2$ can be easily found by considering the next-to-leading order term $\Box\dot{\pi}_c/m^2$ in the derivative expansion. Using the quadratic equation of motion for $\pi_c$, one obtains</footnote>
the effective theory to be valid, we also require

$$m \gg H/\tilde{c}_s. \qquad (8.9)$$

## 8.2 Correlated Sharp Features

Due to the non-linearly realised symmetry, features in the scale-dependence of the power spectrum are also imprinted in the bispectrum scale-dependence. This can be explicitly seen from the mixing Lagrangian in the full theory (8.3) or in the effective theory (8.4). Such correlated features have been extensively studied in the context of the single-field effective theory with a time-varying sound speed [167–170]. Here, we will perform exact calculations in both full and effective theories.

**Time dependence.** Although the mechanisms leading to features can be very model-dependent, we focus on a single form of the quadratic mixing. We choose its time dependence to be parameterised by

$$\rho(t) = \rho_0 (a/a_0)^{-n} \sin[\omega_c(t - t_0)] \times \frac{1}{2} \left[ 1 + \tanh \left( \frac{t - t_0}{\Delta t} \right) \right], \qquad (8.10)$$

where $\mu_c = \omega_c/H$ determines the background oscillation frequency in Hubble units, and $\rho_0$ is the amplitude of the feature. In the following, we will choose $\rho_0$ such that the deviation of the curvature power spectrum from scale invariance satisfies current observational bounds on cosmological scales. The feature is turned on continuously at $t = t_0$ with a step of width $\Delta t$ that dictates its sharpness. We also consider that the feature is diluted by $(a/a_0)^{-n}$ where $a_0 = a(t_0)$. This overall time dependence of the mixing describes well multi-field inflationary scenarios in which the background trajectory undergoes a turn in field space. If the turn is sharp, it will excite negative-frequency modes and produce correlation functions with oscillations linear in $k$ [171, 172]. Those oscillations are characteristic of any sharp feature, be it single field or multi-field, and they are independent of $\omega_c$. They only probe a very short time range located at the sharp feature. Depending on the values of the different parameters, the turn may also transiently excite a background massive field, leading to damped oscillations until the trajectory settles down along the light direction, corresponding to a dilution parameter $n = 3/2$. For $\omega_c \gtrsim H$, these background oscillations resonate with the quantum oscillations of fluctuations inside the Hubble radius, leading to oscillations linear in $\log k$ and of frequency $\omega_c$ [126]. Resonant features are particularly interesting as they can tell apart the inflationary scenario from alternative ones, as exemplified by the primordial standard clocks program [38, 172–179]. Both sharp and resonant features may also be encountered in single-field scenarios with either transient (bump, dip), or ever-present (tiny oscillatory component) deviations from the slow-roll dynamics, see for instance [24, 125, 126, 180, 181].

Before presenting numerical results, we briefly digress to discuss the purpose and limitations of our approach. From an EFT point of view, all couplings of the theory are generic functions of time that need to be inferred from data. However, notice that specific background models result in general in relationships between the time-dependencies of these couplings, whose knowledge can be used to parameterise the EFT couplings. In the following, for simplicity, we adopt a
<figure><img src="image_12.png" /><figcaption>Figure 11: The relative power spectrum to the featureless one $\Delta_{\zeta,0}^2$ and the bispectrum $f_{\text{NL}}$ in the equilateral configuration $k_1 = k_2 = k_3 = k$ as functions of $k/k_0$, computed in the full theory (solid lines) and in the effective single-field theory (dashed lines) both for $\mu = 2$ (left panel) and $\mu = 5$ (right panel). The feature, parameterised by $\rho_0/H = 20$, $\mu_c = 0.1$, $n = 3/2$, and $\Delta t = 0.1H^{-1}$, is localised at $t = t_0$ which corresponds to the scale $k = k_0$.</figcaption></figure>

modelling where only the time dependence of the quadratic mixing is taken into account. First, we assume that the time dependence of the quadratic mixing is not reflected in the evolution of the background spacetime, that is we assume that there is no variation of the Hubble parameter. Moreover, we do not assume any time dependence for the mass of the fluctuation $\sigma$, although in concrete multi-field realisations it is known to contain—in addition to a contribution from the second derivative of the potential—an explicit dependence on the quadratic mixing and the background field-space geometry [182–184].$^{48}$ Therefore, our approach disregards effects such as resonances from the time variation of the mass of the heavy field fluctuation with the same frequency $\omega_c$. Of course, it would be possible to solve the background dynamics and feed CosmoFlow with the actual time evolution of all varying quantities. This is precisely the purpose of the already existing tools PyTransport and CppTransport [122–124]. We instead embrace an effective description of fluctuations, having the advantage of rendering the symmetries manifest, and focusing directly on observables rather than the underlying model building.

**Numerical results.** For $(\omega_c, 1/\Delta t) \ll m$, the adiabaticity conditions (8.7) are satisfied and we expect that the effects coming from the heavy physics are well captured by a reduced speed of sound, provided that (8.9) holds. Whenever $\Delta t \ll m^{-1}$ or $\omega_c \gg m$, the effective theory is expected to break down. Here, we focus on the particularly interesting regime $\Delta t \lesssim m^{-1}$ while maintaining $\mu_c \ll \mu \equiv \sqrt{m^2/H^2 - 9}/4$. Figure 11 shows the numerical results for $\rho_0/H = 20$, $\mu_c = 0.1$, $\Delta t = 0.1H^{-1}$ and $n = 3/2$, both for $\mu = 2$ and $\mu = 5$. These parameters describe a sharp feature where the resonant effects from the oscillations, that occur only for $\mu_c \gtrsim 1$, are absent. We display the power spectrum (blue line) and the bispectrum in the equilateral configuration (red line), computed both in the full (solid line) and the effectively single-field

<footnote>$^{48}$From $m^2(t)\sigma^2 \rightarrow m^2(t+\pi)\sigma^2 \approx (m+\dot{m}\pi)^2\sigma^2$, a non-trivial time variation of the mass of the field $\sigma$ would generate the cubic interaction $\mathcal{L}_{\pi-\sigma}/a^3 \supset (\dot{m}^2)/(2f_{\pi}^2)\pi_c\sigma^2$ (see [142] for the complete cubic action for inflationary fluctuations written in comoving gauge in non-linear two-field sigma model, and [143] for the generalisation to any number of fields, where it is explicitly shown that this operator is present.).</footnote>