Multi-Channel Direct Detection of Light Dark Matter: Theoretical Framework

We present a unified theoretical framework for computing spin-independent direct detection rates via various channels relevant for sub-GeV dark matter -- nuclear recoils, electron transitions and single phonon excitations. Despite the very different physics involved, in each case the rate factorizes into the particle-level matrix element squared, and an integral over a target material- and channel-specific dynamic structure factor. We show how the dynamic structure factor can be derived in all three cases following the same procedure, and extend previous results in the literature in several aspects. For electron transitions, we incorporate directional dependence and point out potential daily modulation signals in anisotropic target materials. For single phonon excitations, we present a new derivation of the rate formula from first principles for generic spin-independent couplings, and include the first calculation of phonon excitation through electron couplings. We also discuss the interplay between single phonon excitations and nuclear recoils, and clarify the role of Umklapp processes, which can dominate the single phonon production rate for dark matter heavier than an MeV. Our results highlight the complementarity between various search channels in probing different kinematic regimes of dark matter scattering, and provide a common reference to connect dark matter theories with ongoing and future direct detection experiments.

To cover a broader mass range, electrons have been considered as an alternate pathway to detecting light DM. A variety of targets have been studied, including noble gas atoms which can be ionized with O(10 eV) energy deposition, semiconductors where electron transitions can happen across O(eV) band gaps [23][24][25][26][27][28][29][30][31][32], as well as systems with O(meV) gaps like superconductors [33][34][35] and Dirac materials [36][37][38]. Electron transitions can potentially extract all of DM's kinetic energy, and thus constitute a more efficient search channel than nuclear recoils. For example, semiconductor targets can probe DM masses down to O(MeV).
When the energy deposition is below the band gap, electron transitions are kinematically forbidden. However, there are condensed matter systems with collective excitations that can couple to the DM. For example, collective excitations in superfluid helium (phonons and rotons) are sensitive to O(meV) energy depositions, especially via phonon pair production [39][40][41][42]. In a crystal target, the active degrees of freedom below the electronic band gap are acoustic and optical phonons -quanta of collective oscillations of atoms/ions. Direct excitation of single phonons in crystals has been recently proposed as a new search channel for light DM [43,44]. Optical phonons typically have energies of O(10-100 meV), and can be excited by DM as light as O(10 keV). Acoustic phonons are gapless and, assuming an O(meV) detector threshold, can also probe DM down to O(10 keV).
All these detection channels do not exist in isolation. Depending on the DM mass and couplings to Standard Model (SM) particles, it may either cause nuclear recoils, or induce electron transitions, or excite phonons in the same target material. Thus, when designing direct detection experiments, an important consideration should be to search for DM across multiple channels in parallel. The kinematic interplay between several channels that we will discuss in detail is illustrated in Fig. 1. with ω ph the phonon energies, above which the rate is suppressed by the Debye-Waller factor. We see that a GeV-mass DM can be probed by all three channels; a 10 MeV DM is out of reach in conventional nuclear recoil searches, but can be searched for via electron transitions in semiconductors and single phonon excitations in crystals; a sub-MeV DM cannot even trigger electron transitions in eV-gap materials, but can still be detected via single phonon excitations.
On the theory side, most of the basic ingredients for the rate calculation are known. However, they have been developed in separate contexts, and at first sight look very different for different detection channels. In our opinion, it would be much more convenient to have a common theoretical framework for all these calculations. This will not only facilitate the comparison of target materials across various existing and proposed search channels, but also provide the necessary calculation tools when new search channels are considered in the future.
It is the purpose of this paper to lay out such a formalism, focusing on spin-independent (SI) DM interactions. 1 As we will see, for each detection channel, the calculation is factorized into a particle physics model-specific part and a target response-specific part. The latter is encoded in a dynamic structure factor, to be computed by quantizing the particle number density operators in the Hilbert space of the excitations under study. We show how this is done in three cases -nuclear recoils, electron transitions and single phonon excitations. While the first two are relatively simple, and our calculation is mostly a formal rederivation of known results, the phonon calculation presented here contains new aspects. Our general framework allows us to derive single phonon excitation rates for arbitrary SI couplings from first principles, such as phonon excitation by coupling to electrons.
In addition to deriving general rate formulae in this unified framework, we also aim to clarify various conceptual and technical issues in direct detection calculations, and present new results that highlight some previously overlooked experimental prospects. For nuclear recoils, we clarify the range of validity of the standard calculation. For electron transitions, we go beyond the commonly made isotropic approximation. In fact, there exist simple materials with large anisotropies. As an example, we consider boron nitride (BN) with a hexagonal crystal structure, and O(eV) band gap, and show that the expected rate can vary by ±(10 -40) % during a day as the DM wind enters from different directions. Such daily modulation signals have been pointed out previously for electron transitions in graphene [46], carbon nanotubes [47] and Dirac materials such as ZrTe 5 and BNQ-TTF [37,38], and for single phonon excitations in sapphire [44] where they help distinguish signal from background. Here we show that also O(eV) band gap three dimensional semiconductors, like BN, can exhibit daily modulation. 2 Finally, for single phonon excitations, we extend the rate calculation to DM heavier than an MeV, where the DM's de Broglie wavelength is shorter than the typical lattice spacing, and Umklapp processes can contribute significantly. We point out an interesting interplay with nuclear recoils, and demonstrate the complementarity between the two 1 The idea of treating various detection channels in a common framework was previously advocated in Ref. [45], where the focus was on DM nuggets. Here we follow the same spirit and develop a formalism for calculating direct detection rates for general DM models, assuming a point-like DM particle. 2 See also Refs. [48][49][50] for proposals that take advantage of direction-dependent threshold effects. channels. We also compute the phonon production rate for generic couplings to the proton, neutron and electron, extending previous results for dark photon mediated interactions.
We focus on the theoretical framework in the present work; in a companion paper [51], we apply the results presented here to carry out a comparative study of many candidate target materials, and discuss strategies to optimize the search across multiple channels. We also note that there are additional detection channels beyond those we discuss in detail here (e.g. excitation of molecular states [52][53][54], multi-excitation production in superfluid helium [39][40][41][42]), which have been pursued and can be studied in the same framework.

