Cutting Cosmological Correlators

The initial conditions of our universe appear to us in the form of a classical probability distribution that we probe with cosmological observations. In the current leading paradigm, this probability distribution arises from a quantum mechanical wavefunction of the universe. Here we ask what the imprint of quantum mechanics is on the late time observables. We show that the requirement of unitary time evolution, colloquially the conservation of probabilities, fixes the analytic structure of the wavefunction and of all the cosmological correlators it encodes. In particular, we derive in perturbation theory an infinite set of single-cut rules that generalize the Cosmological Optical Theorem and relate a certain discontinuity of any tree-level $n$-point function to that of lower-point functions. Our rules are closely related to, but distinct from the recently derived Cosmological Cutting Rules. They follow from the choice of the Bunch-Davies vacuum and a simple property of the (bulk-to-bulk) propagator and are astoundingly general: we prove that they are valid for fields with a linear dispersion relation and any mass, any integer spin and arbitrary local interactions with any number of derivatives. They also apply to general FLRW spacetimes admitting a Bunch-Davies vacuum, including de Sitter, slow-roll inflation, power-law cosmologies and even resonant oscillations in axion monodromy. We verify the single-cut rules in a number of non-trivial examples, including four massless scalars exchanging a massive scalar, as relevant for cosmological collider physics, four gravitons exchanging a graviton, and a scalar five-point function.


Introduction
We have convincing evidence that the origin of cosmological perturbations in our universe can be traced back to a primordial phase preceding the hot big bang. According to the leading paradigm, it is the quantum mechanical vacuum that is responsible for seeding the small inhomogeneities that later seed the formation of structures. It is therefore natural to ask whether the footprint of quantum mechanics can be found anywhere in the classical probability distribution that we use as initial conditions to make predictions for cosmological surveys. To answer this question, it is convenient to work with the quantum mechanical wavefunction of the universe, from which all late time correlators and the related probability distribution can be extracted.
Recently, it has been realized that unitary time evolution during the primordial universe leads to a particular relation at tree-level in perturbation theory between the quartic and cubic wavefunction coefficients for a single scalar field. This has been named the Cosmological Optical Theorem [1] and can be alternatively understood as a conserved quantity under unitary time evolution [2]. One expects that, just as it is the case for the optical theorem in particle physics, this relation has avatars to each order in perturbation theory for any set of fields. In this paper we prove that this is indeed the case and extend previous results to an infinite set of single-cut rules for a general tree-level diagram. A related but distinct result are the Cosmological Cutting Rules, which apply also to any loop diagram. Those are discussed in a separate paper [3].
In the original derivation of the Cosmological Optical Theorem, one starts from the iconic unitarity relation U U † = 1, where U is the time-evolution operator, and then evaluates this operator identity inside correlators in perturbation theory. This derivation makes the role of unitarity very explicit, but is also lengthy and cumbersome to extend beyond the lowest order in the perturbative expansion. Here we take instead a different approach that makes generalizations and extensions much more straightforward, but which leaves the role of unitarity somewhat implicit. In particular, we introduce a series of single-cut rules that cut any given diagram representing a wavefunction coefficient into the product of two lower-order diagrams. This is qualitatively analogous to Cutkosky's cutting rules in flat spacetime [4] (see also [5,6] for alternative derivations) and to the cutting rules in AdS (see e.g. [7,8]).
The main focus of this work is extending the validity of the Cosmological Optical Theorem beyond the simplest case of a single massless scalar field in de Sitter space, which was studied in [1]. Indeed, we find that precisely the same Cosmological Optical Theorem applies much more generally. In addition to the requirement that coupling constants are real, our derivation relies on two properties: (i ) Hermitian analyticity of the (bulk-to-boundary and bulk-to-bulk) propagators, namely K * (−k * ) = K(k), and (ii ) the simple factorization of the imaginary part of the bulk-to-bulk propagator. The latter is guaranteed by the nature of the problem, namely the computation of the wavefunction on some constant-time hypersurface, and reduces to a very analogous factorization of the Feynman propagators in flat spacetime. The first condition is more interesting. Here we prove that it follows under very general assumptions for any set of fields with linear dispersion relation, E ∝ |k| (in the infinite past), on any flat FLRW spacetime on which we can impose a Bunch-Davies initial state. This includes many spacetimes that are relevant for cosmology, such as de Sitter and inflation (quasi-de Sitter), as well as (accelerating) power-law cosmologies. It also applies for arbitrary speeds of sound and masses that can be approximated as constant on sub-Hubble time scales. The presence of different species of particles is just a matter of bookkeeping: every spatial momentum in the Cosmological Optical Theorem can be expanded to include additional quantum numbers such as charges and flavours. Likewise, spinning (boson) fields can be treated on the same footing as scalars by considering helicity as yet another quantum number (and leaving polarization tensors to the vertices).
When using our single-cut rules the derivation of the Cosmological Optical Theorem becomes very straightforward. However, the relations we obtain are highly non-trivial in the general case. In particular, we demonstrate here that our results hold for the scalar four-point functions generated by the exchange of a massive scalar computed in [9]. Verifying this relation is challenging and requires careful treatment of special functions and their branch points/cuts. We also compute a four-graviton diagram from graviton exchange by invoking a particularly simple interaction that is expected to appear in the Effective Field theory of Inflation [10] (see e.g. [11,12]).
The rest of this paper is organized as follows. In Section 2, after reviewing the perturbative computation of the wavefunction of the universe and discussing the associated propagators, we introduce single-cut rules in the simplest case of a single scalar field without time-derivative interactions. Then in Section 3, we discuss generalizations to time-derivative interactions, general FLRW spacetimes and spinning fields. Here we also show that a linear dispersion relation is a necessary property to satisfy our single-cut rules, in their current formulation. In Section 4 we study in detail several non-trivial examples, in particular the scalar four-point functions from the exchange of a conformally coupled and general massive scalar and the graviton four-point function from the exchange of a graviton withγ 3 interactions. Finally we conclude in Section 5 and leave technical details on massive exchange (Appendix A), resonant mode functions (Appendix B) and the WKB approximation (Appendix C) to the appendices.
Notation and conventions Our conventions for the Fourier transform are (1.1) We will use bold letters to refer to vectors, e.g. k, and non-bold letters for their magnitude k ≡ |k|. We will usually call k the "energy" by analogy with on-shell, massless particles in Minkowski. We expand the wavefunction [42] as 5 and use a prime when we omit the momentum-conserving delta function, ψ n (k 1 , . . . , k n ) ≡ ψ n (k 1 , . . . , k n )(2π) 3 δ 3 k a , (1.3) When discussing ψ 4 , we will use the following variables: These are not independent since p 2 s + p 2 t + p 2 u = 4 a=1 k 2 a . We also sometimes use the notation We denote the power spectrum of a φ (α) as where α refers collectively to a set of quantum numbers. In general we will omit the η 0 dependence and, when talking about a single scalar field, we'll omit the quantum numbers α. When discussing general tree-level diagrams we will use k a with a = 1, . . . , n to denote the n external momenta and p m with m = 1, . . . , I to denote the I internal momenta. Since we work at tree level, all the p m are fixed in terms of the k a by momentum (but not energy) conservation at every vertex. In particular, for every p m there is a subset {k} m of external momenta such that p m = a∈{k}m k a . We denote internal and external energies by k a and p m , respectively, and use them as variables together with additional rotational invariant contractions of the external momenta ψ n = ψ n ( external energies ; internal energies ; contractions ) (1.7) = ψ n (k 1 , . . . , k n ; p 1 , . . . , p I ; k a · k b , k a · (k b × k c ), k a · (k b ), . . . ) (1.8) We define the Hermitian analytic image of a function of the energy k and momenta k to be the complex conjugate of the function evaluated at −k * and −k. So, for example, the Hermitian analytic image of the function F (k, k) is F * (−k * , −k). We will refer to a function that is equal to its Hermitian analytic image as Hermitian analytic. When evaluating an expression for real k, we will often drop the complex conjugate in analytically continued energies, −k * → −k. This must be interpreted with the following convention: negative energies −k should be approached from the lower-half complex plane, namely −k = −k − i with > 0 infinitesimal. 5 Notice that our convention differs by a minus sign from that used in [3], ψ here n = −ψ there n . We apologize for this inconvenience.
In this section, we derive a single-cut rule for tree-level wavefunction coefficients, namely a simple expression for the discontinuity of an n-point wavefunction coefficient ψ n in terms of discontinuities of lower-point coefficients. This result follows from the Hermitian analyticity of the bulk-to-boundary propagator, K * −k = K k and the fact that the imaginary part of the bulk-tobulk propagator G factorizes. The Cosmological Optical Theorem (COT) for 4-point exchange diagrams derived in [1] is found as a special case (see also [2] and parallel developments on the AdS side in [7,8]). We start our presentation considering a single canonical massless or conformally coupled scalar field φ in de Sitter with no time-derivative interactions. Generalizations will be discussed later: time-derivatives are discussed in Section 3.1, general FLRW spacetimes in Section 3.2 and multiple fields with arbitrary spin in Section 3.4.

