Higher Spin Supersymmetry at the Cosmological Collider: Sculpting SUSY Rilles in the CMB

We study the imprint of higher spin supermultiplets on cosmological correlators, namely the non-Gaussianity of the cosmic microwave background. Supersymmetry is used as a guide to introduce the contribution of fermionic higher spin particles, which have been neglected thus far in the literature. This necessarily introduces more than just a single additional fermionic superpartner, since the spectrum of massive, higher spin supermultiplets includes two propagating higher spin bosons and two propagating higher spin fermions, which all contribute to the three point function. As an example we consider the half-integer superspin $\textsf{Y}=s+1/2$ supermultiplet, which includes particles of spin values $j=s+1,~j=s+1/2,~j=s+1/2$ and $j=s$. We compute the curvature perturbation 3-point function for higher spin particle exchange and find that the known $P_{s}(\cos \theta)$ angular dependence is accompanied by superpartner contributions that scale as $P_{s+1}(\cos \theta)$ and $\sum_{m}P^{m}_{s} (\cos \theta)$, with $P_{s}$ and $P_{s} ^m$ defined as the Legendre and Associated Legendre polynomials respectively. We also compute the tensor-scalar-scalar 3-point function, and find a complicated angular dependence as an integral over products of Legendre and associated Legendre polynomials.


