Cosmological Correlators Through the Looking Glass: Reality, Parity, and Factorisation

We consider the evolution of quantum fields during inflation, and show that the total-energy singularities appearing in the perturbative expansion of the late-time Wavefunction of the Universe are purely real when the external states are massless scalars and massless gravitons. Our proof relies on the tree-level approximation, Bunch-Davies initial conditions, and exact scale invariance (IR-convergence), but without any assumptions on invariance under de Sitter boosts. We consider all $n$-point functions and allow for the exchange of additional states of any mass and integer spin. Our proof makes use of a decomposition of the inflationary bulk-bulk propagator of massive spinning fields which preserves UV-convergence and ensures that the time-ordered contributions are purely real after we rotate to Euclidean time. We use this reality property to show that the maximally-connected parts of wavefunction coefficients, from which total-energy singularities originate, are purely real. In a theory where all states are in the complementary series, this reality extends to the full wavefunction coefficient. We then use our reality theorem to show that parity-odd correlators (correlators that are mirror asymmetric) are factorised and do not diverge when the total-energy is conserved. We pay special attention to the parity-odd four-point function (trispectrum) of inflationary curvature perturbations and use our reality/factorisation theorems to show that this observable is factorised into a product of cubic diagrams thereby enabling us to derive exact shapes. We present examples of couplings between the inflaton and massive spin-1 and spin-2 fields, with the parity-violation in the trispectrum driven by Chern-Simons corrections to the spinning field two-point function, or from parity-violating cubic interactions which we build within the Effective Field Theory of Inflation.


