# Spectral Representation of Cosmological Correlators

Denis Werth

► To cite this version:

Denis Werth. Spectral Representation of Cosmological Correlators. Journal of High Energy Physics, 2024, 12, pp.017. 10.1007/JHEP12(2024)017. hal-04701154

HAL Id: hal-04701154

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

Submitted on 16 Sep 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
<figure><img src="image_2.png" /></figure>

RECEIVED: September 12, 2024

ACCEPTED: November 5, 2024

PUBLISHED: December 3, 2024

# Spectral representation of cosmological correlators

## Denis Werth

Institut d'Astrophysique de Paris, Sorbonne Université, CNRS,
Paris, F-75014, France

E-mail: werth@iap.fr

**ABSTRACT:** Cosmological correlation functions are significantly more complex than their flat-space analogues, such as tree-level scattering amplitudes. While these amplitudes have simple analytic structure and clear factorisation properties, cosmological correlators often feature branch cuts and lack neat expressions. In this paper, we develop off-shell perturbative methods to study and compute cosmological correlators. We show that such approach not only makes the origin of the correlator singularity structure and factorisation manifest, but also renders practical analytical computations more tractable. Using a spectral representation of massive cosmological propagators that encodes particle production through a suitable $i\epsilon$ prescription, we remove the need to ever perform nested time integrals as they only appear in a factorised form. This approach explicitly shows that complex correlators are constructed by gluing lower-point off-shell correlators, while performing the spectral integral sets the exchanged particles on shell. Notably, in the complex mass plane instead of energy, computing spectral integrals amounts to collecting towers of poles as the simple building blocks are meromorphic functions. We demonstrate this by deriving a new, simple, and partially resummed representation for the four-point function of conformally coupled scalars mediated by tree-level massive scalar exchange in de Sitter. Additionally, we establish cosmological largest-time equations that relate different channels on in-in branches via analytic continuation, analogous to crossing symmetry in flat space. These universal relations provide simple consistency checks and suggest that dispersive methods hold promise for developing cosmological recursion relations, further connecting techniques from modern scattering amplitudes to cosmology.

**KEYWORDS:** de Sitter space, Correlation Functions, Cosmological models, Early Universe Particle Physics

ARXIV EPRINT: 2409.02072
# Contents

<table><tbody><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>From time ordering to contour integrals</td><td>4</td></tr><tr><td>2.1</td><td>Cosmological propagators</td><td>4</td></tr><tr><td>2.2</td><td>Multiple propagators: a matter of contour</td><td>6</td></tr><tr><td>2.3</td><td>Spectral representation of massive cosmological propagators</td><td>8</td></tr><tr><td>3</td><td>Cosmological largest-time equation</td><td>11</td></tr><tr><td>3.1</td><td>Propagator identity from an off-shell perspective</td><td>12</td></tr><tr><td>3.2</td><td>Largest-time equation</td><td>13</td></tr><tr><td>3.3</td><td>Application to simple diagrams</td><td>18</td></tr><tr><td>4</td><td>Massive exchange correlator in de Sitter</td><td>24</td></tr><tr><td>4.1</td><td>Off-shell three-point function</td><td>25</td></tr><tr><td>4.2</td><td>Bootstrapping via the spectral representation</td><td>26</td></tr><tr><td>4.3</td><td>Singularities and analytic structure</td><td>30</td></tr><tr><td>5</td><td>Conclusions</td><td>31</td></tr><tr><td>A</td><td>Feynman $i\epsilon$ prescription from vacuum wave-functional</td><td>33</td></tr><tr><td>B</td><td>Details on massive fields in de Sitter</td><td>34</td></tr><tr><td>B.1</td><td>Useful formulae</td><td>34</td></tr><tr><td>B.2</td><td>Derivation of the spectral representation</td><td>37</td></tr><tr><td>B.3</td><td>Non-analyticity and dispersive integral in the energy domain</td><td>37</td></tr></tbody></table>

## 1 Introduction

The complexity of real physical processes is often made simple when these processes are extended to the complex plane. While thinking about particles propagating with complex momenta and scattering at imaginary angles may seem grotesque, this perspective can offer profound insights into the real world. Even more remarkably, fundamental physical principles themselves can be read off from analytic properties of scattering amplitudes.

The development of off-shell methods for scattering amplitudes is deeply rooted in the power of complex analysis. This approach traces its origins back to the celebrated Kramers-Kronig relations, which link absorption and dissipation in a medium by relating the real and imaginary parts of the response function, thereby establishing the foundational connection between causality and analyticity. Building on this foundation, the complex deformation of momenta has led to significant developments, such as Berends-Giele recursion relations [1], recovering tree-level amplitudes by stitching together lower-point building blocks with one leg off-shell, or crossing symmetry, revealing that different processes are related through
analytic continuation. From this perspective, further rethinking kinematic constraints and symmetries has led to the discovery of fascinating internal mathematical beauty of scattering amplitudes. Methods such as the use of spinor-helicity variables, improved on-shell recursion relations including BCFW [2, 3], or unitarity methods and generalised cuts [4–6] not only deepen our understanding of quantum field theory but also facilitate the computation of complicated physical processes, as evidenced by recent achievements in high-loop calculations (see [7–10] for nice reviews).

In the context of primordial cosmology, the foundations of cosmological correlation functions have not yet achieved the same degree of understanding. Pushing forward the state-of-the-art of computational techniques, even at first order in perturbation theory, demands gruelling mathematical dexterity, for at least three interrelated reasons. First, computing equal-time correlation functions requires integrating over the entire bulk time evolution, as time translation symmetry is lost and an asymptotically free future is absent. Second, the curved geometry of the cosmological background distorts the free propagation of particles. Finally, unlike tree-level scattering amplitudes, which are meromorphic functions of complex kinematics with branch cuts appearing only at the loop level, tree-level cosmological correlators already exhibit branch cuts. This crucial distinction arises from spontaneous particle production.

Despite these challenges, several powerful techniques have recently been developed. At tree level, approaches such as bootstrapping correlators by solving boundary kinematic equations [11–17] and reformulating perturbative diagrammatic rules in Mellin space [18–21] have proven effective. The use of partial Mellin-Barnes representations [15, 22, 23], leveraging the global dilatation symmetry of de Sitter space that replaces time translations, has further improved these methods. Additionally, numerical techniques that trace the time evolution of correlation functions have been developed to explore and carve out the space of most correlators [24–26]. At the loop order, dispersive integrals have been used to tame loop diagrams [27, 28] (also see [29–35] for other techniques).<fnref n="1" />

At the same time, a deeper understanding of the structure of cosmological correlators is gradually emerging as insights from flat-space scattering amplitudes are imported to cosmology [36]. Essentially, flat-space amplitudes are embedded within cosmological correlators as residues on the total energy pole, when energies of the external particles add up to zero [37, 38]. Similarly, on partial energy poles, when the energy of a subdiagram vanishes, correlators factorise into a product of a (shifted) lower-point correlator and a lower-point scattering amplitude. Notably, these singularities can only be reached by analytically continuing (some) external energies to negative values. Perturbative unitarity leads to a set of cutting rules [39–43], while its non-perturbative counterpart is captured through the Källén-Lehmann representation in de Sitter space [20, 44, 45]. This opens up new possibilities to sew together simple building blocks to construct more complex correlators. Furthermore, simple cosmological correlators can be conceptualised as canonical forms of a polytope [46–53] and progress has been made in understanding the differential equations they satisfy, from which “time emerges” [54–57] (with related aspects discussed in [58–60]). As

<footnote n="1">¹The resummation of massless loops in de Sitter space — and the associated secular infrared divergences — has a rich and extensive story that goes beyond our focus in this work.</footnote>
we will argue in this paper, employing an *off-shell* formulation of cosmological correlators not only makes their analytic and factorisation properties more explicit, but also streamlines practical computational techniques.

Most of the analytic properties and the encoding of physical principles are phrased at the level of wavefunction coefficients, which loosely are the cosmological counterparts of scattering amplitudes.<fnref n="2" /> Yet, it remains unclear how these properties can be translated to actual observables, at least beyond leading orders in perturbation theory. Interestingly, cosmological correlators, when pushed to higher order in perturbation theory, are structurally simpler than wavefunction coefficients. For example, they share the same transcendentality as their corresponding flat-space amplitudes [62]. Additionally, these objects are most directly connected to cosmological data as these are the main observables. As such, they can be viewed as the cosmological analogue of cross sections.

In this paper, we propose a systematic off-shell study of cosmological correlation functions. Our approach relies on a spectral representation of massive bulk-to-bulk propagators in de Sitter, as first introduced in [63] (with previous spectral representations in de Sitter proposed in [20, 64–66]).<fnref n="3" /> Essentially, time ordering is replaced by an appropriate contour integral. While spontaneous particle production makes it impossible to use the usual complex energy plane because of branch cuts, the key insight is instead to shift to the complex mass domain, where the dispersive integral takes the form of a split representation of off-shell bulk-to-boundary propagators. In this plane, mode functions are analytic and the integration contour can be safely closed. This approach uses a suitable $i\epsilon$ prescription to account for particle production. By incorporating particle production, unlike the usual Feynman propagator, the number of poles effectively doubles, with different Boltzmann weights that either suppress or enhance their effects. This encodes a crossing symmetry between incoming and outgoing modes separated by a Stokes line, fully capturing particle production. In the flat-space limit, where particle production is turned off, this propagator naturally reduces to the standard Feynman propagator. At the level of correlators, this off-shell propagator trivialises time integrals as they only appear in a factorised manner.

On a practical level, we explicitly compute the spectral integral for the four-point function of conformally coupled scalars mediated by the tree-level exchange of massive scalars, and introduce a new, partially resummed representation for this correlator. This new approach also offers a more intuitive understanding of the underlying physics. For example, the single tower of poles that we find maps directly to the quasi-normal modes of the exchanged massive field. Fundamentally, as we will show, this dispersive representation of cosmological correlators clearly illustrates the analytic structure of individual contributions and their factorisation properties. It makes explicit that exchange correlators are constructed by gluing lower-point off-shell correlators with a sewing kernel that encodes both particle production (from the off-shell leg) and a tower of EFT contributions that emerge after integrating out the exchanged field. By explicitly performing this spectral integral, we set the exchanged particle on-shell, allowing us to reconstruct more complex diagrams.

<footnote n="2">²Still, wavefunction coefficients are sensitive to total derivatives in the action and are not invariant under field redefinitions. To cure this, an *S*-matrix for de Sitter space has been recently proposed in [61].</footnote>

<footnote n="3">³Integral representations of cosmological propagators have been proposed in e.g. [67–69], mostly for the wavefunction coefficients. But all these simple representations do not capture particle production as they were derived either for flat-space objects or for the special case of a conformally coupled field in de Sitter.</footnote>
Finally, we leverage the fundamental properties of in-in Schwinger-Keldysh diagrammatics to derive the largest-time equation satisfied by cosmological correlators. This equation relates processes occurring on different in-in branches through the analytic continuation of external energies, mirroring crossing symmetry in flat space. This relationship is universally applicable to all cosmological correlators, regardless of the order in perturbation theory or whether the theory is unitary, and it can serve as a consistency check when computing complicated diagrams. We conclude by illustrating how the largest-time equation is satisfied through a series of explicit examples.

**Outline.** The outline of the paper is as follows: in section 2, we briefly review the Schwinger-Keldysh formalism and cosmological propagators, with a particular focus on encoding time ordering as a contour integral. We then introduce the spectral representation of massive propagators in de Sitter and discuss their properties. In section 3, we derive the cosmological largest-time equation and prove it for both tree-level and loop-level graphs. For simple correlators, we show how their analytic structure, factorisation properties, and consistency requirements like the largest-time equation become clear when viewed from an off-shell perspective. In section 4, we explicitly perform the spectral integral to obtain a new representation of a massive exchange correlator in de Sitter. We summarise our conclusions in section 5.

**Notation and conventions.** We use the mostly-plus signature for the metric (-, +, +, +). Spatial three-dimensional vectors are written in boldface **k**. We use sans serif letters (a, b, ...) for Schwinger-Keldysh indices. External and internal energies of a graph are denoted by $E_i$ and $K_i$, respectively. Cosmic (physical) time is denoted by $t$ and conformal time, such that $d\tau = dt/a$, by $\tau$. Throughout most of the paper, we work in the Poincaré patch of de Sitter space with the metric $ds^2 = a^2(\tau)(-d\tau^2 + d\boldsymbol{x}^2)$, where $a(\tau) = -(H\tau)^{-1}$ is the scale factor and $H$ is the Hubble parameter. We will often set the Hubble scale to unity $H = 1$. Additional definitions will be introduced as needed in the main text.

## 2 From time ordering to contour integrals

We begin by introducing Schwinger-Keldysh propagators, which play a crucial role in calculating cosmological correlators. We then present the associated spectral representations, which encode time ordering through contour integrals in the complex plane.

### 2.1 Cosmological propagators

We are interested in connected equal-time correlation functions of a bulk scalar field $\varphi(t, \boldsymbol{x})$. The field $\varphi$ should be thought of as the fluctuation around a classical background, as usual in cosmological settings. In the following, we consider a weakly coupled field theory described by a Lagrangian density $\mathcal{L}[\varphi]$.

**Generating functional.** Following the standard Schwinger-Keldysh path integral formalism to compute correlation functions [70], two copies $\varphi_{\pm}$ of the field $\varphi$ are introduced on the two branches of the in-in, with the (+) branch representing the forward evolution from an initial time $t_0$ to $t$ and the (-) branch the backward evolution from $t$ to $t_0$. The condition that both
branches are sewn at the time t imposes that the field configurations $\varphi_{\pm}$ should coincide at t, i.e. $\varphi_{+}(t, \boldsymbol{x}) = \varphi_{-}(t, \boldsymbol{x})$. Equal-time correlation functions $\langle\Omega|\hat{\varphi}_{a_1}(t, \boldsymbol{x}_1)\dots\hat{\varphi}_{a_n}(t, \boldsymbol{x}_n)|\Omega\rangle$, with $a_1, \dots, a_n = \pm$ and where $|\Omega\rangle$ is the vacuum of the fully interacting theory, are computed by taking functional derivatives

$$ \langle\Omega|\hat{\varphi}_{a_1}(t, \boldsymbol{x}_1)\dots\hat{\varphi}_{a_n}(t, \boldsymbol{x}_n)|\Omega\rangle = \frac{1}{i\boldsymbol{a}_1\dots i\boldsymbol{a}_n} \left. \frac{\delta^n Z[J_+, J_-]}{\delta J_{a_1}(t, \boldsymbol{x}_1)\dots\delta J_{a_n}(t, \boldsymbol{x}_n)} \right|_{J_{\pm}=0}, \quad (2.1) $$

where the generating function $Z[J_+, J_-]$ is defined as

