Effective Field Theory for Dark Matter Absorption on Single Phonons

Single phonon excitations, with energies in the $1-100 \, \text{meV}$ range, are a powerful probe of light dark matter (DM). Utilizing effective field theory, we derive a framework to compute DM absorption rates into single phonons starting from general DM-electron, proton, and neutron interactions. We apply the framework to a variety of DM models: Yukawa coupled scalars, axionlike particles (ALPs) with derivative interactions, and vector DM coupling via gauge interactions or Standard Model electric and magnetic dipole moments. We find that GaAs or $\text{Al}_2\text{O}_3$ targets can set powerful constraints on a $U(1)_{B-L}$ model, and targets with electronic spin ordering are similarly sensitive to DM coupling to the electron magnetic dipole moment. Lastly, we make the code, \textsf{PhonoDark-abs} (an extension of the existing \textsf{PhonoDark} code which computes general DM-single phonon scattering rates), publicly available.


I. INTRODUCTION
absorption on single phonons has only been computed for the kinetically mixed dark photon DM model [70,89]. The purpose of this work is to generalize the computation of DM absorption to any DM model which has a Yukawa-like interaction Lagrangian of the form, where g is a perturbatively small coupling constant, ϕ is the DM field, Ψ ∈ {e, p, n} is either an electron, proton, or neutron Standard Model (SM) field, and O is an operator. The simplest example of this interaction Lagrangian is when O = 1, and ϕ is a scalar field coupling to, e.g., electrons; then Eq. (1) is simply L ⊃ g ϕēe. Eq. (1) can also apply to vector DM, V µ , when O has a matching Lorentz index, e.g., L ⊃ g V µē γ µ e and O = γ µ . This can be further extended when the operator O is allowed to contain momentum (derivatives) acting on the DM and SM fields, such as for an axionlike particle (ALP), a, with derivative coupling to electrons: L ⊃ g ∂ µ aēγ µ γ 5 e. In momentum space, this simply corresponds to O = −i q µ γ µ γ 5 , where q µ is the four-momentum of the ALP field.
To compute the DM absorption rate for a general DM model we utilize the self-energy formalism developed in Refs. [67,[90][91][92]  In addition to the general operators encompassed by Eq. (1), we consider targets with and without fermionic spin ordering, e.g., (anti) ferromagnets. We find that targets with different spin orderings can be sensitive to different DM models. For example, in the absence of spin ordering, ALP DM does not couple to phonons via the derivative coupling; however, in a spin ordered target phonons can be excited.
We make the code used to compute the absorption rates for all targets, PhonoDark-abs , publicly available here . This program complements PhonoDark [77,93], which was developed to compute general DM-single phonon scattering rates.
The paper is organized as follows. In Sec. II we provide a theoretical framework to compute the DM absorption rate into single phonons via a general, Yukawa-like interaction in the form of Eq. (1). This derivation will proceed in three steps: first, in Sec. II A we write the absorption rate in terms of self-energy diagrams using the optical theorem. Second, in Sec. II B we compute the phonon contribution to these diagrams, whose imaginary part leads to single phonon absorption, in terms of vertex Feynman rules, which are then derived in detail in Sec. II C. Before computing the DM absorption rate into single phonons, we compute the single phonon contribution to the dielectric function in Sec. III A, which serves both as a cross check of the formalism, and verification of the first principles calculations of the target properties. In Secs. III B -III E we compute the DM absorption rate into single phonons for different DM models and show the projected constraints.
Specifically, we consider scalar DM (Sec. III B), the derivative coupling of ALP DM (Sec. III C), and two models of vector DM, one which couples to the SM vector currents (Sec. III D), and another which couples to the SM particles electric and magnetic dipole moments (Sec. III E).

II. GENERAL SINGLE PHONON ABSORPTION RATE
The purpose of this section is to compute DM single phonon absorption rates due to Yukawalike interactions in the form of Eq. (1). The derivation proceeds in three steps. First, in Sec. II A, we use the optical theorem to write the DM absorption rate in terms of in-medium self-energies.
Second, in Sec. II B, we write the phonon contribution to the self-energies in terms of crystal form factors, F , describing how the DM field couples to phonons. These form factors will depend on the properties of the ions at each lattice site, e.g., the number of protons or the electronic spin.
Third, in Sec. II C, we detail how the form factors are derived using non-relativistic effective field theory (NR EFT) starting from Eq. (1).

A. Absorption Rate In Terms of Self-Energies
Following Refs. [67,[90][91][92], we start by deriving the DM absorption rate in terms of in-medium self-energies, Π. While this formalism was originally developed to compute DM absorption on electrons, it is agnostic about the underlying crystal degrees of freedom and can be similarly applied to the case of phonon absorption. The optical theorem states that the absorption rate of the λ th polarization of the DM field, ϕ, is, where m ϕ is the DM mass, Π λ ϕϕ is the self-energy between two ϕ particles of the λ th polarization. However, when ϕ can mix with the photon field, A, this introduces complications since one must include diagrams mixing ϕ and A. The sum of all the diagrams can be succinctly written in terms of 1PI diagrams by first going to the "in-medium" basis, where the fields are no longer mixed. When ϕ and A are perturbatively coupled, the leading order self-energy of the in-medium "DM-like" state,φ, is given by, where we have introduced the self-energies polarization components defined as, Π λλ ′ ≡ −e λ µ Π µν e λ ′ * ν , where e λ µ are the polarization vectors of the subscripted fields, e.g., Π λλ ′ ϕA = −e λ ϕ,µ Π µν e λ ′ * A,ν . The polarization vectors are defined to diagonalize the Π ϕϕ and Π AA self-energies, i.e., Π λλ ′ AA = Π λ AA δ λλ ′ , and in general will differ from the vacuum longitudinal and transverse polarization vectors. Given Eq. (3), the total DM absorption rate is given by, and the DM-polarization averaged absorption rate, per target mass and exposure time, R, is then, where ρ ϕ is the DM density, taken to be 0.4 GeV/cm 3 , ρ T is target density and n is the number of polarizations of the DM field.
The absorption rate in Eq. (5) can be simplified further if we assume that the photon selfenergy is independent of polarization, i.e., Π λ AA ≃ Π AA , which is true in the isotropic limit. To simplify the following analysis we assume isotropy, and leave a study of anisotropic corrections to future work. In this case, the sum over η can be performed exactly using the completeness relation η e µ A,η e ν, * A,η = −g µν . Moreover, the Ward identity, Q µ Π µν = 0, demands the time components to be q-suppressed relative to the spatial components. Therefore, −g µν Π λµ ϕA Π νλ Aϕ ≈ Π λi ϕA Π iλ Aϕ , and Eq. (5) simplifies to, This leads to our final absorption rate expressions for (pseudo) scalar DM, R S , and vector DM, When deriving the absorption rate for vector DM we have once again used the completeness relation of the polarization vectors to perform the sum over the DM polarizations. Notice how in both the scalar and vector case we were able to remove the dependence on the photon and DM polarizations, such that the problem of computing the absorption rate has now shifted to deriving the spatial components of the in-medium self-energies.
The in-medium self-energies receive contributions from both electronic excitations, Π el , and phonon excitations, Π ph . For example, at one loop, the following graph topologies will contribute: where the first diagram represents the phonon contribution while the last two diagrams encode the contribution from electronic excitations. Here we will assume that the electronic band gap is much larger than the energy of phonon excitations, such that no electron excitation can go on shell at energies relevant for DM absorption into phonons. As a result, Im [Π el ] ≃ 0 and Therefore, in order to compute the absorption rates given in Eq. (7), in addition to the phonon contribution to the self-energies one also has to compute the real contribution from electronic excitations. The electron contribution to in-medium self-energies has been extensively studied in Refs. [57,67,90], therefore, in the following we will focus on the novel phonon contribution and use the values of Re [Π el ] derived in these previous works.

B. Phonon Contribution To Self-Energies
The phonon contributions to the diagrams in Π ϕϕ , Π ϕA and Π AA can all be understood from the same diagram: where Q µ = (ω, q) is the incoming four-momentum, Φ and Φ ′ can be either ϕ or A (the diagram 1 Any directionalq dependence within the rate is averaged over, i.e., we computeR = 1 4π dΩqR(q) for both scalar and vector DM models. Furthermore we take q = ωv0, where v0 = 230 km/s. inherits the Lorentz indices of the field, e.g., Π AA → Π µν AA ), and the central double line is a phonon propagator. To compute the diagram the phonon propagator and vertex rule are needed. The phonon propagator, D νk (ω), is given by [94], where ν, k index the phonon branch and momentum within the first Brillouin zone (1BZ), respectively, ω νk is the phonon energy, and γ νk is the phonon linewidth, or inverse of the phonon lifetime.
Assuming, for now, that the left and right vertex rules are given by, M Φ,νk , M * Φ ′ ,νk , respectively, the self-energy is, where V = N Ω is the volume of the target, N is the number of unit cells, and Ω is the unit cell volume. Analogous to the self-energy, M Φ,νk will inherit the Lorentz indices of the field Φ.
In order to separate the part of the vertex that depends on the structure of the UV Lagrangian from the part that is common among different UV interactions, we parameterize the vertices for scalar and vector fields as respectively, where S is a scalar field, V is a vector field, and the double line indicates a phonon.
We will summarize the meaning of each term here, and provide a detailed derivation of M Φ,νk in Sec. II C. A factor of √ N has been factored out to cancel the 1/N in Eq. (11), as well as the momentum conservation factor δ q,k . j indexes the ions in each unit cell, and the sum indicates that all the ions in the unit cell contribute to the generation of a phonon. T jνk is defined as the phonon transition matrix element, u ℓj is the displacement operator, |νk⟩ = a † νk |0⟩ is a single phonon state at ν, k, x 0 ℓj is the equilibrium position of the j th ion in the ℓ th unit cell, m j is the mass of the j th ion, and ω νk , ϵ jνk are the phonon energies and polarization vectors, respectively. The phonon transition matrix element, T jνk , is coming from a |0⟩ → |νk⟩ transition in the vertex. Lastly, the form factors F S, j and F µ V, j are vectors that contain the information about how the scalar, S, or vector field, V couple to the displacement operator, u ℓj , via macroscopic properties of the ion, e.g., total number of electrons. These form factors may contain contributions from each fermion type, e, p, n at the lattice site j, and therefore can be further decomposed as F j = ψ F jψ , where ψ ∈ {e, p, n}. The detailed derivation of these form factors is given in Sec. II C, and summarized in Table I for the DM models of interest.
Substituting Eqs. (10) and (14) into Eq. (11) gives the final expression for the phonon contribution to the self-energies, where the S, V index on F follows from directly from the Φ, Φ ′ particle type, e.g., when computing Π ϕϕ for scalar ϕ, F S, j should be used. We have also simplified the phonon propagator since for absorption kinematics, q ≪ ω, ω νq ≈ ω ν .

C. Dark Matter-Phonon Interaction Form Factors
In Sec. II A we wrote the single phonon absorption rate in terms of electron and phonon selfenergies, and in Sec. II B we wrote the phonon self-energies, Eq. (16), in terms of some form factors, F j . These results were independent of both the UV Lagrangian and target material, whose dependence manifests within the aforementioned form factors. In this section, we derive these form factors starting from the UV Lagrangian in Eq. (1). Schematically, the derivation proceeds as follows, The first step "NR EFT" (Sec. II C 1) (non-relativistic effective field theory) reduces the UV Lagrangian, L UV , written in terms of the four-component Ψ fields, to the NR Lagrangian, L NR , written in terms of two-component fields, ψ, which describe the electron, proton, and neutrons in the target. The NR expansion is appropriate when the energy and momentum transfers are much smaller than the fermion masses, which is certainly the case for absorption processes. shown in the left column. The leading order form factor is shown for targets with no spin ordering (middle column), and targets with spin ordering (right column). For spin-1 DM we write the i index on the vector F V, j to avoid confusion with the µ index. Explicitly, F 0i V,j (F 0 V,j ) is the left component inside the parentheses, and F ki V,j , the i th component of the vector F k V,j , is the right component. Dashed lines indicate negligible, higher order responses. Note the ψ index on the form factor has been dropped from N, S for simplicity.

No Spin Ordering Spin Ordering
The "0" components of the spin-1 DM form factors are related by the Ward identity, Q µ F µ V, j = 0, where Q µ = (ω, q) is the incoming DM four-momentum.
While the NR Lagrangian is written directly in terms of the particles constituting the target, it does not contain any information about the target itself. In the case of a crystal, the target state is simply a lattice of ions at positions x ℓj = x 0 ℓj + u ℓj , where x 0 ℓj is the equilibrium position of the ion, u ℓj are the displacement operators, and each site indexed by ℓj, where ℓ indexes the unit cell, and j indexes the ion inside the unit cell. This information is added in the second step, labeled "⟨⟩ ℓj " in Eq. (17), or "Target Expectation Value" (Sec. II C 2), which transforms the DM-ψ interaction Lagrangian to DM coupling to the lattice properties at each site, e.g., where n jψ is the number density of ψ particles on the j th site, by summing the expectation values at each lattice site.
In the last step, "|0⟩ → |νk⟩" in Eq. (17), or "Form Factor Calculation" (Sec. II C 3), the form factors, or vertex rules in Eqs. (12) and (13), are derived from the interaction Lagrangian, L(u ℓj ), which is written in terms of the displacement operators. This step is simply computing quantum mechanical matrix elements of the transition between an initial state with no phonon and incoming DM particle, to a final state containing no DM particle and a single phonon indexed by νk.
These three steps are described in the following subsections. While the procedure to connect the UV Lagrangian to the form factors is the same for each operator O in Eq. (1), the details can differ. Therefore to avoid repetition, at each step in the derivation we begin with a general discussion, and then provide an example calculation for vector DM, V , coupling the (spatial part) of a vector current, L UV (Ψ) = gV µΨ γ µ Ψ ⊃ −g V iΨ γ i Ψ. The form factors for all of the DM models considered in Sec. III can be found in Table I.

NR EFT
The purpose of finding the NR limit of a UV Lagrangian is to isolate the dynamics of the two-component field, ψ, which satisfies the Schödigner equation, within the Lagrangian containing two two-component fields inside Ψ. Our starting point is the Dirac Lagrangian, where D µ is the gauge covariant derivative, and m Ψ is the mass of the fermion. The problem becomes more clear after a change of variables, Ψ → e −im Ψ t Ψ which transforms Eq. (18) to, where P ± = (1 ± γ 0 )/2 are projection operators. Eq. (19) describes the dynamics of two twocomponent fields, P ± Ψ, where P + Ψ is massless, and P − Ψ is massive. For this reason, we will refer to P + Ψ as the "light" field and P − Ψ as the "heavy" field. If there were no terms in Eq. (19) which mixed the heavy and light fields then there would be no problem; the dynamics of the two fields are decoupled.
However, the γ i D i term mixes the heavy and light fields, and therefore to isolate the dynamics of the light field, a procedure to remove the heavy field needs to be performed. This is the fundamental problem of NRQED/QCD [95][96][97], and there are many different approaches. We will give a summary of two methods that have been utilized in the context of DM direct detection, Refs. [57] and [90], and refer the reader to these references for more details.
Ref. [57] used the "equation of motion" (EOM) method, which is the most physically intuitive.
One simply solves for the EOM of the heavy field in terms of the light field, and then substitutes the heavy field back into the Lagrangian. This generates a Lagrangian which only depends on the light field. The NR limit of the interaction Lagrangian in Eq. (1) is also readily found; to first order in the DM coupling, one can simply substitute the heavy field which satisfies the EOM when no DM field is present. Including the O(g) dependence in the heavy field EOM only introduces extra O(g 2 ) terms.
While physically straightforward, when integrating out the heavy field extra time derivatives enter, which require careful field redefinitions to keep canonically normalized fields. Another approach, used in Ref. [90] and known more generally as a Foldy-Wouthuysen (FW) transformation [98][99][100][101][102][103], avoids this by removing the mixing with consecutive field redefinitions at each order in 1/m Ψ . That is, one finds n Hermitian operators, {X 0 , X 1 , ..., X n−1 } such that, removes all heavy/light field mixing to O(m −n Ψ ). One can show that the operators, where the Tr is performed over the 2 × 2 block diagonal matrix, and ψ, in the Dirac basis, is the upper two components of Ψ on the right-hand side of Eq. (20). Eq. (22) gives the general form of the first step in Eq. (17). Applying Eq. (22) to the example UV DM Lagrangian, and keeping terms leading order in both 1/m Ψ and absorption kinematics (q ≪ ω), yields, which comes solely from the second term in Eq. (22) using [γ i , γ 0 γ j ] = 2δ ij .

Target Expectation Value
Given the NR Lagrangian in terms of the electron, proton, and neutron fields, ψ, the DM-phonon interaction Lagrangian is simply a sum over the expectation value at each lattice site, These expectation values will then be written in terms of the target properties at each site, e.g., where n jψ is the number density of the ψ field. This step is analogous to the EFT calculation performed in Ref. [77] for DM-single phonon scattering. However, in Ref. [77] it was the scattering potential, V, which was written as a sum of the scattering potential at each lattice site, V = ℓj ⟨ V ⟩ ℓj . Further simplifications were made assuming that the scattering potential only depends on x ℓj in the same way as Eq. (25), i.e., V = ℓj V ℓj (x − x ℓj ). This simplified calculations by allowing the x ℓj dependence to be factored out in the Fourier transform of the scattering potential, Since we are only concerned with single phonon excitations in the q ≪ ω limit here, we perform a different simplification of these expectation values than in Ref. [77] by focusing on the terms that are linear in u ℓj . The other terms will not enter the matrix element calculations of the form factors performed in Sec. II C 3, and avoids the exponential dependence on x ℓj . Additionally, the derivation performed here will keep expectation values that are O(ω) which were subdominant in for the scattering EFT and implicitly dropped when assuming V = ℓj V ℓj (x − x ℓj ).
As an example, the linear order in u ℓj term in Eq. (25) is, An additional simplification can be made when we consider that these expectation values multiply the DM field inside the Lagrangian L NR . Therefore we can integrate by parts and move the derivative acting on the number density to the DM field, and convert to momentum space with Similar simplifications can be performed for the spin density, s i jψ , where the factor of two enters from the definition of spin, S = σ/2.
More complicated operators can be simplified using Ehrenfest's theorem. For example, consider ⟨ψ † k i ψ⟩ ℓj , Ehrenfest's theorem states that, and therefore, where we have used the Schrödinger equation, H 0 ψ = i∂ 0 ψ, and, similar to q previously, ω represents a time derivative acting on the DM field. With one exception that will be discussed later, the operators, 1, σ i and k i are the only operators needed to compute the form factors for all the models discussed here. Furthermore, we assume that there are no background vector gauge fields, With these target expectation values computing the example DM-phonon interaction Lagrangian from Eq. (23) is trivial,

Form Factor Calculation
The last step in the derivation is to identify the form factor from the DM-phonon interaction Lagrangian. This is done by computing the matrix elements from Eqs. (12) and (13), respectively, where the δ/δϕ (δ/δV µ ) simply removes the scalar field, ϕ (vector field, V ) from the interaction vertex, leaving only a function of the phonon operators, u ℓj . It is easiest to understand this formula in practice, and similar simplifications hold for all DM-phonon interaction Lagrangians.
Consider the example L(u ℓj ) in Eq. (31), in this case we have where we have written the displacement operator matrix element in terms of the phonon transition matrix element with Eq. (14). The integral in Eq. (34) can be related to the total particle number, and the sum over ℓ enforces momentum conservation within the crystal, ℓ e i(q−k)·x ℓ = N δ q,k .
After these simplifications, one can isolate the form factor, where the k index corresponds to the spatial part of µ, and i corresponds to the index of the vector, , what gets contracted with T jνk in Eq. (12).

III. APPLICATIONS
In this section, we apply the formalism developed in Sec. II to compute DM absorption rates into single phonon excitations for a variety of targets and DM models. Our focus will be on  Table I. Specifically, for each DM model, we derive the phonon contribution to the in-medium self-energy by substituting the form factors given in Table I into Eq. (16). These self-energies are then substituted into Eq. (7) to compute the total absorption rates.
The phonon transition matrix elements, T jνk , are computed from first principles in two steps.
First, using first principles density functional theory (DFT) [104] calculations within VASP [105][106][107][108][109], the lattice is relaxed to its equilibrium position and the equilibrium positions, x 0 ℓj , are found. Each ion is then displaced from its equilibrium position and the forces on the ion are computed, which generates the spring constants between the ions. VASP is also used to compute the highfrequency dielectric constant, ε ∞ . These three pieces of information are contained in the POSCAR, FORCE SETS, and BORN files output from VASP (or similar DFT software). Second, these files are then input to the phonopy program [110] which diagonalizes the system and calculates the phonon energies, ω νk , and eigenvectors, ϵ jνk .

PhonoDark-abs
is used to compute all absorption rates shown here. PhonoDark-abs performs the second step (calling phonopy) internally, and therefore one simply needs to supply the DFT input files (POSCAR, FORCE SETS, and BORN), for any target material, to compute the absorption rate. PhonoDark-abs is publicly available here .
Before computing DM absorption rates, in Sec. III A, we apply our formalism to compute the long-wavelength dielectric function, ε(ω). This calculation serves a variety of purposes. First, it is needed to compute the screened contribution to the DM absorption rates in Eqs. (7), as the dielectric function is related to the self-energy of the photon via Π AA = ω 2 (1−ε(ω)). Second, it provides a cross-check of our formalism since it allows us to compare our results with previous calculations of the dielectric function in terms of the low-energy electronic and phononic responses [111]. Finally, by comparing our results with experimentally measured values of the dielectric we can tune the phonon widths, γ ν , in the phonon propagator in Eq. (10).
In Secs. III B -III E we compute the DM absorption rates into GaAs, Al 2 O 3 (sapphire), and SiO 2 (quartz) targets. GaAs was the first studied due to the simple structure of its unit cell [70], while Al 2 O 3 is desirable for its large number of resonances and directionality [71] and ready availability in terms of fabrication of ultra-pure single crystals. Both of these targets will be used in the TESSERACT experiment [78]. SiO 2 has been previously identified as an optimal target in terms of reach to light DM scattering off phonons [72]. The DFT input files for GaAs, Al 2 O 3 , and SiO 2 are identical to those used in previous works [72,76,77,81,112]. These targets have no spin ordering, i.e., the fermion spins are not periodically aligned, S ℓj ψ ̸ = S j ψ . Note that spin ordering includes both ferromagnetic and anti-ferromagnetic ordering. Therefore we also consider a magnetically ordered target, FeBr 2 . This magnetic target was chosen because the first-principles calculations of its phonon properties are publicly available [114][115][116][117]; its purpose is to serve as an example calculation, not promote this specific target as a detector concept. The Fe 2+ ion has a magnetic moment of µ ≈ 3.9 µ B [118], where µ B is the Bohr magneton.
Assuming that the 3d electrons are orbitally quenched, as is common for transition metal electrons due to crystal field effects, the spin quantum number is S ≈ 1.8. While FeBr 2 is ferromagnetically ordered within the unit cell [116], it is anti-ferromagnetically ordered in adjacent cells [119]. Since this target is only meant to serve as an example, we will treat it as a ferromagnet and take S e ≈ [0, 0, 1.8] on the Fe 2+ lattice site.

A. Dielectric Function
The long-wavelength (q ≈ 0) dielectric function, ε(ω), receives contributions from both the electron and phonon degrees of freedom in a crystal. This can be understood simply in terms of two contributions to the photon self-energy, Π AA = Π el AA + Π ph AA , where Π el AA is the electronic response, and Π ph AA is the phononic response. Well below the band gap, the electronic contribution is directly related to the "high-frequency" dielectric constant, Π el AA = ω 2 (1 − ε ∞ ), encoding the response of the electrons if the lattice ions were not allowed to move, or "clamped". Using this, and the definition of the total dielectric function defined in terms of the total photon self-energy, Π AA = ω 2 (1 − ε(ω)), we can write the dielectric function as With the help of Table I and Eq. (16), we can then compute Π ph AA to obtain where the sums over ψ have been simplified in terms of the total electric charge of the ion, ψ g ψ N jψ = e N j,p − e N j,e ≡ e Q j , since g p = −g e = e in the QED Lagrangian, L QED ⊃ −eA µē γ µ e+eA µp γ µ p, following the convention in Ref. [120] where e = −|e|. The factor of 3 comes from taking the isotropic limit and averaging over the spatial components, Π AA, ph = Π ii AA, ph /3. This agrees with the standard result in, e.g., Ref. [111], providing a validation of the formalism. 2 In Fig. 1 we compare the imaginary part of the dielectric function (top row) and energy loss function (ELF), Im [−1/ε(ω)] (bottom row), computed from first principles with Eq. (38), to measured data from Ref. [89] for the non-spin ordered targets GaAs, Al 2 O 3 , and SiO 2 . The computed dielectric function is shown for different assumptions about the phonon widths, γ ν ∈ {10 −3 ω ν , 10 −2 ω ν , 10 −1 ω ν }. Smaller widths correspond to a larger resonance peak and smaller off-resonance behavior, and vice versa for larger widths. We find that the measured data can be well reproduced with phonon widths in this range, with slight shifts to the exact locations of the resonances. More sophisticated models of the widths as a function of energy could further improve these fits. For the results shown in Secs. III B -III E we use γ ν = 10 −2 ω ν .
It is known that the absorption rate on phonons of some DM models, e.g., the kinetically mixed dark photon [70,89] and ALPs (in an external magnetic field [113]) can be related to the measured ELF shown in Fig. 1. That first principles calculation can reproduce the measured ELF further validates the first principles approach of computing single phonon absorption rates in these models as studied in, e.g., Ref. [76]. In Fig. 4 we explicitly compare the constraints on the kinetically mixed dark photon model from the measured and calculated ELFs, whose differences are due to the differences shown in Fig. 1. 2 Eq. (38) is derived assuming that the electronic wave functions do not distort under ionic motion. These effects can be incorporated by loosening the assumption, ⟨ψ † γ 0 ψ⟩ ℓj ≈ n ℓj,e (x − x ℓj ) + δn ℓj,e , where δn ℓj,e ≡ iq i δZ ik j u k ℓj , and δZ ik j ≡ Z ik j − Qj, where Z ik j are the "Born effective charges". The spatial component, ⟨ψ † γ i ψ⟩ ℓj , follows a similar simplification by the Ward identity. This adds a form factor to the photon-electron coupling, δF µi je = (q m δZ mi j , ω δZ ki j ), replacing Qj T i jνq → Z ik j T k jνq in Eq. (38). See Ref. [72] for more details. in small band-gap targets, i.e., Al superconductors ("Al-SC", purple) [52,57], and spin-orbit coupled targets ZrTe 5 (turquoise) [67]. Shaded regions correspond to constraints from fifth force experiments (teal) [121,122] and stellar cooling bounds from red giants ("RG", pink) [92] and white dwarfs ("WD", orange) [123].

B. Scalar DM
The first model we consider is scalar DM, ϕ, whose couplings to electrons and nucleons are given by the Lagrangian, Using and Π el ϕA = Π el Aϕ . Since only the imaginary component of Π ϕϕ enters in the absorption rate given in Eq. (7), we have ignored the electron contribution to Π ϕϕ as it is purely real. This will also apply to the self-energies discussed in Secs. III C -III E. The electron contribution to Π ϕA was derived in Refs. [57,67], and is given by the last term in Eq. (41). 3 While the general absorption rate is given by substituting Eqs. (40) and (41) into Eq. (7), it is illuminating to study specific combinations of the coupling constants. For example, if the g Ψ coefficients are "photon-like", i.e., g p = −g e = g, g n = 0, then all the self-energies are proportional to Π AA (assuming an isotropic target), indicating that the total absorption rate can be written in terms of the ELF: Since the absorption rate can be written in terms of the ELF, it can also be related to the dark photon absorption rate, R dp = (ρ ϕ /ρ T ) κ 2 Im [−1/ε(ω)] [89]. Therefore the constraints on g can be related to the constraints on the mixing parameter, κ, of the dark photon model: When the couplings are proportional to the particle masses, g Ψ = g m Ψ , the ions are shaken in-phase, and therefore optical, or out-of-phase, oscillations are not excited. Therefore, as optical phonons are the only ones that can match DM absorption kinematics, this leads to a vanishing absorption rate. This effect, sometimes referred to as the "coupling to mass" effect [70,71,76,80], mathematically corresponds to the statement that in the absorption kinematics limit. While the cancellation is exact for couplings g Ψ ∝ m Ψ , it is also important even when the couplings are approximately proportional to the masses. For example, consider only coupling to the electron number on each site in GaAs, N Ga, e = 28, N As, e = 36.
Parameterically, one would expect F j ∝ N j e , however because the masses are m Ga = 69.7 u, m As = 74.9 u, when one subtracts off the contribution which vanishes due to the coupling to mass effect, , the form factors are roughly a factor of 10 smaller than N j e .
The scalar DM models affected by the coupling to mass effect are fairly generic. This is because both the proton and neutron masses are dominantly dependent on the same quantity, the QCD scale. Therefore scalar DM models which couple to the QCD field strength kinetic term [124], or the benchmark hadrophilic DM model [125] satisfy g p,n = g m p,n , and single phonon excitations will have limited reach due to the coupling to mass effect. Because of this, and since constraints on DM models with photon-like couplings can be trivially related to dark photon constraints, we focus on a DM model with only coupling to electrons, which suffers less from the coupling to mass effect, as discussed previously.
In Fig. 2 we compare the constraints on the electron coupling derived in this work to stellar cooling [92,123] and fifth force [121,122] constraints assuming no backgrounds and a kg · yr exposure. To facilitate the comparison with other conventions for the coupling constant, we show the constraints both on g e , defined in Eq. (39), and the commonly adopted parameterization

C. Axionlike Particle DM
The QCD axion [131][132][133][134] is one of the most theoretically motivated DM candidates since it also provides a solution to the strong-CP problem. The canonical QCD axion DM candidate, with an abundance set by the post-inflationary misalignment mechanism [135][136][137][138][139], is predicted to have a mass in the 10 −6 eV ≲ m a ≲ 10 −5 eV range, well below the scale of gapped phonon excitations in crystal targets. However, non-standard production mechanisms, as well as ALPs that do not necessarily solve the strong-CP problem, can prefer larger values of the axion mass .
Our focus will be on the derivative ALP couplings, 4 L ⊃ Ψ∈{e,p,n} Naïvely, the leading order term in the NR Lagrangian seems to be higher order than the "axion wind" term, ∝ g aΨΨ a q · σ/m Ψ . However, when evaluating the target expectation value of the "wind" term, an additional factor of q enters the form factor via Eq. (28). Therefore, the form factor for the "wind" term is order q 2 /m Ψ , which is much smaller than the form factor for the term in Eq. (45) (see Table I). More generally, when the leading order term in the Lagrangian 4 The Lagrangian in Eq. (45) is equivalent to L ⊃ − Ψ gaΨΨ aΨiγ 5 Ψ. Both forms of the Lagrangians give the same form factor in Table I, but one needs to expand to O(1/m 2 Ψ ) when taking the NR limit of the aΨiγ 5 Ψ Lagrangian. The equivalence can be shown explicitly using the relationship in Eq. (59). See Sec. III E for more details. respectively. Since these are not real targets, their purpose is to give an estimate of a target that does have non-zero spin ordering. As in Fig. 2, dashed lines correspond to projected constraints from absorption on electrons in small band-gap targets, i.e., Al superconductors ("Al-SC", purple) [52,57], and spin-orbit coupled targets ZrTe 5 (turquoise) [67]. In the left panel, the shaded gray region corresponds to constraints from solar axion searches with the XENONnT experiment [41], and the shaded light blue region corresponds to white dwarf ("WD") cooling constraints [123]. The red line corresponds to constraints from red giant ("RG") cooling [126,127], which have come under recent scrutiny [128]. In the middle and right panel, the shaded teal region corresponds to neutron star ("NS") cooling [129]. Tan lines correspond to the prototypical KSVZ and DFSZ QCD axion models [130], assuming 0.28 ≤ tan β ≤ 140 in the DFSZ model.
is O(q/m Ψ ), it is important to check if at next-to-leading order in the O(1/m 2 Ψ ) NR expansion there are terms which dominate. For example, a term of the form ωk/m 2 Ψ , where k is the fermion momentum, is dominant compared to the q/m Ψ term. This feature will also appear when discussing the dipole DM models in Sec. III E.
While the phonon contribution to the self-energies is straightforward to derive using the form factor in Table I, simplifying the electron contribution to Π aA in a spin ordered target requires an additional approximation. This can be understood physically: the high-frequency dielectric corresponds to the electronic response to an electric field, which couples identically to all electrons in the target. Therefore for any effect to be related to the high-frequency dielectric, it must affect all the electrons in the same way, i.e., all the electrons have the same spin such that the spin density is proportional to the number density, s i e =ŝ i e n e . However, while all the electrons must be spin polarized for an exact correspondence, Π aA may be approximately written in terms of ε ∞ as long as the electrons which give the dominant contribution to ε ∞ are spin polarized. Under these approximations the relevant self-energies are, and Π el aA = −Π el Aa . Similar to Sec. III B the absorption rate can be simplified with specific combinations of the coupling constants. For example, if the target is a ferromagnet, and g app = −m p g aee /m e , g ann = 0, then the absorption rate can be expressed in terms of the ELF (and dark photon absorption rate) as, In this case, the constraints on g aee can be related to the constraints on the dark photon coupling κ: g aee ∼ 10 −9 100 meV ω κ 10 −16 (49) Additionally, if ψ g aΨΨ S jψ /m ψ ∝ m j , then the coupling to mass selection rule in Eq. (44) applies and the absorption rate vanishes. While the coupling combinations are more contrived than for the scalar DM models, these two scenarios serve as benchmark points to understand different limits of the theory.
In Fig. 3 we compute the projections (assuming a kg· yr exposure and no backgrounds) for three models, where the only non-zero coupling is the one plotted. We compare these to stellar cooling bounds [123,126,127,129], the canonical DFSZ and KSVZ QCD axion model [130] predictions. be in similar targets with proton or neutron spin ordering, which can be achieved in the presence of a strong magnetic field [161].
We find that single phonon absorption is weaker than stellar cooling bounds for all couplings, especially g app and g ann . We note that the reach on g app and g ann is severely affected by the in GaAs (solid red), Al 2 O 3 (solid blue), and SiO 2 (solid green) targets utilizing single phonon excitations assuming a kg · yr exposure and no backgrounds. Projected constraints on the kinetically mixed dark photon model have also been shown in Ref. [89]; the purpose of the comparison here is to illustrate the good agreement between the first principles calculation performed here, and the data-driven approach (dotted lines) utilizing the ELF [89], also compared in Fig. 1. Dashed lines are projected constraints from targets utilizing electronic absorption: doped Si (orange) [61], Al superconductors ("Al-SC", purple) [52,57], and the spin-orbit coupled target ZrTe 5 (turquoise) [67] Fifth force constraints are from Ref. [162].
coupling to mass effect. If the proton or neutron spins could anti-align on different sites, to avoid the coupling to mass selection rule, the reach would improve. However this seems experimentally unfeasible. While the g aee constraint from Al 2 O 3 * is competitive with the XENONnT bounds [41], and nearly reaches the DFSZ band, the white dwarf [123] and red giant [126,127] cooling bounds are stronger by roughly an order of magnitude on resonance. However, recently there has been some uncertainty surrounding the stellar cooling bounds on g aee [128], which may re-open the parameter space. Absorption into magnons via the wind coupling [76] is still the dominant process in electron spin ordered targets. This is because the magnon response from the wind coupling does not suffer the extra q suppression that the phonon response does, as discussed previously. However, the strict selection rules governing that process [76] severely limits the number of useful modes in simple targets, especially in the absence of an external magnetic field, and single magnon readout is still a developing technology.

D. Vector DM: Gauge Coupling
Next, we consider vector DM, V , coupled to the SM fermion vector currents, generally arising from U (1) gauge theories. The self-energies are straightforwardly computed, and Π el V A = Π el AV . For g p = −g e = κe, and g n = 0, we recover the kinetically mixed dark photon model [163], where κ is the kinetic mixing parameter. As shown in Refs. [70,89], the absorption rate for this model is directly related to ELF shown in Fig. 1 We see that absorption into single phonons in any of the GaAs, Al 2 O 3 or SiO 2 targets can be far superior not only to fifth force constraints in the 30 meV ≲ m V ≲ 100 meV mass window, but also absorption into small gap electronic excitations.

E. Vector DM: Electric and Magnetic Dipole
The last DM models we consider are again vector DM models, but this time with a magnetic or electric dipole coupling to SM fermions. We will refer to these models as MDM and EDM models, respectively. These models were studied in Refs. [90,164,165] in the context of DM with mass m V ≳ keV, where it was shown that production could occur via the standard freeze-in mechanism [164] but is dominated during reheating due to the dimension five nature of the operators. Other production mechanisms could produce the DM nonthermally with the MDM/EDM couplings at much smaller masses [44][45][46][47][48][49], and therefore it is interesting to study the constraints on these models below the eV scale where constraints from electronic excitations begin [90].
The UV Lagrangians of the MDM and EDM models are, The coupling to the spatial part of the vector DM can be further simplified as, and the NR Lagrangians are, Two terms are kept in the NR limit of the MDM model, the first is the leading order response when the target is spin ordered, and the second is the leading order response when the target is not spin ordered. The NR limit of the EDM has multiple terms due to contributions of similar order from the NR limit ofΨσ 0i iγ 5 Ψ andΨσ ij iγ 5 Ψ in Eq. (56). Note that all terms in Eq. (58) will contribute at the same order in the form factors.
Before continuing to the self-energies we comment on the target expectation value of D 0 , D i , which is a bit subtle. We assume that the SM fermions are bound at the lattice site by the temporal component of the gauge fields. For example, the electrons are bound by the potential eA 0 , which is simply the electrostatic potential generated by the ion. Assuming that the protons and neutrons are bound by similar strong forces, then D 0 = ∂ 0 + iV, where V is the binding potential, e.g., the electromagnetic part of V is simply eA 0 . The NR Lagrangian of the SM fermions is then simply the Schrödinger equation with this potential, H = p 2 /2m Ψ + V. Furthermore, assuming that the target expectation value of the spatial part of the gauge fields vanishes (i.e., the spatial part of the gauge fields does not significantly impact binding) then D 0 , D i → H, k i , and therefore, which can then be straightforwardly written in terms of the displacement operator with Eq. (30).
In the case of no spin ordering, only the MDM self-energies are non-zero, where the M subscript denotes the MDM model (the E will subscript denotes the EDM model) and Π el M, V A = Π el M, AV . For targets with spin ordering both MDM and EDM have self-energies, and Π el M/E, V A = −Π el M/E, AV . In Fig. 5 we focus on models with only an electron coupling which were the focus of Ref. [90,164].
Projections are computed assuming a kg · yr exposure with no backgrounds, and the line labeled, "Therm." corresponds to the minimum coupling needed to not thermally produce the vector DM in a universe that reheats right before BBN. That is, we require Γ ≲ H at T = MeV, where H is the Hubble constant, Γ =n V σ E/M is the interaction rate,n V is the equilibrium number density of V particles, and σ E/M ∼ d 2 E/M [164]. In addition to the thermalization bound we also show projected constraints from an Aluminum superconductor, and ZrTe 5 target. Ref. [90] showed that the electronic absorption rate of the MDM model could be related to Im [ε(ω)], and therefore the constraints on g aee can simply be re-scaled accordingly. We see that spin ordering is crucial to be able to probe either of these models, since targets without spin ordering have no EDM response, and an MDM N response only at higher in 1/m Ψ . In these spin ordered targets single phonon excitations are able to probe new parameter space for both the MDM model.

IV. CONCLUSIONS
Single phonon excitations are an exciting avenue for direct detection of light DM with sub-eV thresholds. In Sec. II, using effective field theory (EFT) techniques we provided a framework for computing the DM absorption rate into single phonons starting from a fairly general UV Lagrangian (Eq. (1)). This complements previous work which computed general DM-single phonon scattering rates [70-72, 77, 79-82] and further illustrates the variety of DM models that can excite single phonon excitations. Then in Sec. III we applied this formalism to compute the DM absorption rate of five DM models (Secs. III B -III E) on spin ordered, e.g., FeBr 2 , and non-spin ordered, GaAs, Al 2 O 3 , SiO 2 targets. Additionally, in Sec. III A, we used the formalism to compute the dielectric function which allows for a direct comparison between first principles calculation and experimental data. In Fig. 1 we find good agreement between both the dielectric function, ε(ω), and the ELF, Im [−1/ε(ω)], in GaAs, Al 2 O 3 , and SiO 2 targets, indicating the reliability of the first principles calculations. Moreover, this comparison allows for a data-driven approach to set the only (theoretically) free parameter, the phonon mode widths, γ ν .
In addition to providing a theoretical framework to compute general DM absorption rates into single phonon excitations, we developed PhonoDark-abs . PhonoDark-abs is an extension of PhonoDark [77,93], which computes general DM-single phonon scattering rates (see Refs. [72,77,81,112,166,167] for examples), and numerically computes the DM absorption rate for any target material, given the input DFT files discussed in Sec. III. Currently, PhonoDark-abs can reproduce all the results shown here, and future work will further extend its capabilities and integrate it with PhonoDark completely. PhonoDark-abs is publicly available here .
Using PhonoDark-abs, we find that, assuming a kg · yr exposure and no backgrounds, single phonon excitations in GaAs, Al 2 O 3 , and SiO 2 can probe new parameter space when DM is the gauge boson of a broken U (1) B−L symmetry (Fig. 4), and when DM couples to the electron magnetic dipole moment in spin ordered targets, e.g., the ferromagnetic FeBr 2 (Fig. 5). 5 For the latter projected constraints the spin ordering is crucial; without spin ordering the target response is much higher order and therefore normal targets, e.g., GaAs and Al 2 O 3 , project rather weak constraints. While the projected constraints for the scalar, Fig. 2, and ALP, Fig. 3, DM models coupling to electrons are competitive with targets utilizing small band-gap electronic transitions, e.g., Al superconductors [52,56,57], ZrTe 5 [67], and doped Si [61], strong stellar cooling and fifth force constraints are still superior in this parameter space. Furthermore the projected constraints for ALP DM coupling to protons and neutron spin with hyperpolarized targets are much weaker than the neutron star cooling.
The theoretical framework here may be useful for other other collective excitations, e.g., magnons [73][74][75][76][77]168]. Formulating the absorption rate in terms of self-energies and using NR EFT, as done in Secs. II A and II C 1, has the advantage of being independent of the internal excitations, and therefore may be used to understand general DM absorption into magnon excitations.