Dirac Materials for Sub-MeV Dark Matter Detection: New Targets and Improved Formalism

Because of their tiny band gaps Dirac materials promise to improve the sensitivity for dark matter particles in the sub-MeV mass range by many orders of magnitude. Here we study several candidate materials and calculate the expected rates for dark matter scattering via light and heavy dark photons as well as for dark photon absorption. A particular emphasis is placed on how to distinguish a dark matter signal from background by searching for the characteristic daily modulation of the signal, which arises from the directional sensitivity of anisotropic materials in combination with the rotation of the Earth. We revisit and improve previous calculations and propose two new candidate Dirac materials: BNQ-TTF and Yb$_3$PbO. We perform detailed calculations of the band structures of these materials and of ZrTe$_5$ based on density functional theory and determine the band gap, the Fermi velocities and the dielectric tensor. We show that in both ZrTe$_5$ and BNQ-TTF the amplitude of the daily modulation can be larger than 10% of the total rate, allowing to probe the preferred regions of parameter space even in the presence of sizeable backgrounds. BNQ-TTF is found to be particularly sensitive to small dark matter masses (below 100 keV for scattering and below 50 meV for absorption), while Yb$_3$PbO performs best for heavier particles.


I. INTRODUCTION
The realization that quantum materials, which have been the subject of great attention in recent years, may offer unique opportunities to search for light and very weakly interacting particles has led to a fruitful collaboration between particle physics and condensed matter physics. This development has given new hope to the ongoing search for dark matter (DM) at a time when experimental null results mount increasing pressure on traditional DM models (see, e.g., [1]). Indeed, many novel detection strategies have been developed that promise to probe DM models in regions of parameter space that were previously thought to be experimentally inaccessible [2,3]. This is especially true for DM particles with mass in the keV to MeV range, which would carry so little kinetic energy in the present Universe that their interactions with conventional detectors would be unobservable. While such particles are too light to be produced via the conventional freeze-out mechanism, recent studies have explored many alternative ways to reproduce the observed DM relic abundance, for example via the freezein mechanism [4][5][6][7][8][9][10].
Given the typical velocity of DM particles in the solar neighborhood of v ¼ 10 −3 c, one needs to achieve an energy threshold of less than an eV in order to search for DM particles in the sub-MeV range. Among the proposed materials to achieve this goal are superconductors [11][12][13], superfluids [14][15][16], polar crystals [17][18][19], topological materials [20], and finally Dirac materials [21][22][23], which are the topic of the present work. Dirac materials are defined as materials where the elementary excitations can be effectively described via the Dirac equation [24] with the relativistic flat-metric energy-momentum relation where k denotes the lattice momentum, v F is the Fermi velocity (replacing the speed of light), and 2Δ is the band gap (replacing the rest mass). For jkj ≫ Δ the electrons hence have a linear dispersion relation with coefficient of proportionality given by v F . A crucial advantage of Dirac materials is that the band gap 2Δ, which determines the energy threshold of the material, can be of the order of a few meV. Such small band gaps can arise for example in Dirac semimetals, when a spin degeneracy is lifted by weak spin-orbit coupling or if the underlying symmetry protecting the Dirac node is lifted. A band gap of this magnitude is ideal for the detection of sub-MeV DM particles while at the same time suppressing backgrounds from thermal excitations of electrons. Nevertheless, there are at present no realistic estimates of the expected background level in a Dirac material and existing sensitivity studies in the literature are based on the assumption that backgrounds can be neglected. This might be too optimistic since even in almost perfectly clean samples, states arising in tiny islands of impurity regions can lead to an exponentially small density of states in the mass gap of a Dirac semimetal [25,26]. While this effect is usually negligible, it might play a significant role in rare event searches. As long as one is solely interested in deriving exclusion limits, it may still be justified to ignore backgrounds. However, the question of how the DM nature of a potential signal can be confirmed arises.
In the present work we explore how this question can be answered by searching for a daily modulation in the data. While such a modulation is absent for most backgrounds, it is expected for a DM signal because of the rotation of Earth [18,22,27]. In combination with the motion of the Sun through the Milky Way this rotation leads to a "DM wind" in the laboratory frame that changes its direction over the course of each day. Provided the detector is anisotropic, i.e., that its response depends on the direction of the momentum transfer q, the resulting modulation may allow us to confirm the DM origin of an observed signal.
In Dirac materials such an anisotropy arises from the fact that both the Fermi velocities and the dielectric constants typically differ for the different directions in reciprocal space. It was shown in Ref. [21] that as a result scattering in certain directions may be heavily suppressed or even kinematically forbidden, which makes these materials ideally suited to search for daily modulations. In this work we develop the necessary formalism to calculate the modulation of the DM signal and point out a number of subtleties overlooked in previous studies. We furthermore identify the regions of parameter space of specific models of DM where the modulation is large enough to be detected with statistical significance.
Throughout the paper we will discuss three Dirac materials as potential sensor materials for DM detection. First, ZrTe 5 , which was initially discussed in connection to Dirac materials for DM sensors due to its tiny and well isolated direct gap [21]. Second, we consider the f-electron antiperovskite Yb 3 PbO which was found to exhibit massive Dirac states along certain high-symmetry paths in the Brillouin zone [28]. Third, we follow the outcome of the materials informatics approach to identify potential dark matter sensor materials discussed in Ref. [23] and reveal that one of the three materials mentioned in the study, the quasi-two-dimensional organic molecular crystal bis (naphthoquinone)-tetrathiafulvalene (BNQ-TTF), exhibits various Dirac crossings within the Brillouin zone when spin-orbit coupling is taken into account. These nodes can potentially be gaped by applying stress and as a result breaking some of the crystalline symmetries protecting the Dirac nodes.
In addition to the scattering of sub-MeV DM particles, we also discuss the absorption of bosonic relics with sub-eV masses. We point out that-in contrast to previous claims-the modulation of the signal is absent in this case.
This work is structured as follows. In Sec. II we present the general formalism for the calculation of the expected event rate and its daily modulation, for the cases of both DM scattering and absorption. Section III provides an improved calculation of the polarization tensor for anisotropic Dirac materials. Our numerical calculations of the properties of several candidate Dirac materials are discussed in Sec. IV. In Sec. V we then introduce the statistical method that we employ and present our sensitivity estimates. Finally, we summarize our findings in Sec. VI.

