The phenomenology of nuclear scattering for a WIMP of arbitrary spin

We provide a first systematic and quantitative discussion of the phenomenology of the non-relativistic effective Hamiltonian describing the nuclear scattering process for a Weakly Interacting Massive Particle (WIMP) of arbitrary spin $j_\chi$. To this aim we obtain constraints from a representative sample of present direct detection experiments assuming the WIMP-nucleus scattering process to be driven by each one of the 44 effective couplings that arise for $j_\chi\le$ 2. We find that a high value of the multipolarity $s\le 2 j_\chi$ of the coupling, related to the power of the momentum transfer $q$ appearing in the scattering amplitude, leads to a suppression of the expected rates and pushes the expected differential spectra to large recoil energies $E_R$. For $s\le$ 4 the effective scales probed by direct detection experiments can be suppressed by up to 5 orders of magnitude compared to the case of a standard spin-independent interaction. For operators with large $s$ the expected differential spectra can be pushed to recoil energies in the MeV range, with the largest part of the signal concentrated at $E_R\gtrsim$ 100 keV and a peculiar structure of peaks and minima arising when both the nuclear target and the WIMP are heavy. As a consequence the present bounds on the effective operators can be significantly improved by extending the recoil energy intervals to higher recoil energies. Our analysis assumes effective interaction operators that are irreducible under the rotation group. Such operators drive the interactions of high-multipole dark matter candidates, i.e. states that possess only the highest multipole allowed by their spin. As a consequence our analysis represents also the first phenomenological study of the direct detection of quadrupolar, octupolar, and hexadecapolar dark matter.


Introduction
Weakly Interacting Massive Particles (WIMPs) are the most popular particle candidates to provide the invisible halos of Galaxies, including the Milky Way. WIMPs are expected to have only weak-type interactions and to have a mass m χ falling in the GeV-TeV range. Direct detection (DD) experiments look for the interaction of WIMPs with ordinary matter through the measurement of the recoil energy E R < ∼ 10 keV -100 keV imparted by WIMPs when they scatter elastically off the nuclear targets of a wide range of solid-state and liquid low-background underground detectors.
Model-independent approaches have become increasingly popular to interpret Dark Matter (DM) search experiments  due to the growing tension between the WIMPs arising in popular extensions of the Standard Model (SM) such as Supersymmetry or Large Extra Dimensions and the constraints from the Large Hadron Collider (LHC). In particular, since the DD process is non-relativistic (NR) the WIMP-nucleon interaction can be parameterized with an effective Hamiltonian H that complies with Galilean symmetry.
The effective Hamiltonian H to zero-th order in the WIMP-nucleon relative velocity v and momentum transfer q has been known since at least Ref. [26], and consists of the usual spin-dependent (SD) and spin-independent (SI) terms. To first order in v, the effective Hamiltonian H has been systematically described in [27,28] for WIMPs of spin 0 and 1/2, and less systematically described in [29,30] for WIMPs of spin 1 and in [31] for WIMPs of spin 3/2.
Recently in [32] the NR effective Hamiltonian for WIMP-nucleous scattering has been extended to include WIMPs of arbitrary spin j χ in the approximation of one-nucleon currents. In Ref. [32] H is written in a complete basis of rotationally invariant operators organized according to the rank of the 2j χ + 1 irreducible operator products of up to 2j χ WIMP spin vectors. In particular, for a WIMP of spin j χ a basis of 4 + 20j χ independent operators arise that can be matched to any high-energy model of particle dark matter, including elementary particles and composite states.
So far, only 20 of such operators for a WIMP with j χ ≤1 have been considered in the literature [27][28][29][30]. As a consequence, the introduction in [32] of so many new interaction terms for j χ ≥1 has potentially interesting phenomenological consequences.
In particular, for each new interaction term a different scaling law for the WIMP cross section, needed to compare the results of experiments using different target nuclei, is expected to arise. Moreover, the new high-rank operators that arise at high values of the spin j χ depend on increasing powers of the transferred momentum q ≡ | q|. On the one hand such increasing momentum suppression implies lower expected rates, so that present experimental sensitivities are expected to probe increasingly low values of the effective energy scale of the effective theory. Moreover, for a given choice of the WIMP velocity distribution and m χ > ∼ 100 GeV the increasing powers of q shift the expected rate spectra to high recoil energies E R ≃1 MeV, much higher than those usually expected.
The goal of the present paper is to provide a first systematic and quantitative discussion of the phenomenological aspects outlined above for the effective interaction operators introduced in [32] . To this aim, we will consider each of the 44 operators that can arise in the nuclear scattering of a WIMP with j χ ≤2. Due to the large dimensionality of the parameter space of high-spin DM, starting with analyzing the phenomenology of one operator of such basis at a time appears like a sensible approach. In Section 2.2 we elaborate on the possibility that such basis may have indeed a connection to the physical world, and that DM candidates whose phenomenology is driven by such operators do exist.
Throughout the paper we will assume a standard Maxwellian in the Galactic rest frame cut at the WIMP escape velocity u esc =550 km/s and with reference values for the other relevant parameters: ρ χ =0.3 GeV/cm 3 for the number density of the WIMPs in the neighborhood of the Sun and v rms =270 km/s for the WIMP root-mean-square velocity.
The paper is organized as follows. In Section 2 we discuss the theoretical motivation of DM candidates of arbitrary spin (Section 2.1) and high multipole (Section 2.2). In Section 3 we outline the results of Ref. [32], where the non-relativistic effective operators for j χ ≥1 are introduced and the analytic expressions of the relevant expected rates are calculated. Section 4 is devoted to our phenomenological discussion: in Section 4.1 we summarize the expressions for the scattering rate; Section 4.2 provides a discussion of the values of the energy scale of the effective theory probed by a representative selection of present directdetection experiments; in Section 4.3 we discuss the expected differential rate, showing that when the scattering process is driven by a high-rank (high-spin) operator and both the target and the WIMP are heavy the largest part of the signal is concentrated at E R > ∼ 100 keV with a characteristic pattern of peaks and minima at high recoil energies; in Section 4.4 we show how the present bounds can be significantly improved by extending the recoil energy intervals to higher recoil energies. Finally, we provide our Conclusions in Section 5. The details of the procedure followed to obtain the upper bounds is described in the Appendix.