SCATTERING
In a direct detection event, a non-relativistic DM particle, χ, deposits a certain amount of energy, and triggers a transition |i → |f in the target system. We assume the target system is initially prepared in an energy eigenstate |i (usually the ground state) and, as usual, treat the incoming and outgoing DM particles as momentum eigenstates |p , |p , with p = m χ v, p = p − q.
For a given incoming velocity v and momentum transfer (from the DM to the target) q, the energy deposition is Here and in what follows, we denote q ≡ |q|, where q is the momentum 3-vector. Note that for given DM mass m χ , the energy deposition is bounded by the parabola, ω q ≤ qv max − q 2 /2m χ , as shown in Fig. 1. Applying Fermi's Golden Rule and summing over the final states, we obtain the rate: where δĤ is the interaction Hamiltonian, |p, i = |p ⊗|i , |p , f = |p ⊗|f . We take the quantum states to be unit normalized unless specified otherwise, e.g. p|p = i|i = 1.
The DM part of the matrix element can be evaluated universally at the Born level: where V is the total spatial volume, V(x) is the effective scattering potential felt by the DM, and V is its Fourier transform. We focus on SI couplings in the present work, in which case the scattering potential takes the form 3 Here n p , n n , n e are the proton, neutron and electron number densities in the target, and V p , V n , V e are the respective scattering potentials from a single particle located at the origin. We thus have Note that for SI interactions, V ψ (−q) = V ψ (q) (ψ = p, n, e) are functions of only the magnitude of q. In vacuum, they simply coincide with 2 → 2 scattering matrix elements M χψ (q) familiar from standard quantum field theory calculations. In the target medium, however, they may receive corrections due to screening effects (see Sec. II B). We can define (momentum-dependent) effective in-medium couplings f p , f n , f e to account for screening effects, while the corresponding couplings in the vacuum Lagrangian are denoted by f 0 p , f 0 n , f 0 e . We can write where is the vacuum matrix element for DM scattering off any of the constituent particles (proton, neutron or electron) with unit coupling. The total scattering potential is then Let us rewrite this equation as follows: Depending on the DM model and the process under consideration, we will factor out either M χn or M χe , and define a target form factor, F T (q), composed of contributions from protons, neutrons and electrons, as the quantity in brackets. In other words, we have 3 More generally, DM interactions can be classified by nonrelativistic effective operators [55][56][57][58]. The SI interaction we focus on here is the leading operator if generated without velocity suppression. Other operators result in spin and/or velocity dependence of the scattering potential V(x), and may be probed via additional detection channels beyond those considered in this work. For example, DM coupling to the electron spin can excite magnons in solid state systems with magnetic order [59]. We leave a general effective field theory study of light DM direct detection to future work.
where M stands for M χn or M χe . We can further factor out the q dependence of M, which can only come from the mediator propagator for tree-level scattering: The reference momentum transfer is conventionally chosen to be q 0 = m χ v 0 (with v 0 the DM's velocity dispersion) for DM-neutron scattering, and q 0 = αm e for DM-electron scattering.
The factorization in Eq. (10) is a key component of the formalism. From the target-independent particle-level matrix element M, we define the reference cross sections: where µ denotes the reduced mass. These coincide with the total cross sections of DM-neutron and DM-electron scattering in the heavy mediator case. On the other hand, F T is target specific, from which we define the dynamic structure factor: 4 which encapsulates response of the target to DM couplings to the proton, neutron and electron.
Let us highlight the following regarding the dynamic structure factor S(q, ω).
• S(q, ω) depends on the distribution of constituent particles p, n, e in the target system via n p , n n , n e , which in turn depends on the nucleus types and electron wavefunctions. It is therefore target material specific.
• S(q, ω) also depends on the active degrees of freedom in the target system via the choice of |f , which in turn determines how F T (q) should be quantized. It is therefore excitation (detection channel) specific. 4 Here we adopt a slightly different normalization convention compared to Ref. [45]. The right hand side of Eq. (14) here is identified with 2π Ω S(q, ω) in Ref. [45], where Ω is the primitive cell volume.
• If only one of the constituent particles p, n, e is responsible for the transitions |i → |f , S(q, ω) is DM model independent. Otherwise it depends on ratios (but not the overall strength) of the couplings f 0 p , f 0 n , f 0 e .
• For any given DM mass m χ and incoming velocity v, only a slice in the (q, ω) space, ω = ω q , is probed in the scattering process. The parabolic boundary of kinematic region for each m χ in Fig. 1 is the envelope of these slices for all v directions for fixed magnitude of v.
Finally, to obtain the total rate per target mass, we average over the DM's initial velocity, multiply by the number of DM particles in the detector, and divide by the detector mass, giving where ρ T is the target mass density, ρ χ is the local DM energy density, and f χ is the DM's velocity distribution in the target rest frame. A common choice for f χ is a truncated Maxwell-Boltzmann (MB) distribution boosted by the Earth's velocity with respect to the galactic rest frame, In the calculations presented in this paper, we take ρ χ = 0.4 GeV/cm 3 , v 0 = 230 km/s, v esc = 600 km/s, v e = 240 km/s.
In addition to the total rate, it is often useful to know the differential rate with respect to the energy deposition onto the target ω. This simply requires inserting delta functions into the integrals to pick out the contributions with ω = ω q : To summarize, we have the following algorithm for computing the rate for a given detection channel.
• First, identify the initial and final states |i , |f according to the type of excitation.
• Next, quantize F T (q) in terms of the relevant degrees of freedom such that it acts on the target Hilbert space to induce the transitions |i → |f .
• Then, compute the transition matrix element f |F T (q)|i , and thus the dynamic structure factor S(q, ω) via Eq. (14).
We will carry out this procedure for each detection channel in the next three sections. Before doing so, let us discuss some technical details regarding the phase space integration and in-medium effects.

A. Phase Space Integration
We see from Eqs. (15)-(20) that once the dynamic structure factor S(q, ω) is known, we need to perform a six-dimensional integral over v and q to obtain the event rate R. The integration gives familiar results in the special case of isotropic target response, but is more complicated in the general anisotropic case. We now discuss the two cases in turn.
a) Special case: isotropic target response. If f |F T (q)|i 2 = f |F T (q)|i 2 , as is the case for nuclear recoils, the only dependence on the direction of q is from the δ-function, where θ qv is the angle between q and v. Integrating over the angular variables, we have where The velocity integral then gives where These results are familiar from the standard nuclear recoil calculation [60], and have also been used in previous electron transition calculations, where the target response has been assumed to be isotropic. Note that they hold for any DM velocity distribution f χ (v). In the case of the MB distribution in Eq. (17), the η function can be evaluated analytically, giving We see that five of the six integrals have been done analytically, and we are left only with a one-dimensional integral over q (which can also be done analytically in the case of nuclear recoils).
b) General case: anisotropic target response. Generally, crystal targets are not fully isotropic, as the crystal structures break rotation symmetries. This implies that, for a terrestrial detector, since the DM wind comes in from different directions at different times of the day, there can be daily modulation in the detection rate. While the existence of this effect is well-known [23,27,43], it has been calculated only recently in the contexts of single phonon excitations [44] and electron transitions in Dirac materials [37,38], where the energy deposition is O(meV). In Sec. IV A, we calculate this effect for the first time in electron transitions in an O(eV) gap target.
When f |F T (q)|i 2 depends on the direction of q, the six-dimensional integral generally does not admit a simple analytical solution. To proceed, we first evaluate the velocity integral and define [44,61] The rate can then be written in terms of this g(q, ω) function as For general velocity distributions f χ , we still have to evaluate a six-dimensional integral, which is a numerically intensive task. However, for the commonly assumed MB distribution, Eq. (17), the g(q, ω) function can be evaluated analytically, giving where Thus, only the three-dimensional integral over q needs to be done numerically (in addition to other integrals that may be encountered in the evaluation of the dynamic structure factor).

