Spectral Representation of Cosmological Correlators
Denis Werth




 =15.5pt
empty

0 pt






1818 Spectral Representation of
 
1818 Cosmological Correlators




15pt


1218Denis Werth[werth@iap.frwerth@iap.fr]






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








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ϵ 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.













1.2



1.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. 

4pt
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 <cit.>, 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 <cit.>, or unitarity methods and generalised cuts <cit.> 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 <cit.> for nice reviews).

4pt
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.

4pt
Despite these challenges, several powerful techniques have recently been developed. At tree level, approaches such as bootstrapping correlators by solving boundary kinematic equations <cit.> and reformulating perturbative diagrammatic rules in Mellin space <cit.> have proven effective. The use of partial Mellin-Barnes representations <cit.>, 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 <cit.>. At the loop order, dispersive integrals have been used to tame loop diagrams <cit.> (also see <cit.> for other techniques).[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.] 

4pt
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 <cit.>. 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 <cit.>. 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 <cit.>, while its non-perturbative counterpart is captured through the Källén–Lehmann representation in de Sitter space <cit.>. 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 <cit.> and progress has been made in understanding the differential equations they satisfy, from which “time emerges" <cit.> (with related aspects discussed in <cit.>). As 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.

4pt
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.[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 <cit.>.] 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 <cit.>. 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.

4pt
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 <cit.> (with previous spectral representations in de Sitter proposed in <cit.>).[Integral representations of cosmological propagators have been proposed in e.g. <cit.>, 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.] 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ϵ 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. 

4pt
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.

4pt
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 Sec. <ref>, 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 Sec. <ref>, 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 Sec. <ref>, we explicitly perform the spectral integral to obtain a new representation of a massive exchange correlator in de Sitter. We summarise our conclusions in Sec. <ref>.




  
Notation and Conventions. We use the mostly-plus signature for the metric (-, +, +, +). Spatial three-dimensional vectors are written in boldface k. We use  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τ = dt/a, by τ. Throughout most of the paper, we work in the Poincaré patch of de Sitter space with the metric ṣ^2 = a^2(τ)(-τ̣^2+x^2), where a(τ) = -(Hτ)^-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.







§ 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.




 §.§ Cosmological propagators



We are interested in connected equal-time correlation functions of a bulk scalar field φ(t, x). The field φ 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 ℒ[φ]. 




  
Generating functional. Following the standard Schwinger-Keldysh path integral formalism to compute correlation functions <cit.>, two copies φ_± of the field φ 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 φ_± should coincide at t, i.e. φ_+(t, x) = φ_-(t, x). Equal-time correlation functions ⟨Ω|φ̂_a_1(t, x_1) …φ̂_a_n(t, x_n)|Ω|$⟩, witha_1, …, a_n=±and where|Ω⟩is the vacuum of the fully interacting theory, are computed by taking functional derivatives 

    ⟨Ω|φ̂_a_1(t, x_1) …φ̂_a_n(t, x_n)|Ω|=⟩1/ia_1… ia_n δ^n Z[J_+, J_-]/δ J_a_1(t, x_1) …δ J_a_n(t, x_n)|_J_±=0 ,

where the generating functionZ[J_+, J_-]is defined as

    Z[J_+, J_-] = ∫_𝒞_iϵφ_± exp(i∫_-∞^tṭ' ^̣3x (ℒ[φ_+] + J_+ φ_+ - ℒ[φ_-] - J_- φ_-))  .

The contour_iϵdenotes the usual in-in contour that goes in the upper-half complex plane from-∞^+totand reverts towards-∞^-in the lower-half complex plane, with-∞^±≡-∞(1∓iϵ)andϵ>0an 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. (<ref>), the contour encodes the standardiϵprescription to which we will come back later.




  
Free propagators. Resorting to perturbation theory, the full Lagrangian is split into a free partℒ_0that can be solved exactly, defining the free propagators, and an interacting partℒ_int, i.e.ℒ = ℒ_0 + ℒ_int. The first simplest object that can be defined is the free Wightman function⟨0|φ̂(t_1, x_1) φ̂(t_2, x_1)|0|$⟩ where the expectation value is over the vacuum of the free theory |0⟩. In Fourier space for the spatial coordinates only, after expanding the operator φ̂_k(t) in terms of annihilation and creation operators 

    φ̂_k(t) = u_k(t)  â_k + u_k^*(t) â^†_-k ,

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

    (k; t_1, t_2) = u_k(t_1) u_k^*(t_2) .

We have used invariance under spatial translations and rotations for a homogeneous and isotropic background so that the function  only depends on a single 3-momentum magnitude k ≡ |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 _ab(k; t_1, t_2) where a, b = ±, defined as

    _++(k; t_1, t_2)    = (k; t_1, t_2) Θ(t_1-t_2) + ^*(k; t_1, t_2) Θ(t_2-t_1) ,
    _+-(k; t_1, t_2)    = ^*(k; t_1, t_2) ,
    _-+(k; t_1, t_2)    = (k; t_1, t_2) , 
    _–(k; t_1, t_2)    = ^*(k; t_1, t_2) Θ(t_1-t_2) + (k; t_1, t_2) Θ(t_2-t_1) .

These propagators being not independent is a consequence of the Wightman function satisfying the property (k; t_2, t_1) = ^*(k; t_1, t_2), as can be seen from Eq. (<ref>). From these “bulk-to-bulk" propagators, one can define “bulk-to-boundary" propagators by sending t_2→ t (for t_1<t) and using the condition that the two in-in branches coincide at the external time φ_+(t, x) = φ_-(t, x). They read

    _+(k; t_1) =  ^*(k; t_1, t)  ,   _-(k; t_1) = (k; t_1, t)  .

Following the diagrammatic rules exposed in <cit.>, we represent these propagators with a black dot 0pt[line width=1. pt, scale=2]
[fill=black] (0, 0) circle (.05cm); denoting +, a white dot 0pt[line width=1. pt, scale=2]
[fill=white] (0, 0) circle (.05cm); denoting -, and a square 0pt[line width=1. pt, scale=2]
[fill=white] ([xshift=0pt,yshift=-1.5pt]1, 0) rectangle ++(3pt,3pt); denoting the boundary at external time t where + and - are indistinguishable:

    [c]
            0pt[line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm) node[above=0.5mm] t_1;
    [black] (0.05, 0) – (1, 0);
    [fill=black] (1, 0) circle (.05cm) node[above=0.5mm] t_2;
       =    _++(k; t_1, t_2) , 
    0pt[line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm) node[above=0.5mm] t_1;
    [black] (0.05, 0) – (1, 0);
    [fill=white] (1, 0) circle (.05cm) node[above=0.5mm] t_2;
       =    _+-(k; t_1, t_2) , 
    0pt[line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm) node[above=0.5mm] t_1;
    [black] (0.05, 0) – (1, 0);
    [fill=black] (1, 0) circle (.05cm) node[above=0.5mm] t_2;
       =    _-+(k; t_1, t_2) , 
    0pt[line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm) node[above=0.5mm] t_1;
    [black] (0.05, 0) – (1, 0);
    [fill=white] (1, 0) circle (.05cm) node[above=0.5mm] t_2;
       =    _–(k; t_1, t_2) ,
            [c]
        0pt[line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm) node[above=0.5mm] t_1;
    [black] (0.05, 0) – (1, 0);
    [fill=white] ([xshift=0pt,yshift=-1.5pt]1, 0) rectangle ++(3pt,3pt);
       =    _+(k; t_1) , 
    0pt[line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm) node[above=0.5mm] t_1;
    [black] (0.05, 0) – (1, 0);
    [fill=white] ([xshift=0pt,yshift=-1.5pt]1, 0) rectangle ++(3pt,3pt);
       =    _-(k; t_1) .

A double-coloured dot 0pt[line width=1. pt, scale=2]
 (0, 0) circle (.05cm);
[black] (0,0) – (90:0.05) arc (90:270:0.05) – cycle; means either a black or a white dot. 




 §.§ 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 φ 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) = e^-ik t√(2k). The generalisation to massive fields is straightforward after setting k→ω_k=√(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)Θ(t_1-t_2) + e^+ik(t_1 - t_2)Θ(t_2-t_1) = -2k ∫_-∞^+∞ω̣/2iπe^iω(t_1-t_2)/(ω^2 - k^2)_iϵ ,

where the iϵ prescription is implemented as follows

    1/(ω^2-k^2)_iϵ≡1/2k[1/ω - (k-iϵ) - 1/ω - (-k+iϵ)] = 1/ω^2-k^2+iϵ ,

with ϵ>0 an infinitesimal positive parameter, and where in the second equality we have dropped terms of order 𝒪(ϵ^2) and set 2kϵ = ϵ, which is valid in the limit ϵ→ 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. (<ref>) gives zero and only the second term selects the residue, effectively setting the exponential on shell, i.e. e^iω(t_1-t_2)→ e^-ik(t_1-t_2). One can proceed in the same way for t_1<t_2, this way closing the contour in the lower-half complex plane, to find e^ik(t_1-t_2). Importantly, this requires that the mode function should be analytic, which is valid here for the exponential in the entire complex plane, and that the integrand decays at least faster than 1/ω when ω→∞ so that the semi-circle at infinity does not contribute.

4pt
In the end, the time-ordered propagator _++ can be re-written as

    _++(k; t_1, t_2) = i ∫_-∞^+∞ω̣/2πũ^*_ω(t_1) ũ_ω(t_2)/(ω^2 - k^2)_iϵ ,

where ũ_k(t) = e^-ikt is the mode function after stripping off the overall normalisation. Naturally, we have _–(k; t_1, t_2) = _++^*(k; t_1, t_2). Some comments are in order. First, the use of tilted mode functions ũ_k(t) avoids introducing unnecessary branch cuts in the integrand coming from √(k), thereby preserving its analytic property. Second, both times t_1 and t_2 appear factorised so that time integrals within correlators become trivial. This property will be illustrated with concrete examples in Sec. <ref>. Finally, the iϵ prescription is just a simple way for representing time ordering and specifying a pole prescription. Equivalently, one may deform the contour slightly from -∞^- = -∞(1+iϵ) to +∞^+ = +∞(1+iϵ), while keeping the poles on the real axis. We show in App. <ref> how this iϵ prescription can be recovered directly from the Schwinger-Keldysh path integral.




  
Retarded and advanced propagators. Various Green's functions with different boundary conditions can be recovered from the usual integral representation of the Feynman propagator. This can be easily done by deforming and combining diverse integration contours. For example, using the standard identity satisfied by the Heaviside function Θ(t_1-t_2) + Θ(t_2-t_1) = 1,
the time-ordered propagator _++ can be written as

    _++(k; t_1, t_2) = ^*(k; t_1, t_2) + [(k; t_1, t_2) - ^*(k; t_1, t_2)]Θ(t_1-t_2) .