High-spin and multipolar dark matter
In the present Section we discuss the theoretical motivation for DM candidates of arbitrarily high spin. Moreover, we wish to elaborate on the possibility that it possesses only the highest anomalous moment allowed for its spin, implying that its interaction with ordinary matter is driven by the high-rank irreducible operators introduced in Eq. (3.10) and whose phenomenology is the subject of our analysis.

High-spin DM
While most (but not all) of the theoretical and experimental work on detection of particle dark matter has been focused on dark matter particles that are elementary and have spin 0 or 1/2, there is no compelling reason for dark matter particles to be elementary, or for their spin to be limited to 0 and 1/2. In fact, as is well-known, non-elementary particles and particles of spin higher than 1/2 exist in nature. We include here not only the states of particle physics, such as protons, neutrons, hadrons, gauge bosons, etc., but also larger composite objects like atomic nuclei, atoms, molecules, etc. In particular, the possibility to distinguish between pointlike and non-pointlike dark-matter candidates in direct detection searches through the shape of the nuclear recoil energy spectrum has been discussed in the literature [33].
Sometimes one hears that particles with high spin cannot have gauge interactions, but this is an incorrect generalization of two quite specific statements in the construction of interacting theories: there are difficulties in coupling a massless [34] or massive particle [35] of spin higher than 2 to gravity and there are difficulties in coupling a massless particle to the photon [34]. On the other hand there are no difficulties in coupling any particle to photons at the effective lagrangian level, i.e. in a theory valid at energy scales smaller than the binding energy of the particle. Witness to this is the coupling of complex atoms and molecules, and even macroscopic objects, to the electromagnetic field. Even theoretically, consistent effective Lagrangians coupling towers of particles of higher and higher spin to gravity have been constructed as low energy limits of string theory (see e.g. [36]). And consistent effective couplings of particles to the electromagnetic field have been obtained as four-dimensional reductions of field theories in spacetimes with extra dimensions (see, e.g, [37][38][39]). The non-local character of the fundamental theories involved bypasses the limits of local theories with a single particle of high spin.
Composite dark matter has been introduced for example in technicolor theories (where a new strongly interacting sector is postulated) [40,41] or as another example in mirror dark matter models (where dark copies of the standard model particles are assumed, complete with dark protons, dark neutrons, dark atomic nuclei, and so on; see for example [42] and the review [43]). The continuing interest in composite dark matter can be gathered for example from the references in [44]. In most of these scenarios, the dark matter particles have spin 1/2 and interact with the standard model through kinetic mixing of the dark and visible photon, or the dark and visible Higgs bosons, or the dark and visible neutrinos. Models of this kind have also been considered in the so-called asymmetric dark matter scenarios [45]. Most but not all mirror dark matter candidates are elementary particles of spin 0 or 1/2. An exception we found in the literature is the pangenesis model of [46] in which dark matter is atomic, being a mirror hydrogen atom composed of a mirror proton and a mirror electron.
There are stable light atomic nuclei with spin 3/2 ( 7 Li, 11 Be, 11 B) and spin 3 ( 10 B), and heavier stable nuclei can have spins up to 9/2 (e.g., 73 Ge or 87 Sr). While it may not be easy to find a theoretical mechanism in which heavy dark nuclei would be more abundant that light ones, so little is known about the dynamics in the dark sector that it may well be possible that only dark nuclei with high spin (1, 3/2, 2, ...) interact with ordinary nuclei and are in principle detectable in direct dark matter search experiments. For example, dark heavy elements such as O, Ne, N, C, and Fe (besides H and He) have been considered in the context of direct dark matter detection [47] with a long range interaction mediated by a kinetically-mixed photon-dark photon coupled to the electric charge of the dark nuclei.
Among the particles of spin 1, it is worth mentioning the deuteron and the massive gauge bosons W and Z. The deuteron, which is stable, was produced in the early universe during primordial nucleosynthesis. Its primordial abundance was delicately determined by the particular values of binding energies, lifetimes, and reaction cross sections of hydrogen, deuterium, helium, and the neutron. An analog of the deuteron in the dark sector may have a completely different, and perhaps much larger abundance, thanks to a more favorable dynamics in the dark sector; in fact, in models with mirror dark matter, the dark primordial nucleosynthesis is set to happen at a slightly lower temperature than in the SM sector, to avoid issues with the number of neutrino species and other measures of the abundance of relativistic particles in the early universe, and the lower nucleosynthesis temperature produces a different pattern of primordial abundances, with 4 He being dominant over the other species [48].
One can also imagine a "quasi-mirror" sector that is almost a copy of the SM model but has small differences in the values of some parameters (models with broken mirror symmetry can be found for example in [49][50][51][52]. In a quasi-mirror dark sector the dark neutron may be lighter than the dark proton. Indeed, the experimental fact that the mass of the proton is slightly smaller than the mass of the neutron is not well understood theoretically, as the proton-neutron mass difference is moved one level lower into a mass difference between up and down quarks, which is obtained phenomenologically by fits to experimental data and for which there is no theoretical calculation from first principles. Thus it may well be that a slightly different dynamics in the dark sector produces a dark neutron lighter than the dark proton. Then a free dark proton may decay into a dark neutron, a dark positron, and a dark neutrino, but a dark neutron bound in a dark nucleus may not decay, and the stability of atomic nuclei would be completely changed. This would open the possibility that dark neutrons form dark multineutron states. There is in principle no obstacle for these dark multineutrons to have a large value of their spin. In addition, given enough symmetry, many of the lower electric and magnetic multipoles of the dark multineutron may vanish, making the dominant dark multineutron-photon interaction a multipole of high order in the dark photon-photon kinetic mixing scenario.

High-multipole DM
In Section 4 we discuss the phenomenology of high-spin DM states assuming that their interaction with ordinary matter is driven by one of the irreducible operators introduced in [32] and given in Eq. (3.10). In particular, for a WIMP of spin j χ in Section 4 the phenomenology of the operators of maximal rank s = 2 j χ is shown (although such results can be easily rescaled to lower s values using the coefficient plotted in Fig. 1). As discussed in Ref. [32] the use of the operators of Eq. (3.10) as a basis for the effective Hamiltonian of the WIMP-nucleus interaction has several advantages: for instance, it is complete, it avoids double counting and simplifies the calculation of the cross section, that results in a sum of cleanly separated contributions from irreducible operators of different ranks that do not interfere. From this point of view, however, to use such specific basis of operators instead of any other appears to be a mere technicality, and one may argue that the nonrelativistic limit of a generic Hamiltonian should naturally contain a mix of operators of all ranks, with the contribution of the highest rank operators unavoidably suppressed. In this Section we wish to elaborate on the possibility that, instead, such a basis may have indeed a connection to the physical world, i.e. it is possible to conceive DM candidates whose interaction with ordinary matter is driven by the highest multipole moments connected to the high rank operators whose phenomenology is the subject of our analysis.
There is a general relation between the spin of the particle and the highest nonzero multipole moment it can possess. The Wigner-Eckart theorem and analogous group theoretic arguments in special relativity, allow a particle of spin-j χ to have moments up to order 2j χ + 1: a monopole for spin-0, a dipole for spin-1/2, etc. As a real world example, the neutron has zero net electric charge, nonzero magnetic moment, and all of its higher moments vanish.
Molecules in the dark sector provide an example of both compositeness and high multipole moments. Dipolar molecules (having zero net charge and non-zero permanent electric dipole moment) are well known, while quadrupolar, octupolar, and hexadecapolar molecules may be less known. Which multiple moment a molecule possesses is determined by the kind of symmetry its distribution of charges has. Examples of quadrupolar molecules, which have zero charge and zero dipole moment but non-zero quadrupole moment, are the linear molecules of carbon dioxide CO 2 , CS 2 and the planar molecules of C 6 H 6 . Among the octupolar molecules, whose first nonzero multipole is the octupole, are molecules with tetrahedral symmetry like CH 4 and CF 4 . Hexadecapole molecules, which possess a non-zero hexadecapole but have neither a dipole, quadrupole nor octupole moment, include the octahedral symmetry molecules of SF 6 and UF 6 .
It is possible to envisage a dynamics in the dark sector that may lead to dark molecules of similarly high symmetry. It is important to notice that for what concern the detection of dark matter molecules, it is the charge distribution of the charges coupled to ordinary standard model particles that matters, while the formation and stability of dark molecules may be governed by forces that act within the dark sector exclusively. In addition, one needs to consider the phenomenon of induced polarization. In the case of multipolar dark molecules in the vicinity of ordinary matter one may worry that induced dipole (or other low) moments would come to dominate the interaction at low energies. This may not happen either because of a fortuitous arrangement of charges giving a small polarizability, or more generally because the interaction range of the force between dark and ordinary matter may be so small that there is no time for induced multipoles to develop during the scattering of dark matter molecules with ordinary matter.
The idea that dark matter may carry only multipoles of high order generates the possibility that dark matter may be a neutral particle that possesses only the highest anomalous moment allowed for its spin. The highest multiple moment for a particle of spin j χ is the 2 2jχ multipole, e.g., the dipole for spin-1/2, the quadrupole for spin-1, the octupole for spin-3/2, and the hexadecapole for spin-2. The anomalous moment could be an electric or magnetic moment (examples from the literature are recalled later in this section) or it could be a multipole moment in the interaction of dark matter with new gauge fields in the dark sector, Abelian or non-Abelian. One can thus have "magnetic-dipole dark matter" or "electric-octupole dark matter" or even higher multipole dark matter according to which multipole moment is the lowest non-vanishing moment.
For electromagnetic interactions and multipoles, effective Lagrangians involve the electromagnetic field strength F µν multiplied by an increasing number of derivatives of the particle fields, or of the momentum exchange in the corresponding amplitude, multiplied by an appropriate power of the particle radius or of the inverse of the particle mass. It is important to stress that it is the particle mass, or the particle radius, that accompany the nonrenormalizable operators in a multipolar effective Lagrangian, and not an energy cutoff scale at which the interactions become weak or strong, or at which new degrees of freedom become dynamical. This is exemplified by the Fermi weak theory on one hand and hadrons on the other, for which in the first case the mediator mass is (much) larger than the energy scale involved in the weak processes, and in the second case the scale of strong interactions is lower than the mass of the hadrons.
In absence of CP violation the highest anomalous multipole alternates between electric and magnetic as the spin of the particle increases: a magnetic dipole for spin 1/2, an electric quadrupole for spin 1, a magnetic octupole for spin 3/2, an electric hexadecapole for spin 2, and so on. If F µν is replaced by its dualF µν = ǫ µνλρ F λρ /2, the multipoles in the series switch character: an electric dipole for spin 1/2, a magnetic quadrupole for spin 1, an electric octupole for spin 3/2, a magnetic hexadecapole for spin 2, and so on.
For spin-1/2, a particle can have an electric and a magnetic dipole, the latter with a contribution from the charge (the Dirac magnetic moment) to which an extra contribution may be added (the anomalous magnetic moment). The classic example of a neutral spin-1/2 particle with nonzero anomalous magnetic moment is the neutron. A dark neutron analog for dark matter was suggested already in [53,54], where dipolar dark matter was introduced. In that case, the dark matter was considered to interact with the usual electromagnetic field and have zero electric charge but non-zero anomalous magnetic moment, or non-zero electric dipole moment in case parity-violating interactions were allowed. The interaction lagrangian of a magnetic-dipole dark matter particle χ, of spin 1/2 and zero electric charge, reads where F µν is the electromagnetic field strength tensor, σ µν is the usual combination of Dirac γ matrices, e the elementary unit charge (not the particle charge, which is set to zero), and κ is the anomalous magnetic moment of the particle. The contributions of the anomalous magnetic dipole term in Eq. (2.1) can be distinguished in the three-particle χ(p 1 )-χ(p 2 )-A µ (q) vertex (momentum p 1 incoming, momenta p 2 and q = p 1 − p 2 outgoing) by being those terms that contain the highest allowed power of the momentum q µ carried by the photon ieκ m 2 The same scenario can be set-up for particles of higher spin by keeping only their highest possible multipole. The highest-multipole terms in the three-particle χ-χ-γ vertex can be distinguished as being the terms that carry the 2j χ -th power (the highest power) of the photon momentum q µ (the completely symmetric homogenous combination of 2j χ vectors q µ ).
For spin-1, the W boson has a measured magnetic dipole and electric quadrupole moments, in agreement with the gauge symmetry of the Standard Model [55]. Its magnetic dipole moment is usually parametrized as µ W = e(1 + κ + λ)/2m W , and its electric quadrupole moment as Q W = −e(κ − λ)/m 2 W , where the parameters κ and λ appear in the effective Lagrangian [56]: In the Standard Model, κ = 1 and λ = 0. The matrix element for the elastic scattering of spin-1 quadrupolar dark matter by one-photon exchange can be borrowed from the analogous matrix element for a deuteron [57]: where j µ is the electromagnetic current, p µ 1 , ǫ µ 1 and p µ 2 , ǫ µ 2 are the incoming and outgoing momenta and polarization vectors (of polarization λ 1 and λ 2 ) of the dark matter particle, m χ is the DM particle mass, and q µ = p µ 1 − p µ 2 is the momentum transfer. Perturbative expressions are known explicitly: for a quadrupolar particle χ µ the Feynman rule for the [58]; the momentum p 1 is incoming, the momenta p 2 and q = p 1 − p 2 are outgoing) and as a lagrangian term [59] Here λ is the anomalous quadrupole moment of the spin-1 particle.
In principle, the above considerations on anomalous multipole dark matter are not limited to the exchange of a massless vector mediator. A massive mediator would present the same structure for the interaction lagrangian and the interaction vertex, and would provide short range interactions. A scalar or fermion mediator, whether massive or massless, would also be subjected to a multipole expansion, as the multipole expansion is basically an expansion of a potential in spherical harmonics, or essentially an expansion in (symmetric traceless combinations of) derivatives of the mediator field.
We are not aware of any study of multipolar dark matter beyond the dipole. This paper provides the first phenomenological study of the direct detection of quadrupolar, octupolar, and hexadecapolar dark matter.

