Dark Matter Direct Detection in Materials with Spin-Orbit Coupling

Semiconductors with $\mathcal{O}(\text{meV})$ band gaps have been shown to be promising targets to search for sub-MeV mass dark matter (DM). In this paper we focus on a class of materials where such narrow band gaps arise naturally as a consequence of spin-orbit coupling (SOC). Specifically, we are interested in computing DM-electron scattering and absorption rates in these materials using state-of-the-art density functional theory (DFT) techniques. To do this, we extend the DM interaction rate calculation to include SOC effects which necessitates a generalization to spin-dependent wave functions. We apply our new formalism to calculate limits for several DM benchmark models using an example ZrTe$_{5}$ target and show that the inclusion of SOC can substantially alter projected constraints.

In this work we focus on a specific class of semiconductors for which O(meV) band gaps (as opposed to typical O(eV) band gaps in semiconductors and insulators) arise as a consequence of spin-orbit coupling (SOC) effects. Targets with such small band gaps can probe DM masses down to O(keV) via scattering, and O(meV) via absorption, while still suppressing thermal noise. Moreover, some of these SOC materials have tunable band structures, a property which makes them interesting candidates for direct detection experiments [27].
However, SOC effects introduce some intricacies in the DM-electron interaction rate calculations since the Bloch wave functions are no longer eigenstates of the S z operator, and therefore become two-component objects in spin space. This implies that electron spin sums cannot be trivially performed, and new transition form factors must be computed. For example, spin-dependent vector mediated scattering can no longer be related to its spinindependent counterpart, and must be computed from first principles. We extend the framework in Sec. II and implement the new spin-dependent form factors numerically within EXCEED-DM [30,31], which is publicly available on Github .
To showcase the formalism developed in this paper, we apply it to a target with important SOC effects, ZrTe 5 . This material has been extensively studied in the context of DM direct detection [24][25][26] as a leading candidate for Dirac material targets. Dirac materials are characterized by low-energy excitations which behave like free electrons and satisfy the Dirac equation. The properties of the electronic excitations can then be understood by a simple extension of the standard QED results. An additional consequence is that they have weak electromagnetic screening, even with a small energy gap between the valence and conduction bands, making them a desirable target for sub-MeV dark matter coupled to electrons via a dark photon mediator [25]. In this work, however, we will focus only on ZrTe 5 properties which stem from its SOC nature, and we will not exploit any of the ones deriving from its Dirac nature (which is still debated [32][33][34][35][36][37][38][39][40][41][42][43]).
To illustrate the variety of DM models that an SOC target can probe, we consider several different DM models and processes. Specifically, we will study: • Standard spin-independent (SI) and spindependent (SD) scattering via vector mediators. The fundamental interaction Lagrangians for these models take the form φ µ g χ χγ µ γ 5 χ + g e ψγ µ γ 5 ψ (SD) (1) where ψ and χ are the electron and DM fermion fields, respectively, and φ µ is the dark mediator field.
• Scalar, pseudoscalar, and vector DM absorption. In this case the fundamental interaction Lagrangians take the form g e φψiγ 5 ψ (pseudoscalar DM) g e φ µ ψγ µ ψ (vector DM) . ( The paper is organized as follows. In Sec. II we generalize the DM interaction rate formalism to account for spin-dependent wave functions in general (spin-orbit coupled, anisotropic) targets. Then in Sec. III we apply these results to the candidate material ZrTe 5 and compare the results obtained with and without the inclusion of SOC effects. Further details of DFT calculation are presented in App. A 1, and convergence tests for the results shown in Sec. III can be found in App. A 2.

II. DM INTERACTION RATE FORMALISM
In this section we derive the rates for transitions between electronic energy levels induced by DM absorption and scattering. For the targets of interest here, the electronic energy levels can be labelled by a band index i and a momentum k within the first Brillouin zone (1BZ), which we collectively indicate with an index I = {i, k}. The wave functions of the electronic states can be written in the Bloch form as: where the periodic Bloch wave functions u I are twocomponent vectors in the spin basis, and V is the crystal volume.

A. Absorption
In this subsection, we use the non-relativistic (NR) effective filed theory (EFT) developed in Ref. [19], and summarized in Appendix B, to compute DM absorption rates in materials with sizable SOC.
The absorption rate of a state can be derived from the imaginary part of its self-energy. In a medium, care must be taken due to the possible mixing between the DM, φ, and SM states (in our case the SM photon, A). In the presence of such mixing effects, the DM absorption rate is related to the imaginary part of the self-energy of the "mostly DM" eigenstate, Πφφ: where ω m φ is the energy of the DM state, and is the wave function renormalization which we will approximate as unity in the following. The total absorption rate per unit target mass, R, is given by where n is the number of degrees of freedom of the DM particle (n = 3 for vector DM and n = 1 for scalar and pseudoscalar DM) and we average over the incoming DM polarizations. The DM density, ρ φ , is taken to be 0.4 GeV cm −3 , and ρ T is the target density.
To derive Πφφ we need to diagonalize the in-medium self-energy matrix, which in our case contains a mixing between the DM and the SM photon: where the implicit sum over λ, λ (η, η ) runs over the photon (DM) polarizations, and we have introduced the self-energies polarization components defined as: where λ µ are polarization vectors. In general, the polarization vectors which diagonalize this matrix are not the typical longitudinal and transverse polarization vectors, since mixing can occur (i.e., Π L,T AA = 0). However one can always find an appropriate basis to diagonalize the DM and photon self-energies. In this basis Eq. (6) becomes where Π λ AA and Π η φφ are the eigenvalues of Π λλ AA and Π ηη φφ respectively. The off-diagonal terms in Eq. (8) are perturbatively suppressed by a factor of g e with respect to the Π AA terms. Therefore, working at order O(g 2 e ), we find that the in-medium self-energy for the η polarization of the mostly DM eigenstate is given by: Since vector DM couples to electrons in the same way as the photon, one can derive the relevant self-energies by simply replacing the electromagnetic charge with g e , e.g., Π ηλ φA = −(g e /e) Π λ AA δ ηλ . Doing so allows us to write Eq. (9) in terms of the photon self energy as Scalar and pseudoscalar DM only have one degree of freedom, and therefore Eq. (9) takes the form As usual, the self-energies appearing in the previous equations are computed from the sum of 1PI diagrams. Working at one loop, there are two graph topologies that can contribute where O (1,2) is any operator to which the external field, A or φ (dashed lines), couples. For vector external states these operators carry Lorentz indices that are inherited byΠ andΠ . The full expressions for the self-energies involved in the absorption calculation can be found in Appendix B. However, as we discuss in the same appendix, due to the absorption kinematics and the hierarchy between the DM and electron velocities, a few diagrams dominate these self-energies. Specifically, we find that the diagonalization of the photon in-medium self-energy (and therefore the derivation of Π λ AA ) reduces to diagonalizingΠ v i ,v j , where the velocity operator is defined by From this it follows that the long wavelength limit of the dielectric function, ε(0, ω), which will enter explicitly in the scattering rate calculation, can be derived from The long wavelength dielectric function, ε(0, ω), together with details of its numeric calculation, is reported in Appendix A 2. 1 For scalar and pseudoscalar DM, the leading order terms in the self-energy of the mostly DM eigenstate are found to be where we have introduced the operator 1 Strictly speaking, the dielectric function is a mixed index tensor, as evident from the defining equation, [25] for more details). With a slight abuse of notation we define the matrix ε which has compo- ZrTe5 band structure computed using DFT with SOC (solid lines) and without SOC (dashed lines). The inset highlights the low energy (low E) band dispersion, whose details are sampled using a denser k-point grid. The band gap for the SOC band structure is set to the experimental value of 23.5 meV [44], and the No SOC band structure is shifted accordingly, which gives a larger band gap of 81.6 meV.