B. In-Medium Effects
In the case of a vector or scalar mediator, in-medium effects can cause screening and affect direct detection rates. They must be taken into account when deriving the target response F T (q) (and hence the dynamical structure factor S(q, ω)) when present. While the treatment of in-medium effects has been discussed in various contexts [34,36,37,62], we review it here for completeness.
In particular, we derive the screening factors f ψ (q)/f 0 ψ (ψ = p, n, e) in this subsection. For nonrelativistic systems relevant for direct detection that we focus on here, only electrons can contribute significantly to screening when the energy deposition is above phonon frequencies (ω O(100 meV), corresponding to m χ O(100 keV)), as nuclei are too heavy to respond. At lower frequencies that match energy depositions in phonon excitation processes, there is additional screening in an ionic (polar) crystal due to relative motion of ions. However, as we will see in Sec. V B, the ions' response should be included in the source term in Maxwell's equations in order to be quantized in terms of phonon modes. Thus, also in this case, we consider only electron contributions to in-medium effects. 5 Consider a vector mediator A , and suppose the vacuum Lagrangian takes the form where J µ ψ =ψγ µ ψ (ψ = p, n, e). Here the first line is standard electromagnetism, the second line is the dark sector Lagrangian, and the third line contains A couplings to SM particles. We assume |f 0 ψ | 1, and consistently keep terms only at linear order in these couplings. Because the electron current J µ e couples to the linear combination A µ +κA µ , with κ = −f 0 e /e, as opposed to just A µ , the in-medium photon self-energy Π µν (q) implies the following terms in the momentum space quantum effective action, As in Ref. [37], we can project Π µν onto the three polarizations, (whereq = q/|q|, andê ⊥ is a unit vector perpendicular to q), and diagonalize the 3 × 3 matrix to find the canonical modes. It is worth noting that the polarization vectors satisfy As a result, in the vacuum limit where Π µν = g µν − q µ q ν /(q α q α ) Π and the photon propagator is proportional to In an isotropic medium, and the photon propagators are proportional to 1 q α qα−Π T,L . Generically, for an anisotropic medium, we need to simultaneously rotate A and A into a polarization basis where K is diagonal. In this basis, the quadratic part of the effective action can be diagonalized for each polarization by where Π is an eigenvalue of K. In the A, A basis, the propagators are proportional to Dark matter scattering is mediated by both A and A . Taking both into account, we obtain the following effective interaction: From the last equation, it is clear that the current A couples to contains a screened component and term, which is proportional to the electromagnetic current, gets screened by a factor of q α qα q α qα−Π , whereas the second term is unaffected.
In the special case of a dark photon that kinetically mixes with the SM photon, Eq. (31) follows from diagonalizing the kinetic terms, and κ is equal to the kinetic mixing parameter. In this case, and the coupling to neutrons is not screened. As a final example, a hadrophobic A has f 0 p = f 0 n = 0, resulting in an unscreened DM coupling to protons (which originates from the A-A mixing).
The screening factor q α qα q α qα−Π can be expressed in terms of the dielectric matrix ε(q, ω) by solving the following set of equations for Π µν [34,37]: Note that the three-dimensional quantities are defined by σ = σ i We obtain the following solution: Projecting Π µν onto polarization components, we obtain: We can see explicitly that in the isotropic limit, ε ∝ 1, so K L± = K ∓± = 0, and K LL , K ±± are identified as Π L , Π T , respectively. In this case, K = diag(Π L , Π T , Π T ), and the familiar relations are reproduced. Beyond the isotropic limit, in general one has to diagonalize the K matrix as discussed above. However, assuming anisotropies are not large, the calculation is simplified in the case of nonrelativistic scattering. Here, the currents involved (J µ χ , J µ e , etc.) have velocity suppressed spatial components, so the dominant contribution comes from the polarization that is almost longitudinal, for which Π K LL up to small corrections. As a result, the screening factor in Eq. (40) becomes Now it is straightforward to read off the screening of DM couplings from Eq. (39): In what follows, we will often drop the argument q and just write f p , f n , f e for simplicity. Finally, we note that in the scattering limit, q ω, a scalar mediator has the same coupling as the longitudinal component of a vector mediator, so the same screening factors in Eq. (50) apply.
To close this subsection, we comment that in-medium screening affects different channels differently. Nuclear recoils happen at high enough momentum transfer where ε can be approximated as unity, so f ψ f 0 ψ . For electron transitions, the situation depends on the band gap. For atoms, insulators and semiconductors with O(eV) or larger band gaps, ε approaches unity when q 2π/a ∼ O(keV) [65]; for smaller q, the full ε(q) can be fitted to experimental measurements or calculated using advanced electronic structure techniques. For small-gap systems such as superconductors and Dirac semi-metals, it is important to keep the full energy-momentum dependence in ε(q, ω). For example, in a (super)conductor, ε ∼ λ 2 TF /q 2 at low q, where λ TF ∼ O(keV) is the Thomas-Fermi screening parameter, resulting in significant screening [34]. In contrast, in a Dirac semi-metal, ε approaches a constant at low q, so sensitivity to dark photon mediated scattering (and also dark photon absorption) is much stronger [36,37]. For phonon excitations, screening from electrons should also be accounted for, as we discuss in Sec. V B.