Diagrammatic representation of the wavefunction coefficients
Our starting point is a brief review of the formalism to compute the wavefunction coefficients in perturbation theory (for more details see e.g. App. A of [1] or [43]).
A given (connected) contribution to ψ n is represented as a (connected) diagram with n external lines, each with one end on the boundary at η 0 = 0 and the other end on some vertex at time η A . Vertices are connected pairwise by I internal lines with momenta p m with m = 1, . . . , I completely fixed (at tree level) in terms of the external momenta by momentum conservation at every vertex. Every vertex may involve spatial derivatives and therefore factors of the momenta contracted in a rotational invariant way and time derivatives. To simplify the presentation, for the time being we assume there are no time-derivative interactions. Later, in Section 3.1, we will relax this assumption and arrive at the same results. To capture spatial derivatives we allow for a vertex function F (k a , η), which contains contractions of the momenta ending on a given vertex with each other or with the (3d) Levi-Civita symbol. Later on, in Section 3.4, we will allow F to also include polarization tensors of spinning fields.
Every external line is associated with a bulk-to-boundary propagator K k (η) of momentum k and "energy" k ≡ |k| and every internal line to a bulk-to-bulk propagator G p (η, η ) with internal energy p. These solve the following differential equations subject to the boundary conditions, where O k (η)φ k denotes the linearized equations of motion obtained by demanding that φ is an extremum of the quadratic action. Notice that because of the boundary conditions K and G depend on η 0 , but we will omit writing this dependence explicitly. Also, we will generally have in mind η 0 → 0 − . Both propagators can be written in terms of the positive and negative frequency mode functions φ ± , which solve the same differential equation as K but have different boundary conditions. Namely φ + asymptotes to the positive-frequency Minkowski mode functions at η → −∞ whilst φ − is defined so that the Wronskian of the two solutions, W (φ + , φ − ), is equal to −i and, for real energies, φ − p (η) = φ + p (η) * . Explicitly this gives where P p is the power spectrum of φ defined in (1.6). When we cut internal lines we restrict the momenta of the cut line to be real and so the restriction of this definition to real p is particularly useful. When p is real we can use the relationship φ + p (η) * = φ − p (η) and so this can equivalently be expressed as Notice the overall factor of i in our definition of G. With this definition 6 every diagram has an overall factor of −i and every vertex has no factor of i. For example the vertex corresponding to the interaction λφ 3 /3! is simply λ.
In this section, we consider only massless and conformally coupled fields, with mode functions given by: where we allowed for an arbitrary speed of sound c s (which for a single scalar can be set to unity and included via dimensional analysis). The corresponding action is given in (3.9). Later in Section 3.2.3 we will generalize this discussion to the case of arbitrary masses. The final step in computing ψ n is to integrate over the conformal time of all vertices η A from −∞(1 − i ) to η 0 → 0, where > 0 is a small real number to be taken to zero at the end of the calculation.

Properties of the propagators
In this section, we discuss some remarkable properties displayed by the propagators K and G, which in unitary theories 7 lead to powerful relations among the wavefunction coefficients ψ n . 6 Another prescription would be to remove the overall i from the definition of G as well as the overall factor of i for every diagram and put an i for every vertex, e.g. iλ for λφ 3 /3!. At tree level, this is equivalent to our prescription because V = I + 1. This alternative definition might be more intuitive for some because we are expanding e iS cl in perturbation theory. 7 More precisely, our technical assumption in the derivation is that the coupling constants are real.
At a general level, this discussion is strongly influenced by the cutting rules in Minkowski (see e.g. [4][5][6]44]), but the details are quite different. Let's start by discussing the bulk-to-bulk propagator G. A first thing to notice is that the first two terms in G as written in (2.7) are precisely the standard Feynman propagator ∆ p , namely (for real p) It takes this slightly unusual form because we are writing ∆ p in Fourier space for the spatial coordinates but in position space for the time coordinate. If we used the Minkowski mode functions in (2.9) and Fourier transformed the two times, η and η , to energy (frequency) space E we would find an energy-conserving delta function and the familiar form (E 2 − |p| 2 ) −1 . It is straightforward to see that 13) so, the difference between the bulk-to-bulk propagator and the Feynman propagator is a solution of the homogeneous equation of motion, which reminds us of the presence of a boundary at η 0 where we want to compute the wavefunction. This term is introduced to cancel out with ∆ p in the limit η → η 0 or η → η 0 , so that the boundary conditions in (2.3) are satisfied. In other words, the last term in (2.13) ensures that the interactions are turned off as we approach the boundary 8 . What complicates the calculation of the cosmological wavefunction is the presence of nested time integrals, which is ultimately due to the lack of time-translation invariance. Any way to circumvent or avoid nested integrals enormously simplifies the calculation. In the following we achieve precisely this by taking advantage of the two following properties: the factorized nature of the imaginary part of the bulk-to-bulk propagator G and the Hermitian analyticity of the bulk-to-boundary propagator K.
Factorizing the bulk-to-bulk propagator The "boundary" term of G in (2.13) is already promising because the η and η dependence are factorized and so they cannot give rise to nested integrals. The difficulty sits in the part involving the Feynman propagator. There we can borrow a simple trick from Minkowski spacetime. The key observation is that the imaginary part of G is factorized, where we assumed p ∈ R. The hardest part is to put this observation to good use by finding observable quantities that are computed in terms of Im G as opposed to the full G. This problem was solved in [1], albeit in a different language, and relies on a second key observation. 8 It may happen that the interactions diverge at η0 = 0 faster than G vanishes and the result is IR divergent.
In these cases, we have to evaluate the wavefunction at η0 = 0, where the interactions are finite and G vanishes as (η → η0).
Hermitian analyticity of the bulk-to-boundary propagator Notice that for all the explicit examples of mode functions given in (2.9)-(2.11), the following relation, which we call 9 "Hermitian analyticity", is satisfied where we have allowed for a complex "energy" k to account for the fact that a negative value of k should be thought of as belonging to the analytic continuation of K from real momenta and positive k. To remain compatible with the choice of a Bunch-Davies vacuum, real and negative values of k should always be approached from the lower-half complex plane, k ∈ C − . We will now proceed assuming that this relation holds. Later in Section 3.2, we will show that (2.15) follows very generally from the choice of the Bunch-Davies vacuum. Our general proof will also involve an additional weak technical assumption about the linearized equations of motion (which is satisfied by all the models in the literature for which we have tested it). The property (2.15) implies the following relation for the analytic continuation of the bulk-to-bulk propagator This can be seen straightforwardly for both massless and conformally-coupled scalar fields for which the negative energy solutions, φ − k (η), obey the same Hermitian analyticity condition as the positive energy ones. For massive fields, this property still holds, but one has to be careful when analytically continuing the mode functions to complex energies. We leave this technical point for Section 3.2.3.