The first term is the non-time-ordered propagator _+-(k; t_1, t_2) and therefore is factorised in time, whereas the second one is nothing but the causal retarded Green's function _R(k; t_1, t_2). In terms of integration contours, Eq. (<ref>) reads

    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, -0.1) circle (.02cm);
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0) – (0.4, 0);
    
    =
            [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
    [black, fill = black] (-0.2, 0) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (0.3, 0) arc (0:-360:0.1) – (0.3, 0);
     - 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
    [black, fill = black] (-0.2, 0) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, -0.2) – (0.4, -0.2);
     .

Importantly, non-local processes lead to a vanishing retarded Greens's function. This observation has consequences for the cosmological case, and enables to extract various signals from correlators with different origin. Similarly, one can flip the second Heaviside function in the expression of _++ to obtain

    _++(k; t_1, t_2) = (k; t_1, t_2) + [^*(k; t_1, t_2) - (k; t_1, t_2)]Θ(t_2-t_1) ,

which, in terms of integration contours, reads

    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, -0.1) circle (.02cm);
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0) – (0.4, 0);
    
    =
            [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
    [black, fill = black] (-0.2, 0) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (-0.1, 0) arc (0:360:0.1) – (-0.1, 0);
     - 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
    [black, fill = black] (-0.2, 0) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0.2) – (0.4, 0.2);
     .






 §.§ 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 App. <ref> for more details.[We consider fields in the principal series m/H≥ 3/2 so that μ≡√(m^2/H^2 - 9/4) is real. In particular, we will not deal with late-time secular divergences appearing when exchanging light fields.]




  
Contour prescription. Spontaneous particle production manifests itself as non-analyticity in the complex energy domain k→ω 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 μ plane, allowing the field to acquire any complex mass, as shown in <cit.>. The integral representation of the time-ordered massive cosmological propagator reads[The spectral parameter ν should not be confused with ν = iμ that is commonly used to label light fields in the complementary series.]

    _++(k; τ_1, τ_2) = i ∫_-∞^+∞ν̣ _ν u_k^*(τ_1, ν) u_k^*(τ_2, ν)/(ν^2 - μ^2)_iϵ ,
  
where _ν≡νπsinh(πν) is the de Sitter density of states, u_k^*(τ, μ) is the negative-frequency massive mode function defined in (<ref>), and the contour prescription to select the correct pole structure is given by[The explicit pole structure (and the corresponding residues) can be immediately obtained from the usual flat-space prescription (<ref>) or by noticing that

    1/ν^2 - μ^2 + iϵ = Γ(iν-iμ-ϵ)Γ(-iν-iμ-ϵ)/Γ(iν-iμ-ϵ+1)Γ(-iν-iμ-ϵ+1) .
]

    1/(ν^2-μ^2)_iϵ≡1/2sinh(πμ)[e^+πμ/ν^2-μ^2+iϵ - e^-πμ/ν^2-μ^2-iϵ] .

Naturally, the anti-time ordered propagator is given by _–(k; τ_1, τ_2) = _++^*(k; τ_1, τ_2). The detailed derivation of (<ref>) is presented in App. <ref>. The integral over the spectral parameter ν can be viewed as an integration over the scaling dimension Δ = 32+iν in the principal series that labels infinite-dimensional (scalar) representations of the de Sitter group, which defines eigenvalues of the quadratic Casimir operator Δ (3-Δ), see e.g. <cit.>. In this sense, this integral is closely related to the de Sitter Källén–Lehmann representation albeit in spatial Fourier space <cit.>. Such spectral representation has its analogue in anti-de Sitter space, see e.g. <cit.>, and in other curved geometries <cit.>. 




  
Flat-space limit. The spectral representation (<ref>) extends the familiar flat-space prescription (<ref>) by capturing spontaneous particle production while at the same time having a well-defined flat-space limit. Indeed, in the limit H→0, the Boltzmann suppression is turned off e^-πμ→0 as μ∼ m/H for a fixed mass. The additional unusual poles ν = ±μ± iϵ vanish and the contour prescription reduces to the flat-space one (<ref>)

    1/(ν^2-μ^2)_iϵ1/ν^2-μ^2+iϵ .

Following the derivation of (<ref>) presented in App. <ref>, for τ_1 being in the future of τ_2, and using limiting forms given in App. <ref>, the massive propagator reduces to

    _++(k; τ_1, τ_2)    = i/2∫_-∞^+∞ν̣ ν J_iν(-kτ_1) H_iν^(2)(-kτ_2)/ν^2-μ^2+iϵ = π/2 J_iμ(-kτ_1) H_iμ^(2)(-kτ_2) 
       e^-iμ(t_1-t_2)/√(2μ) ,

where we have selected the single pole at ν = μ-iϵ 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 ω_k ∼μ. 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(τ, μ) = Γ(1+iμ)/√(2μ)(k/2)^-iμ J_iμ(-kτ) e^-iμ t/√(2μ) .

Notice in passing that the late-time limit τ→0 equivalently also corresponds to the large-mass limit μ→∞ as can easily recovered using (<ref>) and the Stirling formula Γ(z) ∼ e^-z z^z √(2π/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(τ, μ) = α_k u_k(τ, μ) + β_k u_k^*(τ, μ) ,

where the Bogolyubov coefficients are

    α_k = Γ(1+iμ)/√(2πμ)(k/2)^-iμ e^+πμ/2-iπ/2 ,   β_k = Γ(1+iμ)/√(2πμ)(k/2)^-iμ e^-πμ/2+iπ/2 ,

from which we obtain the mean particle density of “out" excitations with wavenumber 
k in the “in" vacuum state  |β_k|^2 = 1/(e^2πμ-1). We recover the usual Bose-Einstein distribution, when identifying the energy with the rest mass μ in a thermal bath set by the de Sitter temperature T_dS = 1/2π.







  
Local and non-local parts. The iϵ prescription (<ref>) 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 (<ref>), so that the propagator can be written <cit.>

    _++(k; τ_1, τ_2) = _++^local(k; τ_1, τ_2) + _++^non-local(k; τ_1, τ_2) ,

where both contributions read

    _++^local(k; τ_1, τ_2)    = +iπ/4∫_-∞^+∞ν̣ _ν/sinh^2(πν) J_iν(-kτ_1) J_iν^*(-kτ_2) + (ν↔-ν)/(ν^2-μ^2)_iϵ , 
    _++^non-local(k; τ_1, τ_2)    = -iπ/4∫_-∞^+∞ν̣ _ν/sinh^2(πν) e^-πνJ_iν(-kτ_1) J_iν(-kτ_2) + (ν↔-ν)/(ν^2-μ^2)_iϵ .

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ν(-kτ_1) J_iν^*(-kτ_2) ∼sinh(πν)/πν(τ_1/τ_2)^iν∝ e^-iν(t_1-t_2) .

Notably, the local part is independent of the momentum so that it describes correlation at coincident points in position space. Assuming τ_1>τ_2, we close the contour in the lower-half complex plane for the term J_iν(-kτ_1) J_iν^*(-kτ_2) and in the upper-half complex plane for the complex conjugate piece, which leads to

    _++^local(k; τ_1, τ_2)    = π/4sinh^2(πμ)[e^+πμ J_iμ(-kτ_1)J_iμ^*(-kτ_2) + e^-πμ J_iμ^*(-kτ_1)J_iμ(-kτ_2)] 
       Γ(+iμ)Γ(-iμ)/4π[e^+πμ(τ_1/τ_2)^+iμ + e^-πμ(τ_1/τ_2)^-iμ] .

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 (τ_1τ_2)^± iν→ e^±πμ(τ_1τ_2)^± iμ + e^∓πμ(τ_1τ_2)^∓ iμ, which then combine to reconstruct the correct mode mixing. This is illustrated in Fig. <ref>. On the other hand, the non-local part corresponds to the creation and the propagation of two particles since

    J_iν(-kτ_1) J_iν(-kτ_2) ∼1/Γ(1+iν)^2(k^2τ_1τ_2/4)^iν∝ e^-iν(t_1-t_⋆) e^-iν(t_2-t_⋆) ,

in the late-time limit, where t_⋆=log k is the pair production (cosmic) time. Being non-analytic in the momentum, this piece describes long-range non-local correlations in position space. Similarly to the local contribution, we close the contour in the lower-half complex plane for the term e^-πνJ_iν(-kτ_1) J_iν(-kτ_2) and in the opposite direction for the second one. We obtain

    _++^non-local(k; τ_1, τ_2)    = π/4sinh^2(πμ)[J_iμ(-kτ_1)J_iμ(-kτ_2) + J_iμ^*(-kτ_1)J_iμ^*(-kτ_2)] 
       1/4π[Γ(-iμ)^2 (kτ_1τ_2/4)^+iμ + Γ(+iμ)^2 (kτ_1τ_2/4)^-iμ] .

Importantly, at late times, the choice of closing the contour is completely fixed so that time ordering does not matter. The Feynman propagator _++ effectively reduces to the Wightman function . Non-local processes being real, as can be explicitly seen from the last equation, they do not enter the advanced or retarded Green's functions. This reflects the fact that the endpoints of a soft propagator are in space-like separation. With the spectral representation of the time-ordered propagator, subtracting the commutator of field operators simply amounts to selecting the residues at the poles in the contour prescription

    _++^non-local(k; τ_1, τ_2)    = _++(k; τ_1, τ_2)    -   _R(k; τ_1, τ_2) , 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0.1) circle (.02cm);
    [black, fill = black] (0.2, -0.1) circle (.02cm);
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
    [black, fill = black] (-0.2, -0.1) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (0.3, -0.1) arc (0:-360:0.1) – (0.3, -0.1);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (-0.1, -0.1) arc (0:-360:0.1) – (-0.1, -0.1);
       =
            [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0.1) circle (.02cm);
    [black, fill = black] (0.2, -0.1) circle (.02cm);
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
    [black, fill = black] (-0.2, -0.1) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0) – (0.4, 0);
     - 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0.1) circle (.02cm);
    [black, fill = black] (0.2, -0.1) circle (.02cm);
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
    [black, fill = black] (-0.2, -0.1) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, -0.2) – (0.4, -0.2);
     .

Of course, selecting the poles in upper-half complex plane requires subtracting the advanced Green's function _A. This contour deformation is completely equivalent to the bulk cutting rules that extract non-local signals <cit.>.




  
Conformally coupled field. The special case of a conformally coupled field in de Sitter m^2=2H^2 is related to the flat-space one by a conformal transformation. Physically, particle production of a conformally coupled field in de Sitter is effectively turned off so that its mode function is analytic in the entire complex k plane. Therefore, from the flat-space propagator (<ref>), one immediately obtains

    _++(k; τ_1, τ_2) = i τ_1 τ_2  ∫_-∞^+∞ω̣/2πũ^*_ω(τ_1) ũ_ω(τ_2)/(ω^2 - k^2)_iϵ ,