III. NUCLEAR RECOILS
We now apply the general framework of the previous section to the case of nuclear recoils and reproduce familiar results. For simplicity we shall first assume only one type of nucleus is present, with proton number Z and atomic mass number A, and later generalize to the case of multiple nucleus types with non-degenerate {Z N }, {A N }.
To begin, we assume the nuclei do not interact with each other, so the Hilbert space of the target system, which contains ρ T V /m N nuclei, is a direct product of ρ T V /m N single nucleus Hilbert spaces. We will discuss the validity of this standard assumption in Sec. III A. The target system is prepared in the initial state with k i = 0. In the final state |f , one of the |k i J 's is replaced by |k f J with k f = 0. We can write these states in terms of nucleus creation operators: As usual, we have the canonical commutation relations Now we need to quantize F T (q) = 1 f n f p n p (−q) + f n n n (−q) + f e n e (−q) (53) in terms of nucleus creation and annihilation operatorsb † ,b. Obviously, the electron coupling does not contribute, so we drop the last term. The proton and neutron number densities, on the other hand, can be related to the nucleus number density n N , if we assume elastic scattering (no transition between nuclear states): where n 0 p,n are the proton and neutron number densities around a single nucleus at the origin. Therefore, where f N ≡ f p Z + f n (A − Z), the DM-nucleus coupling in the q → 0 limit (where DM interacts with all nucleons coherently). F N (q) is a nuclear form factor that deviates from unity only for q above the inverse nucleus radius. A commonly used form factor is the Helm form factor [66], where r n 1.14 A 1/3 fm, s 0.9 fm. We can thus write F T (q) in terms ofb † ,b via To obtain the dynamic structure factor, we evaluate the matrix element, and sum over final states, which amounts to summing over the scattered nucleus J (simply multiplying by ρ T V /m N ) and integrating over the final momentum V d 3 k f /(2π) 3 . Therefore, where we have regulated the delta function by We can now reproduce the familiar results for the differential rate. Assuming the nuclear form factor is isotropic, F N (q) = F N (q), as is the case for the Helm form factor in Eq. (56), we can apply Eq. (24) and obtain where η(v min ) is given by Eq. (25) and v min = q 2µ χN in the present case. It is now easy to generalize these results to the case of more than one nucleus type: where N runs over the inequivalent nuclei in the target (e.g. N = Ga, As for GaAs).

A. Validity of the Nuclear Recoil Calculation in Crystal Targets
A key assumption we have made in the derivation above is that the nuclei in the target do not interact with each other (hence the factorization of the Hilbert space). In a crystal target, however, the nuclei are not free, but interact with the neighboring nuclei in the crystal structure.
The justification of treating the nuclei as free particles initially at rest lies in the fact that in the classical limit, the hard scattering process is instantaneous and local. In this case, the nuclei interactions affect only the subsequent secondary processes. For example, secondary phonons can be produced, which allows the energy deposition to be shared by many nuclei.
On the other hand, as detector thresholds are pushed to lower energies, at some point we would get into the quantum regime, where the finite duration and spatial extent of the scattering invalidate the free nuclei assumption. We can make a quick estimate on when this happens from the uncertainty principle. The time scale for the hard scattering to happen is ∼ 1/ω. This should be compared to the intrinsic time scale for atomic vibrations in a crystal, 1/ω ph , with ω ph the phonon energy. The instantaneous interaction approximation in the standard nuclear recoil calculation is valid when the energy deposition is much higher than the energies of all phonon modes, i.e.
ω ω max ph (validity condition for nuclear recoils in crystals) .
An alternative way to reach the same conclusion is the following. Within the length scale 1/q, the DM should see the nucleus as a plane wave for the nuclear recoil calculation to hold. Since the spatial extent of the nucleus wavefunction in a harmonic potential is ∼ (m N ω ph ) −1/2 , we need q (m N ω max ph ) 1/2 . Using the kinematic relation ω = q 2 2m N , we arrive at the same condition as Eq. (63).
To summarize, in crystal targets, the nuclear recoil calculation is valid for energy depositions much higher than the phonon energies, which are typically O(10 − 100) meV. This explains the truncation of the C, Si, Ge, Cs nuclear recoil lines at low ω in Fig. 1. At lower energy depositions, the target Hilbert space does not factorize into individual nuclei, but instead contains single phonon and multi-phonon states as energy eigenstates, and the direct detection rate calculation proceeds differently. We discuss single phonon excitations in Sec. V, which will be the relevant processes when detector thresholds reach the 10-100 meV regime in the future. In the intermediate energy regime -above the single phonon energies yet below the validity range of nuclear recoils -direct multi-phonon production should be considered, which we plan to investigate in future work.