B. Scattering
In this subsection we proceed to derive the DM scattering rate with spin-dependent electronic wave functions. Generalizing the formulas previously derived in Refs. [2,5,10,15,31,45], we can write the DM scattering rate as where the bar indicates a spin average (sum) over the incoming (outgoing) DM states, Ψ I (k) are the Fourier transform of the electronic wave functions defined in Eq. (3), and For the SI and SD models of interest here, we can write the free electron scattering amplitude as where F med q q0 encodes the momentum dependence induced by the mediator propagator, f e (q)/f 0 e is the screening factor introduced by in-medium effects, and σ e is a reference cross section defined by with q 0 = αm e . The total rate per unit detector mass is then where g(q, ω) is the velocity integral defined as with f χ (v) being the DM velocity distribution in the laboratory rest frame, which we take to be a boosted Maxwell-Boltzmann distribution with parameters v 0 = 230 km s −1 , v esc = 600 km s −1 , and v e = 240 km s −1 .
The crystal form factor F II is defined as where the spin operators for the models considered in this work are given by and σ i are the Pauli matrices. Given these expressions the form factors, F II , for the SI and SD models take the form, 2 where we have defined the DM model independent transition form factors, T II and T II , as

III. DETECTION RATES IN ZrTe5
We will now apply the formalism developed in the previous section to our benchmark SOC target: ZrTe 5 . The band structure of ZrTe 5 , with and without the inclusion of SOC effects, is shown in Fig. 1. The details of the DFT calculation can be found in Appendix A 1. The dominant effect of SOC is to shift the valence and conduction bands closer at the Γ point relative to the No SOC calculation.
While in theory the calculation of DM interaction rates is identical for O(meV) and O(eV) gap semiconductors, in practice one must be careful about sampling the 1BZ. This is because these O(meV) energy differences generally only occur in small volumes within the 1BZ. To account for this we sample the 1BZ with a higher k-point density in regions corresponding to the low energy band 2 The absence of the overall factor of two, relative to the SI rate formula given in Ref. [31], can be understood from the sum over the states. If the wave functions are spin independent then IF → IF ss , where s (s ) indexes the initial (final) electron spin state. These spin sums contribute the extra factor of two, bringing Eq. (22) and the rate formula in Ref. [31] into agreement. structure. For ZrTe 5 this occurs near the Γ point, and we split the phase space in to two separate regions, "low E" and "high E", which we describe now. The low E region consists of the highest two (one) valence bands and lowest two (one) conduction bands, for the calculation with (without) SOC, sampled on a "mini-BZ" grid. This mini-BZ grid is a rescaled uniform Monkhorst-Pack grid [46]; each k is scaled by a factor of 1/5, giving a 125× kpoint sampling in that region. This region will give the dominant contribution to absorption of DM with mass 100 meV, as well as DM scattering via a light mediator. The high E region includes all the bands outside the low E region, sampled with a standard Monkhorst-Pack uniform grid. The DM absorption rates in Sec. III A are a combination of the low and high E regions. The DM scattering rates in Sec. III B will be shown for both regions, and it will be clear when one dominates the other.
We compute the Bloch wave functions, Eq. (3), in both regions within the framework of DFT with Quantum ESPRESSO [47][48][49]; details can be found in Appendix A 1. The DM absorption and scattering rates are computed with an extended version of EXCEED-DM [30,31] which includes the formalism developed in Sec. II, and is publicly available on Github .
For each of the models considered in the following subsections, we will show the projected constraints from three different calculations. The curves labelled "SOC" are computed with the inclusion of SOC effects, the curves labelled "No SOC" do not include any SOC effects, and those labelled "Partial SOC" are a combination of the SOC and No SOC calculations, obtained using the energy levels computed with SOC, and the wave functions without SOC. While the Partial SOC results are not a consistent calculation, they aid in understanding how much of the difference between the SOC and No SOC results is due to the changes in the band structure versus the inclusion of the spin dependent wave functions. Generally we find that the changes in the band structure are more influential than the spin dependence in the wave functions, but the latter can still be important.
Lastly, we note that previous works [24][25][26] have derived excitation rates analytically by exploiting the putative Dirac nature of ZrTe 5 . While a direct comparison to assess the validity of the analytic approximations is dubious since we do not observe a conical band structure (see Appendix A 1), previous estimates from Ref. [25] are shown in Figs. (2,3). In App. C we discuss the validity of these analytic approximations in more general Dirac materials.