II. DARK MATTER INTERACTIONS IN DIRAC MATERIALS
While Dirac materials can in principle be used to probe many different models of sub-MeV DM, they are particularly well suited for probing Uð1Þ gauge extensions of the Standard Model. These extensions contain a dark photon A 0 which kinetically mixes with the ordinary photon via L ⊃ − ε 2 F μν F 0μν , where F μν (F 0 μν ) denotes the field strength of the (dark) photon. The dark photon can either be a DM candidate itself or it can mediate the interactions between another DM particle and visible matter. The formalism to calculate the resulting detector signals for Dirac materials has been developed in [21,22]. For the case of anisotropic Dirac materials, however, we find a number of pertinent differences with the expressions provided in these works. We will therefore revisit the derivation of the event rates for DM scattering and absorption in detail and provide improved formulas.

A. Scattering rates in Dirac materials
We first consider a DM particle χ with mass m χ which is charged under the new Uð1Þ gauge group. The total DMelectron scattering rate in a Dirac material with volume V is given by where n e stands for the number of valence band electrons per unit mass and V uc for the volume of the unit cell. The factor g ¼ g s g C is the product of spin degeneracy g s and Dirac cone degeneracy g C [29,30]. The rate for lifting one electron with initial and final lattice momentum k and k 0 from the valence band into the conduction band reads [31] with the four-momentum transfer q μ ¼ ðω; qÞ and q ¼ k 0 − k. The DM density is denoted by ρ χ and the reduced mass of the DM-electron system by μ χe . Furthermore, the fiducial DM-electron cross section is defined asσ It is convenient to evaluate the matrix element M 0 for scattering on a free electron at q 2 0 ¼ α 2 m 2 e , where α and m e stand for the fine structure constant and the electron mass respectively. The momentum dependence of the scattering, which results from the propagator of the exchanged dark photon, is then pulled into the form factor [5] In the main part of this work, we will focus on the case of a very light dark photon and, therefore, neglect m A 0 in this expression. The case of a heavy dark photon mediator will be covered in the Appendix. Next, we turn to the form factor F med ðqÞ which accounts for the optical response of the medium. More specifically, it parametrizes the ratio of the in-medium scattering amplitude M over the free amplitude where j 0 and j denote the DM and the electron current respectively. The difference compared to the vacuum case manifests in the appearance of the in-medium photon propagator D instead of the free propagator D 0 . In the last step, we used the fact that the scattering process is nonrelativistic, which implies j 0 ≫ jjj. 1 The in-medium photon propagator can be derived from the Schwinger-Dyson equation for the electromagnetic field [33,34] where Π stands for the photon polarization tensor. We will explicitly calculate Π for Dirac materials in Sec. III. As we will prove there, the spatial components of Π are negligible in the kinematic regime jqj ≫ ω relevant for DM scattering. Therefore, we obtain The scattering rate, furthermore, depends on the transition form factor f k→k 0 , which results from the electron wave functions in the Dirac material [21,31], where 2Δ is the energy gap between the valence band and the conduction band. The tilde indicates that each three-momentum component is rescaled with the Fermi velocity in the corresponding direction, for examplẽ The last ingredient in Eq. (2.2) is the velocity integral, which arises from an integration over the DM velocity distribution fðvÞ: where the factor 2jqj has been introduced for convenience. The total energy of the initial and final state are denoted by E i and E f . In the so-called standard halo model, the DM velocity distribution is given by where N is a normalization factor, v e is Earth's velocity, v 0 and v esc are the velocity dispersion and the galactic escape velocity, and Θ denotes the Heaviside step function. Note that the velocity distribution only depends on v ¼ jvj and cos θ e ¼v ·v e (the hat indicates unit vectors), i.e., fðvÞ ¼ fðv; cos θ e Þ.
In the nonrelativistic limit, the initial and final energy are given by with q again denoting the momentum transfer and We then find ð2:14Þ where we have introduced cos θ q ¼v ·q and the minimal velocity The velocity integral can hence be written as Without loss of generality, we can choose the coordinate system such that the z axis is aligned with q. Furthermore, we require Earth's velocity vector to reside in the y − z plane. In spherical coordinates ðv; θ; ϕÞ one then finds θ q ¼ θ and cos θ e ¼ sin θ sin ϕ sin ψ þ cos θ cos ψ; ð2:17Þ where ψ denotes the angle between q and v e . The integration over cos θ then yields 0 if v < v min and otherwise sets cos θ ¼ v min =v. We therefore find g ¼gðv min ; ψÞ An important feature of this result is that it does not depend on jqj. Indeedg is entirely determined by the two variables v min and ψ. Calculating this integral numerically and tabulating the results as a function of two variables is straightforward. The result is shown in Fig. 1 and confirms the naive expectation that scattering in the direction of the DM wind (i.e., ψ ≈ 0) is strongly preferred.
We can now transform into the laboratory frame, in which v e is time dependent. We adopt the same coordinate system as proposed in Ref. [18], in which v e points in the z direction at t ¼ 0 days and lies in the y − z plane at t ¼ 0.5 days: where α e ¼ 42°is the angle between Earth's rotation axis and its velocity and β ¼ 2π × t=1 days. Finally, ψ is obtained from cos ψ ¼v e ·q. In order to understand the impact of Earth's rotation on the DM scattering rate, it is instructive to consider collisions with ψ ¼ 0 which dominate the velocity integral. For those, we can derive the inequality The fraction of DM particles which can undergo scattering, hence, strongly depends on the Fermi velocity in the direction of the DM wind. This implies strong daily modulations of the scattering rate in anisotropic Dirac materials with v F;y ≠ v F;z . A final subtlety arises from the fact that the analogy between the electron and a free Dirac fermion only applies for sufficiently small momenta k. For larger momenta, the dispersion relation of the electron will deviate from Eq. (1.1). Of course, electrons with such large momenta may still contribute to the event rate, but the formalism outlined above cannot be applied. To obtain a conservative estimate of the event rate, Ref. [21] introduced a cutoff Λ and considered only scattering processes for which k; k 0 < Λ. For a known band structure the cutoff Λ can be determined by identifying the momentum for which the dispersion relation becomes nonlinear.
In the case of an anisotropic Dirac material the definition of Λ becomes more subtle. Indeed, in this case the cutoff momentum typically depends on the direction, Λ ¼ ðΛ x ; Λ y ; Λ z Þ. While in principle it would be possible to apply different cutoffs in different directions, we will again adopt a simpler and more conservative approach and requirẽ ð2:21Þ Note that this prescription differs from the one proposed in Ref. [21], where the maximum is taken rather than the minimum (presumably because of a typographical mistake).