Effective theory of nuclear scattering for a WIMP of arbitrary spin
In Ref. [32] a systematic approach is introduced to characterize the most general nonrelativistic WIMP-nucleus interaction allowed by Galilean invariance for a WIMP of arbitrary spin j χ in the approximation of one-nucleon currents, i.e. assuming that the WIMP interacts with one nucleon at at time. The procedure consists in organizing the WIMP currents according to the rank of the 2j χ + 1 irreducible operator products of the up to the 2j χ WIMP spin vectors S χ required to mediate transitions where the third component of the WIMP spin changes from ±j χ to ∓j χ . Using index notation S i for the i-th component of the vector S χ (and dropping the subscript χ in S χ,i for more readability), these interaction terms contain no S i or a product of s factors S i up to s = 2j χ , The calculation of the cross section for WIMP-nucleus scattering is greatly simplified if for the products of WIMP spin operators one uses irreducible tensors (i.e., belonging to irreducible representations of the rotation group). Irreducible tensors are completely symmetric under exchange of any two of their indices and have zero trace under contraction of any number of pairs of indices (they are symmetric traceless tensors). In particular irreducible tensor operators of different rank are independent, in the sense that the trace of their product is zero, so that there are no interference terms in the cross section between irreducible operators of different spin. Therefore, the products of Eq. (3.1) are conveniently substituted by the following 2j χ + 1 irreducible spin tensors: where we use an overbracket over an expression containing a set of indices to indicate that the free indices under the bracket are completely symmetrized and all of their contractions are subtracted.
Besides being a convenient basis to describe the WIMP current in spin space the symmetric traceless products of spin operators of Eq. (3.2) can also be seen as the interaction terms appearing in the multipole expansion of the WIMP-nucleus effective potential. For instance, given a scalar potential V (r χN ) (with r χN the WIMP-nucleus distance) the symmetric traceless combinations of l of its derivatives ∂ i 1 ∂ i 2 · · · ∂ is V represents the l-th multipole in its multipole expansion, and when coupled to a generic combination of WIMP spin operators singles out the one with the highest multipole: Analogous expansions in derivatives of vector potentials exist [32]. So operators that are irreducible tensors under the rotation group can drive the effective interaction of the high-multipole DM candidates discussed in Section 2.2.
Once the number of WIMP spin factors is fixed to s, one couples the WIMP currents to one of the five nucleon currents that arise in WIMP-nucleon scattering from the nonrelativistic limit of the free nucleon Dirac bilinears ψ f Γψ i (with Γ any combination of Dirac γ matrices and ψ the Dirac spinor for a relativistic free nucleon) and the most general interaction Hamiltonian depending on the WIMP spin operator and the nucleon spin operator: In the equation above σ N is the vector of Pauli spin matrices acting on the spin states of the nucleon N , while For each operator O X (X=1, Σ, ∆, Φ, Ω) the number of q ≡ q/m N factors is constrained by rotational invariance. In particular, in the case of a scalar nucleon operator O X (X = M, Ω), the WIMP operator o must be a scalar, and all the indices i 1 · · · i s in S i 1 · · · S is must be saturated by terms q i 1 · · · q is . The resulting WIMP operator is S i 1 · · · S is q i 1 · · · q is . On the other hand, in the case of a vector nucleon operator O X (X = Σ, ∆, Φ), a vector WIMP operator o is needed, and the s indices in S i 1 · · · S is must be saturated by an appropriate number of q factors in order to obtain a vector. This can be achieved in three ways: (1) by using s − 1 factors of q to produce S i 1 · · · S is q i 1 · · · q i s−1 with free index i s , (2) by using s factors of q to produce ǫ islm S i 1 · · · S i s−1 S l q i 1 · · · q i s−1 q m , again with free index i s , and (3) by using s + 1 factors of q to produce S i 1 · · · S is q i 1 · · · q is q i s+1 , with free index i s+1 . Using the irreducible spin products in Eq. (3.2) in place of those in Eq. (3.1), one introduces the scalar WIMP operators and the vector WIMP operators so that the following basis of WIMP-nucleon operators O X,s,l , all of which are irreducible in WIMP spin space and Hermitian, is obtained: In terms of the Wilson coefficients c τ X,s,l (q) the WIMP-nucleon effective Hamiltonian H is a linear combination of the WIMP-nucleon operators listed above: with τ an isospin index (0 for isoscalar and 1 for isovector) and t 0 = 1, t 1 = τ 3 are nucleon isospin operators (the 2×2 identity and the third Pauli matrix, respectively). The relations between the isoscalar and isovector coupling constants c 0 X,s,l and c 1 X,s,l and the proton and neutron coupling constants c p X,s,l and c n X,s,l is: In Table 1 we provide a dictionary between the operator basis used in the literature for a WIMP of spin j χ ≤1 [27][28][29][30] and the basis introduced in [32] and summarized in Eq. (3.10). Table 1. Non-relativistic Galilean invariant operators discussed in the literature ( [28][29][30]) for a dark matter particle of spin 0, 1/2 and 1, and their relation with the WIMP-nucleon operators O X,s,l defined in Eqs. (3.10). Notice that the sign convention for the momentum transfer q used in this table and throughout the paper is opposite to that of Refs. [28][29][30].
The unpolarized differential cross section for WIMP-nucleus scattering is given by the expression: where v χT ≡ | v χT | is the WIMP speed in the reference frame of the nucleus center of mass and: with µ χT the reduced WIMP-nucleus mass. The sum in Eq. (3.13) is over 15) and the nuclear response functions W T X are available in the literature for the most common nuclear targets used in direct detection experiments [28,60]. Notice that they do not depend on j χ thanks to the factorization between WIMP and nuclear currents valid in one-nucleon approximation. On the other hand the WIMP response functions R τ τ ′ X are given by [32]: coupling Table 2. The values of X for the nuclear response functions W τ τ ′ T X corresponding to each of the couplings of Eq. (3.10), for the velocity-independent and the velocity-dependent components parts of the WIMP response function, decomposed as in Eq. (3.19). In parenthesis the power of q in the cross section. where: For a WIMP of spin j χ all the operators O X,s,l with s ≤ 2j χ contribute to the cross section. Conversely, for a given operator O X,s,l the dependence on j χ of the squared amplitude is encoded in the B jχ,s factors. A plot of B jχ,s is provided in Fig. 1. Writing: the correspondence between each of the couplings of Eq. (3.10) and the values of X in the nuclear response functions W τ τ ′ T X is provided in Table 2.
4 Experimental sensitivities to the effective operators