A. Absorption
For the models considered, Eq. (2), our results are shown in Fig. 2. For ease of comparison we map the constraints on the g e parameters in Eq. (2) to a more We compare our results with (solid) and without (dotted) SOC for electronic absorption in a ZrTe5 target (red), with the ones for semiconductor silicon (Si, blue) and germanium (Ge, green) targets [19], superconducting aluminum (Al-SC, brown) [19]), phononic absorption in polar materials [50,51] (GaAs in orange and SiO2 in purple), and previous estimates for ZrTe5 (teal) [25]. We also show the projected constraints combining the SOC energy levels with the No SOC wave functions, ("Partial SOC", red, dashed) to explicitly show the effect of the spin dependent wave functions. Constraints are expressed in terms of the commonly adopted parameters shown in Eq. (29). Shaded red bands correspond to different parameterizations of the electron width δ ∈ [10 −1.5 , 10 −0.5 ] ω used in calculating the self-energies (see e.g. Eq. (B17)), with the solid line corresponding to δ = 10 −1 ω. Thin lines indicate results obtained by rescaling the optical data. Also shown are the direct detection limits from XENON10/100 [13], fifth force constraints [52], and stellar cooling constraints from red giants (RG) [53], and white dwarfs (WD) [54]. For the pseudoscalar scenario we also report the couplings corresponding to the QCD axion in KSVZ and DFSZ models, for 0.28 ≤ tan β ≤ 140 [55].
commonly used notation, where M Pl = 1.22 × 10 19 GeV is the Planck mass. For all the benchmark models, the inclusion of SOC effects dominantly impacts the low mass reach where the SOC corrections to the band structure are most relevant. Most notably, the lowest testable DM mass is shifted as a consequence of the different band gaps: 23.5 meV with SOC, and 81.6 meV without SOC. At higher masses the SOC effects are milder and, as expected, the SOC reach approaches the reach without SOC effects. The close agreement between the "Partial SOC" and "SOC" curves indicates that changes to the energy levels are what is mainly driving the difference in the "SOC" and "No SOC" calculations.
For the scalar and vector DM models we find that ZrTe 5 is superior at low DM masses relative to a superconducting aluminum target, another target material with an O(meV) gap (0.6 meV for the Al-SC curves shown here). However for pseudoscalar DM, for m φ eV, Al-SC yields better sensitivity than ZrTe 5 . This can be attributed to the large amount of screening present in the vector DM case (but not in the pseudoscalar DM case) for Al-SC.
Shaded bands correspond to different width parameterizations, δ ∈ [10 −1.5 , 10 −0.5 ] ω. In theory, the absorption rate calculation is independent of the choice of width; however, when sampling the 1BZ discretely this is not the case, and practically the goal is to find results that have a weak dependence on this parameter. The discrepancy in the shaded bands should be viewed as an uncertainty in the calculation. The constraints turn up on the left hand side because of the band gap, and on the right hand side because of the finite number of bands used in the calculation. All bands for which E −E F < 4 eV, where E F is the valence band maximum were included; see Appendix A 2 for more details.