Introduction
Inflationary cosmology provides the initial conditions of standard cosmology, and a mechanism to explain the origin of the large scale structure of the universe. These initial conditions are manifest in the statistical properties of anisotropies in the cosmic microwave background (CMB) radiation, which in addition to being measured to an incredible precision [1], are well described by linearized cosmological perturbation theory. This latter fact means the statistical properties of the CMB are calculable, and combined with the conservation of the primordial curvature perturbation on super-horizon scales [2][3][4][5][6][7][8], makes possible deductions as to the precise particle content of the very early universe.
In particular, while single-field slow-roll inflation predicts adiabatic, Gaussian, nearly-scale invariant perturbations, interactions of the primordial curvature perturbation with new fields can generate deviations from Gaussianity, as encoded at lowest-order by the 3-point function ζζζ . Remarkably, this is sensitive even to fields that are heavier than the Hubble scale during inflation [9][10][11], and thus probes particles at energy scales far above those accessible by terrestrial colliders. The study of ζζζ as a particle detector has been termed 'Cosmological Collider Physics' [12] (see also [13]), and gained significant momentum due to the realization that interactions with higher spin bosons, namely the exchange of a massive spin-s boson, impart a characteristic angular dependence on the non-Gaussianity, ζ(k 1 )ζ(k 2 )ζ(k 3 ) ∝ P s (k 1 ·k 3 ) + k 2 ↔ k 3 , with P s (cos θ) the degree-s Legendre polynomial.
The study of higher spins has a long history dating all the way back to the founding of relativistic field theory 6 . Since then, higher spins have gained fame and attention in large part due to their role in string theory 7 8 , as well as their use in exploring the holographic principle 9 . Additionally, there is the old conjecture [17][18][19][20] that physics beyond Planckian energy scales will have higher symmetries emerging. From this point of view the study of higher spins can be understood as an attempt to classify and realize the various possibilities for these emerging symmetries. The study of manifestly supersymmetric higher spins is a natural extension of the above program, both because supersymmetry is a property of the underlying theory (as in the example of (super)string theory) and, more generally, because it is compatible with the relevant structures (like the symmetries of S-matrix [21]).
For these reasons we are interested in using irreducible representations of the supersymmetric extension of appropriate spacetime symmetry groups which involve higher spin particles. These irreps are classified and labeled by the eigenvalues of the Casimir Operators. In 4D these are the mass (m) and the superspin (Y) which takes either integer (Y = s) or half integer values (Y = s + 1/2) 10 . These representations include multiple representations of the non-supersymmetric spacetime symmetry group: For the massless case (m = 0) a supermultiplet with superspin Y includes massless particles with spins j = Y + 1/2 and j = Y, whereas for the massive case (m = 0) a supermultiplet with superspin Y includes massive particles with spins j = Y + 1/2, j = Y, j = Y and j = Y − 1/2. This implies that supersymmetrizing the cosmological collider does not simply require adding a single higher spin fermion, but rather additional fields as well. 6 First paper by Majorana in 1932 [14] followed by Dirac, Pauli, Fierz, Wigner and others. 7 For example, the UV softness of perturbative string scattering amplitudes originates from the freedom to exchange higher spin particles. 8 More recently, higher spin fields have played a role in constraining the self-consistency of inflation in string theory [15,16]. 9 All available, consistent, fully interacting higher spin theories (such as Vasiliev's or CS in 3D) require an AdS background and a spin two state. Both of these requirements are ingredients of AdS/CFT correspondence. 10 For the purpose of this discussion we will ignore infinite sized representations that go under the name of continuous (super)spin representations [22]. As a step towards combining this with inflationary cosmology, we consider an inflationary sector minimally coupled to a higher spin sector. The inflationary vacuum energy H 2 breaks supersymmetry, generating masses for the inflationary fermionic superpartners, while the higher spin sector, behaving as 'spectator fields' which do not contribute to H 2 and hence do not contribute to supersymmetry breaking, retain their on-shell supersymmetry. Given the candidate bosonic interactions proposed in the literature [12,13], the remnant on-shell supersymmetry of the higher spin sector uniquely determines the interactions of the higher spin fermions with the primordial curvature perturbation. From this one can compute the statistical properties of anisotropies in the CMB, and in this way, use the CMB as a detector for higher spin supersymmetry at the early universe's collider.
Each higher spin particle, as enumerated by the corresponding irreducible representation, induces a contribution to the 3-point function ζζζ , i.e. a signal at the cosmological collider, and in this work we explicitly calculate the non-Gaussianity due to these contributions. Our primary result is the prediction of higher spin supersymmetry for the angular dependence of the non-Gaussianity: we find that the P s contributions to the non-Gaussianity come in a characteristic pattern. Namely, every P s contribution to the non-Gaussianity is accompanied by a P s+1 contribution and a tower of associated Legendre polynomials P m s . The magnitudes, while Boltzmann suppressed, are related by supersymmetric considerations. This paper is organized in the following way. In section 2 we build the effective field theory that will be the framework for our calculations. Section 3, gives a very elementary review of massive, higher superspin supermultiplets focusing on the spectrum of propagating spin particles they include. Furthermore it demonstrates how to construct a supersymmetric extension of the previously studied class of effective field theories. In section 4, using these types of fermionic higher spin terms we consider an effective interacting Lagrangian up to first order in the higher spin fields. The effective Lagrangian is then used to calculate the contribution of higher spin fermions to the three point function ζζζ . The last section 5, gives a summary of our results and a short discussion for future directions, including the tensor-scalar-scalar 3-point function γζζ , which we compute for higher spin fermion exchange.

Setup in Effective Field Theory
In this work we consider a tripartite marriage of 4D, N = 1 supersymmetric higher spins [23][24][25][26][27][28][29][30][31][32][33][34][35][36] with the effective field theory (EFT) of inflation [37] and de Sitter supergravity [38][39][40]. We construct an effective field theory of a supersymmetric theory of higher spins in a quasi-de Sitter spacetime with spontaneously broken supersymmetry and spontaneously broken time-translation invariance. The goal of this construction is to minimally couple the higher spin sector and the inflationary sector, in such a way that the on-shell supersymmetry of the higher spin fields is maintained, despite supersymmetry being broken by the inflationary vacuum energy. The on-shell supersymmetry of the higher spin sector can then be used in conjunction with the effective field theory of inflation to dictate the couplings of higher spin fermions and bosons to primordial perturbations. This is distinct from the supersymmetric EFT of inflation [41] in that we do not focus on the gravity multiplet, but instead on the higher spin supermultiplets.
Our approach allows us to make progress despite not having a full theory of interacting higher spin de Sitter supergravity. While the setup may seem contrived, it bears some similarity with the interplay of supersymmetry and anomaly cancellation in string theory.
The interactions between effective field theory, supersymmetry and anomaly cancellation is not as direct as one might imagine. In some cases a complete superspace formulation or component-level supersymmetrization is known such as in the example of the 4D, N = 1 WZNW-QCD action [42][43][44]. In the development of Heterotic theory, anomaly cancellation required the addition of a new term in the action a la the famous Green-Schwartz mechanism [45], and this term required yet more terms (and years of calculation) in order to restore supersymmetry to the action [46]. Similarly, in Type II/M theories the gravitational anomaly on D/M-branes induced by loops of chiral fermions is canceled via anomaly inflow by a higher-derivative correction to the bulk action [47][48][49], and despite the anomaly not playing any direct role in supersymmetry, as for Green-Schwarz the supersymmetrization of the anomaly-canceling terms requires yet more terms be added, the calculation of which requires a herculean level of technical skill and detail [50]. The same issue applies to the SL(2, Z) symmetry of type IIB: restoring the invariance naively broken by the corrections requires the careful consideration of D-instantons [51] (and again these new terms must be supersymmetrized).
In each of these cases, a seemingly complete theory is found to be anomalous, and cancellation of the anomalies requires new terms. The new terms should respect the symmetries of the action, and generically additional terms must be found to accomplish this task. However, much can be learned even without a complete knowledge of all terms in the theory. For example, in the interim period between [47][48][49] and [50], the AdS/CFT correspondence was discovered [52]. Another example of this is of course the field of String Cosmology [53,54], which makes no recourse to the precise manner in which SL(2, Z) symmetry is maintained. With all this in mind, we construct an effective theory along the lines discussed above.
This approach can be illustrated with simple examples involving chiral superfields in N = 1 supersymmetry. The first non-trivial step is the 'sequestering' of supersymmetry breaking to the inflationary sector, analogous to the Randall-Sundrum scenario [55]. This can be done in a number of ways; one simple example is to take guidance from [56,57] and allow for a non-minimal Kahler potential, as in This exhibits a vacuum at X =X = Y =Ȳ = 0, wherein supersymmetry is broken by X, D X W = M . The scalar potential evaluated for X =X = 0 is given by a constant V = M 2 , leaving the scalar component of Y massless: m 2 Y ≡ ∂ YȲ V = 0 11 . Similarly, the fermion component of of Y remains massless since D Y W = 0. Thus the breaking of supersymmetry is not communicated to the on-shell mass spectra of Y , and the Y superfield retains on-shell N = 1 supersymmetry.
The utility of this approach is the enumeration of interactions and the tree-level couplings, since despite the sequestering of supersymmetry breaking, there are interactions between X and Y , which communicate the SUSY breaking at loop-level. Expanding X → δx and Y → δy, δx, δy ∈ R, one finds the interactions between scalar components, where the ... are higher order terms. Similarly, there are interactions between the fermionic components of X and the fermionic components of Y , and these two sets of interactions will communicate the SUSY breaking to Y . The structure of these interactions is governed by the underlying supersymmetry, which is spontaneously broken by X, and this structure dictates the effect that δy interactions can have on δx correlators.
In our setup, the higher spin sector is analogous to Y while the inflationary sector is analogous to X. It is the above sense in which the HS sector in our setup has on-shell supersymmetry. This can be used to enumerate the interactions and estimate the amplitude of correlation functions. However, this is not the full story: The next puzzle piece is the embedding of supersymmetry and supergravity into cosmological spacetimes.
This can be done via the the framework of de Sitter supergravity [38][39][40]. This theory describes the spontaneous breaking of supersymmetry with no field content other then the gravity multiplet and the goldstino of supersymmetry breaking. The latter can be expressed as a chiral superfield, S, satisfying a constraint equation, This constraint removes the scalar degree of freedom from S, leaving the fermionic component as the only propagating degree of freedom. The most general superpotential is given by, since any additional terms involving S vanish by the nilpotency constraint. Supersymmetry is broken by D S W = M , and the resulting scalar potential, for a minimal Kahler potential K = SS, is a cosmological constant given by which is positive for M > √ 3W 0 , giving a de Sitter spacetime.
Any additional matter sectors in de Sitter supergravity can easily be sequestered from the breaking of supersymmetry. For example, endowing S with a non-trivial Kahler geometry [56,57] and taking W 0 = 0, supersymmetry-breaking is purely in the S-direction provided that D T W = 0 in vacuum, which is guaranteed to be the case since D T W ∝ S = 0, leaving the fermionic component of T massless. Meanwhile, SUSY is broken by S, D S W = M , and the potential is a constant vacuum energy V = Λ = M , leaving the scalar components of T massless. Thus again, T retains on-shell supersymmetry.
To connect this with observational cosmology, and anisotropies in the cosmic microwave background radiation, it is necessary to consider fluctuations. Inflation models can be constructed in de Sitter supergravity along the lines of [56][57][58][59]. Consider a superfield Φ with the real part of the scalar component of Φ identified as the inflaton ϕ. The fluctuations of ϕ in spatially flat gauge are related to the curvature perturbation on uniform density hypersurfaces ζ in uniform field gauge via [60] where ≡ −Ḣ/H 2 is the inflationary slow-roll parameter. This defines the primordial power spectrum, in dimensionless form, where ζ k is a Fourier mode of the field ζ(x, t).
The curvature perturbation ζ can in turn be related to the Goldstone boson of spontaneously broken time-translation invariance, using the machinery of the effective field theory of inflation [37]. This starts from the realization that the time-dependence of the inflaton φ(t) breaks the time diffeomorphisms. The Goldstone boson can be included in the theory by a redefinition of the time-coordinate, with π(t, x) the Goldstone boson. This induces a field fluctuation, and thus corresponds to a curvature perturbation, This also generates a fluctuation to the 00 component of the metric, The interactions of π, and hence ζ, are dictated by the symmetry structure of the action, which is is broken to invariance under spatial rotations. This allows δg 00 to appear explicitly in the action, while δg ij can only appear with all indices contracted.
Similarly, the interactions of π, and hence ζ, with any additional fields are dictated by symmetry considerations. For fields with arbitrary spin, σ µ 1 ...µs , this leads to effective interactions of the form [12,13] where indices i runs over spatial directions: 1, 2, 3, and Λ is a UV scale. These terms descend from higherdimension operators built out of the metric and its derivatives, and have coupling constants that are a priori free parameters of the effective field theory. For example, in the spin-2 case, these terms arise from a coupling of a spin-2 field σ to the extrinsic curvature, √ −gK µν σ µν . The first term in (13) descends from δK µν σ µν while the second term arises from including the metric perturbation δg 00 δK µν σ µν [13]. There can also be additional terms at higher order in π and σ and terms with different distribution of derivatives (up to total derivatives).
We now arrive back at the tripartite marriage: We wish to connect the higher spin interactions in a quasi-dS space to supersymmetry. To do this, one could simply operate along effective field theory lines, and introduce interactions consistent with unbroken spatial rotations. However, an interesting possibility is to consider what we can learn from higher spin supersymmetry, and use this as guidance in constructing our effective field theory describing the interactions of the higher spin fermions. Towards this end, we now develop the machinery of higher spin supersymmetry.

Supersymmetric Higher Spins
The first Lagrangian description of supersymmetric, massless, higher spins in 4D Minkowski space was done in [61,62], using components with on-shell supersymmetry. A natural approach to the off-shell formulation is to use the superspace and superfield methods (see e.g. [63,64]). A superfield description of free supersymmetric massless, higher spin theories was presented for the first time in [23][24][25] for both Minkowski and AdS spaces. This approach has been further explored in [26][27][28][29]. Later studies of free supersymmetric, massless higher spin supermultiplets include [30][31][32].
On the other hand, the Lagrangian description of 4D massive supersymmetric spins for arbitrary values of spin is only known in the component formulation with on-shell supersymmetry [33,34], whereas the off-shell, superspace description has been developed up to superspin Y = 3/2 supermultiplet [35,36]. Nevertheless, independently of what the proper Lagrangian description is, we know that there are two types of such irreps (i) the integer superspin Y = s supermultiplets and (ii) the half integer superspin Y = s + 1/2 supermultiplets. Moreover we know that on-shell they describe two bosonic and two fermionic massive higher spin particles with spin values j = Y + 1/2, j = Y, j = Y and j = Y − 1/2. The half integer superspin Y = s + 1/2 supermultiplet, consisting of components Y = s + 1/2 : (s + 1, s + 1/2, s + 1/2, s) (14) and its on-shell, superspace description is given in terms of a real, bosonic superfield H α(s)α(s) 12 with the following on-shell conditions: where D αs is the superspace covariant derivative. Alternatively, the massive, integer Y = s superspin supermultiplet, comprised of components, has an on-shell superspace description in terms of a fermionic superfield Ψ α(s)α(s−1) with the on-shell equations: We start with the assumption that the higher spin sector respects supersymmetry and therefore can be organized into higher spin supermultiplets. This is extremely useful because supersymmetry will guide us to the introduction of higher spin fermions which have been neglected so far. Once their contribution is better understood, one may choose to drop the assumption of supersymmetry and study these fermionic contributions independently.
The strategy for finding the fermionic higher spin contributions is: (a) Start with the family of effective actions that lead to (13) after breaking the time translation invariance, and elevate them to superspace. This will automatically introduce all fermionic partners. (b) We project back down to a component description to reveal the interactions of the higher spin fermions. (c) Finally, we break supersymmetry appropriately in the inflaton sector.
The first step is to embed the bosonic, massive spin s particle in a massive higher spin supermultiplet described by some higher spin superfield. As we have seen, there are two ways of doing that, we can either use the integer or the half-integer superspin supermultiplet, with components (16) and (14) respectively. For concreteness, we make the latter choice (Y = s + 1/2) which means that our spin s particle will be accompanied by one bosonic higher spin particle j = s + 1 and two more fermionic higher spin particles j = s + 1/2. In this choice the highest propagating spin is s + 1. Similarly we embed the scalar curvature perturbation field in a scalar supermultiplet, which will of course introduce its fermionic superpartner, which we refer to as the inflatino 13 . A simple choice to describe such a scalar supermultiplet is to use a chiral superfield Φ.
Secondly, using superspace, we write quadratic and cubic interaction terms, between H α(s)α(s) and Φ which are linear in the higher spin superfield. The family of such superspace effective Lagrangians takes the form: 12 The notation α(s) signifies a string of s undotted indices α 1 α 2 ...α s which are symmetrized. This type of indices are the spinorial indices of a Weyl spinor of one chirality and take values 1 and 2 in 4D. Similarly forα(s), whereα are the spinorial indices of the opposite chirality Weyl spinor and take also two values1 and2 in 4D. 13 In general the fermionic partner of the curvature perturbation field can be identified as a linear combination of the inflatino and other fermions in the theory.
where I α(s)α(s) is linear in Φ and generates the quadratic interactions part, whereas J α(s)α(s) is quadratic in Φ and generates the cubic 14 part of the interactions. The most general ansatzes for I α(s)α(s) and J α(s)α(s) are, and Using this as a starting point, one can project the superspace Lagrangian to components and find the corresponding field theory (see [63,64] and detailed examples can be found in [29,68]). The result will include the entire spectrum of fields of the supersymmetric theory. In addition to the propagating spins, this includes the set of auxiliary fields required by supersymmetry in order to balance the bosonic and fermionic degrees of freedom and also make the symmetry manifest. However, these auxiliary fields do not have any dynamics and can be integrated out. By doing so, we obtain an effective theory with on-shell supersymmetry which includes two copies of the previously discussed bosonic higher spin interactions. That is because there are two higher spin bosons, one with spin s and one with spin s + 1. Additionally, we obtain terms that depend on the higher spin fermions (ψ α(s+1)α(s) , ξ α(s+1)α(s) ) and the 'inflatino' (χ α ). The interactions are given by, where ... indicates additional and higher-order terms.
As discussed in section 2, on-shell supersymmetry should be preserved only in the higher spin sector, and not in inflaton sector which breaks supersymmetry with the inflationary vacuum energy H 2 . In our effective Lagrangian this information can be entered by hand by removing any correlation between the coupling constants of the inflatino and those of the inflaton. For example we can assign the inflatino χ a mass m χ H. Typically one would then integrate out the inflatino, thereby eliminating all the contributions at linear order in the higher spin fermions (22). However it is the inclusion of such heavy fields that we 14 The reason why we are considering massive higher spin supermultiplets is a consequence of the Higuchi bound [65,66] plus the possibility of higher order mass-like interaction terms for the higher spin superfields, which we do not consider in this work. If that was not the case one should take into account the gauge symmetry of the higher spin (super)fields. The result of that would be that, the spectrum of the half integer supermultiplet will collapse from (s + 1, s + 1/2, s + 1/2, s) to (s + 1, s + 1/2) and more importantly the generator of cubic interactions J α(s)α(s) in (18) will become conserved higher spin supercurrents. Such supercurrents have been found in [67][68][69][70][71][72][73][74] and their structure is consistent with (20). are explicitly after in this work. Indeed, the higher spin fields themselves have mass m > s(s − 1)H by the Higuchi bound [65,66]. The breaking of supersymmetry will also induce differing loop corrections to the on-shell couplings of ζ to bosonic (22) and fermionic (22) higher spin fields. Depending on the precise details of the model, there may also be classical corrections to these parameters, for example, from a quartic interaction involving an additional scalar field that gains a VEV in the SUSY-breaking vacuum. For our analysis, we assume for simplicity that there are no such classical corrections, and that these couplings are equal at tree level. This does not alter the analysis in any way other then the overall prefactor of the result.
To make contact with the framework of effective field theories within we have to work, as presented in section 2, we must break the time translations part of the Poincaré group and write interaction terms which include fermionic higher spin particles up to linear order, and consider terms that contribute to the three point function ζζζ up to one loop corrections. Guided from the discussion of supersymmetry and taking into account contributions coming from the √ −g part of the action, one has to consider the following fermionic interaction Lagrangian: where the coupling toζ enters from the metric perturbation δg 00 , as in equation (13). Armed with this, we can now return to the cosmological collider.

Higher Spin Supersymmetry at the Cosmological Collider
The statistical correlations of temperature fluctuations in the cosmic microwave background descend from the initial conditions prepared for it by inflation. This can be computed via the Schwinger-Keldysh formalism, colloquially called the 'In-In' formalism, in which the choice of integration contour allows for ignorance as to the future evolution of the universe. Introduced to cosmology in [75], there are now many excellent reviews on this topic, see e.g. [76][77][78], and see [79] for a textbook treatment of the field theory aspects.
Correlation functions can be computed in this framework by splitting the Hamiltonian into a 'free' Hamiltonian H 0 and interaction Hamiltonian H int . From this, one can define the 'interaction picture' fields as having propagator determined solely by H 0 , and correlation functions of operators built from the full fields can be computed as contractions of the interaction picture fields with the interaction Hamiltonian.
More precisely, the expectation value of an operator W is given by, where H I int and W I are the interaction Hamiltonian and the operator W built out of interaction picture fields. The simplest quantity one can compute from this is the expectation value given a single insertion of the interaction Hamiltonian. In that case, the above expression reduces to where the ... corresponds to additional insertions of H int . For the case of the curvature perturbation 3-point function, W = ζ 3 , this picks out the intrinsic non-Gaussianity, first computed in [75]. This is the 3-point function induced by the self-interactions of ζ in an inflationary background. The result is given by This is typically expressed in the limit that one of the momenta is much smaller than the other two, in what is referred to as 'squeezed limit'. The result then takes a simplified form which corresponds to 'local shape' non-Gaussianity [80] with amplitude f N L given by The inflationary slow-roll conditions , η 1 thus imply the intrinsic non-Gaussianity in single-field slowroll inflation is extremely small, f N L 1.
Additional insertions of H int capture the effect of particle exchange. Given the slow-roll suppression of the intrinsic non-Gaussianity, this can easily be the dominant effect. It is in this sense that CMB non-Gaussianity is a particle detector, with inflation as the cosmological collider. In this work we are concerned with the non-Gaussianity induced via the exchange of a higher spin particle, as described by the insertion of two interaction Hamiltonians.

Effective Action and Relevant Interactions
We consider an effective action describing the interactions of a scalar ζ, a massive spin-1/2 field χ, and the propagating component fields of the massive, half-integer superspin Y = s + 1/2 supermultiplet, as they have been discussed previously. We consider our action as an expansion in higher spin fields, keeping up to linear order terms. We consider L = L ζ + L χ + L hs + O(hs 2 )...
(a) Exchange of a higher spin Fermion (b) Exchange of a higher spin boson Figure 1: In-In formalism Feynman Diagrams for exchange of a single higher spin particle.
In this work we are particularly interested in the impact of the fermionic higher spin particles, which has thus far been left unstudied. Motivated from the discussion in section 3 and (23), we take the Lagrangian of one fermionic spin-s + 1/2 particle ψ interacting with a dimensionless scalar ζ and the spin-1/2 particle χ, where λ s and g s are dimensionless coupling constants, Λ is a UV cutoff, and fermionic indices are contracted betweenχ and ψ. For simplicity we have taken κ s = 0 in (23), which leads to the same angular dependence for ζζζ as the two terms above. The relevant 'Feynman diagrams', or rather their equivalent in the Schwinger-Keldysh ("in-in" formalism) are shown in Figure 1.

higher spin Fields in de Sitter Space
To evaluate the three-point function in the Schwinger-Keldysh formalism, we first must have expressions for the free fields in de Sitter space.
To begin with, a scalar field φ is quantized in curved space as, The coefficients φ k (k, t) are referred to as "mode functions". The field φ(x) and mode functions φ k are related to two-point correlation functions as follows. The position-space two point function of a free-field φ is given by, The last equality defines the dimensionless power spectrum, ∆ 2 φ (k) ≡ k 3 2π 2 |φ k | 2 . A special case of the above is a scale-invariant spectrum. In this case, ∆ 2 φ (k) is a constant, which owes its name to the implied scaling symmetry of the two-point function, φ(x)φ(y) = φ(λx)ϕ(λy) . During inflation, this scaling symmetry of a massless scalar has its origins in the dilatation symmetry at late times in de Sitter space.
An important example is the curvature perturbation on uniform density hypersurfaces, ζ. This has mode function given by, where = −Ḣ/H 2 1 is the inflationary slow-roll parameter, and η is conformal time dη ≡ a −1 dt, which in de Sitter space is given by η = −1/(aH). This solution is found by explicitly solving the Klein-Gordon equation in de Sitter space. Inflation is a small deviation from de Sitter space, which converts the scaling with k to k −3/2+(ns−1)/2 , where n s defines the spectral index of the power spectrum, Another important example is massless spin-2, e.g. the graviton. Expanded in helicity states λ = ±2, the mode functions are given by [53] γ which, importantly, differs from the curvature perturbation (33) in part by an factor of 1/ √ . The power spectrum of primordial gravitational waves on large scales kη → 0 is then given by [53], which is a direct probe of the energy scale of inflation [81]. The ratio of the tensor power spectrum (36) to scalar power spectrum, termed the 'tensor-to-scalar ratio', is given by where again 1 is the inflationary slow-roll parameter.
In contrast with these two examples, for massive particles the two-point function and hence mode functions are suppressed on large scales and at late times. For a minimally-coupled massive scalar field σ, the two-point function has the exact solution [82] where µ ≡ m/H 2 − 9/4, and H iµ is the Hankel function of the first kind. This corresponds to a mode function σ k ∼ Hη 3/2 e −πµ/2 = a(η) −3/2 e −πµ/2 / √ H, the latter equality using η = −1/aH during inflation.
Now we turn to massive particles with spin. These are constrained by the Higuchi bound to have mass satisfying m 2 ≥ s(s − 1)H 2 . As for scalars, the isometries of dS fixes the scaling of the two-point correlation function of spinning fields [12], which takes the form where all Lorentz indices are contracted with s copies of a null vector, and ∆ is the scaling dimension of the field, For heavy fields, or more precisely the "principle series" [82], one has Re∆ = 3/2, and the mode function is simply The prefactors follow from dimensional analysis, and intuition from solving the Klein-Gordon equation for heavy fields, which leads to an additional e −πµs suppression. Importantly, this applies to general operators with spin, and not to just to bosonic (integer spin) operators, but to fermionic (half-integer spin) as well, and the two-point function of half-integer operators is similarly constrained to scale as k 2∆−3 as in (39).
In this work we will focus on the angular dependence of correlation functions. Given this, for simplicity we ignore the e iµs phase, though we note that this can lead to oscillations in k-space [13], and thus is of potential interest. Dropping this phase, the mode functions for spin-s and spin-(s + 1/2) are at late times given by The scaling with a(η) follows from solving the mode-function equation of motion explicitly 15 , while the See (B.76) and (C.7) of [82]: Using η=1/(aH) and differing factors of H follow from dimensional analysis 16 . This matches with the known result for a heavy spin-1/2 particle in de Sitter space [83]: a spin-1/2 fermion χ with mass m > H has mode function given by, Now we can make this more precise. A massive spin-s boson may be split into helicity components as, and then decomposed into fields of n polarization directions by projecting the spinning field σ µ 1 ...µs onto spatial slices, i.e. via the decomposition where η is the time coordinate. Here, the s index refers to the spin, n refers to the 'spatial spin,' and λ is the helicity of the field. Thus ε λ is...in is a normalized, totally symmetric tensor with spin s and helicity λ. The σ λ n,s satisfy σ λ n,s = 0 for n < |λ| [13]. The quantity that appears in scattering with scalars is λ = 0 and n = s (for more details see [13] or Appendix A), i.e. the quantity σ 0 s,s . Explicitly solving for the mode function, one finds [13] σ 0 s,s (k, η) Moreover, for the above λ = 0 helicity state, one has the important relation, with θ defined as the angle betweenq andk. This follows from more general relations for spin-s polarization vectors, which are detailed in the thesis [82], and given in Appendix A.
Similar to the bosonic case, a massive spin-(s + 1/2) 4-component fermion may be split into helicity components as, where α is a fermionic index and µ 1 ...µ s are bosonic indices, and projected onto spatial slices via the decomposition, where again η is the time coordinate. We construct the spin-(s + 1/2) polarization vectors as a tensor product of spin-1/2 and spin-s. That is, we decompose, with ξ λ a spin-1/2 eigenspinor of helicity λ and ε λ the spin-s polarization vector of helicity λ. The general decomposition of the fermion field can then be written as, Similar to the bosonic case, the fermions can be split into helicity states, and it is the helicity 0 ± 1/2 which contributes to the loop in Figure 1a. More explicitly, the relevant mode function has the form, ψ 0λ s,s a(η) s−3/2 e −πµs/2 .
One can use this, along with the mode functions (44) and (33), and the interaction Lagrangian (30), to compute correlation functions of ζ involving intermediate states of fermions.

Non-Gaussianity from higher spin Particle Exchange
The correlator we wish to compute is of three ζ(k, η) at lates times, η → 0, and in the limit that one of the momenta is small k 1 k 2 , k 3 , i.e. the quantity, given the insertion of interaction vertices between the scalar ζ and higher spin fields. The result for exchange of a higher spin boson are given in [12,13]. While [12] focused on the scaling and angular dependence, [13] explicitly solved the Klein-Gordon equation for higher spin bosonic fields and from this was able to compute all expressions exactly. In our analysis we will follow [12] and focus on the amplitude and angular dependence.
The 3-point function resulting from higher spin boson exchange, diagram 1b, is given by [13], where I (s) (µ s , c π , k 1 , k 3 , k 3 ) is a complicated function of momenta given in the Appendices of [13], and P s the Legendre polynomial. This is characterized by a dimensionless coupling α s , which in the notation of our (13), and taking the Goldstone boson parameters of [13] to be c π = 1 and f π = m pl , is given by, This corresponds to a non-Gaussianity parameter of with shape function and a characteristic angular dependence, with θ the angle between k 1 and k 3 .
We now turn to the fermions. For the Lagrangian (30), to compute the diagram Figure 1a there are two interaction Hamiltonians, given by and with i indices summed over, and where indicates a derivative with respect to conformal time. The 3-point function of three ζ(k, η), at late times η → 0, is given by To proceed we expand H int1 and H int2 in momentum space, which results in 7 momentum integrals, d 3 q 1 ....d 3 q 7 . We define the momenta as follows: let k 1 be the ζ in the left interaction vertex of Figure 1a, and the momentum of χ in the loop and + k 1 the momentum of ψ in the loop. The 'outgoing' ζ momenta are defined as k 2 and k 3 ; under Wick contractions we will need to sum over k 2 → k 3 .
This gives, approximating the mode functions by their forms given in section 4.2, where q ≡ + k 1 . We can factorize this into three pieces: with time-integrals, momentum integrals, and an overall prefactor of .
The loop integral I 2 is UV-divergent, and we apply a cutoff at H, withq now given byq where the latter equality follows working in the limit k 1 → 0. This simplifies the sum over spin-1/2 helicities, as the ξ(k) are normalized to 1. Thus we have, The remaining integral over angles can be performed analytically. Defining k 1 as making angle θ 1 = 0 in the {x, y} plane, and k 3 as making angle θ 13 , such thatk 1 ·k 3 = cos θ 13 ,q ·k 1 = cos θ,q · k 3 = cos(θ − θ 13 ), the integral can be written as, We then use the identity 17 , from which one can evaluate the integral explicitly. The result is with coefficients Finally we can perform the time integration. Using the explicit ζ mode function (33), one can analytically compute these integrals to find, Putting the pieces together, we find for the non-Gaussianity, which can be brought to a canonical form, where ∆ 2 ζ is the dimensionless primordial power spectrum, S(k 1 , k 2 , k 3 ) is a function of the ratios of k i , and A is all remaining prefactors.
The corresponding non-Gaussianity parameter is given by, while the angular dependence is now with the coefficients c m given by (74). The schematic form of the shape function S can be read off from (76), but the exact expression requires solving for the exact mode-functions of the higher spin particles in de Sitter space.

The predictions of Higher Spin Supersymmetry
We can now read-off result for the three-point function ζζζ given the higher spin supermultiplet. We simply add the contributions from the particle content of the half-integer superspin Y = s + 1/2 supermultiplet (14), given the non-Gaussianity from each spin derived in the previous section. The result is where the last indicates that the three terms in (80) have angular dependence given by P s+1 , m P m s , and P s respectively.
The quantitative amplitude of this signal is, as in the non-supersymmetric bosonic case [13], generally small f N L O(1). The primary obstruction making f N L any larger than this is perturbativity of the interaction strength, which at the very least, requires λ s (H/Λ) s−1 1 and g s (H/Λ) s 1, as these are the effective interaction strengths, e.g. appearing in (63). Non-Gaussianity of this size is not what would traditionally be referred to as 'large', but it can be considerably larger than the slow-roll suppressed singlefield slow-roll value, and is within reach for CMB-S4 [81].
The analysis can be repeated for the case of embedding the higher spin particles inside the integer superspin supermultiplet (16) instead of the half integer one we have used as an example. In that case the particle contained is ( s + 1/2, s, s, s − 1/2)therefore one can immediately read off the result. The known P s (cos θ) dependence of spin-s bosons [12,13] is accompanied by two towers of associated Legendre polynomials, m P m s and m P m s−1 , from the s + 1/2 and s − 1/2 fermions respectively.

Discussion
Precision measurements of the cosmic microwave background provide an unprecedented opportunity to search for new physics in the early universe. The 3-point function of primordial curvature perturbations, ζζζ , colloquially referred to as the non-Gaussianity, is sensitive to any new degrees of freedom, including those that are naively too heavy to be excited. One of the most striking results of this research program is the non-Gaussianity due to higher spin particles, and in particular the angular dependence ζζζ ∝ P s (cos θ) due to the exchange of a single spin-s boson [12]. This prompted a flurry of activity, and possibilities for observing [84][85][86][87][88][89][90] the signature of higher spin particles.
higher spin fermions have heretofore been left out of this discussion, but insofar as higher spin theory is understood as a limit of quantum gravity, namely superstring theory, fermions are built into the theory. This is required by the supersymmetric nature of the theory, which is itself a powerful tool for the incorporation of fermions into string theory, e.g. the construction of fermionic D-brane actions is accomplished by relying on the underlying supersymmetry of the theory [91][92][93][94]. Guided by this, and building on recent developments in the construction of supersymmetric higher spin theories [23][24][25][26][27][28][29][30][31][32][33][34][35][36], we have studied the imprint of higher spin supersymmetry at the cosmological collider.
The main result of this paper is a characteristic pattern of the angular dependence of ζζζ due to the exchange of higher spin superpartners. We find the P s (cos θ) signature of higher spin boson exchange, with θ the angle between the short and long wavelength modes, comes along with a P s+1 (cos θ) and a tower of associated Legendre polynomials, arising from a spin-s + 1 boson and a pair of spin-s + 1/2 fermions. For a variant description of higher spin supermultiplet, the partner contributions can be instead two towers of associated Legendre polynomials. The amplitude of the signal is generically not large by comparison to other known sources (e.g. [95,96]), as already known for the non-supersymmetric bosonic case [13], so it is indeed the angular dependence which gives this signal its elevated status. Given this, in this work we have not endeavored to do a rigorous and precise calculation of the shape-function, which requires explicitly solving the mode-functions and computing involved integrals [13]. This latter difficulty motivated the development of the cosmological bootstrap [97], which might be a promising direction to take this work as well.
Remarkably, we have been able to derive these results despite not having a complete theory or model realization of higher spin supergravity inflation. Progress despite incomplete knowledge is a familiar situation in theoretical physics, for example, supersymmetry and the Green-Schwarz mechanism, all work to date pertaining to M-theory [98,99], or in a more recent context, Double Field Theory [100]. To overcome this, we have constructed an effective theory that combines higher spin supersymmetry with de Sitter supergravity and the effective field theory of inflation, to describe a higher spin sector minimally coupled to the inflationary sector such that the higher spin sector retains on-shell supersymmetry. This allowed us to use supersymmetry considerations to deduce the field content and interactions of the higher spin fields with the curvature perturbation.
There are a number of ways forward from here. We have not considered yet the interactions with the graviton, or for that the matter, the gravitino. The former of these, corresponding to primordial gravitational waves γ, itself can lead to an interesting three point function, γζζ , probed by crosscorrelation with CMB B-mode polarization. Starting from an effective field theory guided guess for the relevant interaction, the γζζ computed can be straight-forwardly worked out in a fashion similar to section 4. We provide the calculations for this in Appendix B, and here we give the result: where λ = ±2 is the helicity of the external graviton, ij is the spin-2 polarization tensor,P λ s (x) ≡ (1 − x 2 ) −λ/2 P λ s (x) and E λ 2 (k 1 ·k 3 ) = λ ij (k 1 )k i 3k j 3 as in [13], and we have put all k-dependence in the function S. This result is characterized by an angular dependence that is an integral over Legendre and associated Legendre polynomials, where θ q is the angleq makes in the plane. In the supersymmetric context, this would be joined with contributions from additional interactions. This requires careful consideration of the gravity multiplet, as is the focus of the supersymmetric EFT of inflation [41]. We postpone this analysis to future work.
On the theoretical front, an important next step is to construct the full theory of spontaneously broken supersymmetry (as in de Sitter supergravity and the supersymmetric EFT of inflation) and interacting higher spin fields. From this one can generalize the analysis here to situations where the higher spin fields themselves contribute to the supersymmetry breaking, or perhaps even drive inflation. We leave this possibility, and a host of observational implications, to future work.
As a concluding remark, we would like to express and share our enthusiasm for the work of scientists who are searching for signals of supersymmetry in the cosmos, a sentiment expressed by one of the authors in [101]. As argued for in this work, the non-Gaussianity of the CMB may prove to be a powerful tool of discovery, and with some good fortune, perhaps more and different such tools will later emerge for the SUSY search at the Cosmic Collider.

A Spin-s polarization vectors
This appendix discusses some relevant preliminaries and definitions relating to the free theory of higher spin fields in de Sitter space, which can also be found in [13]. Following [13], vectors will be denoted here in boldface, e.g. k.
It is convenient to project the spinning field, σ µ 1 ...µs onto spatial slices, which we can then write as Here, the s index refers to the spin, n refers to the 'spatial spin,' and λ is the helicity of the field. ε λ is...in is a normalized, totally symmetric tensor with spin s and helicity λ. It must satisfy: corresponding to the symmetric, traceless, and transverse properties. These properties of the polarization tensor imply that we can decompose it into transverse and longitudinal parts as, where ε λ i 1 ...i λ is a maximally transverse polarization tensor, constructed out of polarization vectors ε ± that are perpendicular tok. We must have that ε + = (ε − ) * , so that ε λ i 1 ...i λ can be specified, up to a phase, by a single polarization vector ε. We have also defined f i λ+1 ...is as the longitudinal part of the associated Legendre polynomial, after contraction with momenta.
We then define, The symmetry properties of ε imply that F λ s takes the form [13], in d=3 spatial dimensions, where P λ s is the associated Legendre polynomial. For the special case of λ = 0, which appears in the calculation of 3-point functions after enforcing momentum conservation, the othonormality of differing helicity states λ and λ , and the transverse property (A.2), one has q i 1 ....q is ε 0 i 1 ...is (k) ∝ P s (q ·k), (A.7) with magnitude |q 1 | s , leading to the characteristic angular dependence of the three-point function for spin-s boson exchange. Moreover, the transverse property also implies that the only σ 0 n,s that enters the correlation function is n = s.

B Details of γζζ Calculation
In this appendix, we further explicate the derivation of the tensor-scalar-scalar correlation function. We will limit our analysis to the single-exchange diagram shown in Figure 2.
This diagram has the same form as Figure 1a, however now we have an external graviton carrying momentum k 1 instead of ζ. The relevant interaction Lagrangian we will consider is L =λ s Λ s−2 ∂ i 1 ...i s−2γ i s−1 isχ ψ i 1 ...is + g s Λ sζ ∂ i 1 ...is ζχψ i 1 ...is + h.c., (B.1) which corresponds to two interaction Hamiltonians As in the ζζζ calculation, we would like to expand in Fourier modes. The graviton can be expanded in helicity modes as given in [13]: where the graviton mode function, γ λ k , is given by With the mode functions in hand we can compute the tensor-scalar-scalar three point function γζζ . As before, we would like to expand each interaction Hamiltonian in momentum space and compute the correlator. Much of the calculation remains the same as in the ζζζ case, however, there are some subtleties. The angular dependence due to H int 2 remains largely the same, however now due to the fact that λ = 0, we have H int 2 ∝ E λ 2 (k 3 ·q)P λ s (k 3 ·q), rather than simply P s (k 3 ·q) as in the ζζζ case. This arises from the definition (A.5), and we have defined E λ 2 (k 3 ·q) = ε λ ij (k 3 )q iqj . The angular dependence of H int 1 is similarly complicated, due to the contraction with γ i s−1 ,is .