where ũ_k(τ) = e^-ikτ and where we have set H=1.







§ COSMOLOGICAL LARGEST-TIME EQUATION



Unitarity at the quantum level, often described as the conservation of probability, imposes powerful constraints on observables, even when the full set of states in the Hilbert space is unknown. At a perturbative level, unitarity implies relations between cosmological correlators <cit.> and conserved quantities under unitary time evolution <cit.>. These relations are analogous to the flat-space optical theorem <cit.>. It has been shown that these relations can be organised into a set of cutting rules <cit.>, which stem from the Hermitian analytic nature of propagators and the straightforward factorisation of a combination of bulk-to-bulk propagators, bearing resemblance to the standard S-matrix Cutkosky rules, see e.g. <cit.> for a recent review. However, these cutting rules were derived for the quantum mechanical wavefunction of the universe, as relations satisfied by the corresponding wavefunction coefficients.

4pt
In this section, we derive similar relations directly at the level of cosmological correlators, and relaxing the assumption of unitarity. We show that individual in-in branch components of cosmological correlators, once external legs are properly analytically continued to negative energies on the negative branch, satisfy the largest-time equation <cit.>. We then showcase how this equation is satisfied for relatively simple diagrams.





 §.§ Propagator identity from an off-shell perspective


Using the spectral representation of bulk-to-bulk propagators presented in the previous section, deriving the largest-time equation exactly parallels that for the flat-space S-matrix <cit.>. The essence of this equation lies in a specific trivial identity satisfied by propagators. Using the standard distributional identity

    1/x+iϵ - 1/x-iϵ = -2iπδ(x) ,

in the limit ϵ→ 0, which follows from the Sokhotski-Plemelj identity, the flat-space bulk-to-bulk propagators (<ref>) satisfy the known relation[In addition, we have used the following identity for the delta Dirac function

    δ(ω^2 - k^2) = 1/2k [δ(ω+k) + δ(ω-k)]  ,

for k>0.]

    _++(k; t_1, t_2) + _–(k; t_1, t_2) = _+-(k; t_1, t_2) + _-+(k; t_1, t_2) ,

that is, since _–(k; t_1, t_2) = _++^*(k; t_1, t_2), the real part does not contain time ordering. In terms of integration contour, this identity reads

    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, -0.1) circle (.02cm);
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0) – (0.4, 0);
     - 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0.1) circle (.02cm);
    [black, fill = black] (-0.2, -0.1) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0) – (0.4, 0);
    
    =
            [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
    [black, fill = black] (-0.2, 0) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (0.3, 0) arc (0:-360:0.1) – (0.3, 0);
     - 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
    [black, fill = black] (-0.2, 0) circle (.02cm);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (-0.1, 0) arc (0:-360:0.1) – (-0.1, 0);
     ,

with the additional factor i imposing the correct closed contour orientation. Without branch cuts in the energy domain, the additional identity _–(k; t_1, t_2) = -_++(-k; t_1, t_2)—which can be found directly from the frequency-space representation using 1/(ω^2 - k^2)_iϵ→1/ω^2-k^2-iϵ as k→ -k or from the explicit split representation in terms of mode functions—enables us to interpret the combination of propagators (<ref>) as cutting the internal line of a diagram on the same in-in branch since time integrals factorise. At the level of correlators, the additional input of unitarity relates individual diagrams to their complex conjugate version (with flipped external energies) on the same branch. Consequently, the latest-time equation we derive should be understood as only a direct outcome of locality, as it is fundamentally determined by the structure of propagators. This general property remains valid in the presence of dissipation, fluctuations and non-unitary time evolution. 

4pt
In de Sitter space, we have seen in Sec. <ref> that the effect of particle production amounts to double the usual poles while assigning distinct weights e^±πμ to the associated residues. From the spectral representation (<ref>), it is easy to observe that bulk-to-bulk propagators also satisfy the identity (<ref>), which, in terms of integration contours, is illustrated by

    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0.1) circle (.02cm);
     at (0.35, 0.3) e^-πμ;
    [black, fill = black] (0.2, -0.1) circle (.02cm);
     at (0.35, -0.3) e^+πμ;
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
     at (-0.35, 0.3) e^+πμ;
    [black, fill = black] (-0.2, -0.1) circle (.02cm);
     at (-0.35, -0.3) e^-πμ;
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0) – (0.4, 0);
     - 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0.1) circle (.02cm);
     at (0.35, 0.3) e^+πμ;
    [black, fill = black] (0.2, -0.1) circle (.02cm);
     at (0.35, -0.3) e^-πμ;
    [black, fill = black] (-0.2, 0.1) circle (.02cm);
     at (-0.35, 0.3) e^-πμ;
    [black, fill = black] (-0.2, -0.1) circle (.02cm);
     at (-0.35, -0.3) e^+πμ;
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.7 with [line width=1pt]>] (-0.4, 0) – (0.4, 0);
    
    =
            [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
     at (0.35, 0.3) e^+πμ;
    [black, fill = black] (-0.2, 0) circle (.02cm);
     at (-0.35, 0.3) e^-πμ;
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (0.3, 0) arc (0:-360:0.1) – (0.3, 0);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (-0.1, 0) arc (0:-360:0.1) – (-0.1, 0);
     - 
    [line width=1. pt, scale=2]
    [black, -] (-0.5,0) – (0.5,0) coordinate (xaxis);
    [black, -] (0,-0.4) – (0,0.4) coordinate (yaxis);
    [black, fill = black] (0.2, 0) circle (.02cm);
     at (0.35, 0.3) e^-πμ;
    [black, fill = black] (-0.2, 0) circle (.02cm);
     at (-0.35, 0.3) e^+πμ;
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (0.3, 0) arc (0:-360:0.1) – (0.3, 0);
    [pyred, draw, line width = 1pt, postaction = decorate, decoration=markings,
    			mark=at position 0.4 with [line width=1pt]>] (-0.1, 0) arc (0:-360:0.1) – (-0.1, 0);
     .

Crucially, note that particle production poles change their respective weight under complex conjugation. Sending μ→∞, one immediately recovers the flat-space limit (<ref>).[The combination _K ≡_+++_– is nothing but the standard Keldysh propagator. The advanced, retarded and Keldysh propagators can be recovered after performing a rotation in the φ_± field basis, introducing a new pair of fields φ_r = 1√(2)(φ_+ + φ_-) and φ_a = 1√(2)(φ_+ - φ_-), which is referred to the Keldysh basis <cit.>. The advantage of this basis is that the physical interpretation of various propagators become transparent. Indeed, _R and _A characterise the dynamical properties of the particles in the system under consideration, while _K characterises their statistical distribution and encodes the mean particle density <cit.>.] The identity (<ref>) can of course be recovered from the form of propagators written with the time-ordered usual representation (<ref>) using Θ(t_1-t_2)+Θ(t_2-t_1)=1. Here, we have derived this identity from an off-shell perspective. 





 §.§ Largest-time equation


Our aim is to translate the propagator identity (<ref>) 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 <cit.>.[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:

    
  * Duplicate the Feynman diagram 2^n times by colouring vertices either in black or in white in all possible ways.
    
  * For each vertex, a black-coloured vertex brings a factor i, and a white-coloured one brings a factor (-i).
    
  * The propagator between two black vertices is the usual Feynman propagator Δ_F(x-y) = Θ(x^0-y^0)Δ^+(x-y) + Θ(y^0-x^0)Δ^-(x-y), where Δ^±(x-y) are the usual positive and negative frequency commutation functions, while the propagator between two white vertices is taken to be Δ_F^*(x-y). The propagator between black and white vertices is replaced by Δ^+(x-y), and Δ^-(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. 

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. <cit.>, or cutting rules, see e.g. <cit.>. 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.
]




  
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

    _a(k; t) ≡ 2k _a(k; t) ,   _ab(k; t, t') ≡ 2k _ab(k; t, t') ,

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

    ∏_i=1^n _a(k_i; t) = _a(∑_i=1^n k_i; t) ,

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

    0pt[line width=1. pt, scale=2]
     (0, 0) circle (.05cm) node[above=1mm] E;
    [black] (0,0) – (90:0.05) arc (90:270:0.05) – cycle;
    ≡ (-i) ∑_a=±a∫_-∞^a^0ṭ  e^i aEt ,

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 (<ref>) on a simple example, namely the two-site chain. Even though the precise nature of the field φ, interactions, and the background is unimportant for what follows, we consider polynomial self-interactions in flat space and set coupling constants to unity. 

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

    0pt
    [line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm) node[above=0.5mm] E_1;
    [fill=black] (1, 0) circle (.05cm) node[above=0.5mm] E_2;
    [black] (0.05, 0) – node[above] K (1, 0);
     + 
    0pt
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm) node[above=0.5mm] -E_1;
    [fill=white] (1, 0) circle (.05cm) node[above=0.5mm] -E_2;
    [black] (0.05, 0) – node[above] K (0.95, 0);
     .

It is given by 

    -∫_-∞^+^0 ṭ_1 ∫_-∞^+^0 ṭ_2   e^i E_1 t_1[_++(K; t_1, t_2) + _–(K; t_1, t_2)] e^i E_2 t_2 .

Note that flipping external energies amounts to exchanging both in-in branches. As such, the iϵ prescription that makes the integrals converge needs to be modified. Therefore, we have also deformed the lower bound of the time integrals accordingly so that both graph integrals can be combined. Using the propagator identity (<ref>) and flipping the external energies (and the corresponding lower bound of the time integral) for each integral appearing on the negative branch, we directly obtain

    -∫_-∞^+^0 ṭ_1 ∫_-∞^-^0 ṭ_2   e^i E_1 t_1_+-(K; t_1, t_2) e^-i E_2 t_2 - ∫_-∞^-^0 ṭ_1 ∫_-∞^+^0 ṭ_2   e^-i E_1 t_1_-+(K; t_1, t_2) e^i E_2 t_2 ,

which we recognise as 

    -0pt
    [line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm) node[above=0.5mm] E_1;
    [fill=white] (1, 0) circle (.05cm) node[above=0.5mm] -E_2;
    [black] (0.05, 0) – node[above] K (0.95, 0);
     - 
    0pt
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm) node[above=0.5mm] -E_1;
    [fill=black] (1, 0) circle (.05cm) node[above=0.5mm] E_2;
    [black] (0.05, 0) – node[above] K (0.95, 0);
     .