Single-cut rules
The single-cut rules for the wavefunction coefficients are now easily derived from the two above properties, (2.14) and (2.15). Consider a general tree-level diagram representing a perturbative contribution to ψ n . According to the rules reviewed in the previous subsections, this diagram translates into the following expression Here V is the number of vertices, I = V − 1 the number of internal lines 10 , k a and p m are the external and internal momenta respectively, {k} denote additional non-energy variables obtained from rotation invariant contractions of the external momenta, and we left the time dependence implicit. Now consider choosing as variables for ψ n the norms ("energies") of the external momenta k a , of all the internal momenta p m that appear in the graph, and finally of any additional scalar products k a · k b that might be needed to have a complete set of variables. This should 9 This nomenclature echos that used in the study of amplitudes where the 2-to-2 amplitude enjoys the following property A * (s * ) = A(s), with s the Mandelstam variable. 10 Here we work exclusively at tree level. The cutting rules for loops are derived in [3].
always be possible according to the following counting. Using the graph relation V A n A = 2I + n , (2.18) where n A is the valency of each vertex (the number of internal plus external lines ending on it, matching the number of fields in the corresponding vertex). Restricting to cubic or higher interactions 11 , n A ≥ 3, we find the upper bound I ≤ n − 3. Therefore, the n external plus I internal energies account for variables. For any n ≥ 3, this is always less than or equal to the 3n − 6 variables that are needed to describe a general ψ n (accounting for rotation and translation invariance). We want to consider the following analytic continuation of ψ * n in all external and internal energies except for them-th internal energy, which we denote by S ≡ pm, combined with a parity inversion of all momenta 12 , k → −k, 20) where to simplify the notation we drop the complex conjugate on −k * and state the general convention that all the negative real energies are approached from the lower-half complex plane [1]. Notice that by momentum conservation one also has p m → −p m so all the momenta attached to any vertex get a minus sign. Let's first see what happens to the vertices F A . In real space, the vertex must involve contractions of spatial derivatives either with δ ij (parity even) or with the totally anti-symmetric Levi-Civita symbol ijk (parity odd). In both cases F * (−k, −p) = F (k, p) because the Fourier transform gives ∂ i → ik i and so the minus sign in k cancels with the minus sign from i * = −i.
Using the invariance of the vertices and the Hermitian analyticity of K and G, (2.15) and 11 In particular here we are neglecting "quadratic interactions" which can arise from mixing of fields at linear order. Although such interactions play important roles in many multifield scenarios (e.g. see [45,46]), they bring about a technical difficulty for us: some internal and external energies in their presence become degenerate. We need the distinction between external and internal energies when flipping the sign of some while keeping unchanged the rest. Nevertheless, our results might be naturally generalized to quadratic interactions too if we think of each one of them as the soft limit of a cubic vertex that involves an auxiliary soft conformally coupled field with an independent (external) energy [47]. 12 In [1] only parity invariant scalar theories were considered in which case this parity has no effect and can be neglected. However, both in the presence of parity breaking interactions and spinning fields, this is an essential ingredient to satisfy the COT. This can be easily verified on the parity breaking quartic contact interaction of four different scalars, φ1∂iφ2∂jφ3∂ k φ4 ijk . (2.21) Now we can use the factorization of Im G to re-write the above V nested integrals as a product of two separated sets of (fewer) nested integrals. Then we recognize that, denoting the n-th external leg of a diagram by S and the remaining n − 1 legs by k a , we have where {p} denotes the set of all internal energies. The linear combination in the first line above and in (2.21) is reminiscent of the calculation of discontinuities. Since we will encounter this combination many times in this work, it is convenience to introduce the following notation where f is a generic function. Notice that the spatial momenta always get a minus in the second term, while only the energies k a that do not appear in the argument of Disc are analytically continued to −k a . In other words, the energies appearing in the arguments of Disc remain untouched and are simple spectators.
With this notation, we can re-write the left-hand side of (2.21) using the factorization (2.14) of G and (2.22) into the single-cut rule This relationship is shown diagrammatically in Fig. 1 where we demonstrate the interpretation of the right hand side as a cut. Finally, let's compare these single-cut rules to the recently derived Cosmological Cutting Rules [3], (2.25) First, while both set of rules deal with the discontinuities of wavefunction coefficients, here we only consider single-cuts and so the right-hand side of our expressions contain the product of only two discontinuities. Conversely, the Cosmological Cutting Rules require the sum over all possible ways of cutting internal lines, which in general leads to the product of many discontinuities. In particular, in single-cut rules you get to choose where you want to perform the cut, while there you don't. This is particularly useful in some cases, for example in the derivation of the consequences of manifest locality through the Manifestly Local Test [41]. Second, in the Cosmological Cutting Rules we never analytically continue internal energies {p}. This means that there is no need to arrange so that the variables {p} appear explicitly in the argument of ψ n . In contrast here we need to access each internal energy independently: the uncut internal energies are analytically continued while the cut energy is not. One consequence of this is that for single-cut rules one generally needs to choose different sets of variables for different channels. Third, because of the above considerations it seems challenging to extend our single-cut rules to loop diagrams: in that case there are internal energies that are integrated over and it is not clear how one would analytically continue them by manipulating the kinematical variables.

Generalisations
For the sake of clarity, the derivation in the previous section was written for the simplest case of a massless or conformally coupled scalar field in a de Sitter background with no time-derivative interactions and Bunch-Davies initial conditions. In this section we show that those results can be greatly generalised to general FLRW spacetimes that admit a Bunch-Davies vacuum for any fields with a linear dispersion relation and arbitrary mass, integer spin and speed of sound.

Time Derivatives
When the interaction is allowed to involve time derivatives this introduces derivatives on the Green's function that potentially alter the analysis as it is no longer immediately possible to exploit the result from (2.14). Instead it is necessary to understand the behaviour of terms like This generalisation is simpler than it might seem as the energies are the only complex variables and therefore complex conjugation commutes with the time derivatives and so derivatives of K and G will remain Hermitian analytic, Likewise, we can explore the imaginary part of G for real p, From this it is apparent that we can cut lines involving time derivatives in exactly the same way as non-derivative interactions, except each of the diagrams must include the derivatives previously associated with the bulk-to-bulk propagator on the the external lines that are introduced. We can further clarify this by looking at the full wavefunction coefficient, Where we have introduced the notation whilst suppressing the η dependence. This has an imaginary discontinuity arising from cutting the internal line with momentum S given by This is the same as the expression with non time-derivative interactions except the propagators are now allowed to have derivatives. Therefore, the single-cut rules derived in the previous section apply to any derivative interactions as well. The generalization to time derivatives for spinning fields proceeds similarly as the modefunctions remain Hermitian analtyic and the vertex terms are time independent, see Section 3.4 for more details.

