Shapes of gravity: Tensor non-Gaussianity and massive spin-2 fields

If the graviton is the only high spin particle present during inflation, then the form of the observable tensor three-point function is fixed by de Sitter symmetry at leading order in slow-roll, regardless of the theory, to be a linear combination of two possible shapes. This is because there are only a fixed number of possible on-shell cubic structures through which the graviton can self-interact. If additional massive spin-2 degrees of freedom are present, more cubic interaction structures are possible, including those containing interactions between the new fields and the graviton, and self-interactions of the new fields. We study, in a model-independent way, how these interactions can lead to new shapes for the tensor bispectrum. In general, these shapes cannot be computed analytically, but for the case where the only new field is a partially massless spin-2 field we give simple expressions. It is possible for the contribution from additional spin-2 fields to be larger than the intrinsic Einstein gravity bispectrum and provides a mechanism for enhancing the size of the graviton bispectrum relative to the graviton power spectrum.


Introduction
A generic prediction of inflation is a nearly Gaussian spectrum of primordial tensor perturbations, arising from fluctuations of the graviton about the nearly de Sitter inflationary background. In addition to these two-point correlations, tensor non-Gaussianities are generically produced, for instance three-point correlations between three helicity-2 gravitons. Recently there has been some effort to understand the form of these tensor non-Gaussianities [1][2][3][4][5][6][7][8][9] and upcoming experiments, such as LISA [10], will be able to bound the level of tensor non-Gaussianity in the universe [11,12].
If the only field of spin ≥ 2 active during inflation is the graviton itself-and if its leading interactions are described by General Relativity-then the form of tensor non-Gaussianity is fixed by the cubic graviton self-interaction vertex present in the Einstein-Hilbert action [1,2]. Allowing for higher-derivative corrections, there is another possible shape corresponding to the six-derivative, cubic vertex present in a (Weyl) 3 interaction. These are the only possible 3-point shapes for a massless spin-2 field, regardless of other possible higher-derivative terms in the Lagrangian [2], up to slow roll corrections parametrizing departures from pure de Sitter. 1 This fact follows from de Sitter invariance, which restricts the possible on-shell, cubic vertices (cubic terms in the Lagrangian modulo field re-definitions and total derivatives) to be one of two structures, and the fact that the non-Gaussianity is determined solely by the on-shell cubic vertex. This is equivalent to the statement that conformal symmetry on the boundary of de Sitter space fixes the form of the 3-point correlators up to a finite number of constants.
It is possible that other fields besides the inflaton and graviton are present during inflation with masses of order the Hubble scale. If any of these fields have spin ≥ 2, or if the graviton itself has a small mass, then there will be new possible on-shell three-point vertices which are not possible for the massless graviton alone, and correspondingly new possible shapes of the 3-point function if the new degrees of freedom mix with the graviton tensor modes. For example, a single massive spin-2 field has four possible parity-invariant self-interaction structures, as opposed to the 2 possible for a massless spin-2 field (in four spacetime dimensions). These extra structures arise because the interactions for a massive particle are not required to satisfy the gauge invariance constraints that a massless particle must. If these additional massive particles mix with or otherwise transmit their fluctuations to the graviton, then these new structures can imprint themselves on graviton non-Gaussianities and could be detected by futuristic gravitational wave observations. This would be evidence for the presence of new higher-spin degrees of freedom.
In this paper, we perform a model-independent study of the possible new shapes of tensor non-Gaussianity that can arise due to the presence of additional heavy spin-2 fields. We do this by classifying the possible on-shell cubic vertices through which a massive spin-2 field can interact with itself and with the graviton. Any theory of gravity coupled to additional massive spin-2 1 There is a possible parity-violating cubic interaction L ∼ ǫW W W , which would naively generate another shape, but this term only generates slow-roll suppressed non-Gaussianities [4,5]. We do not consider possible parity-violating interactions in this paper. Throughout the introduction we specialize to four spacetime dimensions. particles (for example, those derived from ghost-free massive gravity [13][14][15], bi-gravity [16,17] and multi-gravity [18] theories) must, on-shell at cubic order, reduce to a linear combination of these vertices. The tensor non-Gaussianity depends only on these vertices and can be computed from them, and so the results apply to any model with these degrees of freedom. A powerful feature of the 3-point function in de Sitter space is that its form is completely fixed by the isometries to be a linear combination of a finite number of shapes. This is true non-perturbatively, and therefore it is not even necessary for the additional spin-2 field to be fundamental: even if it is composite or a resonance, its 3-point interactions and correlators will be fixed to by de Sitter symmetry. 2 The approach we take to evaluating late-time 3-point correlation functions is to compute the cosmological wavefunction as the on-shell action evaluated on a classical solution. This approach makes it manifest that correlation functions are sensitive only to on-shell vertices. In order to perform this computation, we derive the bulk-to-boundary propagator for a general mass spin-2 field, along with the generic on-shell vertices between arbitrary admixtures of massless, partially massless, and massive spin-2 fields, which may be of independent interest for AdS/CFT applications.
For generic spin-2 masses, the integrals required to evaluate the on-shell action and compute the late-time non-Gaussianity cannot be evaluated in closed form. In these cases, the best that can be done is a numerical evaluation of the general expression for particular masses of interest. There are, however, specific mass values for which the non-Gaussianity can be evaluated analytically. One of these is the so-called partially massless value where m 2 = 2H 2 . In these cases which only involve massless and partially massless spin-2 fields, we will give explicit expressions for all possible 3-point correlation functions in exact de Sitter space.
The partially massless point is special for several reasons. It saturates the Higuchi bound m 2 ≥ 2H 2 [19], which gives a minimum mass for stable spin-2 fields on de Sitter space (the massless graviton at m 2 = 0 is the only exception). At the partially massless point, a new scalar gauge symmetry emerges, which removes the longitudinal degree of freedom of the massive spin-2, leaving a field with only tensor and vector modes. There are studies and no-go theorems that would seem to forbid consistent theories of a single interacting partially massless spin-2 particle . However, it may still be possible to have theories of partially massless spin-2 particles interacting with other fields, possibly an infinite number of them (for instance there are known examples of Vasiliev-like higher spin theories with infinite towers of partially massless fields [44][45][46][47][48]). Questions about the non-linear completion of the theory will not affect the systematics of our computations in this paper because to compute bi-spectra we only need consistency on-shell up to cubic order, which amounts to finding the vertices that are invariant on-shell under the linearized partially massless symmetry. These vertices are known [49,50], and we reproduce the partially massless ones as a byproduct of our analysis. In situations where the partially massless symmetry fails to exist beyond cubic order there could be additional contributions to correlation functions beyond those that we consider coming from vertices that do not have the partially massless gauge symmetry. Figure 1. Diagrams giving rise to the T 2 and T 3 wavefunction coefficients, which determine γ 2 and γ 3 to leading order. Wavy lines correspond to graviton bulk-to-boundary propagators. It should be noted that other references (such as [51,52]) use similar diagrams to denote the entire in-in correlator, while for us the above diagrams only correspond to wavefunction coefficients. The three-point interaction vertex in the right diagram can arise from the Einstein-Hilbert term or a W 3 µνρσ higher-derivative interaction (other diffeomorphism invariant interactions are redundant with these [2]).
Once we have computed the wavefunction involving mixtures of massless and partially massless fields, we turn to phenomenology. Wavefunction coefficients involving PM fields on external legs can produce observable signatures in multiple ways. For instance, they can be seen directly if the partially massless spin-2 couples directly to matter. Alternatively, if there exists a linear mixing between the PM spin-2 and the graviton, then these coefficients imprint themselves on the graviton bispectrum. Here we explore a mechanism involving a linear mixing between the partially massless field and the graviton. Such a mixing is not possible in exact de Sitter space, and therefore carries a slow-roll suppression. Mixing of this type has two effects: it can transmit non-Gaussianity into the graviton sector at the end of inflation, so that the wavefunctional computed in exact de Sitter space accurately captures the non-Gaussianities in the tensor sector. Another effect is that the PM spin-2 field can mix into the graviton during inflation. We estimate this latter process, but this requires a more-involved computation to treat fully. Aside from 3-point correlation functions, it is of course possible for additional spin-2 fields to affect higher-point graviton correlation functions, and we comment briefly on possible signatures in the four-point function. Our analysis in this case is somewhat preliminary, but there are some intriguing features.
Our results are also of more formal theoretical interest. There has been great progress recently both in systematizing the computation of correlation functions in cosmology using de Sitter symmetry [2,51,[53][54][55][56][57][58][59] and in the related problem of studying conformal field theories in momentum space [60][61][62][63]. To this point, investigations involving external fields with spin have focused on the massless cases, where gauge invariance (or current conservation) provides numerous simplifications. The results provided here provide the first steps toward an understanding of correlation functions in de Sitter space with more general external states. The building blocks we provide may enable many new computations besides the ones we present here and provide additional data to help further develop our systematic understanding of perturbative field theory in cosmological spacetimes.
Many of our results are rather technical, so to orient the reader we first give a brief overview of the computation we perform. Additionally, we provide a rough bound on the size of non-Gaussianity which can be induced by additional spin-2 fields. For concreteness, we focus on the imprint that partially massless fields can leave on the graviton bispectrum γ 3 .

Estimating the size of non-Gaussianity: A sketch
The approach we follow is to compute graviton non-Gaussianities via the wavefunction of the universe, Ψ[φ k , τ ⋆ ], which when squared gives a probability distribution for fields, collectively denoted by ϕ, to take on a given profile,φ k , at time τ = τ ⋆ . We are interested in momentum space correlators, hence our boundary conditions are phrased in momentum space, as indicated. Equal-time expectation values are computed via the usual quantum mechanics formula: Much of our effort will be devoted to computing the wavefunctional itself. This can be done via the path integral which, at leading order, is approximated by the action evaluated on the classical solution, ϕ cl , which takes on the prescribed boundary valuesφ k at τ = τ ⋆ , integrated up to time τ ⋆ This can be expanded in powers ofφ, schematically as for some functions O n . Using the wavefunctional (1.3), the correlators (1.1) can then be computed perturbatively via standard Gaussian integral formulae. A more detailed discussion of the wavefunction formalism may be found in Appendix A.
The wavefunction coefficients in (1.3) have a convenient diagrammatic representation. For instance, in the case of graviton correlators we use the notation ϕ → γ ij and O → T ij . Considering only self-interactions, the T 2 and T 3 coefficients (indices suppressed) would arise from the diagrams 3 in Fig. 1. Here and throughout, is used to represent graviton lines. The left Figure 3. Wavefunction diagrams arising from interactions between the graviton and a partially massless spin-2 field. When combined with a mixing term T Σ , these terms allow the partially massless field to imprint itself on γ 3 . diagram corresponds to evaluating the quadratic action on the linear classical solution, γ cl ∝γ k , while in the right diagram the same solution is inserted into the cubic interaction. The actions are integrated over all of spacetime up to the τ = τ ⋆ surface where correlators are to be computed. Such calculations are familiar to AdS/CFT practitioners: diagrams such as Fig. 1 are the dS version of Witten diagrams and lines correspond to bulk-to-boundary propagators (or bulk-to-bulk propagators in diagrams involving internal lines).
The graviton power spectrum and bispectrum are related to the wavefunction coefficients T 2 and T 3 via relations of the form We now introduce the partially massless spin-2 fields, for which we use the notation ϕ → σ ij , O → Σ ij , and denote them in graphs by . By including self-interactions and cubic couplings between the graviton and the partially massless field, we can generate non-trivial ΣT 2 , T Σ 2 , and Σ 3 coefficients corresponding to the diagrams in Fig. 3. If there exists a mixing term T Σ in the wavefunction, then these cubic coefficients can imprint upon the graviton bispectrum as in Fig. 4. 5 For instance, the middle diagram in Fig. 4 gives a contribution to γ 3 of the schematic form γ 3 ⊃ + + Figure 4. Partially massless spin-2 contributions to γ 3 . Here, a mixing vertex corresponds to a factor of Re T Σ while a PM line, , is a factor of Re Σ 2 −1 .
k ∼ H by dimensional analysis and diagrammatics. The contributions to γ 3 ∼ E −6 , from various interactions, can be estimated as follows: 6 • Graviton self-interactions: The T 2 coefficient depicted on the left in Fig. 1 is generated by the Einstein-Hilbert operator and hence scales as T 2 ∝ M 2 Pl . The cubic coefficient, T 3 , can be generated by the Einstein-Hilbert vertex, or by a higher derivative ∼ M 2 Weyl-cubed vertex [2], and thus has contributions T 3 ∝ M 2 Pl and T 3 ∝ M 2 Pl /Λ 4 W 3 . Therefore, from (1.4) and dimensional analysis, it follows that on scales k ∼ H the bispectrum is of the form 7 so that for Λ W 3 H the higher-derivative shape can be of the same order as the Einstein-Hilbert contribution.
• Interactions with PM spin-2 fields: For concreteness, we estimate the contribution of the middle diagram in Fig. 4 to γ 3 . We can take Σ 2 ∝ M 2 Pl and further write the mixing coefficient as T Σ ∝ Λ 2 mix . From (1.5), this yields The coefficients T Σ and T Σ 2 determine the size of PM-induced non-Gaussianity. Both coefficients also correct the graviton power spectrum, as shown in Fig. 5. In order for the 6 Strictly speaking, the primed wavefunction coefficient γ 3 ′ which has had its momentum conserving delta function removed scales as E 6 , but in this introduction we omit the primes in order to leave the presentation uncluttered. 7 Here, only factors of H are used to fix the dimensions. This relies on the assumption that γ 3 is time-independent at superhorizon scales (as is suggested by the freeze-out behavior of the graviton mode function) which precludes factors of τ from appearing in the analysis. This assumption does not always hold. For instance, a similar analysis of the parity violating ∼ ǫW W W operator would give an estimate of the form (1.6), while the precise calculation demonstrates that this operator only produces dS non-Gaussianity which decays in time [2,4].
Corrections to the graviton power spectrum due to insertions of T Σ (middle) and T Σ 2 (right). The first diagram is represents the familiar result γ 2 ∼ H 2 cubic wavefunction computation only requires knowledge of the action evaluated on-shell, so in Section 3 we enumerate the possible on-shell cubic vertices involving spin-2 fields. In Section 4 we use these ingredients to compute the cubic wavefunction coefficients depicted in Fig. 3, and in Section 5 we use these to compute the contributions to the graviton non-Gaussianities depicted in Fig. 4. We conclude in Section 6. We collect some useful technical results in a number of appendices. In Appendix A we give a brief orientation to the wavefunctional approach in cosmology. In Appendices B and D we list and give some properties of a set of spatial projection tensors which we use in a number of places in the text and give an explicit basis of helicity-2 polarizations, which we use to evaluate correlation functions. In Sec. C we describe a method for deriving gauge invariant on-shell cubic interactions and use this to derive all consistent cubic interactions between massless and partially massless fields in arbitrary dimensions. The wavefunction coefficients are sensitive to both integrations-by-parts and field redefinitions and we discuss these procedures and their relationship in App. E. Finally in Appendix F we review the structure of CFT 2-point functions for spinning fields.

Conventions:
We use mostly plus signature, work in (d + 1) spacetime dimensions, and use the curvature conventions R ρ σµν = ∂ µ Γ ρ νσ + . . . and R µν = R ρ µρν . On de Sitter space, dS d+1 , the Hubble scale is denoted by H, where R = d(d + 1)H 2 . We work exclusively in the flat slicing of de Sitter where τ ∈ (−∞, 0) is the proper time.
Greek letters µ, ν, ρ, . . . are used for spacetime indices, while lower case Latin letters i, j, k, . . . are reserved for spatial indices. Spacetime indices are always raised and lowered with the full metric g µν . Spatial indices will always be manipulated using the flat δ ij metric, e.g., ∂ i h ij = δ ik ∂ i h kj . Tensors are symmetrized and anti-symmetrized with unit weight, e.g., T (µν) = 1 2! (T µν + T νµ ) and T [µν] = 1 2! (T µν − T νµ ). Spatial vectors are bolded as in x ≡ x or k ≡ k. We use the following Fourier conventions: to minimize the number of explicit 2π factors appearing. In three-point correlators, we denote the sum of magnitudes of momenta by k T ≡ k 1 + k 2 + k 3 .

Free spin-fields on dS d+1
We begin by reviewing the physics of free spin-2 fields with general mass on de Sitter space and then derive their propagators and superhorizon two-point functions.

The quadratic action
The degrees of freedom of a massive spin-2 field of mass m on (d + 1)-dimensional de Sitter space are carried by a symmetric 2-index tensor, h µν , with the following quadratic action, where h ≡ h µ µ and all covariant derivatives and contractions are defined with respect to the background dS d+1 metric g µν . We work in planar inflationary coordinates where the background metric takes the form (1.12). The field, h µν , in (2.1) is dimensionless, and the action is normalized such that the m → 0 limit of (2.1) coincides with the expansion of the Einstein-Hilbert action, with a cosmological constant, about dS d+1 . 9 In (2.1), the derivatives have been organized that the entire second line vanishes when h µν is transverse and traceless, ∇ µ h µν = h = 0, as our solutions will always obey these properties.

Distinguished mass values and equations of motion
The linear equations of motion obtained from the action (2.1) are 3) The number of degrees of freedom described by this equation depends on the value of the mass parameter, m 2 . On de Sitter space, there are two distinguished points: • When m 2 = 0, the spin-2 is massless and h µν only propagates tensor degrees of freedom. The massless quadratic action enjoys the usual linearized diffeomorphism gauge symmetry, with a vector gauge parameter, ξ ν , The two conditions ∇ µ h µν = 0 and h = 0 may be imposed as on-shell gauge conditions [66] and the resulting wave equation for the physical degrees of freedom is • When m 2 = (d−1)H 2 , the spin-2 is partially massless [67] and h µν only propagates tensor and vector degrees of freedom. The partially massless quadratic action enjoys a gauge symmetry with a scalar gauge parameter, χ, It is possible to enforce h = 0 as a gauge choice [68], after which the equations of motion further imply that ∇ ν h µν = 0. In this gauge, the resulting equation of motion is Specifically, the Einstein-Hilbert action which admits (1.12) as a solution is given by Perturbing the metric as ds 2 = (gµν + hµν ) dx µ dx ν , with gµν the background dS metric in (1.12), gives SEH = limm→0 S (2) at O(h 2 ) (after integrations by parts).
For other mass values, the spin-2 is massive and h µν propagates tensor, vector and scalar degrees of freedom. At generic values of the mass, by taking traces and divergences of (2.3) one finds the equations ∇ µ h µν = h = 0, which then simplify (2.3) to With the exception of the massless and partially massless points, the theory is only unitary for m 2 > (d − 1)H 2 , the so-called Higuchi bound [19]. Unitary spin-2 fields on dS must therefore belong to one of three categories: In the language of de Sitter representation theory, massless and partially massless fields belong to the exceptional series (which coincides with the discrete series in d = 3), while massive fields can either belong to the complementary or the principal series, depending on the value of their mass: spin-2 fields in the mass range (d−1) < m 2 H 2 ≤ d 2 4 belong to the complementary series, while fields with masses m 2 H 2 > d 2 4 belong to the principal series. All other mass values correspond to non-unitary representations [69][70][71].