Finally, let us group the factorised contributions in (<ref>) with the nested ones on the same side. The resulting diagrammatic formula reads

    0pt
    [line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm) node[above=0.5mm] E_1;
    [fill=black] (1, 0) circle (.05cm) node[above=0.5mm] E_2;
    [black] (0.05, 0) – node[above] K (1, 0);
     + 
    0pt
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm) node[above=0.5mm] -E_1;
    [fill=white] (1, 0) circle (.05cm) node[above=0.5mm] -E_2;
    [black] (0.05, 0) – node[above] K (0.95, 0);
    
        +0pt
    [line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm) node[above=0.5mm] E_1;
    [fill=white] (1, 0) circle (.05cm) node[above=0.5mm] -E_2;
    [black] (0.05, 0) – node[above] K (0.95, 0);
     +
    0pt
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm) node[above=0.5mm] -E_1;
    [fill=black] (1, 0) circle (.05cm) node[above=0.5mm] E_2;
    [black] (0.05, 0) – node[above] K (0.95, 0);
     = 0 .

Note that in this simple flat-space example, the rotation direction to negative energies does not matter as the mode functions are analytic. However, for massive fields in de Sitter space, one needs to specify the analytic continuation to negative energies so that it does not cross any branch cuts. In particular, energies should be rotated in opposite directions for the two branches

    _+(e^-iπ k; t) = _-(k; t)  ,   _-(e^+iπ k; t) = _+(k; t) .

In most of the cases, we are interested in graphs with analytic external mode functions so that this prescription to evade the branch cut and define negative energies remains implicit. 




  
Largest-time equation. The previous simple example actually generalises to arbitrary graphs. In words, the sum of all individual graph contributions with negative external energies on the negative branch vanishes. More formally and at tree-level, individual graph integral contributions can be written as

    _a_1 …a_V({E_i}, {K_j}) ≡ (-i)^V ∫∏_i=1^V ṭ_i  f_i(t_i)  _a_i(E_i; t_i) ∏_j=1^I _a_j a_j+1(K_j; t_j, t_j+1) ,

where V is the number of vertices, I ≡ V-1 is the number of internal lines (at tree-level), E_i and K_j denote external and internal energies, respectively, and f_i(t_i) are form factors associated to each vertex that account for possible (time-dependent) couplings. We leave the dependence on momenta implicit. The labels a_j and a_j+1 are to be understood as two consecutive Schwinger-Keldysh vertex indices that match those of the bulk-to-boundary propagators a_i (for a_j) and a_i+1 (for a_j+1). Taking the values ±, they indicate on which in-in branch the corresponding vertices are located. With these definitions, the latest-time equation can be simply stated as the following schematic identity


    ∑_a_1, …, a_V = ±_a_1 …a_V({a_i E_i}, {K_j}) = 0  .


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 → -k_i also implies flipping the signs of the corresponding momenta k_i → - k_i.[It should be pointed out that the latest-time equation can be understood as a consequence of the following combinatoric operator identity <cit.>

    ∑_k=0^n (-1)^k ∑_σ∈Π(k, n-k)T̅[Ø(t_σ_1) ⋯Ø(t_σ_k)]  T[Ø(t_σ_k+1) ⋯Ø(t_σ_n)] = 0 ,

where T̅ and T are the (anti-)time-ordering operators, and Π(k, n-k) denotes the set of bi-partitions of {1, …, 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.]




  
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 <cit.> where Cutkosky rules were derived for amputated Feynman diagrams, also see <cit.>.

4pt
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 (<ref>) 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 ± labels. For what follows, it is useful to introduce some diagrammatic notation. We denote by the grey circle

    [line width=1. pt, scale=2]
    [fill=gray] (0, 0) circle (0.3cm);
     (0.3, 0) circle (.05cm) node[right=1mm] E_i;
    [black] (0.3,0) – (0.3, 0.05) arc (90:270:0.05) – cycle;
     ,

the sum of all individual graph contributions with V vertices, where the picked vertex 0pt[line width=1. pt, scale=2]
 (0, 0) circle (.05cm);
[black] (0,0) – (90:0.05) arc (90:270:0.05) – cycle; 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

    [line width=1. pt, scale=2]
    [fill=gray] (0, 0) circle (0.3cm);
    [fill=black] (0.3, 0) circle (.05cm) node[above right=0.1mm] E_i;
    [fill=black] (1.3, 0) circle (.05cm) node[above=0.5mm] E;
    [black] (0.35, 0) – node[below] K (1.35, 0);
    
        +
        
    [line width=1. pt, scale=2]
    [fill=gray] (2.5, 0) circle (0.3cm);
    [fill=white] (2.8, 0) circle (.05cm) node[above right=0.1mm] -E_i;
    [fill=white] (3.8, 0) circle (.05cm) node[above=0.5mm] -E;
    [black] (2.85, 0) – node[below] K (3.75, 0);
     .

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 the iϵ prescription). Internal energies remain untouched. This manipulation isolates the bulk-to-bulk propagators to give

    [line width=1. pt, scale=2]
    [fill=gray] (0, 0) circle (0.3cm);
    ×(
        
    [line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm);
    [fill=black] (1, 0) circle (.05cm);
    [black] (0, 0) – (1, 0);
    
        +
        
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm);
    [fill=white] (1, 0) circle (.05cm);
    [black] (0.05, 0) – (0.95, 0);
    ) .

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

    [line width=1. pt, scale=2]
    [fill=gray] (0, 0) circle (0.3cm);
    ×(
        
    [line width=1. pt, scale=2]
    [fill=black] (0, 0) circle (.05cm);
    [fill=white] (1, 0) circle (.05cm);
    [black] (0.05, 0) – (0.95, 0);
    
        +
        
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (.05cm);
    [fill=black] (1, 0) circle (.05cm);
    [black] (0.05, 0) – (0.95, 0);
    ) .

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

    -
    [line width=1. pt, scale=2]
    [fill=gray] (0, 0) circle (0.3cm);
    [fill=black] (0.3, 0) circle (.05cm) node[above right=0.1mm] E_i;
    [fill=white] (1.3, 0) circle (.05cm) node[above=0.5mm] -E;
    [black] (0.35, 0) – node[below] K (1.25, 0);
    
        -
        
    [line width=1. pt, scale=2]
    [fill=gray] (0, 0) circle (0.3cm);
    [fill=white] (0.3, 0) circle (.05cm) node[above right=0.1mm] -E_i;
    [fill=black] (1.3, 0) circle (.05cm) node[above=0.5mm] E;
    [black] (0.35, 0) – node[below] K (1.35, 0);
     .

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. 

4pt
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 (<ref>), 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.

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

    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
     (0.3, 0) circle (.05cm) node[right=1mm] E;
    [black] (0.3,0) – (0.3, 0.05) arc (90:270:0.05) – cycle;
     .

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

    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
    [fill=black] (0.3, 0) circle (.05cm)node[right=1mm] E;
    
        +
        
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
    [fill=white] (0.3, 0) circle (.05cm)node[right=1mm] -E;
    
        =
        -i∫_-∞^+^0ṭ  e^iE t∫^̣3ℓ/(2π)^3[_++(ℓ; t, t) - _–(ℓ; t, t)] .

This combination of graphs vanishes as both propagators are indistinguishable at equal times, i.e. _++(ℓ; t, t) = _–(ℓ; t, t). Notice that this straightforwardly extends to all one-site graphs to any loop order.

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

    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
     (-0.3, 0) circle (.05cm) node[left=1mm] E_1;
    [black] (-0.3,0) – (-0.3, 0.05) arc (90:270:0.05) – cycle;
     (0.3, 0) circle (.05cm) node[right=1mm] E_2;
    [black] (0.3,0) – (0.3, 0.05) arc (90:270:0.05) – cycle;
     .

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

    ∏_i=1^L+1_++(k_i, t_1, t_2) +∏_i=1^L+1_–(k_i, t_1, t_2) = ∏_i=1^L+1_+-(k_i, t_1, t_2) +∏_i=1^L+1_-+(k_i, t_1, t_2) ,

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 (<ref>) 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 ≠ t_2, one further requires the use of the trivial formula Θ(t_1-t_2)Θ(t_2-t_1) = 0, and at equal times t_1 = t_2, one should notice that the Wightman function is real (k_i, t, t) = ^*(k_i, t, t). The special value Θ(0)=1/2 is used. At one-loop order L=1, we therefore have

    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
    [fill=black] (-0.3, 0) circle (.05cm)node[left=1mm] E_1;
    [fill=black] (0.3, 0) circle (.05cm)node[right=1mm] E_2;
    
        +
        
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
    [fill=white] (-0.3, 0) circle (.05cm)node[left=1mm] -E_1;
    [fill=white] (0.3, 0) circle (.05cm)node[right=1mm] -E_2;
     = -i ∫_-∞^+^0ṭ_1  e^iE_1t_1ṭ_2  e^iE_2t_2
       ×∫^̣3 ℓ/(2π)^3[_++(l; t_1, t_2) _++(|k-l|; t_1, t_2) + _–(l; t_1, t_2) _–(|k-l|; t_1, t_2)]  .

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

    -
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
    [fill=black] (-0.3, 0) circle (.05cm)node[left=1mm] E_1;
    [fill=white] (0.3, 0) circle (.05cm)node[right=1mm] -E_2;
    
        -
        
    [line width=1. pt, scale=2]
    [fill=white] (0, 0) circle (0.3cm);
    [fill=white] (-0.3, 0) circle (.05cm)node[left=1mm] -E_1;
    [fill=black] (0.3, 0) circle (.05cm)node[right=1mm] E_2;
     ,

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





 §.§ 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.