B. Scattering
We now consider DM-electron scattering in ZrTe 5 for the two benchmark models, standard SI and standard SD interactions, shown in Eq. (1). Specifically we consider a light mediator for the SI model and a heavy mediator for the SD model. For the SI model a light mediator was chosen due to its high sensitivity to the lowest energy excitations, as well as for ease of comparison with other proposals which commonly report constraints on this model. The SD model was chosen to highlight the effect of spin dependent wave functions.  [57] (brown), and previous estimates for ZrTe5 (teal) [25]. We also show the projected constraints combining the SOC energy levels with the No SOC wave functions, ("Partial SOC", red, dashed) to explicitly show the effect of the spin dependent wave functions. Stellar constraints (gray) are taken from Ref. [58] and the freeze-in benchmark (orange) is taken from Ref. [59]. Right: SD model with a heavy mediator (F med = 1). Curves labelled "low/high E" include transitions restricted to the low/high E regions discussed in Sec. III.
The results are shown in Fig. 3 and we discuss them in detail here. Constraints computed in this work are shown in red, with shaded bands corresponding to the uncertainty in the calculation of the screening factor/dielectric function from the electron width parameter, discussed previously in Sec. III A.
When considering the SI model with a light mediator we include anisotropic screening effects in the f e (q)/f 0 e = (q · ε(q, ω) ·q) −1 factor. ε(q, ω) is the dielectric tensor, and this screening factor is especially important for the sub-MeV DM masses considered here. Since in this model the scattering rate is dominated by events with small q, we approximate ε(q, ω) ≈ ε(0, ω), such that we replace the dielectric with the anisotropic, long wavelength dielectric function shown in Fig. 6.
For the SI model, we find that the contribution from transitions in the low E region, discussed earlier in Sec. III, dominate the scattering rate. Therefore, in the left panel of Fig. 3, we only show the results derived from transitions within the low E region. For the massive mediator SD model, in the right panel of Fig. 3, we see that the low E contributions dominate at small DM masses. However for m χ 100 keV, when the high E contributions at O(100 meV) become kinematically available, the high E contributions are dominant. This is due to the fact that when scattering via a heavy mediator the rate is no longer dominated by the smallest momentum transfers. While we did not explicitly include transitions between the low and high E regions, we note that these are only expected to be important for masses where the reach is comparable between the regions, and will not affect the conclusions.
We find that, for the SI model with a light mediator, the inclusion of SOC effects significantly alters the reach for the whole DM mass range considered since the rate is dominated by small energy/momentum depositions. For the SD model with a massive mediator the SOC effects are most prominent for low DM masses when the scattering is probing the band structure near the band gap, which is the most affected by SOC effects. We also see that at the lowest masses the "Partial SOC" curve is closer to the "SOC" than the "No SOC" lines. This shows that while the change to the energy levels is the dominant effect when including SOC, the spin dependence of the wave functions can give O(1) variations.
The left hand side of all the constraint curves are determined by the band gap. The smallest kinematically allowed DM mass is m χ = 6 keV for the SOC calculation with E g = 23.5 meV, and m χ = 21 keV for the No SOC with E g = 81.6 meV. As mentioned in Sec. III A, we only consider bands up to 4 eV above the valence band maximum. Kinematically this means that we are only including all contributions for m χ < MeV, and explains why our projections stop there.

