Effective Field Theory of Dark Matter Direct Detection With Collective Excitations

We develop a framework for computing light dark matter direct detection rates through single phonon and magnon excitations via general effective operators. Our work generalizes previous calculations focused on spin-independent interactions involving the total nucleon and electron numbers $N$ (the usual route to excite phonons) and spin-dependent interactions involving the total electron spin $S$ (the usual route to excite magnons), leading us to identify new responses involving the orbital angular momenta $L$, as well as spin-orbit couplings $L\otimes S$ in the target. All four types of responses can excite phonons, while couplings to electron's $S$ and $L$ can also excite magnons. We apply the effective field theory approach to a set of well-motivated relativistic benchmark models, including (pseudo-)scalar mediated interactions, and models where dark matter interacts via a multipole moment, such as a dark electric dipole, magnetic dipole or anapole moment. We find that couplings to point-like degrees of freedom $N$ and $S$ often dominate dark matter detection rates, implying that exotic materials with orbital $L$ order or large spin-orbit couplings $L\otimes S$ are not necessary to have strong reach to a broad class of DM models. We highlight that phonon based crystal experiments in active R&D (such as SPICE) will probe light dark matter models well beyond those having a simple spin-independent interaction, including e.g. models with dipole and anapole interactions. Lastly, we make publicly available a code, PhonoDark, which computes single phonon production rates in a wide variety of materials with the effective field theory framework.

