Spin-dependent dark matter-electron interactions

Detectors with low thresholds for electron recoil open a new window to direct searches of sub-GeV dark matter (DM) candidates. In the past decade, many strong limits on DM-electron interactions have been set, but most on the one which is spin-independent (SI) of both dark matter and electron spins. In this work, we study DM-atom scattering through a spin-dependent (SD) interaction at leading order (LO), using well-benchmarked, state-of-the-art atomic many-body calculations. Exclusion limits on the SD DM-electron cross section are derived with data taken from experiments with xenon and germanium detectors at leading sensitivities. In the DM mass range of 0.1 - 10 GeV, the best limits set by the XENON1T experiment: $\sigma_e^{\textrm{(SD)}}<10^{-41}-10^{-40}\,\textrm{cm}^2$ are comparable to the ones drawn on DM-neutron and DM-proton at slightly bigger DM masses. The detector's responses to the LO SD and SI interactions are analyzed. In nonrelativistic limit, a constant ratio between them leads to an indistinguishability of the SD and SI recoil energy spectra. Relativistic calculations however show the scaling starts to break down at a few hundreds of eV, where spin-orbit effects become sizable. We discuss the prospects of disentangling the SI and SD components in DM-electron interactions via spectral shape measurements, as well as having spin-sensitive experimental signatures without SI background.


INTRODUCTION
Direct searches of dark matter (DM) through its scattering with electrons have being a rapidly growing field in the past decade. With low-threshold capabilities of modern detectors in electron recoil (ER) and new ideas inspired by theoretical studies, the coverage of DM mass, m χ , has gradually been extended from sub-GeV towards increasingly lower reach (for a comprehensive overview, see, e.g., Ref. [1]). So far, most attention is given to the DM-electron interaction which is spin-independent (SI) of both DM and electrons. Various experiments have already set stringent exclusion limits on its cross section, σ (SI) e , in the mass range of MeV to GeV [2][3][4][5][6][7][8][9][10][11][12][13]. The current best limit above 30 MeV is set by XENON1T for their huge exposure mass and time [8,14]; in the range of 1−30 MeV, several experiments capitalizing the condensed phases of materials such as semiconductor silicon [9,11] and germanium [10,12] show potential improvements upon future scale-up. But, what about the spin-dependent (SD) interactions?
While there can be numerous models to motivate studies of SD DM-electron interactions from top down, a more straightforward, bottom-up approach is through nonrelativistic (NR) effective field theory [15][16][17][18][19][20][21][22]. As most DM is believed to be cold, its velocity v χ ∼ 10 −3 provides a natural parameter by which the DM interactions can be formulated as a series expansion. At leading order (LO), i.e., v 0 χ = 1, there are only two terms: 1 χ · 1 e (SI) and S χ · S e (SD), where 1 and S are unity and spin operators, respectively. It is when matching these NR effective operators to their underlying theories that the generality and complementarity of such consideration shows its full power. For example, at the relativistic level, the SI term can come from either a scalar DM with a scalar-scalar coupling, or a fermionic DM with a scalar-scalar or vector-vector coupling to electrons; the SD term can come from a fermionic DM with an pseudovector-pseudovector or tensor-tensor coupling. As a result, null observations of them constrain different parts of the broad DM parameter space. On the other hand, an observation of the SD term would not only be a discovery of DM and a new interaction, but also exclude all scalar DM scenarios -this adds a stronger motivation to search for it.
The LO SD DM-nucleon interactions have been studied, though not as extensively as its SI counterpart, with good recent progress [23][24][25][26][27][28][29][30]. The best limits on the contact interactions with a neutron (n) and a proton (p), in terms of total cross sections, are σ (SD) n < 6.2 × 10 −42 cm 2 with m χ = 30 GeV by XENON1T [28], and σ with m χ = 25 GeV by PICO [29], respectively. Compared with the best limits of order 10 −46 cm 2 on the SI interactions, the huge order-of-magnitude difference is easily understood from the (almost) coherent contributions of A ∼ 100 nucleons for the SI and a single unpaired nucleon for the SD case, respectively, in scattering amplitudes. On the other hand, because a typical atom is of size ∼Å, a cold DM particle can not induce substantial coherent scattering with atomic electrons unless it is as light as m χ MeV. Therefore, in direct searches of DM in the MeV-GeV mass range, one can anticipate constraints of similar orders on both LO SD and SI DM-electron interactions, for the latter has no coherent enhancement in scattering cross sections.
In this work, we study DM-atom scattering through the SD DM-electron interaction at leading order, using wellbenchmarked, state-of-the-art atomic many-body approaches. Exclusion limits are derived from various data of xenon and germanium detectors. The limit on DM-electron cross section σ (SD) e < 10 −41 −10 −40 cm 2 with m χ = 0.1−10 GeV, set by XENON1T data, is comparable to the best ones on σ (SD) n,p mentioned above at slightly bigger m χ . Special attention is on the difference in detector responses to the SI and SD interactions, with spin-orbit effects being found to be the deciding factor. We discuss in the end several strategies for sub-GeV DM searches that can disentangle the SD and SI DM-electron interactions.