IV. CONCLUSIONS
Materials with strong spin-orbit coupling, such as ZrTe 5 , are promising targets in which electronic excitations can be utilized to search for sub-MeV DM. Their O(meV) band gaps lead to sensitivity to new DM parameter space via both absorption and scattering processes, without relying on detecting single collective excitation modes.
However, due to the spin-orbit coupling, in these materials the electron spin is no longer a good quantum number, and the spin sums over electronic states cannot be trivially reduced. This introduces interesting wrinkles in the DM absorption and scattering rate calculations, which we extended to account for these effects. In addition, we updated the EXCEED-DM program [30,31], which computes DM-electron interaction rates from first principles, to be compatible with this input for future study of general targets with spin-orbit coupling.
We considered a wide range of DM models and processes to which materials with SOC are sensitive: absorption of vector, pseudoscalar, and scalar DM in Sec. III A, and scattering via heavy and light mediators via spinindependent and spin-dependent scattering potentials in Sec. III B. We found that for sub-eV vector and scalar DM absorption, ZrTe 5 is a far superior target relative to an aluminum superconductor. We also found more optimistic projections for SI scattering via a light mediator than previous estimates, and computed, for the first time, the projected constraints on an SD model with a heavy mediator. Our projections for ZrTe 5 lay the foundation for further first-principles studies of materials with strong spin-orbit coupling as targets in direct detection experiments. The DFT calculations are carried out within the generalized gradient approximation (GGA) [60] using the Quantum Espresso code [61] with and without spin-orbit coupling (SOC) included. We use the ZrTe 5 experimental lattice constants, a = 3.9797Å, b = 14.470Å, and c = 13.676Å of the orthorhombic crystal structure [62]. We employ fully relativistic pseudopotentials for calculations including SOC, and scalar relativistic pseudopotentials for calculations without SOC, in both cases generated with Pseudo Dojo [63][64][65]. In each case, we use a 3265 eV kinetic energy cutoff on a uniform 4 × 4 × 2 Brillouin zone (BZ) grid to compute the electron density. To systematically converge the absorption and scattering rates, for the high E region we compute the electronic wave functions with 200, 300, and 400 eV cutoffs on 10 × 10 × 10, 12 × 12 × 12, and 14 × 14 × 14 k-grids. For the low E region, we compute the wave functions with 650, 750, and 850 eV cutoffs on 8 × 8 × 8, 9 × 9 × 9, and 10 × 10 × 10 uniform k-grids in a small reciprocal-space volume that includes the low-energy band dispersion. The convergence of these calculations is discussed in Appendix A 2.
The computed band structure of ZrTe 5 is presented in Fig. 1, where we correct the band gap with a scissor shift to match the experimental band gap for the calculation with SOC. The inset shows in detail the dispersion near the band edges, highlighting the linear dispersion along the intralayer directions Γ-Y and Γ-Z. Note that in interlayer directions (not shown in Fig. 1) the dispersion is not linear or conical. This band structure obtained by combining the experiment lattice constant and the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional is consistent with a previous study [66]. While the presence of a Dirac cone in ZrTe 5 is still under debate [32][33][34][35][36][37][38][39][40][41][42][43], pursuing more extensive tests of crystal structure and DFT functionals, or carrying out beyond-DFT band structure calculations, is beyond the scope of this work.