$$ Z[J_+, J_-] = \int_{\mathcal{L}_{i\epsilon}} \mathcal{D}\varphi_{\pm} \exp \left( i \int_{-\infty}^{t} dt' d^3x (\mathcal{L}[\varphi_+] + J_+\varphi_+ - \mathcal{L}[\varphi_-] - J_-\varphi_-) \right). \quad (2.2) $$

The contour $\mathcal{L}_{i\epsilon}$ denotes the usual in-in contour that goes in the upper-half complex plane from $-\infty^+$ to $t$ and reverts towards $-\infty^-$ in the lower-half complex plane, with $-\infty^{\pm} \equiv -\infty(1 \mp i\epsilon)$ and $\epsilon > 0$ an infinitesimal parameter. The vacuum can only be defined in the asymptotic past when fluctuations have wavelengths (in Fourier space) much smaller than the size of the cosmological horizon. As such, in eq. (2.2), the contour encodes the standard $i\epsilon$ prescription to which we will come back later.

**Free propagators.** Resorting to perturbation theory, the full Lagrangian is split into a free part $\mathcal{L}_0$ that can be solved exactly, defining the free propagators, and an interacting part $\mathcal{L}_{\text{int}}$, i.e. $\mathcal{L} = \mathcal{L}_0 + \mathcal{L}_{\text{int}}$. The first simplest object that can be defined is the free Wightman function $\langle 0|\hat{\varphi}(t_1, \boldsymbol{x}_1)\hat{\varphi}(t_2, \boldsymbol{x}_2)|0\rangle$ where the expectation value is over the vacuum of the free theory $|0\rangle$. In Fourier space for the spatial coordinates only, after expanding the operator $\hat{\varphi}_k(t)$ in terms of annihilation and creation operators

$$ \hat{\varphi}_k(t) = u_k(t) \hat{a}_k + u_k^*(t) \hat{a}_{-k}^\dagger, \quad (2.3) $$

where $u_k(t)$ is the mode function, the Wightman function is given by

$$ \mathcal{G}(k; t_1, t_2) = u_k(t_1) u_k^*(t_2). \quad (2.4) $$

We have used invariance under spatial translations and rotations for a homogeneous and isotropic background so that the function $\mathcal{G}$ only depends on a single 3-momentum magnitude $k \equiv |\boldsymbol{k}|$. Although this object does not itself have a natural physical interpretation, it is the building block for the four types of Schwinger-Keldysh propagators $\mathcal{G}_{ab}(k; t_1, t_2)$ where $a, b = \pm$, defined as

$$ \begin{aligned} \mathcal{G}_{++}(k; t_1, t_2) &= \mathcal{G}(k; t_1, t_2)\Theta(t_1 - t_2) + \mathcal{G}^*(k; t_1, t_2)\Theta(t_2 - t_1), \\ \mathcal{G}_{+-}(k; t_1, t_2) &= \mathcal{G}^*(k; t_1, t_2), \\ \mathcal{G}_{-+}(k; t_1, t_2) &= \mathcal{G}(k; t_1, t_2), \\ \mathcal{G}_{--}(k; t_1, t_2) &= \mathcal{G}^*(k; t_1, t_2)\Theta(t_1 - t_2) + \mathcal{G}(k; t_1, t_2)\Theta(t_2 - t_1). \end{aligned} \quad (2.5) $$

These propagators being not independent is a consequence of the Wightman function satisfying the property $\mathcal{G}(k; t_2, t_1) = \mathcal{G}^*(k; t_1, t_2)$, as can be seen from eq. (2.4). From these “bulk-to-bulk” propagators, one can define “bulk-to-boundary” propagators by sending $t_2 \to t$
(for $t_1 < t$) and using the condition that the two in-in branches coincide at the external time $\varphi_+(t, \mathbf{x}) = \varphi_-(t, \mathbf{x})$. They read

$$ \mathcal{K}_+(k; t_1) = \mathcal{G}^*(k; t_1, t), \quad \mathcal{K}_-(k; t_1) = \mathcal{G}(k; t_1, t). \qquad (2.6) $$

Following the diagrammatic rules exposed in [71], we represent these propagators with a black dot ● denoting +, a white dot ○ denoting −, and a square □ denoting the boundary at external time $t$ where + and − are indistinguishable:

$$ \begin{gathered} \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \draw (-.5,0) -- (0,0); \node[above] at (-.5,0) {$t_1$}; \node[above] at (0,0) {$t_2$}; \draw (-.5,0) -- (-.5,1); \node[left] at (-.5,1) {$\bullet$}; \draw (0,0) -- (0,1); \node[right] at (0,1) {$\bullet$}; \end{tikzpicture} = \mathcal{G}_{++}(k; t_1, t_2), \\ \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \draw (-.5,0) -- (0,0); \node[above] at (-.5,0) {$t_1$}; \node[above] at (0,0) {$t_2$}; \draw (-.5,0) -- (-.5,1); \node[left] at (-.5,1) {$\bullet$}; \draw (0,0) -- (0,1); \node[right] at (0,1) {$\bullet$}; \end{tikzpicture} = \mathcal{G}_{+-}(k; t_1, t_2), \quad \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \draw (-.5,0) -- (0,0); \node[above] at (-.5,0) {$t_1$}; \node[above] at (0,0) {$t_2$}; \draw (-.5,0) -- (-.5,1); \node[left] at (-.5,1) {$\bullet$}; \draw (0,0) -- (0,1); \node[right] at (0,1) {$\bullet$}; \node[draw, circle, fill=white, inner sep=1pt] at (0,0) {};} \end{tikzpicture} = \mathcal{K}_+(k; t_1), \\ \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \draw (-.5,0) -- (0,0); \node[above] at (-.5,0) {$t_1$}; \node[above] at (0,0) {$t_2$}; \draw (-.5,0) -- (-.5,1); \node[left] at (-.5,1) {$\bullet$}; \draw (0,0) -- (0,1); \node[right] at (0,1) {$\bullet$}; \end{tikzpicture} = \mathcal{G}_{-+}(k; t_1, t_2), \quad \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \draw (-.5,0) -- (0,0); \node[above] at (-.5,0) {$t_1$}; \node[above] at (0,0) {$t_2$}; \draw (-.5,0) -- (-.5,1); \node[left] at (-.5,1) {$\bullet$}; \draw (0,0) -- (0,1); \node[right] at (0,1) {$\bullet$}; \node[draw, circle, fill=white, inner sep=1pt] at (0,0) {};} \end{tikzpicture} = \mathcal{K}_-(k; t_1). \end{gathered} \qquad (2.7) $$

A double-coloured dot ● means either a black or a white dot.

## 2.2 Multiple propagators: a matter of contour

The previously defined propagators are nothing but Green's functions of the free equation of motion with appropriate boundary conditions, and therefore can be written as contour integrals.

**Feynman propagator.** To illustrate how a frequency-space representation of the propagators can be constructed and for the sake of simplicity, let us first consider a massless scalar field $\varphi$ in flat space. We will upgrade the presented integral representations of propagators to massive fields in de Sitter space in the next section. The corresponding positive-frequency mode function reads $u_k(t) = \frac{e^{-ikt}}{\sqrt{2k}}$. The generalisation to massive fields is straightforward after setting $k \to \omega_k = \sqrt{k^2 + m^2}$. The key insight is to realise that time ordering can be written as a frequency-space integral using the mathematical identity

$$ e^{-ik(t_1-t_2)}\Theta(t_1-t_2) + e^{+ik(t_1-t_2)}\Theta(t_2-t_1) = -2k \int_{-\infty}^{+\infty} \frac{d\omega}{2i\pi} \frac{e^{i\omega(t_1-t_2)}}{(\omega^2 - k^2)_{i\epsilon}}, \qquad (2.8) $$

where the $i\epsilon$ prescription is implemented as follows

$$ \frac{1}{(\omega^2 - k^2)_{i\epsilon}} \equiv \frac{1}{2k} \left[ \frac{1}{\omega - (k - i\epsilon)} - \frac{1}{\omega - (-k + i\epsilon)} \right] = \frac{1}{\omega^2 - k^2 + i\epsilon}, \qquad (2.9) $$

with $\epsilon > 0$ an infinitesimal positive parameter, and where in the second equality we have dropped terms of order $O(\epsilon^2)$ and set $2k\epsilon = \epsilon$, which is valid in the limit $\epsilon \to 0$. This procedure effectively sets the two poles slightly off the real axis in the complex plane. To see that we indeed recover the correct Feynman propagator, we suppose that $t_1 > t_2$ and close the contour in the upper-half complex plane. The first term in eq. (2.9) gives zero and only the second term selects the residue, effectively setting the exponential on shell,
## 2.3 Spectral representation of massive cosmological propagators

We now show how the previous integration contour prescription is modified when taking into account spontaneous particle production. For concreteness in what follows, we consider a real massive scalar field in de Sitter, see appendix B for more details.<fnref n="4" />

**Contour prescription.** Spontaneous particle production manifests itself as non-analyticity in the complex energy domain $k \to \omega$ once the momentum is set off shell. This renders impossible to construct an integral representation of massive cosmological propagators for the reason that implementing time ordering requires crossing the branch cut. This challenge can be alleviated by instead analytically continuing the massive mode function to the complex $\mu$ plane, allowing the field to acquire any complex mass, as shown in [63]. The integral representation of the time-ordered massive cosmological propagator reads<fnref n="5" />

$$ \mathcal{G}_{++}(k; \tau_1, \tau_2) = i \int_{-\infty}^{+\infty} d\nu \mathcal{N}_{\nu} \frac{u_k^*(\tau_1, \nu)u_k^*(\tau_2, \nu)}{(\nu^2 - \mu^2)_{i\epsilon}}, \quad (2.15) $$

where $\mathcal{N}_{\nu} \equiv \frac{\nu}{\pi} \sinh(\pi\nu)$ is the de Sitter density of states, $u_k^*(\tau, \mu)$ is the negative-frequency massive mode function defined in (B.2), and the contour prescription to select the correct pole structure is given by<fnref n="6" />

$$ \frac{1}{(\nu^2 - \mu^2)_{i\epsilon}} \equiv \frac{1}{2\sinh(\pi\mu)} \left[ \frac{e^{+\pi\mu}}{\nu^2 - \mu^2 + i\epsilon} - \frac{e^{-\pi\mu}}{\nu^2 - \mu^2 - i\epsilon} \right]. \quad (2.16) $$

Naturally, the anti-time ordered propagator is given by $\mathcal{G}_{--}(k; \tau_1, \tau_2) = \mathcal{G}_{++}^*(k; \tau_1, \tau_2)$. The detailed derivation of (2.15) is presented in appendix B.2. The integral over the spectral parameter $\nu$ can be viewed as an integration over the scaling dimension $\Delta = \frac{3}{2} + i\nu$ in the principal series that labels infinite-dimensional (scalar) representations of the de Sitter group, which defines eigenvalues of the quadratic Casimir operator $\Delta(3-\Delta)$, see e.g. [72–75]. In this sense, this integral is closely related to the de Sitter Källén-Lehmann representation albeit in spatial Fourier space [44, 45]. Such spectral representation has its analogue in anti-de Sitter space, see e.g. [76], and in other spacetimes [77, 78].

**Flat-space limit.** The spectral representation (2.15) extends the familiar flat-space prescription (2.9) by capturing spontaneous particle production while at the same time having a well-defined flat-space limit. Indeed, in the limit $H \to 0$, the Boltzmann suppression is turned off $e^{-\pi\mu} \to 0$ as $\mu \sim m/H$ for a fixed mass. The additional unusual poles $\nu = \pm\mu \pm i\epsilon$ vanish and the contour prescription reduces to the flat-space one (2.9)

$$ \frac{1}{(\nu^2 - \mu^2)_{i\epsilon}} \xrightarrow{\mu \to \infty} \frac{1}{\nu^2 - \mu^2 + i\epsilon}. \quad (2.17) $$

$$ \frac{1}{\nu^2 - \mu^2 + i\epsilon} = \frac{\Gamma(i\nu - i\mu - \epsilon)\Gamma(-i\nu - i\mu - \epsilon)}{\Gamma(i\nu - i\mu - \epsilon + 1)\Gamma(-i\nu - i\mu - \epsilon + 1)}. $$

<footnote n="4">⁴We consider fields in the principal series $m/H \ge 3/2$ so that $\mu \equiv \sqrt{m^2/H^2 - 9/4}$ is real. In particular, we will not deal with late-time secular divergences appearing when exchanging light fields.</footnote>

<footnote n="5">⁵The spectral parameter $\nu$ should not be confused with $\nu = i\mu$ that is commonly used to label light fields in the complementary series.</footnote>

<footnote n="6">⁶The explicit pole structure (and the corresponding residues) can be immediately obtained from the usual flat-space prescription (2.9) or by noticing that</footnote>
Following the derivation of (2.15) presented in appendix B.2, for $\tau_1$ being in the future of $\tau_2$, and using limiting forms given in appendix B.1, the massive propagator reduces to

$$ \mathcal{G}_{++}(k; \tau_1, \tau_2) = \frac{i}{2} \int_{-\infty}^{+\infty} d\nu \frac{\nu J_{i\nu}(-k\tau_1) H_{i\nu}^{(2)}(-k\tau_2)}{\nu^2 - \mu^2 + i\epsilon} = \frac{\pi}{2} J_{i\mu}(-k\tau_1) H_{i\mu}^{(2)}(-k\tau_2) \quad (2.18) \\ \xrightarrow{\mu \to \infty} \frac{e^{-i\mu(t_1-t_2)}}{\sqrt{2\mu}}, $$

where we have selected the single pole at $\nu = \mu - i\epsilon$ by closing the contour in the lower-half complex plane, and have set $H = 1$. We recover the properly-normalised flat-space propagator composed of plane waves oscillating in cosmic time $t$. In the flat-space limit, this propagator captures the free super-horizon evolution of the massive field when spatial gradients are negligible, leading to a dispersion relation dominated by the mass $\omega_k \sim \mu$. This provides insight into why analytically continuing the mass, rather than the momentum, allows for defining a spectral representation of massive fields.

**Particle production.** Particle production of massive fields manifests itself as the emergence of negative-frequency modes over time, when initialised with a positive-frequency mode as selected by the Bunch-Davies vacuum. This process leads to mode mixing, reflected in a non-zero average particle number observed at later times. The Bessel function defines the mode function in the infinite future when the massive field reaches a steady state, as

$$ v_k(\tau, \mu) = \frac{\Gamma(1 + i\mu)}{\sqrt{2\mu}} \left(\frac{k}{2}\right)^{-i\mu} J_{i\mu}(-k\tau) \xrightarrow{\tau \to 0} \frac{e^{-i\mu t}}{\sqrt{2\mu}}. \quad (2.19) $$

Notice in passing that the late-time limit $\tau \to 0$ equivalently also corresponds to the large-mass limit $\mu \to \infty$ as can easily be recovered using (B.13) and the Stirling formula $\Gamma(z) \sim e^{-z}z^z\sqrt{2\pi/z}$. This is not surprising as both variables are conjugate to each other. Decomposing this outgoing mode function onto the basis of ingoing ones, we obtain

$$ v_k(\tau, \mu) = \alpha_k u_k(\tau, \mu) + \beta_k u_k^*(\tau, \mu), \quad (2.20) $$

where the Bogolyubov coefficients are

$$ \alpha_k = \frac{\Gamma(1 + i\mu)}{\sqrt{2\pi\mu}} \left(\frac{k}{2}\right)^{-i\mu} e^{+\frac{\pi\mu}{2} - \frac{i\pi}{2}}, \quad \beta_k = \frac{\Gamma(1 + i\mu)}{\sqrt{2\pi\mu}} \left(\frac{k}{2}\right)^{-i\mu} e^{-\frac{\pi\mu}{2} + \frac{i\pi}{2}}, \quad (2.21) $$

from which we obtain the mean particle density of “out” excitations with wavenumber $k$ in the “in” vacuum state $|\beta_k|^2 = 1/(e^{2\pi\mu} - 1)$. We recover the usual Bose-Einstein distribution, when identifying the energy with the rest mass $\mu$ in a thermal bath set by the de Sitter temperature $T_{ds} = 1/2\pi$.

**Local and non-local parts.** The $i\epsilon$ prescription (2.16) projects ingoing particle states onto outgoing states, which enables the mixing of positive- and negative-frequency mode functions. This is made manifest by rewriting the off-shell ingoing mode function in terms of off-shell outgoing ones using the connection formula (B.9), so that the propagator can be written [11, 79, 80]

$$ \mathcal{G}_{++}(k; \tau_1, \tau_2) = \mathcal{G}_{++}^{\text{local}}(k; \tau_1, \tau_2) + \mathcal{G}_{++}^{\text{non-local}}(k; \tau_1, \tau_2), \quad (2.22) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of the contour prescription (2.16) in the complex $\nu$ plane. The usual poles $\nu = \pm\mu \mp i\epsilon$ effectively project the off-shell mode onto late-time positive-frequency modes, whereas taking the residues at the Boltzmann suppressed poles $\nu = \pm\mu \pm i\epsilon$ projects the off-shell mode onto the late-time negative-frequency modes. This interchange of ingoing and outgoing modes captures particle production. In the late-time limit, the $\nu \leftrightarrow -\nu$ symmetry of the off-shell modes completely fixes the contour prescription both for the local and non-local contributions, leading to the disappearance of time ordering.</figcaption></figure>

where both contributions read

$$ 
\begin{aligned}
\mathcal{G}_{++}^{\text{local}}(k; \tau_1, \tau_2) &= \frac{+i\pi}{4} \int_{-\infty}^{+\infty} d\nu \frac{\mathcal{N}_{\nu}}{\sinh^2(\pi\nu)} \frac{J_{i\nu}(-k\tau_1)J_{i\nu}^*(-k\tau_2) + (\nu \leftrightarrow -\nu)}{(\nu^2 - \mu^2)_{i\epsilon}}, \\
\mathcal{G}_{++}^{\text{non-local}}(k; \tau_1, \tau_2) &= \frac{-i\pi}{4} \int_{-\infty}^{+\infty} d\nu \frac{\mathcal{N}_{\nu}}{\sinh^2(\pi\nu)} \frac{e^{-\pi\nu}J_{i\nu}(-k\tau_1)J_{i\nu}(-k\tau_2) + (\nu \leftrightarrow -\nu)}{(\nu^2 - \mu^2)_{i\epsilon}}.
\end{aligned}
\quad (2.23) $$

The local part is interpreted as the accumulated dynamical phase of a single particle propagating since in the late-time limit we have

$$ J_{i\nu}(-k\tau_1)J_{i\nu}^*(-k\tau_2) \sim \frac{\sinh(\pi\nu)}{\pi\nu} \left(\frac{\tau_1}{\tau_2}\right)^{i\nu} \propto e^{-i\nu(t_1-t_2)}. \quad (2.24) $$

Notably, the local part is independent of the momentum so that it describes correlation at coincident points in position space. Assuming $\tau_1 > \tau_2$, we close the contour in the lower-half complex plane for the term $J_{i\nu}(-k\tau_1)J_{i\nu}^*(-k\tau_2)$ and in the upper-half complex plane for the complex conjugate piece, which leads to

