{'fraction_non_alphanumeric': 0.04491534372909993, 'fraction_numerical': 0.017506189353139263, 'mean_word_length': 4.110817941952506, 'pattern_counts': {'":': 1, '<': 20, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-Gaussianity as a Particle Detector

We study the imprints of massive particles with spin on cosmological correlators. Using the framework of the effective field theory of inflation, we classify the couplings of these particles to the Goldstone boson of broken time translations and the graviton. We show that it is possible to generate observable non-Gaussianity within the regime of validity of the effective theory, as long as the masses of the particles are close to the Hubble scale and their interactions break the approximate conformal symmetry of the inflationary background. We derive explicit shape functions for the scalar and tensor bispectra that can serve as templates for future observational searches.

Since massive particles decay outside of the horizon during inflation, they cannot be observed directly in late-time correlation functions. Instead, the presence of massive particles has to be inferred from their indirect effects on the correlation functions of ζ = −Hπ and γ ij (see Fig. 1). Some of these effects can be mimicked by adding a local vertex in the low-energy effective Lagrangian, which is the result of integrating out the heavy fields. On the other hand, massive particles may spontaneously be created in an expanding spacetime [3][4][5], an effect which cannot be represented by adding a local vertex to the effective Lagrangian [6]. The role of these nonlocal effects as a means of detecting massive particles during inflation was recently highlighted by Arkani-Hamed and Maldacena (AHM) [6]: the spontaneous particle creation allows us to probe massive fields during inflation, even though we are only observing the late-time expectation values of light fields. The rate of particle production in de Sitter space is exponentially suppressed as a function of mass, e −m/T dS , with T dS ≡ H/2π, so their imprints will only be detectable if their masses are not too far above the Hubble rate H. 2 Since the inflationary scale may be as high as 10 14 GeV, this nevertheless provides an opportunity to probe massive particles far beyond the reach of conventional particle colliders.
Nonlinearities in the decay of the massive particles lead to a non-Gaussianity in the late-time correlation functions of ζ and γ ij . The form of this non-Gaussianity will depend on the masses and the spins of the extra particles. The effects of additional scalar fields during inflation have been explored in many previous works, e.g. in the context of quasi-single-field inflation [8][9][10]. A characteristic signature of these fields are non-analytic scalings in the soft momentum limits of the non-Gaussian correlation functions. These soft limits are particularly clean detection channels, since in single-field inflation their momentum scalings are fixed by the symmetries of the inflationary background [11,12]. The most straightforward interpretation of such non-analyticity in the correlation functions is therefore the presence of extra particles. Scalar fields with masses less than 3 2 H give rise to monotonic scalings in the squeezed limit [8,9], while those with masses greater than 3 2 H lead to oscillatory behavior [6,10,13,14]. The effects of extra massive particles with spin have not been studied in as much detail. Such particles can naturally arise as massive Kaluza-Klein modes or as part of the tower of higher-spin states from string theory [15,16]. As was shown by AHM, the spins of new particles lead to a distinctive angular dependence of the soft limits of the non-Gaussian correlators. The analysis of AHM was restricted to the squeezed limit of the bispectrum and interactions that maintained the approximate conformal invariance of the inflationary background. While this assumption made their analysis particularly well controlled, it also implied that the amplitude of the signal is highly suppressed and only observable in the most optimistic and futuristic scenarios.
We will drop some of the restrictions of the analysis of AHM in our analysis. In particular, we will allow for a large breaking of conformal invariance within the framework of the effective field theory (EFT) of inflation [17]. We will find that the signal due to massive spinning particles can be observable within the regime of validity of the EFT. At the same time, the main spectroscopic features of particles with spin during inflation do not rely on conformal invariance and therefore still apply. On the other hand, couplings to particles with odd spins, which are disallowed in the conformally-invariant case, are permitted in the generic effective theory. We also consider the breaking of special conformal invariance by giving the Goldstone fluctuations a nontrivial sound speed. In that case, we find a reduced exponential suppression in the particle production rate, and thus an enhanced level of non-Gaussianity. Finally, we also study the coupling to an external graviton γ ij . We demonstrate that the soft graviton limit of the correlator γζζ provides an interesting detection channel for extra particles. Like in the case of massive scalar fields, there will be non-analytic scalings of non-Gaussianities close to the soft momentum limit, but this time only for particles with spin greater than or equal to two.
Outline In this paper, we analyze the allowed couplings of massive particles with spin to the Goldstone boson of broken time translations and the graviton, and discuss their observational signatures. In Section 2, we first collect the equations of motion for massive fields with spin in de Sitter space, whose solutions are presented in Appendix A. In Section 3, we then construct the effective action for the leading interactions between the Goldstone boson π, the graviton γ ij , and massive spinning fields σ µ 1 ...µs . We analyze under what conditions the theory is under perturbative control and discuss various constraints on the sizes of the couplings. In Section 4, we compute the correlation functions associated with the interactions of Section 3. We estimate the maximal amount of non-Gaussianity that is consistent with the constraints on the couplings of the effective theory. Details of the in-in computation are relegated to Appendix B, and analytic results for soft limits are given in Appendix C. Our conclusions are presented in Section 5.
Notation and conventions We will use natural units, c = = 1, with reduced Planck mass M 2 pl = 1/8πG. Our metric signature is (−+++). We will use Greek letters for spacetime indices, µ, ν, . . . = 0, 1, 2, 3, and Latin letters for spatial indices, i, j, . . . = 1, 2, 3. Three-dimensional vectors are written in boldface, k, and unit vectors are hatted,k. A shorthand for the symmetrization of tensor indices is a (µ b ν) ≡ 1 2 (a µ b ν + a ν b µ ). Overdots and primes will denote derivatives with respect to physical time t and conformal time η, respectively. The letter π will refer both to 3.141 . . . and the Goldstone boson of broken time translations. The dimensionless power spectrum of a Fourier mode f k is defined as where the prime on the expectation value indicates that the overall momentum-conserving delta function has been dropped.

Spin in de Sitter Space
We begin by reviewing a few elementary facts about massive fields with spin in four-dimensional de Sitter space, dS 4 .

Spin-1
The quadratic action of a massive spin-1 field σ µ in de Sitter space is where m 2 1 ≡ m 2 + 3H 2 , with m being the mass of the field. 3 The structure of the action (2.1) is uniquely fixed by requiring the absence of ghost degrees of freedom. 4 Up to integration by parts, this is equivalent to the Proca action. Variation of the action yields the equation of motion, σ µ − ∇ µ ∇ ν σ ν − m 2 1 σ µ = 0, where ≡ ∇ µ ∇ µ denotes the Laplace-Beltrami operator on dS 4 . Taking the divergence of this equation gives the constraint ∇ µ σ µ = 0. The on-shell equation of motion then becomes − m 2 1 σ µ = 0 . (2.2) In Appendix A, we derive the solutions to this equation for the different helicity components of the field.

Spin-2
The unique ghost-free quadratic action of a massive spin-2 field σ µν in de Sitter space is [19] 3) 3 We define the mass parameter in such a way that the action acquires a gauge invariance in the massless limit, m = 0. This is required in order for massless spinning fields to propagate the right number of degrees of freedom. The mass defined in this way can also be identified as the mass of the field in the flat space limit [18]. 4 The ghost-free structure of the quadratic action will remain valid as long as nonlinear interactions can be treated perturbatively.
Spin-s The Lagrangian for massive fields with arbitrary spin in flat space was constructed in [21,22], and generalized to (A)dS spaces in [23]. For massive fields with spin greater than 2, the action is rather complex and requires introducing auxiliary fields of lower spins. An alternative, which we will follow, is to use a group theoretical approach to find the equations of motion directly [24]. A massive bosonic spin-s field is described by a totally symmetric rank-s tensor, σ µ 1 ···µs , subject to the constraints ∇ µ 1 σ µ 1 ···µs = 0 , σ µ 1 µ 1 ···µs = 0 . where m 2 s ≡ m 2 − (s 2 − 2s − 2)H 2 . The shift in the mass arises from the mismatch between the Casimir and Laplace-Beltrami operators in de Sitter space and is necessary to describe the correct representations for massless fields. Equivalently, it is required by imposing gauge invariance in the massless limit, m = 0. Solutions to equation (2.7) are obtained in Appendix A.
Following Wigner [25], we identify the spectrum of particles by the unitary irreducible representations of the spacetime isometry group. For the de Sitter group SO (1,4), these representations fall into three distinct categories [26,27]: principal series complementary series discrete series for s, t = 0, 1, 2, ..., with t ≤ s. Masses that are not associated with one of the above categories are forbidden and correspond to non-unitary representations. At the specific mass values corresponding to the discrete series, the system gains an additional gauge invariance and some of the lowest helicity modes become pure gauge modes; this phenomenon is called partial masslessness [28].
The spectrum of massive particles is contained in the principal and complementary series. We see that unitarity demands the existence of a lower bound, m 2 > s(s − 1)H 2 , on the masses of fields that belong to this spectrum. For s = 2, this is known as the Higuchi bound [19].
In the late-time limit, the generators of the de Sitter isometries form the 3-dimensional conformal group. The asymptotic scaling of a spin-s field is where the conformal weight of the field is defined as 6 (2.9) In this paper, we will deal mostly with particles belonging to the principal series which covers the largest mass range and corresponds to µ s ≥ 0. For real µ s , the asymptotic scaling is given by a complex-conjugate pair, resulting in a wavefunction that oscillates logarithmically in conformal time. The complementary series has imaginary µ s and corresponds to the interval −iµ s ∈ (0, 1/2). In that case, only the growing mode survives in the late-time limit.

Spin in the Effective Theory of Inflation
In this section, we will construct the leading interactions between the Goldstone boson of broken time translations π, the graviton γ ij , and massive spinning fields σ µ 1 ...µs . We start, in §3.1 and §3.2, by reviewing the effective actions for the Goldstone boson and the graviton. In §3.3, we introduce the couplings to massive particles with spin; first for the special cases s = 1 and 2, and then for arbitrary spin. We close, in §3.4, by discussing how large the mixing interactions can be made while keeping the effective theory under theoretical control.

Goldstone Action
A time-dependent cosmological background induces a "clock", i.e. a preferred time slicingt(t, x) of the spacetime. In the inflationary context, surfaces of constantt may be associated with the homogeneous energy density of the background. The slicing has a timelike gradient, and the unit vector perpendicular to the surface of constantt is The induced spatial metric on the slicing is h µν ≡ g µν +n µ n ν . The metric h µν also serves to project spacetime tensors onto the spatial hypersurfaces. Geometric objects living on the hypersurfaces can be constructed from h µν and n µ . Examples are the intrinsic curvature, (3) R µνρσ [h], and the extrinsic curvature, K µν ≡ h (µ ρ ∇ ρ n ν) . Using the Gauss-Codazzi relation, the intrinsic curvature can be written in terms of (the projection of) the four-dimensional Riemann tensor R µνρσ and the extrinsic curvature K µν . Higher-derivative objects can be constructed using the covariant derivative ∇ µ , defined with respect to h µν .
In unitary gauge, the time coordinate t is chosen to coincide witht. Fluctuations in the clock have been eaten by the metric, and the effective action for adiabatic fluctuations only contains metric perturbations. The action does not have to respect full diffeomorphism invariance, but only has to be invariant under time-dependent spatial diffeomorphisms, x i → x i + ξ i (t, x). Besides terms that are invariant under all diffeomorphisms (such as curvature invariants like R and R µνρσ R µνρσ ), the reduced symmetry of the system now allows many new terms in the action. The normal vector in (3.1) becomes n µ = −δ 0 µ /(−g 00 ) 1/2 in unitary gauge. By contracting covariant tensors with n µ , we produce objects with uncontracted upper 0 indices, such as g 00 and R 00 . It is easy to check that these are scalars under spatial diffeomorphisms. Functions of g 00 , R 00 , etc. are therefore allowed in the effective action. In general, products of any fourdimensional covariant tensors with free upper 0 indices are allowed operators. In addition, we can have operators made out of the three-dimensional quantities describing the geometry of the spatial hypersurfaces (e.g. K µν ). The most general action constructed from these ingredients is [17] where the only free indices entering the functional L are upper 0's. The spacetime indices in (3.2) are contracted with the four-dimensional metric g µν . Terms involving explicit contractions of the induced metric h µν do not lead to new operators.
At leading order in derivatives, the action can be written in terms of g 00 alone, where δg 00 ≡ g 00 + 1. The coefficients of the operators 1 and g 00 have been fixed by the requirement that we are expanding around the correct FRW background with a given expansion rate H(t). This removes all tadpoles and the action starts quadratic in fluctuations. Because time diffeomorphisms are broken, all operators are allowed to have time-dependent coefficients. The limit of slow-roll inflation corresponds to M n → 0.
To make the dynamics of the theory defined by (3.3) more transparent, we introduce the Goldstone boson associated with the spontaneous breaking of time-translation invariance. Through the Stückelberg trick, the field π also restores the full diffeomorphism invariance of the theory. Specifically, we perform a spacetime-dependent time reparameterization, t →t = t + π(t, x). The metric transforms in the usual way: e.g.
Substituting this into (3.3) gives the action for the Goldstone boson. In general, this action contains a complicated mixing between the Goldstone mode and metric fluctuations. However, for most applications of interest, we can take the so-called decoupling limit, and evaluate the Goldstone action in the unperturbed background [17], g µν →ḡ µν . In this case, the transformation (3.4) reduces to g 00 → −1 − 2π + (∂ µ π) 2 , and the Goldstone Lagrangian becomes We see that M 2 = 0 induces a nontrivial sound speed for the Goldstone boson, A small value of c π (large value of M 2 ) is correlated with an enhanced cubic interactionπ(∂ i π) 2 . The Planck constraints on primordial non-Gaussianity imply c π ≥ 0.024 [29]. For purely adiabatic fluctuations, the relationship between the comoving curvature perturbation ζ and the Goldstone boson is ζ = −Hπ + O(π 2 ). The dimensionless power spectrum of ζ is found to be where f 4 π ≡ 2M 2 pl |Ḣ|c π is the symmetry breaking scale [30]. The observed amplitude of the power spectrum is ∆ 2 ζ = (2.14 ± 0.05) × 10 −9 [31].

Graviton Action
The tensor sector of inflation is harder to modify [32]. The leading correction to the Einstein-Hilbert action can be written as where the combination of extrinsic curvature tensors was chosen in a way that doesn't induce a scalar sound speed. Inserting the transverse and traceless tensor perturbation of the metric, g ij = a 2 (δ ij + γ ij ), we find where we have defined the tensor sound speed (3.10) As discussed in detail in [32], the tensor sound speed can be set to unity by a disformal transformation. This transformation makes the tensor sector canonical, and moves all the corrections to the scalar sector. In this paper, we will make the same choice of frame and work with c γ = 1 throughout. The dimensionless power spectrum of γ ij is then given by The current contraint on the tensor-to-scalar ratio, r ≡ ∆ 2 γ /∆ 2 ζ < 0.07 [33], implies that ∆ 2 γ 1.5 × 10 −10 .

Mixing Interactions
Next, we construct the effective action for interactions between the Goldstone boson, the graviton, and massive spinning fields. We will also consider self-interactions of the massive spinning fields, and focus on terms which contribute to the correlation functions ζζζ and γζζ at tree level and at leading order in derivatives. Moreover, we will restrict our presentation to the subset of interactions which give rise to a distinctive angular dependence due to the exchange of the spinning fields.

Couplings to the Goldstone
The construction of the effective action proceeds as above. We first write down all operators consistent with the symmetries. Amongst them will be tadpole terms, which must add up to zero.
In unitary gauge, the basic building blocks involving spinning fields are σ 0···0 and all Lorentzinvariant self-interactions, e.g. σ µ 1 ···µs σ µ 1 ···µs . The latter are invariant under all diffeomorphisms, so they don't lead to a coupling to π after the Stückelberg trick, whereas the former transform as We may also have contractions with the curvature tensors, which appear at higher order in derivatives.
Spin-1 We first analyze the couplings between a massive spin-1 field σ µ and the Goldstone boson π. In unitary gauge, the operators of the effective action involve g 00 and σ 0 . In order not to alter the background solution, these operators have to start at quadratic order in fluctuations.
• At leading order in derivatives and to linear order in σ µ , the mixing Lagrangian is 7 Introducing π using (3.4) and (3.12), we get where we have taken the decoupling limit so that couplings to metric fluctuations become irrelevant. 8 We also used the constraint ∇ µ σ µ = 0, which we assume to hold at the background level, to replaceπσ 0 by ∂ i πσ i . Since only the cubic mixingπ∂ i πσ i will lead to the characteristic angular structure in the resulting correlation functions (see §4.2), we will focus on the bispectrum created created by the combination ofπ∂ i πσ i and ∂ i πσ i . Note that there is a single parameter ω 1 controlling the size of these two interactions. This is a consequence of the nonlinearly-realized time translation symmetry.
• At quadratic order in σ µ , the mixing Lagrangian is (3.16) where in the second line we have introduced the Goldstone and taken the decoupling limit. We see that, this time, the size of the cubic interactionπσ i σ i is independent from the quadratic mixing term.
Combining the above, we can write where π c ≡ f 2 π π is the canonically normalized Goldstone boson, and we defined (3.18) We note that ρ 1 and Λ 1 are correlated, since they are both determined by the parameter ω 1 .
Spin-2 Next, we consider the mixing between a massive spin-2 field and the Goldstone boson.
• At linear order in σ µν , the mixing Lagrangian is L (2) πσ =ω 3 1 δg 00 σ 00 +ω 3 2 (δg 00 ) 2 σ 00 +ω 2 3 δK µν σ µν +ω 2 4 δg 00 δK µν σ µν , (3.19) where it was necessary to include higher-derivative operators to get the relevant interactions for the spatial components σ ij . In the decoupling limit, the mixing with the Goldstone boson is (3.20) We will focus on the traceless part of σ ij , which we denote byσ ij . Only the cubic mixinġ π∂ i ∂ j πσ ij will lead to the characteristic angular structure in the bispectrum. Since the quadratic mixing does not affect the angular structure, we will simply choose ∂ i ∂ j πσ ij as a representative example. Unlike the spin-1 case, the sizes of the quadratic and cubic mixing operators are controlled by two independent parameters,ω 3 andω 4 .
We will study the following mixing Lagrangian where we defined This is similar to the spin-1 mixing Lagrangian (3.17), except that the quadratic and cubic mixing parameters, ρ 2 and Λ 2 , are now independent.
Spin-s Performing the same analysis for a field with arbitrary spin s > 2, we find the following mixing Lagrangian where ∂ i 1 ···is ≡ ∂ i 1 · · · ∂ is . As in the case of spin-2, these interactions generically arise from independent operators, i.e. ρ s , Λ s , and λ s are independent parameters.
The mixing in (3.25) can convert hidden non-Gaussianity in the σ-sector into visible non-Gaussianity in the π-sector. To allow for this possibility, we add cubic self-interactions to the action for σ, which schematically we can write as a 3s L (s) with suitable symmetric contractions of spatial indices.

Couplings to the Graviton
We will also be interested in the couplings between massive spinning fields and the graviton, γ ij . For simplicity, we will only consider linear couplings to γ ij , but the generalization to higher orders will essentially be straightforward.

Spin-1
The leading couplings to the graviton arise from Note that, in our perturbative treatment, the on-shell conditions for σ µ hold at the background level, so thatḡ µν ∇ µ σ ν = 0. The first term in (3.27) is therefore proportional to g µν ∇ µ σ ν = δg µν ∇ µ σ ν and starts at cubic order in fluctuations. In terms of π and γ ij , the mixing Lagrangian becomes where τ 1 ≡ 4ω 2 5 /f 2 π and γ c ij ≡ 1 2 M pl γ ij denotes the canonically normalized graviton. The first term in (3.28) is higher order in derivatives than the operator γ ij ∂ i πσ j . However, the latter only arises from the tadpole σ 0 , and is therefore required to have a vanishing coefficient. Moreover, a quadratic mixing between the spin-1 field and the graviton is forbidden by kinematics: any such mixing will involve spatial gradients and hence must vanish because the graviton is transverse, Spin-2 The couplings between a massive spin-2 field and the graviton follow from Note that we have already encountered the operator δK µν σ µν in (3.19). In our perturbative treatment, the on-shell traceless condition holds at the background level,ḡ µν σ µν = 0, which implies that g µν σ µν = δg µν σ µν , so that the second term in (3.29) starts at cubic order in fluctuations. The cubic operator δg 00 δg µν ∇ µ σ ν0 will not be considered, since its effects are indistinguishable from those of the first term in (3.28). Introducing π and γ ij , the mixing Lagrangian becomes whereρ 2 ≡ −ρ 2 f 2 π and τ 2 ≡ 4ω 3 7 /f 2 π . Note that we have only kept the spatial components in the coupling to the mass term. Unlike in the spin-1 case, there is now a quadratic mixing between the spin-2 field and the graviton, whose size is correlated with the π-σ mixing in (3.23). The other possible form of mixing γ ij σ ij comes from the tadpoleσ and is thus absent.
Spin-s For arbitrary spin s > 2, the leading interactions with the graviton and the Goldstone take the following form whereρ s ≡ −ρ s f 2 π . Again, we have only kept interactions that involve the purely spatial components of the field. In practice, there are other low-dimensional operators that can also contribute to the correlator γζζ with the same angular structure, such asγ ij σ ij0···0 .

Bounds on Mixing Coefficients
It is important to determine how large the mixing interactions of the previous section can be made while keeping the effective theory under theoretical control. In this section, we will discuss bounds arising from i) the requirement that the mixing interactions can be treated perturbatively and ii) the absence of superluminal propagation. 9 Finally, we also consider what range of coefficients yields a technically natural effective field theory, in the sense of stability under radiative corrections. In Section 4, we will consider the implications of these constraints on the size of non-Gaussianities.