4pt
Here, before turning to correlators in a cosmological background which are of primary interest in Sec. <ref>, we first consider correlators in flat space.[In flat-space, we are interested in computing equal-time correlators at a fixed time
t = 0 so that interactions are building up during the time interval -∞ < t < 0. For concreteness, we consider a theory containing a single self-interacting massless scalar field φ. At tree level, the generalisation to (possible different) massive fields is straightforward after setting k→√(k^2+m^2). However, adding a mass to a mediated particle inside a loop is much more complicated.] In this case, the corresponding mode functions are analytic in the time (or energy) domain, and the frequency-space representation of the bulk-to-bulk propagator takes the exact same form as the usual Feynman propagator. Therefore, conventional flat-space methods can be used to compute correlators and reveal their analytic structure. In practice, once the factorised time integrals are evaluated, the remaining frequency integral can be finished with the residue theorem. We illustrate this off-shell factorisation with the simplest of all exchanged diagrams, and explicitly show that the computed correlators consistently satisfy the largest-time equation.




  
Two-site chain. Let us start with the simplest two-chain diagram[Concretely, for the interaction Ł = -g3!φ^3, the s-channel tree-level exchange correlator reads 

    ⟨φ_k_1φ_k_2φ_k_3φ_k_4|'⟩ = g^2/16k_1k_2k_3k_4∑_a, b = ±_ab(E_1=k_12, E_2=k_34, K=s) ,

where s = |k_1+k_2| is the magnitude of the exchange momentum, and k_ij = k_i+k_j.]

    0pt
    [line width=1. pt, scale=2]
     (0, 0) circle (.05cm) node[above=1mm] E_1;
    [black] (0,0) – (0, 0.05) arc (90:270:0.05) – cycle;
     (1, 0) circle (.05cm) node[above=1mm] E_2;
    [black] (1,0) – (1, 0.05) arc (90:270:0.05) – cycle;
    [black] (0.05, 0) – node[above] K (1, 0);
     = ∑_a, b = ±_ab ,  _ab = -ab∫_-∞^a^0 ṭ_1 ∫_-∞^b^0 ṭ_2   e^iaE_1 t_1_ab(K; t_1, t_2) e^ibE_2 t_2 .

In order to compute the time integrals, the iϵ prescription in the in-in contour can be translated to a shift in energy through the formula

    ∫_-∞^±^0 ṭ  e^± izt = ∫_-∞^0 ṭ  e^± i(z ∓ iϵ)t = ∓ i/z∓ iϵ ,

for a real positive energy z>0. This is equivalent to Wick rotating time, as is usually done. However, it gives the correct pole prescription for the later frequency integration when computing the nested contribution. In the limit ϵ→0, the unnested integral is simply given by

    _+- = ∫_-∞^+^0 ṭ_1 e^i(E_1+K)t_1∫_-∞^-^0 ṭ_2 e^-i(E_2+K)t_2 = 1/E_L E_R ,

where E_L ≡ E_1+K and E_R ≡ E_2+K are the sum of energies at the left and right vertices. Note that this contribution has only partial energy poles when E_L, R→0, which emerge from the factorised form. Since the result is purely real, we have _-+ = _+-. Similarly, introducing the frequency-space representation of the time-ordered bulk-to-bulk propagator (<ref>) and performing the time integrals using (<ref>), the time-ordered contribution reads

    _++ = -2K∫_-∞^+∞ω̣/2iπ1/(E_1+ω-iϵ)(ω^2-K^2)_iϵ(E_2-ω-iϵ) .

This representation makes it explicit that both time integrals factorise and that time ordering is encoded in the additional frequency integral. The symmetry under E_1 ↔ E_2 can be restored after changing variables ω→ -ω. The analytic structure of the integrand in the complex ω plane exhibits several interesting features that we show in Fig. <ref>: (i) the pole ω = +E_2-iϵ is the off-shell collinear (folded) limit, (ii) the pole ω = -E_1+iϵ is the off-shell (left) partial energy pole, and (iii) the poles ω = ± (K-iϵ) 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 ω plane. Therefore, selecting for example the physical on-shell condition requires closing the contour in the lower-half complex plane, picking the residues at ω=+K-iϵ and ω=+E_2-iϵ, we obtain

    _++ = 1/E(1/E_L + 1/E_R) ,

where E ≡ 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→0 is made explicit through the identity

    1/(E_1+ω-iϵ)(E_2-ω-iϵ) = 1/E(1/E_1+ω-iϵ + 1/E_2-ω-iϵ) ,

and partial energy poles are reconstructed from the frequency integral. By reality of the result, we have _– = _++. Summing all contributions, the final result reads

    ({E_i}, K) = f({E_i}, K)/E E_L E_R ,   with   f({E_i}, K) = 4(E_1+E_2+K) .


This graph integral is a rational function of the energies and has a singularity when the total energy vanishes E→0, and when partial energies flowing into the left or right vertices add up to zero E_L, R→0, see e.g. <cit.> for more details. The residue at the total energy poles is the corresponding scattering amplitude <cit.>. Indeed, using the identity E_L+E_R=2K in the limit E→0, the graph integral can be written

    lim_E → 0 14K ({E_i}, K) = /E ,

where = 1- is the corresponding s-channel scattering amplitude with =-(k_1^μ + k_2^μ)^2=E_1^2-K^2 being the usual Mandelstam variable. Note that in order to extract the scattering amplitude in Eq. (<ref>), we have imposed the energy conservation condition E_1+E_2=0—valid in the asymptotic past—to write K^2-E_1^2 = E_L E_R. Eventually, the identity (<ref>) makes it explicit that the graph integrals _ab satisfy the largest-time equation

    _++(E_1, E_2, K) + _–(-E_1, -E_2, K) + _+-(E_1, -E_2, K) + _-+(-E_1, E_2, K) = 0 .





  
Three-site chain. We now examine the three-site chain, which presents additional complexities compared to the two-site chain. Notably, this requires evaluating two frequency integrals for the fully nested contribution. The corresponding graph reads

    0pt
    [line width=1. pt, scale=2]
     (0, 0) circle (.05cm) node[above=1mm] E_1;
    [black] (0,0) – (0, 0.05) arc (90:270:0.05) – cycle;
     (1, 0) circle (.05cm) node[above=1mm] E_2;
    [black] (1,0) – (1, 0.05) arc (90:270:0.05) – cycle;
    [black] (0.05, 0) – node[above] K_1 (1, 0);
     (2, 0) circle (.05cm) node[above=1mm] E_3;
    [black] (2,0) – (2, 0.05) arc (90:270:0.05) – cycle;
    [black] (1.05, 0) – node[above] K_2 (2, 0);
       = ∑_a, b, c = ±_abc ,   with
    _abc = i abc∫_-∞^a^0 ṭ_1 ∫_-∞^b^0 ṭ_2 ∫_-∞^c^0 ṭ_3   e^iaE_1t_1   _ab(K_1; t_1, t_2) e^ibE_2t_2_bc(K_2; t_2, t_3) e^icE_3t_3 .

We have in total 2^3=8 integrals to evaluate, among which half of them are related by complex conjugation. Using Eq. (<ref>) and the previous result for evaluating a single frequency integral, the fully factorised and partially nested contributions read

    _+-+ = -1/E_L E_M E_R ,   _++- = -1/E̅_LE_R(1/E_L + 1/E_M) ,   _-++ = -1/E_L E̅_R(1/E_M + 1/E_R) ,

where E_L ≡ E_1+K_1, E_M ≡ E_2+K_1+K_2 and E_R ≡ E_3+K_2 are the sum of energies at the left, middle and right vertices, respectively, and E̅_L ≡ E_1+E_2+K_2 and E̅_R ≡ E_2+E_3+K_1 are the energies of the joined two left and right vertices, respectively. Similar integrals were computed in <cit.> for the corresponding wavefunction coefficient. After using the frequency-space representation for the time-ordered bulk-to-bulk propagators and evaluating the time integrals, the fully nested contribution reads

    _+++ = -(2K_1)(2K_2)    ∫_-∞^+∞ω̣_1/2iπω̣_2/2iπ1/(ω_1^2-K_1^2)_iϵ(ω_2^2-K_2^2)_iϵ
       ×1/(E_1+ω_1-iϵ)(E_2-ω_1+ω_2-iϵ)(E_3-ω_2-iϵ) .

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

    _+++ =    -1/E(1/E_2+E_3+K_1-1/E_2+E_3-K_1)(1/E_3+K_2-1/E_3-K_2) 
       -1/(E_3-K_2)(E_1+E_2+K_2)(1/E_2+K_1+K_2 -1/E_2-K_1+K_2) 
       - 1/(E_1+K_1)(E_2+E_3-K_1)(1/E_3+K_2-1/E_3-K_2) 
       - 1/(E_1+K_1)(E_2-K_1+K_2)(E_3-K_2) .

Notice that it is important to keep track of the correct iϵ 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

    ({E_i}, {K_i}) = f({E_i}, {K_i})/E E_L E̅_L E_M E̅_R E_R , 
    with   f({E_i}, {K_i})    = -8[E̅_LE_M E̅_R + E_M(E_1E̅_L+E_3 E̅_R)+E_1E_3(E̅_L+E̅_R)] .

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

    0    = _+++(E_1, E_2, E_3, K_1, K_2) + _—(-E_1, -E_2, -E_3, K_1, K_2) 
       + _++-(E_1, E_2, -E_3, K_1, K_2) + _–+(-E_1, -E_2, E_3, K_1, K_2) 
       + _-++(-E_1, E_2, E_3, K_1, K_2) + _+–(E_1, -E_2, -E_3, K_1, K_2) 
       + _+-+(E_1, -E_2, E_3, K_1, K_2) + _-+-(-E_1, E_2, -E_3, K_1, K_2) ,

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.[It is possible to map correlators of conformally coupled fields to correlators of massive scalars with half-integer values of ν by applying a set of weight-shifting operators. This procedure follows from the property of the Hankel function that drastically simplifies for ν=n/2 with n integer (n=1 is the conformally coupled case).] The corresponding graph integrals are

    0pt
    [line width=1. pt, scale=2]
     (0, 0) circle (.05cm) node[above=1mm] E_1;
    [black] (0,0) – (0, 0.05) arc (90:270:0.05) – cycle;
     (1, 0) circle (.05cm) node[above=1mm] E_2;
    [black] (1,0) – (1, 0.05) arc (90:270:0.05) – cycle;
    [black] (0.05, 0) – node[above] K (1, 0);
     = ∑_a, b = ±_ab ,  _ab = -ab∫_-∞^a^τ_0τ̣_1/τ_1^2∫_-∞^b^τ_0τ̣_2/τ_2^2  e^iaE_1 τ_1_ab(K; τ_1, τ_2) e^ibE_2 τ_2 ,

where τ_0<0 is a small late-time cutoff. Shifting the energy and neglecting analytic terms in the late-time limit, we use the formula

    ∫_-∞^±^τ_0τ̣/τ e^± i z τ = ∫_-∞^τ_0τ̣/τ e^± i (z∓ i ϵ) τ = log[± i (z∓ iϵ)τ_0] ,

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

    _+- = log(+i E_L τ_0)log(-iE_Rτ_0) .

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

    _++ = 2K ∫_-∞^+∞ω̣/2iπlog[i(E_1+ω-iϵ)]log[i(E_2-ω-iϵ)]/(ω^2-K^2)_iϵ .