Mode functions
In this section we solve the linear equations of motion for generic mass spin-2 fields on dS d+1 in Fourier space. 10 These solutions are typically called mode functions in the cosmology literature. The massless and PM cases can then be obtained as limits of the general solutions. In the following section, these results will be re-packaged into the bulk-to-boundary propagator, which is a solution to the equations of motion with some specified Dirichlet boundary conditions and which plays an important role in the computation of the late-time wavefunctional.

Generic mass solutions
In the generic equations of motion, we trade the mass, m, in favor of a parameter, µ, defined below, after which the equations take on the form The massless and PM cases correspond to iµ = d 2 and iµ = d−2 2 , respectively. As discussed previously, the h = 0 and ∇ µ h µν = 0 conditions appear as equations of motion for generic µ, but correspond to a gauge choice in the massless and PM cases.
We first decompose h µν in ADM-like variables, 10b) 10 A similar calculation in AdS d+1 can be found in [72].
and then further split δN i (τ, x) and ϕ ij (τ, x) into scalar, vector and tensor components which are irreducible with respect to the spatial SO(d) symmetries. This decomposition is most naturally expressed in Fourier space: The various components are transverse and traceless according to 12) and the temporal dependence of all components is being suppressed. The decomposition of ϕ is defined such that the projectors introduced in Appendix B isolate the components shown in (2.11) We now turn to solving the equations of motion. The trace condition, h = 0, fixes the lapse, δN , in terms of S: (2.14) Writing ∇ µ h µν = E ν , the divergence constraint contains the two scalar conditions: E k 0 = k i E k i = 0. Written in terms of the variables in (2.11), these become These equations can be solved for δN k L and Q k : The divergence condition also contains the following vector constraint, 11

18)
11 This is isolated by applying the transverse projector π k ij to E k i , where The constraint (2.18) shows that V T, k i is determined by δN T, k i . It proves fruitful at this point to express the remaining δN T, k i , ϕ T T, k ij , V T, k i , and S k components as time-dependent functions multiplying time-independent tensor structures. Since the vector constraint implies that V T, k i ∝ δN T, k i , these two components also share the same tensor structure, and hence we can write whereφ,V andS are τ -independent. After these replacements, (2.18) reads (2.20) Using the above constraints and definitions, the remaining wave equation can be solved straightforwardly. Only some of its components give non-trivial relations: • The E k 00 component gives an equation for f S alone: • The π k ij E k 0j component gives an equation for f N T alone: • Finally, the transverse, traceless components (Π k T T ) ij lm E k lm give an equation for f T T alone: The general solutions to (2.22), (2.23), and (2.24) are expressions involving Hankel functions of first and second kinds. However, the Bunch-Davies vacuum condition requires solutions which behave as ∼ e +ikτ at early times, which corresponds to only using the Hankel function of the second kind. The desired solutions to the linear equations are then given by iµ (−kτ ), (2.25) where f V was determined through (2.20). The constrained fields are then determined using (2.16) and (2.14). Making the definitions we find the following expressions for the mode functions: Putting all this together, the generic momentum-space solution for h µν can be written as where the (τ, k)-dependence of the f i 's has been suppressed. The generic, massive solution is determined by the boundary data contained in the (d + 1)(d − 2)/2 components of the transverse, traceless tensor,φ T T, k ij , the (d−1) components of the transverse vector,V k i , and the one component ofS k , totaling (d(d + 1) − 2)/2 degrees of freedom, which is the expected counting for a massive spin-2 [14].