Perturbativity
We wish to treat the mixing interactions as perturbative corrections to the free-field actions for the Goldstone boson and the massive spinning fields. Since massive particles decay outside the horizon and oscillate rapidly inside the horizon, the dominant contributions to correlation functions will occur at horizon crossing of the Goldstone boson, corresponding to frequencies of order H. Consistency of the perturbative description therefore requires that the sizes of the mixing interactions at ω ∼ H are smaller than the terms in the free-field actions. This puts constraints on the couplings in the mixing Lagrangians discussed in the previous section. For c π = 1, the criteria for a consistent perturbative treatment require little more than dimensional analysis. The dimensionful couplings of relevant interactions have to be less than H, while those of irrelevant interactions have to be greater than H. The dimensionless couplings of marginal interactions have to be less than unity. For example, for the couplings appearing in (3.23), we require {ρ 2 , λ 2 } < 1 and Λ 2 > H. Similar considerations apply for the couplings in (3.17) and (3.25). For c π = 1, determining the perturbativity constraints on the mixing parameters requires a more careful analysis. Spatial gradients of the Goldstone mode are enhanced and the correlation functions can receive contributions from a second time scale, the time of crossing of the sound horizon. We will return to this complication in Section 4.

Superluminality
The breaking of time diffeomorphism invariance can modify the actions for spinning fields of Section 2 by introducing additional non-Lorentz-invariant interactions. For concreteness, we will confine our discussion in this subsection to the case of spin one, but we expect similar results to hold for higher spins. In unitary gauge, the most general quadratic action for a spin-1 field is where F µν ≡ ∂ µ σ ν − ∂ ν σ µ , and the structure of the kinetic part is enforced by gauge invariance in the massless limit. The departure from the Lorentz-invariant action is characterized by the parameters a 0 and a 1 , which lead to nontrivial sound speeds for the longitudinal mode, c 0 ≡ 1/ √ 1 + a 0 , and for the transverse mode, c 1 ≡ 1/ √ 1 + a 1 . To see this, we consider the on-shell equations of motion in the flat-space limit, where σ t i denotes the transverse mode, ∂ i σ t i = 0. The components of the spin-1 field propagate subluminally with no gradient instability as long as {a 0 , a 1 } ≥ 0. A tachyonic instability is avoided for m 2 > 0.
Even if the spin-1 field propagates subluminally, the mixing with the Goldstone boson π can lead to superluminal propagation in the coupled system. Requiring the absence of superluminality imposes a constraint on the size of the quadratic mixing term in (3.17). To derive this constraint, we consider the on-shell equations of motion for the coupled system, We see that the transverse mode does not mix with π, and hence its dispersion relation is unmodified. After diagonalizing the coupled π-σ system, the dispersion relations obeyed by the normal modes are where δ 2 ≡ c 3 π ρ 2 1 /m 2 . For large k, subluminality implies the following constraint Note that the mixing is required to vanish if either c π or c 0 are equal to 1. (A similar result for the mixing with a scalar field was found in [9].) However, even a relatively small deviation of c π and c 0 from 1 is sufficient to allow ρ 1 to be of order H (i.e. of order the maximal size allowed by perturbativity). For simplicity, we will therefore work with c 0 = c 1 ≈ 1, but incorporating nontrivial sound speeds for the spin-1 field could be done straightforwardly using modifications of the mode functions given in Appendix A. Similarly, we will assume that all spinning fields obey a relativistic dispersion relation.