This integral, besides being formally divergent, is hard to evaluate because of the branch cuts spanning over ω∈ (-E_1+iϵ-i∞, -E_1+iϵ) and ω∈ (E_2-iϵ+i∞, E_2-iϵ) (choosing the log branch cut to be the principal one). These branch cuts notably cross the real axis. As we will see in Sec. <ref>, the presence of branch cuts, rather than poles, indicates particle production. In this scenario, the integrand can be transformed into a meromorphic function by transitioning to the complex mass plane instead of the complex frequency plane. However, for the special case of conformally coupled fields, as first noticed in <cit.>, this integral can be written as an energy integral over the flat-space result

    _++   = 2K ∫_-∞^+∞ω̣/2iπ1/(ω^2-K^2)_iϵ∫_E_1+ω-iϵ^∞∫_E_2-ω-iϵ^∞x̣ỵ/xy
       = 2K ∫_E_1^∞x̣∫_E_2^∞ỵ∫_-∞^+∞ω̣/2iπ1/(x+ω-iϵ)(ω^2-K^2)_iϵ(y-ω-iϵ) ,

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 <cit.> and the result reads[In order to extract the divergence, one can subtract and add the term

    ∫_-∞^+^τ_0τ̣_1/τ_1τ̣_2/τ_2 e^i(E_Lτ_1+E_Rτ_2) = log(iE_L τ_0)log(i E_R τ_0) = ∫_E_1^∞x̣∫_E_2^∞ỵ 1/(x+K-iϵ)(y+K-iϵ) ,

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

    ∫_0^∞x̣ỵ/(x+y+E)(x+E_L)(y+E_R) .

Sophisticated techniques to evaluate such integrals are detailed in <cit.>.]

    _++   = _2(E-E_L/E) + _2(E-E_R/E) + log(E_L/E)log(E_R/E) - π^2/6
       - log(i E_L τ_0)log(i E_Rτ_0) ,

where _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π for z<0. Using the Euler's identity

    _2(z)+_2(1-z)+log(z)log(1-z)=π^2/6 ,

both nested contributions combine to give

    _++(E_1, E_2, K) + _–(-E_1, -E_2, K) =    -log[i(E_1-K)τ_0]log[i(E_2+K)τ_0]
       -log[i(E_1+K)τ_0]log[i(E_2-K)τ_0] ,

which is exactly minus the combination _+-(E_1, -E_2, K)+_-+(-E_1, E_2, K). Consequently, the graph integrals _ab satisfy the largest-time equation. 







§ 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.

4pt
We are interested in computing the de Sitter tree-level exchange graph[Concretely, considering a conformally coupled scalar field φ coupled to a massive scalar field ϕ by Ł_int/a^3(t) = -g2ϕφ^2 in de Sitter, the full correlator reads

    ⟨φ_k_1φ_k_2φ_k_3φ_k_4|'⟩ = g^2 τ_0^4/16k_1k_2k_3k_4  F(E_1 = k_12, E_2 = k_34, K = s) + t + u  channels ,

where we have introduced a late-time cutoff τ_0 and defined F ≡∑_a, b=± F_ab. Notice that we have defined the integrals F_ab without the factor 2 to ease comparison with the literature.]

    0pt
    [line width=1. pt, scale=2]
     (0, 0) circle (.05cm) node[above=1mm] E_1;
    [black] (0,0) – (0, 0.05) arc (90:270:0.05) – cycle;
     (1, 0) circle (.05cm) node[above=1mm] E_2;
    [black] (1,0) – (1, 0.05) arc (90:270:0.05) – cycle;
    [black] (0.05, 0) – node[above] K (1, 0);
       = ∑_a, b = ± F_ab ,
    with   F_ab = -abK ∫_-∞^a^0 τ̣'/(-τ')^1/2e^iaE_1τ'   ∫_-∞^b^0 τ̣”/(-τ”)^1/2e^ibE_2τ”_ab(K; τ', τ”) .

Recall that _ab are the propagators for the canonically normalised field σ(τ, x) ≡ (-τ)^-3/2ϕ(τ, x), hence the unusual power for the conformal time within the integral. We define the usual dimensionless kinematic variables 

    u ≡K/E_1 ,    v ≡K/E_2 ,

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

    F^(3)(z, μ) = 1√(2) |Γ(12+iμ)|^2   P_iμ-1/2(z) ,

where P_iμ-1/2 is the Legendre function defined in Eq. (<ref>). The non-time-ordered contribution, for which both time integrals factorise, is therefore given by

    F_+- = F^(3)(u^-1, μ)   F^(3)(v^-1, μ) .

We have used the useful formula (<ref>) to perform the time integrals. Notice that the Legendre P function is real, i.e. [P_iμ-1/2(z)]^* = P_iμ-1/2(z) for z≥ 0 and μ real. Therefore, the on-shell three-point function F^(3)(z, μ) is real. As a consequence, we obtain the same expression for F_-+. However, anticipating what follows, the off-shell function F^(3)(z, ν) where the mass parameter ν is analytically continued to the complex plane is not real. 





 §.§ 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

    F̅^(3)(z) = ∫_-∞^+∞ν̣ _ν F^(3)(z, ν)/(ν^2-μ^2)_iϵ ,

with z≥ 1, and F^(3) defined in Eq. (<ref>). The object F̅^(3) is typically proportional to F_++ in the folded limit E_2→ K (or E_1→ K by symmetry), equivalently v→ 1 (or u→ 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, ν) in the complex ν plane, as well as its behaviour at infinity, to appropriately determine how to close the contour. 





4pt
The crucial observation is that the off-shell three-point function F^(3) is a meromorphic function in the complex ν plane. At a fixed kinematic configuration z, the function F^(3)(z, ν) has poles at ν = ± i  (12+n) with integer n≥0, which corresponds to the poles of |Γ(12+iν)|^2 = Γ(12+iν) Γ(12-iν), since the Legendre function P is analytic in the entire complex ν plane. However, the Legendre P representation of F^(3) makes the large-ν 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 (<ref>), yielding

    F^(3)(z, ν) = e^-iπ2/√(2)sinh(πν)[Q_-iν-1/2(z) - Q_+iν-1/2(z)] .

The function Q_± iν - 1/2(z) has purely imaginary poles located at ν = ± i  (12+n) (which come from Γ(12± iμ) in Eq. (<ref>)), and should be closed in the lower (upper) complex plane since Q_± iν -1/2(z) ∼ e^∓ iνξ with ξ >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.[The factor 1/sinh(πν) = |Γ(1+iν)|^2/(πν) also brings additional poles to the integrand. However, they are precisely cancelled by the density of states _ν.] We illustrate the analytic structure of the integrand (<ref>) and the contour prescription in Fig <ref>. Eventually, picking up only the particle production poles, we obtain

    F̅^(3)(z) = F^(3)(z e^+iπ, μ) ,

where we have used the analytic continuation of the Legendre function (<ref>) to re-express the result in the Legendre P basis. As already seen in Sec. <ref>, 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 → ze^+iπ. Actually, we show in App. <ref> that the function F^(3)(z, μ) and its analytic continuation F^(3)(ze^+iπ, μ), viewed as functions of z for fixed μ, 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.





 §.§ Bootstrapping via the spectral representation



We now explicitly compute the integral

    F_++ = - ∫_-∞^+∞ν̣ _ν F^(3)(u^-1, ν) F^(3)(v^-1, ν)/(ν^2-μ^2)_iϵ .

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. (<ref>). 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 ν. The integrand takes the form

    _ν    F^(3)(u^-1, ν) F^(3)(v^-1, ν) = -ν/2πsinh(πν)
       ×[Q_-iν-1/2(u^-1)Q_-iν-1/2(v^-1) - Q_-iν-1/2(u^-1)Q_+iν-1/2(v^-1) + (ν↔-ν)] .

First, we immediately observe that both terms related by the shadow transformation ν↔ -ν contribute equally to the integral. Thus, it suffices to evaluate only one of these terms. Second, we notice that that the term Q_-iν-1/2(u^-1)Q_-iν-1/2(v^-1) has a well-defined large-ν 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ν-1/2(u^-1)Q_+iν-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 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_++^<, > ,