Introduction
The fundamental observables of inflationary cosmology are late-time cosmological correlators, namely expectation values of quantum fields evaluated at the end of inflation. In the simplest models of inflation we are interested in correlations between the two massless fields that survive the rapid expansion of the background spacetime: the Goldstone boson of broken time translations, and the graviton. These correlators provide the initial conditions for Hot Big Bang cosmology, and observations of the Cosmic Microwave Background and Large Scale Structure then allow us to in principle distinguish between different models of the very early universe by measuring these spatial correlations. Given the very high energies that could characterise inflation, additional massive particles can be produced from the vacuum and decay into the light states that make it to the end of inflation. Such heavy states leave distinctive imprints on late-time correlators that encode their masses (through time evolution) and spins (through kinematics), thereby leading to the tantalising prospect of using the early universe as a very energetic cosmological collider [1][2][3][4][5].
Traditionally, such correlators are computed in de Sitter space perturbation theory using the inin/Schwinger-Keldysh formalism (see e.g. [6][7][8] for reviews), or the Wavefunction of the Universe. Such computations quickly become very complicated due to the background time dependence which, in contrast to flat-space, means that the time integrals one must compute are very non-trivial. Further complications arise from the fact that the mode functions of massive fields in de Sitter space are usually characterised by Hankel functions (or similar functions) which are harder to integrate compared to the plane wave solutions which are familiar in flat-space. Such complications provide a stumbling block in ones quest to understand the late-time effects of massive fields during inflation, and more generally to understand the fundamental properties of cosmological correlators.
In this work we focus on the Wavefunction of the Universe and the associated wavefunction coefficients which contain the dynamical information about the evolution of quantum fields in de Sitter space. It is now well-known that for fields that satisfy Bunch-Davies initial conditions, i.e. for fields that have Minkowski-like behaviour in the far past, wavefunction coefficients have a restricted set of singularities. Indeed, wavefunction coefficients are only singular when the total-energy of the external states vanishes, or when the total-energy entering a sub-graph vanishes. 4 For physical configurations i.e. for real momenta, such singularities cannot be reached, however much of our understanding of wavefunction coefficients and cosmological correlators stems from analytically continuing away from real momenta in which case energies can become negative and singularities can be probed. For example, the leading total-energy singularity of n-point functions allows us to probe the corresponding flat-space (boost-breaking) scattering amplitude [22,64,65,68] which provides a non-trivial link between the cosmological bootstrap and the S-matrix bootstrap, while at four-points the additional singularities are partial-energy ones on which wavefunction coefficients factorise into a product of a lower-point wavefunction coefficient and a scattering amplitude, see e.g. [18]. When wavefunction coefficients are rational, which is often the case for massless scalars and gravitons [85], the residues of these partial-energy poles can be fixed by unitarity and energy shifts [26,33].
In general, wavefunction coefficients in momentum space are complex functions of the external kinematics and cosmological correlators correspond to taking the real or imaginary parts, at least at tree-level. For parity-even interactions it is the real part that contributes to correlators, while for parity-odd interactions it is the imaginary part [37]. While it is known that the tree-level wavefunction coefficients in simple massless theories are purely real [86] (thereby making parity-odd correlators a neat probe of exotic inflationary physics), in this paper we greatly extend this concept of reality. More precisely, we show that tree-level wavefunction coefficients with massless scalar external states have purely real total-energy poles (both leading and sub-leading) as long as the bulk interactions are IR-convergent. Our proof is valid for all n-point functions at tree-level and we allow for Feynman diagrams corresponding to the exchange of massive spinning fields of any mass and integer spin. More concretely, we show that contributions to wavefunction coefficients from the maximally-connected parts of Feynman diagrams are purely real, where "maximally-connected" corresponds to the contributions to the integrands with the maximal number of θ-functions (and this is where total-energy singularities come from). Our proof makes use of Wick rotations of the time integrals that compute these wavefunction coefficients and a simple property of the bulk-bulk propagator of massive fields in de Sitter space: the time-ordered part is purely real after the time variables are both Wick rotated by 90 • in the complex plane. This property follows as a simple consequence of the differential equation that the bulk-bulk propagator must satisfy given that it is a Green's function. We provide more explicit and detailed proofs within two different set-ups for describing massive spinning fields during inflation: cosmological condensed matter physics (CCM) and cosmological collider physics (CC).
The former was introduced in [87] and requires a sizeable coupling between the new massive degrees of freedom and the time-dependent inflaton. In this set-up fields are classified with respect to how they transform under the unbroken group of symmetries which for cosmology is spatial rotations. The theory should in addition have all of the symmetries of the Effective Field Theory of Inflation (EFToI) [88] and this can be guaranteed thanks to particular couplings with the inflaton which can be built straightforwardly using the building blocks of the EFToI [87]. From the point of view of spontaneous symmetry breaking and the coset construction, such new degrees of freedom are classified as matter fields which can couple to the Goldstone boson of broken time translations. The masses of these fields are not restricted by the Higuchi bound [89], which illustrates that they cannot exist in an exactly de Sitter invariant theory. This set-up is somewhat similar to condensed matter systems where linearly realised Lorentz boosts are not used to construct the effective theory (and hence the name). The latter is perhaps more familiar to the reader and corresponds to describing massive degrees of freedom as representations of the de Sitter group. The masses of these fields must adhere to a lower bound in order for the theory to remain unitary, which is the Higuchi bound [89]. The Lagrangians for such fields are known [90][91][92], but quickly become very complicated due to the need to include auxiliary fields (which ultimately enforce the transverse and traceless conditions on the fields as required by the degrees of freedom counting). One can instead work directly with the equations of motion which take a simpler form [93]. Although the free theories for these new degrees of freedom are de Sitter invariant, de Sitter boosts can be broken when we come to couple these fields to the inflaton. This can again be done within the language of the EFToI as done in [5]. We will review both set-ups in Section 2.
In each case we consider light and heavy fields i.e. those in the complementary and principle series, respectively. We concentrate on the reality properties of the Wick-rotated bulk-bulk propagators which, in cosmology, are composed of the usual time-ordered Feynman propagator and a factorised term necessitated by the boundary conditions. For light fields, we show that the full bulk-bulk propagator is purely real after Wick rotation. For heavy fields, we show that although the full bulk-bulk propagator is complex in general, we are able to add and subtract factorised contributions to the full bulk-bulk propagator in such a way that we cancel the imaginary parts of the Wick-rotated Feynman propagator. The new timeordered part, which enjoys reality after rotation, and which is now a sum of the Feynman propagator and a factorised contribution, is referred to as the connected propagator. This decomposition is diagrammatically represented in (1.1). In addition, this connected propagator enjoys the crucial property that it vanishes in the far past, which ensures that Feynman diagrams constructed using this propagator maintain UV convergence of the associated time integrals. The detailed treatment of the propagator realities differs from one set-up to another. In the CCM set-up, we allow for parity violation in the free theory of the massive spinning field coming from a chemical potential term with a single spatial derivative in the action [38,[94][95][96]. This splits the helicities and changes the mode functions from Hankel functions to Whittaker functions. The propagator realities then require cancellations once we sum over the helicities.
In the limit of a vanishing chemical potential with a parity-conserving propagator, reality holds for each helicity mode separately with the proof for light fields already appearing in [37]. In the CC set-up we maintain parity in the free theory of the massive spinning field (except for spin-1 which is essentially identical to the CCM case) and again show that for light fields the Wick-rotated propagator is purely real, for each helicity mode, while for heavy fields it is again only the connected part that is purely real. In order to arrive at this conclusion we use the fact that the transverse and traceless conditions for these fields relate modes with the same helicity via differential operators that have simple properties under Wick rotation. These proofs make up Section 3.
In Section 4 we then use these general properties of Wick-rotated bulk-bulk propagators to prove that total-energy poles of wavefunction coefficients with external massless scalars are purely real under the assumptions of IR-convergence, scale invariance and the tree-level approximation. Our proof does not rely on de Sitter boosts and is therefore directly applicable to inflationary correlators. We also point out that our proof applies to external gravitons, and to wavefunction coefficients with an even number of conformally coupled scalars. Our approach is to extract the connected part of the full bulk-bulk propagator and show that diagrams that only involve connected propagators (maximally-connected diagrams) are purely real, and indeed total-energy singularities come from such diagrams. This is a general result, but if we restrict ourselves to the exchange of light fields only, then the full wavefunction coefficient is real. In Appendix C, we offer a complementary proof of the reality of total-energy poles using the Hermitian analyticity properties of the external bulk-boundary propagators and the internal bulk-bulk ones. This property combined with exact scale invariance allows us to draw the same conclusions we arrived at using Wick rotations. Hermitian analyticity of bulk-bulk propagators in the context of the CCM scenario has been established in [24], and in Appendix C we extend the analysis to the CC scenario.
As a proof of the usefulness of this observation, we consider the parity-odd scalar trispectrum in Section 5, which has recently gained some attention [37,[97][98][99][100][101][102][103][104]. This correlator is fixed by the imaginary parts of wavefunction coefficients and given that total-energy poles are real, they do not contribute to the parityodd trispectrum (unless there is some IR divergence as in [101], or if they come from loop diagrams as in [42]). This implies that this observable is in fact factorised at tree-level, and can only have partial-energy poles. In computing cosmological correlators the primary difficulties arise when one computes the nested time integrals coming from the time-ordered propagators. Since here these contributions are purely real, computing the parity-odd trispectrum due to the exchange of massive spinning fields reduces to computing lower-order, factorised time integrals which have known closed-form solutions. We present a number of exact parity-odd trispectra for both the CCM and CC descriptions of massive spinning fields focusing primarily on spin-1 and spin-2. We consider different sources of parity-violation: parity-violating bulk interactions and parity violation in the free theory describing the massive field. The latter case is usually studied in the context of cosmological chemical potentials, where it is known that the chemical potential can assist particle production [96,105] and boost the cosmological collider signal [38,41,94,95,[106][107][108][109][110]. In all cases we also allow for a general speed of sound for both the inflaton and the exchanged field which can also boost the signal if the inflaton moves more slowly [34]. We consider both light and heavy fields, and show that the final correlators for these two cases can be converted into each other by analytic continuation. We pay special attention to the spin-1 case with a parity-violating propagator (this corresponds to taking the Proca action and adding a Chern-Simons term that mixes the massive spin-1 field with the inflaton), and compare our exact result for the corresponding parity-odd trispectrum with that recently derived in [104] using a non-local EFT approach. We find that the result from the non-local EFT provides a good approximation to our exact result for small chemical potential, but deviations start to appear as the magnitude of the chemical potential is increased (as expected from numerics [104]). For comparison, the exact result we present here holds for all values of the chemical potential within the perturbative regime.
Our results do not just imply that parity-odd trisepctra are factorised; rather any n-point correlator of massless scalars and massless gravitons that requires one to take the imaginary part of wavefunction coefficients will be factorised and therefore lend itself to a simpler computation. This includes all n-point parity-odd correlators.
We conclude and discuss further directions in Section 6. We also include a number of appendices. In Appendix A we provide full details on how to realise the reality property of the time-ordered part of the Wick-rotated propagator. In Appendix B we provide a proof of our reality theorems using the in-in formalism which might be more familiar to the reader. Appendix C includes a complementary proof of our results using Hermitian analyticity that we have already mentioned above, and in Appendix D we extend our results to other Friedmann-Lemaître-Robertson-Walker (FLRW) spacetimes.

Summary of main results
Before moving to the main body of the paper, let us first summarise our main results for the convenience of the reader.
• Throughout this paper we consider the following decomposition of the wavefunction bulk-bulk propagator G: where we refer to C as the connected part and F as the factorised part. The connected part is defined as containing the time-ordered parts of the full propagator, along with the crucial property that it is purely real after both time variables are Wick rotated by 90 • in the complex plane. As we explain in detail in Section 3, C is not equal to the Feynman propagator; they differ by a factorised contribution. For fields in the complementary series C coincides with G i.e. F = 0, however for fields in the principle series F ̸ = 0. In all cases both C and F vanish in the far past, while only the sum satisfies the future Dirichlet boundary condition with all η 0 dependence (where η 0 is the late-time cut-off) contained in F .
• Using this decomposition and the reality properties of C, we derive reality properties of wavefunction coefficients of massless scalars and gravitons with our main results depicted in Figure 1.
• If we specialise to four-point correlators of inflationary curvature perturbations then the relationship between such a correlator and wavefunction coefficients is depicted in Figure 2.
• We used our factorisation theorem to compute some exact parity-odd trispectra of curvature perturbations. We consider a number of examples in Section 5 but in all cases the final trispectrum takes the following schematic form: where (L, R) correspond to kinematic structures with partial-energy (E L,R ) singularities. There are no total-energy singularities. In all cases the time evolution is characterised by a product of hypergeometric functions.
(4.1) (4.2) (4.3) Figure 1: Scale invariance plus a Bunch-Davies vacuum implies that the maximally-connected parts of wavefunction coefficients are purely real where the maximally-connected part is given by replacing all bulkbulk propagators G with the connected propagators C. Such a maximally-connected part is depicted as a hatched grey blob in this figure. We refer to this result as the k T -reality since all total-energy singularities come from the maximally-connected wavefunction coefficients. This is a general result that applies for the exchange of fields of any mass and integer spin. If the exchanged fields are light i.e. they are in the complementary series, then the full wavefunction coefficient is real (not just the maximally-connected part) since in this case C coincides with G. If we further consider parity-odd correlation functions, then the wavefunction reality implies that such n-point correlators are factorised since only the imaginary part of maximally-connected wavefunction coefficients contribute to these observables and the imaginary parts vanish.
Notations and conventions Throughout this paper, we adopt natural units c = ℏ = 1 and work with the (− + ++) metric sign convention. Fourier transforms are defined by We will always adhere to exact translational and rotational invariance. In fact, we shall work in the Poincaré patch of de Sitter space with the metric as an approximation for an inflationary spacetime. We will mainly work with the wavefunction formalism and define the wavefunction coefficients by the expansion Figure 2: The first line is the usual relationship between (s-channel) four-point correlators and wavefunction coefficients. In the second line we have decomposed the full bulk-bulk propagator into the connected and factorised parts. The k T -reality tells us that the maximally-connected part is purely real which in turn implies that the parity-odd part of the four-point function is factorised (as above).
where φ denotes a general field with indices suppressed. The bulk-boundary and bulk-bulk propagators in the wavefunction formalism are given by with P (k) = φ(η 0 , k)φ * (η 0 , k) denoting the power spectrum at η = η 0 . When computing the wavefunction coefficients, we adopt the amplitude Feynman rules by including a factor of i for every vertex (our bulkbulk propagator therefore differs from that in [22] by a factor of i). Wick rotation turns out to be extremely crucial in this paper. Hence for clarity, we adopt χ > 0 to denote the Wick-rotated conformal time, defined by Under this transformation, the propagators are dressed with tildes to indicate Wick rotation, (1.10) We will pay much attention to the total energy of a diagram with external momenta {k 1 , · · · , k n } given by where k a = |k a |, a = 1, · · · , n are the external energy variables. Note that correlators with a prime denotes the removal of an overall momentum-conserving δ-function, e.g.
In the case of 4-point correlation functions, we evoke the Mandelstam-like variables which satisfy the non-linear relation by momentum conservation. We define a dimensionless curvature trispectrum T by For spinning fields, the CCM and CC set-ups historically chose different conventions for polarisation tensors. We choose to respect these distinct conventions and use different fonts to avoid confusion: (1.16) The different properties of these tensors will be explained in Section 2. Throughout this paper, we make use of a number of mathematical formulae for Hankel and Whittaker functions, which can be found in Chapter 10 and Chapter 13 of NIST [111].

Massive spinning fields during inflation
As we explained in the introduction, in this work we are interested in n-point functions of massless scalar fields (and gravitons) that can be generated at tree-level due to interactions with massive spinning fields. In this section we introduce these massive spinning fields and derive their mode functions. We consider the two different cases of interest, cosmological condensed matter physics and cosmological collider physics, separately. This section does not contain any new results, so readers familiar with these descriptions of massive spinning fields (including parity-violating corrections from the chemical potential) can skip to Section 3.

Cosmological condensed matter physics
We begin with the description of massive spinning fields during inflation advocated in [87]. The idea is to classify states with respect to the unbroken rotational symmetries, rather than as representations of the full de Sitter group. Fields therefore only have spatial indices, and we denote a field of spin S by Σ i1···i S . For this field to carry 2S + 1 degrees of freedom, it should be traceless but not transverse. To ensure that the symmetries of the EFToI are respected, this field is promoted to Σ µ1···µ S where the new temporal components depend on the Goldstone of broken time translation π. For example, for S = 1 we have [87] Σ 0 (π, We then write down a quadratic action for Σ µ1···µ S which will introduce quadratic terms for the spatial components, and interactions between the spatial components and the inflaton with the same coefficients. This illustrates the fact that such a theory cannot exist in the absence of the inflaton. The general action with at most two derivatives is where n µ n µ = −1 is a timelike unit vector that defines a preferred frame in which the spatial rotations remain intact. Notice that in addition to the first four terms which appear in [87], we have also included a fifth term which has a single derivative and is parity-odd. We in principle have five free parameters but we can fix one using our freedom to normalise Σ µ1···µ S which leaves us with four: c, δc, m 2 and κ which respectively correspond to two speed-of-sound parameters, the mass, and the chemical potential. The free theory for Σ i1···i S is then where we have defined σ i1···i S = a −S Σ i1···i S and have converted to conformal time. Here all scale factors are manifest and indices are raised and lowered with the Kronecker symbol δ ij . We see that the kinetic term for this field is the same as that of a canonical scalar in de Sitter. If we had instead tried to directly construct the most general action with at most two derivatives for σ i1···i S that respects rotational invariance and scale invariance, we would also have arrived at (2.4). 5 For Σ i1···i S to be traceless, the trace has to be taken with respect to the induced metric on constant inflaton slices [87]. This implies that in (2.4) we can take σ i1···is to be traceless with respect to δ ij up to terms that are quadratic in σ i1···i S and at least quadratic in π. When discussing the free theory for this massive spinning field we therefore take it to satisfy δ ij σ ijl3···l S = 0. We now convert the action to momentum space and decompose the field in terms of its helicities via where σ h (η, k) are the mode functions and the traceless polarisation tensors satisfying with the first condition following from the reality of the fields in position space, and the second a normalisation choice. For a given helicity the polarisation tensor is a function ofk and two polarisation directionsê ± which are orthogonal tok, and satisfy (2.6) and (2.7). More explicitly, we can construct the polarisation directions asê while for S = 2 we have Further details on these polarisation structures can be found in e.g. [5,24]. Using these properties of the polarisation tensors, (2.4) becomes a decoupled action for each mode function: In deriving this expression we have used iϵ ijk k jê ± k = ±kê ± i . The fact that the parity-violating term splits the helicities is now clear given the helicity factor h in the final term. In addition to this term, we see that each mode has the same mass but different speed of sound which depends on the two original speed of sound parameters, the spin and the helicity. The equation of motion for each helicity mode is then The solution to this equation with Bunch-Davies initial conditions i.e. the solution that has Minkowski-like behaviour in the far past and satisfies the de Sitter Wronskian condition is where W a,b (z) is the Whittaker W -function. This solution is valid for both light (Im ν = 0) and heavy (Re ν = 0) fields. In the limit that the chemical potential vanishes we would expect to recover the solution of [87] corresponding to the solution of a massive scalar field in de Sitter space. This can be verified using the relation where H (1) ν (z) is the Hankel function of the first kind. We will discuss the corresponding power spectra and propagators in Section 3. Notice that in this CCM scenario, the introduction of chemical potential κ is quite natural since it serves as a next-to-leading order correction to the dispersion relation of massive spinning fields in the gradient expansion, and is consistent with spatial rotations and scale invariance. More precisely, it appears as a linear term in momentum in the dispersion relation of a massive field, where k p ≡ k/a(t) and S are the physical momentum and the spin angular momentum of the field mode, respectively. Such a linear correction to the massive field's dispersion relation is not sign-definite, and alters the analytic structure of its equation of motion (2.13), leading to enhanced particle production when m ≳ κ [96]. In the case where m ≲ κ, however, modes with a negative linear term may experience a transient tachyonic phase where w 2 < 0 and grow exponentially, i.e. σ ∼ e −iωt ∼ e |ω|t . Such a tachyonic growth is eventually halted by the finite mass, leading to w 2 ≈ m 2 > 0 in the IR limit k p → 0. Nevertheless, the exponential growth during the tachyonic period may overproduce particles and threaten perturbativity or even the inflationary background. Therefore, in favour of theoretical control, we will require m − κ ≳ −H throughout this paper.

Cosmological collider physics
We now turn to the more familiar description of massive spinning fields in de Sitter/inflation where they are representations of the de Sitter group. In this section we primarily follow [5] and the refer the reader there for further details. Such fields can certainly exist in the absence of the inflaton, in contrast to CCM. A spin-S bosonic field in this CC set-up is described by a symmetric rank-S tensor that satisfies: where the mass parameter is m 2 On-shell we therefore have a transverse and traceless rank-S tensor that satisfies a wave equation. To solve this system we work in a 3+1 decomposition where the components are of the form Φ η···ηi1···in with 0 ⩽ n ⩽ S. We further convert to momentum space and decompose each of these components in terms of helicities via Each mode function is therefore labelled by three numbers corresponding to the spacetime spin (S), spatial spin (n), and helicity component of the spatial spin (h). We have made a distinction between the polarisation tensors used here (e) compared to above in the CCM case (e). This is because a different normalisation is employed in [5] compared to what we used above, and while discussing the CC scenario we will follow the conventions of [5] to (hopefully) avoid confusion for the reader. The polarisation tensors here are still functions ofk and two polarisation directionsê ± , however rather than satisfying (2.6) and (2.7), they are chosen to satisfy where the polarisation tensor with lowest index is defined as In both cases (CCM and CC) the numerical factors that multiply thek andê ± factors can be made purely real or purely imaginary. The magnitudes of these factors are fixed by (2.7) and (2.20) in the two cases, while in the CCM case the condition (2.6) fixes the factors to be imaginary when there are an odd number ofk's, and real when there is an even number, while (2.19) fixes the factors to be always real. This phase difference will ultimately be inconsequential since it can be absorbed into the mode functions which are always only fixed up to a phase.
As an illustration let us spell out the S = 1 case. If we decompose Φ µ into its time and space components Φ η and Φ i , then (□ − m 2 while the transverse constraint is The Φ η field carries only a h = 0 mode, while the Φ i components carry both h = 0 and h = ±1 modes. We then write (2.25) The polarisation vectors are chosen to be The equations of motion then decouple for each mode function and are given by subject to the transverse constraint The equation of motion for the h = ±1 modes does not contain the Hubble friction term since the Maxwell kinetic term is conformally invariant. The solutions to these equations with Bunch-Davies initial conditions are Hk where 34) and the normalisation constants have been fixed by demanding that the commutation relations are the usual ones [5]. We see here a feature that immediately distinguishes this set-up from the CCM one: some of the mode functions in this case are given by a sum of Hankel functions with degenerate order parameters. The dynamics in the two set-ups therefore different.
The story for spin-S is similar. Modes with helicity h can come from all components with n ⩾ |h|, and those with n = |h| satisfy [5] Φ h |h|,S The other mode functions with the same helicity but with n > |h| are then obtained iteratively from the transverse and traceless conditions which fix 6 where . (2.37) One can also solve the recursion relation above and write the mode functions with higher spatial spin as a linear differential operatorD acting on the lowest-spatial-spin one, This form of the general mode function will become useful later in Section 3.2. Note that for a given n, B m,n is only non-zero if n and m differ by an even number. This is because the terms proportional to B m,n come from subtracting traces from Φ η···ηi1···in . The factor of i in this expression comes from converting to momentum space, and no additional factors of i appear since here the coefficients of all polarisation tensors are taken to be real. The solution to (2.35) with Bunch-Davies initial conditions is which is again fixed by demanding that we have the usual commutation relations [5], and As we saw explicitly for S = 1, the mode functions with n > |h| will be given by a sum of Hankel functions. The different solutions fall into different classes depending on the mass of the field. The principle series corresponds to heavy masses with Re ν S = 0, i.e. m 2 ⩾ H 2 (S − 1/2) 2 . The complementary series corresponds to light masses where Im ν S = 0, however, the Higuchi bound sets a lower bound on the mass such that the theory remains unitary. 7 The complementary series is then defined by S(S − 1) < m 2 /H 2 < (S − 1/2) 2 . Finally, we have the discrete series for which m 2 = H 2 [S(S − 1) − T (T − 1)] for S, T = 0, 1, 2, · · · , with T ⩽ S. In these cases there is an additional gauge symmetry that reduces the number of propagating degrees of freedom and corresponds to partial masslessness [112,113]. We will discuss the corresponding power spectra and propagators in the following section.

Mass parameter comparison
As emphasised in [87], in the CCM case the masses of the fields have a wider range as they are not constrained by the Higuchi bound. For comparison, in Figure 3 we show the distribution of the dimensionless mass parameter on the complex plane for light/heavy fields in both the CCM and the CC scenarios. From (2.13) we see that for light fields in the CCM scenario ν cannot exceed the massless boundary of ν = 3/2, while for heavy fields we can keep increasing the mass to keep increasing the value of Im ν. Similarly, in the CC scenario we see from (2.41) that Im ν S continues to increase as we increase the mass while in the principle series. For the complementary series the mass parameter can take values 0 < ν S < 1/2, as dictated by the Higuchi bound. The discrete series, which also lines along the real line in the complex plane of ν S , corresponds to isolated points and with masses that we will also refer to as "light".

General properties of Wick-rotated propagators
We now move on to the propagators that are required to perturbatively compute wavefunction coefficients and cosmological correlators. For each field we have a bulk-boundary propagator for external lines in Feynman diagrams, and a bulk-bulk propagator for internal lines. In this section we will discuss some important properties of the bulk-bulk propagators of massive spinning fields: we will show that after Wick rotation, the time-ordered parts are always purely real regardless of the mass and spin, while for light fields i.e. those in the complementary series, the full bulk-bulk propagator is purely real. To derive these properties we consider the CCM and CC cases separately since the proofs are slightly different in the two cases. For details on the Feynman rules for wavefunction calculations we refer the reader to [22].

Cosmological condensed matter physics
We begin with the power spectrum of massive spinning fields introduced in Section 2.1. The late-time two-point function is where This expression is valid for both light and heavy fields (i.e. for both purely real and purely imaginary ν), and we have used [W a,b (z)] * = W a * ,b * (z * ) and the symmetry property W a,−b (z) = W a,b (z). The power spectrum of each helicity mode is different due to the c h,S andκ dependence. Parity violation is then encoded in the asymmetry between opposite helicities. It can be further shown that the IR expansion of the power spectrum contains enhanced oscillations in η 0 due to particle production assisted by the chemical potential, which leads to lifted cosmological collider signals (see, e.g. [38]). The bulk-boundary propagator of a given helicity mode, which we will denote as K (h) σ (η, k), should satisfy the equation of motion (2.13) subject to the boundary conditions: The first condition simply requires us to add an appropriate normalisation factor, while the second condition requires the bulk-boundary propagator to be fixed by σ * h (η, k) since this is the negative frequency solution that vanishes in the far past and projects onto the Bunch-Davies vacuum. We can therefore write the indexed bulk-boundary propagator as a linear combination of the helical bulk-boundary propagator The factor of e j1···j S (−k) in (3.5) ensures that the different helicity modes remain decoupled during propagation, and the factor of 1/S! follows from the normalisation (2.7). The bulk-bulk propagator for each helicity mode, which we denote as the helical bulk-bulk propagator G and is subjected to the boundary conditions: The second condition again ensures that we project onto the vacuum in the far past while the first ensures that this propagator takes us between two bulk points rather than from a bulk point to the boundary. The solution to this equation is then The manifestly time-ordered parts of this expression correspond to the usual Feynman propagator (timeordered two-point function), while the η 0 -dependent terms are required to satisfy the future boundary condition. We can again dress the helical propagators with polarisation tensors and write the indexed propagator which satisfies where we symmetrise according to e.g. S (ab) = (S ab +S ba )/2. We would now like to consider the properties of this indexed bulk-bulk propagator after we Wick rotate both time variables by 90 • in the complex plane.
We are therefore interested in the properties ofG First consider the LHS of (3.12) where the differential operator that acts on the bulk-bulk propagator is purely real after Wick rotation: the only term that is odd in η 1 comes with a factor of i which itself comes from converting to momentum space. It is scale invariance of the quadratic action that ensures that this differential operator is real after Wick rotation. Indeed, time derivatives are forced to appear as η ∂ ∂η , while spatial derivatives yield factors of iηk. We therefore have a real differential operator acting on the indexed propagator. Moving to the RHS, we first note that the polarisation sum over all helicity modes is purely real which follows from the reality of the fields in position space. 8 After Wick rotation we also do not have the factor of i on the RHS since in this "Euclidean" picture we are computing ψ ∼ e −S rather than ψ ∼ e iS . The RHS is therefore manifestly real after Wick rotation. We therefore conclude that the time-ordered parts of the Wick-rotated bulk-bulk propagator, which are the parts that are fixed by (3.12), are purely real. This does not imply that the full bulk-bulk propagator is purely real after Wick rotation as this discussion still allows for the possibility of adding complex contributions that satisfy the homogeneous equation. This also does not imply that the Feynman propagator i.e. the manifestly timeordered parts in (3.9) are purely real after Wick rotation; we may have to add factorised terms that satisfy the homogeneous equation to make this reality property manifest (this is precisely what happens as we will see below). 9 As long as we are considering fields that are real in position space with scale invariant free theories, this property of the bulk-bulk propagator will hold. Let's see this explicitly by working with (3.9).
Light fields First consider light fields (Im ν = 0) where it was shown in [37] that whenκ = 0 the manifestly time-ordered parts are not real, but once we add the factorised term as required by the future boundary condition, the full bulk-bulk propagator is real. We can check if this remains true whenκ ̸ = 0. We first note that for light fields the bulk-bulk propagator is independent of η 0 . Indeed, we have The bulk-bulk propagator for a given helicity mode can then be written as where we have introduced the Whittaker-M function M a,b (z) which is related to W a,b (z) by In order to arrive at (3.14) we have used the rotation that does not cross the branch cut of W a,b (z) which lies on the negative real axis. We have also used 16) to make sure that the two arguments in (3.14) have the same sign. Note that (3.13) follows from (3.15) since for small arguments W a,b (z) dominates over M a,b (z) (for light fields). The function M a,b (z) is another solution to the Whittaker differential equation, and also satisfies [M a,b (z)] * = M a * ,b * (z * ). We can now consider Wick rotating (3.14), and to do so we rotate both time variables clockwise to avoid the branch cuts on the negative real axis (recall that η 1 , η 2 ⩽ 0). We then havẽ with χ 1 , χ 2 ⩾ 0. The structure of the θ-functions can be understood with a simple flat-space toy example. Consider the bulk-bulk propagator [78] G where the θ-functions ensure convergence in the far past when η = −∞(1−iϵ). Now consider the clockwise rotation we used above. Since we now integrate from 0 to +∞, we need the exponential damping to come from the largest of the two variables. We therefore havẽ as in [15]. This example captures the important features for our bulk-bulk propagator given that the Whittaker functions are also exponentially damped (W )/growing (M ) for large argument.
We can now use the properties of the Whittaker functions and the Γ-functions to conclude that taking the complex conjugate of (3.17) can be compensated for by sendingκ → −κ which is equivalent to h → −h. We therefore conclude that the bulk-bulk propagator is helically real i.e.
This very nice relationship between the bulk-bulk propagator of different helicity modes is not quite enough for us to conclude that the full propagator in (3.11) is real; we also need a relationship between the product of polarisation tensors of different helicities. The reality of this full polarisation sum implies that since contributions with different |h| have a different number of k factors so cannot be related by complex conjugation, while a single contribution from a single helicity has both real and imaginary components. These two relationships therefore allow us to conclude that the full indexed bulk-bulk propagator (3.11), when the field is light, is real after Wick rotation: Heavy fields Let's now consider heavy fields where Re ν = 0. For convenience let us therefore write ν = iµ with µ > 0. The primary difference here compared to the light field case is that the η 0 dependence in the bulk-bulk propagator no longer cancels out. This is perhaps easiest to see in theκ = 0 limit where the mode functions are Hankel functions. As we send η 0 → 0 − , each Hankel function has two comparable oscillating contributions, in contrast to the light field case where they both have one decaying contribution and one growing one, and this ensures that the ratio is time-dependent. Forκ ̸ = 0 the story is the same. Furthermore, one can easily check that the Wick-rotated time-ordered parts of G (h) σ do not satisfy a relation of the form of (3.20), and therefore the Wick-rotated Feynman propagator once we sum over the helicities is not manifestly real. In contrast to the light field case, the factorised terms we add to satisfy the future boundary condition cannot cancel the imaginary parts of the rotated Feynman propagator since they depend on η 0 while the Feynman part does not. However, as we have already alluded to, it is possible to add and subtract factorised contributions in such a way that we can make the time-ordered parts of the bulk-bulk propagator manifestly real after Wick rotation.
To see how this can work let us decompose the full bulk-bulk propagator into two parts which we will refer to as the connected part (C) and the factorised part (F ): where for a given helicity mode these parts of the propagator are where we have added a new contribution to each, ∆G σ , in such a way that the full propagator is unchanged. Diagrammatically, we represent this propagator decomposition as (3.26) This correction term must satisfy the homogeneous equation of motion and so is factorised. We would now like to fix ∆G is manifestly real after Wick rotation. Our argument above suggests that such a ∆G (h) σ exists. There are in principle a number of possible structures that ∆G (h) σ can take, but we would like each term in this decomposition of the bulk-bulk propagator to vanish in the far past such that integrals involving G can be split into integrals involving C and F without losing convergence in the far past. Given that the Feynman propagator vanishes in the far past, thanks to the time ordering, ∆G (h) σ must also vanish. We therefore take it to depend on σ * h rather than σ h . We would also like it to be symmetric under the exchange of η 1 and η 2 such that both C and F are symmetric. We therefore take is a η 0 -independent constant, and so is distinct from the terms we need to add to satisfy the future boundary condition. We now want to fix A h by demanding the connected propagator is helically real after Wick rotation: This condition, along with (3.21), ensures that the connected part of the full bulk-bulk propagator is real after Wick rotation i.e.
To summarise, we define the helical connected bulk-bulk propagator to satisfy the following conditions: satisfies the same equation as the bulk-bulk propagator, i.e. it satisfies (3.7).

C
(h) σ vanishes exponentially fast in the far past under the iϵ-prescription.
In order to fix A h it is wise to first write the connected propagator in terms of M a,b (z) only since its analytic continuation satisfies a simple relation c.f. (3.16), compared to that of W a,b (z) c.f. (3.15). We can eliminate all copies of W a,b (z) using and use (3.16) to make sure that all arguments lie on the positive imaginary axis such that we can rotate each time variable by 90 • clockwise to make all arguments lie on the positive real axis. It is then a simple task to demand (3.28) and fix A h . The most general solution of A h is derived in detail in Appendix A. In short, we find This now gives us a connected bulk-bulk propagator that is manifestly real after we Wick rotate, which we will use in Section 4 to prove the reality of total-energy poles, and a factorised bulk-bulk propagator which we will use in Section 5 to compute exact parity-odd trispectra. 10

Cosmological collider physics
We now turn our attention back to the cosmological collider physics set-up and derive properties of the bulkbulk propagator. We do not find it necessary to discuss the power spectra and bulk-boundary propagators here since they are not required for our needs and details can be found in [5]. The full propagator with covariant indices would naturally be written as G µ1···µ S ν1···ν S (η 1 , η 2 , k) but as we did when we discussed the mode functions, we will consider a 3 + 1 decomposition and consider the properties of propagators with spatial indices: where as before 0 ⩽ n ⩽ S. We can further decompose into helicities and write For a given n, the mode functions with the same |h| are equivalent since here we do not consider parityviolation. Then, given that the combination is real due to (2.19), we only need to consider the reality properties of G h n,S (η 1 , η 2 , k). As we did above, we will consider light and heavy fields separately.
Light fields First consider light fields, i.e. those in the complementary series, and the modes with n = |h| where the solution to the homogeneous equation of motion is given by (2.39). Given that this solution involves a single Hankel function, the properties of G h |h|,S (η 1 , η 2 , k) will be very similar to what we encountered in the CCM case but withκ = 0. Indeed, for light fields we therefore already know that this propagator is purely real after Wick rotation. In any case let us show this explicitly, following [37], since this will allow us to easily see how to extend the proof to the other modes with the same helicity, but with n > |h|. We have where as usual the term on the second line is there to ensure that we satisfy the future boundary condition of the bulk-bulk propagator, and the η 0 dependence drops out of this term since we are considering light fields. It is simple to see that the second term ensures that we satisfy Dirichlet boundary conditions since for light fields H . We can write this expression more compactly as where J ν S (z) is the Bessel function of the first kind. We can now make use of the following integral representations: which are respectively valid for −π < ph z < 0 and Re ν S > − 1 2 , to conclude that where we have used Im ν S = 0. We also have We therefore conclude that the Wick-rotated propagator is purely real i.e.
Given this result we would immediately expect this property to hold for the other modes too since ultimately they form a single multiplet. Let's see this explicitly by considering modes with the same helicity but with n > |h|. The general form of the bulk-bulk propagator is whereD h,n (x, k) are real differential operators in x with k-dependent coefficients. It follows thatD h,n (iη, k) are purely real after Wick rotation. Furthermore, they are either purely even in x and therefore purely real when written in terms of η (if n − |h| is even) or purely odd in x and therefore purely imaginary when written in terms of η (if n − |h| is odd). This final property follows from the fact that the B m,n in (2.36) are real and only non-zero when n and m differ by an even number. We can then write We can then use the integral representations above and the fact thatD h,n (iη, k) is always real after Wick rotation to conclude that this bulk-bulk propagator is real after Wick rotation. We have therefore shown that for light fields 45) and since this holds for all 0 ⩽ n ⩽ S we conclude that It is easy to see that the same reality condition holds for partially-massless fields in the discrete series, since the proof above is valid for arbitrary ν S > 0 as long as the pure gauge modes are excluded.
Heavy fields Now consider heavy fields i.e. those in the principle series, with ν S = iµ S . As we have seen in the CCM scenario, the full bulk-bulk propagator will not be real after Wick rotation, so instead we add and subtract factorised terms such that the connected part of the bulk-bulk propagator is real after rotation. As we did above, we work with G i1···inj1···jn (η 1 , η 2 , k) and work helicity-by-helicity. For each mode we define the following decomposition of the bulk-bulk propagator: As before we take ∆G h n,S (η 1 , η 2 , k) to be factorised, to solve the homogeneous equation of motion, and to vanish in the far past. We therefore write For the n = |h| cases, where the mode function is given by a single Hankel function, we can easily read off the necessary form of A h,h since it is a special case of what we did in the CCM scenario withκ = 0. We derived the constraint that A h,h must satisfy, and the solution, in Appendix A. The result is Note that we could also add any purely real correction to A h,|h| and still satisfy the necessary condition so this choice is the minimal one required to makeC h |h|,S (χ 1 , χ 2 , k) real. With this result in hand, it is simple to deduce what we need to add for modes with n > |h|. We have and by using the fact thatD h,n (iη, k) is purely real for even n − |h|, and purely imaginary for odd n − |h|, we can write Given thatD h,n (iη, k) is always real after Wick rotation, we therefore conclude that the choice ensures that (3.51) is satisfied. Again we can add any purely real correction to A h,n (µ S ), but this is the minimal solution that is sufficient for our purposes. We therefore conclude that 56) and since this holds for all 0 ⩽ n ⩽ S we conclude that We have therefore shown that for the cosmological collider physics set-up, the full bulk-bulk propagator G is real for light fields, while for heavy fields this is a property of the connected partC only. We will use these properties in the next section to prove the reality of total-energy poles in wavefunction coefficients with massless external states, and we useF in Section 5 to compute exact parity-odd trispectra.

Reality and factorisation
In this section we unleash the full power of the reality properties of the indexed propagators for general massive fields that we derived in the previous section, and prove a theorem revealing the universal reality of the total-energy singularities in any tree diagrams with external massless scalars (and massless gravitons) for both the CCM and the CC scenarios. More precisely, we will see that in theories involving light fields of arbitrary mass, spin, couplings and chemical potential, the wavefunction coefficient of external massless scalars is always real. This property also holds for an even number of external conformally coupled scalars, and also external massless gravitons as long as we sum over the two helicities. In more general theories involving heavy fields, despite the fact that the full wavefunction coefficient can be complex, the totalenergy singularities remain real. Indeed, these singularities come from the maximally-connected parts of the wavefunction coefficients (which we define below) and we prove that these parts are purely real. Based on the universal reality of total-energy singularities, we then prove that all parity-odd correlators are necessarily factorised at tree-level, and are free of any total-energy singularities, under the assumption of IR-convergence.

Light fields: the wavefunction reality
Cosmological condensed matter We start by considering the most general theory of interacting light fields in the CCM scenario. The particle spectrum L consists of a set of spinning fields σ f i1···i S f labelled by their flavour f ∈ L , and each field σ f i1···i S f is equipped with an integer spin S f = 0, 1, 2, · · · , a dimensionless mass parameter 0 < ν f ⩽ 3/2, a sound speed c f , and a chemical potential κ f . Among these fields, we will mainly pay attention to the statistics of a "visible" massless scalar field ϕ ≡ σ f =ϕ with S ϕ = 0 and ν ϕ = 3/2, while allowing for exchanges of any of the other fields. For instance, in inflationary cosmology, ϕ is usually associated with the curvature perturbation ζ (or equivalently, the Goldstone π of broken time translation), while σ f i1···i S f , f ̸ = ϕ are associated with various isocurvature perturbations. In position space, the most general interaction Lagrangian at a vertex v schematically reads where D v collectively denotes the scale factors and derivatives at the vertex v, and Here k v , l v ⩾ 0 are integers counting the number of derivatives and p v , q v , N v,f ⩾ 0 are integers counting the spatial indices and number of fields involved in the interaction vertex v. 11 The number of scale factors is fixed by diffeomorphism invariance (or its global subgroup of de Sitter scale invariance). Note also that the coupling coefficients are real constants in time i.e.
as a consequence of unitarity and de Sitter scale invariance.
Let us now try to compute the wavefunction coefficient ψ n for n external massless scalars ϕ. A general tree diagram that contributes to the wavefunction coefficient ψ n is obtained via contracting the fields in adjacent vertices to form bulk-bulk propagators G f i1···i S f j1···j S f and massless bulk-boundary propagators K ϕ , before finally integrating the interaction time η v at the vertices. The bulk-boundary propagator for the massless scalar is which is crucially purely real after Wick rotation. Schematically, for a tree diagram with n external lines, V interaction vertices and I internal lines, we have As dictated by the Feynman rules, we have included a factor of i for each vertex. Now notice that the whole integrand of (4.4) is analytic in the second quadrant of complex η-plane, 12 which allows us to deform the integration contour by performing a Wick rotation under which the propagators become The original integration contour along the negative real axis is thus deformed to that along the positive imaginary axis, together with an arc at infinity. The iϵ-prescription and the Bunch-Davies initial condition for the propagators guarantee that as |η| → ∞, the integrand of (4.4) decays exponentially, leaving a vanishing contribution from the arc at infinity. The only non-vanishing contribution now comes from the integral along the positive imaginary axis where χ > 0. In momentum space, the vertex derivative operators transform as Thus the wavefunction coefficient becomes (4.9) Now we evoke the reality property (3.22) of the Wick-rotated bulk-bulk propagator for light fields, namelỹ together with the reality of the Wick-rotated derivative operator, 12) to see that each individual factor in (4.9) is purely real. Therefore, combining (4.2), (4.12), (4.10) and (4.11), we finally arrive at ψ * n = ψ n . (4.13) More precisely speaking, however, what we have managed to prove is only the reality of the integrand in the perturbative expression of ψ n . In order to make the final logical leap to the reality of the full-fledged ψ n , we need to also ensure the convergence of the integral. Since the propagators and the vertices are well-behaved for any finite Euclidean conformal time 0 < χ < ∞, we only need to check the convergence at the endpoints. In particular, the UV convergence at χ → ∞ is guaranteed by the time-ordering and the Bunch-Davies initial condition, as mentioned above. On the other hand, the IR convergence at χ → 0 is not automatic and is generally model-dependent. For instance, a λσ 4 self-interaction for a massless σ field can bring a logarithmic divergence ln(−η) in the IR limit, which after Wick rotation (4.5) brings an imaginary factor of ln e −i(π/2) = −iπ/2. An odd number of such λσ 4 insertions would lead to a complex ψ n in general. One can also view this as a consequence of scale invariance spontaneously broken by the IR cutoff [37]. Therefore, we will restrict ourselves to the case of IR-convergent interactions and perfect scale invariance. This establishes the reality of ψ n for a general tree diagram expressed by (4.9). Since the total ψ n is a sum of all possible diagrams with n external ϕ lines, the reality property extends to the full ψ n at tree-level for the CCM scenario.
Cosmological collider This proof easily generalises to the CC scenario, where the dimensionless mass of the spinning fields (2.41) takes 0 < ν f < 1/2 for the complementary series and ν f = 1/2 + T, T ∈ N for the discrete series. The bulk-bulk propagators are labelled with covariant indices, (e.g. G f µ1···µ S f ,ν1···ν S f ) and the vertices include contractions of all the covariant indices including that of a time-like unit vector n µ , which allows us to break de Sitter boosts at the level of the interactions as in [5]. The Lagrangian vertices are then Here, ε µνρσ = √ −gϵ µνρσ is the covariant Levi-Civita tensor density. For simplicity, we have introduced rescaled quantities that are related to the original covariant ones byn µ = a −1 n µ andΦ f Note that these rescaled quantities are contracted with the flat metric tensor η µν , and O(∂ lv−1 ) are terms with fewer derivatives that take an analogous form to D v . 13 Henceforth, it suffices to consider the leading operator λ v D vΦ Nv as a general parametrisation of the interactions. In other words, we have expanded all the covariant interaction operators and re-organised the expansion in powers of ordinary derivatives on fields.
The computation of ψ n at tree-level is completely analogous to (4.4), with D v and G e ′ replaced by their covariant cousins. We have shown in (3.46) that the indexed propagator of Φ f satisfies the reality condition Thus the propagator of the rescaled fieldΦ f also satisfies the reality condition after Wick rotation (4. 16) In addition, the vertex derivative operators transforms as with the same reality property as in (4.12) Thus all the ingredients in the perturbative computation of ψ n transform identically as in the CCM scenario, leading to the same conclusion: Theorem 4.1. (ψ n -reality) The tree-level wavefunction coefficient of massless scalar fields is purely real, i.e. Im ψ n = 0, in theories containing an arbitrary number of fields of any light mass, spin, coupling, sound speed and chemical potential, under the assumption of locality, unitarity, scale invariance, IR convergence and a Bunch-Davies vacuum.
Discussion Before we move on to the inclusion of heavy fields, let us make a few remarks: • Although we have chosen a specific massless scalar field f = ϕ as the visible sector, and focused on its wavefunction coefficient ψ n , the same proof straightforwardly generalises to multiple massless scalar fields with different flavours, since their bulk-boundary propagators are all real after Wick rotation, and we did not use any Bose symmetry properties in the proof.
• Apart from massless scalar external lines, the particle spectrum may include massless spinning fields such as the graviton. Since the external polarisation tensors are complex, the helical wavefunction coefficients are in general complex. Thus it is more convenient to work with the indexed wavefunction coefficients where the helicities are added together. The indexed bulk-boundary propagator of the massless graviton is where • We can also consider conformally-coupled scalars φ on the external lines where the bulk-boundary propagators are where the η 0 -dependence ensures that we satisfy the future boundary condition. Clearly this propagator is not real after Wick rotation, instead it is purely imaginary. Our reality theorem then extends to these fields if we have an even number of them on external lines.
• The ψ n -reality automatically implies the reality of the total-energy k T -singularities within. As we shall see in the next subsection, after including heavy fields, the wavefunction coefficient can become complex, but the k T -reality remains true.

Adding heavy fields: the total-energy reality
Now let us move on to the case with heavy fields included. We again label the heavy fields by their flavor f ∈ H , with a dimensionless mass µ f = −iν f > 0. In the CCM scenario, the most general interaction vertex straightforwardly generalises to (4.24) In contrast to light fields, the Wick-rotated bulk-bulk propagator of heavy fields does not enjoy the reality property by itself. However, as we showed in Section 3, we can always achieve reality for the connected part of the propagator by adding appropriate solutions of the homogeneous equation of motion, regardless of the mass of the propagating field. This leads us to a decomposition of the indexed bulk-bulk propagator of any mass, and as shown in (3.29), after Wick rotation, the connected part enjoys the reality property Now the key insight is that the time-ordering θ-functions only appear in the connected part C, and the factorised part F is a sum of products of functions of the vertex times. Therefore, in a general tree diagram, one can extract the maximally-connected contribution by isolating the all-C piece after the decomposition (4.25), where we have explicitly spelled out the k T dependence in ψ C n , with (· · · ) denoting other kinematic variables. For n = 4 (where we only have a single bulk-bulk propagator) the factorised part here is completely factorised in the sense that there are no θ-functions, while for higher-point coefficients the factorised part can still contain θ-functions but crucially fewer than those contained in the all-C (maximally-connected) piece. The total-energy singularities only arise from this maximally-connected part, whereas the factorised piece can only have functional dependence on k T that is analytic at k T → 0. In particular, all the total-energy poles are contained in ψ C n , Now given that the connected part of the propagator C enjoys the reality property for both light and heavy fields regardless of their mass, spin, sound speed and chemical potential, we can go through the same proof in the previous subsection, and conclude that the maximally-connected piece must be real: which immediately implies the reality of all the total-energy poles: Im Res The proof in the CC scenario is analogous to that above after doing the same decomposition and utilizing the reality property of the covariant indexed connected propagator ofΦ f based on (3.57), (4.32) Therefore, we conclude with a reality theorem on the k T -poles of tree-level wavefunction coefficients: The maximally-connected piece of a tree-level wavefunction coefficient for massless scalar fields, along with all the total-energy poles therein, is purely real, i.e. Im ψ C n (k T , · · · ) = Im Res k T →0 k m T ψ n = 0 , m, n ∈ N, in theories containing an arbitrary number of fields of any mass, spin, coupling, sound speed and chemical potential, under the assumption of locality, unitarity, scale invariance, IR convergence and a Bunch-Davies vacuum.
Discussion Notice the k T -reality is concretely established only for k T > 0. For k T -poles inside ψ n , the reality of their residues automatically follows from the reality along the positive real axis. However, in a general tree diagram, the k T → 0 limit may possess singularities other than just poles. For instance, a 4-point exchange diagram with a massive field could contain a logarithmic singularity [16,24], which comes with a branch cut with the branching point k T = 0. More generally, one may expect other types of singularities to occur at k T = 0, but none of them can lie along the positive real axis which is physically accessible. Then what we can say is that if such singularities are part of a function analytic along a strip k T > ϵ > 0, then the coefficient of such a function must be real. To illustrate the idea, suppose we have a branch-cut singularity and an essential singularity at k T = 0, then the k T -reality states that α, β, γ ∈ R and Im c α = Im c βγ = 0. Notice that the full analytic structure at k T → 0 should be completely fixed by the perturbative structure at tree-level and the analytic property of the mode functions, and there could be a constraint on the type of allowed total-energy singularities. Hence some of the singularities (for instance, an essential singularity) may not exist at least at tree-level. However, the study of the analytic structure of singularities in the cosmological wavefunction is still in its infancy, and we leave a more detailed analysis to future work. The reality of the k T -singularities (or equivalently, the maximally-connected component ψ C n ) can be physically understood in a few different ways: • First, from the large-mass EFT perspective, one can choose to integrate out the heavy degrees of freedom in the theory and be left with an EFT of light fields only. This procedure is equivalent to performing a large-mass expansion of the heavy field propagators in a given Feynman diagram. After such an expansion, the heavy field propagators are contracted to contact interactions with a tower of derivative couplings that respect scale invariance. Then by the ψ n -reality theorem (4.1), each contracted diagram is purely real since they only involve light fields. On the other hand, the largemass expansion preserves part of the maximally-connected wavefunction ψ C n . This preserved part is free of any partial-energy singularities involving momenta of the expanded heavy field propagators, but includes part of the original total-energy singularities. Hence the reality of these preserved total-energy singularities can be derived from the ψ n -reality after the large-mass expansion, and is consistent with the full k T -reality.
• Second, from the amplitude perspective, the k T -singularities are generated by the time integral of the connected diagram in the past infinity where η → −∞. Namely, we have In the infinite past, the physical momenta of different modes are much larger than their mass scale as well as the Hubble scale, thus all the fields are effectively massless and interacting in a flat spacetime. The limit k T → 0 therefore probes the energy-conserving scattering processes in the analytically continued sense [78]. In a de Sitter-invariant theory, the residue of the leading k T -pole is expected to recover the on-shell (Lorentz-invariant) scattering amplitude of massless particles A n in flat spacetime [78][79][80]. However, it is known that a Lorentz-invariant tree amplitude of scalars in flat spacetime is manifestly real unless one of the internal lines hit the mass shell and become disconnected. 14 This implies that the leading total-energy pole is real. Here we have also shown that sub-leading ones are also real, which can be argued for given the Manifestly Local Test (MLT) of [33] which states that wavefunction coefficients of massless scalars satisfy ∂ψ n ∂k a ka=0 = 0, (4.36) where the derivative with respect to an external energy is taken while holding all other variables fixed. This equation should be satisfied by both contact and exchange diagrams, and is oblivious to the type of state that is being exchanged. The MLT follows from the simple observation that the bulk-boundary propagator of a massless scalar in de Sitter does not contain a term linear in k: This constraint relates sub-leading total-energy poles to leading ones, and given that it is a real constraint, it implies that sub-leading poles are real once the leading ones are. It can be the case that there are sub-leading total-energy poles that are not tied to the leading ones by this constraint. However, in those cases we expect that the leading pole of such terms also has an amplitude interpretation, coming from an interaction with fewer derivatives, since it has a lower-order pole, and then the argument can be run again. 15 We end this subsection by pointing out the k T -reality also applies to external massless gravitons and an even number of conformally coupled scalars, with the argument mirroring that of light fields which we gave above.

Factorising parity-odd correlators
The universal reality of k T -singularities is not only of theoretical interest, but also provides a powerful tool for the computation of phenomenologically interesting parity-odd correlators. We will show in this subsection that all parity-odd correlators of massless fields must be factorised at tree-level and cannot contain k T -singularities as a consequence of the reality theorems we have just derived. The reason why parity is relevant is simple: for any boundary Hermitian operator ϕ † (x) = ϕ(x), its Hermitian conjugate in momentum space is equivalent to a spatial reversal, 14 To see this fact, one simply counts the factors of i in a tree diagram: i V from vertices, (−i) I from off-shell internal lines, i 2n from 2n vertex derivatives, and an overall i from convention. This leads to i I−V +1 = 1 for a tree topology with 0 = I − V + 1. The imaginary part can only come from Im (p 2 + m 2 − iϵ) −1 = πδ(p 2 + m 2 ), where the diagram is factorised on-shell. 15 We thank Austin Joyce for discussions on these points.
If in addition, ϕ(x) is a parity-even scalar (such as the CMB temperature fluctuations), a spatial reversal is equivalent to a parity transformation. Therefore, n-point correlation functions in momentum space can always be decomposed into a parity-even part and a parity-odd part: The Hermiticity of ϕ(x) implies the parity-even part is always real while the parity-odd part is always imaginary, These boundary correlators are computed by a functional integral of the modulus square of the wavefunction, where σ = {σ f i1,···i S f |f ∈ L ∪ H}| η=η0 collectively denotes all the bulk fields evaluated at the future boundary. The wavefunction exponent is organised as a sum over field products 16 k a ϕ(k 1 ) · · · ϕ(k n ) + (· · · ) , (4. 45) where (· · · ) in the second row denotes terms including fields other than the massless scalar ϕ. Due to the modulus square, the phase information of the wavefunction Ψ[σ, η 0 ] is washed out in correlators, which solely depend on the combination in the probability distribution functional The final correlator receives contributions from various partitions of all possible diagrams at tree-level, where q nm , q nml , · · · are real rational numbers counting the combinatorics of partitions, and momentum conservation is implicit in the individual ρ's. Notice that except the first term, all the other contributions in (4.48) are factorised in the sense that they do not contain any k T -singularities. The only non-factorisable contribution that contains k T -singularities comes from the maximally-connected part of the first term, The k T -reality theorem (4.2) tells us that the maximally-connected massless wavefunction is always real, i.e. ψ C * n = ψ C n . This implies i.e. the maximally-connected density matrix is parity-even. We can then check what the consequence of this is for parity-odd n-point correlators. We have We therefore see that parity-odd correlators of massless scalars are factorised and therefore cannot have k T -singularities, This also ensures that parity-odd correlators admit a well-defined Taylor expansion around k T = 0. The proof straightforwardly generalises to the CC scenario, and provides an understanding of why the final parity-odd trispectrum of [104], computed in a non-local EFT with the massive spinning field integrated out, is manifestly factorised (we will discuss this further in Section 5.1). Thus we conclude with the following theorem, Thus the parity-odd n-point correlator sourced by light fields solely receives contributions from lowerpoint wavefunction coefficients ψ f1···fm with m ⩽ n − 1. In contrast, when heavy fields are involved, ψ n is no longer real by itself, and one has to perform the C-F decomposition to isolate the factorised parts of ψ n , which will also contribute to the final parity-odd n-point correlator. Alternatively, one can say that for light fields, F = 0 and there is nothing to isolate away from ψ n . However, such a computational distinction is an artefact of the wavefunction formalism. Namely, the Dirichlet boundary condition at η = η 0 is introduced as an intermediate tool to organise the perturbative expansion. In the absence of IR divergences, the final correlator does not depend on η 0 , and the η 0 dependence in each contributing piece must cancel out. Yet the Dirichlet boundary condition does not respect the continuity of the mass parameter at ν = iµ = 0, since the IR behaviour of light and heavy fields are different: light fields split into two scaling modes with one dominating over the other in the limit η 0 → 0: Heavy fields split into two oscillatory modes with the same damping power, and are equally important as η 0 → 0: Thus implementing a Dirichlet boundary condition at η 0 is sensitive to the mass of the bulk fields, which should be artificial. This is because there is nothing physically problematic 17 at ν = iµ = 0, and all the physical observables such as the boundary correlator should be continuous across this point. As we shall see in Section 5, this is indeed the case for parity-odd 4-point correlators. Namely, the correlator for exchanging a heavy field can be directly obtained via the analytic continuation ν → iµ in the final result of light field exchange. To complement this discussion and the proofs we have outlined in this section, in Appendix B we show how to understand our results in the in-in/Schwinger-Keldysh formalism where the subtlety of the role of η 0 does not appear.

Exact parity-odd trispectra
In this section we present three examples of parity-odd trispectra, which we are able to compute exactly given our theorem that parity-odd correlators are factorised. Indeed, to arrive at these exact shapes we only need to compute time integrals associated with cubic diagrams, without having to worry about the complicated nested integrals that one usually encounters when computing trispectra. Given these examples, it would be straightforward to extend our methods to more general examples corresponding to other interactions, other spins, etc. Before diving into technical details, we first outline the overall algorithm for computing such parityodd trispectra. Based on the C-F decomposition that isolates the connected propagator C which satisfies helical-reality, , we can compute the s-channel parity-odd trispectrum as a sum of three terms, . (5. 2) The first term on the second line vanishes due to the reality of the maximally-connected wavefunction ψ C 4 . We therefore only need to compute the second and third terms which are products of exactly computable cubic time integrals: • For the exchange of light fields we have F = 0, and so the second term on the second line above vanishes, leaving only the third term from the product of cubic wavefunction coefficients: , (light fields).

(5.3)
• For the exchange of heavy fields, the second term is non-zero and so we need to compute both factorised diagrams: , (heavy fields). (5.4) Both diagrams depend on η 0 , and there is an intricate cancellation between the two which ensures that the final result does not depend on this late-time cut-off.
Again we stress that the final results for heavy fields and light fields are related to each other by analytical continuation (as we will see below).

Example 1: spin-1 exchange in CC with chemical potential
One particularly interesting example is the exchange of a massive spin-1 field with a chemical potential. In this subsection we consider this case within the CC scenario (although things are essentially the same in the CCM as we will explain). The chemical potential in this case is naturally introduced via a dimension-5 Chern-Simons coupling ϕFF that respects the shift symmetry ϕ → ϕ + c. In a recent work [104], it has been demonstrated that with the assistance of the chemical potential and a reduced Goldstone sound speed, a large parity-odd trispectrum can be generated. In [104], the massive vector field is integrated out, resulting in a non-local EFT description which is organized as a time-derivative expansion. From the perspective of the single-field effective theory, the parity-odd trispectrum emerges from the (spatially) non-local feature of the self-interactions of the inflaton and this non-locality allows the evasion of the no-go theorems of [37,86]. Utilizing the non-local EFT, one can derive a simple analytical expression that serves as a good approximation to the parity-odd trispectrum and matches well with numerics in certain parameter regimes. In contrast, in this work, we will set out to obtain the exact result of the parity-odd trispectrum using our factorisation theorem. 18 The action of a massive vector field with a chemical potential is −g is the contravariant Levi-Civita tensor density. The mass term breaks the U (1) gauge symmetry of the spin-1 field. As we will explain in more detail below, all of our results also apply in the massless limit corresponding to axion-U (1) gauge field inflation, see e.g. [100,[114][115][116][117][118][119][120][121][122][123][124], however throughout this subsection we will be more general and keep the mass term.
The inflaton background can be expanded as ϕ = const +φ 0 t, where higher order terms are suppressed by slow-roll parameters. The constant term has no dynamical effects due to the shift symmetry of this dimension-5 operator, and the chemical potential κ is equal toφ 0 /Λ c . The equation of motion of this spin-1 field is then The second term can be eliminate using and taking the divergence of both sides yields the transverse constraint ∇ ν Φ ν = 0. The final equation of motion is then We now convert to momentum space and decompose into the different helicities. We write The equations of motion then decouple for each mode, and only the transverse mode will be affected by the addition of the chemical potential, while the temporal and longitudinal mode remain the same as in Section 2. Since we ultimately care about the parity-odd contributions, which cannot come from the exchange of h = 0 modes, let us only focus on the transverse modes which are subject to 12) and the solution to this equation with Bunch-Davies vacuum conditions is given by the Whittaker-W function: withκ ≡ hκ/H. This is familiar from our discussion of the cosmological condensed matter scenario but now the mass parameter is different, as is the scaling dimension of the field. Here we have Given that only helicity states with ±1 are relevant here (and they have the same speed), we have set the speed of sound of the internal massive field to unity, and incorporated a dependence on the speed of sound of the external Goldstone boson c s . This convention has been adopted in [104] and will make the comparison between the two sets of results more transparent.
We now need to choose interaction vertices of the form ππΦ, and in order to make use of our factorisation theorem the interactions need to be IR-finite. EFToI operators that are quadratic or cubic in building blocks can both yield the desired interactions where here we define a building block as an object that starts at linear order in fluctuations. By using only these operators the tadpole cancellation is guaranteed. For operators that are quadratic in building blocks the presence of ππΦ couplings also induces πΦ couplings. Such couplings have two primary effects: they can contribute to the bispectrum of curvature perturbations through single-exchange diagrams, and yield new trispectrum diagrams which perturbatively capture the corrections to the linear theory. See [34,35] for recent works bootstrapping such single-exchange contributions to the bispectrum. Furthermore, for these interactions the time integrals that we need to compute are not IR-convergent and so we cannot conclude that the parity-odd correlator is factorised. 19 In our quest to write down exact shapes, we will therefore concentrate on EFToI operators for which the leading vertices are ππΦ i.e. those that do not necessarily come with πΦ couplings by symmetry. This essentially tells us that π must appear as π ′ or ∂ i ∂ j π (which come from the EFToI operators δg 00 and δK µν ), and we can add extra derivatives to these objects. By requiring that π always appears in this way we ensure convergence of all time integrals at the conformal boundary. Indeed, the two spatial derivatives come with two factors of η by scale invariance while π ′ yields one power of η by scale invariance but then we have which yields an additional factor. The net contribution from the two factors of π in each vertex is then at least four powers of η which cancels the four inverse powers coming from the integration measure. Adding additional derivatives can only improve convergence thanks to the additional powers of η which are dictated by scale invariance. We can also restrict to at most one time derivative on each of the external bulk-boundary propagators given that higher order ones can be eliminated by the scalar field's equation of motion. The lowest dimension operators which satisfies these properties are dimension-7 [37], and for concreteness we will use 20 which originates from the EFT operator ∇ µ δg 00 δK µ ν Φ ν . The number of scale factors can be understood from scale invariance given that under a scale transformation the vector transforms in the same way as a(η)π c (since it is spin-1). Here π c is the canonically normalized Goldstone boson π c = c −3/2 s f 2 π π with f 4 π = H 4 /(2π∆ ζ ) 2 , and ∆ 2 ζ ≈ 2 × 10 −9 is the observed dimensionless power spectrum. In the following we compute the parity-odd trispectrum due to the exchange of this massive vector due to (5.16). We consider light and heavy fields separately. 19 In general IR-convergence depends on the mass of the exchanged field since its late-time behaviour in time depends on the mass. Here we will assume convergence in the limit that the exchanged field is massless which then guarantees convergence when it is massive. 20 There is a dimension-6 operator, π ′ ∂ i π ′ Φ i , that satisfies our requirements but since only the h = ±1 modes contribute we can take Φ i to be transverse and then this operator is a total spatial derivative. The resulting correlator will then vanish once we impose momentum conservation. The other dimension-7 operator that we could use is ∂ 2 π∂ i π ′ Φ i . The correlator arising from this vertex will only differ from the one we are going to compute in the kinematic factors since the time integrals will be the same.
Light mass case. Let us first consider the light mass case where m < H/2. We remind the reader that the s-channel contribution to the trispectrum is, c.f. (4.48), According to the proof in Section 4, the full ψ 4 is real and ρ 4 does not contribute to the parity-odd trispectrum (see (5.3)). Here we therefore directly compute the factorised contributions i.e. the cubic wavefunction coefficients. We have Here the superscript denotes the helicity of the external vector field, and we have used The dynamical integral can be evaluated exactly using the Laplace transformation of the Whittaker function, where 2F1 is the regularized hypergeometric function: This integral enjoys the helical reality property we have seen many times above, The cubic wavefunction coefficient is therefore given by (5.23) and the cubic density matrix coefficient reads The product of polarisation factors is given by and by putting everything together, projecting onto the parity-odd part of the full trispectrum c.f. (4.43), and transferring π c to the curvature perturbations using the relation ζ = −Hπ = −(2π∆ ζ c 3/2 s /H)π c , we arrive at (s, c s k 12 , ν)I where we have implicitly assumed h = 1 in theκ. The +3 perms yields the correct s-channel trispectrum and then we add the t and u channel permutations to give us the full symmetric trispectrum. To cancel the η 0 factors we used the relation (3.13), and we see that this result is purely imaginary as it should be. Scale invariance tells us that B 4 ∼ k −9 , and we can check that this is indeed the case given that I h n ∼ k 0 . This final result therefore is indeed of the form (1.2), i.e. a sum of terms containing kinematic factors multiplied by a product of hypergeometric functions which come from time evolution. The factor of s · (k 1 × k 3 ) = k 2 · (k 1 × k 3 ) will appear as a factor in all s-channel, parity-odd trispectra.
As a consistency check of this result, and therefore a check of our statement that the quartic wavefunction coefficient is purely real for this light field exchange, we also numerically computed the s-channel trispectrum using the in-in formalism. In Figure 4 we plot the dimensionless trispectrum T s,PO , as defined in (1.15), for both our exact expression and the result from numerics. Clearly the exact solution agrees with the numerics very well as seen from the plot. The mass dependence of T s,PO turns out to be mild and smooth within the complementary series mass range. Interestingly, with a reduced inflaton sound speed c s < 1 and chemical potential κ ≳ H, the trispectrum has a peak roughly lying at a momentum ratio k 1 /s ∼ 2κ/H, which is followed by a dip to zero at k 1 /s ∼ 2κ/(c 1/2 s H) and a second peak at k 1 /s ∼ 2κ/(c s H). The precise location and the relative height of the peaks and the dip depend on the detailed choice of mass and chemical potential, and if probed, can be utilized to break the degeneracy of parameters.
Heavy mass case. We now move to the heavy field case where the mass of the vector Φ is in the principle series. We now have a complex ψ 4 and that will contribute to the final parity-odd trispectrum. However, as we have discussed at length in this paper, the modification of the bulk-bulk propagator by the addition/subtraction of a factorised contribution ensures that the connected component of the wavefunction coefficient is real. Then what we need to take into account are those factorised contributions to the bulk-bulk propagator which we have been denoting by F as a pair of cut dashed lines, as shown diagrammatically by (5.4). For concrete calculations we choose the simplest form of F which is given by with This factorised contribution from F is then given by 29) and the corresponding parity-odd density matrix coefficient is For the light mass case the entire bracket in (5.30) vanishes, which is consistent with ψ 4 being purely real. For the heavy mass case, ρ PO 4 is dependent on η 0 , however this dependence cancels with the contribution to the correlator coming from the cubic wavefunction coefficients thereby rendering the final correlator η 0 -independent (as it should be for IR-finite interactions). The computation of the cubic wavefunction coefficients is identical as above, and is omitted here for brevity. Putting everything together and projecting on to the parity-odd sector, we arrive at × cosh(πκ)A +1,1 I +1 2 (s, c s k 12 , iµ)I +1 2 (s, c s k 34 , iµ) + e πκ Re I +1 2 (s, c s k 12 , iµ)I −1 2 (s, c s k 34 , iµ) − κ → −κ + 3 perms Again this result has the correct momentum scaling, and is purely imaginary. If we compare this result with that of light fields (5.26), we see that they can be converted into each other by replacing iµ ↔ ν. Hence as we expected, there is no discontinuity in the mass parameters. This property extends to other examples too: the heavy field result is always given by an analytic continuation of the light field result. Given that the calculation for light fields is less involved, for the other examples we will concentrate on light fields and then extract the heavy field result via this simple replacement rule. Figure 5: The s-channel dimensionless parity-odd trispectrum Im T s,PO as a function of the momentum ratio k 1 /s. The kinematics is chosen as k 1 = k 3 , k 2 = k 4 = s 2 + k 2 1 and ψ = π/3 being the dihedral angle from the (k 1 , k 2 )-plane to the (k 3 , k 4 )-plane. The parameters are chosen as m = 6H, κ = H (left panel) and m = 6H, κ = 4H (right panel), together with a common sound speed c s = 0.1 and Λ = 30H. The blue and green curves show the leading-order (LO) and the next-to-leading-order (NLO) non-local EFT results. The red curve denotes the exact result (5.31) and the gray dots represent the numerics in the UV theory. We see that the non-local EFT predictions agree with numerics in the small-κ case, yet deviations start to appear at large κ. In principle, such deviations may be cured by adding higher order contributions which are in practice tedious to compute. In contrast, the exact result we have computed in this paper matches the numerics very well for all parameter choices.
Interestingly, for a small inflaton sound speed (i.e. c s ≪ 1), this model of a heavy vector field with chemical potential admits a non-local single-field EFT description in the IR, which well approximates the behaviour of the parity-odd trispectrum in the regime c s κ < c s m < 1 [104]. After partially integrating out the heavy vector, parity violation resurges through emergent non-locality in the effective vertex, and the resulting parity-odd trispectrum is neatly computed as a residue of the non-local pole in the effective vertex. Such a miraculous behaviour as seen from the non-local EFT now becomes understandable from the parity-factorisation perspective. The fact that the parity-odd trispectrum is necessarily factorised in the UV theory is precisely the reason why we only acquire a non-vanishing contribution from the non-local pole in the IR EFT.
To compare our exact result for the parity-odd trispectrum in this model with the non-local EFT prediction, and to check them against numerics, we plot the corresponding dimensionless trispectra in Figure 5. As we can see from the plot, the exact result agrees with numerics very well, while the non-local EFT predictions start to deviate from the exact result when the chemical potential κ is large.
Before moving to some other examples, let us first comment on the massless case with a U (1) gauge symmetry, as promised. In this case we set m = 0 to preserve the gauge symmetry in the free theory of the vector field. We therefore have ν = 1/4. Without adding any additional interactions beyond those in (5.5), a parity-odd trispectrum can be generated at 1-loop due to the πΦΦ coupling. At tree-level we would again need to add interactions of the form ππΦ that preserve the U (1) gauge symmetry. Since the field strength is anti-symmetric, this requires more derivatives than what we have studied so far. Indeed, the first non-zero operator is dimension-8. It would be interesting to study this class of trispectra in more detail.