Naturalness
Finally, we will consider constraints arising from the radiative stability of the effective theory, e.g. we require that the masses of spinning fields do not receive large loop corrections. This is more of a philosophic criterion rather than a strict consistency condition.
• Let us consider the interactionπ c ∂ i 1 ···is π c σ i 1 ···is in (3.25), suppressed by the scale Λ s . At one loop, this term generates the following correction to the non-Lorentz-invariant mass term, Naturalness of the mass of the spinning field requires δm 2 To estimate the size of (3.40), we take the cutoff of the π-loop to be of order the strong coupling scale of the Goldstone sector. For c π = 1, we can, in principle, extend the π-loop up to the symmetry breaking scale, i.e. Λ ∼ f π , while, for c π 1, the effective theory of the Goldstone becomes strongly coupled at Λ ∼ f π c π . We will therefore use Λ ∼ f π c π for all values of c π . The condition for radiative stability then becomes Typically, this constraint requires Λ s to be slightly larger than f π .
• Next, we consider the interaction λ sπc σ 2 i 1 ···is in (3.25). This leads to the following radiative correction to the non-Lorentz-invariant mass term, Cutting off the loop at Λ ∼ f π c π , we obtain the following constraint for radiative stability: (3.43) The interaction λ sπc σ 2 i 1 ···is can also give a correction to the kinetic term for the Goldstone. However, on dimensional grounds, it is easy to see that this interaction only contributes a negligible correction to the sound speed of π.
• Lastly, the radiative correction generated by the cubic self-interaction of the spinning fields is where the couplings ξ s are of dimensions zero and one for odd and even spins, respectively, cf. (3.26). For even spins, we only get a fixed finite correction to the mass term. Since we require ξ s < H for perturbative control, the loop contribution from this interaction is guaranteed to be small. For odd spins, it is natural to take the cutoff for the σ-loop to be the strong coupling scale Λ s . We then get where we have used the naturalness constraint (3.41) on Λ s in the second inequality.

Cosmological Correlators
We will now compute the effects of massive particles with spin on the correlation functions of the Goldstone boson and the graviton. Following [6], we will study separately the contributions from local and non-local processes. Local processes are, by definition, those whose imprint can be mimicked by adding a local operator in the low-energy effective theory of the light fields alone. Non-local processes, on the other hand, capture particle production effects which cannot be mimicked by additional local operators. While the latter are the distinctive signature of extra particles during inflation, the amplitude of such effects is exponentially suppressed for masses above the Hubble scale. We will discover that the sound speed of the Goldstone boson plays a crucial role in controlling the relative size of the local and non-local processes.