General cosmological backgrounds
We know that our single-cut rules cannot be completely generic because when we allow for non Bunch-Davies initial conditions, even within the de Sitter case, our derivation breaks down as the pivotal condition K k (η) = K * −k * (η) is in general no longer valid. Since there are examples in the literature, see e.g. [48,49], of backgrounds that excite negative frequency modes even when they start with Bunch-Davies initial conditions, one might wonder under what conditions our single-cut rules apply to general FLRW spacetimes.
In this section, we demonstrate that, for a fairly generic FLRW background spacetime, the imposition of Bunch-Davies initial conditions is sufficient to ensure that the propagators are Hermitian analytic. We then also present the mode functions for both a massive field in de Sitter and a massless field in an alternative background, where the background excites negative frequency modes from a Bunch-Davies vacuum, and show that the propagators behave in the same way as for massless and conformally coupled fields. Throughout we focus on linear dispersion relations of the form E ∝ |k|, which are generic in relativistic theories to leading order in derivatives (two time or space derivatives in the quadratic action). Conversely, we stress that we don't require interactions to be Lorentz invariant.

Bulk-to-Boundary Propagator
We will consider here scalar fields with a quadratic action of the form This is the most general quadratic action of a real scalar field to leading (quadratic) order in derivatives. We will discuss higher derivatives and more general dispersion relations in Section 3.3. The mode functions φ(k, η) satisfy a second order differential equation of the form We will assume that p and q are analytic functions of k and η over the domain of η. The ordinary differential equation (3.10) has two linearly independent solutions that we call φ ± k (η) from which Equivalently, in words, these two functions are linearly dependent. It is well known [50] that two analytic functions are linearly dependent if their Wronskian, namely vanishes everywhere. Furthermore, if two functions both satisfy the same differential equation of the form in (3.10) and their Wronskian is zero at some point η i then, because the Wronskian is given by it must vanish everywhere (see e.g. [51]) by virtue of the assumption that p and q are analytic on this domain. To fix φ + k (η) we specify Bunch-Davies initial conditions, This solution assumes that c s kη diverges in the infinite past, By imposing this initial condition we have assumed that (3.15) becomes a solution to our differential equation in the infinite past. To understand when this occurs we rewrite the differential equation as We can first of all see that the Bunch-Davies initial condition has no dependence on m or a (except through the prefactor of a −1 ) and so the last two terms multiplying aφ k must be negligible compared to c 2 s k 2 in this limit, We further need that c s is approximately constant. To quantify this condition we insert this asymptotic solution into the differential equation which gives For this to be an asymptotic solution, we generically require that each of the terms in the bracket vanishes individually and so Imposing this initial condition also restricts the asymptotic form of Therefore, provided φ + k (η) and φ + −k * (η) * satisfy the same differential equation, they are guaranteed to be linearly dependent. To determine what differential equation φ + −k * (η) * satisfies, we can take the complex conjugate of (3.10) and replace k → −k * everywhere to give This coincides with (3.10) if p and q are Hermitian analytic. Therefore, the bulk-to-boundary propagator is Hermitian analytic if: • We impose Bunch-Davies initial conditions, which requires that • a a , a(η), m(η) and c s (η) are real, analytic functions in the domain η ∈ (−∞, 0]. Where we have expressed the conditions on p and q in terms of the functions that appear in the quadratic action. One can also show this by making a WKB approximation (see e.g. [51]) of the mode function in a general flat FLRW spacetime, where we can see generically that φ ± = e ikσ ± (k,η) /a, where σ is Hermitian analytic. The details are included in Appendix C.

Bulk-to-Bulk Propagator
In order to extend our results to diagrams with more than one internal line we also need to prove that for a generic background To do this we need to know the Hermitian analytic properties of φ − k (η) in addition to those of φ + k (η), which we have already established. To determine these, consider that φ ± k (η) are defined to have a Wronskian We can treat this as a differential equation in φ − k (η), which can be formally solved: If the above assumptions are valid we can then exploit the Hermitian analytic properties of φ + k (η) to give where A was defined in (3.12). The bulk-to-bulk propagator is and its Hermitian analytic image is Using the relationships in (3.12) and (3.31) we find this is equal to G p (η, η ) and so the bulk-tobulk propagator is Hermitian analytic whenever the bulk-to-boundary propagator is.

Massive Fields in de Sitter
As a concrete example of these results consider the case of a single scalar field of mass m with a constant speed of sound, which we set to 1, in de Sitter. With these assumptions p(k, η) = − 2 η and q(k, η) = m 2 + k 2 which are both analytic functions for all time and, furthermore, are both Hermitian analytic by inspection. Hence in this theory K k (η) = K * −k * (η). This was shown in [1]. We can also see that the bulk-to-bulk propagator is Hermitian analytic by introducing φ ± k for a massive scalar field, When ν is real we will recover the original Hankel functions but when ν is imaginary we pick up a minus sign, this cancels with the sign change in the exponential factor as and so we can just drop the complex conjugates on ν. We wish to express this in terms of the original mode functions and so we replace kη = e iπ (−kη) where this particular choice of argument for −1 is enforced by the fact that Im(k) < 0. We then use the analytic continuation from Section 10.11 of [52], The Hermitian analytic image of the Green's function is Therefore, the Hermitian analyticity of both propagators is manifest in this theory.

Resonant Non-Gaussianity
One might be worried that due to the evolution in the bulk, negative modes may be excited and ruin the Hermitian analyticity of the propagators, even if we start with Bunch-Davies initial conditions. An example to consider is the resonant production of particles in axion monodromy inflation [53,54] and the associated non-Gaussianity [55][56][57]. Here the inflaton potential is modulated by non-perturbative effects which in turn induce small oscillations in the background. Solving for the slow-roll parameters gives where (see [56] for more details) Putting these slow-roll parameters into the Mukhanov-Sasaki equation gives us the following coefficients for the differential equations: Notice that both p(k, η) and q(k, η) are Hermitian analytic (the conditions on the behaviour of the scale factor follow from the fact that this inflationary spacetime is close to de Sitter and m = 0, c s = 1 everywhere). We therefore expect to find that K k (η) is Hermitian analytic. This is not immediately apparent from the results in the paper, however to see that it is, in fact, true we compute the bulk-to-boundary propagator from the modefunctions given in (2.21) of [56] 13 The Hermitian analytic image of this propagator is given by .
Using the definition of c k from [56], we find that, in order for c k to be Hermitian analytic, it is necessary to keep a term that is suppressed by a factor of e is assumed small. Therefore, we have recalculated c k in Appendix B without making any assumptions on the size of this quantity and we find that it remains true that c k is Hermitian analytic. Furthermore, we showed earlier in this section that the Hermitian analyticity of the bulk-to-bulk propagator follows from that of the bulk-toboundary propagator, therefore, even though this background excites negative energy modes the single-cut rules and the Cosmological Optical Theorem still apply exactly as expected.