IV. ELECTRON TRANSITIONS
We next consider electron transitions. The initial state can be written as whereĉ † I are electron creation operators, with I running over all occupied electron states (energy eigenstates). Our normalization convention is such that {ĉ I ,ĉ † I } = δ II , so the electron states are unit-normalized. The final states are labeled by I 1 , I 2 , where one of the electrons has transitioned from I 1 to an unoccupied state I 2 : |f =ĉ † I 2ĉ I 1 |i .
The relevant piece in F T (q) is simply where the creation and annihilation operators are for momentum eigenstates, and satisfy {ĉ k ,ĉ † k } = (2π) 3 δ 3 (k − k ), etc. As discussed in Sec. II B, the screening factor is for a vector or scalar mediator.
The dynamic structure factor is therefore where we have usedĉ † I 1 |i =ĉ I 2 |i = 0, and that the anticommutators are just numbers. To evaluate the anticommutators, we expand the energy eigenstates in terms of momentum eigenstates: where ψ I (k) is the momentum space wavefunction, which satisfies the orthonormality condition The dynamic structure factor in Eq. (70) applies for any target system where DM scattering can trigger electron transitions -atoms, crystals, superconductors, Dirac materials, etc. -once the energy levels and wavefunctions are known. In what follows, we examine the case of periodic crystals in more detail. Here, the energy eigenstates of an electron are Bloch waves labeled by a band index and a wavevector within the first Brillouin zone (1BZ), e.g.
where G 1 runs over all reciprocal lattice vectors. Note that the state labeled by i 1 , k 1 has Fourier components of k 1 plus any reciprocal lattice vector. The coefficients u i 1 (k 1 + G 1 ) are normalized as G 1 |u i 1 (k 1 + G 1 )| 2 = 1. The dynamic structure factor now becomes where the prefactor 2 comes from summing over contributions from degenerate spin states, and the sums over the final state quantum numbers k 1,2 have been replaced by integrals in the continuum limit. As in Ref. [27], we define a crystal form factor for the transition i 1 k 1 → i 2 k 2 with an Umklapp G. This simply encodes the wavefunction overlap, summed over all Fourier components consistent with momentum conservation. The dynamic structure factor can now be written more concisely as Note that we have again used the identity (2π) 3 δ 3 (0) = d 3 x e i0·x = V . The material-specific quantities appearing in S(q, ω) are the electron band structures (energy eigenvalues E i,k ) and Bloch wavefunction coefficients u i (k + G). They can be computed by density functional theory (DFT) methods which we discuss more in our companion paper [51].
Finally, performing the phase space integration, we obtain the total rate per target mass: where The g(q, ω) function, the mediator form factor F med , the screening factor f e /f 0 e and the crystal form factor f [i 1 k 1 ,i 2 k 2 ,G] are given by Eqs. (27), (12), (67) and (74), respectively. This generalizes the formula derived in Ref. [27] to account for possible anisotropies in the target response.

A. Target Anisotropies and Daily Modulation
The simplest crystal targets that have been considered for direct detection via electron transitions, like silicon and germanium, are quite isotropic. As a result, the rate is essentially independent of the direction of the incoming DM's velocity. However, this is not the case for materials with large anisotropies in the electron band structures or wavefunctions. For terrestrial experiments, as the target rotates with the Earth, the DM wind comes in from different directions at different times of the day, resulting in a daily modulation of the rate. This is on top of the annual modulation signal expected due to the variation of the average DM velocity as the Earth orbits around the Sun [26,27]. If observed, it would be a smoking-gun signature of DM that is distinct from possible backgrounds. Our rate formula Eq. (76) incorporates directional information, and is well-suited for calculating the daily modulation signal.
As an example target, we consider hexagonal boron nitride (BN), shown in Fig. 2. The numerical calculation of electron band structures and wavefunction coefficients, as well as direct detection rates, proceeds in the same way as in our companion paper [51]. We include the calculation details specific for BN in Appendix A. As a result of the layered crystal structure, the rate is strongly dependent on the angle between the DM wind and the layers. We note, however, that BN has a three-dimensional crystal structure with the layers of BN repeating in the out-of-plane direction, in contrast to single-layer graphene previously considered in Ref. [46].
To show this directional dependence, we consider the same experimental setup as in Refs. [37,44], where the crystal c-axis is aligned with the Earth's velocity v e at time t = 0. With this choice, daily modulation signal is independent of the location of the laboratory. In Fig. 3, we pick three DM masses m χ = 5, 10, 100 MeV to show how the expected detection rates -both total (left panel) and differential (right panel) -change during a sidereal day, assuming a light mediator and negligible in-medium effects. For all three masses, we see that the rate is maximized at t = 12 hours when the DM wind is roughly aligned with the crystal a-b plane, and minimized at t = 0 when the DM wind is aligned with the crystal c-axis. This can be understood from the fact that electron wavefunctions are more localized in the c direction and thus have smaller low-momentum components, whereas the DM scattering matrix element peaks at low q for a light mediator. We also observe that modulation is stronger for lighter DM. Generically, with a smaller energy deposition, the rate is more strongly affected by band structure anisotropies near the band gap; far from the band gap, the electron band structures and wavefunctions approach those for individual, isotropic ions. For DM heavier than 100 MeV, we find roughly the same amount of daily modulation as the m χ = 100 MeV case. This is again because the momentum integral is dominated by small q, which corresponds to the same kinematic region ω q q · v in the large m χ limit. On the other hand, once we go below m χ = 5 MeV, the total rate quickly approaches zero, as the DM does not carry sufficient kinetic energy to trigger a transition across the band gap, which is ∼ 6 eV in BN.