ζζ
Before discussing a potentially richer structure in the bispectra, we will gain some useful insights by first examining the effect of massive particles on the power spectrum ζζ (see Fig. 2). We will separate the correlation function into distinct contributions coming from local and non-local processes. Spin doesn't play a big role in the correction to the power spectrum, so for simplicity we will consider a minimally-coupled massive scalar field σ, whose two-point function in de Sitter space is where iµ is the Hankel function of the first kind and µ ≡ m 2 /H 2 − 9/4. We will focus on massive particles belonging to the principal series, so that µ ≥ 0. The local part of the two-point function has support only at coincident points in position space, while the non-local part describes correlations over long distances. In Fourier space, the local and non-local parts of the two-point function are analytic and non-analytic in the momentum k, respectively. In the late-time limit, we have Away from the late-time limit, we use a series expansion of the Hankel function, to decompose the two-point function (4.1) into its local and non-local pieces. Summing over the set of local and non-local contributions, the two-point function can be split into 10 where J iµ denotes the Bessel function of the first kind.
To illustrate the distinct roles played by local and non-local parts, let us consider a coupling between π and σ of the form d 4 x a 3 ρπ c σ [37,38]. At tree level, this produces the following correction to the power spectrum of ζ: diagram corresponding to c π = 1 (c π < µ −1 ). The Hubble radius is denoted by r H ≡ H −1 . We see that for c π < µ −1 the horizon crossing of the Goldstone boson occurs before the turning point of the massive particles, while for c π = 1 it occurs after. where The integral in (4.8) can be evaluated analytically to give where 2 F 1 is the hypergeometric function. It is instructive to consider the c π → 1 and c π → 0 limits of the result (4.10): • For c π = 1, the hypergeometric function becomes unity, and (4.10) scales as e −2πµ for large µ, as expected for the pair-production of massive particles.
• In the limit c π → 0, we instead get which scales as e −πµ for large µ instead of the usual Boltzmann factor e −2πµ .
To see why the exponential suppression of C 1 changes for c π 1, we need to consider the change in the dynamics of σ and π. There are two relevant timescales in the problem: i) at the turning point, |kη| ∼ µ, the mode function of the massive particle starts to decay, ii) at the sound horizon crossing, |kη| ∼ c −1 π , the Goldstone boson freezes. For c π = 1, event i) occurs before ii), while for c π < µ −1 , the order is reversed (see Fig. 3). As a consequence, the integral in (4.8) is dominated at the horizon crossing of π for c π = 1, while it is 0.01 0.1 1 10 100 dominated by the turning point of σ for c π < µ −1 . This is illustrated in Fig. 4, where we show the Wick-rotated integrand of the integral in (4.8) as a function of x = |kη|. A notable feature is the peak at x ∼ µ, which increases for small c π . For c π = 1, the turning point occurs before horizon crossing and the overlap between π and σ is suppressed. For c π < µ −1 , on the other hand, the turning point occurs after the freeze-out of the Goldstone, which enhances the feature at x ∼ µ. This qualitatively explains the boost in the amplitude of C 1 for small c π .
Let us now consider the time-ordered integral C 2 in (4.9). For general c π , it cannot be evaluated analytically, but some insights can be obtained by taking the limits c π → 1 and c π → 0: • For c π = 1, the above decomposition of the σ-propagator into local and non-local pieces leads to [37] where ψ (1) (z) = ∂ 2 z ln Γ(z) is the polygamma function of order 1. 11 For large µ, the first term in (4.12) scales as µ −2 , which has a simple interpretation: a heavy field contributes to non-renormalizable interactions in the low-energy effective theory of the light fields with coefficients given by inverse powers of the mass of the heavy field. The second term is instead suppressed by e −2πµ , describing an effect which cannot be captured by a local Lagrangian of the light fields alone. Finally, we see that the non-local part of the σ-propagator does not contribute to the correction to the power spectrum.
• In the limit c π → 0, we find We wish to highlight several features of this result. First, the local contribution to C 2 vanishes. This follows from the simple fact that the Goldstone bosons become non-propagating when c π = 0; hence, they can only communicate to each other through non-local effects. Second, the non-local contributions to C 1 and C 2 precisely cancel each other, implying that the correction to the two-point function (4.7) vanishes faster than c 2 π in the limit c π → 0. This is the result of the cancellation between the contributions from the forward and backward branches of the integration contour. A way to see this is to drop the exponentials in c π in (4.8) and (4.9), and notice that C 1 + C 2 is now proportional to the sum of all Schwinger-Keldysh propagators for the σ field; these propagators add up to zero. Of course, for small (but finite) c π , we do not expect this cancellation to be exact. To understand how the result for general c π interpolates between these two limiting behaviours, we evaluate C 2 numerically. Figure 5 shows the analytical result (4.10) for C 1 and a numerical computation of C 2 , both as functions of c π . As c π is lowered, the exponential dependence on µ for both of the integrals changes. For C 1 , this happens relatively quickly when c π µ −1 , agreeing with the intuition that reversing the ordering of the turning point of σ and the horizon exit of π changes the solution qualitatively. On the other hand, the transition in the exponential behavior for C 2 only occurs for very small c π , typically much smaller than the lower limit required for perturbative control of the non-renormalizable interactionπ(∂ i π) 2 associated with c π . This implies that, while for C 2 the dependence on µ > 1 will not change much within the allowed range of c π > 10 −2 , the exponential suppression e −2πµ of C 1 can be reduced to e −πµ when c π < µ −1 .

ζζζ
Next, we consider the imprints of massive spinning particles on the three-point function ζζζ . In single-field inflation, a long-wavelength curvature perturbation locally corresponds to a rescaling of the background experienced by short-wavelength fluctuations. As a result, the bispectrum ζζζ satisfies a consistency relation for the squeezed limit [11,12,39]. In particular, we can write a Taylor expansion around the squeezed limit, where the leading coefficient is determined by the tilt of the scalar power spectrum, b 0 = −(n s −1). The consistency condition furthermore fixes the coefficient of the linear term, b 1 , and partially constrains higher-order coefficients [40][41][42][43][44]. Since the contributions coming from b 0 and b 1 cannot be measured by a local observer [45,46], any physical effect will only appear at order (k 1 /k 3 ) 2 [47]. A crucial consequence of the consistency relation is the existence of the Taylor expansion (4. 16) with only integer powers of k 1 /k 3 . Interesting non-analytic deviations from (4.16), however, are known to arise in the presence of additional fields. For example, fractional powers (k 1 /k 3 ) ν can be present in quasi-single-field inflation [8], with scaling 0 < ν ≤ 3/2 in between the fully constrained (k 1 /k 3 ) 0 term and the physical (k 1 /k 3 ) 2 term. In this section, we will study such deviations for additional fields that carry spin. Figure 6 shows all possible tree-level contributions to ζζζ . The three diagrams share many qualitative features, so to avoid repetition we will mostly concentrate on the analysis of the single-exchange diagram [(a)], and only highlight the differences that arise for the other two diagrams [(b,c)]. We will split the contributions to the bispectrum into its local and non-local parts. To avoid confusion with the alternative usage of "local non-Gaussianity", we will refer to these contributions as analytic and non-analytic, respectively. (This terminology highlights the distinctive scaling behavior in the squeezed limit.) Although we will ultimately be interested in the behavior of the latter, the observability of the signal will depend on the full bispectrum, so we will present the results for both types of contributions. As before, we will mostly restrict our analysis to particles in the principal series, with µ s ≥ 0. Single-exchange diagram We will first compute the bispectrum associated with the exchange of a single spinning field (Fig. 6a). The relevant interaction Lagrangian is [cf. eq. (3.25)] We obtain the following bispectrum where an integral representation of the function I (s) is given in Appendix B. The dimensionless parameters α s are where we have included powers of c π in α s , so that the function I (s) does not scale parametrically with c π . By this we mean that I (s) saturates to a constant value in the limit of small c π , similar to the behavior of the integrals (4.8) and (4.9). The requirement of a perturbative treatment of non-Gaussianity implies that α s < 1 . (4.20) Notice that we have a stronger perturbativity condition on the bare parameters ρ s and Λ s for subluminal c π , which takes into account the fact that the dispersion relation, ω = c π k, is nonrelativistic.
Size of NG.-It is customary to quantify the size of non-Gaussianity by the parameter where the bispectrum is evaluated in the equilateral configuration, k 1 = k 2 = k 3 ≡ k. The overall size of the non-Gaussianity can only partially be read off from the prefactor in (4.18), since there is a hidden dependence on µ s in the function I (s) . An estimate for the size of non-Gaussian signal is where f (µ s ) gives the appropriate mass suppressions for the analytic and non-analytic parts 12, 13 analytic, e −πµs non-analytic, c π = 1, e −πµs/2 non-analytic, c π < µ −1 s .

(4.23)
We see that there are two sources of suppression in the signal: the mass suppression as a function of µ s and the mixing efficiency parameterized by α s . At the same time, there is a ∆ −1 ζ ≈ 10 5 enhancement in the signal. It is this large factor that can, in principle, allow for observable non-Gaussianity even in the presence of the above suppressions. The size of the analytic part is only power-law suppressed and thus dominates for large mass, whereas the non-analytic part is always accompanied by an exponential Boltzmann suppression. For c π = 1, the dominant non-analytic term is suppressed by e −πµs . As explained in [6,13], this arises from the quantum interference of two wavefunctions: Ψ[2σ] ∝ e −πµs for pair-produced massive particles and Ψ[0σ] for the wavefunction involving no spontaneously created massive particles. This interference contribution is larger than the probability of pair-producing massive particles, which is |Ψ[2σ]| 2 ∝ e −2πµs . For c π < µ −1 s 1, the exponential suppression of the non-analytic part changes to e −πµs/2 . We have already encountered this phenomenon in §4.1: for small c π the horizon crossing of the Goldstone boson occurs before the turning point in the mode function of the massive particle. In this case, we are picking out the contribution of the wavefunction for a pair of massive particles not in the late-time limit, but at the turning point, which comes with a different exponential factor.
In §3.4.3, we derived naturalness constraints on the mixing parameters of the effective theory. For the parameter α s in (4.19), the radiative stability of the mass (3.41) implies