DM Interaction Constraint Convergence and Dielectric Function
In this appendix we will discuss some details of the DM scattering and absorption rate calculations, as well as the long wavelength, anisotropic dielectric function, ε(0, ω). Since the main focus of this paper is the effect of SOC, only the electronic wave functions near the Fermi surface are needed. This is because, in ZrTe 5 , SOC effects are approximately O(10 meV), and therefore a very small perturbation for states > eV away from the Fermi surface. We are therefore safely within the "valence to conduction" regime, discussed in more detail in Ref. [31], and do not need to study deeper, core electronic levels, or larger energy states where the electrons are close to free. DFT is the preferred tool for studying these transitions, and the two main convergence parameters are the number of k-points in the 1BZ sampling, and the plane wave expansion cutoff, E cut . In both the low E and high E regions we sample k points uniformly with a Monkhorst-Pack grid. The only difference is that the low E points are scaled by 1/5 relative to the high E region. Convergence of the DM absorption and scattering constraints with respect to the k point sampling and E cut parameters are shown in Fig. 4 and Fig. 5 respectively. The constraints in the main text are identical to the most converged constraints shown in Figs. (4,5). Generally we see faster convergence with respect to E cut than the k point density, and slightly faster convergence for the DFT calculation which omits SOC effects than those which include them. We also note that all-electron reconstruction effects were omitted here since we are focusing on very small DM masses, and therefore kinematically limited to small q transitions. However these effects could be important for studies of DM scattering in ZrTe 5 at higher masses, or for larger experimental thresholds.
The dielectric function in the long wavelength limit is shown in Fig. 6, was used as an intermediate to compute a few different constraints. Specifically it was used to screen the SI scattering rate, and it can be shown that the vector DM absorption rate, as well as the pseudoscalar DM absorption rate when wave functions are spin independent, can be related to the dielectric function. Moreover this calculation serves as a useful benchmark to compare future DFT calculations. Convergence of the constraints on DM scattering, for the models discussed in Sec. III B, with respect to the k point sampling (k-grid) and plane wave energy cutoff, Ecut. The first row includes SOC effects while the second row does not. The collection of constraints dominant at the lowest masses corresponds to the low E transitions, and the other set corresponds to the high E transitions. Similar to Fig. 4, the Ecut parameters in the legend correspond to the values used for the low/high E regions.