FORMALISM
At leading order (LO), the effective SD interaction between DM (χ) and an electron (e) is where the coupling constants c 4 and d 4 (following the convention of Refs. [15,16]) are for the short-and long-range interactions, respectively; and q the magnitude of the three-momentum transfer q.
The unpolarized differential scattering cross section in the laboratory frame can be derived straightforwardly (see, e.g., Ref. [31]), and we focus only on the ionization processes that yield ER signals dσ dT Note that for later convenience we redefine the coupling constants (c 4 ,d 4 ) = s χ /4(c 4 , d 4 ) by absorbing (i) the spin factors χ = s χ (s χ + 1)/3, which results from an average of the initial, and a sum of the final DM spin states, with s χ = 1/2 or 1 applying to the case of a fermionic or vector DM particle, and (ii) the square of the electron spin s e = 1/2. The SD response function involves a statistical average of the initial states |I and a sum of all final states |F allowed by energy conservation (imposed by the delta function). To incorporate relativistic effects in atoms, the electron spin operator is expressed by the 4 × 4 Dirac spin matrices σ D e = 2 S e , with σ D = γ 0 γ = σ 0 0 σ where σ are the normal 2 × 2 Pauli matrices.
The transition operator contains all Z electrons (the summation of i), and spin operators of orthogonal j-directions contribute incoherently (the summation of j). Our evaluation of R (ion) SD proceeds similarly to our previous work on the SI case [14]. All the essential formulae, including those new to the SD case, are given in Appendix A . To emphasize once again the importance of relativistic and many-body physics, we also highlight therein the differences from NR, independent-particle approaches widely adopted in literature.
While the computation of R having an additional spin operator σ D and the associated summation of three directions. Therefore, before presenting our full numerical results, it is instructive to study the possible relationship between them. For the purpose, we define a scaling factor ξ = R (ion) First in Ref. [32], it was shown in DM scattering off free electrons that ξ (0) = 3. Later in Ref. [31], the same conclusion was reached for hydrogen-like atoms that ξ (H) = 3. Recently, in Ref. [33], the same conclusion was generalized to many-electron atoms with two conditions: (i) the electron basis wave functions are labeled by the principle (for bound states) or energy (for continuum states), orbital, and magnetic quantum numbers: n or k (with E = k 2 /2m e ), l, m l , and (ii) the final states are purely one-particle-one-hole excitations from the ground state. In such a NR, purely mean-field framework, it also yields that ξ (NR) = 3. Because a NR electron's orbital and spin wave functions are factorized, the scaling factor of 3 can be heuristic deducted from the multiplicity of the σ (for the SD case) over the unity (for the SI case) operator, when the transition involves only one independent electron.
However, as there is spin-orbit interaction (SOI) in atoms, this scaling relation can only be approximate at best. Its effects show up in all parts of the response function. First, the SOI causes the energy shift of two electron states of same l but different total angular momenta j = l ± 1/2, so the energy conserving delta functions are different. Second, the electron basis wave functions (which are best solved by relativistic mean field methods), their associated one-electron matrix elements, and the density matrix (which contains the weights of all atomic configurations that can be mixed by the residual interaction) all have explicit dependence on j. Therefore, one chief goal of this paper is to quantify to what extent that ξ ≈ 3 is a good approximation to apply in DM-atom scattering.
It should be pointed out that the analogous behavior for DM-nucleus scattering is rather different. Unlike the atomic SOI, which is electromagnetic with its strength characterized by the fine structure constant α EM 1/137, the nuclear SOI is a component of the nuclear strong interaction, and its strength is of the order of the strong coupling constant α S ∼ 1. Consequently, there is no such simple scaling in DM-nucleus scattering through the SD and SI DM-nucleon interactions (see, e.g., Ref. [16]).