The smallest DM mass that can be probed is limited by the band gap in these materials, typically For sensitivity to smaller energy deposits, and optimal reach to light DM and mediating particles, we look toward excitations at sub-eV energies. Such excitations exist and are derived from collective behaviors of atoms, ions or electrons in condensed matter systems. Phonons were proposed in Ref. [45] and further studied in Refs. [46][47][48][49][50] for direct detection in superfluid helium (where maxon and roton excitations also contribute), and were also discussed in the context of bosonic DM absorption in superconductors [51] and semiconductors [52], though ultimately, acoustic and optical phonons in (polar) crystals were advanced [53] and shown to have the best experimental prospects and sensitivity to light dark matter [33,34,[54][55][56]. Magnons -quanta of collective spin excitations -were also proposed in Ref. [57]. Both phonons and magnons in crystal targets have typical energies up to O(100 meV). To date, the work in the literature has focused on demonstrating the sensitivity of phonons and magnons to simple DM models. Only spin-independent (SI) interactions, via couplings to linear combinations of the proton, neutron and electron numbers, have been considered for phonon excitations, while a few benchmark models have been studied for magnon excitations.
The goal of this paper is to extend these results to general types of DM interactions. Effective field theory (EFT) is well-suited for this purpose: we can match a relativistic theory of DM onto a nonrelativistic (NR) EFT, and then compute the target response to the EFT operators. Within this framework, starting from a UV model consisting of relativistic operators coupling the DM to the proton, neutron, and/or electron, we can systematically compute direct detection rates via single phonon and magnon excitations in various target materials. The idea is along the lines of previous works on EFT calculations of nuclear recoils [58][59][60][61][62][63] (which extend earlier studies focused on standard SI and spin-dependent (SD) DM-nucleon interactions), and, more recently, of electron excitations in atoms [27] and crystals [37] (which extends earlier studies focused on SI DM-electron interactions), but technically there are important differences. Specifically, our EFT approach to DM-induced single phonon and magnon excitations consists of the following steps: 1. Matching of a relativistic theory of DM interactions onto the NR EFT (DM model-specific).
2. Matching of NR operators onto DM couplings to lattice degrees of freedom (universal).
3. Calculation of phonon or magnon excitation matrix elements (target-and excitation-specific).
We explain each of these steps in the three subsections of Sec. II. The first step -matching relativistic DM theories to the NR EFT -follows a similar procedure as previous works [27,37,[58][59][60][61][62][63], but involves a larger set of independent operators due to the absence of Galilean invariance in a medium. For nuclear recoils, one then derives the nuclear responses to the EFT operators.
Analogously, the key quantities in the present case are crystal responses which determines how DM couples to the collective excitations. (We emphasize, however, that despite the similar choice of terminology, collective excitations are associated with a different kinematic regime and degrees of freedom than nuclear recoils and therefore require a distinct EFT calculation.) Technically, for both phonon and magnon excitations in crystal targets, the second step listed above involves matching the NR EFT of DM-nucleon and DM-electron interactions onto an effective scattering potential that involves ionic degrees of freedom in the crystal lattice -in the long wavelength (low momentum transfer) limit relevant for light DM scattering, these (as we will highlight throughout) are quantities that characterize an ion as a whole, including the total particle numbers N ψ for the proton, neutron and electron (ψ = p, n, e), total spins S ψ , orbital angular momenta L ψ , as well as spin-orbit couplings L ψ ⊗ S ψ (a tensor with components L i ψ S k ψ summed over the constituent nucleons/electrons). Finally, in the third step, we quantize the scattering potential to obtain the phonon and magnon modes in a specific target material and compute the matrix elements for exciting them. All four types of crystal responses highlighted above can lead to phonon excitation in appropriately chosen targets, while S e and L e can also lead to magnon excitation.
Our new results significantly extend the searchable DM model space via phonon and magnon excitations, which we showcase in Sec. III with a variety of well-motivated benchmark models. We present full numerical calculations for several representative target materials, and apply simple analytic estimates to understand the results. We compare for which operators and interactions one expects phonon versus magnon excitations to dominate the rate, quantifying and generalizing the discussion in Ref. [57]. These calculations highlight the complementarity between phonon and magnon excitations, and between different targets, in probing the light DM theory space.
Our code for computing single phonon excitation rates, PhonoDark, is publicly available [64] and will be explained in detail in a forthcoming publication; it integrates the open-source phonon eigensystem solver phonopy [65], and takes general NR EFT operator coefficients, together with density functional theory (DFT) calculations of material properties, as input. Our magnon code, based on the Toth-Lake algorithm [66] for solving the magnon eigensystem, is also available upon request.

II. EFFECTIVE FIELD THEORY CALCULATION OF DARK MATTER INDUCED COLLECTIVE EXCITATIONS
Our goal is to present a framework for computing direct detection rates for general DM models, for the process where a DM particle scatters off a crystal target and induces a quasiparticle excitation in the crystal. This quantum mechanical process follows Fermi's golden rule which, when the incoming and outgoing DM particles are momentum eigenstates in free space, takes the form where v is the incoming DM's velocity, V is the total target volume, |i and |f are the initial and final states of the target system (defined with NR normalization: i|i = f |f = 1), and V is the Fourier transform of the scattering potential. The momentum transfer from the DM to the target, q, is integrated over, while the energy deposition onto the target is constrained to be See Ref. [33] for a review of the general formalism.
We now need to specify the type of transitions |i → |f in the target system to calculate the matrix element f | V(−q, v) |i . Here we focus on excitation of single phonon or magnon in a crystal target at zero temperature. We therefore take |i to be the ground state |0 , and |f to be the one-phonon or one-magnon states |ν, k , labeled by branch ν and momentum k within the first Brillouin zone (1BZ). For a crystal target, we write the scattering potential as a sum of contributions from individual ions: 1 where l = 1, . . . , N labels the primitive cells, j = 1, . . . , n labels the ions within each primitive cell, and x lj is the position of the ion labeled by l, j. Therefore, and we obtain The central quantity for the rate calculation is then the lattice potential V lj which the DM senses. This will depend on both the specific DM model and on the lattice degrees of freedom (e.g. the nucleon/electron number or total electronic spin of the ions) available to scatter from.
We will determine the lattice potential V lj in two steps previously mentioned in the introduction: first, in Sec. II A, we review the procedure of matching relativistic DM models onto NR effective operators; next, in Sec. II B, we further reduce these effective operators to DM couplings to the lattice degrees of freedom. In the simplest case of SI interactions, there is only one effective operator, O 1 = 1, and V lj is a linear combination of N p lj , N n lj and N e lj (proton, neutron and electron numbers of the ions, respectively) [33]. More generally, a DM model can generate a larger set of effective operators that involve the spins, momentum transfer, and velocities. The resulting lattice potential V lj depends on lattice degrees of freedom that include, in addition to the particle numbers N ψ lj (ψ = p, n, e), also their spins S ψ lj , orbital angular momenta L ψ lj , as well as spin-orbit couplings L ψ ⊗ S ψ lj (a tensor with components L i ψ S k ψ lj , see Eq. (24) below). The last step in computing the scattering rate is to quantize V lj in terms of phonon/magnon creation and annihilation operators; we carry out this exercise in Sec. II C. The framework in this section will provide the basis for concrete calculations of direct detection rates via single phonon and magnon excitations, and will be applied to a set of benchmark models in Sec. III.

A. From Dark Matter Models to Nonrelativistic Effective Operators
In this subsection, we take a top-down approach in deriving the EFT, focusing on how the effective operators arise from NR matching of well-motivated relativistic models. While one can also 1 For simplicity, we will refer to either atoms or ions on lattice sites as ions. take a bottom-up approach as in e.g. Ref. [58], and construct the EFT by enumerating operators consistent with rotation and translation invariance, we find it useful to have a set of benchmark UV models to develop intuition on how realistic theories of DM, which often predict correlations between EFT operators [61,67], can be probed by experiment.
Let us start from a relativistic model of a DM particle χ interacting with the proton (p), neutron (n) and electron (e); 2 we denote these Standard Model (SM) particles collectively by ψ in the following. To compute the NR EFT, we take the NR limit of the relativistic theory and map it on to the appropriate NR degrees of freedom. The EFT consists of the NR fields χ ± , ψ ± , generally defined by (using the SM fermion ψ for example): Here the sum is over energy eigenstates, ε I = E I − m ψ are the energy eigenvalues minus the rest mass, Ψ I (x) are the wavefunctions (which are two-component for spin-1 2 fermions) andb I are the annihilation operators. In the familiar case of a fermion in free space, the energy eigenstates are labeled by momentum k and spin s = ±, with eigenvalues ε k, where ξ + = ( 1 0 ), ξ − = ( 0 1 ). For a spin-1 2 fermion, the relation between the relativistic field ψ and NR field ψ + is (see at leading order in m −1 ψ , where k, ε are operators acting on ψ + . For a fermion in free space, we have k = −i∇, ε = i∂ t , which become simply numbers in momentum space. In the presence of an external potential (Φ, A) (e.g. electromagnetic fields from the ions), k = −i∇ − A is the kinematical momentum, while ε = i∂ t − Φ. Eq. (8) applies for the SM fermions ψ = p, n, e. If the DM χ is a spin-1 2 fermion, it also applies for the DM, with ψ replaced by χ. For a spin-0 DM, on the other hand, χ = e −imχt χ + , with χ + given by Eq. (7) without the ξ s factor. 2 The DM-proton and DM-neutron couplings follow from the DM-quark and DM-gluon couplings in the fundamental Lagrangian by standard methods, see e.g. Ref. [59]. 3 In this and the next subsection, we shall use k to denote a SM fermion's momentum while deriving the lattice potential, which should not be confused with the phonon momentum in Eq. (5). Afterward, starting from Sec. II C, we will no longer need to deal with fermion momenta, and the notation k will be recycled for phonon momentum.

Lagrangian Term
Coupling Type (Effective) Current → NR Limit Types of couplings between a spin-1 2 fermion ψ and a scalar (vector) mediator φ (V µ ). The (effective) currents are defined by L ⊃ g X φJ X (X = S, P ) or g X V µ J µ X (X = V, A, edm, mdm, ana, V 2), upon integration by parts in the last four cases. The expressions following the arrows are the leading operators in the NR reduction of the currents (assuming scattering kinematics), which appear between the nonrelativistic fields ψ − and ψ + -see e.g. Eq. (9). These will be used to derive the NR operators generated by specific DM models involving tree-level exchange of a scalar or vector mediator in Table II. To demonstrate the procedure of matching a relativistic model onto the NR EFT, we focus on tree level DM scattering mediated by a spin-0 or abelian spin-1 particle, denoted by φ and V µ respectively. While it should be kept in mind that the EFT is capable of describing a broader class of models, including e.g. loop-mediated scattering, we find it useful to organize our thinking by categorizing mediator couplings to fermion bilinears. In Table I, we list the commonly considered types of couplings at the level of the relativistic Lagrangian, and their NR limits. We explain the table in detail in the following two paragraphs.
For a spin-0 mediator φ, we consider its couplings to the scalar and pseudoscalar currents J S , J P . For a spin-1 mediator V µ , we consider both minimal coupling to the vector and axial-vector currents J µ V , J µ A , and non-minimal couplings to the field strength V µν . 4 The latter include a series of higher dimensional operators. At dimension five, we have the electric dipole moment (edm) and magnetic dipole moment (mdm) couplings. Upon integration by parts, they can be cast in the same form, V µ J µ , as the minimal coupling case, with effective currents J µ edm , J µ mdm listed in the last column of Table I. Next, at dimension six, we consider ∂ ν V µν coupling to the axial-vector and vector currents. The former represents a new type of coupling known as the anapole [67,69,70], and the corresponding effective current is denoted by J µ ana . On the other hand, ∂ ν V µν coupling to the vector current gives an O(q 2 ) contribution to the same form factors that J µ V induces (i.e. the familiar charge and magnetic dipole in quantum electrodynamics), so we denote the effective current by J µ V 2 . It is useful to note that all the (effective) currents that couple to a spin-1 mediator, except J µ A , are conserved: q µ J µ X = 0 (X = V, edm, mdm, ana, V 2). In the NR limit, we can substitute Eq. (8) for the relativistic fermion field ψ into the expressions for the (effective) currents in Table I, and expand in powers of k m ψ and ε m ψ . For example, for J µ V = (J 0 V , J V ), we find, at leading order, where S ψ = σ 2 is the spin operator, and with k defined as acting on the ψ − field on the left, k = i ← − ∇ −A, giving the kinematical momentum of the final state ψ. We can carry out the same exercise for the other (effective) currents. The results, up to the first nonvanishing order, are listed after the arrows in the last column of Table I, with ψ − on the left and ψ + on the right implicit. We see that all currents reduce to operators involving S ψ , K and iq; in the case of the electric dipole coupling, ω ≡ ε − ε also appears, with ε defined as acting on ψ − on the left.
With Table I, it is straightforward to derive the NR effective operators generated by tree-level exchange of a spin-0 or spin-1 mediator between a DM current and a SM current. Concretely, let us consider a set of benchmark models of spin-1 2 DM [61], listed in Table II. In each model, the DM χ and a SM fermion ψ each couple to the mediator via a linear combination of the currents in the last column of Table I, whose NR limits can be directly read off. Integrating out the mediator, 4 Other operators, such as those with derivatives acting on ψ and those involving the dual field strength Vµν , are not independent -see e.g. Ref. [68].
we then arrive at a NR EFT for DM scattering of the form For convenience, we reserve k and k for the momentum operators acting on ψ ± , and write the same operators as p and p when they act on χ ± . We normalize the operators by powers of m ψ so that O  Table II (to be discussed in detail shortly). These coefficients contain all the information for constructing the lattice potential V lj for a given DM model, and will be exploited below for computing the DM detection rate.
For kinematic conventions, we take to denote the momentum transfer from the DM to the target, which agrees with Refs. [59,63] but has an opposite sign compared to the definitions in Refs. [58,[60][61][62]. There are two other independent combinations of momenta: where P = p + p, K = k + k. Note that v χ should not be confused with the incoming DM's velocity, which we denote by v = p mχ ; the two are related by v χ = v − q 2mχ . The list of NR operators O (ψ) i up to linear order in v χ , v ψ is presented in Table III (grouped into four categories to be explained below). These encompass all the operators generated at leading order in the benchmark models in Table II. Our operator basis here is an extension of the familiar one from previous works on the EFT for direct detection via nuclear recoils [58][59][60][61][62][63]. In the latter case, due to Galilean invariance, NR effective operators involve only the linear combination v ⊥ ≡ v χ − v ψ . In contrast, for collective excitations considered in this work, in-medium effects, which break Galilean invariance, can be important, so v χ and v ψ must be treated separately. We adopt the operator numbering convention of Ref. [63], and split each v ⊥ -dependent operator into two; one term dependent on v χ and the other v ψ . For example, O  5 Note that the standard SD interaction cannot be realized with a light mediator. In that case the leading interaction

Model UV Lagrangian NR EFT Responses
Standard SI φ gχJS,χ + g ψ J S,ψ or c a Heavy mediator only. in the second to last column. g eff ψ are the screened couplings defined in Eq. (14), and µ ψ = 1 + δ µ ψ is half the Landé g-factor of ψ ( µ p 2.8, µ n −1.9, µ e 1). The last column lists the lattice degrees of freedom which enter the scattering potential, Eq. (27). All models can excite phonons, and models with S or L response generated by DM-electron coupling can also excite magnons. motivated class of (hidden sector) models contain DM particles coupling to a vector mediator via a multipole moment, which in turn kinetically mixes with the photon (see e.g. Refs. [61,67,[69][70][71][72][73][74]).
We consider the electric dipole, magnetic dipole and anapole DM models, which generate O 11 , O 1,4,5a,5b,6 and O 8a,8b,9 , respectively. Finally, Table II includes a model where a vector mediator couples to the SM fermion's magnetic dipole moment J µ mdm , and as a result generates O 3b . Among other things, this leads to a coupling to the SM fermion's spin-orbit coupling, which can be the leading interaction if one simultaneously introduces a coupling to the "O(q 2 ) vector current" J µ V 2 (see Table I), with a coefficient (relative to J µ mdm ) tuned to κ = −1 to cancel the standard SI interaction O 1 .
We also note that, in the case of a vector mediator coupling to the electron's vector current J µ V,e , in-medium screening effects modify the effective couplings to the proton and electron [33,39,42,43,75]. For NR scattering, screening is negligible for transverse photon exchange, but can be significant for longitudinal photon exchange, which generates the first category of operators (O 1,5a,8a,11 ) in Table III. As shown in Refs. [33,75], this amounts to replacing where ε is the dielectric tensor, and g p,e are the tree-level (unscreened) couplings. The same is true for a scalar mediator coupling to the electron's scalar current J S,e [76]. For single phonon and magnon excitations below the electronic band gap that we focus on in this work, one can use the high-frequency dielectric ε ∞ , which captures the screening due to fast electron responses [33,54,77].
We will study the reach phonon and magnon detectors have to these benchmark models in Sec. III, after developing the formalism of rate calculations within the EFT in the rest of this section.

B. Matching Effective Operators Onto Lattice Degrees of Freedom
We now match the effective operators O (ψ) i onto lattice degrees of freedom (highlighted for clarity) that appear in the DM-ion scattering potentials V lj . In Table III, we have organized the operators into four categories, according to whether O ∝ S ψ ("coupling to spin"), and whether the operator involves v ψ . Since our focus is light DM that evades conventional searches via nuclear recoils and electronic excitations, we will work in the long wavelength limit, where the momentum transfer is small compared to the inverse ionic radius (corresponding to m χ 10 MeV), so at leading order, the only relevant degrees of freedom are those that characterize the ion as a whole. Intuitively, we expect couplings to charge and spin of a constituent particle ψ = p, n, e to match onto couplings to the total number N ψ and spin S ψ of that particle, respectively. These are point-like degrees of freedom that do not involve the internal motions of the ion constituents; they are the only degrees of freedom to which DM couples if the operator is velocity-independent. On the other hand, v ψ -dependent operators are expected to couple DM to the motion of ψ particles inside an ion, manifest as the total orbital angular momenta L ψ and spin-orbit couplings L ψ ⊗ S ψ , which are "composite" degrees of freedom.
In the rest of this subsection, we will see concretely how these intuitive expectations are borne out. The final result of this calculation is the lattice potential in terms of the NR EFT operator coefficients c (ψ) i , given below in Eq. (27). Since the calculation proceeds in much the same way for all operators in the same category, to avoid tedious repetition we pick one operator from each category to explain the procedure: 3b , with ψ taken to be one of p, n, e. To obtain the DM-ion scattering potentials V lj , we need to compute the matrix elements of these operators between the incoming and outgoing states of the DM-ion system. Since the initial and final DM states are plane waves, the DM part of the matrix element simply yields a phase factor, so where α runs over all the ψ fermions associated with the ion labeled by l, j, and · represents the ionic expectation value (assuming the ionic state is unchanged for the low energy depositions of interest). Computing these expectation values in full generality is a tedious task that involves numerical integration over nuclear and electronic wavefunctions. However, the calculation is dramatically simplified in the long wavelength limit of interest here, where we can expand e iq·xα = 1 + iq · x α + . . . and keep just the leading nonvanishing terms. In the following two paragraphs, we discuss in turn the v ψ -independent operators O (ψ) (first line of Eq. (15)) and 4 . For these, it is sufficient to set e iq·xα → 1: So we obtain, respectively, the expectation values of the number and total spin of ψ particles for ion l, j, as one would expect for the lowest order "coupling to charge" (O (ψ) 1 ) and "coupling to spin" (O (ψ) 4 ) operators. Note that S ψ lj should not be confused with the total nuclear or ionic spin, which may also contain orbital angular momentum components. We will see in the next subsection that the total ionic spin (from electrons) is relevant for magnon excitations, and we will need to work out its decomposition into spin and orbital components (see Eq. (32) below); the total nuclear spin, on the other hand, does not enter the calculation of phonon or magnon excitations.
probability current, and its treatment is analogous to the nuclear recoil calculation [58]. Assuming the ionic states are energy eigenstates implies that the probability density is constant in time, and therefore by the continuity equation, ∂ i v i ψ,α lj = 0. This means that v i ψ,α can be written as a total derivative, v i ψ,α = ∂ k x i α v k ψ,α , and therefore has vanishing expectation value. The leading contribution then comes from expanding the e iq·xα to the next order in q: To go further, we note that x i α v k ψ,α lj is anti-symmetric in i ↔ k since the symmetric part can be written as a total derivative, α and therefore has vanishing expectation value. Expanding the anti-symmetric part gives which after integration by parts can be simplified to where L α is the angular momentum operator. We therefore have where are Cartesian components of the spin-orbit coupling tensor. So we finally obtain As alluded to previously, the v ψ -dependent operators O 3b induce DM couplings to the ψ particles' total orbital angular momentum and spin-orbit coupling.
We can carry out the same calculation for the other operators in Table III. The result is where v χ = v − q 2mχ (with the incoming DM's velocity v and momentum transfer q to be integrated over when calculating detection rates), and summation over repeated Cartesian indices is implicit.
To summarize, in the long wavelength limit, the DM-ion scattering potential V lj involves a set of quantities that characterize properties of the ion: the total numbers N ψ , spins S ψ and orbital angular momenta L ψ of the constituent particles ψ = p, n, e, as well as the spin-orbit coupling tensors L ψ ⊗ S ψ . We will refer to these as different types of crystal responses, as DM couplings to these quantities drive collective excitations in the crystal; they play a similar role to the nuclear responses in nuclear recoil calculations (which similarly reduce to the total nucleon numbers, spins, etc. in the long wavelength limit [58,[60][61][62]). We emphasize, however, that in contrast to standard nuclear recoil where nuclei are treated as free -a valid approximation at energy depositions 500 meV [33] -collective excitations arise in a lower energy regime where inter-ionic interactions become important; the EFT therefore involves different degrees of freedom and the calculation proceeds differently. 6 We will sometimes abbreviate the crystal responses v ψ -independent operators (the first two categories in Table III) generate point-like responses, while v ψ -dependent operators (the last two categories in Table III)  although from the top-down point of view, it seems difficult to construct well-motivated simple models that dominantly generate a composite response (L or L ⊗ S) without being accompanied by a point-like response (N or S) of at least comparable size, similar to the case of nuclear recoil as highlighted in Ref. [61]. We will elaborate on this in Sec. III.

C. Quantization of Lattice Potential for Phonons and Magnons
Now that we have obtained V lj in terms of the lattice degrees of freedom, Eq. (27), it remains to compute the matrix elements by quantizing the lattice potential in terms of phonon or magnon modes. The simplest cases, where phonon excitations in a crystal proceed through N ψ (via the SI operator O 1 = 1) and magnon excitations proceed through S e were considered previously in Refs. [33,34,53,54] and Ref. [57], respectively. Here we extend those calculations to include all four crystal responses ( N ψ , S ψ , L ψ , L ψ ⊗ S ψ ) identified in the previous subsection, which can be generated by the full set of effective operators.
Phonons arise from the ions' displacements with respect to their equilibrium positions x 0 lj : Recall that N (without subscript, not to be confused with N ψ ) is the total number of primitive cells in the crystal lattice, to be sent to infinity at the end of the calculation. The phonon creation and annihilation operators satisfy the canonical commutation relations, [â ν,k ,â † ν ,k ] = δ νν δ k,k with all others vanishing. The eigenenergies ω ν,k and eigenvectors ν,k,j (normalized such that j | ν,k,j | 2 = 1) are solved for by diagonalizing the quadratic crystal potential. The quadratic crystal potential, and equilibrium positions, are computed with DFT [78] (see Refs. [34,54] for details) and the diagonalization is performed with phonopy [65]. At leading order, dependence of the matrix element in Eq. (28) on u lj comes only from the phase factor e iq·x lj ; we assume the DM-ion scattering potentials V lj (−q, v) are not significantly affected by ionic displacements and can thus be pulled out of the matrix element. 7 Then, evaluating the matrix element of 7 If V lj receives contributions from DM-electron couplings, the scattering potential can depend on u lj directly, as ionic displacements distort the electron wavefunctions. This correction can be taken into account via the Born effective charges in the case of SI interactions in the long wavelength limit, as discussed in Ref. [33]. the phase factor, ν, k| e iq·x lj |0 , follows the standard procedure of expanding x lj as in Eq. (29) and applying the Baker-Campbell-Hausdorff formula to normal-order the phonon creation and annihilation operators [33]. As a result, where is the Debye-Waller factor. Crucially, the 1 √ N factor (which originates from Eq. (29) and is to be squared when computing the rate), together with the prefactor 1 V in the rate formula Eq. (5), indicates that the rate Γ would scale as 1 N 2 → 0 unless the l sum in Eq. (30) scales with N . This in turn requires the N terms in the l sum to add up coherently, which is possible only when i) the phonon momentum k matches the momentum transfer q up to reciprocal lattice vectors, which is the statement of lattice momentum conservation, and ii) l V lj ∼ N , i.e. the DM couples coherently across the crystal lattice. The second requirement is trivially satisfied for DM couplings to the scalar quantities N ψ , tr( L ψ ⊗ S ψ ). For couplings to the vector and tensor quantities S ψ , L ψ , L ψ ⊗ S ψ (modulo the trace part), on the other hand, coherence is possible only when they are ordered (or polarized), so that they point in the same directions in all primitive cells; in the case of S ψ , this can be achieved by spontaneous magnetic ordering for ψ = e, or by applying an external magnetic field for ψ = p, n.
Up to possible small corrections due to the presence of different isotopes, we can set V lj = V j , which is independent of l. We then obtain the single phonon excitation rate: where Ω is the volume of the primitive cell, x 0 j is the equilibrium position of the jth ion with respect to the cell center, and it is implicit that q = k + G where G is a reciprocal lattice vector.
To map q to a vector k within the 1BZ, we first write q = 3 i=1 a i b i , with b i the basis vectors of the reciprocal lattice, then construct a set of eight candidate G vectors whose components in reduced coordinates take the floor and ceiling integer values of a i , and finally choose the correct G vector to be the one that minimizes |q − G|.
The DM-ion scattering potential V j that enters Eq. (31) is simply given by Eq. (27) above, with the l subscripts dropped, assuming S ψ , L ψ , L ψ ⊗ S ψ are ordered, as explained above; in the absence of ordering, the corresponding terms should be dropped (with L ψ ⊗ S ψ set to its scalar component 1 3 tr( L ψ ⊗ S ψ ) 1 = 1 3 L ψ · S ψ 1). In the special case of SI interactions, one has only c (ψ) 1 , so V j = ψ c (ψ) 1 N ψ j , reproducing the results in Ref. [33], whereas in the full EFT, all four crystal responses can contribute to phonon excitations.
Next we move on to magnons. They are collective spin excitations in a magnetically ordered phase, and can thus respond to DM scattering only if the potentials V lj depend on the magnetic ions' effective spins S lj . Generally, S lj can come from the electrons' spin and orbital angular momenta, S e lj and L e lj , respectively. When projected onto the Hilbert space spanned by S lj , they become where λ S,j , λ L,j are numbers (which we will say more about shortly). Therefore, from Eq. (27) we obtain the matrix element for exciting a magnon mode |ν, k : where iq m e q × S χ .
As in Eq. (27), we have defined q ≡ |q|,q ≡ q/q, and v χ = v − q 2mχ . Now we need to compute ν, k|S lj |0 . The calculation follows Ref. [57], which we encourage the reader to consult for more details. The magnetic order is captured by a set of rotation matrices R j that take each S lj to a local coordinate system where it points in the +z direction: We restrict ourselves to commensurate orders, in which case the rotation matrices R j do not depend on the primitive cell label l. We then apply the Holstein-Primakoff transformation and expand S lj around the ground state in terms of bosonic creation and annihilation operators: Magnon eigenstates are obtained by diagonalizing the spin Hamiltonian, which is specific to the target material; in the simplest cases, the target can be modeled by Heisenberg exchange interactions S lj · S l j between neighboring sites, while more complicated model descriptions are needed in other cases. For a general spin Hamiltonian, the diagonalization can be achieved by a Bogoliubov transformation in momentum space: where U, V are n × n matrices (with n the number of magnetic ions per cell), andb † j,k ,b j,k are the creation and annihilation operators for the magnon eigenstates satisfying canonical commutation relations, [b ν,k ,b † ν ,k ] = δ νν δ k,k with all others vanishing. An efficient algorithm for the diagonalization can be found in Ref. [66] (see also Refs. [57,79]). Now computing the magnon excitation matrix element ν, k|S lj |0 , and hence the DM scattering rate, is reduced to standard algebra. We obtain [57,79] where r j = (R xx j , R yx j , R zx j ) + i (R xy j , R yy j , R zy j ). As in the phonon case, it is implicit that k matches q up to a reciprocal lattice vector, q = k + G, due to lattice momentum conservation.
A comment is in order about the target choice. In the case where the total S lj involve only spin degrees of freedom (as is the case for yttrium iron garnet (YIG) discussed in Ref. [57]), λ S,j = 1, λ L,j = 0, and only the first two lines of Eq. (34) are relevant. Targets for which λ L,j = 0 are more exotic. One class of materials with λ L,j = 0 is spin-orbit-entangled Mott insulators [80][81][82], where the combined effect of crystal fields and spin-orbit coupling results in effective spins S j = 1 2 , and we can show that λ S,j = − 1 3 , λ L,j = − 4 3 (see Appendix B for details, and Refs. [81][82][83][84] for related discussion), so the magnetic ions' effective spins are in fact dominated by their orbital components.
Perovskite irridates such as Sr 2 IrO 4 [80,83] and Kitaev materials Na 2 IrO 3 , α-RuCl 3 [82,[84][85][86] are among the materials with this feature that have been actively studied recently by the condensed matter physics community. While perhaps futuristic as DM detectors, such materials have the novel feature of being sensitive to DM couplings with electrons' orbital angular momenta.
As a final remark, we note from the derivation above that when the same crystal response, S e or L e , excites both phonons and magnons, the phonon excitation rate is parametrically suppressed by q 2 m ion ω ∼ 10 −2 q  Table III with ψ = e, which generates S e response, single magnon excitation is expected to achieve better sensitivity than single phonon excitation for the same exposure and detector efficiency. On the other hand, since phonons can be excited also by other crystal responses, they have a broader coverage of the DM theory space. We will investigate the interplay between single phonon and magnon excitations in the context of our benchmark models in the next section.

III. APPLICATION TO BENCHMARK MODELS
We now apply the general results of the previous section to the set of benchmark models in In order to present the results in a concise way, let us introduce the following definitions. For single phonon excitation, we define (cf. Eq. (31)) where X represents one of the crystal responses, X = N , S, L, L ⊗ S; note that F (ψ) X,ν are vector (tensor) quantities when X = S, L (X = L ⊗ S), and will be written as F where X = S, L. These are formally analogous to polarization vectors of a vector field. In both Eqs. (39) and (40), k is the phonon momentum inside the 1BZ that satisfies q = k + G for some reciprocal lattice vector G; as emphasized below Eq. (31), k is uniquely determined by mapping q into the 1BZ through reciprocal lattice vectors. We further define a set of quantities Σ ν (q), for both single phonon and single magnon excitations, by (cf. Eq. (5)) We will refer to Σ ν (q), which have mass dimension −4, as "differential rates." Practically, Σ ν (q) are obtained simply by taking V lj in Eq. (27), substituting X ψ lj by F (ψ) X,ν (for ψ = p, n, e and X = N , S, L, L ⊗ S) or E X,ν (for ψ = e only, and X = S, L), squaring it and averaging over the DM's spin (which amounts to replacing S i χ S k χ → 1 4 δ ik ). As we will see, written in terms of the dimensionless quantities F (ψ) X,ν and E X,ν defined above, Σ ν (q) can be expressed in a concise form for each benchmark model. This will be convenient when we compare the rates between different models, and between phonon and magnon excitations.
Our final results will be presented in terms of the rate per unit target mass, where ρ T is the target's mass density that we take from Ref. [87], ρ χ = 0.4 GeV/cm 3 is the local DM mass density, and f χ (v) is the DM's velocity distribution, taken to be a boosted and truncated Maxwell-Boltzmann distribution -see Appendix C for technical details of evaluating the velocity integrals. For the projected reach, we assume 3 events per kilogram-year exposure, corresponding to 95% C.L. exclusion in a background-free counting experiment, and assume a detector energy threshold of 1 meV. While we will present full numerical results, the main features can usually be understood by simple parametric estimates. Generally, noting that the velocity integral over the energy conserving delta function δ(ω ν,k − ω q ) yields a function that scales as q −1 (see Appendix C), we have from Eqs. (41) and (42), parametrically, where m cell = ρ T Ω is the mass of the primitive cell, and as before, q = |q|. Then, from the formulas for Σ ν (q) presented below for each model in terms of F (ψ) X,ν and E X,ν defined in Eqs. (39) and (40), we can estimate the rate R by In the case of single phonon excitations, we should further note that ω, which appears in F On the target side, we will consider the following representative set of materials: • GaAs [phonons, subject of R&D]. As the first-studied target for DM detection via phonons, GaAs is already in R&D as a target for both electron excitations and phonon excitations [88].
Phonons in GaAs form 3 acoustic and 3 optical branches, and have energies up to ∼35 meV.
• SiO 2 (quartz) [phonons, optimal sensitivity]. Based on our previous theoretical study comparing the phonon reach of a variety of target materials [34], we have advocated quartz as having good sensitivity to DM couplings to both acoustic and optical phonons. Also, quartz has complementary features compared to GaAs: while GaAs has a simple crystal structure and relatively low phonon energies, quartz has a large number of phonon branches (3 acoustic, 24 optical), with energies up to ∼ 150 meV.
• Y 3 Fe 5 O 12 (YIG) [mostly magnons, also phonons for comparison]. YIG is a well studied material with ferrimagnetic order, and is already used in an axion DM detection experiment QUAX [89][90][91][92][93]. The effective spin Hamiltonian is a Heisenberg model, with S j = 5 2 for the magnetic Fe 3+ ions coming entirely from electron spins S e (i.e. λ S,j = 1, λ L,j = 0 in Eq. (40)). We take the antiferromagnetic exchange coupling parameters from Ref. [94], together with the crystal parameters from Ref. [87], to compute the magnon spectrum and rotation matrices. YIG has 20 magnon branches, one of which is gapless and has a quadratic dispersion at small k. The gapped magnons have energies up to ∼ 90 meV. We will mostly consider YIG as a candidate material for DM detection via magnon excitations, but will also consider phonon excitations in YIG via DM couplings to the ordered electron spins in Sec. III A for comparison; in this case the scattering potential is determined by S e lj of the Fe 3+ ions, which have magnitude 5 2 and directions set by the ferrimagnetic order. YIG has 80 ions in total in the primitive cell and therefore 240 phonon branches (3 acoustic, 237 optical), with energies up to ∼ 120 meV.
• α-RuCl 3 [small-gap magnons with orbital component]. As discussed below Eq. (34), α-RuCl 3 is one of the materials where the effective ionic spins involve orbital degrees of freedom, and is therefore sensitive to DM couplings to the electrons' orbital angular momenta.
The magnetic ions Ru 3+ have S lj = 1 2 , coming from both S e and L e with λ S,j = − 1 3 , λ L,j = − 4 3 , as discussed in Appendix B. The effective spin Hamiltonian features Kitaevtype bond-directional exchange couplings. We use the Hamiltonian parameters derived from neutron scattering data in Ref. [95], which also includes an antiferromagnetic Heisenberg exchange; see Ref. [86] for a summary of some alternative model parameterizations derived from a variety of experimental and numerical techniques. The resulting magnetic order is zig-zag antiferromagnetic. Magnons in α-RuCl 3 , of which there are 4 branches, are at very low energy, below 7 meV, and can thus probe lighter DM than YIG. Also, since all magnon branches are gapped at zero momentum, the sensitivity is not significantly affected by the finite detector threshold. This is in contrast with YIG, where the assumed 1 meV energy threshold limits the momentum transfer to be greater than ∼ 80 eV in order to excite magnons on the gapless branch. Therefore, even though the experimental prospects of α-RuCl 3 itself are unclear, it can be regarded as a useful benchmark which highlights the generic advantage of small-gap targets.
Our main results are Figs. 1-4. We give a brief summary here and discuss them in more detail in the following subsections. A major issue of interest is the comparison of sensitivity to various types of DM interactions, via single phonon and magnon excitations induced by various crystal responses. First, we consider the standard SD interaction in Fig. 1, where we see that magnons outperform phonons, typically, by more than an order of magnitude in terms of the coupling reach.
Next, in Fig. 2, we compare the phonon and magnon rates for the four combinations of scalar mediator couplings; the phonon production rate is larger, if the scalar and pseudoscalar couplings are of the same order, while magnons allow access to the models where the mediator dominantly couples to the pseudoscalar currents of SM fermions. Next, we compare the reach of phonons and magnons to multipole models in Fig. 3; for the magnetic dipole and anapole models we expect the magnon reach to be better, and indeed it is. However, the phonon reach from quartz is sufficiently strong that, given the greater experimental challenges currently associated with magnon read-out, quartz should be considered a competitor for these models. Lastly, in Fig. 4, we compare theoretical reach in the (L · S)-interacting model, where magnons outperform phonons for sub-MeV DM with the same exposure; however, the (L · S)-interacting model is difficult to UV complete, and our calculation is perhaps somewhat an academic exercise that demonstrate aspects of the EFT.
We now discuss each benchmark model in turn.

A. Standard Spin-Dependent Interaction
For the standard SD interaction there is only one operator, O 4 , which generates the S response, and can excite both phonons and magnons in a magnetically ordered target. Here, only couplings to electrons (whose spins are ordered) are relevant, and we obtain, for the differential rates, In Fig. 1, we compare the phonon and magnon reach with YIG. As a technical note, in the absence of a DFT calculation for the crystal potential in YIG which is necessary for computing the  Table II from single magnon (red) and phonon (blue) excitations in YIG. The phonon rate is estimated in two ways, as discussed in the text, which lead to the solid and dashed curves, respectively. Since this model generates only the S response, magnons are seen to have better sensitivity than phonons. phonon eigenmodes, we estimate the rate in two ways. First, we carry out an approximate analytic calculation taking into account long-wavelength acoustic phonons, as explained in Appendix D.
This results in the dashed reach curve in Fig. 1, which is truncated at the DM mass for which the maximum momentum transfer reaches the edge of the 1BZ, so that the approximations we make cease to hold. Second, we borrow the crystal potential of Y 3 Ga 5 O 12 (YGG) which is publicly available [96]. YGG has the same crystal structure as YIG, with Fe replaced by Ga, and the phonon dispersions we obtain for YGG are very similar to those of YIG [97]. The resulting reach is shown by the solid blue curve in Fig. 1. We see from the figure that both estimates are in good agreement near m χ ∼ 10 −2 MeV, where acoustic phonons dominate, while including optical phonon contributions in the second approach improves the reach at lower and higher m χ .
We can understand these curves by estimating the rates using Eqs. (43) and (44). The q integrals are dominated by q max ∼ m χ v. As a result, Fixing R, the coupling plotted in Fig. 1 , and the advantage becomes more significant at smaller m χ (though the magnon curve hits the kinematic threshold at higher m χ due to the dispersion being quadratic).

B. Scalar Mediator Models
We next consider scalar mediator models with both scalar and pseudoscalar couplings. We take the mediator couplings to SM fermions to be proportional to their masses, g ψ ∝ m ψ (motivated by Higgs-portal hidden sector theories, see Ref. [98] for a recent review), and consider each of the four combinations of currents, which we denote by S × S, P × S, S × P and P × P . Among them, S × S (i.e. standard SI considered previously in Refs. [33,34,53,54]) and P × S can excite phonons via the N response, 8 while S × P and P × P can excite both phonons and magnons in a magnetically ordered target via the S response. However, similar to the standard SD interaction in Sec. III A, the phonon excitation rate will be suppressed relative to the magnon excitation rate, so we focus on the latter here. We obtain the following expressions for the differential rates defined in Eq. (41): Note that for the S × S and P × S models, screening effects have been taken into account by using g eff ψ in place of g ψ , as discussed around Eq. (14); the dielectric tensors ε ∞ of the phonon targets GaAs and SiO 2 are obtained from DFT calculations [56]. scalar mediator. The couplings to SM fermions are taken proportional to their masses, g p = g n = mp me g e , and we fix g χ g e = 10 −13 . Each curve is labeled with the model type as in Table II and the excitation type (phonon or magnon) that can probe each model. The phonon curves assume SiO 2 (solid) and GaAs (dashed) targets, and the magnon curves assume YIG (solid) and α-RuCl 3 (dashed) targets.
In Fig. 2, we plot the expected rate for each of the four coupling combinations, for a common value for the product of couplings, to illustrate the hierarchy between the rate from the different interactions. We have chosen to show the rate instead of projected reach here so that the general case where more than one types of interactions are present, it would be straightforward to rescale the curves to see which one is dominant. For example, if g ψ , we have the highest rate from phonon excitations via the S × S coupling, i.e. the standard SI interaction, as expected. On the other hand, if the couplings to SM fermions are dominantly pseudoscalar, ψ /g (S) ψ 10 7 , magnon excitations have better sensitivity than phonon excitations for the same exposure; this is one of the benchmark models considered previously in Ref. [57]. The hierarchy seen in Fig. 2, and also some main features of the curves, can be understood following Eqs. (43) and (44), as we now explain.
First consider the light mediator case, m V q (left panel of Fig. 2). For phonon excitations in the S × S and P × S models, since the couplings to all ions have the same sign, the rate is dominated by acoustic phonons. For q within the 1BZ, setting ω ∼ c s q, we obtain where ω min = c s q min . These are consistent with the m −1 χ and m −2 χ scaling of the green and purple curves for m χ up to ∼ MeV. Also, consistent with the figure, the ratio between them is For magnon excitations in the S × P and P × P models, we have again consistent with the m −1 χ scaling of the corresponding curves in Fig. 2 (though the YIG curves have a bump near MeV due to the gapped magnons starting to contribute, as discussed in Ref. [57], which slightly obscures the overall scaling with m χ ). Comparing the two models, we see that The heavy mediator case, m V q (right panel of Fig. 2), follows a similar analysis. All the q integrals are now peaked at q max ∼ m χ v, and we find  Projected reach on the multipole DM models listed in Table II, assuming dark photon-like couplings to SM particles: g p = −g e , g n = 0. The left panel shows the hierarchy of sensitivities of single phonon excitations, in GaAs and in SiO 2 , to the three multipole DM models, together with the SI interaction model for comparison. The center and right panels focus on the magnetic dipole and anapole DM models, respectively, and compare the phonon reach of GaAs and SiO 2 (via the N response), and the magnon reach of YIG (via the S response) and α-RuCl 3 (via both S and L responses); these models are best probed by magnons, though the phonon sensitivity with an optimal target like SiO 2 may be competitive.

C. Multipole Dark Matter Models
We now turn to the electric dipole, magnetic dipole, and anapole DM models in Table II. For comparison, we also include the SI interaction model with a vector mediator. Motivated by the kinetic mixing benchmark, we assume the mediator couples to electric charge, g p = −g e , g n = 0, and is much lighter than the smallest momentum transfer, m V eV. The SI and electric dipole DM models generate O 1 and O 11 at leading order, respectively, both of which induce only the N response, which can be probed by single phonon excitation. The differential rates are Eq. (61) is in agreement with previous results in Refs. [33,34,54]. The magnetic dipole and the anapole DM models generate, in addition to N , also S and L responses, and can therefore be probed by both phonons and magnons. For single phonon excitation, we have Note that for an unordered/unpolarized target, F (ψ) L,ν = 0. For single magnon excitation, we have which extend the results in Ref. [57].
A comparison of the phonon reach in these models is shown in the left panel of Fig. 3. The center and right panels of Fig. 3 zoom in on the magnetic dipole and anapole DM models, respectively, and compare the reach of phonon and magnon excitations.
We can carry out a similar analysis as in the previous subsections to understand the main features in Fig. 3. For single phonon excitation in GaAs and SiO 2 , we keep only the F (ψ) N,ν terms in the Σ ν (q) formulae above, and note that, as in the SI case discussed previously in Refs. [33,34,53,54], the DM-ion couplings, being proportional to N p − N e = Q ion , have opposite signs for oppositely charged ions, so the optical phonon modes with ω ∼ q 0 give the dominant contributions.
Next, let us compare the reach of different target materials, and via phonons versus magnons.
For phonon excitations, the factor in parentheses in Eq. (67) reproduces the "quality factor" identified in Ref. [34], up to O(1) factors we have dropped here. It captures the material properties that determine the sensitivity to the SI model with a dark photon mediator, and is the quantity to maximize in order to optimize target choice. For example, SiO 2 has a quality factor that is about 80 times that of GaAs, which explains its significantly better reach, by almost an order of magnitude on the coupling g χ g e , as seen in Fig. 3 (and also previously in Ref. [34]).
For magnon excitations for the magnetic dipole and anapole DM models, we have considered YIG, which probes only the S response, and α-RuCl 3 , which probes both S and L. Since for these models, DM couples to the linear combination 2S e + L e -the spin of an elementary particle has a Landé g-factor of 2 -the additional L response that α-RuCl 3 has does not qualitatively improve the sensitivity. Indeed, we see from ω , which is the phonon quality factor with the dimensionful parameters normalized in a way close to Ref. [34].
Its values are typically O(10 −7 -10 −5 ), with GaAs and SiO 2 residing on the lower and higher ends of the interval, respectively. We find where m cell is for the target for magnon excitations, and we have substituted the numbers for YIG in the last equation. We see that, for the magnetic dipole and anapole DM models, magnons are indeed more sensitive than phonons, though choosing high phonon quality factor targets, such as SiO 2 with Q ∼ 10 −5 can approach the magnon sensitivity. Up to O(1) factors, this is consistent with the center and right panels of Fig. 3.

D. (L · S)-Interacting Dark Matter
We finally consider the (L · S)-interacting DM model, which induces N , S and L ⊗ S responses.
Taking the mediator to couple only to electrons for simplicity, we obtain the differential rates: (1 + κ)F In the absence of magnetic order, F S,ν = 0. Also, unless κ is tuned to be very close to −1, we do not expect the F (e) L⊗S,ν term in Eq. (72) to dominate -the total spin-orbit coupling vanishes for full shells, and is otherwise often suppressed by crystal fields, especially for lighter elements. Thus, while an interesting feature of this model, the coupling to L · S does not suggest a better probe than searching for phonon excitations via the N response with already proposed target materials.
In Fig. 4, we include only the F for a heavy mediator (m V q). These equations explain the scaling of the phonon curves in The magnon reach curves for YIG and α-RuCl 3 can be understood in a similar way. We have for a light mediator (m V q), and for a heavy mediator (m V q). In contrast to the phonon case, the reach on g χ g e (g χ g e

IV. CONCLUSIONS
We have formulated an EFT framework for computing direct detection rates via single phonon and magnon excitations for general DM interactions, and illustrated its application with a set of benchmark models, listed in Table II, that cover a wide range of possibilities for a spin-1 2 DM particle interacting with SM fermions ψ = p, n, e (proton, neutron and electron). The procedure consists of first matching a relativistic DM model onto a set of NR effective operators, listed in Table III, and then matching these operators onto lattice degrees of freedom, including particle numbers N ψ , spins S ψ , orbital angular momenta L ψ and spin-orbit couplings L ψ ⊗ S ψ for the ψ = p, n, e particles in an ion. These define the four types of crystal responses and enter the rate formula for single phonon excitation, while a subset of them -S e and L e -also enter the rate formula for single magnon excitation.  Table II for our benchmark models), one simply replaces the ionic expectation values X ψ lj by F (ψ) X,ν defined in Eq. (39) (for ψ = p, n, e and X = N , S, L, L ⊗ S) or E X,ν defined in Eq. (40) (for ψ = e and X = S, L), squares the expression and takes the DM spin average. This gives the differential rates Σ ν (q), which are then substituted into Eqs. (41) and (42) for the total rate of single phonon or magnon excitation.
The set of crystal responses that we have identified point to various possibilities of optimizing detector target choice. However, a general observation from our calculations in Sec. III is that, among the four types of crystal responses, N ψ and S ψ , which are associated with point-like degrees of freedom, tend to dominate the rate, compared to the composite responses L ψ and L ψ ⊗ S ψ . This implies that, purely from the point of view of maximizing the rate, exotic materials with orbital orders or strong spin-orbit couplings are not necessary for improving the reach to a broad class of DM models.

Meanwhile, as phonon DM experiments focused on crystal targets, such as SPICE (Sub-eV Polar
Interactions Cryogenic Experiment), which is part of the TESSERACT (Transition Edge Sensors with Sub-EV Resolution And Cryogenic Targets) project [88], move into R&D, it is important to note that their discovery potential extends well beyond the simplest models with spin-independent interactions studied previously. As we showed in Sec. III, phonon excitations have broad sensitivity to light DM models. Perhaps surprisingly, with judicious choice of target material, phonon excitations may even be competitive with magnon excitations for some DM models where the latter is expected to have a parametrically higher rate, such as the magnetic dipole and anapole DM models. Given the greater challenges associated with single magnon detection relative to phonons, this is encouraging for phonon-based experiments in the near term.
In this appendix, we review the procedure of decomposing a Dirac fermion field ψ in the NR limit. Consider the following unperturbed relativistic Lagrangian: In free space, we would expand the ψ field in plane waves multiplied by the usual u, v spinors satisfying the free particle Dirac equation. Here, we allow the presence of an external gauge potential A µ = (Φ, A), which may not be a small perturbation. For example, if ψ is an electron in a crystal, it is bound by the electromagnetic potential from the ions, and the bound state wavefunctions are very different from plane waves. Generally, we can expand the ψ field in the basis of energy eigenstate wavefunctions. Dropping the antiparticle part, we have where the c-number u I spinors satisfy Writing with Ψ I , Θ I two-component wavefunctions, we see that Eq. (A3) is solved by This immediately leads to Eq. (8), repeated here for easy reference: where k = −i∇ − A, ε = i∂ t − Φ, and ψ + (x, t) = I e −iε I t Ψ I (x)b I with ε I = E I − m ψ . The prefactor has been chosen such that the NR field ψ ± 's kinetic term is normalized at leading order as in Eq. (11).
In the NR limit, |Θ I | |Φ I |. The large component Ψ I satisfies At leading order, we replace 1 2m ψ +ε I −Φ(x) → 1 2m ψ , and recover the NR Schrödinger equation: Corrections to this equation can be incorporated order by order if needed.

Appendix B: Projection of Angular Momentum Operators
In this appendix, we detail the steps that lead to the numbers λ S,j = − 1 3 , λ L,j = − 4 3 in the case of α-RuCl 3 , following the projection of angular momentum operators S e , L e in Eq. (32).
The formation of effective ionic spins S j = 1 2 is due to the combined effect of crystal fields and spin-orbit coupling [82]. First, octahedral crystal fields split the five degenerate 3d orbitals ( = 2) of Ru 3+ into two higher-energy e g orbitals and three lower-energy t 2g orbitals with an effective orbital moment eff = 1. The energy difference between the e g and t 2g orbitals is O(eV), rendering the (unoccupied) e g orbitals irrelevant for the discussion. For the t 2g orbitals, spin-orbit coupling further splits these eff = 1 states into j eff = 3 2 and 1 2 . With five 3d electrons, the lower-energy j eff = 3 2 states are fully occupied, while the higher-energy j eff = 1 2 Kramers doublet is occupied by a single electron -it is this electron that contributes to the magnetic order. Therefore, the goal is to project the angular momentum operators S, L (dropping subscript e from here on for simplicity) onto the j eff = 1 2 subspace. The first step is to project L onto the t 2g subspace. The t 2g states are denoted by d yz , d zx , d xy .
The angular part of their wavefunctions are linear combinations of spherical harmonics Y m =2 (θ, φ) (see e.g. Ref. [99]); equivalently, these t 2g states are linear combinations of | , m states with = 2: To compute P t 2g L P t 2g , with the projection operator we make use of the familiar formulae where L ± = L x ± iL y , and obtain, for the matrix representation in the |d yz , |d zx , |d xy basis: These might not look familiar, but they are nothing but = 1 angular momentum operators in the |p x , |p y , |p z basis, which is related to the | , m basis with = 1 by [99] The angular momentum operators in this basis read Comparing Eq. (B4) and (B6), we see that L acts as an effective angular momentum with = 1 on the t 2g subspace: The second step is to combine this effective orbital angular momentum eff = 1 with the electron's spin s = 1 2 . This follows the standard angular momentum addition, and we obtain, for the j eff = 1 2 states: where the coefficients are Clebsch-Gordan coefficients. It is now straightforward to project L eff and S onto the j eff = 1 2 subspace: We see that both L eff and S are proportional to J eff = σ 2 (identified as the total ionic spin as discussed above) when acting on the j eff = 1 2 subspace. So finally, we obtain When the velocity dependent rate Γ(v), given by Eq. (41), is convoluted with the incoming DM's velocity distribution f χ (v) to yield the total rate, Eq. (42), we encounter the following scalar, vector and tensor velocity integrals: where v χ = v − q 2mχ , and ω q = q · v − q 2 2mχ . From the expressions of differential rates Σ ν (q) throughout Sec. III, it should be easy to see how these integrals emerge. Note that for velocityindependent interactions, only the scalar integral g 0 appears [33,43,54].
As we now show, all three velocity integrals above can be evaluated analytically for a boosted and truncated Maxwell-Boltzmann distribution, which we assume in this work: where and we take v 0 = 230 km/s, v esc = 600 km/s, v e = 240 km/s. For all the target materials considered in Sec. III, the rates are insensitive to the direction of v e . The analytic results obtained here are key to efficient rate calculations, as they reduce the six-dimensional integral d 3 v d 3 q to just a three-dimensional integral d 3 q, which we then compute numerically.
We then obtain Next, the vector integral g 1 can be decomposed as The first term can be computed by shifting v → v − v e as before, but this time the integrand also depends on the azimuthal angle φ: dv v e −v 2 /v 2 0 = v * q g 0 (q, ω) , wheren 1 ,n 2 are orthogonal unit vectors in the plane perpendicular to q. Plugging in the definition of v * in Eq. (C7), we obtain g 1 (q, ω) = ω qq − (1 −qq) · v e g 0 (q, ω) .
They follow from q · v χ = ω q , and can be easily checked using the explicit expressions above.

Appendix D: Estimation of Single Phonon Excitation Rate in YIG
In this appendix, we explain the analytic estimation that results in the dashed curve in Fig. 1.
For the standard SD interaction considered in Sec. III A, the single phonon excitation rate is where F (e) S,ν (q) = j=Fe 3+ e −W j (q) e iG·x 0 j q · * ν,k,j 2m j ω ν,k S e j .
See Eqs. (41), (45) and (39). For YIG, ν runs from 1 to 240. However, since DM has same-sign couplings to all the Fe 3+ ions (and zero couplings to the other ions), we expect acoustic phonons to give an O(1) contribution to the total rate at low momentum transfer. Further, the dot product q · * ν,k,j in F (e) S,ν singles out the longitudinal acoustic branch, ν = 3, which has the following general properties at low momentum [100]: where c s is the longitudinal acoustic sound speed. Also, we can set G = 0, k = q, and W j 0 at low q. Therefore, where v min (q) ≡ q 2m χ + c s .
Now we can write the total rate per unit target mass in terms of the commonly used η function, defined by The result is This is the formula we use to estimate the single phonon excitation rate in YIG in Sec. III A.
The material parameters are c s = 7.2 km/s [101], S cell = 10, m cell = ρ T Ω, with ρ T = 4.95 g/cm 3 , Ω = 990.683Å 3 . The analytic expression for the η(v min ) function for the Maxwell-Boltzmann distribution of Eq. (C4) can be found in e.g. Ref. [33]. Since the η function has support up to q max 2m χ (v e + v esc ), we cut off the dashed curve in Fig. 1 at the m χ value for which q max reaches π Ω 1/3 , roughly the edge of the 1BZ.