Appendix B: Generalized Self-Energies
In this appendix we provide the expressions for the self-energies used in the main text, namely Π AA , Π Aφ , and Π φφ . Since electrons in the target are non-relativistic, we will work in the framework of NR EFT closely following Ref. [19], generalizing the results to anisotropic materials with sizable SOC.
At leading order in the NR EFT, it can be shown [19] that the electron-photon coupling reads:  where Σ = diag(σ, σ), and ψ + = 1 2 (1 + γ 0 )ψ NR with ψ NR being the NR electron field defined as For vector DM, by simply replacing eA µ → eA µ − g e φ µ in eq. (B1), we obtain: In deriving the effective interaction Lagrangian for scalar and pseudoscalar DM, we have to keep some NLO terms in the NR expansion. This is because, as discussed in [19], the LO order terms contain factors of the momentum transfer, q, which in the absorption limit induces a larger suppression compared to the electron velocity. Therefore, keeping all the NLO order terms that do not contain factors of q we obtain With these effective interactions, we are now ready to derive the expressions for the self-energies. By using the photon-electron coupling given in eq. (B1), we obtain the expression for Π µν AA in terms of the loop diagrams Π O1O2 and Π O defined in eq. (12) and (13): where ω p = nee 2 me is the plasma frequency, and we have highlighted in red terms that vanish in absence of sizable spin-orbit coupling (in this specific case, they vanish because tr[σ i ] = 0 in absence of SOC). Since vector DM couples e TABLE I. Self-energies scaling with the DM and electron velocities in the absorption limit. Notice that each insertion of the identity operator induces a suppression of order vev φ due to the wave-function orthogonality, and that parity odd self-energies receive an additional suppression of order q/k.
to electrons in the same way of the photon but with a rescaled coupling, κ = g e /e, we have For scalar DM, by using the interactions given in eq. (B4), we get where we have introduced the operatorṽ ij ≡ σ i v j , and as before highlighted in red the terms that vanish in absence of sizable SOC. The terms highlighted in blue, instead, vanish in isotropic materials without SOC. Similarly, by using the couplings given in eq. (B4), we derive the expression for the self-energies of pseudoscalar DM: Due to the absorption kinematics (q ∼ m φ v φ ω ∼ m φ ), and the hierarchy that exists between the DM velocity, v φ ∼ 10 −3 , and the electrons' typical velocity in a crystal, v e ∼ 10 −2 ; only a few terms are actually relevant in the self-energy expressions given above. To facilitate the following discussion, in Table I we summarize the velocity scaling of all the terms appearing in the self-energy expression given above. By using these scaling relations it is easy to see that the photon self-energy (and therefore also the DM self-energy) is dominated by its spatial components, specifically by theΠ v i v j term. For scalar DM, Π φφ is dominated by the termΠv2v2, and the mixing self-energies Π Aφ are suppressed by one power of v φ . Finally, for pseudoscalar DM, Π φφ is dominated byΠṽṽ and the mixing self-energies are again suppressed.
So far we have ignored the tadpole termsΠ O . They can be written in terms of the electronic wave functions as [19]: and are usually related to macroscopic quantities of the material. Specifically,Π vi andΠ σi are related to the current and spin densities of the material (which both vanish for the case of ZrTe 5 ). For the case of ZrTe 5 , the only nonvanishing tadpole term isΠ 1 = n e , which enters in the expression for the vector self-energy. However, we never explicitly compute this term. Instead, we exploit the relation q j me . Indeed, as discussed in [67], a direct numerical derivation ofΠ v i v j − δ ij meΠ 1 would be affected by numerical errors in the ω → 0 limit.
Let's conclude this section by discussing more in detail the scaling relations given in Table I. Indeed, while some of them are trivial, others require some explanation. The expression for the loop diagramsΠ O1,O2 is given by [19]: where V is the total volume, ω I I ≡ E I − E I , δ I I ≡ δ sgn(ω I I ), and f I , f I are the occupation numbers (which, at zero temperature, equal one for states below the Fermi surface, and zero for states above it). From this expression we can see that self-energies involving the identity operators contain the matrix element i , k | e iq·x |i, k which vanishes in the q → 0 limit since |i , k and |i, k are distinct energy eigenstates and therefore orthogonal. At O(q), we have i , k | e iq·x |i, k iq · i , k | x |i, k . One way to compute this matrix element is to trade the position operator for the momentum operator via its commutator with the Hamiltonian. Here we will assume that the Hamiltonian has the form H = p 2 2me + V (x), ignoring the possibility of momentum-dependent or non-local terms in the potential. While these terms can introduce mild corrections (O(10%)) we do not expect them to change the overall scaling of the self-energies so we can ignore them in this context. With this assumption in mind, we can write the matrix element involving the position operator as Writing the wave functions in the Bloch form, we find: where ω i i, k ≡ E i ,k −E i,k . Therefore, in the absorption limit, each identity operator entering in a self-energy diagram induces a suppression of order v e v φ . Parity-odd self energies also vanish in the q → 0 limit. Let's show this explicitly for the case ofΠ v iv2 . By rewriting the electronic wave function in the Bloch form, we can writeΠ v iv2 as By parity invariance the Bloch coefficients satisfy the relation u s ikG = u s i−k−G , therefore, at order q 0 we haveΠ v iv2 = −Π v iv2 = 0. The first non-vanishing contribution arises at order q and is given bȳ Therefore, instead of the naive v 3 e scaling,Πv2 v i scales as (m φ /m e )v 2 e v φ in the absorption limit. By following an analogous derivation, we can conclude that any parity-odd operator receives an additional m φ v φ meve ∼ q k suppression respect to its naive scaling.
where 2∆ is a band gap between the two cones, and with v F the directionally dependent Fermi velocity. The solutions to Eq. (C2) can be found analytically, and most previous works [24][25][26] studying DMelectron interactions in 3D Dirac materials used these analytic solutions as the Bloch wave functions in Eq. (3). Specifically, they used these analytic wave functions to derive scattering and absorption rates.
However, the subtlety is that solutions to Eq. (C2) cannot be the electronic Bloch wave functions since they are not eigenstates of the crystal Hamiltonian, H = p 2 /2m e + V . Therefore while the excitations which satisfy the rescaled Dirac equation, Eq. (C2), are certainly related to the electronic Bloch wavefunctions, they are, generally, not the appropriate wave functions to use when computing DM interaction rates.
To further illustrate this point we will briefly discuss the most well known Dirac material, graphene. Even though it is only two dimensional it will serve as a good example to illustrate the difference between the electronic Bloch wave functions and those which satisfy the Dirac equation. Our discussion here will closely follow Ref. [68], to which we refer the reader for further details.
Graphene has two carbon atoms within a unit cell which form a hexagonal lattice structure. The Bloch wave functions, satisfying the crystal Hamiltonian, are typically found using the "tight-binding" method, which assumes that the Bloch wave functions are a linear combination of the atomic wave wave functions of each of the carbon atoms, where the"A" and "B" indexes refer to the individual carbon atoms (equivalently the individual carbon atom sublattices), ψ j,k are some coefficient functions and X j,k are the specific linear combination of the atomic wave functions which forms a Bloch state, Here r is a lattice vector, r 0 j is the equilibrium position of the carbon atom on the j th sublattice, and N is the number of unit cells in the lattice. The Bloch nature of the X j,k functions can be seen explicitly by noticing that X j,k (x + r) = e ik·r X j,k (x). Assuming that the ψ j,k 's are lattice periodic implies that Ψ i,k is also a valid Bloch state. The idea behind this decomposition is that the ψ j,k 's are slowly varying functions, or envelope functions, in position space, while the atomic wave functions contain the high frequency behavior, being very localized to the atomic sites. Using this intuition we can simplify the full Schödinger equation near the Dirac point This equation can now be "coarse-grained" by integrating out the pieces close to the center of the atoms with the operator, Ω l d 3 x X * l,k for both l ∈ {A, B} sublattices. Assuming that ψ varies slowly over these regions, we can pull ψ j,k out of these integrals and Eq. (C6) becomes two equations, 0 = j=A,B − 1 m e X l,k |∇|X j,k · ∇ − δ l,j E i,k ψ j,k for each l = A, B, where the expectation value of −∇ 2 /2m e + V with respect to X i,k vanishes since we are implicitly assuming k is close to the Dirac point, i.e. at the peak of the conical band.. From symmetry arguments it can be shown that X A,k |∇|X B,k ∝x − iŷ and therefore Eq. (C7) can be further simplified to, which is exactly the rescaled Dirac equation, with v F the Fermi velocity parameter. Therefore we see that the ψ i,k components of the total Bloch wave functions in Eq. (C3) are what satisfies the Dirac equation, not the Ψ i,k which should be used in the excitation rate calculations. Moreover note that the σ operator does not act in spin-space but rather in "sublattice" space, and therefore for spin-dependent excitation rates the spin dependence follows from the X j,k functions.
There are circumstances where the analytic expressions can be used as as approximation. If the tight-binding approximation is valid, and the Bloch wave functions can be cleanly separated in to high and low momentum components (as was just done for graphene), then for q much smaller than typical momentum scale of the X functions the spin independent transition form factors, e.g., Eq. (27) if Ψ is spin-independent, can reduce to the previously used analytic expressions. In these targets the agreement between an analytic and numeric approach is then indicative of how good the tight-binding approximation is. However not all Dirac cones necessarily appear from the same tight-binding approximation as in graphene, and a detailed study of the Bloch wave functions, along with the band structure, should be done to understand whether any analytic approximations will be valid.