The scattering rate
In this Section we will assume that the WIMP-nucleus scattering process is driven by each one of the 44 effective couplings in Eq. (3.10) for j χ ≤ 2, and obtain constraints from a representative sample of the present DD experiments by comparing the expected rate in each model to the corresponding experimental upper bound. In particular the expected number of events in a WIMP direct detection experiment in the interval of visible energy E ′ 1 ≤ E ′ ≤ E ′ 2 is given by: where T indicates the target nuclei present in the detector. For a given target T : In the equations above E R is the recoil energy deposited in the scattering process (indicated in keVnr), while E ee =q T (E R )E R (indicated in keVee) is the fraction of E R that goes into ionization and scintillation (wth q T (E R ) the target quenching factor) and G(E ′ , E ee ) is the effect of the energy resolution (so that E ′ is measured instead of E ee through a calibration procedure). Finally ǫ(E ′ ) is the measured experimental acceptance. The theoretical differential rate is given by: 3) with f ( v, t) the WIMP velocity distribution in the Earth's rest frame, while: with m T the nuclear target and µ χT the WIMP-target reduced mass.

Present constraints
For each of the 44 couplings we parameterize the corresponding Wilson coefficient as: where g represents an effective coupling constant while M a mediator mass. In the central plot of Fig. 2 we show the upper bound on the coupling g as a function of M for the effective operator O Φ,4,5 , m χ =1 TeV and j χ =2. In our analysis we consider 5 representative direct detection experiments (XENON1T [61], XENON100 [62], SuperCDMS [63], PICO-60 [64] and COSINE-100 [65] that use 4 different targets: Xe, Ge, C 3 F 8 and N aI (see Appendix A for a summary of our procedure). In particular, in Fig. 2 the bound is driven by XENON100, which provides the stronger constraint. The central plot of Fig. 2 shows a clear transition from the regime M ≪ q, for which c τ Φ,4,5 ≃ g 2 /q 2 (long-range interaction) to regime for which M ≫ q and c τ Φ,4,5 ≃ g 2 /M 2 (contact-interaction). In particular the square (green) markers, that show the result of an accurate evaluation of the bound using the full expression (4.5), follow closely the solid line that represents the curve g 2 /(M 2 + q 2 0 )=K, where the two parameters q 0 and K are fixed so that the bound for M ≪ q coincides to the constraint g < g lim shown in the left plot for a long-range interaction (i.e. assuming c τ Φ,4,5 ≃ g 2 /q 2 ) and for M ≫ q coincides to the constraint M/g < (M/g) lim shown in the right plot for a short-range interaction (i.e. assuming c τ Φ,4,5 ≃ g 2 /M 2 ). The transition between the two asymptotic regimes, represented in the central plot of Fig. 2 by the horizontal (red) dotted line and by the dot-dashed (black) line corresponds to M = q 0 = g lim (M/g) lim .
In Figs. 3 and 4 we show the result of a systematic analysis on the most constraining lower bound (M/g) lim for m χ =100 GeV and m χ =1 TeV, respectively for all the 44 models and for a contact interaction (c τ X,s,l = g 2 /M 2 ). In each figure the left-hand plot assumes that the WIMP couples to protons only, and the right-hand plot to neutrons only. For each operator O X,s,l the bound assumes j χ =s/2. The constraint for j χ >s/2 can be obtained from that for j χ = s/2 by multiplication times the factor (B s/2,s /B jχ,s ) 1/4 (with B jχ,s given in Eq. (3.17) and plotted in Fig. 1). In all the figures each filled marker represents the most constraining present bound from one among the five experiments that we analyze while, if present, an open marker indicates an estimation of the improvement that could be obtained in the limit by extending the experimental energy ranges beyond the present ones (see Section 4.4). Moreover, for each effective model, a vertical line represents the minimal value of M/g compatible to the assumption of a contact interaction, obtained by combining M > q 0 (following for each model the procedure outlined in Fig. 2) and the perturbativity requirement g 2 /(4π) <1.
A first conclusion one can draw from Figs. 3 and 4 is that for all the couplings the bound on M/g is compatible to the assumption of a contact interaction. Moreover, one observes an anticorrelation between the value of (M/g) lim and the s parameter. This is due to the fact that a larger value of s corresponds to stronger momentum suppression in the expected rate and, as a consequence, to a weaker bound. More specifically, as shown in Table 2 for a given value of s the power of the transferred momentum q ranges from 2s − 2 and 2s + 4.
In Figs   depends on the nuclear response function W τ τ ′ Φ ′′ . Such nuclear response function favors heavier elements with large nuclear shell and scales with the nuclear target similarly to the SI interaction, so that for most isotopes it is the most sizeable among the W τ τ ′ T X with the exception of X = M [27,28] (this holds for all the nuclear targets that we consider with the exception of the semi-magic isotope 72 Ge, for which W τ τ ′ Φ ′′ vanishes). As far as the other operators are concerned, for X= ∆, Σ the constraint gets systematically less stringent at growing l, as one expects due to the enhanced momentum suppression, i.e. (M/g) lim (∆, 1, 0) > (M/g) lim (∆, 1, 1) > (M/g) lim (∆, 1, 2) and (M/g) lim (Σ, 1, 0) > (M/g) lim (Σ, 1, 1) > (M/g) lim (Σ, 1, 2) An exception to this pattern is provided by X = Φ for which one observes instead (M/g) lim (Φ, 1, 0) > (M/g) lim (Φ, 1, 2) > (M/g) lim (Φ, 1, 1). One understands this inversion between (M/g) lim (Φ, 1, 1) and (M/g) lim (Φ, 1, 2) by noticing that, as again shown in Table 2 , as already pointed out. The same pattern among the 10 corresponding constraints on M/g repeats also for s =2, 3,4. Another feature arising from Figs. 3 and 4 is that in most cases the strongest constraint   is provided by xenon detectors (XENON1T and XENON100) with the few exceptions of fluorine (PICO-60) and N aI (COSINE-100), but only for a WIMP-proton interaction. Indeed in most cases xenon experiments are the most competitive due to the large exposure. In order to understand why there are exceptions to this one needs again to consider the mapping provided in Table 2 between each effective coupling and the nuclear response functions W τ τ ′ X . In particular the 10 operators arising at a given value of s can be divided into 4 broad classes: i) operators whose cross section is driven by the W τ τ ′ In particular we notice that the operators O ∆,s,s−1 and O ∆,s,s couple to W τ τ ′ M through a velocity-suppressed term (i.e. R τ τ ′ 1X in Eq. (3.19)). However, in spite of the v 2 T /c 2 ≃ 10 −6 suppression, for heavy nuclei (Xe, I) the scattering amplitude of such operators is dominated by the R τ τ ′ 1X terms thanks to the huge hierarchy between the two response functions W τ τ ′ M ≫ W τ τ ′ ∆ . On the other hand, for lighter targets (Ge, N a, F ) the contributions from R τ τ ′ 0X and R τ τ ′ 1X are of the same order. This feature was already observed in the context of spin-1/2 WIMPs for O 5 =-O ∆,1,1 and O 8 = O ∆,1,0 [24].
The mapping summarized above allows to understand that xenon experiments are the most sensitive to the couplings driven by W τ τ ′ M and W τ τ ′ Φ ′′ , which favor heavy elements. On the other hand, W τ τ ′ Σ ′′ and W τ τ ′ Σ ′ correspond to a coupling of the WIMP to the nucleon spins. Since inside nuclei the nucleon spins tend to cancel each other the contribution from even-numbered nucleons to such response functions is strongly suppressed. As a consequence of this for such interactions neutron-odd targets (such as xenon and germanium) are mostly sensitive to the WIMP-neutron coupling, while proton-odd targets (such as fluorine, sodium and iodine) are mostly sensitive to the WIMP-proton coupling. This implies that in the left-hand plots of Figs. 3 and 4 the rate in xenon detectors is strongly suppressed for X = Σ, Ω. In such cases the strongest constraint is provided by fluorine in PICO-60 which, due to the large exposure, is the most competitive among proton-odd detectors, albeit only for moderate values of s. Indeed, for s = 3, 4 the strongest constraint can instead be provided by xenon or iodine. We postpone the explanation of this fact to Section 4.3, where the expected spectral shape of the differential rate in our scenarios will be discussed. Finally, the response function W τ τ ′ Φ ′ requires a nuclear spin ≥ 1 [27] and vanishes in fluorine, so that also in this case xenon turns out to be the most competitive target.
The results shown in Figs. 2, 3 and 4 are also valid in the case of the multipolar DM discussed in Section 2.2. Notice that the interaction driving the scattering process and the binding force within the DM composite state need not be the same, and that the latter is not accessible experimentally. If, on the other hand, one assumes the same interaction, and that during the scattering process polarization does not induce lower multipoles in the DM state, the nonrenormalizable operators in the multipolar effective Lagrangian depend on the particle mass or the particle radius and not on the energy cutoff scale at which the interactions become weak or strong, or at which new degrees of freedom become dynamical. For instance, the analysis in [66] concludes that the generic energy cutoff for effective theories of electromagnetic interactions of a particle with mass m, spin j and coupling g is ∼ m/g 1/(2j−1) for j ≥ 3/2.
We conclude this Section by pointing out that within the context of the non-relativistic effective theory used in our analysis it is not possible to assess the compatibility of the bounds obtained in Figs. 3 and 4 for the mediator mass M with other constraints, such as those from accelerator physics. Instead, in the next Section we will use such bounds to calculate the expected maximal signals in future direct detection experiments.