V. SINGLE PHONON EXCITATIONS
Finally, we derive single phonon production rates following the same procedure. Assuming zero temperature, the initial state is the ground state with no phonons, and the final state contains one phonon: where the canonical commutation relations read â ν,k ,â † ν ,k = δ νν δ kk , etc. Note that phonons are labeled by a branch index ν = 1, . . . , 3n, where n is the number of atoms/ions in each primitive To see how F T (q) should be quantized in the phonon Hilbert space, we note that phonons arise from atom/ion displacements: where x lj is the position of the jth atom/ion in the lth primitive cell, x 0 lj is the equilibrium position, m j are the atom/ion masses, ω ν,k are the phonon energies, and ν,k,j are the phonon polarization vectors, normalized such that j | ν,k,j | 2 = 1. The task is thus to find how F T (q) depends on the atom/ion positions x lj and displacements u lj .
To do so, let us revisit the scattering potential in Eq. (4). For a periodic crystal, it can be written as a sum over contributions from individual atoms/ions: where Ω lj is a volume surrounding the lattice site l, j. Within each site volume, we have changed the integration variable to r = x − x lj , the position relative to the center of the site, and defined n lj p (r) ≡ n p (x lj + r), etc. For protons and neutrons, n lj p,n here coincides with n 0 p,n introduced in Sec. III for the nucleus at site l, j. Also, displacing an atom/ion does not change the nucleon distributions inside of a nucleus. Thus, we can write n lj p,n (r) = n j p,n (r) , where the last expression assumes the effect of electron redistribution following an atom/ion displacement is weak and local. This is usually a good approximation for ionic crystals such as gallium e iq·x lj f p n j p (−q) + f n n j n (−q) + f e n j e (−q) + f e δ n lj e (−q) δu lj · u lj where f j = f p Z j + f n (A j − Z j ), and F N j (q) is the nuclear form factor (introduced in Sec. III) for the nucleus occupying site j in each primitive cell. We therefore obtain with where f 0 ψ = f 0 n (f 0 e ) if the rate is written in terms of σ n (σ e ). Note that ∆ j is independent of l due to lattice translation symmetries. From Eq. (84) we see that F T (q) depends on u lj -which are quantized in terms of phonon modes as in Eq. (79) -via both the phase factor e iq·x lj = e iq·(x 0 lj +u lj ) and the ∆ j (q) · u lj term.
With F T (q) quantized in the phonon Hilbert space, we now move on to calculate the matrix element ν, k|F T (q)|0 . We first apply the Baker-Campbell-Hausdorff (BCH) formula to the phase factor e iq·x lj to move annihilation operators to the right: where we have used the fact that the commutator between creation and annihilation operators is a classical number so the BCH series terminates. In the last equation, is the Debye-Waller factor (in the continuum limit 3 with Ω the volume of the primitive cell). The physical meaning of this factor is that a transition |i → |f can be accompanied by additional phonons' creation out of the vacuum followed by their annihilation, and all these processes are resummed into the exponential. The matrix element thus becomes The l sum can be eliminated via the identity where x 0 lj = x l +x 0 j with x l being the position of the lth primitive cell and x 0 j being the equilibrium position of the jth atom/ion within the primitive cell, and G runs over the reciprocal lattice vectors.
In fact, at most one term in the G sum is picked out for given q and k, since k ∈ 1BZ. We will thus drop the G sum in what follows. On each phonon branch, as we sum over k, only the mode that satisfies q = k + G can give a nonzero contribution to the dynamic structure factor, as a result of lattice momentum conservation.
It is worth emphasizing that the notion of momentum conservation here differs from the one familiar in particle physics, due to the spontaneous breaking of continuous translation symmetries.
While each phonon can be thought of as carrying a momentum k within the 1BZ, it can be excited even when the momentum transfer q is outside the 1BZ via Umklapp scattering, in which case G = 0. For DM heavier than ∼ MeV, the momentum transfer can exceed ∼ keV, the typical size of the 1BZ. In this case, Umklapp processes can contribute significantly if the matrix element has support at high q (which is the case for a heavy mediator). We will see an example of this in Sec. V A. Note that momentum is still conserved at the fundamental level: the extra momentum G leads to a recoil of the entire crystal, which becomes unobservable in the limit N → ∞. On the other hand, the notion of energy conservation is the same, as continuous time translation symmetry remains unbroken. As a result, the energy deposition has to match the phonon energy for a phonon mode to be excited.
With the equations above, we obtain the dynamic structure factor: where We have made it implicit in the last line of Eq. (90) that the k vector is the one inside the first Brillouin zone that satisfies q = k + G.
Finally, integrating over the DM velocity distribution, we obtain the rate per target mass: where m cell = ρ T Ω is the mass contained in a primitive cell. The mediator form factor F med , the Debye-Waller factor W j (q) and the g(q, ω) function are given by Eqs. (12), (87) and (27), respectively. The DM couplings are encoded in the Y j vectors given in Eq. (91), with F 0 j , ∆ j defined in Eq. (85). Meanwhile, the material specific quantities -phonon dispersions ω ν,k and polarization vectors ν,k,j -can be numerically computed using DFT methods detailed in our companion paper [51].
In the following subsections, we discuss the phonon excitation calculation in more detail. It is clear from the discussion above that Y j are the key quantities to compute for any specific DM model. In Sec. V A, we consider the simpler case where DM couples only to nucleons but not electrons, and point out an interesting complementarity with nuclear recoils. We also discuss the relevance of Umklapp processes for DM heavier than an MeV, for both heavy and light mediators.
Including DM-electron couplings introduces complications, but we show in Sec. V B that Y j take a simple form in the low q limit for general couplings f p,n,e . Note that the dark photon mediator benchmark (f 0 p = f 0 e , f 0 n = 0) has been studied in Refs. [43,44] based on the Fröhlich Hamiltonian. Our calculation here reproduces previous results, and helps clarify their range of validity.