with (we have accounted for the factor 2 coming from the ν↔ -ν symmetry)

    F_++^0    = 1/π∫_-∞^+∞ν̣ ν/sinh(πν)Q_-iν-1/2(u^-1) Q_-iν-1/2(v^-1)/(ν^2-μ^2)_iϵ , 
    
            F_++^<, >   = -1/π∫_-∞^+∞ν̣ ν/sinh(πν)Q_-iν-1/2(u^-1) Q_+iν-1/2(v^-1)/(ν^2-μ^2)_iϵ .

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

    
  * A first tower of poles comes from the factor ν/sinh(πν) = Γ(1+iν) Γ(1-iν)/π. They are located on the entire imaginary axis ν = ± i  (1+n) with n≥ 0, and cannot be evaded by the contour prescription.
    
  * The factor Q_-iν-1/2(u^-1) Q_-iν-1/2(v^-1) also generates an infinite tower of (second-order) poles located on the imaginary axis ν=-i  (12+n) with n≥0. However, at large ν, the integrand behaves as ∼ e^+iν(ξ_u+ξ_v) with ξ_u = cosh^-1(u^-1) ≥ 0 for u^-1≥1 (and similarly for ξ_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 Sec. <ref>. 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ν-1/2(u^-1) Q_+iν-1/2(v^-1) are located at ν=± i  (12+n) with n≥0. At large ν, it behaves as ∼ e^+iν(ξ_u-ξ_v). For |u|<|v|, we have ξ_u-ξ_v≥ 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ϵ prescription.

Let us now collect these various poles in turn. 




  
Particle production residues. The poles coming from the iϵ prescription are the easiest to evaluate. We obtain

    F_++^0 ⊃-i/2sinh^2(πμ)[e^+πμ Q_+iμ-1/2(u^-1)Q_+iμ-1/2(v^-1) + (μ↔ -μ)] .

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

    F_++^< ⊃i/2sinh^2(πμ)[e^+πμ Q_+iμ-1/2(u^-1)Q_-iμ-1/2(v^-1) + (μ↔ -μ)] .

Combining both contributions, and projecting back onto the Legendre P basis using (<ref>) and (<ref>) yields

    F_++^collider = - F^(3)(-u^-1, μ) F^(3)(v^-1, μ) ,

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 Sec. <ref>. The result (<ref>) can be directly found by performing the time integrals after substituting _++→_-+.




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

    F_++^0 ⊃ -2/π ∑_n=0^+∞ (-1)^n  (n+1)/(n+1)^2+μ^2  Q_n+1/2(u^-1) Q_n+1/2(v^-1) ,

where we have used that the residue of Γ(1+iν) at ν = +i(1+n) is (-1)^nn!. The terms in 1/μ^2, characteristic of an EFT expansion, come from 

    1/(ν^2-μ^2)_iϵ-1/(n+1)^2 - μ^2 .

The large n behaviour of the series coefficients is given by

    |(n+1)/(n+1)^2+μ^2  Q_n+1/2(u^-1) Q_n+1/2(v^-1)| ∼1/(n+1)^2+μ^2e^-(n+1)(ξ_u+ξ_v)/√(sinhξ_u sinhξ_v) ,

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_++^< ⊃2/π ∑_n=0^+∞ (-1)^n  (n+1)/(n+1)^2+μ^2  Q_n+1/2(u^-1) Q_-n-3/2(v^-1) .

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

    Q_-n-3/2(v^-1) - Q_n+1/2(v^-1) ∝sin[π(n+1)] .

This effect, which generally occurs in odd spatial dimensions, was previously noted in <cit.> 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 ν= ± i (12+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_++^EFT = u ∑_n=0^+∞ (-1)^n/(n+ 12)^2+μ^2 (u/v)^n  211+n2, 1+n232+nu^2 211-n2, -n212-nv^2 .

The Legendre Q function being not defined for negative integer order, the above expression is found after evaluating the residue of the Γ function entering (<ref>).




  
Full result. The full result is found by summing all contributions, namely the non-time ordered pieces (<ref>) and the time-ordered pieces combining the collider signal (<ref>) with the EFT series (<ref>). 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 <cit.>. Fig. <ref> 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 <cit.>. Second, compared to the bootstrap result of <cit.>, 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 ,

from their explicit forms in (<ref>), (<ref>) and (<ref>). Indeed, both EFT series cancel against each other thanks to the overall factor u, i.e. F_–^EFT(-u, -v) = - F_++^EFT(u, v). Similarly, F_++^collider(u, v) cancels against F_-+(-u, v), and F_–^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. 





 §.§ 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ϵ prescription encoding particle production—, the spectral integral reduces to a contact four-point function generated by the leading bulk interaction φ^4 in a gradient expansion. This simplification is revealed through the generalised Mehler-Fock transformation (Eq. (14.20.14) of <cit.>)

    ∫_-∞^+∞ν̣ _ν  F^(3)(u^-1, ν) F^(3)(v^-1, ν) = uv/u+v .

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→0, i.e. internal soft momentum K→0 while keeping all external momenta hard, and (ii) the partial-energy pole u→ -1 while keeping v fixed.




  
Collapsed limit. In the limit u, v→ 0, we use the asymptotic behaviour of the Legendre P function (<ref>) for large argument to easily obtain

    lim_u, v→0 F_+- = (uv/4)^12+iμ Γ(12+iμ)^2Γ(-iμ)^2/2π .

Taking the limit inside the spectral integral, we obtain

    lim_u, v→0 F_++ = -1/2π(uv/4)^12∫_-∞^+∞ν̣/(ν^2-μ^2)_iϵ(uv/4)^+iνΓ(12+iν)^2 Γ(-iν)/Γ(+iν) .

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

    lim_u, v→0 F_++⊃i/2[e^+πμ(uv/4)^12+iμΓ(12+iμ)^2Γ(-iμ)^2/2π + (μ↔-μ)] .

Notice that two infinite towers of poles arise: one from Γ(12+iν) located at ν = +i  (12+n) that we do not collect as they lie in the upper-half complex plane, and another from Γ(-iν) located at ν=-i n with n≥ 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→0 limit can be directly recovered from its full expression (<ref>). Eventually, keeping only the terms which are non-analytic in uv as u, v→0 (and therefore also non-analytic in K), we obtain

    lim_u, v→0 F = (uv/4)^12+iμ (1+isinhπμ)  Γ(12+iμ)^2Γ(-iμ)^2/2π + c.c. ,

which exactly reproduces the expression found in <cit.>.




  
Partial-energy pole. At fixed v, let us probe the limit u→ -1. The three-point function F^(3)(z, μ) 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→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, μ) ∼1/√(2)lim_z→0∫_z^+∞x̣/x e^-ix(1+u^-1) = 1/√(2)lim_z→0E_1[(1+u^-1)iz] ,

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

    lim_u→-1 F^(3)(u^-1, μ) = -1√(2)log(1+u) .

The non-time-ordered contributions give

    lim_u→-1(F_+- + F_-+) = -√(2)log(1+u) × F^(3)(v^-1, μ) .

Similarly, the time-ordered contribution is found to be

    lim_u→-1 F_++ = 1/√(2)log(1+u) ∫_-∞^+∞ν̣ _ν F^(3)(v^-1, ν)/(ν^2-μ^2)_iϵ = 1/√(2)log(1+u) × F^(3)(-v^-1, μ) ,

where we have used (<ref>). Summing all contributions, we eventually obtain

    lim_u→-1F = -√(2)log(1+u) ×[F^(3)(v^-1, μ) - F^(3)(-v^-1, μ)] .

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. <cit.>. Similarly, the correlator has a singularity in the flat-space limit u→ -v, i.e. E = E_1+E_2→0 or equivalently μ→∞ as seen in Sec. <ref>, with a coefficient that is related to the flat-space scattering amplitude <cit.>.







§ 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.

4pt
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.





  
Acknowledgements. 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). 
This article is distributed under the Creative Commons Attribution International Licence (https://creativecommons.org/licenses/by/4.0/CC-BY 4.0).







§ FEYNMAN IΕ PRESCRIPTION FROM VACUUM WAVE-FUNCTIONAL



In this appendix, we come back to the iϵ prescription defining the in-in generating function (<ref>). As is well known, deforming the time integration contour around the infinite past explicitly breaks unitarity <cit.>. 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ϵ prescription in the Feynman propagator (<ref>).

4pt
Reaching convergence in the infinite past is reminiscent of adiabatically projecting the vacuum of the fully interacting theory |Ω⟩ onto the free one |0⟩. Within the in-in path integral (<ref>), this procedure is achieved with the additional term Ψ[φ_+](t_0)Ψ^*[φ_-](t_0) where Ψ[φ](t) ≡⟨φ(t)|Ω|$⟩ is the vacuum wave-functional that describes the transition amplitude from the vacuum|Ω⟩to a specific field configurationφ(x)at some timet_0<cit.>. Let us therefore compute this object and take the limitt_0 →-∞. The key insight is to use the defining property of the annihilation operatorâ_k|Ω⟩ = 0. After inverting the canonical quantisation of the fieldφ̂_kand its conjugate momentump̂_k, and using the Wronskian conditionu_k p_k^* - u_k^*p_k = i, valid at all time and fixed by the commutation relation (p_kis the mode function associated top̂_k), one obtains

    â_k = 1/i (p_k^* φ̂_k - u_k^* p̂_k) .

Using the field-space representation of the momentum operatorp̂_k →-iδ/δφ_k, withδφ_q/δφ_k = (2π)^3 δ^(3)(q+k), the conditionâ_k|Ω⟩ = 0projected onto⟨φ(t_0)|yields the functional differential equation

    (u_k^* δ/δφ_k - i p_k^* φ_k)Ψ[φ](t_0) = 0 .

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

    Ψ[φ](t_0) = exp(-1/2∫^̣3k/(2π)^3φ_k(t_0)ω_k(t_0)φ_-k(t_0)) ,

whereω_k(t_0) = -ip_k^*/u_k^*, with the mode functions evaluated att_0. The overall normalisation is fixed to unity due to the correct vacuum normalisation⟨Ω|Ω|=⟩ 1. We now take the limitt_0→-∞. Reaching the Bunch-Davies state imposes that mode functions behave as simple plane-wavesu_k(t) ∼e^-iktso thatω_k = kis nothing but the dispersion relation. As for the remaining terms, we introduce the following regularisation scheme

    f(-∞) = ϵ→ 0lim  ϵ∫_-∞^t ṭ' f(t') e^ϵ t' ,

valid for any smooth-enough functionf(t), that can be recovered after integration by parts, to rewrite the contribution from vacuum wave-functionals entering (<ref>) as

    t_0 → -∞lim Ψ[φ_+](t_0) Ψ^*[φ_-](t_0) = ϵ→ 0lim  exp(-ϵ/2∫_-∞^tṭ' ∫^̣3k/(2π)^3 k[φ_k, +(t') φ_-k, +(t') + (+↔ -)]) .

We have sete^ϵt' = 1since we only care about the leading term asϵ→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φ_±written in Fourier spaceℒ[φ_k, ±] = -1/2φ_k, ±(t)(-∂_t^2 - k^2)φ_-k, ±(t). Adding the term (<ref>) amounts to modifying the weight in the exponential of the generating function (<ref>) to

    ±i/2∫^̣3k/(2π)^3∫_-∞^t ṭ'  φ_k, ±(t')[-∂_t'^2 - (k± iϵ)^2]φ_-k, ±(t') ,

where we have used(k±iϵ)^2 ≈k^2 ±2iϵand relabelled2ϵ→ϵ. In the end, performing the converging Gaussian path integral over the fieldsφ_±leads to the correctiϵprescription for the Feynman propagator (<ref>).







§ DETAILS ON MASSIVE FIELDS IN DE SITTER


Consider a real scalar fieldϕ(τ, x)of massm, where-∞<τ≤0is the conformal time, evolving in de Sitter space with three spatial
dimensions. We will focus on fields in the principal series so thatμ≡√(m^2/H^2 - 9/4)is real. In what follows, we setH=1for convenience. Defining the canonically normalised fieldσ(τ, x) ≡(-τ)^-3/2 ϕ(τ, x), its positive-frequency mode function is given by

    u_k(τ, μ) = i√(π)/2  e^-πμ/2  H_iμ^(1)(-kτ) e^-ikτ/√(2kτ) ,

whereH_iμ^(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_iϵin (<ref>) requires analytically continuing the positive-frequency mode function (<ref>) to the lower-halfkτ-complex plane to ensure reaching the asymptotic vacuum. This precisely avoids crossing the branch cut of the Hankel functionH_iμ^(1)(z), that is chosen to lie on the real negative axis, away from the physical regime. Therefore, using the complex conjugate relation (<ref>), the negative-frequency mode can be written

    u_k^*(τ, μ) = -i√(π)/2 e^+πμ/2 H_iμ^(2)(-kτ) .



 §.§ 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μ^(1)(z)]^* = e^+πμH_iμ^(2)(z^*) ,   [H_iμ^(2)(z)]^* = e^-πμH_iμ^(1)(z^*) ,

forz ∈ℂ92{ℝ^-}. Notice that these relations are not verified when the argumentzlies on the branch cut.




  
Symmetries. The Hankel functions obey the following relations

    H_-iμ^(1)(z) = e^-πμH_iμ^(1)(z) ,    H_-iμ^(2)(z) = e^+πμH_iμ^(2)(z) ,

forz ∈ℂ. These identities render manifest an underlying shadowμ↔-μsymmetry for the positive- and negative-frequency mode functions

    u_k(τ, -μ) = u_k(τ, μ) ,    u_k^*(τ, -μ) = u_k^*(τ, μ) .