Energy spectra for high-spin WIMPs
Besides leading to a suppression of the expected rates, the larger powers of the momentum transfer q = √ 2m T E R that arise in high-multipole operators lead to another important effect: they push the expected differential spectra to larger recoil energies E R . An example of this is shown in Fig 5, where the differential rate of Eq. (4.2) is plotted for each of the 44 couplings fixing m χ =1 TeV, assuming an isoscalar interaction (c p X,s,l = c n X,s,l ) and fixing M/g to the present upper bound, calculated as in Figs. 3 and 4.
Indeed, in some cases the WIMP signal not only extends well beyond E R = 1 MeV (for high multipoles s and for large-enough values of the WIMP mass m χ and of the target mass m T ) but its largest part is concentrated to E R > ∼ 100 keV, showing there a structure of peaks and minima. The peculiar features of such plots can be understood as a combination of three ingredients: the nuclear structure function W τ τ ′ X , the WIMP velocity distribution f ( v) and the power of q that appears in the cross section of each operator (as summarized in Table 2). In particular the peaks appearing in Fig. 5 descend directly from the diffractive nature of the W τ τ ′ T X nuclear structure response functions, which are calculated as Fourier transforms of nuclear current densities. An example of the W τ τ ′ M and W τ τ ′ Φ ′′ functions for xenon as calculated in [28] is shown in Fig. 6 for an isoscalar interaction (τ = τ ′ = 0). In the differential rate the W τ τ ′ T X function is convoluted with the velocity distribution f ( v), which in the standard halo model we adopt is exponentially suppressed when v T,min (E R ) approaches the escape velocity. For a standard SI or SD interaction with no explicit momentum dependence such suppression prevents the peak structure of the W τ τ ′ T X to emerge, and leads to an energy spectrum that falls monotonically with the recoil energy. However, in the general case of the O X,s,l operators when s is large enough the power of q appearing in the cross section enhances the peaks at high recoil energies, so that they become visible, if, at the same energies, the spectrum is not cut by the velocity distribution. For heavy nuclei (Xe, I) this is indeed the case, so that at large s the diffractive peaks become visible (notice that, in the most extreme case, the cross section for O Φ,4,5 is enhanced at high recoil energy by a q 12 ≃ E 6 R factor). This effect is also enhanced for heavier WIMP masses, as illustrated in Fig. 7, where the expected differential rates, calculated with the same assumptions of Fig. 5, are shown for O Φ,4,5 and m χ = 100 GeV, 300 GeV and 1 TeV: in particular as m χ is increased the overall signal rate at high energy gets larger while the spectrum becomes harder, with a growing relative size of the peak with the highest energy. A different situation is observed however for lighter nuclei (F , Ge, N a). In this case, irrespective on m χ , the escape velocity cut on the recoil energy is in the 100-200 keV range, and the diffractive peaks correspond to values of v T beyond the escape velocity. As a consequence, as far as the overall size of the expected signal for O X,s,l is concerned, larger l values (i.e. larger powers of q in the operator) increasingly favour larger-mass targets compared to lighter ones. Another way to see this is that, due to the cut from the escape velocity, the average momentum transfer q light of the expected spectrum of a light target is smaller than the corresponding q heavy of the spectrum of a heavy target, so that the overall rate of a heavy target is enhanced compared to that of a light one by a factor (q heavy /q light ) n ≫ 1 if the scattering amplitudes depend on q n with n ≫1. This effect is observed in the left-hand plots of Figs. 3 and 4, where large-enough values of l suppress the sensitivity of F in favour of the heavier Xe and I targets, that pick up and become the more competitive ones. In particular the latter turn out to have a comparable sensitivity because Xe is neutron odd and is disfavoured by the coupling, although it is favoured by the large exposure and low residual background of the XENON1T and XENON100 experiments, while I is proton-odd, so is favoured by the coupling, although it is disfavoured by the much lower exposure and the higher residual background of the COSINE-100 experiment.
In light of the above discussion, heavy targets are the most suitable to search for largespin WIMPs, if their expected rate is driven by a large-multipole operator 1 Figure 6. Nuclear structure functions W 00 M and W 00 Φ ′′ for xenon as calculated in [28] as a function of the nuclear recoil energy E R for an isoscalar interaction (τ = τ ′ = 0). recoil energy behaviour of the spectra in Fig. 5 is at strong variance with the usual WIMP DD paradigm, where the search for new physics is focused on the lowest-energy part of the spectrum. Indeed, among present experiments only PICO-60, that is a threshold detector, includes the full range of recoil energies in its analysis (although, as already explained, F targets are only sensitive to E R < ∼ 100-200 keVnr because of the cut from the velocity distribution). For instance, the Region of Interest (ROI) analyzed by XENON1T [61] is limited to E R < ∼ 27 keV, that of SuperCDMS [63] stops at E R = 100 keV, while that from COSINE-100 [65] is limited to E ee ≤ 6 keVee, corresponding to E R ≤ 18 keVnr for sodium targets and E R ≤ 60 keVnr for iodine targets. Clearly, all such experimental WIMPs. We did not include CRESST in our analysis because the nuclear response functions W τ τ ′ T X for tungsten are not available in the literature.
searches are not optimized to look for some of the signals shown in Fig. 5. In addition, COSINE-100 achieved an exclusion plot at the level of the DAMA effect (using the same N aI target material) in spite of measuring a residual background a factor ≃ 3 higher than that of DAMA, thanks to an aggressive background subtraction. However such fit was performed using the energy part of the spectrum for E ee ≥ 6 keVee assuming it signal free, an assumption that is clearly not valid for high-s effective operators. Actually, the role of the high-energy part of the spectrum in momentum-dependent effective models was already partially understood for j χ =1/2 in [17,68]. In particular in [17] the data from XENON100 were reanalyzed extending the recoil energy interval up to 240 keVnr. For this reason to estimate the present bounds for a xenon target in Figs. 3 and 4 we included the result of Ref. [17] besides that from XENON1T [61] 2 .