A. Dark Matter Coupling Only to Nucleons
Setting f e = 0 and f 0 ψ = f 0 n = f n in Eq. (85), we have In this case, Y j is simply F 0 j q, and the rate Eq. (92) becomes It is interesting to compare to the nuclear recoils case. If there is only one atom in the primitive cell, we have m cell = m j = m N , and The differential rate reads On the other hand, we can rewrite Eq. (60) for nuclear recoils in terms of the g(q, ω) function via q dq η(v min ) → 2 d 3 q (2π) 3 g(q, ω), and multiply the integrand by 1 = q 2 2m N ω : One can clearly see the similarity between Eqs. (96) and (97). However, a key difference between nuclear recoils and phonon excitations is the way in which contributions from different atoms add up in the case of more than one atoms in the primitive cell. Comparing Eq. (94) against Eq. (62), we see that, in contrast to the nuclear recoils case where we add up the rates from inequivalent nuclei, for phonon excitations the sum over j is taken at the amplitude level. It is worth noting, however, that this apparent coherence does not result in a more favorable scaling of the detection rate. In fact, the total rate per target mass scales with neither the number of nuclei in the primitive cell, nor the total number of atoms/ions in the crystal. The former can be seen from the fact that phonon polarization vectors scale as ν,k,j ∼ m j /m cell , which, together with the prefactor, means the denominator of Eq. (94) scales as m 2 cell . The latter is because of the 1/ √ N normalization factor when expanding u lj in terms of phonon creation and annihilation operators (see Eq. (79)). The intuition here is that, despite the collective nature of phonon excitations, we have to project the motion of each atom onto the phonon modes that match the energy-momentum transfer. As a result, coherence between more atoms comes with a price of a smaller overlap with phonon modes.
Another key difference between nuclear recoils and phonon excitations, alluded to in Fig. 1 and Sec. III A, is the kinematic regimes probed. In the phonon case, the Debye-Waller factor e −W j cuts off the momentum integral for q √ m N ω ph , the inverse spatial extent of the nucleus wavefunction.
The ω ph here should be thought of as an average phonon energy over the entire 1BZ, which is of the same order as ω max ph . As discussed in Sec. III A, this high q regime is exactly where the nuclear recoil calculation becomes valid. In addition, nuclear recoils happen at much higher energy depositions ω = q 2 /2m N ω max ph than phonon excitations. A multi-channel search can exploit this complementarity between nuclear recoils and phonon excitations. Let us consider, as a benchmark model, a hadrophilic scalar mediator coupling identically to protons and neutrons (f p = f n , f e = 0). In Fig. 4, we compare the reach of the two channels, using GaAs as an example target material. 6 For a heavy mediator (left panel), we see that with sub-eV energy thresholds, nuclear recoils can probe DM masses above ∼ 100 MeV -this is the mass regime where the single phonon excitation rate suffers from Debye-Waller suppression.
Below ∼ 100 MeV where nuclear recoils lose sensitivity, single phonon excitations can probe a few more orders of magnitude of m χ , depending on the energy threshold. For a light mediator (right panel), on the other hand, single phonon excitations outperform nuclear recoils for all m χ . This is because the momentum integral is dominated by the lowest q, which only depends on the energy currently unconstrained parameter space. The heavy mediator case is free from stellar constraints for m φ 400 MeV [62], and the neutrino floor is taken from Ref. [68]. Currently, the best experimental nuclear recoil constraints in this region of parameter space are from DarkSide-50 [8] (assuming binomial fluctuations), and XENON1T (combined limits from [21,22]). We also show the constraint from CRESST-II [4], which is stronger than the DarkSide-50 constraint at low masses assuming no fluctuation in energy quenching. A more complete collection of nuclear recoil constraints can be found in Refs. [8,13,22]. For a light mediator with m φ = 1 eV, fifth force experiments provide the dominant constraint on mediator-nucleon couplings [62].
Meanwhile, the mediator-χ coupling is constrained by DM self interactions (SIDM) if χ makes up all the DM [62], or just by perturbativity (Pert.) if χ is a DM subcomponent (in which case the projected reach can be easily rescaled). 6 threshold, q min ω min /v max . The mass scaling of the curves in Fig. 4 can be understood with a close examination of phase space integrals; we reserve a detailed discussion, including how the various features of the curves depend on material properties, for the companion paper [51].
It is also worth noting that while direct production of single phonons has been proposed mainly as a channel to search for sub-MeV DM, we see from Fig. 4 that its sensitivity extends well above MeV DM masses, which is important for covering the parameter space out of reach in nuclear recoils. A DM particle heavier than ∼ MeV carries a momentum larger than the typical size of the 1BZ (or equivalently, the inverse lattice spacing). However, as explained below Eq. (89), a crystal target is able to absorb a momentum transfer beyond the 1BZ while still producing a phonon, provided the energy deposition matches that of the phonon energy. Such Umklapp processes can contribute significantly to the rate. In Fig. 5, we examine the role of Umklapp scattering by Umklapp processes account for the differences between solid and dashed histograms.
comparing the full rate (solid) vs. contributions from q ∈ 1BZ (dashed), for three DM masses.
We show the differential distribution up to 34 meV, the highest phonon energy in GaAs. For m χ = 0.1 MeV, the maximum momentum transfer q max 2m χ v max 0.56 keV is within the 1BZ, so the solid and dashed histograms coincide. Also, only acoustic phonons with energies below c s q max 9 meV (where c s is the speed of sound) and optical phonons are kinematically accessible; contributions from optical phonons are suppressed at low q [69], so the total rate is dominated by the low energy acoustic phonons. For m χ = 1 MeV and 10 MeV, Umklapp processes dominate the rate in the heavy mediator case, since the momentum integral is dominated by large q. In the light mediator case, the matrix element peaks at small q, so the total rate is well approximated by the 1BZ contribution for sufficiently low energy thresholds (e.g. 1 meV). However, Umklapp scattering can still contribute significantly in the highest energy bins, and dominate the rate if the energy threshold is higher (e.g. 30 meV).