Partially massless solutions
We now discuss simplifications which can be made in the partially massless limit, where iµ = d−2 2 . After choosing the gauge such that h = ∇ µ h µν = 0, there remain residual gauge-transformations of the form δh µν = ∇ µ ∇ ν + H 2 g µν χ , We can solve this equation in Fourier space in the same way as in the previous Section. We find that the solution for a residual χ is where all the mode functions should be understood to be evaluated at iµ = d−2 2 . The PM solution is determined by the (d + 1)(d − 2)/2 components of the transverse, tracelessφ T T, k ij and the (d − 1) components of the transverseV T, k i , totaling (d(d + 1) − 4)/2 degrees of freedom.

Massless solutions
Even more drastic simplifications can be made in the massless limit, where iµ = d 2 . After choosing the gauge such that h = ∇ µ h µν = 0, there remain residual gauge-transformations of the form Using the same techniques, we can find the solution for a residual gauge transformation as a function of boundary data: (2.36b) The massless solution is determined by the (d + 1)(d − 2)/2 components of the transverse, traceless ϕ T T, k ij .

Bulk-to-boundary propagators
The perturbative approach to computing the wavefunction we are taking amounts to computing the on-shell action as a function of boundary data on some τ = τ ⋆ surface. The solutions to the linear equations of motion which satisfy some specified set of Dirichlet boundary conditions are known as bulk-to-boundary propagators. It is straightforward to re-interpret the solutions from the previous Section as these objects. In particular, we will write this relation in the form whereφ k ij is the boundary value of the ϕ ij field appearing in (2.10) and K k µν ij (τ ) is the bulk-toboundary propagator. Since we requireφ ij to be the boundary value of the field h µν , we therefore require the spatial parts of the bulk-to-boundary propagator to satisfy a relation of the schematic form where 1 is a type of identity matrix whose precise form will depend on the case at hand.
We now build the massless and PM propagators, which are relatively simple, and end with the massive case, which is slightly more cumbersome. In doing this, it is extremely useful to use the spatial projection tensors defined in Appendix B to organize the bulk-to-boundary propagator into its irreducible components.

The massless bulk-to-boundary propagator
We first construct the massless bulk-to-boundary propagator. From (2.36), the massless propagator is purely spatial, K k 0ν ij = 0. Its components are given by At the boundary, K k ij lm (τ ) reduces to the identity matrix for transverse, traceless tensors by construction, i.e., K(τ ) τ →τ⋆ −−−→ Π T T . In the massless case, we define σ ij → γ ij in order to conform to standard notation for massless spin-2's in the literature and to distinguish the graviton from the PM and massive spin-2 fields. Hence, the bulk solution will be written as rather than as (2.37), andγ k ij is required to be transverse and traceless. 12

The partially massless bulk-to-boundary propagator
Next we consider the bulk-to-boundary propagator of a partially massless field. From (2.33), the PM propagator has both spatial and temporal components. These are given by 12 A sidenote: since K is transverse and traceless, so too will the wavefunctional coefficients be. However, we know that these coefficients should have non-transverse-traceless pieces dictated by the stress tensor Ward identities. In fact, these additional pieces can be reconstructed using the Ward identities [60]. Alternatively, these pieces can be computed directly in a different gauge where the bulk-to-boundary propagator is not transverse-traceless.
At the boundary, K k ij lm (τ ) reduces to the identity matrix for symmetric, two-index tensors which only carry tensor and vector components by construction, i.e., K(τ ) In the PM case, we will redefine ϕ ij → σ ij in order to distinguish the PM from the massless and massive spin-2 fields. As previously stated, for PM fields, the boundary dataσ k ij is required to be free of scalar components, because they can always be gauged to zero: Π k S ·σ k = Π k Q ·σ k = 0, in condensed notation. Hence, the bulk solution will be written as

The massive bulk-to-boundary propagator
We now turn to the fully massive bulk-to-boundary propagator. From (2.30), we see that the generic mass propagator has both spatial and temporal components. In this case, all of the components of the boundary dataφ ij are generically non-zero, but they cannot all be independently chosenthey have to be consistent with the constraint structure of the theory-as was seen in Sec. 2.2.1. Specifically, the Π k Q ·φ k component of the boundary data is constrained in terms of the Π k S · ϕ k component. This makes the construction of the bulk-to-boundary propagator slightly more complicated than the massless and PM cases. Due to this fact, we will defineφ k ij , which contains only the unconstrained components ofφ k ij , i.e.,φ k ij obeys Π Q ·φ = 0 but coincides withφ ij otherwise, 13 and the bulk solution will be expressed as The components of the bulk-to-boundary propagator are then given by 13 Explicitly, using the projectors in App. B we can writē where f Q (τ, k) is the lengthy expression in (2.29). We have chosen K k ij lm (τ ) such that its scalar, vector, and tensor components reduce to their respective identity matrices when τ = τ ⋆ , i.e.,

Quadratic on-shell action and two-point functions
In this section we evaluate the quadratic on-shell action and derive the superhorizon power spectrum for the transverse, traceless components of spin-2 fields in general dimensions.
Conforming to the standard variables used in cosmology, we will compute the correlators and wavefunctions corresponding to the field ϕ k and the notation ϕ → γ and ϕ → σ will be used for the massless and PM cases, respectively. Because we are only focusing on the transverse, traceless components of correlators, our on-shell solutions for ϕ k ij (τ ) are related to the boundary dataφ k ij (τ ) via For later convenience, we explicitly write the form of the transverse, traceless d = 3 solutions for the massless and partially massless fields:

The on-shell action and generic two-point functions
We can obtain the transverse, traceless two-point function for generic spin-2 modes by evaluating the quadratic action on the solution (2.48).
Starting from the quadratic action (2.1), we integrate by parts to obtain Hereḡ ij is the induced metric on the τ = τ ⋆ boundary, n µ is the outward-pointing normal vector to the boundary, and E µν ρσ is the differential operator appearing in the free equation of motion: E ρσ µν h ρσ = 0. When evaluated on-shell, the second term vanishes and we are left with When we further use the ∇ µ h µν = h = 0 restriction, all terms but the first vanish and the on-shell action reduces to Finally, if we convert from h µν to ϕ ij via (2.47), the on-shell action for the transverse, traceless modes becomes The action is obtained as a function of boundary data by using (2.48) for ϕ ij : The wavefunctional is written in terms of the on-shell action as Ψ ∼ e iS and we parameterize its quadratic component as where the Gaussian wavefunctional coefficient is given by (2.55) and where Π k T T is the projector onto symmetric, transverse, traceless tensors (B.2). Only the real part of the wavefunctional coefficient (respectively the imaginary component of the action) affects correlation functions of ϕ ij , as they are governed by the probability distribution |Ψ| 2 . The real part of the coefficient is 14 When iµ ∈ R + , as in the cases of a massless or partially massless spin-2 field, the superhorizon kτ ⋆ → 0 limit of (2.56) is given by We now specialize to the cases of a massless spin-2 field and a partially massless spin-2 field.

The massless spin-2 two-point function
In the massless limit, we write the quadratic wavefunction as 15 where the second expression follows from evaluating (2.57) with iµ = d 2 . By squaring and integrating the wavefunctional overγ, and using (2.58) and (A.14), we find the superhorizon two-point function is given by We are primarily interested in the d → 3 limit, where we recover the well-known result (2.60) 14 Here, we have used H iµ (−kτ ) which holds for both real and imaginary µ [51]. 15 We write the "dual operator" to the massless γ k ij field as T ij k and that of the partially massless σ k ij field as Σ ij k , using the language of holography.

The partially massless spin-2 two-point function
In the partially massless limit, we write the quadratic wavefunction as where the second expression follows from evaluating (2.57) with iµ = d−2 2 . We have only written the part of the wavefunctional corresponding to tensor polarizations. This result is sufficient for the interests of this paper, but in general a partially massless spin-2 solution will carry both tensor and vector modes which propagate. 16 Using (2.58) and (A.14), we can compute the transverse, traceless parts of the superhorizon two-point function for a PM field We are primarily interested in the d → 3 limit where we obtain:

On-shell cubic vertices
We now turn to interactions. Our focus is on three-point correlation functions involving mixtures of spin-2 fields. There is a certain universality in these objects; in the de Sitter limit, their form is completely fixed by conformal invariance up to a finite number of constants [73], much as the structure of three-point on-shell scattering amplitudes is fixed by Lorentz invariance [74]. In order to compute three-point correlators, only the cubic component of the on-shell action is required. For instance, the GR contribution to the massless tensor bispectrum γ 3 comes from the T 3 wavefunction coefficient, which is related to the cubic on-shell action via The universality of three-point correlation functions manifests itself in the fact that there are only a finite number of independent on-shell cubic interactions for a given set of fields of various 16 That is, this calculation has only included the ΠT T parts of (2.41). Had we also included the ΠV vector mode components, then Σ 2 would also contain a piece ∝ ΠV . However, because we are ultimately interested in the induced non-Gaussianity for γij , all such vector components are projected out since the final answer can only be proportional to ΠT T tensors. For this reason, vector modes are ignored throughout this paper. If desired, the ΠV components of the PM two-point function can be restored using the results of Appendix F, as they are tied to the ΠT T terms by conformal symmetry.
spins. Therefore, we can compute the most general, model-independent, three-point correlation function by finding a basis for these cubic structures. In this Section, we describe such a basis-the detailed construction of these interactions is discussed in Appendix C.

Generic spin-2 on-shell interactions
We begin by describing a basis for generic mass cubic spin-2 operators in dS d+1 . Throughout, the vertices are on-shell, meaning we are imposing the conditions ∇ µ h µν = h = 0 and h µν ∝ h µν when identifying independent operators.
• Basis of three different spin-2 fields: We begin with the most general parity preserving cubic interactions amongst three distinguishable spin-2 fields h (i)µν , i ∈ {1, 2, 3}, which are each taken to have unique and generic masses. There are 11 parity preserving CFT structures involving three spin-2 primaries with generic weights [73,75]. Hence, we expect there to be 11 independent on-shell cubic interactions. A basis of eleven operators which are independent is given by: Any cubic term in the action can be written as a linear combination of these vertices after integrations by parts and using the conditions on h µν stated at the start of this section.
• Basis when two fields are identical: When two out of the three spin-2 fields are identical, the action is symmetric under their interchange and some of the operators in (3.2) degenerate.
Only eight independent interactions remain in this case. Again, this is the expected counting from CFT or S-matrix considerations [73,75]. Our basis for one h (1) interacting with two h (2) fields is • Basis for a single spin-2 field: Finally, we consider the self-interactions of a single spin-2, h (i) µν = h µν . In this case, the interactions have to be totally symmetric, which further reduces the number of independent operators to five. Again, this is the expected counting from CFT or S-matrix considerations [73,75,76], and the basis we use for self-interactions is The operator bases (3.2), (3.3), and (3.4) are written for generic mass spin-2 fields. When one or more of the h (i)µν fields are massless or partially massless, gauge invariance places further constraints on the operator combinations. This reduces the number of independent cubic interaction terms. In Appendix C we describe how to implement on-shell gauge invariance at the level of the action. In what follows, we merely catalog the results of this procedure.