(4.24)
For c π = 1, this naturalness constraint is rather strong, implying that large non-Gaussianity, f NL > 1, is only possible if additional physics, such as supersymmetry, stabilizes the mass of the spinning particle, or if the mass term is fine-tuned. For c π = 1, the current observational constraint c π ≥ 0.024 [29] still allows for naturally large non-Gaussianity, although within a rather narrow range in the small c π regime.
Some comments are in order concerning the observability of particles with odd spins. In [6], it was shown that the diagram due to the exchange of an odd-spin particle vanish exactly at leading order in the weak breaking of conformal symmetry. At subleading orders, however, there are non-zero contributions from odd-spin particles. 14 When conformal symmetry is strongly broken, these terms become as important as the leading ones, and odd-spin particles can leave 14 When the approximate conformal invariance is valid, we can think of this in terms of correlation functions of the inflaton Φ(t, x) = φ(t) + ϕ(t, x), whereφ = 0 characterizes the weak breaking of conformal symmetry. The leading three-point function for the inflaton perturbation ϕ will be given by the four-point function of Φ with one external leg set toφ: where · · · inf denotes an inflationary correlation function which breaks conformal symmetry [48]. However, in the conformally symmetric case, ϕϕσ vanishes when σ has odd spin [49]. The next-to-leading order result is given by the six-point function with three insertions ofφ, This is suppressed by an additional factor ofφ 2 , but notice that the correlator ϕϕσ inf , not being constrained by conformal symmetry, does not have to vanish for odd-spin σ.
an equally relevant imprint on the correlation function ζζζ . Nevertheless, the amplitude of the bispectrum with an intermediate spin-1 particle is As long as the mixing is perturbative, ρ 1 < H, this non-Gaussianity is constrained to be less than unity. We see that a spin-1 particle cannot lead to large non-Gaussianity because the size of the cubic vertex in (3.17) is tied to the quadratic mixing coefficient. In fact, the same reasoning applies to the coupling to scalar fields, which is why the single-exchange diagram has been neglected in the context of quasi-single-field inflation [8,9]. This fact, however, is only tied to spins zero and one, and the bispectrum does not have to be suppressed for higher odd-spin particles. Moreover, we will see that the diagrams involving more than a single exchange can allow for observable non-Gaussianity, even for spin one.
Shape of NG.-Before considering the general shape of the bispectrum, we will first analyze the singular behavior of the bispectrum in the squeezed limit, mainly concentrating on particles with even spins. We will quote results whose derivations can be found in Appendix C.
• For the analytic part of the bispectrum, we get (4.28) We see that the local effects of massive particles lead to the same squeezed limit behavior as for single-field inflation, cf. (4.16). This is expected, since the massive particle can be integrated out for large µ s , producing an effective cubic vertex of the formπ(∂ i 1 ···is π) 2 . The presence of extra particles therefore cannot be inferred from this part of the signal. Although the analytic part of the non-Gaussianity is itself interesting and more information can be gained by analyzing its shape for general momentum configurations, we have to treat it as an effective noise in the squeezed limit as far as the detection of extra particles is concerned.
• For the non-analytic part, we find where the phase φ s is uniquely fixed in terms of µ s and c π (see Appendix C). The suppression factor (k 1 /k 3 ) 3/2 represents the dilution of the physical particle number density due to the volume expansion. This non-analytic scaling in the squeezed limit, corresponding to an intrinsically non-local process, cannot be mimicked by a local interaction within the effective theory of a single field. The signal contains oscillations in ln(k 1 /k 3 ), with a frequency set by the mass of the spinning particle. This is due to the fact that the wavefunctions of massive particles oscillate logarithmically in time on superhorizon scales. The spin of the extra particle is reflected in the angular dependence, which is given by a Legendre polynomial of the angle between the short and long momenta.
The above behavior applies for particles in the principal series, for which µ s ≥ 0. For particles in the complementary series, µ s becomes imaginary and the scaling of the squeezed bispectrum changes to with ν s ≡ −iµ s real. For s ≥ 2, unitarity implies ν s ∈ [0, 1/2), and the singular behavior in the squeezed limit is suppressed by at least k 1 /k 3 compared to the leading term in the consistency relation (4.16).
The fact that the polarization tensors corresponding to odd-spin particles are odd under the exchange of two short momenta, together with momentum conservation, implies that the signal will gain an extra suppression factor of k 1 /k 3 in the squeezed limit compared to the case of even spin. This means that the non-analytic part due to odd-spin particles scales as (k 1 /k 3 ) 5/2 in the squeezed limit, which is more suppressed than the analytic part that scales as (k 1 /k 3 ) 2 . The latter, however, have an analytic dependence on momenta and correspond to local correlations in position space. Thus, the presence of odd-spin particles could still be inferred from long-distance correlations, although it might be subdominant compared to other non-local effects.
It is possible to understand the different behaviors in the squeezed limit intuitively. For concreteness, let us consider the exchange of a spin-2 field involving the interactions ∂ i ∂ j πσ ij anḋ π∂ i ∂ j πσ ij . The bispectrum in the isosceles-triangle configuration, k 2 = k 3 , consists of three different permutations of the external legs: , (4.31) where π n ≡ π(k n ), σ n ≡ σ ij (k n ) and I(k 1 , k 2 , k 3 ) ∝ P 2 (k 1 ·k 3 ) I (2) (µ 2 , c π , k 1 , k 2 , k 3 ). The nonanalytic squeezed limit (4.29) arises if the massive exchange particle carries the soft momentum, corresponding to the contribution I 1 in (4.31). This describes a non-local conversion process between the massive particle and the Goldstone boson between the horizon crossing times of the long and short modes. However, when the mass of the extra particle becomes large, it can be integrated out and the same effect will be captured by a local vertex. In that case, the bispectrum should become indistinguishable from that produced by a self-interaction of π, namelyπ(∂ ij π) 2 . Note, in particular, that this interaction is symmetric under the exchange of the momenta associated with the two external legs with spatial gradients. This allows us to gauge how well the interaction is approximated by a local vertex by looking at how similar the terms I 1 and I 2 are. Both I 2 and I 3 will lead to analytic scalings in the squeezed limit, where the latter produces (4.28).
To analyze the shape of the bispectrum for general momentum configurations, we proceed numerically. For this purpose, it is convenient to define a dimensionless shape function (4.32) Figure 7 shows two-dimensional projections of the shape function for spin 2 with µ 2 = 3, 5, 7 and c π = 1, 0.1 in the isosceles-triangle configuration, k 2 = k 3 . For the reasons explained in the previous paragraph, in Fig. 7 we have shown separately the shape functions corresponding to the contributions I 1 and I 2 in (4.31). 15 As anticipated, these contributions exhibit different scalings in the squeezed limit. The plots show (k 3 /k 1 ) × S, so that the analytic part is expected to approach a constant in the squeezed limit, while the non-analytic part grows as (k 1 /k 3 ) −1/2 for small k 1 . We see that the shape of the bispectrum is mostly governed by the non-analytic part for small mass, giving almost pure oscillations. The amplitude of this effect, however, goes as e −πµ 2 for large µ 2 . The analytic part, being power-law suppressed, therefore takes over in size as the mass increases, and the shape approaches the equilateral form in the limit of large mass. For large mass, it is clear that the non-Gaussianity is dominated by the analytic piece, with small oscillations coming from the non-analytic piece indicating the presence of a heavy mode. For c π = 1, the contributions I 1 and I 2 lead to the same shape of the bispectrum for µ 2 = 7, indicating that the π-σ conversion process has become local. Indeed, in this case the bispectrum precisely overlaps with that of the local interactionπ(∂ ij π) 2 . For small c π , we have argued that the exponential suppression is instead e −πµ 2 /2 . The fact that we see more pronounced oscillations for c π = 0.1 is a consequence of this. Moreover, for small c π , the shapes of the contributions I 1 and I 2 are no longer identical. Note that, in order for the massive particle to be integrated out, the time of its turning point should be much earlier than the time at which the Goldstone boson crosses its sound horizon, which translates into the condition c π > µ −1 2 . For c π = 0.1, this condition is not satisfied for the list of mass parameters used in the figure, which is the reason why we do not see the convergence to the local behavior. We have checked that the convergence does indeed happen for sufficiently large µ 2 > c −1 π . Another characteristic of the signal due to spinning particles is its angular dependence. Figure 8 shows the shape function of the total signal as a function of the angle between the long and short momenta, θ ≡ cos −1 (k 1 ·k 3 ), for a range of momentum configurations with fixed k 1 /k 3 . For visualization purposes, the plot has been rescaled so that it can be compared more easily to the Legendre polynomial P 2 (cos θ). As expected, the angular dependence converges to the pure Legendre behavior as the triangle becomes squeezed, k 1 /k 3 1. The non-zero offset is due to the analytic part which doesn't carry any angular dependence. We also see that the angular dependence deviates from the pure Legendre behavior as the triangle approaches the equilateral shape. Still, the peak around the flat triangle (θ = 180 • ) remains prominent regardless of the momentum configuration. This suggests that the information about a particle's spin can still be inferred without necessarily going to very squeezed momentum configurations, since the width of the peak is still fixed by the polarization tensor of the spinning particle. This property can  Figure 7: Shape functions (in units of α 2 ∆ −1 ζ ) for the spin-2 single-exchange diagram in the isoscelestriangle configuration, k 2 = k 3 , with µ 2 = 3 (top), µ 2 = 5 (middle), and µ 2 = 7 (bottom) for c π = 1 (left) and c π = 0.1 (right). The solid and dashed lines correspond to the numerical results for the parts of the signal corresponding to the terms I 1 and I 2 in (4.31), respectively. Not shown in the figure is the term I 3 , which produces an analytic scaling in the squeezed limit and is needed to obtain the full bispectrum. Convergence of the solid and dashed lines indicates that the same effect can be captured by a local verteẋ π(∂ ij π) 2 in the single-field EFT. The dotted lines show the analytical predictions for the non-analytic part. serve as an important tool for detecting odd-spin particles, whose signal in the squeezed limit necessarily gains an extra suppression in the soft momentum.