Non-linear dispersion relations: a non Hermitian-analytic example
So far we showed that Hermitian analyticity is very general for linear dispersion relations, as derived from the leading-derivative quadratic action in (3.9). Here we want to point out that this is not the case for general non-linear dispersion relations. As a toy model, consider the following quadratic action for a real scalar field 55) 13 In this paper they define the modefunctions, which they label as R k , to be where c n is some real, possibly time-dependent parameter, n is a positive integer and for concreteness we take a to be the scale factor in de Sitter 14 . This action has equations of motion that are second order in time, but involves an arbitrary number, 2n, of spatial derivatives. For n = 1, this theory reduces to what we studied previously, (3.9), while for n = 2, it reduces to the ghost condensate studied in [58]. For any n = 1 we have a non-linear dispersion relation E 2 = c 2 n k 2n . The equation of motion for φ is which satisfies the requirement that p and q are Hermitian analytic for real a and c 2 . The positive energy initial conditions must now be implemented as (3.57) To ensure that this is finite in the infinite past we require that Im(k n ) < 0. If n is even, then the initial condition for the Hermitian analytic image of φ + k is We see that (3.57) and (3.58) are genuinely different solutions, and so the Wronskian will not vanish in the infinite past. Therefore these functions are linearly independent and the bulkto-boundary propagator will not be Hermitian analytic. This can also be seen explicitly by considering the concrete example of n = 2, namely the ghost condensate. The mode function and bulk-to-boundary propagator are found to be (3.60) We can then compute the Hermitian analytic image of K: This means that the single-cut rules as written in this paper to do not apply to this case. Because the dispersion relation for the ghost condensate is still a simple monomial, E ∼ k 2 , we can find an ad hoc modification of our single-cut rules that does work. In particular, consider the transformation k →k = ±ik * . Using this to define a modified Hermitian analytic image, we find Therefore, in this a theory, we expect it to be possible to derive similar results to those presented in the rest of the paper, with the replacement −k * →k. We don't pursue this further here.
Instead, we notice that, for more general, non-monomial dispersion relations, even this modified Hermitian analyticity does not exist. For example consider the theory for constant, real c s , c 2 . In this case we require The higher derivative term will come to dominate in the infinite past, so, the positive energy initial condition is lim This is equal to the same limit of φ + will ensure that K * k (η) = K k (η). (3.70) The problem with (3.69) is that nowk depends on time. If we defined a Disc with this Hermitian analytic image, the Disc would not commute with the time integrals appearing in the Feynman rules for the calculation of the wavefunction and our derivation would not work. We don't pursue this further here, but notice that because of this it seems unlikely that single-cut or general cutting rules can be derived for general non-linear dispersion relations.

Spinning fields
In this section we discuss the generalization of our single-cut rules to integer spin fields. For concreteness we focus on the very general class of free theories for such fields that was developed in [59] where Φ i 1 ...is is a totally-symmetric, traceless tensor with only spatial indices, i 1 = 1, 2, 3. This theory arises in generic models of inflation where the background of the inflaton selects a preferred time foliation of spacetime into spatial hypersurfaces. The above expression can be written in a covariant way by using the Goldstone boson π of time translations to upgrade the spatial tensor Φ i 1 ...is to a covariant spacetime tensor. The coupling of Φ i 1 ...is to π is also dictated by this constructions but we will not need this here. Notice that Φ i 1 ...is has (2s + 1) components, which each create states ("particles") with helicities 0, ±1, . . . , ±s, respectively.

Hermitian anaylicity of the propagators
The equation of motion for the field Φ i 1 ...is is given by: The field can be separated into two parts: ..is is the transverse part of the field, obeying and Φ R i 1 ...is is the remainder. It is straightforward to see that Φ T has 2 degrees of freedom and represents the components with helicity ±s, while Φ R has 2s−1 components with lower helicities.
For Φ T , the penultimate term in (3.72) vanishes, and the equation of motion becomes: This equation is in the same form as (3.10), therefore we can directly apply the analysis in Section 3.2 to show that the propagators of Φ T are Hermitian analytic. For Φ R we can take the divergence of (3.72), which gives us: Once again the equation is in the same form as (3.10), but with c 2 s replaced with c 2 s + δc 2 s . We can again directly apply the analysis in Section 3.2. Working in Fourier space, this tells us that ik j Φ R j...is is Hermitian analytic. Since ik is Hermitian analytic, and ik j Φ R j...is has exactly 2s − 1 degrees of freedom, we deduce that the propagators of Φ R are also Hermitian analytic. We conclude that the propagator of the full field Φ must be Hermitian analytic, which establishes the crucial property of our derivation of single-cut rules for free fields of any integer spin (in the spontaneously boost-breaking theories of [59]).

Helicity basis and the diagonalization of propagators
For practical calculations, we would like to work in a basis where the propagators have a simple form. This can be achieved by looking at the helicity basis of the fields. These are irreps of ISO(3), the isometry group of a flat FLRW spacetime. As we show below, fields of different helicities decouple from each other and the corresponding propagators in this basis become diagonal.
The only non-diagonal term in the action is We are guaranteed to be able to diagonalise this term because it is real and symmetric in i's and j's. To understand this diagonalisation procedure, let's start by looking at the vector case, s = 1, for which the tensor equation can be understood as a matrix multiplication, and so is diagonalised by finding the eigenvalues, λ h , and eigenvectors, h , of M , The eigenvalues of M are We define the eigenvectors so that they satisfy the inversion relationship h (−k) = h (k) * 15 , for some normal vectorn that is perpendicular to k. These eigenvectors are orthogonal to each other, where C h (k 2 ) is a polynomial in k 2 (so is guaranteed to be Hermitian analytic) that comes from the normalisation of the eigenvectors. We can therefore express M and the identity in terms of these eigenvectors as We can then see that in the so called "helicity basis", the action diagonalises, 15 Naively it may appear that this is sufficient to ensure that the Disc commutes with the vertex contributions.
However, as is discussed in Section 3.4.3, the explicit dependence of ± on the energies will ruin this relationship, because not all energies are analytically continued.
This also makes it clear why we imposed that h (−k) = h (k) * as it ensures that Now that we have our eigenvectors for the spin-1 case this procedure can be generalised to arbitrary spin. To keep the symmetries of our field manifest we define a symmetric, traceless basis containing 2s + 1 tensors which are constructed from the symmetrised direct product of s copies of the vector h , . . .
The tracelessness of these terms relies on the relationship + = − * which ensures that any contractions like ± i ± i vanish by orthogonality. It can be shown that these tensors inherit the orthogonality of the vectors so (3.95) where C h is given, as for the vector case, by The action is therefore exactly that given in (3.87) but with h running from −s to s. As this is diagonal, it gives a separate differential equation for each helicity mode 16 , which ensures that the propagators are diagonal in this basis, Here, K h k and G h p are constructed from the positive energy modefunctions that satisfy (3.97) subject to the Bunch-Davies initial condition (3.15) with the substitution This is k independent because λ h ∝ k 2 for all h. The proof that these propagators are Hermitian analytic therefore follows similarly to the scalar case.

Interaction vertices
In order to derive the single-cut rules, we need the interaction vertices to be Hermitian analytic. This is indeed the case if we follow the prescription that the all polarization tensors commute with the Disc operation. In the helicity basis, a generic interaction vertex has the following form: The interaction vertex is constructed by taking various contractions of ha i 1 ...is (k a ) and ik a (the latter comes from a spatial derivative). Clearly ik is Hermitian analytic since (i(−k)) * = ik. Next, we look at the Hermitian analytic image h i 1 ...is (−k) * . Naively one would like to use the property of polarization tensors h i 1 ...is (−k) * = h i 1 ...is (k), but this is subtle because the Disc actually analytically continues some energies (all uncut lines) and not others. To avoid confusion, we provide a clear cut prescription for which our single cut rules are valid: all polarization tensors h i 1 ...is should be factorized outside all the Disc's. In other words, we should substitute h i 1 ...is (−k) * in the second term inside each Disc with h i 1 ...is (k) 17 , which is precisely what appears in the first term. With this prescription, we conclude that the vertex function is Hermitian analytic.
Due to the form of the propagator, polarization tensors associated with bulk-to-bulk propagators must also come with a sum over helicities. Aside from this, the proof of the single-cut rules proceeds in the same way as the scalar case. Therefore, we have the following single-cut rules: Here P h Φ is the power spectrum of the exchanged field, where C h is defined in (3.96).