B. Absorption of dark photon dark matter
Let us now consider the case that the dark photon itself constitutes the DM. It can then be absorbed in a Dirac material in analogy to the photoelectric effect. Specifically, we are interested in the absorption of nonrelativistic dark photons with rest mass comparable to the band gap, which implies that ω ≃ m A 0 ≫ jqj. In this regime-as we will show in the next section-the spatial components of the inmedium photon propagator can be approximated as We now want to determine the effective in-medium mixing angle ε med between the dark and the ordinary photon. For anisotropic materials ε med depends on the polarization. Since the dark photons are nonrelativistic, we can conveniently choose the polarization vectors ðϵ x Þ μ ¼ ð0; 1; 0; 0Þ; ðϵ y Þ μ ¼ ð0; 0; 1; 0Þ; The in-medium mixing angle for an x-polarized dark photon is then obtained from the relation 2 ð2:25Þ In the above expression we explicitly indicate that the polarization tensor has to be evaluated at q 2 ¼ m 2 A 0 . The x-polarized dark photon absorption rate is determined from the optical theorem (see, e.g., [36]) Absorption rates for the other two polarizations are obtained in complete analogy. One simply has to replace Π 11 by Π 22 (Π 33 ) for y-polarized (z-polarized) dark photons.
In principle the incoming dark photon polarization needs to be evaluated in the laboratory frame. This complication is, however, usually irrelevant since the dark photons in the solar neighborhood are expected to be unpolarized. Therefore, the rate is simply given by the average The total absorption rate in the detector per unit mass is obtained as The lowest dark photon mass which can be probed by a Dirac material is set by the band gap. Furthermore, the rate has to be cut off when the largest energy deposit consistent with the linear dispersion relation is reached at m A 0 ¼ 2Λ [21]. We emphasize that the absorption of unpolarized dark photons is time independent. This is because the spatial components of the polarization tensor Π ii (with i ¼ 1, 2, 3) are independent of the three-momentum transfer (in the relevant limit ω ≫ jqj). This statement disagrees with Ref. [22], which found a large daily modulation in anisotropic Dirac materials. The discrepancy arises because the photon polarization tensor employed in Ref. [22] carries a residualq 2 dependence which would favor scattering in the direction of the largest Fermi velocity. We will show in the next section that such a momentum dependence is absent and that the absorption rate remains constant with time.