Double-exchange diagram
The bispectrum for the double-exchange diagram (Fig. 6b) is where the function J (s) is given explicitly in Appendix B, and the dimensionless parametersα s areα with λ s and ρ s defined in (3.25).
The size of the non-Gaussianity associated with the double-exchange diagram can be read off from (4.22) after replacing α s byα s , but with an extra suppression of µ −2 s , because this diagram involves another particle exchange. The condition for radiative stability (3.43) imposes the following upper limit on the size of the mixing parameter: Notice that this is a much weaker constraint than the corresponding constraint for the singleexchange diagram (4.24). Depending on the values of c π , this may or may not be stronger than the requirement for perturbativity, f NL < ∆ −1 ζ . This diagram can thus naturally produce detectable levels of non-Gaussianity, even for c π = 1.
Note that this diagram involves two π-σ conversion processes. When one of these processes becomes local, the double-exchange diagram becomes essentially equivalent to the single-exchange diagram. This can be seen by replacing one of theσ i 1 ···is legs in the cubic vertexπσ 2 i 1 ···is by ∂ i 1 ···is π, after which the interaction becomes the same as the cubic vertex for the single-exchange diagram. As a result, the squeezed-limit behavior for this diagram is essentially the same as that of the single-exchange diagram. Hence, the analysis we have presented for the single-exchange diagram applies also to the double-exchange diagram.
Triple-exchange diagram As indicated in (3.25), there is a slight difference between the form of the cubic self-interaction of spinning fields for even and odd spins. For concreteness, we will present the results for the former. The bispectrum for the triple-exchange diagram (Fig. 6c) is where P (k 1 ,k 2 ,k 3 ) ≡ ε 0 (k 1 )·ε 0 (k 2 )·ε 0 (k 3 ) is a symmetric contraction of the longitudinal polarization tensors ε 0 i 1 ···is (see Appendix A for the precise definition of the polarization tensor) that reduces to P s (k 1 ·k 3 ) in the squeezed limit. The couplingsα s arê where ξ s was introduced in (3.25). The function K (s) can be found in Appendix B.
The size of the non-Gaussianity associated with this diagram can, again, be read off from (4.22), with α s replaced byα s , and taking into account an extra suppression of µ −4 s . Although the qualitative features of the non-analytic signal will be similar to that of the other diagrams, there are some relevant differences. First, as shown in §3.4.3, naturalness does not constrain the size of the coupling ξ s , so the triple-exchange diagram allows for a naturally large non-Gaussianity. This is to be contrasted especially with the single-exchange diagram, where the naturalness criterion imposed a strong constraint on the size of the corresponding non-Gaussianity. Second, when the mass of the particle becomes large, the bispectrum is well-captured by a local vertex, namely (∂ i 1 ···is π) 3 with symmetric contraction of indices. Notice that, due to the number of spatial gradients, for s > 2 the squeezed-limit bispectrum is suppressed by more than (k 1 /k 3 ) 2 for small k 1 . This makes the non-analytic part, scaling as (k 1 /k 3 ) 3/2 , a rather clean signal in the squeezed limit.
Summary All diagrams in Fig. 6, except for the single-exchange diagram for spin one, can yield sizable non-Gaussianities within the perturbative regime. In order for this to be natural, the single-exchange diagram requires new physics or fine-tuning to stabilize the mass of the spinning particle, whereas both the double-and triple-exchange diagrams can naturally produce large non-Gaussianities. The non-analytic part of the bispectrum is suppressed by e −πµs for c π = 1, but only by e −πµs/2 for small c π . Typically, we find that f NL O(1) from the non-analytic part is possible if µ s 5 for c π = 1 and µ s 10 for c π 1.

γζζ
Lastly, we consider the tensor-scalar-scalar correlation function γζζ . In single-field inflation, a long-wavelength tensor fluctuation is locally equivalent to a spatially anisotropic coordinate transformation. Again, we can Taylor expand the expectation value around the squeezed limit, thus obtaining where γ λ , with λ = ±2, denotes the positive or negative helicity components of the graviton. As in the case of the scalar bispectrum, the leading coefficients are determined by the single-field consistency relation [11] (see also [40,42]). In particular, d 0 in (4.38) is given by When the consistency relation holds, it also completely fixes the linear term d 1 in (4.38), and physical effects appear at order (k 1 /k 3 ) 2 . The presence of new particles during inflation invalidates the Taylor expansion and leads to nonanalytic scalings in (4.38). Our goal in this section is to study these characteristic signatures of massive spinning particles. Figure 9: Tree-level diagrams contributing to γζζ . The solid, dashed, and wavy lines represent the curvature perturbation ζ, a spinning field σ i1···is , and the graviton γ ij , respectively.
All tree-level diagrams contributing to γζζ are shown in Fig. 9. Not all of these diagrams can lead to a nontrivial deviation from the consistency relation. For the diagrams [(a-c)] the same symmetry that generates the tensor consistency relation enforces corrections to the power spectrum, so that the relation in (4.38) and (4.39) is preserved [50]. Only the diagrams [(d-f)], which involve a quadratic mixing between the graviton and the intermediate particle, can lead to such a deviation. These diagrams have the same structure as those in Fig. 6, except that one of the legs in the quadratic mixing is replaced by an external graviton, so that the exchanging particle must carry the same helicity as the graviton. In the following, we will present results for the diagrams [(d-f)], mostly focusing on the single-exchange diagram [(d)] to avoid repetition. The quadratic γ-σ mixing vanishes for spins 0 and 1, so only particles with s ≥ 2 will contribute.
Single-exchange diagram We first consider the single-exchange diagram (Fig. 9d). The relevant interaction Lagrangian is [c.f. eqs. (3.25) and (3.31)] Using (3.7) and (3.11), we can write the coefficient of the quadratic mixing term as −ρ s r/8H. The perturbativity condition on the π-σ mixing, ρ s < 1, implies that the γ-σ mixing carries an extra suppression factor of r/8. The bispectrum corresponding to the single-exchange diagram is whereP λ s ≡ (1 − x 2 ) −λ/2 P λ s , with P λ s the associated Legendre polynomial. The function B (s) is given explicitly in Appendix B.
Size of NG.-We quantify the size of the tensor-scalar-scalar bispectrum by where the bispectrum is evaluated in the equilateral configuration, k 1 = k 2 = k 3 ≡ k, with vectors maximally aligned with the polarization tensor. This choice of normalization agrees with that adopted in [51] and implies f γζζ NL = √ r/16 for single-field slow-roll inflation [11]. An estimate of the size of the non-Gaussianity from the single-exchange diagram is where g(µ s ) denotes the appropriate mass suppressions for the analytic and non-analytic parts, which in the large µ s limit scale as 16 analytic, e −πµs non-analytic. (4.44) The enhancement of f γζζ NL by the large factor ∆ −1 ζ means that, in principle, the signal could be significantly larger than the one predicted from single-field slow-roll inflation, f γζζ NL √ r/16, even in the perturbative regime. As in the scalar case, the condition for radiative stability gives a rather strong constraint on the naturally allowed size of the bispectrum associated with the single-exchange diagram [cf. (4.24)]. While the size of the single-exchange diagram is strongly constrained by naturalness, both the diagrams [(e,f)] can lead to naturally large non-Gaussianity, as in the case of the scalar bispectrum. Future constraints on f γζζ NL from observations of the BT T correlator of CMB anisotropies were discussed in [51]. The proposed CMB Stage IV experiments [52] will have the sensitivity to reach σ( √ rf γζζ NL ) ∼ 0.1, which suggests that the tensor non-Gaussianity due to massive spinning particles would be detectable for r 10 −4 [g(µ s )α s ] −1 . 17 Shape of NG.-In the squeezed limit, γζζ behaves in the following ways: • The analytic part scales as Notice that the suppression of the analytic part in the squeezed limit increases with spin. This can be understood by looking at the form of the local vertex after integrating out the massive particle, which becomesπ∂ i 1 ···is π∂ i 3 ···isγi 1 i 2 . As we will see below, this means that the analytic part of the signal will be subdominant compared to its non-analytic counterpart in the soft graviton limit.
• For µ s ≥ 0, the squeezed limit of the non-analytic part of the bispectrum scales as where the phaseφ s is a function of µ s and c π (see Appendix C). Coupling to a particle with spin greater than two induces an extra angular structure. For imaginary µ s , we instead have with ν s ≡ −iµ s ∈ [0, 1/2). This gives a non-analytic (k 1 /k 3 ) 3/2−νs correction to the leading term of the consistency relation (4.38). Since unitarity implies ν s < 1/2, the squeezed-limit bispectrum due to massive spinning particles will be suppressed by at least k 1 /k 3 compared to the leading term in the tensor consistency relation. 18 Other diagrams The extensions to the diagrams [(e,f)] are completely analogous to the scalar case. Similar to the scalar three-point function, these diagrams have the advantage that they are less constrained by naturalness considerations. 17 Producing a large tensor contribution while keeping the scalar contribution small may require some fine-tuned cancellation between interactions in the scalar sector. This is because the interaction vertices in (4.40) and (4.17) arise from the same operators in unitary gauge. Suppressing the effects of the interactions in (4.17) would require balancing them against additional interactions such asπσ0···0. 18 A deviation from the leading term of the consistency relation due to spinning particles can arise in a number of ways: First, the unitarity bound can be evaded if the de Sitter isometries are not fully respected in the quadratic action of the spinning field [53,54]. Another possibility involves partially massless fields with spin greater than two, since the late-time behavior of these fields does not obey the same restrictions as for the massive case. It would be interesting to explore these possibilities further.