B. Dark Matter With Couplings to Electrons
In the presence of electron couplings f e = 0, information about electron distributions is needed for the rate calculation. We focus on ionic crystals in this subsection, for which Eq. (82) is a good approximation, and the rate formula Eq. (92) directly applies. In this case, we need n j e and δ n lj e /δu lj as input. While n j e can be derived from the same electron wavefunctions as those used in electron transition calculations in Sec. IV, δ n lj e /δu lj is challenging to compute numerically for general q and u lj .
However, the calculation simplifies in the limit q r −1 ion , the inverse ionic radii. As in classical electromagnetism, we can make a multipole expansion, where N e,j is the number of electrons associated with site l, j, and P e,j is the electron contribution to the polarization in the volume Ω lj . Consider the response of the total polarization of the volume to a lattice displacement u lj : where Q j = Z j − N e,j is the total charge. This defines the Born effective charge tensor: 7 Thus, transition |0 → |ν, k . As such, it enters the source term rather than the dielectric matrix ε in Maxwell's equations. In the low q limit, electron contributions to ε below the electronic band gap approach a constant ε ∞ , referred to as the high-frequency dielectric constant.
In the special case of a dark photon mediator that kinetically mixes with the SM photon, (104) and (50), and setting ε → ε ∞ , we obtain By Eq. (92), the rate is therefore Note that since Eq. (104) for Y j is derived in the limit q r −1 ion ∼ O(keV), Eq. (106) holds only when the integral is dominated by this region. This is the case for a light dark photon mediator for any DM mass, since the integrand peaks at small q. In this case, Eq. (106) is in agreement with the result obtain in Ref. [44] based on the Fröhlich Hamiltonian. For a heavy mediator, on the other hand, the integrand peaks at q max = 2m χ v max , so Eq. (106) holds only for Beyond the previously studied dark photon mediator case, our first-principle rate derivation here allows us to compute the reach for other DM models with couplings to electrons. As examples, we consider two benchmark models from Ref. m χ 100 MeV and ω min = 1 meV, the reach extends slightly past the astrophysical constraints, but the rest of the parameter space is constrained. Therefore, as in Ref. [62], we consider the case where χ is a 5% subcomponent of DM, in which case SIDM constraints are absent and single phonon excitations can probe currently unconstrained parameter space. The projected reach for both benchmark models is shown in Fig. 6, where a mediator mass of 1 eV is assumed for definiteness.

VI. CONCLUSIONS
Dark matter direct detection has entered an era in which not only the mass coverage is extending beyond the classic WIMP window -especially into the sub-GeV regime -but also multi-channel target response is becoming an important consideration when designing new experiments. In this paper, we detailed a theoretical framework for calculating spin-independent direct detection rates that can be applied across multiple search channels. Starting from generic DM couplings to the proton, neutron and electron, we factored out material and channel dependent target response into the dynamic structure factor, and derived a procedure to compute this factor which involves quantizing number density operators in the appropriate Hilbert space. We focused on O(eV)gap crystal targets where existing and proposed search channels include nuclear recoils, electron transitions and single phonon excitations, each probing a different kinematic regime (see Fig. 1).
Despite the apparently very different physics involved, the calculation proceeds analogously for all three channels.
While part of this paper has been devoted to rederiving known results in this unified framework, we also obtained several new results, which we summarize in the following: • We have clarified the range of validity of the standard nuclear recoils calculation (Sec. III A).
For energy depositions lower than O(100 meV) in a crystal target, the picture of scattering off single nuclei breaks down. Collective motions of all nuclei have to be considered, with phonons being the appropriate degrees of freedom. The situation is analogous in fluids, though the energy cutoff can be lower (e.g. O(meV) for superfluid helium).
• We have extended the electron transition calculation to account for anisotropic target response, and pointed out the resulting daily modulation can be significant (Sec. IV A). As an example, we considered hexagonal boron nitride, a semiconductor with a 6 eV gap and layered crystal structure, and showed that ±(10 -40)% daily modulation can be expected, depending on the DM mass (Fig. 3).
• As a major new result, we have presented a first-principle derivation of single phonon excitation rates for generic SI couplings. The final result is Eq. (92), where dependence on the relative couplings to the proton, neutron, and electron is fully captured by the quantities Y j .
Computing Y j is straightforward for DM coupling only to nucleons (Sec. V A), but nontrivial in the presence of coupling to electrons (Sec. V B). In the latter case, we have shown that Y j are related to the Born effective charges in an ionic crystal for a general light mediator (not necessarily a dark photon) -see Eq. (104). As examples, we computed the reach for DM scattering via a light hadrophobic scalar or U (1) B−L vector mediator (Fig. 6), where single phonon excitations offer a complementary search channel with competitive sensitivities to previous proposals [62].
• We have pointed out that sensitivity of the single phonon excitation channel is not restricted to sub-MeV DM. For heavier DM, Umklapp contribution can be significant (Fig. 5), and single phonon excitations and nuclear recoils play complementary roles in probing the DM parameter space (Fig. 4).
In addition to shedding light on the connection and complementarity between various existing and proposed direct detection channels, the theoretical framework presented here also makes clear that there is a common algorithm one can follow to study yet unexplored novel detection channels in the future. Some of them will require extending our present formalism beyond SI interactions, a task we plan to take on in future work.