RESULTS
To get a prediction for the event rate at a detector with N T atoms the differential cross section is weighted and averaged by the standard Maxwell-Boltzmann velocity distribution of DM [34], f ( v χ ), with conventional choices of DM parameters (same as in Ref. [14]): the local density ρ χ = 0.4 GeV/cm 3 , circular velocity v 0 = 220 km/s, escape velocity v esc = 544 km/s, averaged Earth relative velocity v E = 232 km/s so that  [24], XENON10 [2], XENON100 [3], XENON1T [8], and PandaX-II [13] data. The SD and SI limits of XENON1T at 400 eV threshold are depicted as solid and dotted red lines, respectively, to illustrate their expected recoil spectral shape are different.
v max = v esc + v E , and minimum v min = 2T /m χ . We follow the same efficient velocity average scheme as in Refs. [32,35] [34,36]. This simplification has been explicitly verified in Ref. [14] for the same kinematic region.
In Fig. 1, the results of d σv χ /dT are shown for the cases with m χ = 1 GeV, c 1 =c 4 = 1/GeV −2 , and d 1 = d 4 = 10 −9 . The fully many-body calculations through (multiconfiguration) relativistic random phase approximation, (MC)RRPA, are performed at a few points (for them being computer-time-consuming) labeled by circles (open squares) for xenon (germanium) atoms as benchmarks. The main results in solid curves are obtained from a more efficient approach: the frozen core approximation (FCA), which is in the spirit of an independent particle approximation but with a more realistic mean field to compute the ionized electron wave function (details in Ref. [14]). In general, FCA shows good agreement with (MC)RRPA benchmark points except near edge energies. Considering the 5% accuracy of (MC)RRPA achieved in photoabsorption cross sections with T ≥12 and 80 eV for xenon [37] and germanium [38,39], respectively, we conservatively assign an overall theory error of 20%. (The theory error can be reduced to 5%, if we perform (MC)RRPA calculations throughout the work without the limitation of computer resources.) Complete data of d σv χ /dT with different m χ are provided in the folder of ancillary files.
Data from CDMSlite [24], XENON10 [2], XENON100 [3], XENON1T [8], and PandaX-II [13] are analyzed against the theory predictions of dR/dT in the same way as in Ref. [14]). In Fig. 2, the exclusion limits on the coupling constantsc 4 andd 4 , and equivalently the SD DM-electron total cross sections σ n,p at slightly bigger m χ . This is an important complement to mainstream searches of weakly-interacting massive particles through DM-nucleus scattering: not only because of some overlapping in mass range (1-10 GeV), but also for classes of DM which can be both hadroand lepto-phillic. For smaller m χ down to 10 MeV, PandaX-II instead set better limits because of its lower ER threshold at 80 eV. It defines the exclusion boundary at low m χ for having lower background ∼ 2.64 × 10 −4 cpd/kg/keV compared to ∼ 0.2 of XENON10, despite the latter has an even lower threshold at 13.8 eV.
We emphasize here the derived exclusion limits depend critically on the theory predictions for the DM scattering event rates. Substantial differences exist between our relativistic many-body approach and several NR mean-field approaches Refs. [2-4, 33, 35] including the package QEdark that was used in the data analyses of Refs. [8,13]. Detailed comparisons, reasons that cause the differences, and justifications of our approach are presented in Ref. [14]. Based on similar arguments, one would generically anticipate weaker limits on the short-range interaction in the high-mass region, and tighter limits on the long-range interaction in the low-mass region, when above-mentioned approaches are applied. The main reasons are: (i) the relativistic effects substantially increase the scattering event rate at high T , and (ii) the many-body effects play an important role at all distances, including the Coulomb screening at long distances which substantially decrease the scattering rate at very low T .
Not included in our data analyses are several low-threshold semiconductor experiments, which already set limits on σ (SI) e below 10 MeV [5][6][7][9][10][11] and have the potential to compete with liquid noble-gas detectors above the 10 MeV region, particularly for the long-range interactions, in the future. The associated many-body physics of condensed matter is beyond the scope of this work, but has been taken in Refs. [40][41][42][43][44][45][46] for the SI interaction, and recently in Refs. [22,[45][46][47] that also include the SD interactions. Currently, there are discrepancies among predictions from different theoretical schemes.

DISENTANGLEMENT OF SD AND SI INTERACTIONS
Conventional DM detectors only measure recoil energy. Detector responses to SI and SD signatures must be different in order to disentangle them. Therefore, if the NR scaling ξ (NR) = 3 is valid, the SD and SI signals are indistinguishable in the measurable recoil energy spectrum shapes and the exclusion limits of σ In the bottom panel of Fig. 1, this statement is examined through theξ parameter, which is the ratio of dT . At T 200 eV, the scaling relationξ ≈ 3 works well for both xenon and germanium. The scaling deviation starts to grow as T increases, and the larger deviations in xenon than germanium demonstrates the effects of a stronger SOI in an atom of higher Z. We note that the SOI is only part of the relativistic corrections in atoms. Therefore, even though it has been shown that relativistic corrections become sizable for T 200 eV in DM-xenon scattering [14], the departure ofξ from the NR predictionξ (NR) = 3 is comparatively modest.
It follows that precision spectral shape measurements at high energy allow SD and SI to be distinguished in the case of positive discovery signals. As SD limits are driven by detector thresholds lower than 200 eV, the exclusion curves of Fig. 2 differ only slightly from the ones for SI given in Fig. 2 of Ref. [14]. If a higher cutoff is applied, as the example shown in Fig. 2 with 400-eV cutoff to XENON1T data, noticeable difference shows up.
There are intriguing possibilities of experimental signatures unique to SD interactions without SI background. For example, in a spin polarizable target and with known spin states of the ionized electrons, only part of the SD interaction λ=±1 (−1) λ S χ,λ S e,−λ (where λ denotes the spherical component) can generate signals as spin-flips. There are intense recent interests to seek solutions in condensed phases of matter which have natural spin, or magnetic, order. For instance, a ferrimagnetic material: yttrium iron garnet has been proposed as a candidate target [22,47] for light DM searches with m χ ≤ 10 MeV. The SD DM-electron interaction can perturb this spin system and generate magnons, which are the low-energy quanta of spin waves and experimentally measurable.
Motivated by the role of SOI, another approach is to realize experimental configurations where spins couple strongly to orbital angular momenta. For atoms, the fine structure constant is fixed, so that SOI enhancement can only be found in the inner-shell electrons, which have higher ionization energies. The novel Dirac materials do have the desired property that the effective fine structure constants can be larger than α EM . Take graphene as an example: a naive dimensional analysis yields α g = α EM /v F ≈ 2.2, adopting a generic Fermi velocity at Dirac cone v F ≈ 10 6 m/s. While recent calculations [48] and measurements [49,50] reported smaller values in between 0.1 and 1 (mostly due to the suppression by Coulomb screening), they are still significantly bigger than α EM . The smallness of v F compared to the speed of light is the most important factor that lowers the bar to reach the relativistic limit. Unfortunately, the intrinsic SOI in graphene is found to be small [51], so whether it could generate a sufficiently different SD response requires more study. More promising candidates are topological insulators [52,53] and transition-metal dichalcogenides [54], whose structures critically depend on SOIs. In addition, among a few Dirac materials being proposed to search for SI DM-electron interactions, the semimetal ZrTe 5 [55][56][57] and organic crystal (BEDT-TTF) · Br [58] show corrections from SOIs in band structure calculations.

SUMMARY
The SD DM-electron interactions is an important, but less-attended subject in direct DM searches. Its studies complement the very active research on the SI interactions, and together they can provide a more comprehensive understanding about the nature of DM and its interactions with matter. We derive in this paper, for the first time, limits on the SD DM-electron cross sections at leading order with state-of-the-arts atomic many-body calculations and current best experiment data. As a result of spin-orbit interaction, the spectral shapes of SD and SI recoil spectra differ at high energies, which provides a way of differentiating them experimentally. In addition, signals unique to SD interactions and systems of enhanced SOIs that can be found in novel phases of matter may offer more direct experimental probes for SD from sub-GeV DM searches in the future.

Appendix A: Atomic Many-body Response Functions
It is a well-acknowledged fact that a detector's response to scattering of nonrelativistic sub-GeV DM involves many-body physics at atomic, molecular, or condensed matter scale, depending on the configuration of a detector. In this Appendix, we summarize the important steps and formulae in a truly many-body calculation of the atomic response functions, and illustrate the differences to approaches based solely on independent particle pictures, either relativistically or non-relativistically. We shall follow the notations and terminology of Ref. [16] (a seminal work on the nuclear responses to DM-nucleus scattering) as much as possible. However, a necessary addition is using Dirac spinors for electrons, as relativistic effects are important and manifest in atomic physics. A further note is while a formulation in molecular or condensed matter systems proceeds similarly in spirit, there are substantial changes of basic elements and notations which we leave for future work.
For brevity, consider only the short-range DM-electron interactions with both leading-order SI (c 1 ) and SD (c 4 ) terms included. The unpolarized differential scattering cross section in the laboratory frame is given by dσ dT wherec 4 = s χ (s χ + 1)/12 c 4 that absorbs the factors due to DM spin, s χ , and electron spin s e = 1/2. The SI and SD detector responses are encoded in two functions where the delta function imposes energy conservation E F = E I + T .
A substantial simplification to the problem can usually be achieved by multipole expansion. First, because total angular momentum and its z-projection, are good quantum numbers, so we label the initial and final states as |I = |I, J I M J I and |F = |F, J F M J F , respectively, where I and F are collective labels for other quantum numbers. Second, each multipole operator has its definite spherical rank, J, component, M J , and parity, so the selection rules help to organize all allowed transitions in a systematic way. Third, by the Wigner-Eckart theorem, the summation over M J F , M J I , and M J can be carried out easily, and what remain to be calculated are the reduced matrix elements. Most important of all, for low energy processes, the expansion in spherical rank converges efficiently.
Following the notation of Ref. [16], the SI and SD response functions can now be decomposed as where j J (qr) is the spherical Bessel function, Y M J J (Ω r ) the spherical harmonics of solid angle Ω r , and Y M J JL (Ω r ) the vector spherical harmonics formed by recoupling of Y M L L (Ω r ) and the unit vectorr, whose spherical projection is proportional to Y λ 1 . A general feature of all these matrix elements to be computed is they involve many-body states |I and |F , but the operators are one-body, i.e, only one electron makes transition at one time. Such matrix elements of a generic one-body operatorÔ (1) = iô (i) can be compactly written by second quantization that where the one-body density matrix gives the statistical weight of annihilating an electron with a quantum state labeled by β and creating an electron with a quantum state labeled by α in this transition, and the total transition amplitude is the weighted sum of each single-particle transition with an amplitude given by α|ô|β (the reference to the ith electron is dropped because electrons are identical. Implementing the above scheme to a spherical basis with spherical operators is standard but too lengthy to reproduce here, but we should mention two basic ingredients. First, all single-particle states have good quantum labels of total angular momentum, so one can label them as |α = |aj a m ja and |β = |bj b m j b . Second, by recoupling the creation and annihilation operatorĉ † α andĉ β , a spherical tensor operator of rank J and component M J can be formed where the braket is the standard Clebsch-Gordan coefficient. Consequently, the reduced matrix element of a spherical operator can be recast as and the resulting density matrix is also reduced in the sense that it no longer depends on total angular momentum substates m ja and m j b , and the spherical multipole component M J . At this point we should note that as long as many-body wave functions can be obtained exactly, the choice of a complete set of single-particle basis states {|α , |β , . . .} is irrelevant, because the density matrix is also known exactly. However, this does not happen in practical cases, and the arts of many-body methods is mostly about choosing good sets of basis states and solving the density matrix reliably. The approach by Ref. [16], the nuclear shell model, is built upon a simple basis spanned by the harmonic oscillator eigenstates, which have nice properties that yield algebraic results of single-particle matrix elements. This simplistic mean field is compensated by diagonalizing the residual interaction in a large, but still truncated, model space with proper care in using the effective operators. Our approach, the (multiconfiguration) relativistic random phase approximation ((MC)RRPA), proceeds in a rather different direction. The mean field is first solved self-consistently by the Dirac-Hartree-Fock equation, so the singleparticle basis states are optimized, i.e, they already include a large part of many-body physics. Then the effects due to the residual interactions are accounted for by the random phase approximation. The resulting density matrices, though not as dense as the ones of shell model approaches, are certainly not sparse when residual correlation induces substantial configuration mixing.
which depend on the single-particle orbital wave functions. For the electrons in the core, as the mean field methods based on variational principles are designed to get a best approximation for the ground state energy, the resulting wave functions can be taken with certain confidence. However, it is more problematic for the ionized electron, and such a highly excited state goes beyond the reach of conventional mean field design. There are many methods applied to DM-atom scattering, including plane wave, hydrogen-like, and frozen core approximations. We refer readers to Ref. [14] for more discussions and comparison of these methods, and only point out here that even though our main results in this paper are obtained with a frozen core approximation (FCA), the latter approach is justified by a few benchmark calculations using (MC)RRPA.

Nonrelativistic Limit
By further taking the nonrelativistic (NR) limit to the above independent particle approximation, we now show explicitly how our formulation can be reduced to the conventional form widely used in the literature, and recover the NR scaling factor ξ (NR) = 3.
In numerical computations, the NR limit can be implemented by taking the speed of light to infinity. Theoretically, the most important features are (i) the Dirac spinor now has a vanishing small component and collapses to a Pauli spinor, i.e., all functions f εlj (r) can be taken to zero, and (ii) the radial wave functions now only depends on ε and l, since there is no longer spin-orbit interaction to break the j = l ± 1/2 degeneracy. As a result, the only change in our formulation is in the radial integral where u εlj (r)'s are the NR radial wave functions.
In the NR-IPA scheme, the SI response function