III. POLARIZATION TENSOR IN DIRAC MATERIALS
In this section, we will derive the photon polarization tensor for Dirac materials. The Lagrangian describing photons and electronic excitations in Dirac materials reads For convenience, we introduced the rescaled gamma matricesγ Compared to the Lagrangian of quantum electrodynamics, the speed of light is replaced by the Fermi velocity in the corresponding spatial direction. Furthermore, the role of the electron mass term is played by Δ which is half the band gap. Notice that the structures of the electron kinetic term and the electron-photon vertex coincide as required by gauge invariance. At first order in perturbation theory, the photon polarization tensor is obtained from the diagram shown in Fig. 2. The corresponding amplitude reads where the rescaled four-momentaq μ andk μ are defined analogous toγ μ in Eq. (3.2). The case of a single Dirac fermion in the loop corresponds to a single Dirac cone (g C ¼ 1) with spin degeneracy g s ¼ 2. We keep the factor g in the above expression in order to allow for generic Dirac cone degeneracy. The background dielectric constant κ is taken to be isotropic for the moment (the case of anisotropic κ will be discussed below). In order to employ the well-known result for the vacuum polarization tensor in QED (see, e.g., Ref. [32]), it is convenient to transform the integration measure from k tok. Furthermore, we need to regularize the integral. Choosing the dimensional regularization scheme, we perform the following replacement, ð3:4Þ whereμ denotes the renormalization scale. The resulting polarization tensor can be written in the form where we neglected quartic terms in the v F;i . This is justified since the Fermi velocities are much smaller than the speed of light. For later convenience, we have not included κ in the definition of πðq 2 Þ. As a consistency check, one can easily verify that the polarization tensor fulfills the Ward identities Π μ ν q μ ¼ Π μ ν q ν ¼ 0. This implies that the photon remains massless within the Dirac material. Let us now turn to the polarization function. We find In the modified minimal subtraction scheme (MS), one replaces ð3:7Þ and hence ð3:8Þ In the following, we set the renormalization scale to the cutoffΛ. This choice is motivated by the matching condition for the effective electron charge which is given by e= ffiffi ffi κ p at the cutoff (where electrons should behave as in an insulator). In order to recover the standard expression of Π 00 for an isotropic Dirac material with vanishing band gap (see, e.g., [37]), the precise replacement is The imaginary part of the polarization function, which arises from a negative argument in the logarithm, can be evaluated analytically. One finds FIG. 2. Feynman diagram for the polarization tensor.
For convenience we also state the result for vanishing band gap, We finally want to generalize the photon polarization tensor to the case of an anisotropic background dielectric tensor. Along the principal axes, the latter can be chosen diagonal such that we have Given this form, the spatial components of the polarization tensor can be obtained by the replacement Π ii =κ → Π ii =κ ii in Eq. (3.5) [21]. The remaining components are fixed by the Ward identities. The most general expression for the polarization tensor thus reads Notice that in the kinematic regime relevant for DM scattering jqj ≫ ω, the polarization tensor is strongly dominated by the Π 00 component. With this simplification, the Schwinger-Dyson equation leads to the photon propagator of Eq. (2.7) and therefore This expression improves the corresponding expression in Ref. [21], where the geometric mean of the components of κ is taken instead of including them individually.
In the opposite regime ω ≫ jqj which is relevant for dark photon absorption, the spatial components of Π dominate and one obtains the photon propagator of Eq. (2.22). We make the important observation that for ω ≫ jqj, Π ij becomes independent of the three-momentum q [since the function πðq 2 Þ in Eq. (3.13) only depends on ω 2 in this regime]. As stated earlier, this implies that the dark photon absorption rate in Dirac materials does not depend on the direction of the momentum transfer.

IV. CANDIDATE DIRAC MATERIALS
For our study we consider three potential candidates for Dirac materials based DM sensors: ZrTe 5 , Yb 3 PbO, and BNQ-TTF. In this section we present calculations of their respective band structures and determine the relevant properties. The ab initio calculations were performed in the framework of the density functional theory (DFT) using a pseudopotential projector augmented-wave method [38][39][40][41], as implemented in the Vienna ab initio simulation package (VASP) [42,43]. We compare results for the experimental crystal structure (NR) with results from structurally optimized crystal structures, which where obtained by allowing the unit cell volume to change, but keeping the unit cell shape and the atomic positions unchanged (ISIF7). For the structural optimization and the band structure calculations, we have used the semilocal meta-generalized gradient approximation (GGA) functional (SCAN) [44,45]. To get reliable optimized structural ground states we added van der Waals corrections according to Tkatchenko and Scheffler [46] for the calculations concerning ZrTe 5 and BNQ-TTF.
For the ⃗ k-space integration, we chose a Γ-centered mesh [47] with 14 × 4 × 4 points for ZrTe 5 , 10 × 10 × 10 points for Yb 3 PbO, and 14 × 8 × 2 points for BNQ-TTF. The cutoff energy was set to 600 eV. The calculation of the dielectric tensor was performed using the generalized gradient approximation according to Perdew, Burke, and Ernzerhof [48] and density functional perturbation theory.
The calculations for the band structure and dielectric tensor were performed with spin-orbit coupling, the structural optimization was done without spin-orbit coupling. For Yb 3 PbO the f electrons are considered to be occupied. To push related electronic bands occurring at the Fermi level into the valence band we applied the GGA þ Hubbard U correction using a value of U ¼ 10 eV for the Yb-f orbitals as suggested in Ref. [28].
The unit cells and obtained lattice parameters from the structural optimization in comparison with the reported experimental lattice constants are shown in Fig. 3. We observe that the overall unit cell volume for the computational ground state is slightly decreased for ZrTe 5 and Yb 3 PbO and slightly increased for BNQ-TTF. The increase of the unit cell volume for the structural ground state for organic materials is common and can be traced back to a slightly increased bond length occurring in the DFT calculations.
The obtained band structures for ZrTe 5 , Yb 3 PbO, and BNQ-TTF are shown in Fig. 4.  A slightly strained sample of BNQ-TTF breaking the space group symmetries of the material is therefore likely to introduce a tiny gap. For our sensitivity estimates in the following section, we will consider a band gap of 10 meV (Δ ¼ 5 meV).
We furthermore performed additional band structure calculations for all three materials to fit the occurring Fermi velocities and determine the cutoff radii Λ i where the Dirac dispersion approximately holds. All values are summarized in Table I. The highest Fermi velocities are found for ZrTe 5 , which are of the order of 10 −3 . In contrast, the flatbands of BNQ-TTF lead to very small Fermi velocities in the order of 10 −4 . Due to the low symmetry of ZrTe 5 and BNQ-TTF all three Fermi velocities come with different values. Furthermore, BNQ-TTF is a quasitwo-dimensional material were the dispersion in the k z direction of the Brillouin zone is extremely flat and the corresponding Fermi velocity vanishes. This effect is related to the weak hopping of electrons in the c direction of the crystal stemming from the layered structure of the material. Applying pressure on the sample along the crystallographic c direction will decrease the distance of molecules in the c direction and therefore increase the hoping amplitudes between the molecules. As a result, an increased hopping amplitude will lead to a stronger dispersion of bands opening the opportunity to lift the flatness of the band and tune the Fermi velocity. In the following, we assume that in a sufficiently strained sample the flat direction will take a value of v F;z ¼ 10 −4 for our sensitivity estimates. 3 In comparison to ZrTe 5 and BNQ-TTF, Yb 3 PbO crystallizes in a high-symmetry space group Pm3m. As the Dirac point is observed, e.g., along the path ΓX the little group of ⃗ k is given by C 4v [54,55]. Hence, the rotational symmetry enforces the two Fermi velocities corresponding to the directions orthogonal to ΓX to be degenerated, i.e., v F;y ¼ v F;z ≠ v F;x . For Yb 3 PbO we observe slightly different values for v F;i in the conduction and valence bands. In the conduction (valence) band, the value for v F;x is about v F;x ≈ 2v F;y (v F;x ≈ 1 2 v F;y ). Hence the averaged values for v F;x given in Table I do not reflect this anisotropy. In the following we will use these averaged values to estimate the sensitivity of Yb 3 PbO, but we will not attempt to calculate TABLE I. Calculated Fermi velocities, band gaps, cutoff radii, and Dirac point positions in the Brillouin zone. We compare values obtained using optimized structures (ISIF7) and experimental structures (NR). The values implemented for our sensitivity estimates are highlighted in bold.