Explicit examples: general relativity and massive gravity
As a demonstration of the various general properties discussed above, let us look at explicit examples involving massive gravity and general relativity. The simplest case is that of general relativity. In this theory the (massless) graviton has the same (positive energy) mode functions as a scalar field, (2.10), 104) 17 To see how this is compatible with our explicit expression for ± note that we are free to specify how to send k → −k under the Disc. This is distinct from in the inversion relationship where we must rotate k and its normal vectors with it. We therefore choose to invert k by reflecting it in a plane perpendicular to itself which leaves any vectors perpendicular to it unchanged. This leads to the desired result for all uncut lines whilst for the cut line the helicities are reversed. This is not an issue as all helicities are summed over for internal lines and there is a symmetry between the plus and minus helicities that ensures their propagators are the same.
where now h = ±2 since the lower-helicty modes are removed by diff invariance. As we have seen for the scalar field, the propagators corresponding to this mode function must be Hermitian analytic.
As an example of an interaction, consider the cubic graviton interaction induced by the spatial Ricci scalar R (3) [25]: The two spatial derivatives are Hermitian analytic thanks to factor of i 2 in front, and the polarization tensors can be taken to obey Hermitian analyticity because of the prescription outlined in Section 3.4.3. As a more interesting example, we also look at the propagators in a theory of massive gravity (see [43] for more details): As shown in the paper, the mode function for the helicity mode +2 and +1 are found to be: We immediately notice that the +2 helicity has a Hermitian analytic propagator. For the +1 helicity mode, we use the recurrence relation of Hankel function to rearrange the mode function as: from which we see that this mode function will give rise to a Hermitian analytic propagator.

Explicit examples
In this section we show how the single-cut rules are satisfied is a few non-trivial examples. Due to the similarity of the results to those discussed in [1] we will not reproduce the checks performed there but instead discuss more complex cases. First, we consider the four point function for 4 external conformally coupled scalars for non-derivative cubic interactions with both a conformally coupled and a massive exchanged field. Then we discuss a novel case involving 4 gravitons with a graviton exchange. Finally, we explore the 5-point function with derivative interactions to demonstrate how our results generalise for higher point derivative diagrams.

Conformally coupled exchange
We consider the case of a conformally coupled field with a cubic polynomial interaction λφ 3 , as discussed in [13] (see also [60]). There the in-in correlator was computed, but here we are interested in the wavefunction coefficient so we redo the calculation. It is necessary to keep the dependence on the time at which the future boundary is taken because the propagators diverge as this is taken to zero. For this interaction the wavefunction coefficients are given by where the propagators are defined to be For simplicity we have set H = 1, since it can be recovered in the final result using dimensional analysis. We are specifically looking to relate the s channel of ψ 4 , shown in Fig. 2, to ψ 3 and so we focus on where k ij = k i + k j . This can be approached as in [13] by rewriting the divergent 1 η terms as integrals over momentum. Furthermore, as the integrals are finite in the limit η 0 → 0 we find .
This integral can be performed and the resulting logarithms can be combined under the assumption that all energies have negative imaginary parts to give which has been simplified using the fact that the dilogarithm satisfies We wish to compare this to the same limit of the product of the discontinuities in the two three point functions that make up the s channel diagram, From this it is clear that, in this limit, Disc ps iψ s 4 = −iP ps Disc ps iψ 3 (k 1 , k 2 , p s ) Disc ps iψ 3 (k 3 , k 4 , p s ) . (4.14) Therefore, our results hold for this interaction involving conformally coupled fields.

Massive exchange
It has been observed that the four point function of a conformally coupled scalar field serves as a seed solution from which various correlators can be obtained [9,13,35,39,61]. Therefore, it will be useful to demonstrate that the single-cut rule is satisfied by this primary building block of the bootstrap program 18 . In this section, we check the single-cut rule for the 4pt exchange diagram with an intermediate heavy scalar field. In Appendix A we present the case for an arbitrary mass particle as well as some additional details. The three point function arising from the ϕ 2 σ interaction is given by [9,13,62] ψ ϕϕσ (u, 18 We thank Sébastien Renaux-Petel for stimulating discussions about this section.
where we have set H = 1, defined u = ps k 12 and introduced the associated Legendre Function of the first kind, see Section 14 of [52], which is defined in terms of hypergeometric functions as It is also helpful to introduce the power spectrum for the massive field which we will keep implicit as P σ ps = σ + ps (η 0 )σ − ps (η 0 ). (4.17) and we note that p s will be considered real throughout this section in accordance with the procedure for cutting internal lines. The four point correlator for this interaction was computed in [9], however, the corresponding wavefunction coefficient has not been calculated in the literature as far as we are aware. Fortunately, the techniques used in [9] can be naturally extended to the wavefunction of the universe coefficients which are defined by The Green's function here obeys exactly the same differential equation as the Green's functions for the correlators, up to some conventional factors of i, we therefore define F = − psη 4 0 4λ 2 ψ 4 so that it satisfies the differential equation studied in [9], where v = ps k 34 . The ansatz in [9] uses, as a basis, the functions However, the inclusion of ν in the prefactor of these functions requires coefficients which, when we analytically continue in u, lead to complications in the analysis. Furthermore, the Hermitian analytic image of these functions depends on the mass of the particles through the behaviour of the complex conjugate of ν. Therefore, we wish to use an alternative basis. The basis we chose isF Where P and Q are the associated Legendre functions of the first and second type 19 . We use the Legendre functions rather than generic hypergeometric functions because it simplifies the notation and unifies the two hypergeometric functions that appear in the three and four point correlators. We could have chosen any two linearly independent combinations of these solutions but these particular ones are chosen because • The branch cut in P ν (z) is along z ∈ (−∞, −1], the branch cut in Q ν (z) is along z ∈ (−∞, 1] so neither of these solutions have a branch cut with u in its physical range u ∈ (0, 1) and one has no branch cut for u > 0.
• The Wronskian for these two functions is the same as for F ± (u) and so the matching condition that arises from the particular integral will take the same form in terms ofF ± (u) as in [9] therefore, the most general solution that is symmetric under exchange of u, v is 24) and c mn are real coefficients that are already completely fixed by the inhomogeneous part of the differential equation.
• The Hermitian analytic image of these functions does not depend on whether or not ν is real, • F + is finite in the limit u → 1, and so in this limit F (u, v) is given by, We want this to be finite and so we require that β − = 0, β 0 = −1. This is slightly simpler than the form in [9] where all four possible terms must be kept.
• Both functions are real for u ∈ (0, 1) for both real and imaginary ν.
These are all nice properties for our solutions to have but it is possible that some other choice of basis may work. In [9] they fix β 0 using the u → −1 limit. Here we will instead find a β 0 that satisfies the single-cut rule and then show that, for real values of u our solution agrees with that in [9]. We will start by calculating the discontinuities of the cut diagrams, introducing the notation that (4.30) To equate this to the left hand side it is easiest to express each term on the left hand side in terms of P ± . We start by considering the case that |u| ≤ |v| and note that c mn are real so the sum is guaranteed to vanish on the left hand side 20 therefore note that the reality of c mn ensures that the sum cancels and this expression is symmetric in u, v, ensuring that we will find an identical result for |u| > |v|. By comparison with (4.30) it is straightforward to see that the single-cut rule holds if . (4.32)