Upon rotating time, the Hankel functions satisfy

    H_-iμ^(1)(e^+iπ z) = -H_iμ^(2)(z) ,    H_-iμ^(2)(e^-iπz) = -H_iμ^(1)(z) .

The first relation is valid forzin the lower-half complex plane excluding the negative real and imaginary axis, and the 
second relation is valid forzin 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, μ) = u_k^*(e^-iπz, -μ) ,    u_k^*(z, μ) = u_k(e^+iπz, -μ) .

Additional analytic continuations of the Hankel function are given by

    H_iμ^(1)(e^-iπz)    = e^+πμH_iμ^(2)(z) + 2cosh(πμ) H_iμ^(1)(z) , 
     H_iμ^(2)(e^+iπz)    = e^-πμH_iμ^(1)(z) + 2cosh(πμ) H_iμ^(2)(z) .


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

    H_iμ^(1)(z) = e^+πμJ_+iμ(z) - J_-iμ(z)/sinh(πμ) ,   
            H_iμ^(2)(z) = J_-iμ(z) - e^-πμJ_+iμ(z)/sinh(πμ) .

The Bessel function of the first kind satisfies 

    J_iμ^*(z) = J_-iμ(z^*) ,

forz ∈ℂ92{ℝ^-}, and has the following power-law asymptotic form for small argumentz →0J_iμ(z) ∼1/Γ(1+iμ)(z/2)^iμ ,

whereas the Hankel function has the following plane-wave asymptotic form for large argumentz→∞H_iμ^(1)(z) ∼√(2/π z)  e^i(z - π4) + πμ2 .


  
Large-order asymptotic forms. For largeμwithz(≠0) fixed, we have

    J_iμ(z) ∼1/√(2iπμ)(ez/2iμ)^iμ ,    H_iμ^(1)(z) ∼ - H_iμ^(2)(z) ∼ -i √(2/iπμ)(ez/2iμ)^-iμ .


  
Legendre functions. The Legendre functions are equivalent to Bessel (and Hankel) functions with higher transcendentality. Explicitly, they are given by

    P_iμ-1/2(z)    ≡2112+iμ, 12-iμ11-z/2 , 
    
            Q_iμ-1/2(z)    ≡√(π/2) Γ(12+iμ)/Γ(1 + iμ)  2^-iμ z^-12-iμ 213/2+iμ/2, 1/2+iμ21+iμ1/z^2 ,

in terms of the hypergeometric function_2 F_1that is defined by the following formal series 

    21a, bcz = ∑_n=0^∞(a)
        _n (b)_n/(c)_nz^n/n! ,

on the disk|z|<1and by analytic continuation elsewhere, where(a)_n ≡Γ(a+n)/Γ(a)is the Pochhammer symbol. These functions areP_α(z)=  andQ_α(z)=  in Mathematica. The optional argument “3" selects the correct branch cut structure as it is falsely implemented by default.[The Legendre functions P and Q generalise to associated Legendre functions (also called Ferrers function of the first and second kind) P_iμ-1/2(z)→ P_iμ-1/2^1/2-j(z) and Q_iμ-1/2(z)→ Q_iμ-1/2^1/2-j(z), see Chap. 14 of <cit.>, 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.] At fixed argument, the LegendrePfunction is analytic in the entire complexμplane. However, the LegendreQfunction is meromorphic as it has poles coming from the factorΓ(12+iμ).




  
Integral formula. The Legendre functionPis the solution of the following integrals

    -i √(π)/2e^+πν2∫_-∞^+^0 τ̣  e^i kτ/(-τ)^1/2   H_iν^(2)(-sτ)    = e^-iπ4/√(2s) |Γ(12+iν)|^2   P_iν-1/2(k/s) , 
    
            +i √(π)/2e^-πν2∫_-∞^-^0 τ̣  e^-i kτ/(-τ)^1/2  H_iν^(1)(-sτ)    = e^+iπ4/√(2s) |Γ(12+iν)|^2   P_iν-1/2(k/s) .


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

    |Γ(12+iμ)|^2   P_iμ-1/2(z) = e^-iπ2/sinh(πμ)[Q_-iμ-1/2(z) - Q_+iμ-1/2(z)] .

Its analytic continuation is defined by

    |Γ(12+iμ)|^2   P_iμ-1/2(z e^+iπ) = 1/sinh(πμ)[e^+πμQ_+iμ-1/2(z) - e^-πμ Q_-iμ-1/2(z)] .


  
Large-order asymptotic forms. The LegendrePfunction has a simple behaviour atz→0(similar to BesselJfunction), while the LegendreQfunction has a simple behaviour atz→∞(similar to HankelHfunction). Indeed, they are given by

    P_iν-1/2(coshξ)    ∼(ξ/sinhξ)^1/2 , 
    
            Q_iν-1/2(coshξ)    ∼π/2i(2ξ/πνξsinhξ)^1/2  e^-i(νξ - π4) .

These expressions equivalently give the large-order asymptotic forms. The asymptotic behaviour of the LegendrePfunction reads

    P_iν-1/2(z) ∼1/√(π)Γ(-iν)/Γ(12-iν)(1/2z)^12+iν ,

forz→∞.





 §.§ Derivation of the spectral representation


For completeness, we present the proof of the spectral representation of massive scalar field in de Sitter, closely following <cit.>. Let us first consider the case where the timeτ_1is in the future ofτ_2, i.e.-kτ_1<-kτ_2, so that the time-ordered propagator reduces to_++(k; τ_1, τ_2) = u_k(τ_1, μ) u_k^*(τ_2, μ). Expanding the first Hankel functionH_iν^(2)(-kτ_1)in (<ref>) in terms of Bessel functions using (<ref>), using Eq. (<ref>), and the symmetry of the pole prescription with respect toν↔-ν, the integral can be written

    _++(k; τ_1, τ_2) = i/2∫_-∞^+∞ν̣ ν J_iν(-kτ_1) H_iν^(2)(-kτ_2)/(ν^2-μ^2)_iϵ .

The crucial stage is to observe that the numeratorνJ_iν(-kτ_1) H_iν^(2)(-kτ_2)is analytic in the entire complexνplane, so that we choose to close the integration contour in the lower-half complex plane <cit.>. Since the large-νexpansion (<ref>) gives

    ν J_iν(z_1) H_iν^(2)(z_2) ∼e^iνlog(z_1/z_2)/π ,

the arc at infinity does not contribute forz_1<z_2andIm(ν) →-∞. By Cauchy's residue theorem, only the two poles located atν= ±μ-iϵare selected by the closed contour integral, which yields

    _++(k; τ_1, τ_2)    = π/4sinh(πμ)(e^+πμJ_+iμ(-kτ_1) H_+iμ^(2)(-kτ_2) - e^-πμ J_-iμ(-kτ_1) H_-iμ^(2)(-kτ_2))
       = u_k(τ_1, μ) u_k^*(τ_2, μ) ,

where we have again used Eqs. (<ref>) and (<ref>). For the case where the timeτ_2is in the future ofτ_1, i.e.-kτ_2<-kτ_1, one needs to expand the second Hankel functionH_iν^(2)(-kτ_2)in (<ref>) in terms of Bessel functions and then close the contour, picking up the corresponding residues. 





 §.§ Non-analyticity and dispersive integral in the energy domain


In Sec. <ref>, 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 energyF^(3)(z, μ) →F^(3)(ze^+iπ, μ). Similarly, in Sec. <ref> for the four-point function, we have observed that both the factorised contributionF_-+and the collider contribution of the nested channelF_++^colliderare 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 <cit.>.




  
Complex energy domain. Let us define the off-shell three-point function by its integral representation only

    F^(3)(z, μ) ≡ -i √(π)/2  e^πμ2+ iπ4√(K) ∫_-∞^0 τ̣  e^iE τ/(-τ)^1/2   H_iμ^(2)(-Kτ) ,

wherez ≡E/K, withEandKdenoting the external and internal energies, respectively. The physical kinematic region is given byz ≥1. We now want to analytically continue this function in the entire complexzplane, and in particular outside the physical region, without actually evaluating the integral explicitly (the result is given in (<ref>)). First, from the asymptotic form of the Hankel function, we notice that the UV convergence asτ→-∞is controlled by the phase ofE+K, i.e. the sum of the energies entering this three-point function. The integral converges only whenE+Khas a negative imaginary part. As well explained in <cit.>, the integral can be regularised by deforming the contour-∞→-∞^+(E+K)^-1, with-∞^+ ≡-∞(1-iϵ), so that the lower boundτ=-∞is approached from a direction that depends on the argument ofE+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τ∈(0, +∞), as usual. This gives rise to a discontinuity in the energy domain alongE+K ∈(-∞, 0), equivalently alongz ∈(-∞, -1).




  
Discontinuity. We now show that the full functionF^(3)(z, μ)can be recovered from the knowledge of its discontinuity only. By slightly deforming the time integration, the discontinuity of the integral along the branch cutz ∈(-∞, -1)is given by the discontinuity of the integrand itself

    z[F^(3)(z, μ)]    = ϵ→0lim[F^(3)(ze^-iϵ, μ) - F^(3)(ze^+iϵ, μ)] 
       = i √(π)/2  e^πμ2+ iπ4√(K) ∫_∞^0τ̣ τ[ e^iE τ/τ^1/2   H_iμ^(2)(Kτ)] .

Using the known analytic continuation of the Hankel function (<ref>) and (<ref>), the discontinuity of the integrand is given by

    τ[ e^iE τ/τ^1/2   H_iμ^(2)(Kτ)] = -2icosh(πμ)  e^iE τ/τ^1/2   H_iμ^(2)(Kτ) .

Eventually, after changing variables, we obtain

    z[F^(3)(z, μ)] = -2icosh(πμ)   F^(3)(-z, μ) .


  
Dispersive integral. At fixedμ, since the functionF^(3)(z, μ)is analytic everywhere except on the branch cut, we can use Cauchy's integral formula to write

    F^(3)(z, μ) = ∮_ẓ'/2iπF^(3)(z', μ)/z'-z = ∫_-∞^-1ẓ'/2iπ_z'[F^(3)(z', μ)]/z'-z ,

where the integral contourencircleszcounterclockwise. The second equality is found by deforming the contour as shown in Fig. <ref>, as the integral along the large arc vanishes. Using the found discontinuity of the three-point function (<ref>), we finally get

    F^(3)(z, μ) = -2icosh(πμ)  ∫_-∞^-1ẓ'/2iπF^(3)(-z', μ)/z'-z .

Inversely,F^(3)(-z, μ)is obtained fromF^(3)(z, μ)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 <cit.>.






tocsectionReferencesutphys