Conclusions
In this paper, we have studied the imprints of massive particles with spin on cosmological correlators using the framework of the effective field theory of inflation [17]. This generalizes the work of Arkani-Hamed and Maldacena (AHM) [6] to cases where conformal symmetry is strongly broken. Let us summarize our results and contrast them with the conclusions of AHM: • In AHM's more conservative analysis, the overall size of non-Gaussianity was too small to be observable even in the most optimistic experimental scenarios. Our results are cautiously more optimistic. Within the regime of validity of the effective field theory, we can accommodate observable non-Gaussianity as long as the masses of the new particles aren't too far above the Hubble scale during inflation.
• The key spectroscopic features of massive particles with spin do not rely on conformal invariance and therefore continue to hold in our analysis. As explained in [6], the masses and spins of extra particles during inflation can be extracted by measuring the momentum dependence in the squeezed limit.
• Our systematic effective field theory treatment of massive spinning particles during inflation allows for a complete characterization of their effects on non-Gaussian cosmological correlators, including their imprints beyond the squeezed limit. We showed that the characteristic angular dependence resulting from the presence of particles with spin persists even for more general momentum configurations. Having access to the complete correlation functions will be valuable for future data analysis.
• We also studied the effects of an explicit breaking of special conformal symmetry by introducing a sound speed c π for the Goldstone fluctuations. We found that, for c π < µ −1 s , the exponential suppression in the production of the massive particles, e −πµs , is changed to e −πµs/2 . For a given mass, the size of non-Gaussianity is therefore enhanced (or less suppressed) for small c π .
• Finally, we showed that particles with spin greater than or equal to two lead to a signature in the squeezed limit of γζζ . This signal may be observable in the BT T correlator of CMB anisotropies [51]. Figure 10 is a schematic illustration of current and future constraints on (scale-invariant) primordial non-Gaussianities. We see that the perturbatively interesting regime spans about seven orders of magnitude in f NL . Of this regime, three orders of magnitude have been ruled out by current CMB observations, leaving a window of opportunity of about four orders of magnitude. Accessing these low levels of non-Gaussianity will be challenging. Even optimistic projections for future CMB observations won't reduce the constraints by more than an order of magnitude. Digging deeper will require new cosmological probes, such as observations of the largescale structure (LSS) of the universe [55] and the tomography of the 21cm transition of neutral hydrogen gas [56]. Our results, together with [6][7][8][9][10], will help to find optimal observational strategies for extracting the subtle imprints of extra particles during the inflationary era. 21cm? Figure 10: Schematic illustration of current and future constraints on (scale-invariant) primordial non-Gaussianity. The "gravitational floor" denotes the minimal level of non-Gaussianity created by purely gravitational interactions during inflation [11].
A More on Spin in de Sitter Space In this appendix, we will derive various mathematical results that have been used in this work. In §A.1, we obtain the mode functions for massive spinning fields in de Sitter space by solving their equations of motion. We then derive the formula for the two-point function in §A.2.
Preliminaries We will work with the components of the spinning field σ µ 1 ···µs projected onto spatial slices, i.e. σ i 1 ···inη···η . We will find it convenient to write these as where ε λ i 1 ···in is a suitably normalized polarization tensor (see insert below). The sub/superscripts on the mode functions σ λ n,s label three "quantum numbers": s is the spin (or the rank) of the spacetime tensor field, n is its "spatial" spin, and λ is the helicity component of the spatial spin.
Polarization tensors.-In this insert, we will derive explicit expressions for the polarization tensors of arbitrary spin and helicity. The longitudinal polarization tensors are functions ofk, while the transverse polarization tensors in addition depend on two polarization directionsε ± , withk ·ε ± = 0. Sinceε + andε − are related to each other by the reality conditionε + = (ε − ) * , let us denote one of them byε. The polarization tensors of helicity λ satisfy the following conditions: ii) traceless: ε λ iii3···is = 0.
iii) transverse:k i1 · · ·k in ε λ i1···is = 0, when n > s − λ. The last condition implies that the polarization tensor is of the form wherek i1 ε λ i1···i λ (ε) = 0 and f i1···i s−λ is some tensor. Let us contract with vectors q and define where we have defined x ≡ q 2 , y ≡ q ·k, and z ≡ q i1 · · · q i λ ε λ i1···i λ . The function F λ s is a homogeneous polynomial in q, so that The transverse and traceless conditions translate into where d is the number of spatial dimensions. Taking derivatives of (A.4) and (A.5), and substituting into (A.6), we get Without loss of generality, we now set x = q 2 ≡ 1. The solution to (A.5) and (A.7) is F λ s (y, z) ∝ zP β λ βs (y) , (A. 8) whereP β λ βs is part of the associated Legendre polynomial P β λ βs of degree β s ≡ 1 2 (2s + d − 3) and order β λ ≡ 1 2 (2λ + d − 3), defined by P β λ βs (y) = (1 − y 2 ) β λ /2P β λ βs . We will set P β λ s ≡ P |β λ | s and distinguish the opposite helicities only by the phase. For d = 3, this reduces to This result includes longitudinal polarization tensors for λ = 0 and z = 1. It is straightforward to obtain explicit expressions for the polarization tensors by stripping off the contractions with q in (A.9) and symmetrizing the indices: The self-contraction of the polarization tensors can be written as When choosing the orthogonal direction to be in, say, the z-direction, there will be in total of 2 s non-zero components for the polarized tensor ε s i1···is , which are ±1 or ±i up to a phase. This means that ε s i1···is ε s * i1···is = 2 s with some overall normalization which we set to unity for convenience.

A.1 Mode Functions
In this section, we will derive the de Sitter mode functions for fields with spin. We will explicitly derive the mode functions for fields with spins 1 and 2, and present the results for arbitrary spin at the end.

Spin-1
The equation of motion of a massive spin-1 field σ µ is with ∇ µ σ µ = 0 and m 2 1 = m 2 + 3H 2 . The components σ η and σ i then satisfy where a prime denotes a derivative with respect to conformal time, and To decouple equations (A.14) and (A.15), we expand the field σ µ into its different helicity components, We demand that the polarization vectors ε λ i (k) satisfŷ and the transverse condition (A.16) becomes The solutions to these equations with the Bunch-Davies initial condition are where A 1 ≡ e iπ/4 e −πµ 1 /2 and Z ±1 1 denotes the normalization constant for the helicity-±1 mode of the spin-1 field. We have also suppressed the argument −kη of the Hankel functions H iµ 1 ≡ H (1) iµ 1 for brevity.
A few comments are in order. First, note that for m = 0 equation (A.23) for the transverse mode becomes the flat space wave equation, whose solutions are simply plane waves. This is because the action of a massless spin-1 field is conformally invariant, so the mode in de Sitter space behaves as if it were in flat space. On the other hand, we do not see this behavior for the longitudinal mode. In particular, the longitudinal mode blows up relative to the transverse mode as we go to the infinite past η → −∞. We can understand this as follows. The mass term m 2 /H 2 η 2 in the action (2.1) is time dependent, so the spin-1 field is effectively massless in the infinite past, in which case the longitudinal mode turns into a pure gauge mode.
We still need to determine the normalization constants N 1 and Z ±1 1 . This is done by imposing orthonormality of mode functions under the inner product This orthonormality condition guarantees that we get the standard equal-time commutation relation upon canonical quantization. We have where W denotes the Wronskian. Substituting (A.25) and (A.27), we obtain Note that the time dependences in (A.30) and (A.31) cancel in (A.29). Imposing (A.28), we then get The normalization for the transverse mode can be determined in a similar way. We get Notice that the normalization for the longitudinal mode blows up when m = 0, which, again, does not signal any pathologies, since the longitudinal mode becomes a pure gauge mode in this limit.