$$ 
\begin{aligned}
\mathcal{G}_{++}^{\text{local}}(k; \tau_1, \tau_2) &= \frac{\pi}{4 \sinh^2(\pi \mu)} \left[ e^{+\pi \mu} J_{i\mu}(-k\tau_1) J_{i\mu}^*(-k\tau_2) + e^{-\pi \mu} J_{i\mu}^*(-k\tau_1) J_{i\mu}(-k\tau_2) \right] \\
&\xrightarrow[\tau_1, \tau_2 \to 0]{\Gamma(+i\mu)\Gamma(-i\mu)} \frac{\Gamma(+i\mu)\Gamma(-i\mu)}{4\pi} \left[ e^{+\pi\mu} \left(\frac{\tau_1}{\tau_2}\right)^{+i\mu} + e^{-\pi\mu} \left(\frac{\tau_1}{\tau_2}\right)^{-i\mu} \right].
\end{aligned}
\quad (2.25) $$

For the reverse time ordering, the contours should be closed in the opposite directions for both contributions. The doubling of poles projects the off-shell positive-frequency mode into positive- and negative-frequency ones as $(\tau_1/\tau_2)^{\pm i\nu} \rightarrow e^{\pm\pi\mu}(\tau_1/\tau_2)^{\pm i\mu} + e^{\mp\pi\mu}(\tau_1/\tau_2)^{\mp i\mu}$, which then
by the structure of propagators. This general property remains valid in the presence of dissipation, fluctuations and non-unitary time evolution.

In de Sitter space, we have seen in section 2.3 that the effect of particle production amounts to double the usual poles while assigning distinct weights $e^{\pm\pi\mu}$ to the associated residues. From the spectral representation (2.15), it is easy to observe that bulk-to-bulk propagators also satisfy the identity (3.2), which, in terms of integration contours, is illustrated by

$$ \frac{e^{+\pi\mu}}{e^{-\pi\mu}} \begin{array}{c} \bullet \\ \bullet \end{array} - \frac{e^{-\pi\mu}}{e^{+\pi\mu}} \begin{array}{c} \bullet \\ \bullet \end{array} = \frac{e^{-\pi\mu}}{e^{-\pi\mu}} \begin{array}{c} \bullet \\ \bullet \end{array} - \frac{e^{+\pi\mu}}{e^{+\pi\mu}} \begin{array}{c} \bullet \\ \bullet \end{array} \quad (3.4) $$

Crucially, note that particle production poles change their respective weight under complex conjugation. Sending $\mu \to \infty$, one immediately recovers the flat-space limit (3.3).<fnref n="8" /> The identity (3.2) can of course be recovered from the form of propagators written with the time-ordered usual representation (2.5) using $\Theta(t_1 - t_2) + \Theta(t_2 - t_1) = 1$. Here, we have derived this identity from an off-shell perspective.

## 3.2 Largest-time equation

Our aim is to translate the propagator identity (3.2) at the level of correlators, resulting in an identity satisfied by individual branch contributions of cosmological correlators. In the context of flat-space scattering amplitudes, this identity is known as the largest-time equation [90].<fnref n="9" />

The flat-space largest-time equation states that the sum of all these contributions must vanish

$$ F(\{x_i\}) + F^*(\{x_i\}) + F(\{x_i\}) = 0, $$

where $F$ denotes the usual Feynman diagram, $F^*$ its complex conjugate, and $F$ all other $2^n - 2$ contributions. This is often taken as the starting point to derive recursion relations, see e.g. [1, 2, 6], or cutting rules, see e.g. [86, 90, 91]. From the above colouring rules, it comes at no surprise that the Schwinger-Keldysh diagrammatic rules, and therefore cosmological correlators, are the most suited objects to phrase the largest-time equation.

<footnote n="8">⁸The combination $G_K \equiv G_{++} + G_{--}$ is nothing but the standard Keldysh propagator. The advanced, retarded and Keldysh propagators can be recovered after performing a rotation in the $\phi_\pm$ field basis, introducing a new pair of fields $\phi_r = \frac{1}{\sqrt{2}}(\phi_+ + \phi_-)$ and $\phi_a = \frac{1}{\sqrt{2}}(\phi_+ - \phi_-)$, which is referred to the Keldysh basis [88]. The advantage of this basis is that the physical interpretation of various propagators become transparent. Indeed, $G_R$ and $G_A$ characterise the dynamical properties of the particles in the system under consideration, while $G_K$ characterises their statistical distribution and encodes the mean particle density [89].</footnote>

<footnote n="9">⁹Consider $F(\{x_i\})$ an amputated (real-space) Feynman diagram with $n$ vertices in flat space. Let us perform the following manipulations on this diagram:</footnote>

<footnote>* Duplicate the Feynman diagram $2^n$ times by colouring vertices either in black or in white in all possible ways.</footnote>

<footnote>* For each vertex, a black-coloured vertex brings a factor $i$, and a white-coloured one brings a factor $(-i)$.</footnote>

<footnote>* The propagator between two black vertices is the usual Feynman propagator $\Delta_F(x - y) = \Theta(x^0 - y^0)\Delta^+(x - y) + \Theta(y^0 - x^0)\Delta^-(x - y)$, where $\Delta^\pm(x - y)$ are the usual positive and negative frequency commutation functions, while the propagator between two white vertices is taken to be $\Delta_F^*(x - y)$. The propagator between black and white vertices is replaced by $\Delta^+(x - y)$, and $\Delta^-(x - y)$ for white and black vertices. Clearly, the all-black diagram is the usual Feynman diagram while the all-white one is its complex conjugate.</footnote>
**Graph.** Computing correlators boils down to evaluating a set of multi-nested time integrals, whose number equals the number of vertices. Thus, instead of working with correlators, we define a graph with V vertices and N internal lines by its collection of external energies {$E_i$} ($i = 1, ..., V$) and internal energies {$K_j$} ($j = 1, ..., N$). The resulting graph integrals are found after removing all overall factors of $1/2k$ from bulk-to-boundary and bulk-to-bulk propagators, hence defining the rescaled propagators

$$ \tilde{\mathcal{K}}_a(k; t) \equiv 2k \mathcal{K}_a(k; t), \quad \tilde{\mathcal{G}}_{ab}(k; t, t') \equiv 2k \mathcal{G}_{ab}(k; t, t'), \qquad (3.5) $$

and overall coupling constants. We also multiply each graph integral by $(-i)^V$. Of course, in practice, the external energies are function of external momentum magnitudes. The main motivation is that the number of bulk-to-boundary propagators attached to a vertex, which we take to be massless in flat space and conformally coupled in de Sitter, does not matter, as they satisfy the trivial relation

$$ \prod_{i=1}^{n} \tilde{\mathcal{K}}_a(k_i; t) = \tilde{\mathcal{K}}_a \left( \sum_{i=1}^{n} k_i; t \right), \qquad (3.6) $$

where $n$ is the number of external legs of a vertex. The sum of all external momentum magnitudes is defined to be the energy of the vertex, i.e. $k_1 + ... + k_n = E$. Diagrammatically, we represent a graph the same way as a correlator (with each vertex carrying a Schwinger-Keldysh index) albeit with external legs removed, and instead labelling each vertex with its energy. For example, a one-site graph integral in flat space reads

$$ \bullet \quad E \quad \equiv (-i) \sum_{a=\pm} a \int_{-\infty^a}^{0} dt \, e^{iaEt}, \qquad (3.7) $$

which represents any contact correlator with arbitrary number of external legs. We also label internal lines with their corresponding internal energies.

**A simple illustrative example.** Before stating the general largest-time equation, we illustrate the consequence of the identity (3.2) on a simple example, namely the two-site chain. Even though the precise nature of the field $\phi$, interactions, and the background is unimportant for what follows, we consider polynomial self-interactions in flat space and set coupling constants to unity.

First, let us consider the sum of both time-ordered contributions with flipped external energies on the negative branch

$$ \bullet \quad E_1 \quad K \quad E_2 \quad \bullet + \quad \bullet \quad -E_1 \quad K \quad -E_2 \quad \bullet . \qquad (3.8) $$

It is given by

$$ - \int_{-\infty^+}^0 dt_1 \int_{-\infty^+}^0 dt_2 e^{iE_1 t_1} [\tilde{\mathcal{G}}_{++}(K; t_1, t_2) + \tilde{\mathcal{G}}_{--}(K; t_1, t_2)] e^{iE_2 t_2}. \qquad (3.9) $$

Note that flipping external energies amounts to exchanging both in-in branches. As such, the $i\epsilon$ prescription that makes the integrals converge needs to be modified. Therefore, we
The above sum contains $2^V$ terms. In the presence of spatial derivative interactions, parity-violating interactions, or with external spinning fields, correlators can explicitly depend on external momenta $k_i$, that can be contracted with e.g. polarisation tensors. In such cases, analytically continuing the energies $k_i \to -k_i$ also implies flipping the signs of the corresponding momenta $k_i \to -k_i$.<fnref n="10" />

**Proof at tree-level.** To prove this identity for tree-level graphs, we use an inductive approach, as any tree graph can be constructed iteratively by successively adding internal lines. We follow the same lines as in the original work [90] where Cutkosky rules were derived for amputated Feynman diagrams, also see [41, 83, 95].

As seen above, the largest-time equation is satisfied for the simple two-site graph (the even more simple case of a one-site graph is trivial). We now assume that the latest-time equation (3.15) is valid for a tree-level graph with $V$ vertices and prove that it still remains valid when we add an additional vertex to the graph. Therefore, given a graph with $V$ vertices, we choose a vertex (say the farthest right one) and denote its corresponding external energy by $E_i$. Of course, this vertex carries $\pm$ labels. For what follows, it is useful to introduce some diagrammatic notation. We denote by the grey circle

$$ \bullet \quad E_i \quad , \qquad (3.16) $$

the sum of all individual graph contributions with $V$ vertices, where the picked vertex $\bullet$ has been brought forward. We now attach to it an additional vertex with external energy $E$ and we label the internal energy flowing through the bulk-to-bulk propagator by $K$. This increases the number of individual contributions by a factor 2, i.e. from $2^V$ to $2^{V+1}$. Similarly to the simple illustrative example that we have treated above, we now consider the special combination of attaching both time-ordered bulk-to-bulk propagators with negative energies on the negative branch, which pictorially is represented by

$$ \bullet \quad E_i \quad \bullet + \quad \bullet \quad -E_i \quad \bullet \quad . \qquad (3.17) $$

Note that we have added an additional factor of $i$ (which we include in the grey circle as it does not spoil the argument). Adding this vertex also brings an additional time integral. Flipping the external energies on the negative branch allows us to factor out all the information about external lines (and group both terms under same time integrals after also flipping

$$ \sum_{k=0}^{n} (-1)^k \sum_{\sigma \in \Pi(k, n-k)} \bar{T}[\mathcal{O}(t_{\sigma_1}) \cdots \mathcal{O}(t_{\sigma_k})] T[\mathcal{O}(t_{\sigma_{k+1}}) \cdots \mathcal{O}(t_{\sigma_n})] = 0, $$

where $\bar{T}$ and $T$ are the (anti-)time-ordering operators, and $\Pi(k, n-k)$ denotes the set of bi-partitions of $\{1, \dots, n\}$ with size $k$ and $n-k$. This “operator optical theorem” can be proved by induction. This means that the identity relating various in-in branches holds more generally at the operator level.

<footnote n="10">¹⁰It should be pointed out that the latest-time equation can be understood as a consequence of the following combinatoric operator identity [83, 92–94]</footnote>
the $i\epsilon$ prescription). Internal energies remain untouched. This manipulation isolates the bulk-to-bulk propagators to give

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E_i$}; \node[above] at (left) {$-E$}; \node[above] at (right) {$E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} \times \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E_i$}; \node[above] at (left) {$-E$}; \node[above] at (right) {$E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} . \quad (3.18) $$

This way, we can use the propagator identity (3.2) to replace all time-ordered propagators with non-time-ordered ones, yielding

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E_i$}; \node[above] at (left) {$-E$}; \node[above] at (right) {$E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} \times \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E_i$}; \node[above] at (left) {$-E$}; \node[above] at (right) {$E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} . \quad (3.19) $$

These terms can then be redistributed to generate two individual graphs. At this stage, we recognise the two factorised contributions with flipped external energies on the negative branch, up to an overall minus sign coming from the Schwinger-Keldysh vertex indices. We obtain

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E_i$}; \node[above] at (left) {$-E$}; \node[above] at (right) {$E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} \times \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E_i$}; \node[above] at (left) {$-E$}; \node[above] at (right) {$E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} . \quad (3.20) $$

Combining all terms on the same side, we finally obtain the desired diagrammatic identity. Since the simplest graph explicitly verifies this identity, the largest-time equation is proved for all tree-level graphs.

**Loop-level diagrams.** Cuts, in the broad sense, are at the heart of modern flat-space scattering amplitude methods, especially at the loop level. For cosmological correlators, we will see here that the largest-time equation naturally generalises to loop graphs, as this identity stems from general properties of bulk-to-bulk propagators.

A loop diagram can essentially be generated from a tree-level diagram by merging two vertices. However, this process involves addressing two key points. First, the propagator identity (3.2), which lies at the core of the largest-time equation at tree-level, must be extended to accommodate products of propagators. Second, loop diagrams contain unconstrained momenta flowing in the loops that need to be integrated over, which often results in divergences. As we are not primarily interested in loops, we will not explicitly perform loop integrals nor paying attention to regularising them. Instead, we will show how the largest-time equation holds at the loop level for a few simple graphs.

**One-site graph.** Let us consider the following simplest one-site graph at one-loop level with the vertex energy $E$

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E$}; \node[above] at (left) {}; \node[above] at (right) {}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} . \quad (3.21) $$