Massless self-interactions
We first consider the self-interactions of a single massless field. Before imposing gauge invariance, the possible interactions are of the form (3.4). In general dimension there exists a three-parameter family of on-shell gauge-invariant combinations of these interactions, where a 1 and a 2 are fixed to be and a 3 , a 4 , a 5 are free. It can be verified that the set of interactions described by (3.5) is on-shell equivalent to a sum of the Einstein-Hilbert term (2.2), the Weyl tensor squared, and the Weyl tensor cubed, up to integrations-by-parts. Furthermore, this counting matches the well-known counting of massless spin-2 scattering amplitudes or of stress tensor 3-point functions in general dimensions.
In d = 3, it is instead fruitful to map (3.5) to a sum of the Einstein-Hilbert term (2.2), the Gauss-Bonnet combination, and the Weyl tensor cubed. In particular, the d = 3 Gauss-Bonnet with a 2 = a 4 = a 5 = 0, which can be shown by direct computation. The Gauss-Bonnet combination is a total derivative in d = 3.

Partially massless self-interactions
In this section, we present the cubic self-interactions of partially massless fields in arbitrary dimensions. Again the fields are identical, so the initial basis of operators is of the form (3.4) and we find two branches of solutions for the a i : one branch which is gauge invariant for arbitrary d and one which only exists in d = 3.
Arbitrary d Solution: On one branch of solutions, d can be kept arbitrary. In this case, we find a two-parameter family of solutions, three of the coefficients are fixed to be where a 4 , a 5 are arbitrary. d = 3 Solution: On the other branch of solutions we are forced to set d = 3, in which case we find a three-parameter family of solutions where we fix with free coefficients, a 3 , a 4 , a 5 .
Our d-dimensional result (3.7) matches the counting found in [49]: there is one operator whose highest derivative components are O(∇ 6 ) and one whose are O(∇ 4 ). In the d → 3 limit of (3.7), the four-derivative a 4 interaction degenerates to the Gauss-Bonnet total derivative which is invariant under the partially massless gauge symmetry only in d = 3: On the d = 3 branch of solutions (3.8), we have one additional operator whose highest derivative components are only O(∇ 2 ). One combination of the a i 's reproduces the six-derivative interaction (one six-derivative combination is on-shell equivalent to the cubic interactions coming from a W 3 µνρσ term) of (3.7), one combination of the a 4 and a 3 terms corresponds to L GB , and another combination of a 4 and a 3 reproduces the cubic self-interaction found in [26].

PM-massless-massless interactions
Next we consider cubic interactions between two massless fields and one partially massless spin-2.
The basis of operators is of the form (3.3), where h (1)µν is partially massless and the two h (2)µν 's are massless. We find two branches of solutions for the b i 's: one which holds for arbitrary d and another which requires d = 3.
Arbitrary d solution: On the first branch of solutions, d can be kept arbitrary. In this case, we find a one-parameter family of solutions for the b i , where On the other branch of solutions we are forced to set d = 3, in which case we find a two-parameter family of solutions for the b i : where the free coefficients are b 4b and b 5 .
Our d-dimensional result (3.10) matches the counting found in [49]: there is one operator whose highest derivative components are O(∇ 6 ). On the d = 3 branch of solutions (3.11), we have one additional operator whose highest derivative components are only O(∇ 4 ). One combination of b i in (3.11) reproduces the six-derivative interaction of (3.10).

Massless-PM-PM interactions
In this section we consider the cubic interactions between two partially massless fields and one massless spin-2. The basis of operators is of the form (3.3), where h (1)µν is massless and the two h (2)µν 's are partially massless. Imposing gauge invariance, we find a single branch of solutions which holds in arbitrary dimensions. The family of solutions depend on five free parameters, the coefficients are fixed, while the 5 parameters b 3a , b 3b , b 4a , b 4b and b 5 are free. The above results match the counting of [49]: there is one operator whose highest derivative components are O(∇ 6 ), two whose are O(∇ 4 ), and two whose are O(∇ 2 ).

Massive-massless-massless interactions
In this section we consider the cubic interactions between two massless fields and one massive spin-2. The basis of operators is of the form (3.3), where h (1)µν is massive and the two h (2)µν 's are massless. Imposing gauge invariance, we find a single branch of solutions which holds in arbitrary dimensions. The family of solutions depend on two free parameters, the coefficients are fixed, while the 2 parameters b 4b and b 5 are free. This matches the counting of 3-point scattering amplitudes found, for example, in [77], and the counting of CFT correlators [73,75].

Massless-massive-massive interactions
In this section we consider the cubic interactions between two massive fields and one massless spin-2. The basis of operators is of the form (3.3), where h (1)µν is massless and the two h (2)µν 's are massive. Imposing gauge invariance, we find a single branch of solutions which holds in arbitrary dimensions. The family of solutions depend on six free parameters, the coefficients are fixed, while the 6 parameters b 2a , b 3a , b 3b , b 4a , b 4b and b 5 are free. This again matches the counting of scattering amplitudes [77] and the counting of CFT correlators [73,75].

Massive-PM-PM interactions
In this section we consider the cubic interactions between two partially massless fields and one massive spin-2. The basis of operators is of the form (3.3), where h (1)µν is massive and the two h (2)µν 's are partially massless. Imposing gauge invariance, we find a single branch of solutions which holds in arbitrary dimensions. The family of solutions depend on five free parameters, the coefficients are fixed, while the 5 parameters b 3a , b 3b , b 4a , b 4b and b 5 are free. (Here [49] find instead a 6 parameter family.)

PM-massive-massive interactions
In this section we consider the cubic interactions between two massive fields and one partially massless spin-2. The basis of operators is of the form (3.3), where h (1)µν is partially massless and the two h (2)µν 's are massive. Imposing gauge invariance, we find a single branch of solutions which holds in arbitrary dimensions. The family of solutions depend on five free parameters, the coefficients are fixed, while the 5 parameters b 3a , b 3b , b 4a , b 4b and b 5 are free. (The case of two different masses is studied in [49] and they find instead an 8 parameter family.)

Massless-PM-massive interactions
In this section we consider the cubic interactions between two massless fields and one massive spin-2. The basis of operators is of the form (3.2), where h (1)µν is massless, h (2)µν is partially massless, and h (1)µν is massive. Imposing gauge invariance, we find a single branch of solutions which holds in arbitrary dimensions. The family of solutions depend on two free parameters, the coefficients

Cubic wavefunction coefficients
We now have all of the elements required to compute three-point correlation functions involving arbitrary admixtures of massless and massive fields. Through the use of the bulk-to-boundary propagators listed in Section 2.3 and the on-shell cubic vertices enumerated in Section 3, we can, in principle, compute three-point correlation functions in full generality. However, we will we make a number of simplifications in our concrete calculations: • We only consider cases involving massless and partially massless fields. 17 • We work in d = 3.
• We only consider helicity-2 polarizations for the fields (equivalently, we restrict to transverse, traceless boundary data for all fields).
The first two restrictions are chosen for technical convenience in order to obtain closed form solutions. In general dimensions, or for general graviton masses, the τ integrals involved in computing the wavefunctional do not evaluate to closed-form expressions, but there is no in-principle obstruction to evaluating them numerically for any given case. Additionally, partially massless spin-2 fields are the lightest unitary massive spin-2 representations; we therefore may expect them to be of the most phenomenological interest. The final restriction arises as our primary interest is to understand how the presence of additional fields can modify the correlation functions of the massless graviton, this can be relaxed straightforwardly. Since the massless graviton only carries helicity-2 modes, only the helicity-2 modes of other particles can affect the graviton bispectrum.

Presentation of results
Even in the simplified cases that we consider, the resulting correlation functions are rather complex, so it is worthwhile to briefly describe how the cubic wavefunction coefficients will be presented.
On de Sitter space, spin-2 degrees of freedom can always be diagonalized so that the late-time wavefunctional takes the following form: (4.1) In the above, we have specialized to the case with one massless and one PM field, withγ ij andσ ij the transverse, traceless boundary values of these respective fields on the τ = τ ⋆ time slice, but similar expressions would hold in the presence of additional spin-2 particles. The variables T and Σ transform under the de Sitter isometries as conformal primaries with weights ∆ T = d, ∆ Σ = d − 1.
In (4.1), theγ ij andσ ij tensors play a role analogous to that of polarization tensors in scattering amplitudes. In d = 3, it can sometimes be inconvenient to keepγ ij andσ ij arbitrary because of the presence of Gram/Schouten identities-examples of which are discussed in Appendix D. It is therefore convenient to choose an explicit basis of polarizations and compute the wavefunctional coefficients in this basis. There is no loss of information in doing this. In principle the wavefunctional coefficient for arbitrary boundary data can be reconstructed as a linear combination of the values in an explicit basis.
We use the same basis of polarizations considered in [2] and choose theγ k ij 's to be one of the explicit polarization tensors ǫ kP ij and ǫ kX ij , which are defined in Appendix D. The polarization ǫ kP ij is parity even and ǫ kX ij is parity odd. Since the interactions we consider preserve parity, only contractions with three ǫ P 's or one ǫ P and two ǫ X are non-vanishing. Hence we define, for example and similar for the Σ 3 , ΣT 2 , and T Σ 2 terms. It is these combinations which we calculate. Also, we will report only the real part of these wavefunctional coefficients, as these are the only parts which contribute to |Ψ| 2 and hence to correlation functions for the graviton. This is merely to make the final expressions manageable.
Another subtlety to address is that our bases for interactions (3.2), (3.3), and (3.4) are ambiguous up to integrations-by-parts. Other, equally valid, operator bases hence will differ by boundary terms and we must therefore keep track of possible contributions to correlation functions coming from such operators. Remarkably, all possible boundary terms seem to produce the same shape in correlation functions when evaluated in an explicit basis of polarizations-at least in the case where all fields are the same, or two fields are the same and one is different. 18 We may therefore take as boundary term representatives the interaction in the case where two different fields interact or in the case of self-interactions, where λ is an arbitrary coefficient. As discussed in Appendix E, the resulting shapes are also those associated to local redefinitions of the fields.
When presenting results for wavefunction coefficients, we therefore give both the result which follows from our choice of basis as well as the shape produced by (4.4) or (4.5), which represent the ambiguous parts of the coefficients. It turns out that in our special case of interest where d = 3, the boundary term shapes are redundant with those arising from our bulk interactions due to the existence of dimension-dependent Gauss-Bonnet total derivatives. However, we still give the form of the boundary shapes in order to demonstrate this fact. 19 When non-Gaussianities are computed in Sec. 5, we will similarly compute the shapes that follow from our choice of basis as well as the shapes produced by boundary terms.
One way to fix the boundary term ambiguity in general is to demand that the wavefunction transform correctly under gauge transformations on the τ = τ ⋆ surface [78]. Indeed, one can think of Maldacena's consistency relation [1] as fixing the contact term ambiguity in the cosmological wavefunction, as f local NL corresponds precisely to a contact term. For pure GR, it was demonstrated in [79] that demanding the gauge-invariance of the wavefunctional fully fixes the ambiguous terms, which in turn fixes the squeezed limit of γ 3 . However, our cubic on-shell analysis only requires the lowest-order gauge transformations of γ ij and σ ij , and performing a similar analysis in our cases would require knowing how the gauge transformations are deformed at next order-such an off-shell analysis is beyond the scope of the present paper. 20