Spin-2
The equations of motion and the constraints satisfied by a massive spin-2 field σ µν are In terms of components, these are As before, we expand the Fourier modes into helicity eigenstates Let us denote the traceless part of the spatial tensor byσ ij , so that σ ij =σ ij + 1 3 σ ηη δ ij , and decompose the mode functions into different helicities: Demanding that the polarization tensors satisfŷ and fixes ε ±2 ij up to a phase. Fork along the z-direction, this can be chosen to be The equations satisfied by the different helicity modes are subject to the transverse conditions The solutions with Bunch-Davies initial conditions are for the longitudinal modes, and for the higher-helicity modes.
To fix the normalization, we again impose orthonormality of the mode functions We have The condition (A.60) then sets the normalization constant to be We see that this diverges at m 2 = 0 and m 2 = 2H 2 . This is again to be expected. For m = 0, the action gains gauge invariance, in which case only the helicity-±2 modes are physical. For m 2 = 2H 2 , the field becomes partially massless, and the number of propagating degrees of freedom becomes four. In both cases, the longitudinal mode becomes a pure gauge mode. Finally, determining the normalizations of the transverse modes in an analogous way, we get In the massless limit, Z ±1 2 diverges and only Z ±2 2 remains finite.
Spin-s For spins higher than two, we need to solve the on-shell equations (2.7). In order to decouple these equations, we expand the field σ µ 1 ···µs into its different helicity components A mode of helicity λ and n polarization directions can be written as where σ λ n,s = 0 for n < |λ|. The helicity-λ mode function with n = |λ| number of polarization directions satisfies whose solution is given by The other mode functions can then be obtained iteratively using the following recursion relation: . (A.71) Having obtained the formula that enables us to compute the mode functions of arbitrary spin and helicity, let us now fix their normalization constants. In order to do so, we first define an inner product between two mode functions. Note that if f µ 1 ···µs and h ν 1 ···µs are two solutions to (2.7), then the current is conserved, ∇ µ J µ = 0. This means that we can define an inner product of two solutions where Σ denotes a spacelike hypersurface,ĝ is the determinant of the spatial metric, and n µ is the timelike unit vector orthogonal to Σ. The conservation of the current (A.72) implies that the inner product is time independent. For the FRW metric, the above inner product reduces to The normalization in (A.69) is then determined by imposing orthonormality under the inner product (A.74): Since the inner product is time independent, it does not matter which time slice we choose to evaluate the integral on. We will therefore evaluate the integral on the future boundary by taking the limit η → 0. From (A.70), we note that σ λ n 1 ,s is subleading compared to σ λ n 2 ,s in the limit η → 0 for all n 1 < n 2 , so we simply need to compute the Wronskian of the mode with the highest number of polarization directions, σ λ s,s . If we had kept all the Wronskians, then the subleading time-dependent terms would cancel. Note also that the trace terms in (A.70) become subleading in the limit η → 0, so we will drop these terms. Since (A.74) is a constant, the leading term in the Wronskian must scale as η 2(1−s) to cancel off the factor a 2(1−s) . In the insert below, we will show that the orthonormality condition fixes the normalization constant to be Note that the normalization constant has poles at µ 2 s = {−(n+ 1 2 ) 2 } s n=λ , at which the spinning field becomes (partially) massless and some of the helicity modes become unphysical. For convenience, we will denote the normalization of the longitudinal mode by N s ≡ Z 0 s (N s ≡ Z 0 s ).
Derivation of (A.76).-First, note that the n-th mode function can be cast in the form σ λ n,s = A s Z λ s (−kη) 3/2−n (x n + iy n )H iµs + (w n + iz n )kηH iµs+1 , (A.78) by use of the recursion relation H iµs+1 (x) + H iµs−1 (x) = (2iµ s /x)H iµs (x). The coefficients x n , y n , w n , and z n can in general depend on time, but are constant in the limit η → 0. The Wronskian is W σ λ n,s , σ λ * n,s = 4ik(Z λ s ) 2 π(kη) 2(n−1) X n − 2µ s Y n (coth πµ s − 1) , Let us show that in fact Y n = 0 for any n-th order mode function. We do this by induction. First, it is trivial to check that this is satisfied by the mode (A.69). Now, assume that Y n = 0 is satisfied at some n-th order. Using the recursion relation (A.70), and taking the limit η → 0, we get where 2x n+1 = −2µ s x n + (2n + 1)y n , 2y n+1 = −(2n + 1)x n − 2µ s y n , 2w n+1 = −2y n + 2µ s w n + (2n + 1)z n , 2z n+1 = 2x n − (2n + 1)w n + 2µ s z n . (A.82) These coefficients then give Hence, Y n+1 = 0. Since n was arbitrary, we conclude that Y n = 0 for all n. Next, we show that the Wronskian of the n-th longitudinal mode function has the form W σ λ n,s , σ λ * n,s = 4ik(Z λ s ) 2 π(kη) 2(n−1) The Wronskian of the mode function (A.69) is and hence satisfies (A.84). Assuming that (A.84) is true at n-th order and using (A.83), we get where in the last line we have use the fact that Thus, we have proven (A.84). Finally, the inner product (A.74) is given by Note that our final normalization depends on the normalization of the polarization tensors. This does not affect correlation functions, however, as we show in the next section. Plugging (A.12) into (A.88) and imposing (A.75), we obtain (A.76).

A.2 Two-Point Function
In this section, we will compute the two-point functions of spinning fields. For this purpose, it will be convenient to contract free indices of the spinning fields with auxiliary vectors. In other words, we will compute (n · σ) 2 s ≡ (n i 1 · · · n is σ i 1 ···is (η)) ñ j 1 · · ·ñ js σ j 1 ···js (η ) , (A.89) where the prime on the expectation value indicates the removal of the momentum-conserving delta function, and n ≡ (cos α, sin α, i) andñ ≡ (cos β, sin β, −i) are null vectors. For generic η and η , the two-point function is where χ ≡ α − β. In the late-time limit (or the long-wavelength limit), the two-point function simplifies considerably. We get This late-time expectation value matches the two-point function of a spin-s field of a conformal field theory living on the future boundary, which have been computed in [6]. (n i1 · · · n is ε λ i1···is )(ñ j1 · · ·ñ js ε λ * j1···js ) σ λ s,s σ λ * s,s . (A.93) Let us compute σ λ s,s σ λ * s,s in the late-time limit. First, recall that we can cast the mode function in the form σ λ n,s = A s Z λ s (−kη) 3/2−n (x n + iy n )H iµs + (w n + iz n )kηH iµs+1 . (A.94) Taking the asymptotic limits of the Hankel functions, we get σ λ n,s σ λ * n,s η,η →0 where W n ≡ x 2 n + y 2 n + 2µ s (x n + iy n )(iw n + z n ) . (A.96) Using (A.82), we obtain the recursion relation Following similar arguments as in the previous section, it can then be shown that where we have dropped the local terms and defined To obtain an expression for I λ s , let us first recall that the structure of the polarization tensors are given by the (associated) Legendre polynomials. Contracting with null vectors, only the term with the leading power in k survives (with no Kronecker delta's), whose coefficient is (2s−1)!!/[(2λ−1)!!(s−λ)!]. This means that where we used the fact that we get one factor of e iα for each contraction with a null vector, i.e. n i1 · · · n is ε s i1···is = e isα . Combining (A.102) and (A.12), we get

B In-In Results
In this appendix, we present details of the in-in computations of Section 4. In particular, we will give explicit expressions for the shape functions introduced in (4.18), (4.33), (4.36) and (4.41).
Results In Section 4, the results for the bispectra were defined in terms of a number of momentum-dependent functions. In the following, we give explicit integral expressions for these functions: • For s ≥ 2, the functions I (s) in (4.18) are given by iµs (c π , y) , (B.12) iµs (c π , y) , (B.13) iµs (c π , k 1 , k 2 , k 3 , y) , (B.14) where κ ij ≡ k i /k j and N s ≡ Z 0 s is the normalization constant defined in (A.77). The integrands are represented by the functions where G iµ 1 is in fact IR divergent. To avoid this issue, we integrate by parts and work with F (1) • The functions J (s) in (4.33) are given by iµs (c π , k 1 , k 2 , k 3 , y) iµs (c π , k 1 , k 2 , k 3 , z) , (B.28) • The functions K (s) in (4.36) are given by iµs (k 1 , k 2 , k 3 , y)

C Soft Limits
In this appendix, we will derive analytic formulas for the soft limits of the non-analytic parts of all correlation functions that we considered in this work.

C.1 ζζζ
We will focus on the squeezed limit of the scalar three-point function for the single-exchange diagram (cf. Fig. 6a), and consider even spins first. This leads to a non-analytic behavior if the quadratic mixing leg is taken to be soft. In the squeezed limit, k 1 k 2 ≈ k 3 , this contribution is given by iµs introduced in (B.7). The integrals in (C.11) and (C.12) can be computed analytically for arbitrary c π . To derive the results below, we will use the formula e −πµs/2 ∞ 0 dx x n H iµs (bx)e iax = (i/2) n √ π b n+1 F 21 (n + 3/2, µ s , (b − a)/2b) , (C. 13) where F 21 (a, µ s , z) ≡ Γ(a − 1 2 − iµ s )Γ(a − 1 2 + iµ s ) Γ(a) 2 F 1 a − 1 2 − iµ s , a − 1 2 + iµ s , a + s, z . (C.14) In the squeezed limit, k 1 c π k 3 , the result for (C.11) is Since we cannot take a soft limit of the integral (C.12), its general expression is rather complicated. For simplicity, let us display the results for the two limiting cases, c π = 1 and c π 1, for which (C.12) reduces to Notice that the mixing integral becomes independent of c π in the small c π limit. The function f (s) (c π ) is precisely the difference between evaluating the integral (C.12) with the mode function G for even spins, whereas the result for odd spins is given by replacing 1 + isinhπµ s → i cosh πµ s . The final answer is then obtained by summing the permutations (k 2 ↔ k 3 ) in (C.1). Momentum conservation impliesk Writing the spin as s = 2 + 1 for odd spins, with an integer, we get For odd spins, the leading terms cancel in the sum over the two permutations, and the squeezed limit scales as (k 1 /k 3 ) 5/2 . Note that the right-hand side of (C.25) is an even-degree polynomial of the angle. For even spin s = 2 , we have instead (see also [57]) where the leading terms add up and thus scale as (k 1 /k 3 ) 3/2 . The next-to-leading term at order (k 1 /k 3 ) 5/2 is an odd-degree polynomial of the angle.

C.2 γζζ
We also studied the tensor-scalar-scalar bispectrum γζζ . Its squeezed limit can be written as The integral (C.33) is given bỹ wheref (s) a polynomial of µ s that encodes the difference between evaluating the integral with G φ s ≡ arg B s − µ s ln 4c π , (C. 38) for even spins. The result for odd spins requires the replacement 1 + isinhπµ s → cosh πµ s . The final bispectrum is then obtained by summing over the permutations (k 2 ↔ k 3 ).