Comparison with the Correlator
We have used the single-cut rule as a condition to fix the wavefunction coefficient. In order to use this to check that this rule holds for the result given in [9] we need to compare them. We start by relating the functions Q ± (u) to the F ± basis (4.20), where α ± are defined as in [9] by The ambiguity in the definition of these functions discussed under (4.20) can be ignored here because we are only interested in evaluating the correlator on the positive real axis and so the 20 The reality of cmn's is rooted in the fact that the m n . . . part of the correlator is equivalent to a sume over an infinite tower of quartic contact diagrams that emerges after integrating out the massive particle σ. Unitarity (in the form of the COT) and scale invariance (for conformally coupled fields) then demands the reality of each contact diagram, hence the reality of the coefficients cmn.
power of ν in α ± will appropriately cancel with that in F ± . The correlator is related to the wavefunction coefficient by As mentioned previouslyF ± are real for u > 0 and so taking the real part of ψ 4 reduces to taking the real part of the coefficients. We can then express allF ± in terms of Q ± and use (4.33) to relate these to F ± , details in Appendix A, so the s channel of the correlator is where we have introduced the function and defined, for the sake of comparison, β = − 1 sin(πν) . This is exactly the result found in [9] (except they use µ = iν). Therefore, by using the single-cut rule as a condition on the wavefunction of the universe coefficients, we find a result that is consistent with the results in the literature for a four point interaction involving the exchange of a massive scalar.

Four graviton exchange
We consider the case of graviton exchange with a vertex of the form [11]: where we have set the coupling constant to unity to simplify our notation. Setting M Pl = 1 for this subsection, the mode function and bulk-to-boundary propagator of the graviton are With this, we can easily obtain the wavefunction coefficient ψ λ 1 λ 2 λ 3 3 as ψ λ 1 λ 2 λ 3 bulk-to-bulk propagator Then we find s e ips(η+η ) .

(4.45)
Evaluating the integral gives us the following: (4.46) Here E L = k 1 + k 2 + p s , E R = k 3 + k 4 + p s , and k T = k 1 + k 2 + k 3 + k 4 . Adding the Hermitian analytic image reads: Figure 4: Diagram showing the geometry for the 5 point interaction which is being considered in (4.51). There are other diagrams that involve permutations of the labeling of the momenta that must be considered so that the wavefunction coefficient is symmetric in the momenta.
The power spectrum is given by: Therefore, we have: With this, it is straightforward to verify the single-cut rule for this interaction.

A five-point function
The simplest case in which we can't remove all derivatives from all the Green's functions by integration by parts is the five point function with a single three particle interactionφ 3 so we will consider this interaction here. Unfortunately, the five point function is very complicated and so writing it down here is rather unhelpful. Instead, we include a Mathematica file that demonstrates that the single-cut rules do hold in this case. We do this also for a second five point interaction that involves a combination ofφ 3 andφ∂ i φ∂ i φ vertices. Despite the complications with the explicit form of the five point function it was noted by Pajer and Hillman [63] that it is possible to express de Sitter wavefunction coefficients in terms of derivatives of flat space ones.
In the case at hand with onlyφ 3 interactions, this gives (4.50) Whereψ f 5 is a flat space wavefunction coefficient The geometry for this diagram is shown in Fig. 4. Here K f k (t) is the flat space bulk-to-boundary propagator, To calculate the discontinuities of the cut diagrams, first notice that the derivatives over k T in ψ 3 can equivalently be written as derivatives over k 12 so that the ψ 3 discontinuity is Disc p 1 iψ 3 (k 1 , k 2 , p 1 ) = Hk 2 1 k 2 2 p 2 1 ∂ 2 ∂k 12 ∂k 12 Disc p 1 iψ f 3 = iHk 2 1 k 2 2 p 2 1 ∂ 2 ∂k 12 ∂k 12 (4.61) Likewise, the k 12 derivatives in ψ 4 can be expressed as k 2 derivatives so that .

(4.62)
It is then possible to combine all of this together to give the single-cut rule for this 3 site chain, (4.63)

Conclusion
In this work, we have derived single-cut rules for the coefficients of the wavefunction of the universe that vastly extend the validity of the Cosmological Optical Theorem [1]. Our derivation leverages some simple analytic properties of the bulk-to-bulk and bulk-to-boundary propagators. Just like cutting rules in flat space, our results should be regarded as a consequence of unitarity.
In particular, our main achievements are summarized as follows: • We generalised the Cosmological Optical Theorem to an arbitrary number of spinning bosonic fields with a linear dispersion relation and arbitrary mass, arbitrary speed of sound and general local interactions at tree level (cutting rules for loops are discussed in [3]).
In particular, we explicitly checked that our relations are obeyed by the four-point scalar correlators from conformally coupled and general massive scalar exchange derived in [13,15]. We also discussed a four-graviton correlator from graviton exchange to demonstrate our treatment of spinning fields.
• We proved that the Cosmological Optical Theorem applies to all FLRW spacetimes where a Bunch-Davies initial state can be consistently chosen. This includes most spacetimes relevant for cosmology, such as de Sitter, slow and fast roll inflation and accelerating powerlaw cosmologies. We also checked that it applies to axion-monodromy inflation, where oscillations in the inflaton potential lead to a resonant particle creation and characteristic non-Gaussianities.
There are several directions for future investigation: • While valid to all orders in perturbation theory at tree level, our results do not give a non-perturbative statement of the Cosmological Optical Theorem. In analogy with flat spacetime, such a non-perturbative formulation would be highly desirable and could provide an important piece of the puzzle to derive positivity constraints on cosmological observables (see [64][65][66] for developments in that direction) and perhaps numerically bootstrap nonperturbative correlators, along the lines of [67,68] (see also [69] for recent non-perturbative results).
• It has recently been shown how to bootstrap cosmological correlators for massless scalars and tensors using the Cosmological Optical Theorem and a set of Bootstrap Rules [36,40,41]. Given our results here, it would be interesting to see if one can extend this derivation to the case of exchanged massive and possibly spinning fields, with potential applications to the cosmological collider phenomenology [13].
• It would be interesting to investigate what is the holographic interpretation of our single-cut rules in term of a hypothetical boundary field theory. Around de Sitter space one would expect the boundary theory to be a non-unitary CFT [70,71], but it is not clear what additional property needs to be satisfied to ensure that the bulk time evolution is unitary.
The fundamental and general nature of our results in this work strongly suggests that there are still basic and very general facts about quantum field theory on cosmological spacetimes that are awaiting to be discovered. Because of the ever growing body of cosmological dataset, advancements on the theory side are likely to have important repercussion on the phenomenology and ultimately make a long standing contribution to our understanding of the very early universe.
also interested in the proportionality factor that is neglected there. Therefore, in this appendix we calculate the three point function from the explicit bulk time integral, As in the main text we set H = 1 and expand this in terms of the conformally coupled modefunctions defined in (2.10) which gives where we have taken out the factor of k 3 2 3 from σ + k (η) so that it is a function of k 3 η only. We then define x = k 3 η so that Now we introduce the variable U = k 12 k 3 and consider the action of a slightly suggestive differential operator, (A.4) It is possible to express this in terms of x derivatives of the exponential, (A.5) which we integrate by parts We can exploit the differential equation satisfied by σ + (x), to give We are interested in the η 0 → 0 limit and so we define The homogeneous equation is identical to the one considered in [13] from the symmetries of the theory and is also the associated Legendre equation whilst the inhomogeneous part is independent of U and so the general solution is We want to avoid terms with spurious singularities and so B = 0. In order to fix A we explore the U → −1 limit Comparing this to the same limit in (A.11), . Therefore, For contact with the four point function we then define u = 1 U and take k 3 = p s so that 21 For later convinience we define this as This additional, inhomogeneous, term can be ignored if where α and β are some ν dependent constants. For imaginary ν this is always true because the term in brackets is bounded. For real ν this is true provided where this final condition also includes the case for imaginary ν.