Example 2: spin-2 exchange in CCM
We now move to a second example where we consider the exchange of a spin-2 field with its dynamics described by the cosmological condensed matter physics scenario. In this case we take the bulk-bulk propagator to be parity-evenκ = 0, and source the parity-violation via the interaction vertices. We therefore need one vertex to have an even number of spatial derivatives, and for the other to have an odd number. As before, we need to ensure IR-convergence. We will therefore work with the following interactions: where the first term is dimension-6 while the second is dimension-7. These are the only operators with those mass dimensions and are the leading ones which are IR-finite. The corresponding EFToI operators are δg 00 δK µν Σ µν and n µ ε µναβ ∇ ν δg 00 δK αγ Σ γ β . In the CCM scenario the conformal weight of the massive field is the same as that of a massless scalar so the counting of the scale factors is simply 4 − (total number of derivatives). As with Section 2 we write 33) and the polarisation tensors are chosen to satisfy conditions (2.6) and (2.7), and are given by The mode functions for each helicity are given by where as we mentioned before we take the chemical potential to vanishκ = 0. In this limit the Whittaker function recovers the Hankel function of the first kind, c.f. (2.15). When the mass of the spin-2 field is light, the only contributions to the parity-odd trispectrum come from the cubic wavefunction coefficients which are given by In the absence of the chemical potential, the I h n integral is identical for each helicity (except for the sound speeds) and purely real. The helicity dependence then resides in the kinematic factors rather than the dynamical ones. The density matrix then reads To make the angular dependence transparent, let us decompose the s-channel trispectrum into two separate parts arising from the exchange of different helicity modes: where we have dropped the contribution from h = 0 since scalar exchanges cannot yield a parity-odd contribution. For the higher helicity modes, we need to fix the polarisation sums. For spin-1 the form of h=±1ê i (k)ê j (−k) can be easily fixed without choosing any particular basis. The result should be parity-even and real given the properties of the polarisation vectors. Scale invariance further constrains it to only depend on δ ij andk i = k i /k. The free parameters can be then fixed by requiring the result to be transverse, and appropriately normalised:ê ± i (k)ê ± i (−k) = 1. We then have For spin-2, where polarisation tensors are combinations ofk i andê ± i , we can proceed in a similar way using (5.33) and (5.39). We have h=±1 e h ij (k)e h mn (−k) =k ikm π jn +k jkm π in +k ikn π jm +k jkn π im , and h=±2 e h ij (k)e h mn (−k) = π im π jn + π in π jm − π ij π mn . (5.41) By combining these polarisation sums evaluated at k = −s and the density matrix coefficients we can extract the contribution from individual helicity exchanges to the final parity-odd trispectrum, I 0 1 (c 1,2 s, k 12 , ν)I 0 2 (c 1,2 s, k 34 , ν) + 7 perms I 0 1 (c 2,2 s, k 12 , ν)I 0 2 (c 2,2 s, k 34 , ν) + 7 perms As a consistency check, we see that the terms proportional to (k 2 · s) (k 4 · s) in (5.42) and (5.43) differ only by a sign and the sound speeds of the different modes. These two contributions would therefore cancel once added together if the sound speeds were identical (c 1,2 = c 2,2 ). In that case the total trispectrum would be independent of s i , which is to be expected since in that case the three polarisation sums add up to an object that is independent of s i . We also see that the result is purely imaginary, and has the correct momentum scaling. For the special cases of ν = 3/2 and ν = 1/2, where the mode functions simplify to exponentials, the trispectrum vanishes which is consistent with the no-go theorem of [37]. Note that here we add +7 perms rather than +3 perms which we had in example 1 since here the two vertices on either side of a diagram are different. It is noteworthy that in de-Sitter/inflationary four-point functions, spin-1 exchange is typically characterised by linear factors of t 2 − u 2 in s-channel diagrams, which originate from contractions between momenta and the polarisation sum. Exchanges of higher spin are then non-linear in t 2 − u 2 . However, here things are slightly different due to the Levi-Civita ϵ-tensor, and it is easy to check no such factor arises for spin-1 exchange as indicated by (5.42). For the exchange of h = ±2 modes, this dependence appears from the k 2 · k 4 factor and its corresponding permutation. Indeed we have For general spin-S, it is simple to see that helicity ±h exchange will introduce contributions to the parityodd trispectrum with a factor of t 2 − u 2 |h|−1 , from which we can read off the spin of the exchanged field.
Here we have computed the trispectrum for light field exchange, and as we discussed above, from this result we can extract that of heavy field exchange by sending ν → iµ. We have also checked by explicit calculation that this analytic continuation yields the correct result.