This diagram consists of two distinct contributions that are complex conjugate to each other. As usual, we sum these contributions after flipping the energy on the negative branch, resulting in

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$E$}; \node[above] at (left) {$-E$}; \node[above] at (right) {$E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} + \begin{tikzpicture}[baseline=(current bounding box.center)] \node[draw, circle, fill=gray!20, minimum size=1.5cm] (center) at (0,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (left) at (-1.5,0) {}; \node[draw, circle, fill=gray!20, minimum size=1.5cm] (right) at (1.5,0) {}; \draw[->] (center) -- (left); \draw[->] (center) -- (right); \node[above] at (center) {$-E$}; \node[above] at (left) {$E$}; \node[above] at (right) {$-E$}; \node[below] at (center) {}; \node[below] at (left) {}; \node[below] at (right) {}; \node[above right] at (center) {$K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \node[above right] at (center) {$-K$}; \node[below right] at (center) {}; \node[above left] at (center) {}; \node[below left] at (center) {}; \end{tikzpicture} . \quad (3.22) $$
This combination of graphs vanishes as both propagators are indistinguishable at equal times, i.e. $\tilde{G}_{++}(\ell; t, t) = \tilde{G}_{--}(\ell; t, t)$. Notice that this straightforwardly extends to all one-site graphs to any loop order.

Two-site graph. We then consider the two-site graph at one-loop level

$$ E_1 \bullet \bullet E_2 . \qquad (3.23) $$

Each contribution contains a product of bulk-to-bulk propagators. The key insight is to realise that bulk-to-bulk propagators satisfy

$$ \prod_{i=1}^{L+1} \mathcal{G}_{++}(k_i, t_1, t_2) + \prod_{i=1}^{L+1} \mathcal{G}_{--}(k_i, t_1, t_2) = \prod_{i=1}^{L+1} \mathcal{G}_{+-}(k_i, t_1, t_2) + \prod_{i=1}^{L+1} \mathcal{G}_{-+}(k_i, t_1, t_2), \quad (3.24) $$

where $L$ denotes the number of loops. Of course, this relation is also valid for the rescaled (tilted) propagators. For $L = 0$, we recover the identity (3.2) that was used in the tree-level case. Proving this identity is simply a matter of manipulating and rearranging time orderings. For unequal times $t_1 \neq t_2$, one further requires the use of the trivial formula $\Theta(t_1 - t_2)\Theta(t_2 - t_1) = 0$, and at equal times $t_1 = t_2$, one should notice that the Wightman function is real $\mathcal{G}(k_i, t, t) = \mathcal{G}^*(k_i, t, t)$. The special value $\Theta(0) = 1/2$ is used. At one-loop order $L = 1$, we therefore have

$$ E_1 \bullet \bullet E_2 + -E_1 \bullet \bullet -E_2 = -i \int_{-\infty+}^{0} dt_1 e^{iE_1 t_1} dt_2 e^{iE_2 t_2} \times \int \frac{d^3 \ell}{(2\pi)^3} [\tilde{G}_{++}(l; t_1, t_2) \tilde{G}_{++}(|\mathbf{k} - \mathbf{l}|; t_1, t_2) + \tilde{G}_{--}(l; t_1, t_2) \tilde{G}_{--}(|\mathbf{k} - \mathbf{l}|; t_1, t_2)] . \qquad (3.25) $$

Using (3.24), we recognise the sum of non-time-ordered contributions with flipped energies on the negative branch

$$ -E_1 \bullet \bullet -E_2 - -E_1 \bullet \bullet E_2 , \qquad (3.26) $$

albeit with an overall minus sign. Putting all terms on the same side, we recover the largest-time equation. Using (3.24) for $L > 2$, this naturally generalises to any loop order.

### 3.3 Application to simple diagrams

The primary challenge in performing perturbative calculations of cosmological correlators lies in the presence of nested time integrals, which are generally difficult to evaluate. However, as we have seen previously, time-ordering can be interpreted as a specific choice of contour to set the exchanged fields on shell. We can do this by considering a suitable integral representation for the bulk-to-bulk propagators, so that the off-shell propagators appear explicitly factorised.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of the analytic structure of the integrand of $\mathcal{I}_{++}$ in the complex $\omega$ plane. Different choices of closing the contour, or equivalently different time orderings, relate physical poles in the lower-half plane to unphysical poles in the upper-half plane.</figcaption></figure>

restored after changing variables $\omega \to -\omega$. The analytic structure of the integrand in the complex $\omega$ plane exhibits several interesting features that we show in figure 2: (i) the pole $\omega = +E_2 - i\epsilon$ is the off-shell collinear (folded) limit, (ii) the pole $\omega = -E_1 + i\epsilon$ is the off-shell (left) partial energy pole, and (iii) the poles $\omega = \pm(K - i\epsilon)$ are the on-shell conditions. Notice that poles that can only be reached after analytically continuing the energies to the unphysical kinematic configuration (after setting the exchanged particle on shell) are located in the upper-half complex $\omega$ plane. Therefore, selecting for example the physical on-shell condition requires closing the contour in the lower-half complex plane, picking the residues at $\omega = +K - i\epsilon$ and $\omega = +E_2 - i\epsilon$, we obtain

$$ \mathcal{I}_{++} = \frac{1}{E} \left( \frac{1}{E_L} + \frac{1}{E_R} \right), \qquad (3.31) $$

where $E \equiv E_1 + E_2$ is the total energy of the graph. Of course, closing the contour upwards yields the same result. This shows that partial energy singularities reached for unphysical kinematic configurations are related to folded singularities reached for physical ones by different time ordering choices. In this representation, the total energy pole $E \to 0$ is made explicit through the identity

$$ \frac{1}{(E_1 + \omega - i\epsilon)(E_2 - \omega - i\epsilon)} = \frac{1}{E} \left( \frac{1}{E_1 + \omega - i\epsilon} + \frac{1}{E_2 - \omega - i\epsilon} \right), \qquad (3.32) $$

and partial energy poles are reconstructed from the frequency integral. By reality of the result, we have $\mathcal{I}_{--} = \mathcal{I}_{++}$. Summing all contributions, the final result reads

$$ \mathcal{I}(\{E_i\}, K) = \frac{f(\{E_i\}, K)}{EE_L E_R}, \quad \text{with} \quad f(\{E_i\}, K) = 4(E_1 + E_2 + K). \qquad (3.33) $$

This graph integral is a rational function of the energies and has a singularity when the total energy vanishes $E \to 0$, and when partial energies flowing into the left or right vertices
This double integral can be evaluated after expanding the integrand into partial fractions, rendering the off-shell poles manifest, and then picking up residues. We obtain

$$ \mathcal{I}_{+++} = -\frac{1}{E} \left( \frac{1}{E_2 + E_3 + K_1} - \frac{1}{E_2 + E_3 - K_1} \right) \left( \frac{1}{E_3 + K_2} - \frac{1}{E_3 - K_2} \right) \\ - \frac{1}{(E_3 - K_2)(E_1 + E_2 + K_2)} \left( \frac{1}{E_2 + K_1 + K_2} - \frac{1}{E_2 - K_1 + K_2} \right) \\ - \frac{1}{(E_1 + K_1)(E_2 + E_3 - K_1)} \left( \frac{1}{E_3 + K_2} - \frac{1}{E_3 - K_2} \right) \\ - \frac{1}{(E_1 + K_1)(E_2 - K_1 + K_2)(E_3 - K_2)}. \qquad (3.39) $$

Notice that it is important to keep track of the correct $i\epsilon$ pole prescription to collect the correct residues. Eventually, since all contributions are real, the remaining ones are identical. Summing all contributions, the graph integral is

$$ \mathcal{I}(\{E_i\}, \{K_i\}) = \frac{f(\{E_i\}, \{K_i\})}{E E_L \bar{E}_L E_M \bar{E}_R E_R}, \qquad (3.40) \\ \text{with } f(\{E_i\}, \{K_i\}) = -8[\bar{E}_L E_M \bar{E}_R + E_M(E_1 \bar{E}_L + E_3 \bar{E}_R) + E_1 E_3 (\bar{E}_L + \bar{E}_R)]. $$

The fairly non-trivial form factor $f$ is well symmetric under $(E_1 \leftrightarrow E_3, K_1 \leftrightarrow K_2)$, and the final result displays the correct $E_L, E_M, E_R, \bar{E}_L, \bar{E}_R$ and $E$ poles, as expected. One can explicitly verify that individual contributions satisfy the largest-time equation

$$ \begin{aligned} 0 = \mathcal{I}_{+++}(E_1, E_2, E_3, K_1, K_2) + \mathcal{I}_{---}(-E_1, -E_2, -E_3, K_1, K_2) \\ + \mathcal{I}_{++}(-E_1, E_2, -E_3, K_1, K_2) + \mathcal{I}_{--}(-E_1, -E_2, E_3, K_1, K_2) \\ + \mathcal{I}_{-++}(-E_1, E_2, E_3, K_1, K_2) + \mathcal{I}_{+-}(-E_1, -E_2, -E_3, K_1, K_2) \\ + \mathcal{I}_{+-}(-E_1, -E_2, E_3, K_1, K_2) + \mathcal{I}_{-+-}(-E_1, E_2, -E_3, K_1, K_2), \end{aligned} \qquad (3.41) $$

providing a non-trivial check of the final result. In principle, this procedure can be generalised to arbitrary tree-level diagrams. However, the number of required integrals and the complexity of multi-dimensional frequency integrals rapidly become intractable.

**Conformal two-site chain in de Sitter.** As a last example, we consider the two-site chain of conformally coupled scalar fields (with mass $m^2 = 2H^2$) in de Sitter, whose computation conceptually follows the same line as in flat space, as the corresponding time-ordered bulk-to-bulk propagators can be recasted as frequency-space integrals.<fnref n="13" /> The corresponding graph integrals are

$$ E_1 \bullet K \bullet E_2 = \sum_{a,b=\pm} \mathcal{I}_{ab}, \quad \mathcal{I}_{ab} = -ab \int_{-\infty^a}^{\tau_0} \frac{d\tau_1}{\tau_1^2} \int_{-\infty^b}^{\tau_0} \frac{d\tau_2}{\tau_2^2} e^{i a E_1 \tau_1} \tilde{G}_{ab}(K; \tau_1, \tau_2) e^{i b E_2 \tau_2}, \qquad (3.42) $$

<footnote n="13">¹³It is possible to map correlators of conformally coupled fields to correlators of massive scalars with half-integer values of $\nu$ by applying a set of weight-shifting operators. This procedure follows from the property of the Hankel function that drastically simplifies for $\nu = n/2$ with $n$ integer ($n=1$ is the conformally coupled case).</footnote>
where $\tau_0 < 0$ is a small late-time cutoff. Shifting the energy and neglecting analytic terms in the late-time limit, we use the formula

$$ \int_{-\infty}^{\tau_0} \frac{d\tau}{\tau} e^{\pm i z \tau} = \int_{-\infty}^{\tau_0} \frac{d\tau}{\tau} e^{\pm i(z \mp i\epsilon)\tau} = \log[\pm i(z \mp i\epsilon)\tau_0], \qquad (3.43) $$

for a real positive energy $z > 0$, to write the factorised contribution as

$$ \mathcal{I}_{+-} = \log(+iE_L\tau_0) \log(-iE_R\tau_0). \qquad (3.44) $$

Similarly, after performing the time integrals for the nested contribution, we obtain

$$ \mathcal{I}_{++} = 2K \int_{-\infty}^{+\infty} \frac{d\omega}{2i\pi} \frac{\log[i(E_1 + \omega - i\epsilon)] \log[i(E_2 - \omega - i\epsilon)]}{(\omega^2 - K^2)_{i\epsilon}}. \qquad (3.45) $$

This integral, besides being formally divergent, is hard to evaluate because of the branch cuts spanning over $\omega \in (-E_1 + i\epsilon - i\infty, -E_1 + i\epsilon)$ and $\omega \in (E_2 - i\epsilon + i\infty, E_2 - i\epsilon)$ (choosing the log branch cut to be the principal one). These branch cuts notably cross the real axis, and arise from the cancellation of would-be IR divergences, as explicitly shown in footnote 14, that are otherwise absent in the presence of derivative interactions. As first noticed in [11], this integral can be written as an energy integral over the flat-space result

$$ \begin{aligned} \mathcal{I}_{++} &= 2K \int_{-\infty}^{+\infty} \frac{d\omega}{2i\pi} \frac{1}{(\omega^2 - K^2)_{i\epsilon}} \int_{E_1+\omega-i\epsilon}^{\infty} \int_{E_2-\omega-i\epsilon}^{\infty} \frac{dxdy}{xy} \\ &= 2K \int_{E_1}^{\infty} dx \int_{E_2}^{\infty} dy \int_{-\infty}^{+\infty} \frac{d\omega}{2i\pi} \frac{1}{(x + \omega - i\epsilon)(\omega^2 - K^2)_{i\epsilon}(y - \omega - i\epsilon)}, \end{aligned} \qquad (3.46) $$

where we have shifted the kinematic integrals so that their limits do not depend on the off-shell frequency. We recognise the flat-space frequency integral that we previously computed. The energy integrals were first computed in [11] and the result reads<fnref n="14" />

$$ \begin{aligned} \mathcal{I}_{++} ={}& \mathrm{Li}_2\left(\frac{E - E_L}{E}\right) + \mathrm{Li}_2\left(\frac{E - E_R}{E}\right) + \log\left(\frac{E_L}{E}\right)\log\left(\frac{E_R}{E}\right) - \frac{\pi^2}{6} \\ & - \log(iE_L\tau_0)\log(iE_R\tau_0), \end{aligned} \qquad (3.47) $$

where $\mathrm{Li}_2$ is the dilogarithm function. When summing all contributions, the kinematic parts of divergent pieces cancel against each other thanks to the property $\log(iz) - \log(-iz) = -i\pi$ for $z < 0$. Using the Euler's identity

$$ \mathrm{Li}_2(z) + \mathrm{Li}_2(1-z) + \log(z)\log(1-z) = \frac{\pi^2}{6}, \qquad (3.48) $$

$$ \int_{-\infty}^{\tau_0} \frac{d\tau_1}{\tau_1} \frac{d\tau_2}{\tau_2} e^{i(E_L \tau_1 + E_R \tau_2)} = \log(iE_L \tau_0) \log(iE_R \tau_0) = \int_{E_1}^{\infty} dx \int_{E_2}^{\infty} dy \frac{1}{(x+K-i\epsilon)(y+K-i\epsilon)}, $$

which can either be written as a product of logarithms (hence isolating the divergence) or as energy integrals that can be combined with the original ones. This effectively mimics the corresponding wavefunction coefficient calculation. The remaining convergent integral to be computed is

$$ \int_0^\infty \frac{dxdy}{(x+y+E)(x+E_L)(y+E_R)}. $$

Sophisticated techniques to evaluate such integrals are detailed in [54–57].

<footnote n="14">¹⁴In order to extract the divergence, one can subtract and add the term</footnote>
both nested contributions combine to give

$$ \mathcal{I}_{++}(E_1, E_2, K) + \mathcal{I}_{--}(-E_1, -E_2, K) = -\log[i(E_1 - K)\tau_0]\log[i(E_2 + K)\tau_0] \\ -\log[i(E_1 + K)\tau_0]\log[i(E_2 - K)\tau_0], \qquad (3.49) $$

which is exactly minus the combination $\mathcal{I}_{+-}(E_1, -E_2, K) + \mathcal{I}_{-+}(-E_1, E_2, K)$. Consequently, the graph integrals $\mathcal{I}_{ab}$ satisfy the largest-time equation.

## 4 Massive exchange correlator in de Sitter

We now turn to the case of the two-site chain of conformally coupled scalars mediated by the tree-level exchange of massive scalars in de Sitter. By explicitly performing the spectral integral, we derive a new closed-form expression for this correlator. As we will show, the branch cut in the energy domain due to particle production translates into a tower of poles in the complex mass domain. Summing over these residues leads to a new simple and partially resummed series solution.

We are interested in computing the de Sitter tree-level exchange graph<fnref n="15" />

$$ \begin{gathered} E_1 \xrightarrow{K} E_2 = \sum_{a,b=\pm} F_{ab}, \qquad (4.1) \\ \text{with } F_{ab} = -abK \int_{-\infty}^{0} \frac{d\tau'}{(-\tau')^{1/2}} e^{iaE_1\tau'} \int_{-\infty}^{0} \frac{d\tau''}{(-\tau'')^{1/2}} e^{ibE_2\tau''} \mathcal{G}_{ab}(K; \tau', \tau''). \end{gathered} $$

Recall that $\mathcal{G}_{ab}$ are the propagators for the canonically normalised field $\sigma(\tau, x) \equiv (-\tau)^{-3/2}\phi(\tau, x)$, hence the unusual power for the conformal time within the integral. We define the usual dimensionless kinematic variables

$$ u \equiv \frac{K}{E_1}, \quad v \equiv \frac{K}{E_2}, \qquad (4.2) $$

which lie in the unit disk, i.e. $|u|, |v| \le 1$, and the three-point function

$$ F^{(3)}(z, \mu) = \frac{1}{\sqrt{2}} |\Gamma(\frac{1}{2} + i\mu)|^2 P_{i\mu-1/2}(z), \qquad (4.3) $$

where $P_{i\mu-1/2}$ is the Legendre function defined in eq. (B.14). The non-time-ordered contribution, for which both time integrals factorise, is therefore given by

$$ F_{+-} = F^{(3)}(u^{-1}, \mu) F^{(3)}(v^{-1}, \mu). \qquad (4.4) $$

We have used the useful formula (B.16) to perform the time integrals. Notice that the Legendre $P$ function is real, i.e. $[P_{i\mu-1/2}(z)]^* = P_{i\mu-1/2}(z)$ for $z \ge 0$ and $\mu$ real. Therefore, the on-shell three-point function $F^{(3)}(z, \mu)$ is real. As a consequence, we obtain the same expression for $F_{-+}$. However, anticipating what follows, the off-shell function $F^{(3)}(z, \nu)$ where the mass parameter $\nu$ is analytically continued to the complex plane is not real.

$$ \langle \varphi_{k_1} \varphi_{k_2} \varphi_{k_3} \varphi_{k_4} \rangle' = \frac{g^2 \tau_0^4}{16k_1 k_2 k_3 k_4} F(E_1 = k_{12}, E_2 = k_{34}, K = s) + t + u \text{ channels,} $$

<footnote n="15">¹⁵Concretely, considering a conformally coupled scalar field $\varphi$ coupled to a massive scalar field $\phi$ by $\mathcal{L}_{\text{int}}/a^3(t) = -\frac{g^2}{2}\phi\varphi^2$ in de Sitter, the full correlator reads</footnote>

<footnote>where we have introduced a late-time cutoff $\tau_0$ and defined $F \equiv \sum_{a,b=\pm} F_{ab}$. Notice that we have defined the integrals $F_{ab}$ without the factor 2 to ease comparison with the literature.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Illustration of the analytic structure of the off-shell three-point function integrand in (4.5) in the complex $\nu$ plane, written in the Legendre $Q$ basis. The contour prescription avoids the infinite tower of poles lying on the imaginary axis and only picks up particle production poles.</figcaption></figure>

## 4.1 Off-shell three-point function

For pedagogical reasons, before computing the spectral integral for the time-ordered contribution $F_{++}$, we first study the following simpler integral

$$ \bar{F}^{(3)}(z) = \int_{-\infty}^{+\infty} d\nu \mathcal{N}_{\nu} \frac{F^{(3)}(z, \nu)}{(\nu^2 - \mu^2)_{i\epsilon}}, \qquad (4.5) $$

with $z \ge 1$, and $F^{(3)}$ defined in eq. (4.3). The object $\bar{F}^{(3)}$ is typically proportional to $F_{++}$ in the folded limit $E_2 \to K$ (or $E_1 \to K$ by symmetry), equivalently $v \to 1$ (or $u \to 1$), and already exhibits an interesting analytic structure. To be able to perform the spectral integral, it is essential to understand the analytic structure of the function $F^{(3)}(z, \nu)$ in the complex $\nu$ plane, as well as its behaviour at infinity, to appropriately determine how to close the contour.

The crucial observation is that the off-shell three-point function $F^{(3)}$ is a meromorphic function in the complex $\nu$ plane. At a fixed kinematic configuration $z$, the function $F^{(3)}(z, \nu)$ has poles at $\nu = \pm i(\frac{1}{2} + n)$ with integer $n \ge 0$, which corresponds to the poles of $|\Gamma(\frac{1}{2} + i\nu)|^2 = \Gamma(\frac{1}{2} + i\nu)\Gamma(\frac{1}{2} - i\nu)$, since the Legendre function $P$ is analytic in the entire complex $\nu$ plane. However, the Legendre $P$ representation of $F^{(3)}$ makes the large-$\nu$ asymptotic form not neat, i.e. it depends on the phase of $z$. In order to disentangle both behaviours at infinity, we project the off-shell three point function onto the Legendre $Q$ basis using the connection formula (B.17), yielding

$$ F^{(3)}(z, \nu) = \frac{e^{-\frac{i\pi}{2}}}{\sqrt{2} \sinh(\pi\nu)} \left[ Q_{-i\nu-1/2}(z) - Q_{+i\nu-1/2}(z) \right]. \qquad (4.6) $$

The function $Q_{\pm i\nu - 1/2}(z)$ has purely imaginary poles located at $\nu = \pm i(\frac{1}{2} + n)$ (which come from $\Gamma(\frac{1}{2} \pm i\mu)$ in eq. (B.14)), and should be closed in the lower (upper) complex plane since $Q_{\pm i\nu - 1/2}(z) \sim e^{\mp i\nu\xi}$ with $\xi > 0$. Therefore, when performing the spectral integral, we never
need to pick up these poles because the contour is closed to avoid the infinite tower.<fnref n="16" /> We illustrate the analytic structure of the integrand (4.5) and the contour prescription in figure 3. Eventually, picking up only the particle production poles, we obtain

$$ \bar{F}^{(3)}(z) = F^{(3)}(ze^{+i\pi}, \mu), \qquad (4.7) $$

where we have used the analytic continuation of the Legendre function (B.18) to re-express the result in the Legendre P basis. As already seen in section 2.3, the particle production pole structure projects ingoing modes on outgoing ones and vice versa. At the level of the three-point function, this amounts to rotate the kinematic configuration $z \to ze^{+i\pi}$. Actually, we show in appendix B.3 that the function $F^{(3)}(z, \mu)$ and its analytic continuation $F^{(3)}(ze^{+i\pi}, \mu)$, viewed as functions of $z$ for fixed $\mu$, are related by a dispersive integral in the complex energy domain. As we will now see, for the time-ordered contribution $F_{++}$, the infinite tower of poles cannot be evaded by the integration contour, and will result in series solution for the correlator.

## 4.2 Bootstrapping via the spectral representation

We now explicitly compute the integral

$$ F_{++} = - \int_{-\infty}^{+\infty} d\nu \mathcal{N}_{\nu} \frac{F^{(3)}(u^{-1}, \nu)F^{(3)}(v^{-1}, \nu)}{(\nu^2 - \mu^2)_{i\epsilon}}. \qquad (4.8) $$

This spectral representation makes it explicit that the time-ordered contribution $F_{++}$ is the spectral integral of the off-shell function $F_{+-}$ given in eq. (4.4). In practice, it means that solutions for the time-ordered contributions are of higher-transcendentality than the factorised ones. We expect that the result contains a local and non-local contributions which take a simple resummed form, as well as an EFT series contribution coming from infinite towers of poles. Proceeding the same way as previously for the three-point function, we decompose the functions $F^{(3)}$ in the Legendre Q basis to isolate well-behaved asymptotic forms at large $\nu$. The integrand takes the form

$$ \mathcal{N}_{\nu} F^{(3)}(u^{-1}, \nu) F^{(3)}(v^{-1}, \nu) = \frac{-\nu}{2\pi \sinh(\pi\nu)} \times [Q_{-i\nu-1/2}(u^{-1})Q_{-i\nu-1/2}(v^{-1}) - Q_{-i\nu-1/2}(u^{-1})Q_{+i\nu-1/2}(v^{-1}) + (\nu \leftrightarrow -\nu)] \qquad (4.9) $$

First, we immediately observe that both terms related by the shadow transformation $\nu \leftrightarrow -\nu$ contribute equally to the integral. Thus, it suffices to evaluate only one of these terms. Second, we notice that the term $Q_{-i\nu-1/2}(u^{-1})Q_{-i\nu-1/2}(v^{-1})$ has a well-defined large-$\nu$ behaviour that is independent of the hierarchy between the kinematic variables $u$ and $v$. However, the choice of whether to close the contour upwards of downwards for the term $Q_{-i\nu-1/2}(u^{-1})Q_{+i\nu-1/2}(v^{-1})$ depends on the relative magnitude between $u$ and $v$. As such, we will obtain a solution for $|u| < |v|$ and another one for $|u| > |v|$. The matching condition

<footnote n="16">¹⁶The factor $1/\sinh(\pi\nu) = |\Gamma(1 + i\nu)|^2/(\pi\nu)$ also brings additional poles to the integrand. However, they are precisely cancelled by the density of states $\mathcal{N}_{\nu}$.</footnote>
at $u = v$ is inherently satisfied, as it is encoded in the spectral integral. Therefore, we decompose the spectral integral into two pieces

$$F_{++} = F_{++}^{0} + F_{++}^{<, >}, \qquad (4.10)$$

with (we have accounted for the factor 2 coming from the $\nu \leftrightarrow -\nu$ symmetry)

$$F_{++}^{0} = \frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{d\nu \, \nu}{\sinh(\pi\nu)} \frac{Q_{-i\nu-1/2}(u^{-1})Q_{-i\nu-1/2}(v^{-1})}{(\nu^2 - \mu^2)_{i\epsilon}}, \qquad (4.11) \\ F_{++}^{<, >} = \frac{-1}{\pi} \int_{-\infty}^{+\infty} \frac{d\nu \, \nu}{\sinh(\pi\nu)} \frac{Q_{-i\nu-1/2}(u^{-1})Q_{+i\nu-1/2}(v^{-1})}{(\nu^2 - \mu^2)_{i\epsilon}}.$$

Despite the appearance, the contribution $F_{++}^{<, >}$ is well symmetric under $u \leftrightarrow v$, which can be explicitly checked after changing variables $\nu \leftrightarrow -\nu$. The analytic structure of these integrands reveals the presence of several poles with different origins:

* A first tower of poles comes from the factor $\nu/\sinh(\pi\nu) = \Gamma(1+i\nu)\Gamma(1-i\nu)/\pi$. They are located on the entire imaginary axis $\nu = \pm i(1+n)$ with $n \ge 0$, and cannot be evaded by the contour prescription.

* The factor $Q_{-i\nu-1/2}(u^{-1})Q_{-i\nu-1/2}(v^{-1})$ also generates an infinite tower of (second-order) poles located on the imaginary axis $\nu = -i(\frac{1}{2} + n)$ with $n \ge 0$. However, at large $\nu$, the integrand behaves as $\sim e^{+i\nu(\xi_u + \xi_v)}$ with $\xi_u = \cosh^{-1}(u^{-1}) \ge 0$ for $u^{-1} \ge 1$ (and similarly for $\xi_v$), so the contour should be closed in the upper-half complex plane. As a result, we do not collect these poles. This case is similar to the off-shell three-point function case seen in section 4.1. Here, both off-shell exchanged modes interfere constructively so that their spectrum spans only over half the imaginary axis, with residues doubling their weight.

* Poles of the term $Q_{-i\nu-1/2}(u^{-1})Q_{+i\nu-1/2}(v^{-1})$ are located at $\nu = \pm i(\frac{1}{2} + n)$ with $n \ge 0$. At large $\nu$, it behaves as $\sim e^{+i\nu(\xi_u - \xi_v)}$. For $|u| < |v|$, we have $\xi_u - \xi_v \ge 0$ and we close the contour in the upper-half complex plane, picking the tower of (first-order) poles. This case reflects the destructive interference between off-shell exchanged modes, resulting in the spectrum spanning over the entire imaginary axis.

* Eventually, we also have the usual particle production poles coming from the $i\epsilon$ prescription.

Let us now collect these various poles in turn.

**Particle production residues.** The poles coming from the $i\epsilon$ prescription are the easiest to evaluate. We obtain

$$F_{++}^{0} \supset \frac{-i}{2\sinh^2(\pi\mu)} \left[ e^{+\pi\mu} Q_{+i\mu-1/2}(u^{-1}) Q_{+i\mu-1/2}(v^{-1}) + (\mu \leftrightarrow -\mu) \right]. \qquad (4.12)$$

For the contribution $F_{++}^{<}$ for which $|u| < |v|$, collecting both poles on the upper-half complex plane results in

$$F_{++}^{<} \supset \frac{i}{2\sinh^2(\pi\mu)} \left[ e^{+\pi\mu} Q_{+i\mu-1/2}(u^{-1}) Q_{-i\mu-1/2}(v^{-1}) + (\mu \leftrightarrow -\mu) \right]. \qquad (4.13)$$
Combining both contributions, and projecting back onto the Legendre P basis using (B.17) and (B.18) yields

$$F_{++}^{\text{collider}} = -F^{(3)}(-u^{-1}, \mu)F^{(3)}(v^{-1}, \mu), \quad (4.14)$$

for $|u| < |v|$. The contribution for which $|u| > |v|$ is simply found after swapping $u$ and $v$. The found factorised form is not surprising as in the late-time limit time ordering vanishes, as seen in section 2.3. The result (4.14) can be directly found by performing the time integrals after substituting $G_{++} \rightarrow G_{-+}$.

**EFT residues.** Contributions coming from EFT residues are analytic in both $u$ and $v$ as $u, v \rightarrow 0$, and therefore admits series representations. We first sum over the residues from the factor $1/\sinh(\pi\nu)$. The integration contour picks up only half of this tower of first-order poles on the imaginary axis. Explicitly, we obtain

$$F_{++}^{0} \supset -\frac{2}{\pi} \sum_{n=0}^{+\infty} (-1)^n \frac{(n+1)}{(n+1)^2 + \mu^2} Q_{n+1/2}(u^{-1})Q_{n+1/2}(v^{-1}), \quad (4.15)$$

where we have used that the residue of $\Gamma(1 + i\nu)$ at $\nu = +i(1 + n)$ is $\frac{(-1)^n}{n!}$. The terms in $1/\mu^2$, characteristic of an EFT expansion, come from

$$\frac{1}{(\nu^2 - \mu^2)_{i\epsilon}} \underset{\nu \to +i(n+1)}{\longrightarrow} \frac{-1}{(n+1)^2 - \mu^2}. \quad (4.16)$$

The large $n$ behaviour of the series coefficients is given by

$$\left| \frac{(n+1)}{(n+1)^2 + \mu^2} Q_{n+1/2}(u^{-1}) Q_{n+1/2}(v^{-1}) \right| \sim \frac{1}{(n+1)^2 + \mu^2} \frac{e^{-(n+1)(\xi_u + \xi_v)}}{\sqrt{\sinh \xi_u \sinh \xi_v}}, \quad (4.17)$$

which is well convergent. Notice that close to $u, v = 1$, the exponential damping becomes less efficient and the series converges slower, i.e. as $1/n^2$. Similarly, the contribution of these residues to $F_{++}^<$ is found to be

$$F_{++}^{<} \supset \frac{2}{\pi} \sum_{n=0}^{+\infty} (-1)^n \frac{(n+1)}{(n+1)^2 + \mu^2} Q_{n+1/2}(u^{-1}) Q_{-n-3/2}(v^{-1}). \quad (4.18)$$

Quite interestingly, the two towers of residues (4.15) and (4.18) exactly cancel each other. Indeed, both contributions can be combined in a single series, and using (B.17), each series coefficient is proportional to

$$Q_{-n-3/2}(v^{-1}) - Q_{n+1/2}(v^{-1}) \propto \sin[\pi(n+1)]. \quad (4.19)$$

This effect, which generally occurs in odd spatial dimensions, was previously noted in [18] and is applicable to contact diagrams as well. Notably, while this phenomenon was initially understood as interference between the two branches of the in-in contour, we show here that this cancellation happens within a single branch. We finish with the poles located at $\nu = \pm i(\frac{1}{2} + n)$ which only enter the contribution $F_{++}^{<,>}$. They are the only ones contributing to the final expression. For $|u| < |v|$, we obtain the rather simple expression

$$F_{++}^{\text{EFT}} = u \sum_{n=0}^{+\infty} \frac{(-1)^n}{(n + \frac{1}{2})^2 + \mu^2} \left(\frac{u}{v}\right)^n {}_2F_1\left[\frac{1+n}{2} + \frac{1-n}{2}; u^2\right] {}_2F_1\left[\frac{1-n}{2} - \frac{n}{2}; v^2\right]. \quad (4.20)$$

The Legendre $Q$ function being not defined for negative integer order, the above expression is found after evaluating the residue of the $\Gamma$ function entering (B.14).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The four-point correlation function $F(u, v) = \sum_{a,b=\pm} F_{ab}$, as defined in eq. (4.1), as function of the dimensionless kinematic variable $u$ with fixed $v = 0.5$, for $\mu = 1$ (left panel) and $\mu = 2$ (right panel). We include terms in the series from $n = 1, \dots, N$. These results reproduce those presented in [12] (figure 6) and in [15] (figure 6). The Mathematica notebook with the analytical expression for $F$ is available on the Github repository [https://github.com/deniswerth/Spectral-Representation].</figcaption></figure>

**Full result.** The full result is found by summing all contributions, namely the non-time ordered pieces (4.4) and the time-ordered pieces combining the collider signal (4.14) with the EFT series (4.20). A few comments are in order. First, this method yields a convergent result in all kinematic configurations. Although not obvious from the analytical expression, we have checked numerically that the found result perfectly matches the one originally found in [12]. Figure 4 shows the full closed-form result when keeping only a few terms in the series. We observe that the series convergence rate is the same as the bootstrap result but slower than the series solution found using the partial Mellin-Barnes method [15]. Second, compared to the bootstrap result of [12], the spectral integration naturally resums one of the two series, which in practice makes the evaluation much faster. Lastly, since the spectral representation is solution to the bootstrap equation, there is a one-to-one correspondence between the basis of functions onto which the homogeneous solution of the bootstrap equation is expanded and the Legendre $P$ or $Q$ basis. Additionally, when solving the boundary differential equation, fixing the free coefficients requires two boundary conditions for specific values (or limits) of $u$ and $v$. Here, these coefficients are fully encoded in the spectral integral. Eventually, the spectral method does not requires matching different solutions at $u = v$ because the transition of closing the contour, whether upwards or downwards, is continuous. However, the first derivative of the result is discontinuous at the transition, resulting in a noticeable kink. In short, rather than integrating the boundary differential equation, we directly perform the spectral integral, with the bounds inherently determining the free integration constants.

**Largest-time equation.** Assuming $|u| < |v|$, it is trivial to observe that

$$F_{++}(u, v) + F_{--}(-u, -v) + F_{+-}(-u, -v) + F_{-+}(-u, v) = 0, \quad (4.21)$$

from their explicit forms in (4.4), (4.14) and (4.20). Indeed, both EFT series cancel against each other thanks to the overall factor $u$, i.e. $F^{\text{EFT}}(-u, -v) = -F^{\text{EFT}}(u, v)$. Similarly,
$F_{++}^{\text{collider}}(u, v)$ cancels against $F_{+-}(-u, v)$, and $F_{--}^{\text{collider}}(-u, -v)$ against $F_{+-}(u, -v)$. Of course, swapping $u$ and $v$ yields the same result. This provides an additional consistency check for the derived formula.

## 4.3 Singularities and analytic structure

The spectral representation of the four-point function enables straightforward examination of various limits by directly inspecting the spectral integral. A first observation is that removing the dynamical part of the propagating exchanged field — that is, discarding the $i\epsilon$ prescription encoding particle production—, the spectral integral reduces to a contact four-point function generated by the leading bulk interaction $\phi^4$ in a gradient expansion. This simplification is revealed through the generalised Mehler-Fock transformation (eq. (14.20.14) of [96])

$$ \int_{-\infty}^{+\infty} d\nu \mathcal{N}_{\nu} F^{(3)}(u^{-1}, \nu) F^{(3)}(v^{-1}, \nu) = \frac{uv}{u+v}. \qquad (4.22) $$

This identity is analogous to setting the differential operator of the bootstrap equation to unity and therefore recovering the usual source term. Let us now probe the following limits: (i) collapsed limit $u, v \to 0$, i.e. internal soft momentum $K \to 0$ while keeping all external momenta hard, and (ii) the partial-energy pole $u \to -1$ while keeping $v$ fixed.

** Collapsed limit.** In the limit $u, v \to 0$, we use the asymptotic behaviour of the Legendre $P$ function (B.20) for large argument to easily obtain

$$ \lim_{u,v \to 0} F_{+-} = \left(\frac{uv}{4}\right)^{\frac{1}{2}+i\mu} \frac{\Gamma(\frac{1}{2}+i\mu)^2 \Gamma(-i\mu)^2}{2\pi}. \qquad (4.23) $$

Taking the limit inside the spectral integral, we obtain

$$ \lim_{u,v \to 0} F_{++} = -\frac{1}{2\pi} \left(\frac{uv}{4}\right)^{\frac{1}{2}} \int_{-\infty}^{+\infty} \frac{d\nu}{(\nu^2 - \mu^2)_{i\epsilon}} \left(\frac{uv}{4}\right)^{+i\nu} \Gamma(\frac{1}{2} + i\nu)^2 \frac{\Gamma(-i\nu)}{\Gamma(+i\nu)}. \qquad (4.24) $$

At large $\nu$, the integrand scales as $\sim e^{i\nu \log(uv/4)}$ (up to a phase) so that the contour is closed in the lower-half complex plane (recall $\log(uv/4) \le 0$). The usual particle production poles give

$$ \lim_{u,v \to 0} F_{++} \supset \frac{i}{2} \left[ e^{+\pi\mu} \left(\frac{uv}{4}\right)^{\frac{1}{2}+i\mu} \frac{\Gamma(\frac{1}{2}+i\mu)^2 \Gamma(-i\mu)^2}{2\pi} + (\mu \leftrightarrow -\mu) \right]. \qquad (4.25) $$

Notice that two infinite towers of poles arise: one from $\Gamma(\frac{1}{2} + i\nu)$ located at $\nu = +i(\frac{1}{2} + n)$ that we do not collect as they lie in the upper-half complex plane, and another from $\Gamma(-i\nu)$ located at $\nu = -in$ with $n \ge 0$ integer, which we do collect. By definition, this EFT signal is analytic in both $u$ and $v$ as it admits a series representation. Its $u, v \to 0$ limit can be directly recovered from its full expression (4.20). Eventually, keeping only the terms which are non-analytic in $uv$ as $u, v \to 0$ (and therefore also non-analytic in $K$), we obtain

$$ \lim_{u,v \to 0} F = \left(\frac{uv}{4}\right)^{\frac{1}{2}+i\mu} (1 + i \sinh \pi\mu) \frac{\Gamma(\frac{1}{2}+i\mu)^2 \Gamma(-i\mu)^2}{2\pi} + \text{c.c.}, \qquad (4.26) $$

which exactly reproduces the expression found in [11].
**Partial-energy pole.** At fixed $v$, let us probe the limit $u \to -1$. The three-point function $F^{(3)}(z, \mu)$ has a branch point at $z = -1$ so we analyse its behaviour directly at the level of the time integral. The limit $E_1 + K \to 0$ probes the early time-limit of the left vertex, so we can take the early-time limit of the Hankel function inside the integral

$$F^{(3)}(u^{-1}, \mu) \sim \frac{1}{\sqrt{2}} \lim_{z \to 0} \int_z^{+\infty} \frac{dx}{x} e^{-ix(1+u^{-1})} = \frac{1}{\sqrt{2}} \lim_{z \to 0} E_1 \left[ (1+u^{-1}) iz \right], \quad (4.27)$$

where $E_1$ is the exponential integral, which has a logarithmic branch. Isolating the leading term, we obtain

$$\lim_{u \to -1} F^{(3)}(u^{-1}, \mu) = -\frac{1}{\sqrt{2}} \log(1+u). \quad (4.28)$$

The non-time-ordered contributions give

$$\lim_{u \to -1} (F_{+-} + F_{-+}) = -\sqrt{2} \log(1+u) \times F^{(3)}(v^{-1}, \mu). \quad (4.29)$$

Similarly, the time-ordered contribution is found to be

$$\lim_{u \to -1} F_{++} = \frac{1}{\sqrt{2}} \log(1+u) \int_{-\infty}^{+\infty} d\nu N_\nu \frac{F^{(3)}(v^{-1}, \nu)}{(\nu^2 - \mu^2)_{i\epsilon}} = \frac{1}{\sqrt{2}} \log(1+u) \times F^{(3)}(-v^{-1}, \mu), \quad (4.30)$$

where we have used (4.7). Summing all contributions, we eventually obtain

$$\lim_{u \to -1} F = -\sqrt{2} \log(1+u) \times [F^{(3)}(v^{-1}, \mu) - F^{(3)}(-v^{-1}, \mu)]. \quad (4.31)$$

In this limit, the residue of the partial-energy pole is proportional to the shifted three-point function. This behaviour generalises to any exchanged process where the sum of all “energies” entering a vertex vanishes, see e.g. [46–48]. Similarly, the correlator has a singularity in the flat-space limit $u \to -v$, i.e. $E = E_1 + E_2 \to 0$ or equivalently $\mu \to \infty$ as seen in section 2.3, with a coefficient that is related to the flat-space scattering amplitude [37, 38].

## 5 Conclusions

In this paper, we have explored an off-shell approach to cosmological correlators. For massive field experiencing particle production, going to the complex mass plane allows one to write a dispersive integral for the propagator that naturally encodes time ordering. Using this object, we have shown that exchanged correlators can be obtained by gluing off-shell lower-point correlators. In this procedure, performing the resulting spectral integral, which amounts to collecting poles as the integrand is meromorphic, effectively sets the exchanged particle on shell. We argued that this representation not only clarifies the analytic structure and factorisation properties of cosmological correlators but also simplifies explicit calculations. As a specific example, we derived a new, simple closed-form formula for the four-point correlator of conformally coupled fields exchanging a massive field in de Sitter space and examined its analytic properties. Additionally, we derived cosmological largest-time equations, which relate individual correlator in-in branch channels through the analytic continuation of certain external energies. These relations, explicitly illustrated for simple diagrams, can serve as consistency checks for more complex correlators.
Eventually, our work opens several interesting avenues for future investigation:

* First of all, the spectral representation of massive propagators used in this work does not crucially rely on de Sitter isometries. Meanwhile, a large level of non-Gaussianity can be reached in cosmological processes that strongly break de Sitter boosts. It would therefore be natural to use this spectral representation to obtain clean and simpler closed-form solutions for correlators of particles featuring reduced sound speeds or in the presence of a chemical potential. The latter case is known to boost the particle production rate, which can lead to enhanced cosmological collider signals. However, since mode functions are described by Whittaker functions, the spectral representation should be upgraded accordingly.

* In the context of modern scattering amplitude techniques, complex amplitudes can be efficiently constructed by sewing together lower-point amplitudes using sophisticated recursion relations. We have shown a glimpse of how such procedure can be applied to cosmological correlators in simple cases. Since time integrals are trivialised, computing higher-point correlation functions primarily involves performing a series of spectral integrals, which effectively reduces to summing over several towers of poles, given that off-shell lower-point correlators are meromorphic in the complex mass plane. It will be interesting to explore how this off-shell approach can render computations of more complex correlation functions more tractable. Ultimately, we believe this approach lays the groundwork for deriving more general cosmological recursion relations beyond rational correlators.

* Finally, as in flat space, states in a unitary quantum field theory in de Sitter are classified by unitary irreducible representations of the isometries, which are parametrised by the conformal dimension (related to the mass) and spin. Imposing unitarity as an additional consistency condition imposes a set of positivity constraints. While previous such bounds have been largely theoretical, exploring their phenomenological implications could yield valuable insights. The spectral representation employed in this work can be viewed as the leading-order perturbative Källén-Lehmann representation in spatial Fourier space. By perturbatively computing massive self-energy corrections, it may be possible to extend this approach further in perturbation theory, potentially deriving useful bounds from the positivity and meromorphic properties of the spectral density. Ultimately, one would hope to use this spectral representation to “dress” massive propagators in de Sitter space by resumming bubble diagrams.

## Acknowledgments

We thank Guillaume Faye, Jean-Baptiste Fouvry, Mang Hei Gordon Lee, Scott Melville, Enrico Pajer, Sébastien Renaux-Petel, Xi Tong and Zhong-Zhi Xianyu for helpful discussions. We also thank Lucas Pinol, Arthur Poisson and Sébastien Renaux-Petel for comments on the draft. DW is supported by the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement No 758792, Starting Grant project GEODESI).
## A Feynman $i\epsilon$ prescription from vacuum wave-functional

In this appendix, we come back to the $i\epsilon$ prescription defining the in-in generating function (2.2). As is well known, deforming the time integration contour around the infinite past explicitly breaks unitarity [97–100]. This is a consequence of both in-in branches not being invariant under time reversal, producing a commutator of time evolution operators and not identity. This exact same prescription, when thought of as analytically continuing energy instead of time, not only preserves unitarity, but leads to the correct $i\epsilon$ prescription in the Feynman propagator (2.9).

Reaching convergence in the infinite past is reminiscent of adiabatically projecting the vacuum of the fully interacting theory $|\Omega\rangle$ onto the free one $|0\rangle$. Within the in-in path integral (2.2), this procedure is achieved with the additional term $\Psi[\varphi_+](t_0)\Psi^*[\varphi_-](t_0)$ where $\Psi[\varphi](t) \equiv \langle\varphi(t)|\Omega\rangle$ is the vacuum wave-functional that describes the transition amplitude from the vacuum $|\Omega\rangle$ to a specific field configuration $\varphi(x)$ at some time $t_0$ [71, 81, 101, 102]. Let us therefore compute this object and take the limit $t_0 \to -\infty$. The key insight is to use the defining property of the annihilation operator $\hat{a}_k |\Omega\rangle = 0$. After inverting the canonical quantisation of the field $\hat{\varphi}_k$ and its conjugate momentum $\hat{p}_k$, and using the Wronskian condition $u_k p_k^* - u_k^* p_k = i$, valid at all time and fixed by the commutation relation $(p_k)$ is the mode function associated to $\hat{p}_k$), one obtains