Prospects of improvement of present constraints by extending the experimental energy range
In light of what pointed out above, with the exception of PICO-60 the present bounds on high-multipole effective operators are expected to be substantially improved by extending the experimental energy windows beyond the present ones. In order to estimate such improvement, for each model we have recalculated the bounds on M/g for a xenon, germanium and sodium iodide target by requiring that the corresponding expected differential rate in the full energy range where it is non-vanishing stays below the same residual background levels (in events/kg/day/keV) achieved by present experiments at lower energies. Notice that experimentally the large energy part of the spectrum is devoid of all the uncertainties that arise close to threshold, where the efficiency of applied cuts, the energy scale from light-yields or quenching factors or the energy resolution are sometimes challenging to determine. In particular such effects can be safely neglected in the calculation of the differential rate at high recoil energies. there is no compelling reason for dark matter particles to be elementary, or for their spin to be limited to 0 and 1/2.
In the present paper we have provided a first systematic and quantitative discussion of the phenomenology of the non-relativistic effective Hamiltonian introduced in Ref. [32] to describe the nuclear scattering process for a WIMP of arbitrary spin j χ . To this aim we obtained constraints from a representative sample of present direct detection experiments assuming the WIMP-nucleus scattering process to be driven by each one of the 44 effective couplings O X,s,l that arise for j χ ≤2.
We found that high values of the multipole parameter s, related to powers of the momentum transfer q appearing in the scattering amplitude, can push the expected differential spectra to recoil energies E R much larger than usually assumed, with the largest part of the signal concentrated at E R > ∼ 100 keV and a peculiar structure of peaks and minima arising when both the nuclear target and the WIMP are heavy. This phenomenology is at strong variance with the usual WIMP DD paradigm, where the search for new physics is focused on the lowest-energy part of the spectrum. In particular we have shown that the present bounds on the effective operators can be significantly improved by extending the recoil energy intervals to recoil energies up to ≃ 1 MeV.
A large multipolarity s leads also to a suppression of the expected rates. In particular we found quantitatively that for s ≤ 4 the effective scales probed by direct detection experiments can be suppressed by up to 5 orders of magnitude compared to the case s=0. It is possible to conceive DM candidates whose interaction with ordinary matter is driven by the highest multipole moments connected to the high-rank operators whose phenomenology is the subject of our analysis. An example is provided by molecules in the dark sector, where a particularly high symmetry cancels all lower multipoles except the highest one. In this sense this paper provides the first phenomenological study of the direct detection of quadrupolar, octupolar, and hexadecapolar dark matter.