Material
Mode  3 We note that for small Fermi velocities the effective coupling strength α eff ¼ α=ðκv F Þ increases and the material becomes increasingly strongly coupled. The considered Fermi velocities of BNQ-TTF imply α eff ∼ 10. It is conceivable that perturbation theory still applies to systems with α eff in this range (see the discussion in [30]). Indeed, this has experimentally been verified for the case of graphene [52]. Nevertheless, we wish to point out that our one-loop calculation of the polarization tensor should only be seen as a qualitative estimate for the case of BNQ-TTF. the modulation signal, which would require an extended formalism allowing for different Fermi velocities in the valence and conduction bands.
We furthermore calculated the values for the dielectric tensor by using density functional perturbation theory as implemented in the code VASP. The values are given in Table II. Due to the tiny gaps present in these materials these calculations are very sensitive to the gap size. However, we observe that for all three materials the diagonal elements κ xx , κ yy , and κ zz dominate over the off-diagonal components. The largest values are found for ZrTe 5 , which is highly anisotropic with κ xx ≈ 308, but κ yy ≈ 21. The smallest values are seen for BNQ-TTF with κ xx ≈ 19 and κ yy ≈ 6. The cubic symmetry in Yb 3 PbO enforces κ xx ¼ κ yy ¼ κ zz ≈ 43.
Finally, we need to determine the optimum orientation of the three Dirac materials in the laboratory. The coordinate system that we introduced above implies that the DM wind points in the z direction at t ¼ 0 days and approximately in the y direction at t ¼ 0.5 days. We hence want to align the materials in such a way that the largest anisotropy is observed in the y-z plane. In the following, we will always align the materials such that the smallest Fermi velocity points in the y direction, while the largest Fermi velocity points in the z direction. 4 Since the event rate is largest when the DM wind is aligned with the smallest Fermi velocity, we expect a daily modulation that peaks at t ¼ 0.5 days.

V. SENSITIVITY ESTIMATES
We are now in a position to calculate the predicted DM signal as a function of time in the Dirac materials that we consider and to estimate their sensitivity. Before presenting our results, we first introduce the statistical approach that we employ.

A. Statistical method for daily modulation
We will consider two possible outcomes for the experiments under consideration. First, we consider the case that the DM hypothesis is incorrect and that the experiments do not observe any DM signal. For example, if no events are observed at all, any parameter point predicting three or more events can be excluded at 95% confidence level. If the experiment observes a number N b of background events, it can still exclude all parameter points for which the probability to observe at most N b signal events is less than 5%. 5 The second outcome we consider is that the experiments do observe a DM signal. In this case it will be essential to confirm the DM nature of the excess by performing a test for daily modulation. Whether or not the daily modulation will be observable depends on both the amplitude of the modulation and the total (i.e., unmodulated) rate. We use the following approach to quantify the significance of the daily modulation.
Each day of observation is divided into the 12 hours around the expected maximum of the modulation and the 12 hours around the minimum of the modulation. Let N max be the total number of events that fall into the former window and N min the remaining events. Assuming N max , N min ≫ 1 the event numbers are expected to follow a normal distribution with estimated standard deviation ffiffiffiffiffiffiffiffiffiffi N max p and ffiffiffiffiffiffiffiffiffi ffi N min p respectively. Hence, the difference N max − N min should follow a normal distribution with standard deviation ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N max þ N min p . In the absence of a daily modulation, the expectation value of this quantity vanishes. To test the hypothesis that there is no modulation, we can hence define the test statistic which we have confirmed to follow a χ 2 distribution with 1 degree of freedom under the null hypothesis using explicit Monte Carlo simulation. If χ 2 s ≫ 1 there is positive evidence for a daily modulation and the hypothesis of no modulation can be rejected. For example, to reject the null hypothesis at 95% confidence level, one would require χ 2 s > 3.84. More generally, the significance of the modulation is given by ffiffiffiffi ffi χ 2 s p standard deviations. PbO, and BNQ-TTF calculated using density functional perturbation theory. In the final column we also specify the respective degeneracy g ¼ g s g C .