$$ \hat{a}_k = \frac{1}{i} (p_k^* \hat{\varphi}_k - u_k^* \hat{p}_k). \qquad (A.1) $$

Using the field-space representation of the momentum operator $\hat{p}_k \to -i\delta/\delta\varphi_k$, with $\delta\varphi_q/\delta\varphi_k = (2\pi)^3\delta^{(3)}(q+k)$, the condition $\hat{a}_k |\Omega\rangle = 0$ projected onto $\langle\varphi(t_0)|$ yields the functional differential equation

$$ \left( u_k^* \frac{\delta}{\delta \varphi_k} - i p_k^* \varphi_k \right) \Psi[\varphi](t_0) = 0. \qquad (A.2) $$

Solving this equation gives the following Gaussian solution for the vacuum wave-functional

$$ \Psi[\varphi](t_0) = \exp \left( -\frac{1}{2} \int \frac{d^3 k}{(2\pi)^3} \varphi_k(t_0) \omega_k(t_0) \varphi_{-k}(t_0) \right), \qquad (A.3) $$

where $\omega_k(t_0) = -ip_k^*/u_k^*$, with the mode functions evaluated at $t_0$. The overall normalisation is fixed to unity due to the correct vacuum normalisation $\langle\Omega|\Omega\rangle = 1$. We now take the limit $t_0 \to -\infty$. Reaching the Bunch-Davies state imposes that mode functions behave as simple plane-waves $u_k(t) \sim e^{-ikt}$ so that $\omega_k = k$ is nothing but the dispersion relation. As for the remaining terms, we introduce the following regularisation scheme