Example 3: spin-2 exchange in CC
In this final example we again consider spin-2 exchange with parity-violation arising from the interactions vertices, however now we describe the dynamics of the spin-2 field in the cosmological collider physics setup. As we discussed in Section 2, the mode functions in the CC and CCM scenarios can differ significantly, leading to a distinct parity-odd trispectrum compared to what we have just computed in the previous subsection. We will denote the massive spin-2 field by Φ µν . Given that scalar modes do not contribute to the final trispectrum, our focus will be on the components Φ 0j and Φ ij . Again we need IR-finite interaction vertices and we choose which arise from the following EFToI operators: Again the scale factors are fixed by scale invariance (note that here Φ scales in the same way as a 2 (η)π c ). We now decompose the field into the helicity basis: and from now on we will ignore the longitudinal modes Φ 0 1,2 and Φ 0 2,2 since they will not contribute to the parity-odd trispectrum. This leaves us with three modes: Φ ±1 1,2 , Φ ±1 2,2 , and Φ ±2 2,2 . The polarization tensors for these modes, which satisfy (2.19) and (2.20), are The factor √ 2 arises because we adhere to the convention of [5], where e (±1) i e * (±1) i = 2. Following our discussion in Section 2, we can derive the following mode functions where Z |h| s is given by (2.40), and in the last line of Φ ±1 2,2 we have used the recursion relation of Hankel function to simplify the expression. As in the previous example, we consider three contributions to the final parity-odd trispectrum arising from the exchange of these three modes. We denote the different components as B ζ (n,h) , and therefore the full trispectrum is The computation of each component is the same as what we have been through above. Once again, considering light fields buys us the privilege of limiting the focus on the cubic wavefunction coefficients only (i.e. (5.3)), before obtaining the density matrix coefficients and summing over helicities. The result for heavy fields then follows from analytic continuation. In short, we find I 0 3 (s, k 12 , ν)I 0 4 (s, k 34 , ν) + 7 perms Again we see that each component is purely imaginary and has the correct scaling with momenta (in checking this point we use the fact that Z h 2 (s) ∼ s 1/2 ). We also see that for the special case of ν = 1/2 the trispectrum vanishes. This limit corresponds to the partially-massless limit where the massive spin-2 field has only four propagating degrees of freedom (since the h = 0 modes do not contribute here we don't need to worry about the divergence of the corresponding two-point function in this PM limit). In this limit the mode functions simplify to exponentials and again the no-go theorem of [37] dictates that the result should vanish.

Conclusions and outlook
In this paper we have derived some reality theorems relating to cosmological wavefunction coefficients of massless scalars, massless gravitons and conformally coupled scalars. We have shown that the maximallyconnected part of these coefficients, which is i) the most difficult to compute since it contains the maximum number of nested time integrals and ii) the part that can be singular as the total-energy goes to zero k T → 0, is a purely real function of the kinematics. Our results allow for the exchange of states with any mass and integer spin, and in deriving our results we considered two distinct descriptions for the dynamics of massive spinning fields during inflation: cosmological condensed matter physics (where states are representations of the group of rotations) and cosmological collider physics (where states are representations of the de Sitter group). Furthermore, if all exchanged fields are in the complementary series i.e. they have light masses, then our reality theorem extends beyond the maximally-connected part to the full wavefunction coefficient. Our results apply under the following assumptions: • Tree-level approximation: we considered tree-level Feynman diagrams which allowed us to avoid having to analytically continue the spacetime dimensions. This meant that we could work with fixed external propagators which have simple properties in D = 3+1, namely that they are real after Wick rotation. In Appendix C we offer an alternative proof of our theorems using Hermitian analyticity of all propagators, and scale invariance. Here we use the relation V = I + 1 where V is the number of vertices and I is the number of internal propagators. This relation holds only at tree-level. By relaxing the tree-level approximation the maximally-connected parts of the wavefunction coefficients can have imaginary parts [42].
• Bunch-Davies vacuum conditions: this assumption enabled us to rotate all time variables by 90 • in the complex plane as a tool for computing wavefunction coefficients. The fact that we are computing the vacuum wavefunction for which fields vanish exponentially fast in the far past allowed us to close the contour and drop any contributions from the arc such that integration along the real line could be replaced by integration along the imaginary line. This assumption is relevant for our proof in Appendix C since Hermitian analyticity is closely tied to having Bunch-Davies vacuum conditions. If we relax this assumption, for example with the Ghost Condensate, then the maximallyconnected parts of the wavefunction coefficients can have imaginary parts [37].
• Scale invariance: this assumption ensured that the vertex operators were real after Wick rotation. Indeed, scale invariance ensures that time derivatives enter as η∂ η and spatial momenta enter as iηk. This could then be combined with the reality properties of the propagators to prove that various integrands are purely real. If we relax this assumption, for example by allowing for time-dependent couplings or going to general FLRW spacetimes (except when the scale factor is an odd-power-law function of the conformal time), then the maximally-connected parts of the wavefunction coefficients can have imaginary parts [37,125].
• IR-convergence: this assumption enabled us to make the leap from proving the reality of integrands to the reality of the integrated result. Indeed, IR-convergence meant that the final result was independent of η 0 and so we did not need to worry about rotating this cut-off. In the presence of an IR-divergence we would need to also rotate η 0 and this can yield imaginary parts with a logarithmicdivergence, for example. In Appendix C we combined scale invariance of the bulk interactions with IR-convergence to use the fact that the wavefunction coefficients have a fixed scaling with momenta yielding simple transformation properties as all momenta and energies flip sign. If we relax this assumption, for example by allowing for IR-divergent bulk interactions as occurs for the minimal coupling between the inflaton and the massless graviton, then the maximally-connected parts of the wavefunction coefficients can have imaginary parts [101].
The reality of the maximally-connected part of wavefunction coefficients is not just of theoretical interest, rather it can make the computation of phenomenologically relevant ( [97][98][99]102,103]) cosmological correlators a far simpler task than naively expected. We have shown this in Section 5 by considering the parity-odd trispectra of curvature perturbations due to the coupling of the inflaton with another sector with massive spinning fields and parity-violation. Since parity-odd correlators depend on the imaginary part of the maximally-connected wavefunction coefficients, these trispectra are factorised and computed by considering cubic diagrams only. We presented a number examples considering both the CCM and CC scenarios, both light and heavy fields (with the final answers related by analytic continuation), and both parity-violation arising from the free theory of the massive spinning fields and from the bulk interactions. In particular, we considered a parity-violating correction to the action of a massive vector field in Section 5.1 and compared our result with that computed in [104] using a non-local EFT arising from integrating out the massive vector field. Our result recovers the EFT result in the appropriate limit, but also gives an exact result in the regime where the EFT breaks down. This example also includes axion-U (1) gauge field inflation. We also considered examples with massive spin-2 fields in Sections 5.2 and 5.3. For all spins we allowed for a chemical potential in our analysis, and for hierarchies between the speed of the inflaton and the exchanged fields, such that the size of trispectra can be enhanced.
There are many avenues for future research directions and here we outline a few: • Moving on to loops: it would be interesting to see if any of our reality theorems hold at loop level. As we mentioned above, total-energy poles coming from loops can be imaginary, but perhaps the structure of such imaginary terms can be constrained given that the reality properties of bulk-bulk propagators still hold. In fact, in the original D = 3 + 1 spacetime dimensions, the loop integrand for ψ n of light fields still appears to be purely real after Wick rotation, but the dimensional regulator demands an evaluation in D = 4 − ϵ dimensions. The Wick rotation in such a case is expected to bring factors of (−η) ϵ = (−iχ) ϵ = (1−iπϵ/2)χ ϵ . The O(ϵ) imaginary part turns out to be cancellable by the 1/ϵ divergences, giving finite imaginary contributions to ψ n . This has been demonstrated for massless and conformally coupled theories in [42]. It would be interesting to see if such a phenomenon persists for general massive fields and at higher loops.
• Distinguishing between massive spinning field set-ups: we have considered two different descriptions of massive spinning fields during inflation (CC and CCM). It would be interesting to investigate if these two set-ups could be distinguished from each other at the level of massless scalar and massless graviton correlation functions.
• Kramers-Kronig for correlators: it would be interesting to investigate if the parity-even part of a correlator can be constrained given the parity-odd part. This is conceivable for examples where the parity-violation is driven by a chemical potential correction to the free theory. Perhaps consistency of higher-point functions could be used to constrain lower-point ones in this regard. Since the parityodd part is always imaginary and the parity-even part is always real, such reconstruction of a full correlator from its imaginary part, if possible, will be an interesting cosmological analogy of the well-celebrated Kramers-Kronig relations in electromagnetism.
• EAdS perspective: the Wick rotation used extensively throughout this work is in fact a contour deformation, i.e. a change of integrated bulk time which are dummy variables. However it is tempting to conceive a further Wick rotation for the boundary time η 0 = iχ 0 which is explicit. In doing so, one arrives at a quantity different from the wavefunction of the universe. The reality of this quantity is especially transparent since the whole spacetime becomes Euclidean. In fact, a further rotation of the Hubble parameter H = −i/L AdS yields a theory defined in Euclidean Anti-de Sitter (EAdS) space. The wavefunction then becomes the partition function of the boundary CFT of the EAdS bulk [71]. The difficulty, however, is in the fact that the continuations of η 0 and H do not seem to commute with the recipe of obtaining correlators from wavefunction coefficients, which relies on the unitarity of the wavefunction. It would be helpful to understand how to correctly perform this continuation, i.e. how to understand in-in correlators in Euclidean field theories.
• Classifying singularities: we have concentrated on the total-energy singularities in this work, yet they are by no means the only singularities of the wavefunction. Notably, there are also partialenergy singularities when the total energy flowing into a sub-component of a Feynman diagram approaches zero. For IR-convergent massless or conformally coupled theories in de Sitter and flat spacetime, these singularities are always poles at tree-level since the wavefunction coefficients are rational functions of momenta [28,78]. However, in more general theories with arbitrary mass, spin, sound speed and chemical potential, these singularities have not yet been classified even at treelevel. The classification of these singularities would crucially serve as a first step toward a complete understanding of the analytic structure of wavefunction of the universe. is helically real after Wick rotation. We allow for a chemical potential κ in the massive field's dispersion relation in the CCM scenario such that the relevant mode function is given by a Whittaker function, c.f. (2.14). Recall that we need to decompose the bulk-bulk propagator into two parts, In order to ensure that the connected part still satisfies the propagator equation (3.7), the added term must satisfy the homogeneous equation The UV convergence of bulk time integrals requires the Bunch-Davies initial condition at η → −∞, namely Note that in contrast to (3.8), we do not impose any boundary condition at η 0 for ∆G (h) . Upon symmetrizing over η 1 ↔ η 2 , we are left with ) is a helicity-dependent constant to be determined. Now we demand that this connected propagator is helically real after rotation, i.e.
Due to the symmetry between χ 1 ↔ χ 2 , we are free to pick χ 1 < χ 2 without loss of generality. Thus the Wick-rotated connected propagator becomes Notice that Whittaker functions have a branch cut along the negative real axis, which is why we have kept the e iϵ factor to ensure that the branch cut is never crossed. However, the complex conjugation of Whittaker functions is most transparent when their argument lies along the positive real axis, To inspect the complex conjugation [C (h) σ ] * , we can first expand the Whittaker W -functions in terms of Whittaker M -functions, and then rotate away arguments lying below the branch cut using This concludes the construction of the desired connected propagator in the CCM scenario.
In the CC scenario, we turn off the chemical potential, and the general constraint (A.15) for the mode with n = |h| reduces to The simplest choice is to demand A −h,|h| = A h,|h| = −A * h,|h| , which gives Note that this solution establishes the reality of C h |h|,S , while that of C h n,S , n > |h| follows trivially from acting with real derivative operators on C h |h|,S , as we explained in the main text. Finally, we point out that although the inclusion of heavy fields motivated the ∆G piece that achieves the reality of C, the same procedure works equally well for light fields except when ν is a half-integer. It is easy to see that all the derivations above straightforwardly go through with the replacement of µ → −iν. In fact, after choosing the specific choice (A.17) and replacing µ by −iν, the connected propagator becomes identical to the original bulk-bulk propagator in the η 0 → 0 limit (c.f., (3.13)), i.e. C = G and F = 0. For the case of ν = n/2, n ∈ N + , the equation (A.9) reaches a singularity, rendering the proof invalid. However, it is easy to check the validity of (A.17) by inserting it into (A.6), using the same derivation that appeared in Section 3.1. This crucially demonstrates that the proof of k T -reality does not make any artificial distinction between light fields and heavy fields; we can extract the connected propagator part of all the internal lines regardless of their masses, and make use of their reality property to show the total-energy poles are always real.