Material
Mode κ xx κ yy κ zz κ xy κ xz κ yz κ yx κ zx κ zy g In the following it will be useful to define the total number of signal events N s ¼N max þN min and the modulation fraction A¼ðN max −N min Þ=ðN max þN min Þ. With this definition, the test statistic can simply be written as χ 2 s ¼ A 2 N s . Hence, for a modulation fraction of A ¼ 20% it is necessary to observe N s ≈ 225 events to detect 3σ evidence for a modulation, while for A ¼ 50% fewer than 40 events may be sufficient. Note that given actual data, more sophisticated methods, such as a Lomb-Scargle [56,57] analysis, may reveal even higher significance for a modulation (see, e.g., [58]).
Our approach is easily extended to include a number N b of background events. Assuming that the background does not modulate, it will cancel in the numerator but contribute to the denominator of Eq. (5.1), giving : ð5:2Þ We emphasize that this expression corresponds to the most conservative case without background subtraction and does not require any model of the expected background.
Let us consider the example of an unknown background which has a rate of one event per day. The total exposure is assumed to be 1 kg yr. Based on the total number of observed events alone one can exclude all parameter points that would predict more than ∼400 signal events. Nevertheless, provided the modulation amplitude is sufficiently large, a substantially smaller number of signal events may be sufficient to identify a daily modulation. Indeed, given a modulation fraction of 50% (30%) it would only require about 135 (250) signal events to obtain 3σ evidence for daily modulation.

B. Results for dark matter scattering
We present our main results in Fig. 5 for three different Dirac materials. The two panels in the top row show the expected sensitivity for ZrTe 5 , assuming the Fermi velocities and the band gap obtained from our calculations (left panels) and from experimental measurements (right panels). The two panels in the bottom row correspond to FIG. 5. Expected sensitivity for the Dirac materials ZrTe 5 , BNQ-TTF, and Yb 3 PbO. The two panels in the top row correspond to different assumed properties for ZrTe 5 (see text for details). In the case of a null result, all parameter points above the dashed lines (corresponding to three expected events) can be excluded. In the shaded parameter region it will be possible to identify a daily modulation with 3σ significance in the case in which a DM signal is observed.
BNQ-TTF and Yb 3 PbO, respectively. In each panel the dashed line indicates the exclusion bound from a null result, and the shaded region in the first three panels indicates the parameter space where a daily modulation can be identified with 3σ significance. For the moment we assume that experimental backgrounds are negligible.
For comparison we show in each panel the combination of parameters for which the observed DM relic abundance can be reproduced via the freeze-in mechanism in a model with a massless dark photon. We include the contribution from plasmon decays, recently studied in Refs. [9,10,59]. In the top row we furthermore indicate two benchmark points, corresponding to m χ ¼ 20 keV, σ e ¼ 2 × 10 −41 (orange) and m χ ¼ 50 keV, σ e ¼ 3 × 10 −41 . The predicted event rate for these two points as a function of time is shown in Fig. 6.
We make the surprising observation that the modulation signal is extremely sensitive to the assumed properties of the Dirac material. For the case of ZrTe 5 both the total rate and the modulation amplitude differ substantially for the different values of the Fermi velocities and the band gap. This is investigated more closely in Fig. 7, which in the left panel shows the derivative of the total rate with respect to the cosine of the angle ψ between the momentum transfer q and the velocity of Earth v e . We can see that for t ¼ 0.5 days (i.e., close to the maximum of the rate) the differential event rate looks similar in the two cases and is strongly peaked toward cos ψ ¼ 1, such that the momentum transfer is aligned with the direction of the DM wind.
For t ¼ 0 days on the other hand, there are decisive differences between the two cases. While for the theoretically calculated Fermi velocities and band gap the differential rate still peaks at cos ψ ¼ 1, for the experimental values scattering with cos ψ ≈ 1 is strongly suppressed. This can be traced back to the fact that in this case the Fermi velocity pointing in the direction of the DM wind is  v F ¼ 1.6 × 10 −3 and hence close to the maximum velocity of DM particles in the galactic halo. As a result, only very few DM particles have sufficient kinetic energy to induce scattering with cos ψ ≈ 1 and the event rate is suppressed.
As a result, for the theoretically calculated properties of ZrTe 5 we find a larger total rate but a smaller modulation amplitude than for the experimentally measured properties. This is illustrated in the right panel in Fig. 7, which shows the modulation amplitude (blue lines) and the significance for a daily modulation for σ e ¼ 10 −41 cm 2 (orange lines) in the two cases. We can see that for the experimentally measured properties the modulation amplitude is substantially larger and hence the significance of a daily modulation is increased in spite of the smaller total rate.
Finally, we note that for the theoretical properties of ZrTe 5 the amplitude of the modulation vanishes for m χ ∼ 500 keV and becomes negative for larger DM masses. This is a result of two competing effects. The velocity integral gives the largest contribution if the DM wind points in the direction of the smallest Fermi velocity. At the same time, the combination of form factors F DM and f k→k 0 favors small q but largeq. It, hence, prefers scattering in the direction of the larger Fermi velocities. For small DM masses, the former effect dominates and leads to a modulation peaked at t ¼ 0.5 days while for larger DM masses the second effect can be comparable or even dominant. 6 This can lead to a vanishing modulation amplitude for specific values of the DM mass or even an antimodulation peaked at t ¼ 0 days. Since our definition of χ 2 is symmetric in N max and N min the case of antimodulation is automatically included in our test for daily modulation.
An interesting side remark concerns the dark matter form factor. Since we focused on the exchange of a light dark photon mediator, the latter was taken to scale as F DM ∝ q −2 with the four-momentum transfer. This behavior changes if we consider a heavy mediator exchange for which F DM approaches a constant. We find that the momentum scaling of F DM has profound implications on the modulation of the DM scattering rate. For illustration, we depict the sensitivity of Dirac materials for the heavy mediator case in the Appendix. Most remarkably, the modulation fraction is increased and the flip in the modulation amplitude at m χ ∼ 500 keV completely disappears for ZrTe 5 .
Thanks to its tiny band gap the organic Dirac material BNQ-TTF can probe significantly smaller DM masses than ZrTe 5 . For the assumed value Δ ¼ 5 keV the sensitivity extends down to m χ > 4 keV. Close to the threshold the modulation amplitude is found to be quite large, but it decreases rapidly for heavier DM particles and switches sign for m χ > 100 keV.
For the last material Yb 3 PbO we only show the sensitivity based on the absolute rate. The modulation signal is suppressed due to the very symmetric nature of this material.
Finally, we consider the case where backgrounds are non-negligible and assume for concreteness a background rate of one event per kg day (corresponding to 365 events from background in the assumed exposure of 1 kg yr). The estimated sensitivities in this case are shown in Fig. 8. In this case, only parameter points predicting more than 400 signal events can be excluded based on the absolute rate, and the resulting bounds are therefore much weaker than in Fig. 5. However, the parameter region where a DM signal can be identified based on its daily modulation remains almost unchanged. In fact, for DM masses close to the kinematic threshold, the modulation fraction can be so large that a DM signal can be identified even if the DM signal is significantly smaller than the number of background events.
To conclude this discussion, let us briefly comment on the dependence of our results on the cutoff Λ. It is clear that introducing such a cutoff leads to conservative results for the total rate, since only a part of the Brillouin zone is included in the calculation. However, for the modulation amplitude it could in principle happen that the region excluded from the calculation contributes to the modulation with the opposite phase and would hence reduce rather than increase the modulation amplitude. We have therefore explicitly confirmed that variations in the cutoff Λ do not significantly modify any of the results presented in this section. The reason is that the dominant contribution to DM scattering stems from collisions with small momentum transfer, for which the precise value of the cutoff is irrelevant. Only if the modulation amplitude nearly vanishes (i.e., close to the transition from modulation to antimodulation) can the dependence on Λ be sizable. Since the search for a daily modulation is essentially insensitive in this particular region, this dependence on Λ is of little practical importance.