A.2 General form of the four-point function
In this appendix we present the details of the calculation of the four-point function for a field of arbitrary mass. Just as for the three point function we will consider the bulk integral representation to ensure that we do not ignore any boundary terms that might be present, We can expand K ϕ to give We then use the differential equation solved by G σ , We can integrate this by parts to move all derivatives from G σ ps onto e ik 12 where they can be performed so that We can recognise this first term as ψ 3 and express each of the remaining terms as derivatives of ψ 4 so that This is once again the associated Legendre equation however the particular integral arising from the first term on the right hand side is not obvious. We will therefore employ the same tools that were developed in [9] by defining F = − η 4 0 s 4λ 2 ψ 4 , u = ps k 12 and v = ps k 34 so that For ∆ u defined as in [9], As we identified in Section 4.2, it is preferable to take the two linearly independent solutions to the homogeneous differential equation to bẽ We fix the particular integral arising from uv u+v just as in [9] by first considering a series expansion for |u| < |v|,F the series coefficients are c mn = (1) n (n + 1)...(n + 2m) Because these solutions have an identical Wronskian to those used in [9] the particular integral takes an identical form, (A.31) We also have a second particular integral that arises from the boundary term, We expect this to be symmetric in u, v so to find this particular integral we first look at the differential equation solved by ψ ϕϕσ So the homogeneous part of ψ ϕϕσ is also a homogeneous solution to the differential equation satisfied by F P therefore, to ensure symmetry in u, v we must add the homogeneous solution in u to this particular integral, This also ensures that it properly satisfies the differential equation in v. The only thing that remains is to find the complementary function which is defined so that the total solution is symmetric in u, v and free from unphysical singularities. The first of these conditions is implemented identically to [9] so that However, the removal of unphysical singularities as u, v → 1 occurs slightly differently. ψ 3 is already free from such singularities and the convergence of the sum is ensured by the piece-wise definition. The only term that can diverge is thereforeg. Physically u, v < 1 and so we must take the u → 1 limit with u > v, and we need to ensure that is finite. Therefore, we need β − = 0 and β 0 = −1 sog(u, v) becomes We will fix the remaining free coefficient, β + using the single-cut rule and so need to consider the Hermitian analytic image of this function. For this we need to use that Using this we find that To explore the Hermitian analytic image of F − we need that, [72], where the final equality holds for Im(u) > 0, thereforẽ It is convinient at this point to introduce ψ 3 includes only P + and so to prove that the single-cut rule holds it is convenient to re-express F − (u) in terms of P ± and we find that, for Im(u) > 0, Using this we can express theF − terms as The sum term is trivially equal to minus its Hermitian analytic image so, for |u| ≤ |v| the discontinuity of the four point wavefunction coefficient is where we have used that for p s > 0, Noting that this is symmetric when we exchange u and v we can conclude that this is also valid for |v| ≤ |u| which is required for the single-cut rule to hold. We now need to compare this to the discontinuity of the cut diagram which depends on .
From this we find the right hand side . (A.54) So this satisfies the single-cut rule if .

A.3 Comparison to the correlator
The results presented in [9] were expressed in terms of This is defined explicitly for real u so that Q ± can also be expressed in terms of hyper geometric functions, Therefore, so, to compare our results to those in [9] we should re-express everything in terms of Q ± . As we noted in Section 4.2,F ± are real for physical u, v and so we first take the real part and then use (A.47) to give π 2 cos(πν) Reg(u, v) = π 2 cos(πν)F Similarly, we need to find the real part of the three point function, (A.62) Where we have used that the Wronskian for the plus and minus solutions is −i. We can then calculate the trispectrum from this using that So, expressing the s channel in terms of F gives Therefore, for |u| ≤ |v|, we find We can also see that Therefore, this term cancels with the inhomogeneous contribution to the three point function and so We simplify this using that Re σ − ps (η 0 ) 2 − 2 Re σ − ps (η 0 ) 2 P σ ps = − Re(σ − ps (η 0 )) 2 + Im(σ − ps (η 0 ) 2 Re(σ − ps (η 0 )) 2 + Im(σ − ps (η 0 )) 2 = −1, where we have introduced β = − 1 sin(πν) . We now rewrite this in terms of F ± using (A.60) so that Finally, we simplify this using Euler's reflection formula [51], 4π sin(πν) cos (πν) = − 4π β cos(πν) .
(A.71) Therefore, (A.72) The expression on the final line is exactlyĝ(u, v) from [9]. We can find the expression for |v| ≤ |u| by exchanging u ↔ v and so we have This is exactly the same expression for B s 4 as (B.29) from [9] and so we find agreement with the results in [9].

A.4 The total discontinuity
We have shown that it is possible to cut an internal massive line in a way that is consistent with the massless case however, our more general results also rely on the ability to leave massive internal lines uncut. To see that it is valid to do this whilst also cutting an internal line would require the calculation of a new, complicated, diagram and so it is convenient to instead consider the total discontinuity of this diagram which, by the Hermitian analyticity of the propagators, is expected to vanish. As before we first restrict to the case where |u| ≤ |v|, We consider each of the terms in order Where we have used (3.42) as well as the fact that, for Im(p s ) < 0 to cancel the second line. Note that this same relationship also gives the total discontinuity in ψ ϕϕσ , which is what has previously been referred to as the contact COT [1]. The sum trivially cancels as c mn are real and the final term is g(u, v) −g(u * , v * ) = β + − β * + F + (u)F + (v) = 2iπ − π cos(πν) which is symmetric in u, v and so is valid for all u, v. This is exactly the result that was predicted by our consideration of the general case and so we have confirmed that uncut massive lines behave in a way consistent with our cutting rules.

B Resonant non-Gaussianity
In this appendix, we present some technical details of the proof that the propagators of perturbations in axion monodromy inflation are Hermitial analytic even in the presence of background oscillations that lead to the resonant production of perturbations. In order to make direct contact with the bulk-to-boundary propagator we will calculate φ + (as opposed to R k (−kη) ∝ φ − which was considered in [56]). To this end, we consider the ansatz φ + k (η) = C (−kη) ν H (2) ν (−kη) + b * c k (−kη)(−kη) where C is a k-independent constant and the factor of b * has been taken out of c k to make the relative size of the terms clear. The second term is of order 3 2 because ν = 3 2 + 2 0 + δ 0 and 0 , δ 0 = O(b * ) so to linear order it is enough to keep just this term. This is important because we don't have a convenient expression for the Hermitian analytic image of H (1) ν (x) for generic ν and Im(x) < 0. The differential equation that φ + satisfies is where A and B are x independent constants. Since c k must vanish in the infinite past we choose B = 0. Generically, for k with a negative imaginary part and real v we have that the Hermitian analytic image of the Hankel function is Since we are interested in the case where the mode function approaches e ikη in the far past, we make the following ansatz: σ(k, η) = ±σ 0 (η) + 1 k σ 1 (η) + 1 k 2 σ 2 (η) + This means that σ 1 is constant unless c s vanishes somewhere in the bulk. This constant can be absorbed into the normalization of the mode function, so we will ignore its contribution.
At O(k −2 ) we have: −2c s (η)σ 2 + m 2 a 2 − a a = 0, (C.7) Since everything within the integral is real, we expect σ 2 to be real as well. By induction, we see that σ r must be real for even r and pure imaginary for odd r. Therefore, as long as this series expansion converges, we conclude that σ(k, η) is Hermitian analytic.
As an example of how this WKB expansion gives us the mode function, let us consider the case of a massless scalar in de Sitter space with c s = 1. For de Sitter, a = 1 Hη , therefore a a = 2 η 2 . Since m 2 = 0, we have: Taking this result, we can obtain the other σ r order by order: (C.14) With the help of combinatorics and using induction, we can see that: (C.15) Therefore, we have: (C. 16) Setting C = −ik/ √ 2k 3 , and remembering that φ = f /a, we have which is the usual de Sitter mode function of a massless scalar.