The T 3 coefficient
We begin by computing the cubic wavefunctional coefficient for three massless gravitons, T ij which arises from the interactions (3.4) with coefficients (3.5).

Bulk interactions
Here we list the contribution to the wavefunctional coefficients proportional to each of the free parameters-in this case there are 3 free parameters, each of which multiply a particular linear combination of the interactions (3.4).
• The shape proportional to a 3 is: (4.6) • The shape proportional to a 4 is (4.7) • The shape proportional to a 5 is: For instance, due to the interactions with the PM spin-2 particle, the graviton's linear transformation law is not expected to simply be of the usual schematic γ → γ + ∂ǫ + ǫ∂γ diffeomorphism form, but will also include O(σ) terms. The analysis of the present paper is not sensitive to these corrections.
All other non-trivial correlators follow from permutations of the above results. In Sec. 5.1 we will verify explicitly that the above results reproduce the results of [2]. Note that each of these coefficients is singular in the k T → 0 limit, where the total energy of the bulk interaction is conserved. The residue of this singularity is precisely the flat space scattering amplitude. We can think of this as a signature of the fact that these correlation functions came from local interactions in the bulk [2,59,80,81].
As expected, each of these coefficients is time-independent, reflecting the fact that graviton perturbations freeze out at long wavelengths. Something worth noting is that if we compute the wavefunctional for the wrong linear combinations of operators, i.e., if combinations other than those in (3.5) are used, then the wavefunction will have ∼ log kτ ⋆ factors. However, these all cancel once the appropriate gauge-invariant combinations are used. Similar results hold for the wavefunction coefficients calculated in later sections: logarithms appear when non-gauge invariant values for the coefficients are used, but cancel once gauge-invariance is imposed on the interactions.

Boundary terms
Using the boundary term (4.5), converting from h µν → γ ij , and evaluating the above result on the solution (2.49), the result is where we have only kept the leading imaginary terms in the superhorizon limit. This result is in agreement with the expectations of Appendix E, as (4.9) is the shape that arises from a local field-redefinition of the form γ ij (x) → γ ij (x) + c γ il (x)γ l j (x) for constant c. Evaluating (4.9) on the explicit polarization tensors gives: (4.10) Note that these terms are completely regular in the k T → 0 limit, which is consistent with them having an intrinsically boundary origin.

Comments on the Gauss-Bonnet term
In d = 3, the Gauss-Bonnet combination is a total derivative. This manifests itself in terms of the wavefunctional coefficients as the fact that the following sum vanishes when evaluated on the following parameter values This can be explicitly seen by substituting the above values of the a i s into (3.5), which gives a result proportional to (3.6). Similar Gauss-Bonnet-like combinations exist for the other cubic interactions of interest; 21 hence in the following we will find additional combinations of correlators that are degenerate. These identities make our results somewhat insensitive to the ambiguities regarding choosing a basis of interactions and integrations by parts, since the effect of boundary terms can be traded for shifts in the values of the bulk interaction coefficients.
Note that the vanishing of (4.11) requires a contribution from a boundary term. A direct calculation of the cubic wavefunction coefficient induced by the Gauss-Bonnet term, L GB , proceeds similarly. If we compute the cubic coefficient for the variable γ ij introduced as h ij = 1 (Hτ ) 2 γ ij without performing any integrations-by-parts in L GB , the result is non-zero and of the form (4.10). This non-zero result can be understood from the fact when a manifold has a boundary, such as the τ = τ ⋆ surface, integrating L GB over the manifold only produces a topological invariant if an appropriate boundary term is added to the action [82]. If the boundary term is added to the action, then the bulk result is completely cancelled and the wavefunction coefficient induced by the Gauss-Bonnet term vanishes, as it should. 22 21 For instance, the Gauss-Bonnet-like cubic interaction for an interaction involving one h (1)µν field with two h (2)µν fields can be derived from taking the O(h 3 µν ) terms in LGB, replacing hµν → h (1)µν +h (2)µν everywhere and extracting the O(h (1)µν h 2 (2)µν ) terms. 22 Alternatively, if we instead introduce γij via hij = 1 (Hτ ) 2 exp[γ]ij , then the bulk and boundary contributions to the wavefunction coefficient both separately vanish because the field redefinition introduces the correct boundary term to cancel the bulk contribution.

The Σ 3 coefficient
In this Section, we compute the cubic wavefunctional coefficients Σ ij ′ for a self-interacting partially massless field in d = 3, corresponding to the interactions (3.4) with parameter values (3.8).

Bulk interactions
We find the following contributions to Σ 3 , there is a three-parameter family of shapes • The shape proportional to a 3 is: (4.13) • The shape proportional to a 4 is: (4.14) • The shape proportional to a 5 is: All other non-trivial correlators can be obtained from permutations of these.
Note that each of these coefficients scales as ∼ τ −2 ⋆ . This is not the scaling that one would naively expect from the fact that the partially massless mode functions decay as ∼ τ ⋆ on superhorizon scales. Instead we would expect the partially massless three-point correlator to scale as σ 3 ∼ τ 3 ⋆ . We can translate this into a scaling for the wavefunctional coefficient using and we can see from (2.63) that Re Σ 2 ∝ τ −2 ⋆ . Therefore, the natural expectation is that Re Σ 3 ∝ τ −3 ⋆ , while the explicit answer scales as Re Σ 3 ∝ τ −2 ⋆ , resulting in a faster-than-expected decay of correlation functions. 23 This also indicates that we should not expect self-interactions of partially massless fields to imprint themselves on the cosmologically interesting part of the graviton bispectrum, as will be seen in Sec. 5.2.2.

Boundary terms
Using the boundary term (4.5), converting from h µν → σ ij , and evaluating the above result on the solution (2.49), the result is where we have only kept the leading imaginary terms in the superhorizon limit. The result (4.17) again corresponds to the shape arising from a local field-redefinition of the form σ ij (x) → σ ij (x) + c σ il (x)σ l j (x). Evaluating (4.17) for explicit polarization tensors gives: It can be checked that there is another Gauss-Bonnet relation of the type (4.11) when the results (4.18) are added to the above coefficients (4.13), (4.14), and (4.15).

The ΣT 2 coefficient
Next we consider the interaction between two massless spin-2 fields and a partially massless spin-2 field, Σ ij ′ . This corresponds to substituting the parameter values (3.11) into the interactions (3.3).

Bulk interactions
In this case, the entire leading late-time answer is subleading terms in the calculation, it is found that all O(τ −1 ⋆ ) pieces vanish entirely. 24 This implies that the ΣT 2 coefficient does not source the time-independent component of the superhorizon graviton bispectrum; see Sec. 5.2.3.

Boundary terms
Using the boundary term (4.4), converting from h (1)µν → σ ij and h (2)µν → γ ij , and evaluating the above result on the solutions (2.49) and (2.49), the result is where we have only kept the leading imaginary terms in the superhorizon limit. This is the same form as (4.19), as discussed. The result (4.20) corresponds to the shape arising from a local field-redefinition of the form 25 . Evaluating (4.20) on the explicit polarization tensors gives: (4.21)

The T Σ 2 coefficient
Finally, we consider the situation with two partially massless fields interacting with a single massless field, leading to the cubic wavefunctional coefficient T ij

Bulk interactions
The parameter choices (3.12) leave 5 free parameters, corresponding to the following shapes • The shape proportional to b 3a is: The O(τ −1 ⋆ ) terms in Im ΣT 2 are not vanishing, even when evaluated on (3.11), which are again terms which would be interesting in the AdS/CFT analysis of the same bulk fields. 25 In (4.20) there are also subleading ∼ k 3 σγ 2 contributions which correspond to redefining γ → γ + cσγ.
• The shape proportional to b 3b is: (4.23) • The shape proportional to b 4a is: (4.24) • The shape proportional to b 4b is: • The shape proportional to b 5 is: All other non-trivial correlators follow from permutations of the above results.

Boundary terms
Using the boundary term (4.4), converting from h (1)µν → γ ij and h (2)µν → σ ij , and evaluating the above result on the solutions (2.49) and (2.49), the result is where we have only kept the leading imaginary terms in the superhorizon limit. The result (4.27) corresponds to the shape arising from a local field-redefinition of the form 26 σ ij (x) → σ ij (x) + c σ il (x)γ l j (x). Evaluating (4.27) on the explicit polarization tensors gives: Once again, it can be verified that there is a Gauss-Bonnet relation of the type (4.11) when the results (4.21) are added to the sum of T Σ 2 coefficients above.

Graviton non-Gaussianities
The wavefunctional coefficients computed in Section 4 involve partially massless spin-2 fields on the external lines, and therefore will lead to mixed massless-PM 3-point correlation functions. Though these may be of independent phenomenological interest if the partially massless spin-2 itself couples to baryonic matter, in this section we wish to consider the effects of these wavefunctional coefficients on pure graviton correlation functions. This requires a linear mixing that can convert the PM spin-2 external lines into massless graviton legs. Such a mixing is not possible in pure de Sitter space, so we expect that such mixings will be slow-roll suppressed. Nevertheless, this suppression can be compensated by the size of the 3-point couplings, as was demonstrated in Sec. 1.1.
In this section we discuss two different effects of mixing. In the first case, we imagine that the conversion to massless gravitons happens at the end of inflation, so that the late-time wavefunctional computed in exact de Sitter space gives the leading contribution, and that the mixing occurs due to the off-diagonal 2-point function between the PM field and the graviton in the wavefunction. Here we focus on the 3-point interactions involving 2 PM fields and one massless graviton, since it happens that these are the only vertices that lead to a constant late-time graviton 3-point function. We also consider the effect of linear mixings directly in quasi-de Sitter space, which allows the conversion to take place during inflation, which can also affect graviton correlation functions at late times, particularly in soft limits. In technical terms, the former processes depend on the cubic coefficients calculated in the previous section, while the latter stem from contributions of partially massless fields to T 3 via bulk-mixing diagrams which have not yet been discussed.
We give the contributions of the computed coefficients to the bispectrum γ 3 , and comment on their contributions to the trispectrum γ 4 .