C. Summary for dark matter scattering and absorption
Our sensitivity studies for the three considered Dirac materials are summarized in Fig. 9. The left panel covers DM scattering, while the right panel refers to dark photon absorption. Intriguingly, the two new materials suggested in this work, BNQ-TTF and Yb 3 PbO, reach a very competitive sensitivity compared to ZrTe 5 for both cases. The smallness of the band gap and Fermi velocities make BNQ-TTF the ideal target to search for scattering (absorption) of DM particles with masses down to a few keV (meV). Yb 3 PbO, on the other hand, has a relatively large band gap and is therefore only sensitive to m χ ≳ 20 keV and m A 0 ≳ 40 meV respectively. While the large amount of symmetry makes this material unsuitable to search for a daily modulation, the small Fermi velocities combined with the large cutoff scale imply the best sensitivity to DM particles with m χ > 100 keV (m A 0 ≳ 50 meV) based on the total rate alone. white dwarfs and red giants [60]. For larger DM masses, on the other hand, the leading constraint comes from SN1987A, which is compatible with σ e ≲ 10 −35 cm 2 [62]. Dirac materials may therefore improve on these constraints by up to 8 orders of magnitude. A similar picture emerges for dark photon absorption, where Dirac material detectors could improve existing limits on the kinetic mixing by 4-6 orders of magnitude in the range m A 0 ¼ 10 meV-1 eV. In order to highlight this impressive sensitivity, let us note that even the extremely tiny kinetic mixing induced by gravity at six-loop order [63] which has ϵ ≲ 10 −13 is within reach for Dirac materials.