$$ f(-\infty) = \lim_{\epsilon \to 0} \epsilon \int_{-\infty}^{t} dt' f(t') e^{\epsilon t'}, \qquad (A.4) $$

valid for any smooth-enough function $f(t)$, that can be recovered after integration by parts, to rewrite the contribution from vacuum wave-functionals entering (2.2) as

$$ \lim_{t_0 \to -\infty} \Psi[\varphi_+](t_0) \Psi^*[\varphi_-](t_0) = \lim_{\epsilon \to 0} \exp \left( -\frac{\epsilon}{2} \int_{-\infty}^{t} dt' \int \frac{d^3 k}{(2\pi)^3} k [\varphi_{k,+}(t')\varphi_{-k,+}(t') + (+\leftrightarrow -)] \right). \qquad (A.5) $$
We have set $e^{\epsilon t'} = 1$ since we only care about the leading term as $\epsilon \to 0$. Obtaining a Gaussian weight in the fields with negative real values, as opposed to purely imaginary ones in the standard action, turns out to be crucial to ensure convergence of the path integral in the infinite past. Indeed, let us consider the free quadratic Lagrangian for the fields $\varphi_{\pm}$ written in Fourier space $\mathcal{L}[\varphi_{k,\pm}] = -\frac{1}{2}\varphi_{k,\pm}(t)(-\partial_t^2 - k^2)\varphi_{-k,\pm}(t)$. Adding the term (A.5) amounts to modifying the weight in the exponential of the generating function (2.2) to

$$ \pm \frac{i}{2} \int \frac{d^3 k}{(2\pi)^3} \int_{-\infty}^{t} dt' \varphi_{k,\pm}(t') [-\partial_{t'}^2 - (k \pm i\epsilon)^2] \varphi_{-k,\pm}(t'), \quad (A.6) $$

where we have used $(k \pm i\epsilon)^2 \approx k^2 \pm 2i\epsilon$ and relabelled $2\epsilon \to \epsilon$. In the end, performing the converging Gaussian path integral over the fields $\varphi_{\pm}$ leads to the correct $i\epsilon$ prescription for the Feynman propagator (2.9).

## B Details on massive fields in de Sitter

Consider a real scalar field $\phi(\tau, x)$ of mass $m$, where $-\infty < \tau \le 0$ is the conformal time, evolving in de Sitter space with three spatial dimensions. We will focus on fields in the principal series so that $\mu \equiv \sqrt{m^2/H^2 - 9/4}$ is real. In what follows, we set $H = 1$ for convenience. Defining the canonically normalised field $\sigma(\tau, x) \equiv (-\tau)^{-3/2}\phi(\tau, x)$, its positive-frequency mode function is given by

$$ u_k(\tau, \mu) = \frac{i\sqrt{\pi}}{2} e^{-\frac{\pi\mu}{2}} H_{i\mu}^{(1)}(-k\tau) \underset{|k\tau|\to\infty}{\longrightarrow} \frac{e^{-ik\tau}}{\sqrt{2k\tau}}, \quad (B.1) $$

where $H_{i\mu}^{(1)}(z)$ is the Hankel function of the first kind, such that it reduces to the usual Bunch-Davies vacuum in the far past, up to an overall unimportant phase. The in-in contour $C_{i\epsilon}$ in (2.2) requires analytically continuing the positive-frequency mode function (B.1) to the lower-half $k\tau$-complex plane to ensure reaching the asymptotic vacuum. This precisely avoids crossing the branch cut of the Hankel function $H_{i\mu}^{(1)}(z)$, that is chosen to lie on the real negative axis, away from the physical regime. Therefore, using the complex conjugate relation (B.3), the negative-frequency mode can be written

$$ u_k^*(\tau, \mu) = -\frac{i\sqrt{\pi}}{2} e^{+\frac{\pi\mu}{2}} H_{i\mu}^{(2)}(-k\tau). \quad (B.2) $$

### B.1 Useful formulae

The following mathematical identities are useful for manipulating mode functions of a massive scalar field in de Sitter.

**Complex conjugate relations.** The Hankel functions of the first and second kind are related by the complex conjugate relations

$$ [H_{i\mu}^{(1)}(z)]^* = e^{+\pi\mu} H_{i\mu}^{(2)}(z^*), \quad [H_{i\mu}^{(2)}(z)]^* = e^{-\pi\mu} H_{i\mu}^{(1)}(z^*), \quad (B.3) $$

for $z \in \mathbb{C} \setminus \{\mathbb{R}^-\}$. Notice that these relations are not verified when the argument $z$ lies on the branch cut.
**Symmetries.** The Hankel functions obey the following relations

$$H_{-i\mu}^{(1)}(z) = e^{-\pi\mu} H_{i\mu}^{(1)}(z), \quad H_{-i\mu}^{(2)}(z) = e^{+\pi\mu} H_{i\mu}^{(2)}(z), \qquad (B.4)$$

for $z \in \mathbb{C}$. These identities render manifest an underlying shadow $\mu \leftrightarrow -\mu$ symmetry for the positive- and negative-frequency mode functions

$$u_k(\tau, -\mu) = u_k(\tau, \mu), \quad u_k^*(\tau, -\mu) = u_k^*(\tau, \mu). \qquad (B.5)$$

Upon rotating time, the Hankel functions satisfy

$$H_{-i\mu}^{(1)}(e^{+i\pi}z) = -H_{i\mu}^{(2)}(z), \quad H_{-i\mu}^{(2)}(e^{-i\pi}z) = -H_{i\mu}^{(1)}(z). \qquad (B.6)$$

The first relation is valid for $z$ in the lower-half complex plane excluding the negative real and imaginary axis, and the second relation is valid for $z$ in the upper-half complex plane excluding the positive real and imaginary axis. In the physical time domain, these identities relate the mode functions by a CPT symmetry

$$u_k(z, \mu) = u_k^*(e^{-i\pi}z, -\mu), \quad u_k^*(z, \mu) = u_k(e^{+i\pi}z, -\mu). \qquad (B.7)$$

Additional analytic continuations of the Hankel function are given by

$$\begin{aligned}
H_{i\mu}^{(1)}(e^{-i\pi}z) &= e^{+\pi\mu}H_{i\mu}^{(2)}(z) + 2\cosh(\pi\mu)H_{i\mu}^{(1)}(z), \\
H_{i\mu}^{(2)}(e^{+i\pi}z) &= e^{-\pi\mu}H_{i\mu}^{(1)}(z) + 2\cosh(\pi\mu)H_{i\mu}^{(2)}(z).
\end{aligned} \qquad (B.8)$$

**Connection formulae.** Hankel functions can be expanded in terms of Bessel functions with the following connection formulae

$$H_{i\mu}^{(1)}(z) = \frac{e^{+\pi\mu}J_{+i\mu}(z) - J_{-i\mu}(z)}{\sinh(\pi\mu)}, \quad H_{i\mu}^{(2)}(z) = \frac{J_{-i\mu}(z) - e^{-\pi\mu}J_{+i\mu}(z)}{\sinh(\pi\mu)}. \qquad (B.9)$$

The Bessel function of the first kind satisfies

$$J_{i\mu}^*(z) = J_{-i\mu}(z^*), \qquad (B.10)$$

for $z \in \mathbb{C} \setminus \{\mathbb{R}^{-}\}$, and has the following power-law asymptotic form for small argument $z \to 0$

$$J_{i\mu}(z) \sim \frac{1}{\Gamma(1 + i\mu)} \left(\frac{z}{2}\right)^{i\mu}, \qquad (B.11)$$

whereas the Hankel function has the following plane-wave asymptotic form for large argument $z \to \infty$

$$H_{i\mu}^{(1)}(z) \sim \sqrt{\frac{2}{\pi z}} e^{i(z-\frac{\pi}{4})+\frac{\pi\mu}{2}}. \qquad (B.12)$$

**Large-order asymptotic forms.** For large $\mu$ with $z \neq 0$ fixed, we have

$$J_{i\mu}(z) \sim \frac{1}{\sqrt{2i\pi\mu}} \left(\frac{ez}{2i\mu}\right)^{i\mu}, \quad H_{i\mu}^{(1)}(z) \sim -H_{i\mu}^{(2)}(z) \sim -i\sqrt{\frac{2}{i\pi\mu}} \left(\frac{ez}{2i\mu}\right)^{-i\mu}. \qquad (B.13)$$
**Legendre functions.** The Legendre functions are equivalent to Bessel (and Hankel) functions with higher transcendentality. Explicitly, they are given by

$$ P_{i\mu-1/2}(z) \equiv {}_2F_1\left[\begin{array}{c}\frac{1}{2}+i\mu, \frac{1}{2}-i\mu, \frac{1-z}{2} \\ 1\end{array}\right], \qquad (B.14) \\ Q_{i\mu-1/2}(z) \equiv \sqrt{\frac{\pi}{2}} \frac{\Gamma(\frac{1}{2}+i\mu)}{\Gamma(1+i\mu)} 2^{-i\mu}z^{-\frac{1}{2}-i\mu} {}_2F_1\left[\begin{array}{c}\frac{3/2+i\mu}{2}, \frac{1/2+i\mu}{2} \\ 1+i\mu \end{array}\middle|\frac{1}{z^2}\right], $$

in terms of the hypergeometric function ${}_2F_1$ that is defined by the following formal series

$$ {}_2F_1\left[\begin{array}{c}a \ b \\ c \end{array}; z\right] = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n}{(c)_n} \frac{z^n}{n!}, \qquad (B.15) $$

on the disk $|z| < 1$ and by analytic continuation elsewhere, where $(a)_n \equiv \Gamma(a+n)/\Gamma(a)$ is the Pochhammer symbol. These functions are $P_\alpha(z) = \text{LegendreP}[\alpha, 0, 3, z]$ and $Q_\alpha(z) = \text{LegendreQ}[\alpha, 0, 3, z]$ in Mathematica. The optional argument “3” selects the correct branch cut structure as it is falsely implemented by default.<fnref n="17" /> At fixed argument, the Legendre $P$ function is analytic in the entire complex $\mu$ plane. However, the Legendre $Q$ function is meromorphic as it has poles coming from the factor $\Gamma(\frac{1}{2} + i\mu)$.

**Integral formula.** The Legendre function $P$ is the solution of the following integrals

$$ -i \frac{\sqrt{\pi}}{2} e^{+\frac{\pi \nu}{2}} \int_{-\infty^+}^0 \frac{d\tau e^{ik\tau}}{(-\tau)^{1/2}} H_{i\nu}^{(2)}(-s\tau) = \frac{e^{-\frac{i\pi}{4}}}{\sqrt{2s}} |\Gamma(\frac{1}{2} + i\nu)|^2 P_{i\nu-1/2}(k/s), \qquad (B.16) \\ +i \frac{\sqrt{\pi}}{2} e^{-\frac{\pi \nu}{2}} \int_{-\infty^-}^0 \frac{d\tau e^{-ik\tau}}{(-\tau)^{1/2}} H_{i\nu}^{(1)}(-s\tau) = \frac{e^{+\frac{i\pi}{4}}}{\sqrt{2s}} |\Gamma(\frac{1}{2} + i\nu)|^2 P_{i\nu-1/2}(k/s). $$

**Connection formula.** Analogously to Bessel and Hankel functions, the Legendre functions satisfy the connection formula

$$ |\Gamma(\frac{1}{2} + i\mu)|^2 P_{i\mu-1/2}(z) = \frac{e^{-\frac{i\pi}{2}}}{\sinh(\pi\mu)} \left[Q_{-i\mu-1/2}(z) - Q_{+i\mu-1/2}(z)\right]. \qquad (B.17) $$

Its analytic continuation is defined by

$$ |\Gamma(\frac{1}{2} + i\mu)|^2 P_{i\mu-1/2}(ze^{+i\pi}) = \frac{1}{\sinh(\pi\mu)} \left[ e^{+\pi\mu} Q_{+i\mu-1/2}(z) - e^{-\pi\mu} Q_{-i\mu-1/2}(z) \right]. \qquad (B.18) $$

**Large-order asymptotic forms.** The Legendre $P$ function has a simple behaviour at $z \to 0$ (similar to Bessel $J$ function), while the Legendre $Q$ function has a simple behaviour at $z \to \infty$ (similar to Hankel $H$ function). Indeed, they are given by

$$ P_{i\nu-1/2}(\cosh \xi) \sim \left( \frac{\xi}{\sinh \xi} \right)^{1/2}, \qquad (B.19) \\ Q_{i\nu-1/2}(\cosh \xi) \sim \frac{\pi}{2i} \left( \frac{2\xi}{\pi\nu\xi\sinh\xi} \right)^{1/2} e^{-i(\nu\xi - \frac{\pi}{4})}. $$

<footnote n="17">¹⁷The Legendre functions $P$ and $Q$ generalise to associated Legendre functions (also called Ferrers function of the first and second kind) $P_{i\mu-1/2}(z) \to P_{i\mu-1/2}^{1/2-j}(z)$ and $Q_{i\mu-1/2}(z) \to Q_{i\mu-1/2}^{1/2-j}(z)$, see Chap. 14 of [96], that are useful when dealing with spatial derivative interactions, i.e. for different powers of conformal time in the correlator time integrals. The non-derivative interaction case simply corresponds to $j = 1/2$.</footnote>
These expressions equivalently give the large-order asymptotic forms. The asymptotic behaviour of the Legendre P function reads

$$ P_{i\nu-1/2}(z) \sim \frac{1}{\sqrt{\pi}} \frac{\Gamma(-i\nu)}{\Gamma(\frac{1}{2}-i\nu)} \left(\frac{1}{2z}\right)^{\frac{1}{2}+i\nu}, \qquad (B.20) $$

for $z \to \infty$.

## B.2 Derivation of the spectral representation

For completeness, we present the proof of the spectral representation of massive scalar field in de Sitter, closely following [63]. Let us first consider the case where the time $\tau_1$ is in the future of $\tau_2$, i.e. $-k\tau_1 < -k\tau_2$, so that the time-ordered propagator reduces to $G_{++}(k; \tau_1, \tau_2) = u_k(\tau_1, \mu)u_k^*(\tau_2, \mu)$. Expanding the first Hankel function $H_{i\nu}^{(2)}(-k\tau_1)$ in (2.15) in terms of Bessel functions using (B.9), using eq. (B.4), and the symmetry of the pole prescription with respect to $\nu \leftrightarrow -\nu$, the integral can be written

$$ G_{++}(k; \tau_1, \tau_2) = \frac{i}{2} \int_{-\infty}^{+\infty} d\nu \frac{\nu J_{i\nu}(-k\tau_1) H_{i\nu}^{(2)}(-k\tau_2)}{(\nu^2 - \mu^2)_{i\epsilon}}. \qquad (B.21) $$

The crucial stage is to observe that the numerator $\nu J_{i\nu}(-k\tau_1) H_{i\nu}^{(2)}(-k\tau_2)$ is analytic in the entire complex $\nu$ plane, so that we choose to close the integration contour in the lower-half complex plane [103]. Since the large-$\nu$ expansion (B.13) gives

$$ \nu J_{i\nu}(z_1) H_{i\nu}^{(2)}(z_2) \sim \frac{e^{i\nu \log(z_1/z_2)}}{\pi}, \qquad (B.22) $$

the arc at infinity does not contribute for $z_1 < z_2$ and $\mathrm{Im}(\nu) \to -\infty$. By Cauchy's residue theorem, only the two poles located at $\nu = \pm\mu - i\epsilon$ are selected by the closed contour integral, which yields

$$ \begin{aligned} G_{++}(k; \tau_1, \tau_2) &= \frac{\pi}{4 \sinh(\pi \mu)} \left( e^{+\pi \mu} J_{+i\mu}(-k\tau_1) H_{+i\mu}^{(2)}(-k\tau_2) - e^{-\pi \mu} J_{-i\mu}(-k\tau_1) H_{-i\mu}^{(2)}(-k\tau_2) \right) \\ &= u_k(\tau_1, \mu) u_k^*(\tau_2, \mu), \end{aligned} \qquad (B.23) $$

where we have again used eqs. (B.4) and (B.9). For the case where the time $\tau_2$ is in the future of $\tau_1$, i.e. $-k\tau_2 < -k\tau_1$, one needs to expand the second Hankel function $H_{i\nu}^{(2)}(-k\tau_2)$ in (2.15) in terms of Bessel functions and then close the contour, picking up the corresponding residues.

## B.3 Non-analyticity and dispersive integral in the energy domain

In section 4.1, by performing the spectral integral in the mass domain, we showed that the effect of the particle production poles on the off-shell three-point function is to rotate its external energy $F^{(3)}(z, \mu) \to F^{(3)}(ze^{+i\pi}, \mu)$. Similarly, in section 4.2 for the four-point function, we have observed that both the factorised contribution $F_{-+}$ and the collider contribution of the nested channel $F_{++}^{\text{collider}}$ are connected through the analytic continuation of external energies. In the following, we show that this connection is manifested through a dispersive integral in the energy domain [28].
**Complex energy domain.** Let us define the off-shell three-point function by its integral representation only

$$F^{(3)}(z, \mu) = -i \frac{\sqrt{\pi}}{2} e^{\frac{\pi \mu}{2} + \frac{i \pi}{4}} \sqrt{K} \int_{-\infty}^{0} \frac{d\tau e^{iE\tau}}{(-\tau)^{1/2}} H_{i\mu}^{(2)}(-K\tau), \quad (B.24)$$

where $z \equiv E/K$, with $E$ and $K$ denoting the external and internal energies, respectively. The physical kinematic region is given by $z \ge 1$. We now want to analytically continue this function in the entire complex $z$ plane, and in particular outside the physical region, without actually evaluating the integral explicitly (the result is given in (4.3)). First, from the asymptotic form of the Hankel function, we notice that the UV convergence as $\tau \to -\infty$ is controlled by the phase of $E+K$, i.e. the sum of the energies entering this three-point function. The integral converges only when $E+K$ has a negative imaginary part. As well explained in [28], the integral can be regularised by deforming the contour $-\infty \to -\infty+(E+K)^{-1}$, with $-\infty^+ \equiv -\infty(1 - i\epsilon)$, so that the lower bound $\tau = -\infty$ is approached from a direction that depends on the argument of $E+K$. Then, we note that the integrand exhibits a branch cut, from the Hankel function and the square root, that we choose to be on the positive real axis $\tau \in (0, +\infty)$, as usual. This gives rise to a discontinuity in the energy domain along $E+K \in (-\infty, 0)$, equivalently along $z \in (-\infty, -1)$.

**Discontinuity.** We now show that the full function $F^{(3)}(z, \mu)$ can be recovered from the knowledge of its discontinuity only. By slightly deforming the time integration, the discontinuity of the integral along the branch cut $z \in (-\infty, -1)$ is given by the discontinuity of the integrand itself

$$\begin{aligned} \mathrm{Disc}_z [F^{(3)}(z, \mu)] &= \lim_{\epsilon \to 0} [F^{(3)}(ze^{-i\epsilon}, \mu) - F^{(3)}(ze^{+i\epsilon}, \mu)] \\ &= i \frac{\sqrt{\pi}}{2} e^{\frac{\pi \mu}{2} + \frac{i \pi}{4}} \sqrt{K} \int_{\infty}^{0} d\tau \mathrm{Disc}_{\tau} \left[ \frac{e^{iE\tau}}{\tau^{1/2}} H_{i\mu}^{(2)}(K\tau) \right]. \end{aligned} \quad (B.25)$$

Using the known analytic continuation of the Hankel function (B.6) and (B.8), the discontinuity of the integrand is given by

$$\mathrm{Disc}_{\tau} \left[ \frac{e^{iE\tau}}{\tau^{1/2}} H_{i\mu}^{(2)}(K\tau) \right] = -2i \cosh(\pi\mu) \frac{e^{iE\tau}}{\tau^{1/2}} H_{i\mu}^{(2)}(K\tau). \quad (B.26)$$

Eventually, after changing variables, we obtain

$$\mathrm{Disc}_z [F^{(3)}(z, \mu)] = -2i \cosh(\pi\mu) F^{(3)}(-z, \mu). \quad (B.27)$$

**Dispersive integral.** At fixed $\mu$, since the function $F^{(3)}(z, \mu)$ is analytic everywhere except on the branch cut, we can use Cauchy's integral formula to write

$$F^{(3)}(z, \mu) = \oint_C \frac{dz' F^{(3)}(z', \mu)}{2i\pi (z' - z)} = \int_{-\infty}^{-1} \frac{dz' \mathrm{Disc}_{z'} [F^{(3)}(z', \mu)]}{2i\pi (z' - z)}, \quad (B.28)$$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Illustration of the contour deformation used in (B.28) in the complex $z'$ plane. At fixed mass $\mu$, the off-shell three-point function $F^{(3)}(z, \mu)$ in the physical kinematic region $z \ge 1$ is related to its discontinuity along the branch cut $z \le -1$ through a dispersive integral. At a fixed internal energy $K$ flowing in the off-shell leg, the point $z = +1$ is the folded limit $E \to K$ whereas the branch point $z = -1$ is the standard early-time energy singularity $E + K \to 0$.</figcaption></figure>

where the integral contour $C$ encircles $z$ counterclockwise. The second equality is found by deforming the contour as shown in figure 5, as the integral along the large arc vanishes. Using the found discontinuity of the three-point function (B.27), we finally get

$$F^{(3)}(z, \mu) = -2i \cosh(\pi\mu) \int_{-\infty}^{-1} \frac{dz'}{2i\pi} \frac{F^{(3)}(-z', \mu)}{z' - z} . \qquad (B.29)$$

Inversely, $F^{(3)}(-z, \mu)$ is obtained from $F^{(3)}(z, \mu)$ in the same way. The same approach can be employed to reconstruct the complete four-point correlator, including the tower of EFT poles, using only the non-local signal. This powerful method was used and further explored in [28].

**Data Availability Statement.** This article has no associated data or the data will not be deposited.

**Code Availability Statement.** This article has no associated code or the code will not be deposited.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

References

[1] F.A. Berends and W.T. Giele, *Recursive Calculations for Processes with n Gluons*, Nucl. Phys. B 306 (1988) 759 [INSPIRE].

[2] R. Britto, F. Cachazo and B. Feng, *New recursion relations for tree amplitudes of gluons*, Nucl. Phys. B 715 (2005) 499 [hep-th/0412308] [INSPIRE].
[3] R. Britto, F. Cachazo, B. Feng and E. Witten, *Direct proof of tree-level recursion relation in Yang-Mills theory*, Phys. Rev. Lett. **94** (2005) 181602 [hep-th/0501052] [INSPIRE].

[4] Z. Bern, L.J. Dixon, D.C. Dunbar and D.A. Kosower, *One loop n point gauge theory amplitudes, unitarity and collinear limits*, Nucl. Phys. **B** **425** (1994) 217 [hep-ph/9403226] [INSPIRE].

[5] Z. Bern, L.J. Dixon, D.C. Dunbar and D.A. Kosower, *Fusing gauge theory tree amplitudes into loop amplitudes*, Nucl. Phys. **B** **435** (1995) 59 [hep-ph/9409265] [INSPIRE].

[6] Z. Bern and Y.-T. Huang, *Basics of Generalized Unitarity*, J. Phys. **A** **44** (2011) 454003 [arXiv:1103.1869] [INSPIRE].

[7] P. Benincasa, *New structures in scattering amplitudes: a review*, Int. J. Mod. Phys. **A** **29** (2014) 1430005 [arXiv:1312.5583] [INSPIRE].

[8] H. Elvang and Y.-T. Huang, *Scattering Amplitudes*, arXiv:1308.1697 [INSPIRE].

[9] S. Weinzierl, *Tales of 1001 Gluons*, Phys. Rept. **676** (2017) 1 [arXiv:1610.05318] [INSPIRE].

[10] C. Cheung, *TASI lectures on scattering amplitudes*, in the proceedings of the *Theoretical Advanced Study Institute in Elementary Particle Physics: Anticipating the Next Discoveries in Particle Physics*, Boulder, U.S.A., June 6 – July 1 (2016). [DOI:10.1142/9789813233348_0008] [arXiv:1708.03872] [INSPIRE].

[11] N. Arkani-Hamed and J. Maldacena, *Cosmological Collider Physics*, arXiv:1503.08043 [INSPIRE].

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

[13] G.L. Pimentel and D.-G. Wang, *Boostless cosmological collider bootstrap*, JHEP **10** (2022) 177 [arXiv:2205.00013] [INSPIRE].

[14] S. Jazayeri and S. Renaux-Petel, *Cosmological bootstrap in slow motion*, JHEP **12** (2022) 137 [arXiv:2205.10340] [INSPIRE].

[15] Z. Qin and Z.-Z. Xianyu, *Helical inflation correlators: partial Mellin-Barnes and bootstrap equations*, JHEP **04** (2023) 059 [arXiv:2208.13790] [INSPIRE].

[16] Z. Qin and Z.-Z. Xianyu, *Closed-form formulae for inflation correlators*, JHEP **07** (2023) 001 [arXiv:2301.07047] [INSPIRE].

[17] S. Aoki et al., *Cosmological correlators with double massive exchanges: bootstrap equation and phenomenology*, JHEP **09** (2024) 176 [arXiv:2404.09547] [INSPIRE].

[18] C. Sleight, *A Mellin Space Approach to Cosmological Correlators*, JHEP **01** (2020) 090 [arXiv:1906.12302] [INSPIRE].

[19] C. Sleight and M. Taronna, *Bootstrapping Inflationary Correlators in Mellin Space*, JHEP **02** (2020) 098 [arXiv:1907.01143] [INSPIRE].

[20] C. Sleight and M. Taronna, *From AdS to dS exchanges: Spectral representation, Mellin amplitudes, and crossing*, Phys. Rev. **D** **104** (2021) L081902 [arXiv:2007.09993] [INSPIRE].

[21] C. Sleight and M. Taronna, *From dS to AdS and back*, JHEP **12** (2021) 074 [arXiv:2109.02725] [INSPIRE].

[22] Z. Qin and Z.-Z. Xianyu, *Phase information in cosmological collider signals*, JHEP **10** (2022) 192 [arXiv:2205.01692] [INSPIRE].
[23] Z.-Z. Xianyu and J. Zang, Inflation correlators with multiple massive exchanges, JHEP 03 (2024) 070 [arXiv:2309.10849] [INSPIRE].

[24] D. Werth, L. Pinol and S. Renaux-Petel, Cosmological Flow of Primordial Correlators, Phys. Rev. Lett. 133 (2024) 141002 [arXiv:2302.00655] [INSPIRE].

[25] L. Pinol, S. Renaux-Petel and D. Werth, The Cosmological Flow: A Systematic Approach to Primordial Correlators, arXiv:2312.06559 [INSPIRE].

[26] D. Werth, L. Pinol and S. Renaux-Petel, CosmoFlow: Python Package for Cosmological Correlators, Class. Quant. Grav. 41 (2024) 175015 [arXiv:2402.03693] [INSPIRE].

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

[28] H. Liu, Z. Qin and Z.-Z. Xianyu, Dispersive Bootstrap of Massive Inflation Correlators, arXiv:2407.12299 [INSPIRE].

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

[30] X. Chen, Y. Wang and Z.-Z. Xianyu, Standard Model Mass Spectrum in Inflationary Universe, JHEP 04 (2017) 058 [arXiv:1612.08122] [INSPIRE].

[31] X. Chen, Y. Wang and Z.-Z. Xianyu, Neutrino Signatures in Primordial Non-Gaussianities, JHEP 09 (2018) 022 [arXiv:1805.02656] [INSPIRE].

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

[33] T. Heckelbacher, I. Sachs, E. Skvortsov and P. Vanhove, Analytical evaluation of cosmological correlation functions, JHEP 08 (2022) 139 [arXiv:2204.07217] [INSPIRE].

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

[35] Z. Qin and Z.-Z. Xianyu, Nonanalyticity and on-shell factorization of inflation correlators at all loop orders, JHEP 01 (2024) 168 [arXiv:2308.14802] [INSPIRE].

[36] D. Baumann et al., Snowmass White Paper: The Cosmological Bootstrap, in the proceedings of the Snowmass 2021, Seattle, U.S.A., July 17–26 (2022) [arXiv:2203.08121] [INSPIRE].

[37] J.M. Maldacena and G.L. Pimentel, On graviton non-Gaussianities during inflation, JHEP 09 (2011) 045 [arXiv:1104.2846] [INSPIRE].

[38] S. Raju, New Recursion Relations and a Flat Space Limit for AdS/CFT Correlators, Phys. Rev. D 85 (2012) 126009 [arXiv:1201.6449] [INSPIRE].

[39] H. Goodhew, S. Jazayeri and E. Pajer, The Cosmological Optical Theorem, JCAP 04 (2021) 021 [arXiv:2009.02898] [INSPIRE].

[40] S. Céspedes, A.-C. Davis and S. Melville, On the time evolution of cosmological correlators, JHEP 02 (2021) 012 [arXiv:2009.07874] [INSPIRE].

[41] D. Baumann et al., Linking the singularities of cosmological correlators, JHEP 09 (2022) 010 [arXiv:2106.05294] [INSPIRE].

[42] S. Melville and E. Pajer, Cosmological Cutting Rules, JHEP 05 (2021) 249 [arXiv:2103.09832] [INSPIRE].

[43] H. Goodhew, S. Jazayeri, M.H.G. Lee and E. Pajer, Cutting cosmological correlators, JCAP 08 (2021) 003 [arXiv:2104.06587] [INSPIRE].
[44] M. Hogervorst, J. Penedones and K.S. Vaziri, Towards the non-perturbative cosmological bootstrap, JHEP 02 (2023) 162 [arXiv:2107.13871] [INSPIRE].

[45] M. Loparco, J. Penedones, K. Salehi Vaziri and Z. Sun, The Källén-Lehmann representation in de Sitter spacetime, JHEP 12 (2023) 159 [arXiv:2306.00090] [INSPIRE].

[46] N. Arkani-Hamed, P. Benincasa and A. Postnikov, Cosmological Polytopes and the Wavefunction of the Universe, arXiv:1709.02813 [INSPIRE].

[47] N. Arkani-Hamed and P. Benincasa, On the Emergence of Lorentz Invariance and Unitarity from the Scattering Facet of Cosmological Polytopes, arXiv:1811.01125 [INSPIRE].

[48] P. Benincasa, From the flat-space S-matrix to the Wavefunction of the Universe, arXiv:1811.02515 [INSPIRE].

[49] P. Benincasa, Cosmological Polytopes and the Wavefunction of the Universe for Light States, arXiv:1909.02517 [INSPIRE].

[50] P. Benincasa, A.J. McLeod and C. Vergu, Steinmann Relations and the Wavefunction of the Universe, Phys. Rev. D 102 (2020) 125004 [arXiv:2009.03047] [INSPIRE].

[51] P. Benincasa and W.J.T. Bobadilla, Physical representations for scattering amplitudes and the wavefunction of the universe, SciPost Phys. 12 (2022) 192 [arXiv:2112.09028] [INSPIRE].

[52] P. Benincasa and G. Dian, The Geometry of Cosmological Correlators, arXiv:2401.05207 [INSPIRE].

[53] P. Benincasa and F. Vazão, The Asymptotic Structure of Cosmological Integrals, arXiv:2402.06558 [INSPIRE].

[54] A. Hillman, Symbol Recursion for the dS Wave Function, arXiv:1912.09450 [INSPIRE].

[55] S. De and A. Pokraka, Cosmology meets cohomology, JHEP 03 (2024) 156 [arXiv:2308.03753] [INSPIRE].

[56] N. Arkani-Hamed et al., Differential Equations for Cosmological Correlators, arXiv:2312.05303 [INSPIRE].

[57] B. Fan and Z.-Z. Xianyu, Cosmological Amplitudes in Power-Law FRW Universe, arXiv:2403.07050 [INSPIRE].

[58] H. Gomez, R.L. Jusinskas and A. Lipstein, Cosmological Scattering Equations, Phys. Rev. Lett. 127 (2021) 251604 [arXiv:2106.11903] [INSPIRE].

[59] H. Gomez, R. Lipinski Jusinskas and A. Lipstein, Cosmological scattering equations at tree-level and one-loop, JHEP 07 (2022) 004 [arXiv:2112.12695] [INSPIRE].

[60] C. Armstrong et al., Effective field theories and cosmological scattering equations, JHEP 08 (2022) 054 [arXiv:2204.08931] [INSPIRE].

[61] S. Melville and G.L. Pimentel, A de Sitter S-matrix for the masses, arXiv:2309.07092 [INSPIRE].

[62] C. Chowdhury et al., The Subtle Simplicity of Cosmological Correlators, arXiv:2312.13803 [INSPIRE].

[63] S. Melville and G.L. Pimentel, A de Sitter S-matrix from amputated cosmological correlators, JHEP 08 (2024) 211 [arXiv:2404.05712] [INSPIRE].

[64] D. Marolf, I.A. Morrison and M. Srednicki, Perturbative S-matrix for massive scalar fields in global de Sitter space, Class. Quant. Grav. 30 (2013) 155023 [arXiv:1209.6039] [INSPIRE].
[65] L. Di Pietro, V. Gorbenko and S. Komatsu, Analyticity and unitarity for cosmological correlators, JHEP 03 (2022) 023 [arXiv:2108.01695] [INSPIRE].

[66] L. Di Pietro, V. Gorbenko and S. Komatsu, Cosmological Correlators at Finite Coupling, arXiv:2312.17195 [INSPIRE].

[67] D. Meltzer, The inflationary wavefunction from analyticity and factorization, JCAP 12 (2021) 018 [arXiv:2107.10266] [INSPIRE].

[68] M.H.G. Lee, From amplitudes to analytic wavefunctions, JHEP 03 (2024) 058 [arXiv:2310.01525] [INSPIRE].

[69] Y. Donath and E. Pajer, The in-out formalism for in-in correlators, JHEP 07 (2024) 064 [arXiv:2402.05999] [INSPIRE].

[70] S. Weinberg, Quantum contributions to cosmological correlations, Phys. Rev. D 72 (2005) 043514 [hep-th/0506236] [INSPIRE].

[71] X. Chen, Y. Wang and Z.-Z. Xianyu, Schwinger-Keldysh Diagrammatics for Primordial Perturbations, JCAP 12 (2017) 006 [arXiv:1703.10166] [INSPIRE].

[72] T. Basile, X. Bekaert and N. Boulanger, Mixed-symmetry fields in de Sitter space: a group theoretical glance, JHEP 05 (2017) 081 [arXiv:1612.08166] [INSPIRE].

[73] D. Karateev, P. Kravchuk and D. Simmons-Duffin, Harmonic Analysis and Mean Field Theory, JHEP 10 (2019) 217 [arXiv:1809.05111] [INSPIRE].

[74] G. Sengör and C. Skordis, Unitarity at the Late time Boundary of de Sitter, JHEP 06 (2020) 041 [arXiv:1912.09885] [INSPIRE].

[75] Z. Sun, A note on the representations of SO(1, d + 1), arXiv:2111.04591 [INSPIRE].

[76] J. Penedones, Writing CFT correlation functions as AdS scattering amplitudes, JHEP 03 (2011) 025 [arXiv:1011.1485] [INSPIRE].

[77] U. Moschella and R. Schaeffer, Quantum theory on Lobatchevski spaces, Class. Quant. Grav. 24 (2007) 3571 [arXiv:0709.2795] [INSPIRE].

[78] S. Duary and S. Maji, Spectral representation in Klein space: simplifying celestial leaf amplitudes, JHEP 08 (2024) 079 [arXiv:2406.02342] [INSPIRE].

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

[80] X. Tong, Y. Wang and Y. Zhu, Cutting rule for cosmological collider signals: a bulk evolution perspective, JHEP 03 (2022) 181 [arXiv:2112.03448] [INSPIRE].

[81] S. Weinberg, The quantum theory of fields. Volume 2: Modern applications, Cambridge University Press (2013) [DOI:10.1017/CBO9781139644174] [INSPIRE].

[82] M.D. Schwartz, Quantum Field Theory and the Standard Model, Cambridge University Press (2014) [DOI:10.1017/9781139540940].

[83] D. Meltzer and A. Sivaramakrishnan, CFT unitarity and the AdS Cutkosky rules, JHEP 11 (2020) 073 [arXiv:2008.11730] [INSPIRE].

[84] R. Britto, C. Duhr, H.S. Hannesdottir and S. Mizera, Cutting-Edge Tools for Cutting Edges, arXiv:2402.19415 [INSPIRE].

[85] M.J.G. Veltman, Diagrammatica: The Path to Feynman rules, Cambridge University Press (2012) [DOI:10.1017/cbo9780511564079] [INSPIRE].
[86] R.E. Cutkosky, *Singularities and discontinuities of Feynman amplitudes*, J. Math. Phys. 1 (1960) 429 [INSPIRE].

[87] G. 't Hooft and M.J.G. Veltman, *Diagrammar*, NATO Sci. Ser. B 4 (1974) 177 [INSPIRE].

[88] L.V. Keldysh, *Diagram technique for nonequilibrium processes*, Zh. Eksp. Teor. Fiz. 47 (1964) 1515 [INSPIRE].

[89] A. Kamenev, *Field Theory of Non-Equilibrium Systems*, Cambridge University Press (2011) [DOI:10.1017/cbo9781139003667].

[90] M.J.G. Veltman, *Unitarity and causality in a renormalizable field theory with unstable particles*, Physica 29 (1963) 186 [INSPIRE].

[91] R.J. Eden, P.V. Landshoff, D.I. Olive and J.C. Polkinghorne, *The analytic S-matrix*, Cambridge University Press, Cambridge (1966).

[92] M. Gillioz, X. Lu and M.A. Luty, *Scale Anomalies, States, and Rates in Conformal Field Theory*, JHEP 04 (2017) 171 [arXiv:1612.07800] [INSPIRE].

[93] M. Gillioz, X. Lu and M.A. Luty, *Graviton Scattering and a Sum Rule for the c Anomaly in 4D CFT*, JHEP 09 (2018) 025 [arXiv:1801.05807] [INSPIRE].

[94] M. Gillioz, *Momentum-space conformal blocks on the light cone*, JHEP 10 (2018) 125 [arXiv:1807.07003] [INSPIRE].

[95] J.L. Bourjaily et al., *Sequential Discontinuities of Feynman Integrals and the Monodromy Group*, JHEP 01 (2021) 205 [arXiv:2007.13747] [INSPIRE].

[96] I. Thompson, *NIST Handbook of Mathematical Functions*, edited by Frank W.J. Olver, Daniel W. Lozier, Ronald F. Boisvert, Charles W. Clark, *Contemp. Phys.* 52 (2011) 497.

[97] R.P. Feynman and F.L. Vernon Jr., *The theory of a general quantum system interacting with a linear dissipative system*, Annals Phys. 24 (1963) 118 [INSPIRE].

[98] A. Kaya, *On iɛ Prescription in Cosmology*, JCAP 04 (2019) 002 [arXiv:1810.12324] [INSPIRE].

[99] M. Baumgart and R. Sundrum, *Manifestly Causal In-In Perturbation Theory about the Interacting Vacuum*, JHEP 03 (2021) 080 [arXiv:2010.10785] [INSPIRE].

[100] S. Albayrak, P. Benincasa and C. Duaso Pueyo, *Perturbative unitarity and the wavefunction of the Universe*, SciPost Phys. 16 (2024) 157 [arXiv:2305.19686] [INSPIRE].

[101] H.P. Breuer and F. Petruccione, *The theory of open quantum systems*, Oxford University Press (2002) [DOI:10.1093/acprof:oso/9780199213900.001.0001].

[102] E.A. Calzetta and B.-L.B. Hu, *Nonequilibrium Quantum Field Theory*, Cambridge University Press (2008) [DOI:10.1017/cbo9780511535123].

[103] Y.E. Gutiérrez-Tovar and J.M.R. Méndez-Pérez, *The Kontorovich-Lebedev integral transformation with a Hankel function kernel in a space of generalized functions of doubly exponential descent*, J. Math. Anal. Appl. 328 (2007) 359.