Pure graviton bispectrum
We begin by computing the graviton bispectrum induced by the T 3 wavefunction coefficient. These are the contributions to γ 3 from General Relativity and higher order curvature terms, which were originally computed in [2]. We verify here that our calculations reproduce those results.
The relation between γ 3 and the wavefunction coefficients follows from a straightforward generalization of (A.14). Completing the computation, we find that the shapes due to the various components are: • The shape proportional to a 3 is: • The shape proportional to a 4 is • The shape proportional to a 5 is: We also calculate the contribution to γ 3 from the purely boundary terms (4.9): Comparison to [2]: The results in (5.1) reproduce the Einstein-Hilbert shape found in [2] if we set a 3 = M 2 Pl /4 and a 4 = a 5 = 0. The W 3 µνρσ /Λ 2 W shape of [2] is reproduced by setting . Due to the relation (4.11), the values of the a i 's can be shifted around in these matchings at the cost of also adding boundary term shapes to the correlator.

Graviton bispectrum from interactions involving partially massless fields
All remaining cubic wavefunction coefficients involve at least one factor of Σ, which encodes the effects of partially massless fields. In order for the these components to affect the observable graviton bispectrum γ 3 , we require a way to convert the PM spin-2 into the massless graviton: this requires the presence of a mixing term ∼ T Σ γσ in the wavefunction. In the following sections, we discuss the generation of such a mixing term and use the result to compute the imprint of Σ-dependent coefficients on γ 3 .

A Spin-2 mixing interaction
In pure de Sitter space, linear interactions between two spin-2 fields of different mass can always be diagonalized, so a mixing of the kind desired is impossible. 27 However, inflationary spacetimes are only quasi-de Sitter, with the departure characterized by the slow-roll parameters, the first two of which are given by 27 This is most straightforward to see from the dual perspective: the two-point function of spin-2 fields is constrained by conformal invariance. In particular, invariance under special conformal transformations requires that the spin-2 fields have the same conformal weight (equivalently the same superhorizon time dependence). This requires the fields to have the same mass.
with t the proper time. By relaxing the constraint of special conformal invariance-but retaining all other de Sitter symmetry requirements-the leading order form of the mixing coefficient is restricted to take the form for some f (ε, η) that vanishes as the slow-roll parameters approach zero. 28 This form of mixing is realized by using, for instance, the following explicit interaction between a massless spin-2 (h (γ) µν ), a partially massless spin-2 (h (σ) µν ), and the inflaton (φ): 29 Putting φ on the inflationary background by replacing φ(x µ ) →φ(τ ), converting from h µν → σ ij , and finally going on-shell, one finds Pl and the standard slow-roll approximations ε ≈ constant and a(τ ) ≈ 1 H(−τ ) were used. The result (5.8) corresponds to (5.6) with f (ε, η) = πεM 2 Pl /Λ 2 . In the calculations to follow, we will take the form of T Σ to be where all slow-roll parameters and energy scales apart from Hubble have been absorbed into Λ mix . 30 In the context of the example (5.7), we have Λ 2 mix = πεM 2 Pl H 2 /Λ 2 . The mixing term (5.9) corrects the graviton power spectrum through diagrams of the form shown in Fig. 6. At leading order, the corrections come in the form Pl . (5.10) In order for our computations to be reliable, we want these corrections to the graviton propagator to be small, so we demand that Λ mix ≪ M Pl . 28 The k dependence follows from dilation invariance with ∆T = 3, ∆Σ = 2. The τ⋆-dependence in (5.6) follows from the late time behavior of the solutions (2.49) and the form of the other quadratic wavefunction coefficients. At late times γ ∼ τ 0 and σ ∼ τ and so we expect γσ ∝ τ⋆. This correlator is built from wavefunction coefficients via γσ ∝ Re T 2 −1 Re T Σ Re Σ 2 −1 and from (2.58) and (2.61) it follows that γσ ∝ τ 2 ⋆ Re T Σ , thus we expect Re T Σ ∝ τ −1 ⋆ . 29 We are assuming here that this interaction can be made gauge invariant under both the massless and partially massless spin-2 gauge transformations (perhaps by including other interaction terms, adjusting the overall coefficient, and/or having φ transform non-trivially under these gauge symmetries) and that the resulting theory admits a slowroll inflationary solution. Verifying this explicitly is beyond the scope of the present paper, however. See [64] for a different discussion of couplings between partially massless field and the inflaton. 30 The ∼ H −2 form of the prefactor was chosen to mimic the H −2 factors in T 2 and Σ 2 and is a convenient choice for later expressions.
. . Figure 6. The T Σ term also corrects the graviton power spectrum through diagrams of the above type.
These corrections are small in the limit Λ mix /M Pl ≪ 1.

Graviton bispectrum from Σ 3
We first consider the contribution to the graviton bispectrum coming from pure PM spin-2 interactions, generated by the Σ 3 wavefunction coefficient. It turns out that the induced bispectrum decays as τ ⋆ , since the contribution is of the schematic form (see Fig. 4) where in the last relation we have used (2.60), (2.63), and (5.9) to extract the time dependence of the two-point coefficients. As noted in Sec. 4.3, the leading component of Re Σ 3 is only O(τ −2 ⋆ ) and hence this contribution decays as τ ⋆ at late times, and will not produce an interesting signal in graviton non-Gaussianities.

Graviton bispectrum from ΣT 2
Next, we consider the contribution to γ 3 from interactions involving a single partially massless spin-2 and two massless gravitons, induced by the ΣT 2 wavefunction coefficient. The associated non-Gaussianity is rather strange: it grows as ∝ τ −1 ⋆ , but it also has purely the same shape as the ambiguous boundary term.
The contribution from ΣT 2 to γ 3 is of the schematic form (see Fig. 4) as follows from an analysis similar to that of the preceding section. As noted in Sec. 4.4, the leading component of Re ΣT 2 at late times grows as O(τ −2 ⋆ ), hence the Re ΣT 2 term would seem to produce non-Gaussianity which diverges as τ −1 ⋆ . Specifically, the following shapes are generated: As also noted in Sec. 4.4, these correspond to the local shapes which can be reproduced (or removed) by field redefinitions, so it is somewhat unclear what-if any-significance we should ascribe to the late-time divergence.

Graviton bispectrum from T Σ 2
In this section, we discuss the graviton bispectrum arising from 3-point interactions involving two partially massless fields and one massless graviton. This corresponds to the contribution of the T Σ 2 wavefunction coefficient to γ 3 . This is the only coefficient which sources the τ ⋆ -independent component of γ 3 , and is therefore the most phenomenologically-interesting situation.
The contribution from T Σ 2 to γ 3 is of the schematic form (see Fig. 4) as follows from an analysis similar to that of the preceding sections. The leading component of Re T Σ 2 scales with conformal time as O(τ −2 ⋆ ), hence T Σ 2 sources long-lived O(τ 0 ⋆ ) graviton non-Gaussianities. We find the following shapes: • The shape produced by the b 3a coefficient is: • The shape produced by the b 3b coefficient is: • The shape produced by the b 4a coefficient is: • The shape produced by the b 4b coefficient is: We also calculate the contribution to γ 3 from the purely boundary terms (4.28): The precise details of the shapes enumerated in this section are perhaps not particularly important at this juncture. What we would like to emphasize, however, is that they are different from the pure graviton shapes considered in Section 5.1. In this sense, a measurement of one of these shapes would be a sharp indication of the presence of an additional degree of freedom during inflation.