A Experiments
In this Appendix we summarize the procedure that we adopted to obtain the upper bounds discussed in Section 4, categorizing them according to the target: Xe, Ge, F and N aI.

A.1 Xenon target: XENON1T & XENON100
For XENON1T we have assumed 7 WIMP candidate events in the range of 3PE ≤ S 1 ≤ 70PE, as shown in Fig. 3 of Ref. [61] for the primary scintillation signal S1 (directly in Photo Electrons, PE), with an exposure of 278.8 days and a fiducial volume of 1.3 ton of xenon (corresponding to a residual background b res,XEN ON 1T ≃ 7.7×10 −7 events/kg/day/keV). We have used the efficiency taken from Fig. 1 of [61] and employed a light collection efficiency g 1 =0.055; for the light yield L y we have extracted the best estimation curve for photon yields n ph /E from Fig. 7 in [70] with an electric field of 90 V/cm (with these assumptions the energy range analyzed by XENON1T corresponds to 2 keVnr < ∼ E R < ∼ 27 keVnr). The energy resolution was modeled combining a Poisson fluctuation for the observed primary signal S 1 and a Gaussian response of the photomultiplier with σ P M T = 0.5 [71].
Given the relevance of high recoil energies to constrain the effective models discussed in the present paper, in our analysis we have included the study of Ref. [62], where the data from run II of XENON100 (34 kg × 224.6 live days) were reanalyzed in the increased recoil energy interval (6.6-240) keV nr . To calculate the bounds we combined the number of observed events n k and the expected background rates b k for the 9 bins (k=1...9) listed in Table I of [62] with the corresponding expected rates r k (m χ , M/g) in the likelihood: L(m χ , M/g) = 2 k [r k (m χ , M/g) + b k − n k log(r k + b k )] (A.1) and found the upper bounds on M/g imposing the condition L(m χ , M/g) − L min ≤ 1.64 2 at 90% C.L., with L min the minimum of L. To calculate the expected rates r k we have directly convoluted the signal model detector response tables provided for each of the 9 analysis bins in numerical form in [62] with the differential rate dR/dE R calculated in each effective model (see Eq.(B5) of [62]). In Section 4.4 we estimate the improvements to the present bounds by extending the experimental energy windows beyond the present ones. To do so for a xenon target we require the corresponding expected differential rate as given in Eq. (4.3) to be below the same residual background b res,XEN ON 1T ≃ 7.7×10 −7 events/kg/day/keV presently achieved by XENON1T at low energy in the full energy range where it is non-vanishing. In the calculation of the differential rate at high recoil energies we assume efficiency equal to unity and neglect the effects from the energy resolution.

