Higher Spin Dark Matter

Little is known about dark matter beyond the fact that it does not interact with the standard model or itself, or else does so incredibly weakly. A natural candidate, given the history of no-go theorems against their interactions, are higher spin fields. Here we develop the scenario of higher spin (spin $s>2$) dark matter. We show that the gravitational production of superheavy bosonic higher spin fields during inflation can provide all the dark matter we observe today. We consider the observable signatures, and find a potential characteristic signature of bosonic higher spin dark matter in directional direct detection; we find that there are distinct spin-dependent contributions to the double differential recoil rate, which complement the oscillatory imprint of higher spin fields in the cosmic microwave background. We consider the extension to higher spin fermions and supersymmetric higher spins.


I. INTRODUCTION
The precise identity of dark matter (DM) remains a mystery, despite decades of theorizing and detection efforts. Observations suggest that the dark matter does not interact with the standard model, or does so extremely weakly, creating significant room for model builders. Taken at face value, the observational evidence is at odds with the conventional origin story of dark matter, namely a thermal history, wherein the dark matter was initially in a state of thermal equilibrium with the standard model, sustained by interactions.
There are now many alternative dark matter origin stories. A particularly compelling possibility, by virtue of its simplicity, is the genesis of dark matter via gravitational particle production (GPP) in the early universe [1][2][3][4][5], e.g., during cosmic inflation. A generic feature of inflation is that the exponential expansion acts as a gravitational amplifier for particle production. While many of these particles would be redshifted, some, depending on their intrinsic properties such as mass and spin, can survive as a relic after inflation ends. This idea was introduced in the context of superheavy 'WIMPzilla' dark matter, characterized by dark matter masses greater than the Hubble scale at the end of inflation (m DM > H) [1][2][3][4][5][6], and has since been explored in a variety of cosmological contexts (see e.g. [7][8][9][10][11][12]).
On the other hand, inflation is known to exhibit a UV sensitivity [13], motivating the search for a UV completion of inflation in theories beyond the standard model, such as string theory. In this context we may ask which states could be generically produced in an inflationary model that is connected to string theory. That massive higher spin particles are a natural consequence of string theory, and that generic cosmological inflation models induce particle production, is suggestive of the potential implications for connections to other physics, namely the dark matter problem. Moreover, the Higuchi bound on the mass of higher spin fields, m 2 ≥ s(s − 1)H 2 , which must be satisfied at all times during an early universe genesis mechanism, naturally suggests dark matter in the superheavy regime, and hence, in light of [1][2][3][4][5], a gravitational origin of higher spin dark matter.
In this work we consider inflationary production and cosmological implications of higher spin particles and find that they can naturally serve as 100 percent of the dark matter: Higher Spin Dark Matter (HSDM). We consider the implications of a small interaction of HSDM with the standard model, and find a characteristic angular dependence of nuclear recoil events for direct detection experiments. This mirrors the angular dependence in the cosmic microwave background non-Gaussianity that is predicted due to the the production of higher spins during inflation [24,37,39].
The structure of the paper is as follows: in Section II we introduce the relevant higher spin formalism. In Section III we calculate the gravitational production of higher spin dark matter (HSDM) and show that there is a parameter space such that higher spin particles can account for all the dark matter. In Section IV we discuss the possibilities for directional direct detection and show that there is a spin dependent contribution to the double differential recoil rate. Lastly, in Section V we speculate on other possible observable avenues and conclude with a discussion in Section VI.
There are well known 'no-go' theorems that significantly limit the interactions of HS particles in a selfconsistent quantum field theory. Generally, such theorems make it the case that in flat space, massless HS particles are forbidden from interacting with electromagnetism or gravity 1 [74][75][76][77][78][79]. Two notable 'no-go' theorems are Weinberg's theorem [74], which necessitates that, in flat space, there are no long range interactions with spin greater than two, and the Coleman-Mandula theorem [75], which demonstrates that, assuming an Smatrix and finite degrees of freedom, there can be no conserved higher spin charges associated with particles of s > 2.
A caveat to these arguments is massive higher spin theories. The mass term explicitly breaks the higher spin gauge invariance, such that there is no conserved current, and hence no conflict with the Coleman-Mandula theorem. As such, massive higher spins are not plagued by the same restrictions due to no-go theorems [80]. Indeed, massive higher spin excitations are intrinsic to string theory, and comprise the Regge trajectories. Higher spin fields have been considered in studies of inflation in string theory [81][82][83], and in the AdS/CFT correspondence [46,48,84]. It is thought that the tensionless limit of string theory is a higher spin field theory [85][86][87][88][89], and it has been suggested that string theory itself is a symmetry broken phase of a HS field theory [46,[90][91][92].
While the full theory of higher spins is not known, progress can be made by enumerating the irreducible representations of the spacetime symmetry group, thereby identifying the building the blocks of the theory. Although the representation theory of HS fields in a general Freidmann-Robertson-Walker spacetime is not known, the representations are known for flat space and (A)dS.
For cosmological purposes, in particular during inflation, we may make use of the results for de Sitter space.
In this context, a lower bound on the higher spin mass is given by the Higuchi bound: m 2 ≥ s(s − 1)H 2 [93,94]. Beyond this, fields can be organized into three categories of unitary, irreducible representations of the spacetime isometry group [95,96]. They are the principal series: the complementary series: and the discrete series: In addition, In this work, we focus on the complementary series and principal series representations. The evolution of HS fields in de Sitter space was derived in [37], which we summarize below. The spin-s generalization of Klein-Gordon equation is the Casimir eigenvalue equation of the de Sitter group [37], This is supplemented by constraint equations corresponding to transverse and traceless conditions on σ. To solve this equation, 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 polarization vector ε λ i1···in is symmetric, transverse, and traceless; for details, see [37].
The helicity-λ mode function with n = |λ| number of polarization directions satisfies, This admits an exact solution, given by [37], where µ s is defined as The normalization coefficients are given by, and with The other mode functions can then be obtained iteratively from the recursion relation: . (14) Care should be taken when considering these mode functions, since they are normalized with respect to σ with all lower indices. The quantity of physical interest is the two point function of two contracted σ, i.e., More generally, for σ λ n,s , we are interested in the twopoint correlation function, The remaining contraction of polarization vectors can be computed from [37] where we have used the normalization ε s i1···is ε s * i1···is = 2 s [37].

III. GRAVITATIONAL PRODUCTION OF HIGHER SPIN DARK MATTER
The Higuchi bound m 2 > s(s − 1)H 2 [93,94] suggests that higher spin fields as realized in nature, insofar as they can be described by a single effective field theory in both the very early universe and in the late universe, should be cosmologically heavy. Guided by past literature [1][2][3][4][5], it is logical then to consider gravitational particle production as a genesis mechanism for higher spin dark matter. With this in mind, our first goal is to make a conservative estimate of the gravitational production of higher spin particles in the very early universe. This will serve as a proof of principle of gravitational particle production (GPP) as a genesis mechanism for higher spin dark matter (HSDM).
We will consider only the gravitational production during inflation, and not the transition between inflation and the radiation dominated phase that is the usual focus of works on GPP of dark matter [1-5, 9, 10]. This simplification is not made for convenience, but rather due to the limited knowledge of higher spin field theories as discussed previously. Our calculation provides a lower bound on the production, suitable for a demonstration that early universe GPP can provide enough higher spin particles to explain the observed DM abundance. We expect a more detailed analysis (e.g., directly from string theory) to change slightly the quantitative relationship between s, m, and H, that leads to the correct relic density, but not the qualitative result.
We focus on HS fields that during inflation are in either the complementary or principal series, defined by Eqs. (2) and (1) respectively. We make use of the fact that cosmic history is thought to be book-ended by de Sitter phases: cosmic inflation in the first moments and dark energy in the present. The full structure of the higher spin theory is not known in for the intervening time period. However, post-inflation we are left with a collection of non-relativistic particles that simply redshift as matter, and hence we do not need to consider the detailed field theory dynamics. This negates the need to have complete knowledge of higher spin theories. In general, gravitational particle production occurs when the field mass, including all contributions from quantum and gravitational effects, changes nonadiabatically. The canonical example is the primordial curvature, ζ, which obeys the equation of motion (in ex- Adiabaticity is violated when k = √ 2/|η| aH, i.e., when a given mode exits the horizon. The resulting particle production can be thought of as Hawking radiation emitted by the de Sitter horizon [97], see e.g. the discussion in [98], and indeed computed using conventional particle production methods [99,100]. Alternatively, the equation of motion can be solved exactly at all times as a function of kη, see e.g. [101].
For a massive scalar field the effective mass is similarly given by [101] For m/H 1, adiabaticity is violated at k µ/|η| = µaH where µ ≡ m/H [10]. The scale k * = µaH defines an effective horizon, and all modes which have k > k * at the end of inflation, i.e. which have exited this effective horizon, have undergone particle production during inflation. This is the principle behind superheavy dark matter. The total energy density in a massive scalar can be computed as corresponding to a particle number, n = ρ ϕ /m ϕ , The phenomenology of superheavy dark matter is thus determined by the dark matter mass and the energy scale of inflation. Due to the decay outside the horizon of the heavy field, along with the phase space suppression in the limit k → 0, the integral is dominated by the contribution from the upper bound. The number density scales with H 3 , where H is the Hubble constant during inflation. The relic density can be tuned to the observed value by tuning the mass m DM as a function of H; for the canonical super-heavy dark matter, one finds m DM H leads to the observed density [1][2][3][4][5].
The generalization of this to the energy density of a bosonic higher spin field σ produced gravitationally during inflation is given by, where the sum is over all values of n and λ, and the contraction of polarization vectors is given by Eq. (17).
To evaluate this we begin from the equation of motion Eq. (8). To put this in a more familiar form, we rescale σ in Eq. (8) asσ = η 1−λ σ. This removes the first-derivative term, leading to, which describes a simple harmonic oscillator with a timedependent mass. We can see that the λ dependence has cancelled identically, and the effective mass is given by, where we define the quantity µ * ≡ s(s − 1) − µ 2 s , with µ s given by Eq. (9).
In analogy with a massive scalar, adiabaticity is maximally violated (|ω k |/ω 2 k is peaked) when k = √ 2µ * aH µ * aH. Thus, we take the UV cutoff of the energy density integral to be this scale of adiabaticity violation. Incidentally, this cutoff is also the dividing line between relativistic and non-relativistic particles. Thus, independent of discussions of adiabaticity violation, this approach is equivalent to considering only those perturbations which become non-relativistic already during inflation and thus can be treated as non-relativistic at all times following inflation.
In principle, Eq. (22) involves the sum over all excitations of the spin-s field σ. However, the dominant excitation of a spin-s field near the Higuchi bound is σ 0 s,s . This can be seen qualitatively from Eqs. (8) and (13). Eq. (8) will have a maximum value at λ = 0 and all other modes with λ > 0 will be suppressed in comparison. It can then be seen from Eq. (13) that the dominant contribution will be the n = s state, due to the contribution from B m,n . Curiously, this is the same mode which is considered in the 'cosmological collider physics' program [24][25][26], leading to the characteristic Legendre polynomial angular dependence of 3-point functions [24,37,39].
Given the subdominance of all other modes, we approximate the full density as that which comes from σ 0 s,s . We use the recursion relation Eq. (13) to explicitly numerically calculate the mode functions, σ 0 n,s . We then calculate the density of excitations at the end of inflation using where t i denotes the end of inflation, and µ * is given by Eq. (24). The factorial prefactors come from evaluating the contraction of λ = 0, n = s, polarization vectors. The integral is again dominated by the contribution from the upper bound, and in the numerics that follow we make use of this approximation. Further, since µ * s > 1, by the Higuchi bound, the DM density is dominated by particles that are still sub-horizon at the end of inflation (k/a < H) and those which re-enter the horizon immediately following inflation k/a ∼ H. Since m/H is greater than 1 for all subsequent times, we can approximate these particles as non-relativistic for all times post-inflation. From this we define the present DM density as, where ρ inflation is defined by the integral Eq. (25), and a(t i ) is the scale factor of the universe at the end of inflation (we normalize a = 1 today). The acceptable parameter space for our HSDM model will be that for which the relic density Eq. (26) matches the observed dark matter abundance. The observed dark matter density is given by ρ DM0 = 3m 2 pl H 2 0 Ω CDM where Ω CDM is the abundance observed to be Ω CDM h 2 0.12 [102], and H 0 ≡ 100h km/s/Mpc = 2.13h × 10 −33 eV with h 0.7. From this we find ρ DM0 evaluates to ρ DM0 = 3.95 × 10 −11 eV 4 . Meanwhile, the redshift factor in Eq. (26) can be simplified by expand a(t i ) as a ratio of redshifts, where eq refers to matter-radiation equality. We have z eq 3400 and (1 + z) ∝ T for z z eq . Taking T eq 0.8eV 1eV and instant reheating T re (g * π 2 /30) −1/4 Hm pl , with g * ∼ 100, the above becomes a(t i ) = 1.43 × 10 −18 eV H . Putting things together, we find that the DM density after inflation must satisfy, with ρ(t i ) given by Eq. (25).
To gain intuition as to the range of masses that can provide the correct relic density, we note the peculiar case of the complementary series, which occupies a narrow range of masses just above the Higuchi bound. At the saturation limit of the Higuchi bound, we approach partial masslessness and states become gauge redundancies [38]. To avoid this we deform away from the Higuchi bound by a small amount 2δ and consider m 2 /H 2 = s(s − 1) + 2δ. In this limit, µ s , Eq. (9), becomes One can appreciate from the above that µ s is purely imaginary. This is a feature of the complementary series Eq. (2), which corresponds to 0 < δ < 1/4. It follows that the exponential suppression, which one might anticipate for excitations of a heavy field, and is encoded in Eq. (10), becomes a phase factor -i.e., is not a suppression at all. In Fig. 1 we illustrate the values of s and H for which the correct relic density of particles is produced. By varying the mass, points in the colored regions can achieve the correct relic density. The colors pink and blue denote masses in the complementary series and principle series respectively.
The lower and upper dashed lines of Fig. 1 denote the lower edge of the complementary series (the Higuchi bound) and upper edge of the complementary series. The parameter space below the lower edge is ruled out: the DM particle mass cannot be lower than the Higuchi bound, and decreasing H while leaving m/H fixed will lead to an underproduction of DM during inflation. The upper bound of this band represents the upper limit on particles with masses in the complementary series. If m/H and s remain fixed, increasing H will lead to an overproduction of the DM. A simple way out is to increase m/H, leaving the complementary series, and thereby generating an exponential suppression of the amount of DM, which is denoted by the light blue portion of Figure 1.
By considering masses slightly in the principal series we are able to obtain the correct relic density over the entire blue region of Fig. 1. This imposes a relation between m and H, at fixed values of s, as shown in Fig. 2. From left to right, fixed curves from s = 2 − 8, respectively show that for any spin, there is a continuous range of allowed H values with increasing mass. The allowed parameter space is bounded by the Planck mass, above which is forbidden and denoted by the grey shaded region. Although we have only explicitly shown up to s = 8, one can appreciate that there is a much larger space where solutions will be present, up to m/H = 10 18 and correspondingly increasing spin. This allows for a wide range of H and s values, which makes HSDM amenable to a variety of inflation models.

IV. DIRECTIONAL DIRECT DETECTION
Direct detection is a prominent detection strategy for dark matter. In this approach, one hopes to observe nuclear recoil events generated by scattering of incoming dark matter particles. As realized early on [103], a signature prediction of the motion of the earth through the enveloping dark matter halo is a preferred direction of nuclear recoil events, suggesting an approach known as directional direct detection [104]. More recently, it has been demonstrated [105][106][107] that the angular dependence of the directional direct detection signal can distinguish between spin-0, and spin-1/2, and spin-1 dark matter. It is logical, therefore, to consider the directional direct detection signature of higher spin dark matter.
Any such direct detection signal is premised upon an interaction of dark matter with the standard model. As discussed in Sec. II, while higher spin interactions are naively strongly constrained, this is relaxed for massive higher spins. The interactions of massive higher spin fields can be understood as a low energy effective field theory, with a UV completion given by string theory. The possible interactions have been studied in detail and enumerated in e.g. [62,63,66,67].
We will consider a simple interaction between a higher spin boson and a standard model nucleus, eg., Xenon, that we model as a Dirac fermion. We construct an interaction through a derivative coupling of the higher spin boson to the fermion vector current. We consider the low energy effective interaction Lagrangian, where Λ is a UV scale that corresponds to the cutoff of the massive higher spin theory, and g s is a coupling constant. Note that here σ refers to a spin-(s + 1) field, rather than spin-s. We compute the nuclear recoil scattering cross section in App. A. We obtain, where p and p and k and k refer to the ingoing and outgoing momenta of the SM nucleus and DM, respectively, v is the relative velocity of the incoming DM particle, and P s are the Legendre polynomials. The details of the computation of the cross section can be found in App. A. The quantity relevant to direct directional detection is the double differential rate of nuclear recoil events [105,106]. This is given in standard notation by (32) Here, v is the relative velocity of the incoming DM particle relative to the target nucleus, κ DM is related to the local halo density of DM near the earth (ρ DM 0.3 GeV/cm 3 ), f is the velocity distribution of DM in the galactic halo, which is dependent on v as well as v ⊕ , the time dependent Earth velocity in the galactic rest frame. We also define w as a unit vector pointing in the direction of nuclear recoil, w q = q/(2µ). We will follow the methods in [106], outlined below, to obtain an expression for our particular case.
Following [106], we assume a Maxwell-Boltzmann velocity distribution f (v), truncated at the galactic escape velocity, which we take to be v esc = 544 km/s. The most probable speed is v 0 = 220 km/s, which is the circular speed of the local standard of rest. This can be written as where Here v e is the Earth velocity in the galactic rest frame, which is generally a function of time, but for simplicity can be taken to be its constant magnitude, v e = 232 km/s. The overall normalization N is a constant given by We can then break up the velocity integral into two parts: For simplicity here we focus on the first term; the resulting spin and angular dependence is the same in both cases. We make use of the identity, wherev is a unit vector in the direction of the incoming dark matter velocity,q is a unit vector in the recoil direction, and µ is the reduced mass of nucleon-DM system, µ = m N mDM m N +mDM . Then, performing the integral over the Dirac delta function we obtain, where we have taken the dark matter mass to be m. Now let us define the angles α, β, θ, and φ, such that q= (sinαcosβ, sinαsinβ, cosα), (40) v= (sinθcosφ, sinθsinφ, cosθ) wherev is the direction of the incoming DM particle,q is the direction of nuclear recoil,v is the direction of the outgoing HS DM particle, andp is the unit vector in the direction of the incoming nucleus, which we will take to be in theẑ direction. One can deduce v from conservation of momentum. Substituting into Eq. (39) the cross section Eq. (31), we find the double differential recoil rate, where dΩ is the integration over incoming momenta, and f SI (v) is the spin-independent contribution to the scattering rate, given by, Finally, in the limit that the dark matter mass is larger than the momentum transfer, i.e., m q, the term proportional to mv in the second line will be dominant over the first.
The final result is then given by, where the prefactor N is defined as, Here, note the distinct dependence on the Legendre polynomials, P s . The angular dependence in f SI (v) is spinindependent and is present in conventional WIMP models.
The angular dependence on the Legendre polynomials, one the other hand, is a direct consequence of considering higher spins. This generalizes previous works (see e.g. [105][106][107][108]) that considered the imprints of spin-0, 1/2, and 1 DM in the double differential recoil rate,leading to ring-like features. Our HSDM model modulates the bosonic signal further, with an added angular dependence on Legendre polynomials, P s , due to scattering of a spins + 1 boson. This avenue for direct detection is complementary to work that has been done within the 'cosmological collider physics' program, which predicts that higher spin particles produced during inflation will leave behind a distinctive signature in the cosmic microwave background, proportional to Legendre polynomials [25,37,39]. An implication of our results is the potential for a dual 'smoking gun' signature and relationship between CMB experiments and dark matter direct detection experiments. We see that the angular dependence on Legendre polynomials in the CMB non-gaussianity is mirrored by the appearance of Legendre polynomials in the double differential recoil rate.

V. OTHER OBSERVABLE WINDOWS
Aside from direct detection, one might wonder what signature higher spin dark matter may leave at the other pillars of dark matter detection: collider production, and indirect detection (e.g., at the galactic center [109]). To this end, we consider the other possible interactions that could couple higher spin dark matter candidate to the Standard Model.
The simplest possibility is to directly couple a higher spin boson to the Standard Model Higgs boson. For example, an interaction of the form, where ∂ (s) denotes s number of derivatives. From this interaction, a high energy Higgs boson, e.g. generated at a collider, could radiate a spin-s boson σ. At quadratic order in σ, one could have a σ 2 |H| 2 interaction, allowing 2 ↔ 2 scattering of the Higgs and HS boson. The field content of the standard model is not restricted to bosons, and nor are higher spin field theories. Higher spin fermions are interesting in their own right, and may serve as their own dark matter candidate. Unfortunately, a relic density computation as has been done here is not readily repeated for fermions, since the exact solution of HS fermions in dS space is not known. Nonetheless, one may expect a higher spin fermon dark matter candidate to couple to the standard model fermions via 4-fermion interactions, e.g., of the form, where Ψ is a spin-s + 1/2 fermion, and f is a standard model fermion. Such an interaction allows for signatures at precision electroweak experiments, and via annihilation of HS fermions into standard model fermions, a signature in the galactic center [109]. Finally, there is the possibility that the higher spin bosons and higher spin fermions could be organized into multiplets, that is, into super-multiplets of a supersymmetric theory of higher spin fields [42][43][44] (for recent work, see e.g. [54-57, 59, 61-64, 66-69]). Supersymmetry constrains both the spectrum of the theory and the interactions, as discussed in a cosmological context in [39]. Even if supersymmetry is broken at a high scale, one would expect some remnant of this structure to remain at energies accessible by terrestrial experiments.
The cosmological collider analysis of higher spin supersymmetry [39] revealed correlated signals, with the usual P s (cos θ) angular dependence accompanied by superpartner contributions that scale as P s+1 (cos θ) and m P m s (cos θ). It will be interesting to compute the dual signal in directional direct detection.

VI. DISCUSSION
In this work we have considered gravitationally produced massive higher spin particles as a model of dark matter. We have shown that there is a wide range of parameter space for superheavy particles with s > 2 for which the correct relic density of dark matter is produced. We have also explored a potential directional direct detection signature, showing that there is distinctive spin dependent angular dependence in the double differential recoil rate. This enters in the form of Legendre polynomials and is complementary to the 'cosmological collider' signature in the cosmic microwave background.
This opens up opportunities to explore a wide range of new models and parameter spaces to aid in the search for dark matter. We have mentioned several possibilities, but there is certainly much that is still unknown. In future work we will perform a more rigorous numerical exploration of our results in the context of directional direct detection, in order to show more explicitly the impact of the higher spin angular dependence in the double differential recoil rate. It would also be of interest to build a similar model for fermionic HSDM and build connections to HS supersymmetry. We leave these explorations and others to future work.
Comment added: During the final preparation of this manuscript we became aware of recent work [110] with thematic overlap to this paper. The effective field theory approach taken there is promising and may yield further directional direct detection signatures in addition to that considered here. There is no overlap in the content of the papers.

ACKNOWLEDGMENTS
We thank Sylvester James Gates, Jr., Rocky Kolb, Konstantinos Koutrolikos, and Andrew Long, for useful discussions.

Appendix A: Matrix Element Calculation
We here show technical details of the differential cross section used in Section IV. To calculate the differential cross section for a HSDM particle scattering off a SM fermion, there are several subtleties that must be considered. The vertex factor for scattering with σ 0 s,s is given by where we have considered all possible combinations of derivatives onψψ, and k 1 and k 2 correspond to the momenta of the standard model nucleons. We assume that the incoming nucleon is at rest, therefore the only terms that contributes to the vertex will be those which have either k s 1 or k s 2 . The helicity state of the higher spin particle is λ. Thus, this expression can be simplified as Each higher spin external leg carries a factor of the spins + 1 polarization tensor [λ ]i1...is µ (k 3 ). This can be decomposed into a spin-1 component and a spin-s component as follows: where λ = −1, 0, 1 and λ = −s...s are the possible helicity states of the spin-1 and spin-s components, respectively. Lastly, note that generally, working in an expanding background will lead to additional factors of the scale factor, a(t). However, for the remainder of the calculation we normalize a(t) = 1, to account for the insensitivity of particle physics experiments today to the previous expansion of the universe. The matrix element can easily be found in analogy with standard QED computations, and with the use of the relation q i1 ...q is λ i1...is ≡ E λ λ (θ, φ)P λ s (cos θ), where cos θ =q ·k, cos φ =q · and E λ λ and P λ s are the transverse and longitudinal parts of the spherical harmonics, respectively [37,38], where for the λ = 0 modes we simply haveq i1 ...q is λ i1...is = P s (cos θ). Then, we find the matrix element to be |M| 2 = g 4 s Λ 4s p s P s (k ·p ) + (p + k) s P s (k ·(p + k)) 2 · p s P s (k ·p) + (p + k) s P s (k ·(p + k)) 2 · (2m) 2 (2p · k) 2 (16p 2 + 64p · k + 32k 2 ), where k and k refer to the incoming and outgoing HSDM particle, respectively,p and p refer to the standard model nucleus, and we note(p + k) ≡ ( p + k)/| p + k|. From the matrix element, one can find the differential scattering cross section given by dσ dE R = 2m πv 2 1 (2J + 1)(2s χ + 1) We will simplify this by noting that in our construction, we assume that the standard model nucleus is stationary and the incoming HS particle has a much larger momentum, k p. In this limit, we can also take k · (p + k) = k · k. Thus, keeping only the relevant dominant terms in k, we obtain for the cross section dσ dE R m 3 πv 2 1 (2s + 1) g 4 s Λ 4s p s P s (k ·p ) + k s P s (k ·k) 2 · p s P s (k ·p) 2 · 32k 2 (p · k) 2 . (A7)