VI. CONCLUSIONS
Detectors built from Dirac materials with sub-eV band gap are one of the most promising strategies to search for sub-MeV DM particles interacting with electrons via the exchange of a dark photon. At the same time they can search directly for the absorption of dark photons with sub-eV masses. In the present work we have studied the properties of several different Dirac materials in order to answer the question of how a potential DM signal in such a material can be distinguished from backgrounds. The central observation is that in anisotropic materials the DM signal is predicted to exhibit a daily modulation due to the rotation of Earth relative to the incoming DM wind, which can be used to reject the background hypothesis.
In the first part of this work (Secs. II and III) we have revisited the formalism to calculate experimental event rates for the scattering or absorption of DM particles in anisotropic Dirac materials and have provided a number of improvements to previous results: (a) In Eq. (2.18) we have introduced a simple way to include the anisotropy of the DM velocity distribution by calculating the velocity integral in terms of the minimum velocity v min and the angle ψ between the momentum transfer and the velocity of Earth. (b) In Eq. (2.21) we have proposed an improved way of defining the cutoffΛ that determines the region of reciprocal space where the electrons behave like a free Dirac fermion. (c) We have shown that in the case of (unpolarized) dark photon absorption there is no daily modulation. This conclusion has been drawn from the fact that the absorption rate is determined by the spatial components of the polarization tensor, see Eq. (2.27). The latter has been found to be independent of the threemomentum transfer in the kinematic regime relevant for absorption (ω ≫ jqj), see Eq. (3.13). (d) Equation (3.14) has provided the correct expression for the in-medium form factor F med in the case in which the dielectric tensor is anisotropic. (e) We point out that the daily modulation depends sensitively on the assumed DM form factor and that the modulation amplitude is significantly larger for the case of a heavy mediator than for a light mediator (see the Appendix). In the second part of this work (Sec. IV) we have presented a number of candidate Dirac materials that possess the required properties to detect DM particles in the sub-MeV range. We have performed an improved calculation of the band structure of ZrTe 5 and have determined the band gap, Fermi velocities, and the dielectric tensor (see Tables I and II). In particular, we have confirmed the finding that this material exhibits a sizable anisotropy, which makes it particularly well suited to search for a daily modulation.
We furthermore have proposed two new Dirac materials, BNQ-TTF and Yb 3 PbO, which have not been previously considered in the context of DM physics. Both materials have significantly smaller Fermi velocities and therefore potentially much larger sensitivity to DM scattering than ZrTe 5 . While Yb 3 PbO crystallizes in a cubic lattice and therefore exhibits little anisotropy, BNQ-TTF is found to be highly anisotropic and furthermore exhibits a tiny band gap, making this material extremely attractive for further investigations. As reported in Ref. [51], a macroscopic sample of BNQ-TTF has already been synthesized, feasible for usage in devices. It will be exciting to see whether the properties that we predict can be confirmed in the laboratory.
Finally, in Sec. V we have provided our sensitivity estimates for the three Dirac materials (see Fig. 5). We have identified the parameter regions that can be excluded by a null result as well as the parameter regions where the daily modulation is large enough to provide a way to confirm the DM nature of an observed signal. The statistical method that we use to search for daily modulations can easily be extended to include a nonmodulating background contribution and we have found that anisotropic Dirac materials retain an impressive sensitivity to DM scattering even in the presence of sizable backgrounds (see Fig. 8). However, we have also concluded that the modulation signal depends sensitively on the properties of the Dirac material, in particular the Fermi velocities, making a precise determination of these properties essential.
Clearly, there is still a long way to go before the first DM detector based on a Dirac material will be built. Nevertheless, as the interest for DM models in the sub-MeV mass range grows rapidly, there will be increasing incentive to exploit the great potential of this technology. Both improved calculations and experimental measurements will be essential in order to identify the materials most suited for exploring this uncharted territory of DM physics.  [22] appeared, which also considers the daily modulation of DM signals in Dirac materials. For the case of DM scattering in ZrTe 5 our results are in qualitative agreement, but there are important differences for the case of DM absorption, as discussed in detail in the text.

APPENDIX: SENSITIVITY OF DIRAC MATERIALS FOR DM SCATTERING WITH HEAVY MEDIATORS
In this appendix we consider the case of a heavy dark photon with m A 0 ≫ jqj. In this case the DM-electron scattering cross section becomes independent of the momentum transfer and the DM form factor simply becomes F DM ðqÞ ¼ 1. Compared to the case of a light dark photon, this leads to a strong suppression of the scattering rate for small DM masses (i.e., small momentum transfer) but a much milder suppression for large DM masses (large momentum transfer). This can be seen in Fig. 10, which shows the parameter regions that could be excluded by an experimental null result (corresponding to three expected events) as well as the parameter regions where one could detect a daily modulation at 3σ significance.
We observe that for the case of a heavy mediator the modulation fraction is generally enhanced. This is because the DM form factor for a heavy mediator leads to a much weaker preference for scattering in the direction of large Fermi velocities than the DM form factor for a light mediator. For example, in ZrTe 5 we find the modulation fraction to be greater than 10% for all DM masses and no change of sign for large DM masses. Unfortunately, we find that for large DM masses our sensitivity estimates depend on the adopted value of the cutoff Λ, because scattering with momenta close to the cutoff gives an important contribution. This explains in particular why the sensitivity of ZrTe 5 is much worse when using the properties from Ref. [21] (including Λ ¼ 0.07 Å) than for the properties determined from our calculations (for which Λ is significantly larger). The sensitivities shown in Fig. 10 should therefore be interpreted as conservative estimates. 7 Compared to m A 0 ¼ 0, the heavy mediator case is more strongly constrained by cosmological and astrophysical observations like the effective neutrino number and supernovae (see, e.g., Ref. [64]). However, in the heavy mediator case, the cosmological bounds onσ e are sensitive to the choice of the dark sector gauge coupling and the precise dark photon mass. Since heavy mediators are not our main focus, a more detailed investigation of the corresponding astrophysical and cosmological constraints is beyond the scope of this work.