A.2 Germanium target: SuperCDMS
The SuperCDMS analysis [63] in 2017 observed 1 event between 4 and 100 keVnr with an exposure of 1690 kg days (corresponding to a residual background b res,SuperCDM S ≃ 6.2×10 −6 events/kg/day/keV). We have taken the efficiency from Fig.1 of [63] and the energy resolution σ = √ 0.293 2 + 0.056 2 E ee from [72]. To analyze the observed spectrum we apply the optimal interval method [73].
In Section 4.4 for a germanium target we follow the same procedure of Section A.1 using b res,SuperCDM S ≃ 6.2×10 −6 events/kg/day/keV.

A.3 Fluorine target: PICO-60
PICO-60 is a threshold experiment utilizing a bubble chamber. We analyzed the data obtained with a C 3 F 8 target [64] using two thresholds: an exposure of 1404 kg day at threshold E th =2.45 (with 3 observed candidate events and 1 event from the expected background, implying an upper bound of 6.42 events at 90%C.L. [74]) and an exposure of 1167 kg day keV at threshold E th =3.3 keV (with zero observed candidate events and negligible expected background, implying a 90% C.L. upper bound of 2.3 events). For the two runs we have assumed the nucleation probabilities in Fig. 3 of [64].

A.4 Sodium Iodide target: COSINE-100
The exclusion plot for COSINE-100 [65] relies on a Montecarlo to subtract the different backgrounds of each of the eight crystals used in the analysis. In Ref. [65] the amount of residual background after subtraction is not provided, so we have assumed a constant background b at low energy (2 keVee< E ee < 8 keVee), and estimated b by tuning it to reproduce the exclusion plot in Fig.4 of Ref. [65] for the isoscalar spin-independent elastic case. The result of our procedure yields b res,COSIN E−100 ≃0.13 events/kg/day/keVee, which implies a subtraction of about 95% of the background. We take the energy resolution σ/keV = 0.3171 E ee /keVee + 0.008189E ee /keVee averaged over the COSINE-100 crystals [75] and the efficiency for nuclear recoils from Fig.1 of Ref. [65]. Quenching factors for sodium and iodine are assumed to be equal to 0.3 and 0.09 respectively, the same values used by DAMA. In Section 4.4 for a sodium iodide target we follow the same procedure of Section A.1 using b res,COSIN E−100 = 0.13 events/kg/day/keV.