B Proofs via the in-in/Schwinger-Keldysh formalism
In this appendix we streamline the derivation of the reality and parity-odd factorisation theorems in the language of the more conventional in-in and Schwinger-Keldysh formalisms. Both of them focus directly on the observable, i.e. the n-point correlation function itself, without introducing intermediate quantities such as the wavefunction representation of the quantum state of the universe. Therefore, we shall translate the k T -reality theorem for wavefunction coefficients (4.2) to the correlator level, and show that consequently all parity-odd correlators are factorised at tree-level. Here we will only deal with the CCM case, and the proof straightforwardly extends to the CC case as shown in the main text.
The perturbative computation of correlators can be organized using diagrammatics with slightly different Feynman rules from those of the cosmological wavefunction (see, for example, [8]). In short, every n-point correlator is computed by a set of diagrams with coloured vertices indicating whether they are time-ordered (black or "+") or anti-time-ordered (white or "−"). Thus in general, a diagram of V vertices will comprise of 2 V coloured copies which need to be summed. A change from a black vertex to a white one (i.e. from + to −) corresponds to a local complex conjugation plus a flip of the direction of all the 3-momenta flowing into the vertex. The internal propagators connecting these vertices are thus classified into four types according to the colour of their vertices: The external propagators are simply obtained from sending one of the vertices to the boundary in the internal propagators, All of these propagators satisfy the Bunch-Davies (or anti-Bunch-Davies for anti-time-ordered vertices) initial condition in the far past, while no boundary condition at η 0 is posed for them. Instead, they satisfy the conjugation rule Notice that the flip of vertex colour is accompanied by a flip of momentum in the kinematic structure in addition to the complex conjugation. We start by noticing that the total-energy singularities can only reside in monochromatic diagrams where the vertices are either all-black or all-white. This is simply a consequence of the factorised nature of the Wightman functions G ±∓ (η 1 , η 2 , k), i.e. any polychromatic diagram is necessarily disconnected in time at the internal line with opposite colours. In the wavefunction formalism language, they correspond to contributions from the factorised third term in the bulk-bulk propagator (1.7) together with sewing disconnected ρ n in (4.48). Since the all-white diagram is just the complex conjugation plus momentum inversion of the all-black diagram, we will focus on the all-black diagram without loss of generality, where D v is given by (4.1) for the CCM scenario. As with the derivation for the cosmological wavefunction ψ n , we perform a Wick rotation η = iχ, The solution of ∆G is identical to that in the wavefunction approach (see Appendix A). Thus after expanding the internal propagators of the all-black diagram, we deduce that the maximally-connected contribution is purely real (see the diagrammatic illustration in Figure 6). Consequently, all the total-energy poles inside the full tree-level correlator must also be real, Other total-energy singularities are also real when understood as analytically continued from the positive- Figure 6: A five-point illustration of the k T -reality in the in-in/Schwinger-Keldysh diagrammatics. The all-black vertices indicate the diagram is a fully time-ordered diagram ⟨φ 5 ⟩ + . We then expand the internal Schwinger-Keldysh propagators (solid lines) into the connected (double lines) and factorized parts (dashed lines), and use the reality of the connected propagator to conclude the reality of the maximally-connected first term and the total-energy singularities therein.
k T direction, see the discussions in Section 4.2. The factorisation of parity-odd correlators then directly follows from the k T -reality, where we have applied (B.5) and (B.13). Finally, we comment that in the in-in/Schwinger-Keldysh formalism, there is no explicit reference to the boundary time η 0 except for the external propagator K ± which can be trivialised by sending η 0 → 0 first, and consequently the parity-odd correlators factorise in the same fashion for light fields and heavy fields. This generalises the proof of [37], which applied to massless and conformally-coupled mode functions, allowing for general massive mode functions for the internal lines.