Soft limits
In the context of scalar non-Gaussianity, soft limits of correlation functions are a particularly sensitive probe of the presence of additional heavy particles [51,52,[83][84][85]. We are therefore motivated to see how the presence of additional spin-2 fields can show up in tensor 3-point functions.
In this section we discuss these soft limits-in particular the dominant contribution comes from bulk mixing diagrams that have been neglected up to this point, where the PM spin-2 converts into the massless graviton during inflation. In order to connect to familiar formulae in soft limits, we present results expressed in terms of arbitrary polarization tensors, rather than expressing equations in terms of ǫ X and ǫ P .

Heretofore neglected bulk diagrams
The mixing term produced by the interaction (5.7) also allows partially massless fields to imprint themselves on the T 3 coefficient via thus-far neglected diagrams of the form shown in Fig. 7. These diagrams are expected to contribute to non-Gaussianity at the same order as the contributions which have been presented so far. Unfortunately, the calculation of such contributions is difficult, as they require the use of bulk-to-bulk propagators for the internal lines. We only consider the single diagram in Fig. 7 for what follows, and leave the calculation of more complicated bulk-mixing calculations for future work.
For computations such as that shown in Fig. 7, only the spatial, transverse, traceless parts of the bulk-to-bulk propagator for the PM field are required. These components take on the form The function G(τ, τ ′ , k) can be derived as a Green's function for the equation of motion E τ,k G(τ, τ ′ , k) = −δ(τ − τ ′ ), where E is the differential operator appearing in the EOM: δS 2 [γ] δγ(τ,x) ij ≡ E τ,x γ ij (τ, x), subject to the following boundary conditions: Figure 7. The simplest bulk mixing diagram in which a PM field contributes to T 3 .
along with a jump condition on the first derivative that follows from integrating the equation of motion. The solution for G(τ, τ ′ , k) is given by Given (5.22), we can calculate the contribution shown in Fig. 7 to T 3 via (see Appendix A.2 for details) This is just the explicit form of the bulk interaction which generates the mixing term in (5.9), though there could be more general possibilities.

Soft-Limits from bulk mixing
We now compute (5.23) in the limit where one of the momenta is soft. Until now we have been symmetrizing over the three momenta appearing in T 3 . For the mixing diagram in Fig. 7 it is convenient to instead only symmetrize over k 2 and k 3 , and to keep k 1 as the leg where the mixing occurs, as indicated in Fig. 7. When exploring soft limits, this is useful since it lets us isolate the effect of taking the leg where the mixing occurs to be soft, versus the effect of taking one of the other two legs soft. Both calculations can be done analytically.
First, we take one of the non-mixing legs to be soft by sending k 3 → q, k 2 → k, and k 1 → −k−q. In the q → 0 limit, the leading terms are, 31 Tr (γ k ·γ q ·γ −k−q ) .
(5.25) where we have used matrix product notation to suppress the indices.
Next, we can take the mixing leg to be soft by sending k 1 → q, k 2 → k, and k 3 → −k − q. In this case, the leading terms are We can then convert these wavefunctional coefficients to squeezed bispectra. The wavefunctional coefficient (5.25) gives a contribution 32 while (5.26) contributes a logarithmic piece in the squeezed limit: The result (5.28) scales with a higher power of q than both (5.27) and subleading contributions to (5.27) which have been left unwritten. We have separated out the sub-leading terms (5.28) in order to highlight the logarithmic scaling in the soft limit, which is reminiscent of findings in "cosmological collider" studies [51,52]. In those cases, however, the logarithm appears as an oscillatory factor ∼ cos ln(q/k).

Other soft limits
In this section, we compute the graviton bispectrum soft-limit induced by the T 3 and T Σ 2 coefficients, which were the unique terms that sourced the time-independent component of γ 3 . 31 Explicit calculation demonstrates that both Tr (ǫ k · ǫ q · ǫ −k−q ) and Tr (ǫ k · ǫ −k−q )(k · ǫ q ·k) are O(1) while all other contractions are at least O(q/k). When wavefunction coefficients are converted into γ 3 correlators we effectively replaceγ → ǫ and we have used this fact to simplify expressions when deriving leading contributions. 32 We use the abbreviated notation γ k 1 γ k 2 γ k 3 T 3 Soft limits: There is a single soft limit that can be taken: (5.29) The soft limit from the boundary term shape (4.9) is: The above terms give the following contribution to the squeezed bispectrum: The structure on the first line is what was found in [1] and the structure in the second line is produced by the ambiguous boundary term operator (4.9), as shown. As discussed in Sec. 5.1, the Einstein-Hilbert shape corresponds to a 3 = M 2 Pl /4 and a 4 = a 5 = λ = 0, in which case we find the standard result which agrees with [1], after accounting for notational differences. 33 T Σ 2 Soft limits: There are two distinct soft limits which can be taken: • First, we can take the T leg to be soft: • Second, we can take one of the Σ legs to be soft: 33 Namely, factors of two arise when converting from our γij to the γ s of [1] due to the relations γij = s ǫ s ij γ s and γ s = 1 2 ǫ ij,s γij which follows from the standard polarization tensor conventions given in Appendix D. Figure 8. Contributions to T 4 from gravitons and partially massless fields.
The same analysis can be performed on the boundary term (4.27): • First, we can take the T leg to be soft: • Second, we can take one of the Σ legs to be soft: The above terms give the following contribution to the squeezed bispectrum: The first line of (5.37) is of the same form as the leading GR soft-limit result (5.32), though with a different coefficient, while the second line is produced by the ambiguous boundary term operator (4.27), as shown. Note that the dominant contributions to γ 3 arise from contributions where T is soft. We cannot conclude from (5.37) that this violates the usual soft limit result [86], since (5.37) is missing contributions from bulk-mixing diagrams.

Comments on the graviton trispectrum
Finally, we make some brief comments on the contribution of the coefficients we have calculated to the graviton trispectrum γ 4 . While the graviton trispectrum is expected to be even less observationally relevant than the bispectrum, for our purposes it has the advantage that partially γ 4 ⊃ + + Figure 9. Contributions to the graviton trispectrum γ 4 from gravitons and partially massless fields. The three-point vertices come from the T 3 and ΣT 2 coefficients we have already calculated. The four-point vertex depends on the T 4 coefficient which we have not calculated and whose associated diagrams are shown in Fig. 8.
massless fields can imprint themselves on γ 4 even in pure de Sitter, without a mixing term of the type considered in Sec. 5.2.1. A primary finding is that if the time-dependence of the ΣT 2 coefficients found in Sec. 4.4 is taken seriously, then these provide a contribution to γ 4 which diverges as ∝ τ −2 ⋆ . The leading-order contributions are due to interactions between one partially massless spin-2 particle and two gravitons. This interaction both produces the ΣT 2 coefficient considered in Sec. 4.4 and contributes to T 4 through the middle diagram in Fig. 8. These coefficients, along with T 3 , determine γ 4 at leading order through the diagrams in Fig. 9.
In order to understand the effects of partially massless spin-2 exchange, it is useful to compare and contrast the contributions of T 3 and ΣT 2 to γ 4 : • The leftmost diagram in Fig. 9 is schematically given by using the scaling of these wavefunction coefficients, we see that this is a τ ⋆ -independent expression scaling as γ 4 ∼ k −9 .
• The middle diagram in Fig. 9 is schematically given by appear due to the presence of new possible on-shell cubic vertices involving the massive spin-2 fields. This can also be understood from the dual perspective: massive spin-2 fields correspond to nonconserved spin-2 currents as conformal representations, for which there are more possible 3-point structures. This provides evidence for the presence of new particles during inflation: measurement of a shape not allowed for a massless graviton would point toward other spinning degrees of freedom.
We have derived the general bulk-to-boundary propagators and on-shell cubic vertices required to compute arbitrary tensor non-Gaussianities involving spin-2 fields of any mass. In order to obtain closed-form results, we have focused on the case where the additional spin-2 field is partially massless. These interactions also provide a mechanism for increasing the size of tensor non-Gaussianity while only negligibly affecting the tensor power spectrum.
The tensor bi-spectrum is expected to be very small, and even the tensor spectrum itself has not been detected yet, so any observation is necessarily futuristic. However, it is possible that information about the tensor bi-spectrum may be gained from the future LISA gravitational wave experiment [12]. Measurement of a non-Einstein-Hilbert 3-point correlation function for the graviton would be exciting for a number of reasons. Observation of one of the shapes we enumerate in this paper (or their general mass versions) would provide evidence for the presence of an additional heavy field during inflation, but it would also strongly suggest that there should be additional heavy states to discover. In [87], it was shown that large higher-derivative cubic interactions for a massless graviton in flat space lead to asymptotic causality problems. These problems can be cured by introducing a tower of higher-spin states-as in string theory. If we assume that something similar happens on de Sitter space (though this is not proven), we should expect that large higher-derivative non-Gaussianity for the graviton would also imply the presence of a tower of higher-spin states [2,87]. The shapes that we are interested in here are not pure gravity shapes, but rather involve other spin-2 fields: however, it was shown in [76,77]-again in flat space-that theories involving massive spin-2 fields with large cubic couplings have the same superluminality issues, unless their couplings are of the Einstein-Hilbert form. Therefore, a natural expectation is that non-Einstein-Hilbert non-Gaussian shapes in general require a tower of higher-spin fields to be consistent.
It would be very interesting to understand more fully to what extent the signatures we have computed here are degenerate with either modifications of the initial vacuum or to slow-roll suppressed corrections to graviton correlation functions. Our expectation is that at least the non-analytic exchange contribution cannot be mimicked by any local interaction or non-singular vacuum state, but it would be interesting to understand these issues systematically. We leave such a study for the future.
In this paper we have focused on the "bulk" version of the wavefunctional calculation, by directly computing the on-shell action in de Sitter space. However, our final results are controlled by the isometries of de Sitter space, which act as conformal transformations at late times. Conformal symmetry is known to fix the 3-point correlation function of any operators, up to a finite number of constants. Therefore, it should be possible to re-derive our results from the "boundary" perspective using conformal symmetry. Though this seems to be a formidable task, given that the constraints of conformal symmetry for spinning operators in momentum space are quite complex, the final results we obtain are rather simple: they are rational functions of momenta, so it may be more tractable than it naively seems. It would be very interesting to carry out such an analysis, as it may provide a systematic way to compute to higher-order correlation functions involving massive and massless spinning operators in cosmology and in AdS/CFT applications.

A The wavefunction
One of the central objects of our interest is the late-time cosmological wavefunction. For the convenience of the reader, in this Appendix we briefly review the wavefunctional approach to cosmological correctors. This is a straightforward extension of flat space Schrödinger field theory. For some other applications, see [81,[89][90][91] A.1 Introduction to Ψ Consider a theory of a scalar field ϕ, for concreteness, whose action is of the form where L int stands for any interaction terms we may be interested in. The generalization to other types of fields is straightforward.
The wavefunction of the universe is simply the amplitude for finding ϕ(x µ ) in a particular spatial configurationφ(x) at time τ ⋆ : Ψ[φ, τ ⋆ ] = φ|Ψ(τ ⋆ ) . This is the field-basis representation of a state in the Hilbert space. Given Ψ[φ, τ ⋆ ], arbitrary equal time expectation values can be calculated via the usual quantum mechanics formula There are two ways in which Ψ[φ, τ ⋆ ] may be practically calculated. First, from S one can build the hamiltonian H[ϕ, Π ϕ ], where Π ϕ is the momentum conjugate to ϕ, and solve the functional Schrödinger equation This quickly becomes infeasible, so a more tractable approach is to use the path integral where we have suppressed the dependence on the initial state; typically we project onto the vacuum state of the theory by employing an iǫ prescription as we send τ → −∞.
A.2 The path integral calculation of Ψ In order to compute the late-time wavefunction in de Sitter space, we utilize the path integral method and restrict all of our calculations to the tree-level, classical approximation in which we can approximate the path integral by its saddle point so that where ϕ cl is a solution to the classical equations of motion and we have omitted the overall normalization factor (it is immaterial). The particular classical solution, ϕ cl (τ, x), that we are interested in is the one satisfying the boundary conditions The choice of ϕ vac will be discussed momentarily.
Using these ingredients, an implicit expression for the classical solution is given by .
This can be solved iteratively to find ϕ cl to any desired order inφ.
The on-shell action S cl [φ cl ] is efficiently calculated by first writing the quadratic action as a boundary contribution plus the equation of motion: where n µ is the outward pointing normal vector to the τ = τ ⋆ surface 34 andḡ ij is the induced metric on this hypersurface. Plugging in the expression (A.8) for ϕ cl into the above and doing some algebra, the end result is 35 Here (and throughout this paper), all temporal integrals range from −∞ to τ ⋆ , while spatial integrals range from −∞ to ∞. In the present paper, we will primarily focus on three-point correlation functions, in which case the final line of (A.11) does not contribute.

A.3 Correlators from Ψ
Once we have the wavefunction, there is yet another step required to extract from it correlation functions. The wavefunction, Ψ[φ, τ ⋆ ], can be written as an expansion inφ. We are interested in momentum space correlators, in which case we write where the O k 1 . . . O kn 's are simply functions of the indicated arguments. These coefficients contain the information that specifies correlation functions of the fields ϕ. Due to homogeneity and isotropy of τ = const. slices, these wavefunction coefficients are proportional to delta functions which enforce d-momentum conservation. We use a prime to indicate that a delta function and the associated 2π factors have been dropped from a quantity, i.e., Expectation values are related to wavefunction coefficients via standard perturbative Gaussian integral formulas, such as Fall-off conditions are chosen to ensure that ϕ vanishes on all other (spatial) boundaries. 35 Deriving this requires the relation lim τ ′ →τ⋆ √ḡ n µ ∇ (τ ′ ,y) µ G(τ, x; τ ′ , y) = K(τ, x; y) , (A. 10) which follows from the integral identity M K G − G K = ∂M n 2 (Kn · ∇G − Gn · ∇K) where n µ the outward pointing normal vector. This property is responsible for important cancellations in the derivation of (A.11). (A.14) where k ij ≡ k i + k j . When generalized to tensor fields, inverse wavefunction factors become matrix inverses. Note that the real components of the O n wavefunction coefficients arise from the imaginary component of the action because Ψ ∼ e iS . Though the action is defined in terms of real fields, it develops imaginary terms because it is being evaluated on complex field configurations. This is essentially due to the early time vacuum condition, which is an imaginary constraint.

B Projectors
In several places it is useful to be able to project spatial tensors onto their irreducible components. We therefore introduce a basis of projection tensors which accomplishes this task. It is most straightforwardly built out of the transverse tensor Using this, we form the tensors (dropping the k label on π ij to simplify notation) These tensors obey the usual orthonormality and completeness relations: in condensed notation. 36 36 For instance, the identity matrix 1 is explicitly given by (1) ij lm = δ

C On-shell, gauge invariant cubic interactions
In this Appendix, we describe the derivation of the on-shell cubic vertices involving massless or partially massless fields, which require the imposition of on-shell gauge invariance. A similar problem was treated in [49,92] making use of embedding space techniques, but for our purposes it will be more convenient to derive expressions directly in the physical (d + 1)-dimensional de Sitter space.
The bases of operators introduced in Sec. 3 (see (3.2), (3.3) and (3.4)) are written for spin-2 fields with generic masses. When a h (i)µν field is massless or partially massless, gauge invariance places further constraints on the operator combinations. In this Appendix, we discuss how gauge invariance is enforced in the on-shell action. We produce the counting of operators found in [49] for cases where d is arbitrary and find additional invariant operators for special choices of d, namely d = 3.

C.1 Off-shell vs. on-shell gauge invariance
Given a theory involving a spin-2 field, h µν , which enjoys a gauge symmetry, we can expand both the action and gauge transformation in powers of h: ξ h µν are O(h n ) and ξ represents the dependence on the gauge parameter. The off-shell statement of gauge invariance is that S[h] = S[h ′ ] which can be written order-by-order in the fields as where we have defined δ All of the conditions (C.2) mean equality up to total derivatives in the action.
However, to compute 3-point correlation functions, we only require the on-shell action. In this case, the constraints simplify. In particular, when the fields are put on-shell, the second condition in (C.2) becomes where ∼ = indicates equality up to terms which vanish on-shell. This happens because δS 2 /δh is the linear equation of motion, which is zero on-shell to this order in the fields. The on-shell gauge invariance constraints are in some ways simpler than the off-shell ones since the knowledge of the first-order gauge transformation δ (1) ξ is not required. However, the requirement of keeping all fields on-shell also brings in complications which are discussed in the following section using a concrete example.
C.1.1 Example: the Einstein-Hilbert interaction in d = 3 There are two main subtleties in solving (C.4): the form of the gauge parameter is restricted, since it must preserve ∇ µ h µν = h = 0, 37 and there exist terms which are total derivatives only when the EOM are used (and hence are not annihilated by the naive variational derivative). We demonstrate both of these subtleties and how to deal with them in the concrete case of a massless spin-2 in d = 3, where we derive the on-shell cubic vertex corresponding to the Einstein-Hilbert action (2.2).
On-Shell Gauge Transformations: Generic gauge transformations do not preserve the conditions ∇ µ h µν = h = 0. Under the massless gauge transformation h µν → h ′ µν = h µν + ∇ (µ ξ ν) (2.4) we instead have in d = 3: Therefore, after imposing these conditions, we can only demand on-shell gauge invariance under residual gauge transformations that obey the equations On-Shell Total Derivatives: The condition δ is a sum of terms which are total derivatives after accounting for equations of motion. We can parametrize the most general possible total derivative as and O µν built from appropriate powers of ξ µ and h µν . If we were solving for the off-shell gauge invariance constraints, only the explicit total derivative term, ∇ µ J µ , would appear, but due to the on-shell conditions (2.5), (C.7), and ∇ µ h µν = 0, the additional terms above are also total derivatives.
The appearance of the extra terms in (C.8) makes the process of solving for the gauge-invariant interaction somewhat different than in the off-shell scenario. A practical method for solving (C.8) is to integrate all derivatives off of the gauge parameter ξ µ . This is accomplished by acting on both sides with the operator d 4 x ξ µ (x) δ δξµ(x) , yielding 38 for some operators O ′ , O ′ µ . In contrast, the analogous off-shell gauge invariance constraint at this order would read with nothing appearing on the RHS.
Finding the Einstein-Hilbert Interaction: Accounting for the previous points, we can now solve for the d = 3 Einstein-Hilbert interaction. From derivative counting, the on-shell action must take on the form (3.4) with a 4 = a 5 = 0: Calculating the LHS of (C.9) for this action yields and we need to find the operators O and O ν for which The relation (C.12) is then solved at each order in derivatives, writing O = n O (n) where O (n) is O(∇ n ) and similar for O ν . Starting with the O(∇ 3 ) terms, we can systematically include all on-shell independent terms in O and O ν in (C.12) which contribute at this order with arbitrary coefficients. These are explicitly given by Solving (C.12) at this order places constraints both on the O and O ν coefficients and the a i 's: The remaining terms in (C.12) are O(∇) and read Choosing O (0) = (a 3 − 5A 3 )H 2 h µν h µν cancels the first term above, while the second cannot be canceled so we are forced to set a 3 = a 1 /(3H 2 ).
The interaction has therefore been fixed up to an overall factor: Expanding the Einstein-Hilbert action (2.2) to cubic order, integrating by parts, and using the onshell conditions, it can be explicitly shown that the cubic interactions agree if we set a 1 = 3 4 M 2 Pl H 2 . In the following section we describe the generalization of this procedure to general dimension and to include partially massless fields.

D Explicit d = 3 transverse-traceless polarization tensors
In this Appendix, we discuss the explicit transverse, traceless polarization tensors which are used when evaluating tensor non-Gaussianity in d = 3. We write our basis as ǫ kaP ij = z i z j − u a i u a j , ǫ kaX ij = u a i z j + z i u a j , (D.1) following [2], though changing the overall normalization. In (D.1), z i is a unit vector orthogonal to the plane spanned by the three k a 's and u a is a unit vector in the plane of this triangle which is orthogonal to k a ; see Fig. 10. The action of parity is reflection across the plane in which the momenta lie. The u a are invariant under this transformation, while z flips sign. Hence, ǫ P is even under parity, while ǫ X is odd.
The normalization in (D.1) is chosen such that s ǫ * , s kaij ǫ s kalm = 2(Π ka T T ) ijlm , ǫ sij ka ǫ * , s ′ kaij = 2δ ss ′ , (D .2) where s ∈ {X, P } and Π T T is the transverse, traceless projector defined in (B.2). The factors of 2 are conventional in the literature. For instance, the dimensionless tensor power spectrum ∆ 2 γ is naturally defined as and using (2.60), which reads γ k Pl k 3 (Π k T T ) ijlm and (Π k T T ) ij ij = 2, it follows that Pl which is the standard expression.

D.1 Dimensionally-dependent identities
A motivation for introducing an explicit basis of polarizations is that in d = 3, there are many accidental identities involving contractions between transverse, traceless polarization tensors and 3-momenta, which can make it difficult to identify the independent structures appearing in contractions such as T ij no . As an illustration of this difficulty we present some highly non-trivial identities relating various contractions of momenta and polarizations. After using the transverse property and 3-momentum conservation to canonicalize all contractions between k a 's and ǫ's by eliminating all contractions of the form k i a ǫ k b ij with a < b, cyclically understood, we find the following identities via direct calculation. First, we find ǫ i k 1 j ǫ j k 2 l ǫ l k 3 i = −4 k T a (k T − 2k a ) k i 1 k j 1 ǫ k 3 ij k l 2 ǫ k 1 lm ǫ k 2 m n k n 3 + 2 perms. .

E Integrations by parts and the wavefunction
Perturbative computation of the wavefunction in terms of the on-shell action appears to have an inherent ambiguity: the addition of boundary terms supported on the time slice τ = τ ⋆ . 40 Actions related by integrations-by-parts differ precisely by such boundary terms, so it behooves us to understand their effect on the wavefunctional. In this Appendix, we discuss the effect that integrations-by-parts have on wavefunction coefficients. These operations do not leave the wavefunction invariant. Instead, they change Ψ in a very specific way: integration-by-parts has the same effect as local field redefinitions. We illustrate these concepts using the example of a massless scalar field. See [94][95][96] for discussions of boundary terms in the canonical in-in formalism and [97] for a discussion from a dS/CFT and wavefunctional perspective.

E.1 Example: massless scalar field
Consider the action Each of the a i operators in (E.1) can be integrated-by-parts to be proportional to φ. Hence, when evaluated on the classical linear solution obeying ϕ cl = 0, the entire effect of these operators is captured by a boundary term: d d x √ḡ n µ a 1 2 ϕ 2 cl ∇ µ ϕ cl + a 2 2 (∇ϕ cl ) 2 ∇ µ ϕ cl .
(E.2) The boundary terms above affect Ψ non-trivially since ϕ cl is non-vanishing on the τ = τ ⋆ surface. However, their contributions to the wavefunction have the same effect as redefiningφ in Ψ[φ, τ ⋆ ].
We illustrate this by explicitly calculating the wavefunction corresponding to (E.1) in the d = 3, super-horizon limit. The linear solution is Substituting this solution into (E.2), the result is (E.4) 40 More specifically, the only boundary terms which can have interesting effects on the wavefunction are those involving derivatives normal to the boundary. Operators which only depend on tangential derivatives (i.e. the analogues of holographic counterterms [93]) can only affect the phase of the wavefunction, due to the reality condition on fields at the τ = τ⋆ surface.
The boundary contribution to (E.4) is non-zero, but it can also be removed by a local field redefinition, due to the ∼ k 3 form of the quadratic terms in (E.4). Specifically, making the following local, position space field redefinition precisely cancels the contribution of boundary terms to the wavefunction:φ (x) →φ(x) + a 1 2φ (x) 2 . (E.5) A systematic construction of all possible cubic boundary terms involving up to five derivatives confirms that total derivatives only generate the ∼ i k 3 i φ k 1φ k 2φ k 3 structure. We have found that the same pattern holds for spin-2 wavefunctions: integration-by-parts ambiguities in the wavefunction are of precisely the same form as field redefinition ambiguities.
E.2 A general cubic argument: field redefinitions correspond to boundary terms In fact, it is possible to show more generally that local field redefinitions of the type ϕ(x µ ) → ϕ(x µ ) + a ϕ(x µ ) 2 , with constant a, have the same effect on the cubic order wavefunction as adding boundary terms to the action, so that the ambiguities associated with boundary terms are the same contact term ambiguities associated with making a choice of field variables.

Consider a generic cubic action
where ϕ stands for one or more fields of arbitrary type. Making the above field redefinition yields explicitly demonstrating that the result of a local field redefinition is to add a boundary term to the on-shell action. The terms on the RHS of (E.11) only affect three-point correlators by contact terms, where two of the three operators are brought to coincident points, i.e. their contribution is proportional to a delta function. Such semi-local contributions, in the language of [60], are often ignored in position space, but can be relevant for cosmological correlators.

F CFT two-point functions
In this Appendix, we review the structure of momentum-space CFT two-point functions for arbitrary symmetric, traceless primary operators. For fields on de Sitter, both the quadratic wavefunction coefficients and the two-point equal-time correlators take on the form of CFT correlators.
An efficient method for deriving the momentum space two-point functions of generic spin, symmetric, traceless primary fields ϕ i 1 ...is is to first contract all indices with auxiliary null-vectors y and z to create the index-free correlator The constraints of scale and special conformal invariance can then be efficiently imposed on ϕ k s ϕ −k s ′ which is forced to take the form [51] ϕ k s ϕ −k s ′ ∝ k 2∆−d−2s (y · k) s (z · k) s 2 F 1 −s, −1 + ∆, 1 − d 2 − s + ∆, k 2 2 (y · z) (y · k)(z · k) (F.2) with 2 F 1 a hypergeometric function. 41 The y and z factors can then be stripped from the above to produce a traceless correlator by acting s-times with the operator [73] and s-times with the similar D (z) i operator.
The above construction also reproduces the quadratic wavefunction coefficients and two-point functions derived in Sec. 2.4. In the present paper we are primarily concerned with the transverse, traceless parts of correlators, but the other components can be restored using the preceding results. 41 The 2F1 convention we use is that of functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1.