C Reality from Hermitian analyticity
In [37], the Cosmological Optical Theorem (COT) of [22] was used to deduce that contact diagrams of massless scalars arising from IR-finite interactions are purely real, which in turn implies that such diagrams do not contribute to parity-odd trispectra. This result itself suggests that such a trispectrum is a very nice probe of exotic inflationary physics. The same result was also derived in [86] using Wick rotations. It is therefore tempting to wonder if the more general results we have derived in this paper, namely the inclusion of exchange processes and massive fields, can also be understood using the COT (or more generally without invoking Wick rotations). In other words, can we deduce from the COT that the imaginary part of wavefunction coefficients is factorised? Let us first review the argument for contact diagrams. The COT states that where in the second term we flip both the energies and the momenta. Note that the COT relies on our ability to analytically continue the momenta away from the physical region. The COT follows from having real coupling constants (by unitarity) and Hermitian analyticity of the massless bulk-boundary propagators, K(η, k) = K * (η, −k) c.f. (4.3), and the spatial momenta which enter as ik. We refer the reader to [22][23][24]29] for full details. If we have exact scale invariance then ψ n ∼ k 3 , where the cubic scaling is there to cancel the scaling of the momentum-conserving delta function in three spatial dimensions, and therefore (C.1) becomes Since it is the imaginary part of the wavefunction coefficient that contributes to the parity-odd correlator, c.f. (4.46), this tells us that contact diagrams do not contribute. Here we have used exact scale invariance of the wavefunction coefficients which of course requires scale invariance of the bulk interactions, but also IR-convergence of the time integrals otherwise scale invariance is broken by the IR cut-off η 0 . In that case the wavefunction can indeed have an imaginary part [37]. Now consider a four-point exchange diagram. For an s-channel diagram the COT reads [22]  where the factorised RHS depends on the three-point sub-diagrams that contribute to this four-point coefficient. In addition to unitarity and Hermitian analyticity of the bulk-boundary propagators and spatial momenta, this expression holds since the real part of the bulk-bulk propagator is factorised [24] (we remind the reader that in this paper our bulk-bulk propagator differs by a factor of i from that of [24] which is why the real part rather than the imaginary part is factorised). We can see this explicitly. Indeed, G(η 1 , η 2 , k) = σ(η 1 , k)σ * (η 2 , k)θ(η 1 − η 2 ) + σ(η 2 , k)σ * (η 1 , k)θ(η 2 − η 1 ) + factorised , (C.4) and therefore G(η 1 , η 2 , k) + G * (η 1 , η 2 , k) = [σ(η 1 , k)σ * (η 2 , k) + σ * (η 1 , k)σ(η 2 , k)][θ(η 1 − η 2 ) + θ(η 2 − η 1 )] + factorised = factorised . G e ′ (p e ′ ) , (C.14) where we have used the fact that at tree-level we have V = I + 1. For wavefunction coefficients of massless scalars that adhere to exact scale invariance, i.e. have no η 0 dependence and therefore scale as ψ n ∼ k 3 , we can then write which establishes the reality of ψ n . Here it was crucial that there was no η 0 dependence. If there is an η 0 dependence then flipping the signs of all energies and momenta does not simply yield an overall minus sign since η 0 itself carries a conformal weight. As we have discussed a number times in this paper, the bulk-bulk propagator is independent of η 0 only for light fields, whereas for heavy fields it does indeed depend on η 0 . This proof therefore applies for light fields only, and offers a complementary proof of the result we derived in Section 4 using Wick rotations. For heavy fields this argument does not hold (and indeed we wouldn't expect it to hold since we have already seen that for heavy fields ψ n is not real), but the discussion and the C-F decomposition we made in Section 3 suggests a clear way forward. Indeed, consider the connected bulk-bulk propagator c.f. (3.24), where we have taken the minimal solution of A h which we derived in Appendix A. The relative minus sign between the two terms in the brackets comes from the fact we have written this propagator in terms of σ − and σ + rather than σ and σ * . Here we have written the solution for A h that is valid for both light and heavy fields. We can then use the Hermitian analytic properties of the mode functions to conclude that this connected bulk-bulk propagator is anti-Hermitian analytic: which holds for both light and heavy fields. In Appendix A we saw that we could add any real term to A h while maintaining the reality of the connected bulk-bulk propagator after Wick rotation. Here we can also add any real term to A h and still realise the anti-Hermitian analyticity given (C.10). We can now run the same argument as above but with the full bulk-bulk propagator replaced by the connected one, and with the crucial difference that the connected propagator is independent of η 0 such that we can use exact scale invariance, 22 to conclude that ψ C n is real which complements the proof we derived in Section 4 using Wick rotations.
Cosmological collider physics As expected, things are a little more involved for the CC scenario since the mode functions are more complicated, however here we prove anti-Hermitian analyticity of the full bulk-bulk propagator, and the connected part, for each helicity mode. To the best our of knowledge this has not been shown before.
(C. 19) The main observation that allows us to make progress is that the differential operatorsD h,n (iη, k) are Hermitian analytic. This follows straightforwardly from the fact that the differential operator in (2.36) is Hermitian analytic (and so acting with it iteratively also yields something Hermitian analytic). The same Hermitian analyticity relations we used above for CCM can then be used to infer that this bulk-bulk propagator is anti-Hermitian analytic: G h n,S (η 1 , η 2 , −k * ) * = −G h n,S (η 1 , η 2 , k) . (C.20) Note that in arriving at this conclusion we also had to make use of the fact thatD h,n (iη, k) are either purely real, for even n − |h|, or purely imaginary, for odd n − |h|, as we discussed in Section 3. This leads toD This is used to cancel the cos(πν) pieces that come from (C.11). We can now use a similar argument as above for the CCM scenario, using the general wavefunction coefficients we discussed in Section 4 for the CC scenario, to conclude that the wavefunction coefficients of massless scalars exchanging such massive spinning fields are purely real, as long as the η 0 dependence cancels out. As with the CCM case, this is only the case for light fields c.f. (3.44). This complements the proof we detailed in Section 4 using Wick rotations.
The situation for heavy fields now follows in the same way as for the CCM scenario: the connected bulkbulk propagator for all modes is Hermitian analytic from which we can easily deduce that the connected part of wavefunction coefficients is purely real. In this CC scenario the connected bulk-bulk propagator is given by (3.48) with reality of this propagator after Wick rotation fixing A h,n (µ) = (−1) n−|h| i cos(πν S ) , (C. 22) which we have again written in a way that is valid for both light and heavy fields. We can then again use the Hermitian analyticity ofD h,n (iη, k), and the fact that these differential operators are purely real for even n − |h| and purely imaginary for odd n − |h|, to conclude that this connected bulk-bulk propagator is anti-Hermitian analytic: C h n,S (η 1 , η 2 , −k * ) * = −C h n,S (η 1 , η 2 , k) .

(C.23)
A caution on locality and the anomaly of Hermitian analyticity In all of the above proofs, locality of the general form of the interactions, meaning that the vertex operator D v is composed of derivatives but not inverse derivatives, is an implicit assumption. In momentum space, locality tells us that D v ∼ (ik) n is a polynomial in the energy variable k = |k| with n ∈ N. Such a local interaction vertex trivially satisfies Hermitian analyticity, i.e. D * v (−k) = D v (k), leading to (C.14) and the reality theorems. However, if the assumption of locality is dropped, and D v (k) is allowed to have a non-polynomial dependence on k, there can be an intriguing "anomaly" of Hermitian analyticity after performing the time integrals, thereby invalidating the conclusions about reality of wavefunction coefficients.
To demonstrate the essential idea, consider the following toy model of a scale-invariant non-local interaction: where ∇ 2 = δ ij ∂ i ∂ j is the three-dimensional Laplacian in flat space. This non-local theory can be understood as describing massless ϕ particles interacting via a Yukawa-like force. It can be derived from a UV theory of ϕ and a massive scalar σ, by integrating out σ and taking the leading order contribution from the time-derivative expansion [34]. The resulting s-channel contact wavefunction is 23 as expected from (C.14). We might then be tempted to use scale invariance and conclude that such a wavefunction coefficient is purely imaginary. However, this is not the case. Indeed, we can carry on and compute the time integral to obtain which, under Hermitian-analytic conjugation, becomes Consequently, the wavefunction coefficient ψ 4 is complex in general (rather than being purely imaginary). Such an anomaly of Hermitian analyticity stems from the non-analytic behaviour of the vertex function with respect to the energy variable s: at any finite time η, there are poles at s = ±i/η which affect the definition of the Hermitian analytic image. One can choose to continue path-wise in the s-plane from either side of the poles, but since the integration time η ranges from the origin all the way to infinity, there is no uniform way to perform the continuation throughout time (see Figure 7). The Hermitian analytic properties of the integrand therefore do not imply some simple relations for the final integrated result when there is some element of non-locality in the interactions. Figure 7: Left panel: At any fixed time η, one can choose to continue from s to −s * by either passing above or below the singularity at i/η. Right panel: After finishing the time integral, the singularities merge into a branch cut that goes all the way from zero to infinity, preventing a uniform definition of the analytic continuation.
Things are somewhat clearer from the Wick rotation perspective which we have primarily used in this paper: the pole at η c = i/s on the complex η-plane prevents the deformation of the integration contour into the Wick-rotated one. Instead, we must include half of the residue at η c to account for the half-circle touring around the pole, which gives exactly (C.31). In fact, this is precisely how parity violation is generated in the non-local single-field EFT model [104].

D Beyond scale invariance: reality in other FLRW spacetimes
The discussion in Appendix C suggests a strong connection between the reality properties we have derived in this work, namely that of wavefunction propagators after Wick rotation, and their Hermitian analyticity properties that have been discussed in the literature. Indeed, with the exact scale invariance of de Sitter space, these two properties are equivalent at the level of equations of motion. To see this equivalence more clearly, recall the equation of motion for a free bosonic field without the chemical potential term: Our reality theorems rely on the fact that after Wick rotation, η = iχ, the equation of motion remains real: Ê (iχ, k) * =Ê(iχ, k) , orÊ * (−iχ, k) =Ê(iχ, k) .