Determining Dark Matter-Electron Scattering Rates from the Dielectric Function

We show that the rate for dark matter-electron scattering in an arbitrary material is determined by an experimentally measurable quantity, the complex dielectric function, for any dark matter interaction that couples to electron density. This formulation automatically includes many-body effects, eliminates all systematic theoretical uncertainties on the electronic wavefunctions, and allows a direct calibration of the spectrum by electromagnetic probes such as infrared spectroscopy, X-ray scattering, and electron energy-loss spectroscopy (EELS). Our formalism applies for several common benchmark models, including spin-independent interactions through scalar and vector mediators of arbitrary mass. We discuss the consequences for standard semiconductor and superconductor targets, and find that the true reach of superconductor detectors for light mediators exceeds previous estimates by several orders of magnitude, with further enhancements possible due to the low-energy tail of the plasmon. Using a heavy-fermion superconductor as an example, we show how our formulation allows a rapid and systematic investigation of novel electron scattering targets.

Dark matter (DM)-electron scattering was first proposed for sub-GeV DM detection less than a decade ago [1], and there has been enormous theoretical  and experimental [38][39][40][41][42][43][44][45][46][47][48][49][50][51] progress since then. Since electrons are not free particles, but are bound in atoms or delocalized across solids, they have favorable kinematics for light DM scattering. However, the rich complexity of condensed matter systems complicates the calculation of scattering rates. Not only do bound electrons have different wavefunctions than their free-particle counterparts [52], many condensed matter systems exhibit collective electronic modes such as plasmons [53]. A formalism describing DM scattering with a single electronic state [3,25] can potentially miss important electron interaction and correlation effects, and must carefully account for 'screening' where the electron density rearranges itself to partially cancel out DM-induced perturbations [6].
In this Letter we propose to bypass the single-particle formulation entirely, and frame the problem of DMelectron scattering in terms of matrix elements of the many-body electron density operator. This perspective is inspired by a classic paper on collective energy loss in solids [54], and since it does not rely on a particular choice of eigenstates, it is equally applicable to all systems: atoms, molecules, metals, insulators, or more exotic materials. Moreover, it intrinsically accounts for all electron interactions and correlations in the target by relating the scattering rate to an experimentally-measurable quantity, the complex dielectric function (q, ω). Crucially, since (q, ω) is defined as a linear response function, the response of the target to a momentum transfer q and energy deposit ω is determined by density matrix elements which are the same whether measured by DM-electron scattering or by an electromagnetic probe [55,56]. The assumption of linear response applies as long as DM interactions are weaker than electromagnetism.
The key result of this Letter is that the total scattering rate for DM with mass m χ and velocity v χ in an arbitrary target is given by where ω q = q·v χ − q 2 2mχ , q = |q|, e is the electron charge, and V (q) is the non-relativistic DM-electron potential. The full derivation can be found in the Supplemental Material (SM), and follows mainly from the arguments made in Ref. [54]. The target-dependent object which appears in the integrand, is known as the loss function. The only assumptions we have made about the DM interactions in deriving Eq. (1) FIG. 1. Schematic depiction of the relevant kinematics for sub-GeV DM. The shaded purple parabolas represent the kinematically-allowed region of q and ω for the labeled DM masses, as in Ref. [25], for a fixed DM speed vχ = 10 −3 , with upper boundary ω = qvχ independent of mχ. The blue and orange shaded regions represent the support of the plasmon part of the loss function. The tail extends into the DM region for conventional materials such as Al and Si, and for heavyfermion materials such as URu2Si2, the plasmon peak lies in the DM region. The range of support for the free electron gas (FEG) loss function is shown in shaded gray, and can be used to approximate the rate in both superconductors and semiconductors over a limited range of ω. The dot-dashed vertical line indicates the size of the Brillouin zone (q ≈ 2.3 keV) of Si, while the horizontal dashed line indicates the band gap above which electron scattering can produce ionization.
are (i ) that the non-relativistic Hamiltonian coupling DM to electrons takes the formĤ int = i V (r χ −r i ), depending only on the electron position operatorsr i and no other operators such as spin or momentum, and (ii ) that H int can be treated perturbatively. The consequences of Eq. (1) are of immediate importance for DM-electron scattering. Spin-independent Hamiltonians arise in many common benchmark models, including those for scattering through scalar and vector mediators. The presence of W implies that all of these interactions are screened. The importance of screening was first noted for a kinetically-mixed dark photon mediator in a solid-state target [6], and later for a scalar mediator [35]. Our results show that a scalar force which couples equally and oppositely to electrons and protons, whether short-or long-ranged, is screened exactly like a kinetically-mixed dark photon. Furthermore, as long as ion contributions to the loss function are negligible (as in semiconductors well above the gap), forces that couple differently to nucleons and electrons are still screened identically. All such screening effects are invisible in a single-particle picture.
Since W(q, ω) is directly measurable through electro-magnetic scattering, DM-electron scattering experiments can be calibrated experimentally, exactly as was done for DM absorption [57][58][59] using the measured real conductivity σ 1 (ω) = (ω/4π)| (0, ω)| 2 W(0, ω). The advantage of our approach is that the loss function can also be modeled semi-analytically in certain relevant energy and momentum regimes, and such models can be compared directly to data. This enables rapid assessment of candidate experimental targets, and potentially bypasses the need for numerical electron wavefunctions to determine the reach of novel detector materials. As shown in the SM, the loss function contains a sum over all possible final states of the target, and thus Eq. (1) represents the maximum possible scattering rate which could be observed at any experiment sensitive to a particular subset of excitations, for example, electron-hole pairs. In the following sections, we show that in a material with free carriers, the loss function scales as W(q, ω) ∝ q at small ω, which can be interpreted as the familiar screening which partially suppresses the 1/q 4 enhancement characteristic of a light mediator. We then show that if W(0, ω) is nonvanishing, a rate enhancement at small q remains whenever ω is kinematically accessible. This behavior of the loss function can arise in two qualitatively different ways: interband transitions in insulators, and long-range plasmons which are generically present in all materials. As we will show, the low-energy plasmon tail may improve the sensitivity of superconducting detectors to light DM by several orders of magnitude, and materials with Fermi velocities slower than v χ may allow DM to access the bulk of the loss function rather than the tail. We illustrate these kinematic regimes in Fig. 1.
In this Letter, we adopt a generic form for the potential, V (q) = V (q) = gχge q 2 +m 2 φ,V , which is valid for DM coupling through a scalar mediator φ or vector V . We compute scattering rates by integrating Eq. (1) over the DM velocity distribution, for which we take the Standard Halo Model (see SM for details). We frame our results in terms of a reference cross section σ e = (µ 2 eχ /π)|V (q 0 )| 2 where µ eχ is the electron-DM reduced mass and q 0 = αm e 3.7 keV is a reference momentum. We show results for a light mediator m 2 φ,V q 2 , with heavy mediator results given in the SM (Fig. S6).

CONVENTIONAL SUPERCONDUCTORS
Ref. [5] first proposed using superconducting metals such as aluminum (Al) as targets for DM-electron scattering. Ref. [6] soon pointed out that long-range Coulomb forces among electrons would screen DM interactions if mediated by a kinetically-mixed dark photon. This effect was incorporated by multiplying the freeparticle matrix element by 1/| RPA (q, ω)| 2 , where RPA is the dielectric function of a free electron gas (FEG) in For Al, the solid line uses W from Ref. [60], and the top of the shaded region uses the FEG model, both with ω ∈ [1 meV, 1 eV]. The URu2Si2 loss function is taken from Ref. [61] with ω ∈ [1 meV, 74 meV], and the shaded region spans W measured along two crystal axes. Si is treated as a FEG with a 2e − threshold, using the ionization model of Ref. [3]. We also show the reach for a Dirac material with density 10 g/cm 3 , gap 2∆ = 20 meV, Fermi velocity vF = 4 × 10 −4 , background dielectric constant κ = 40, and Dirac band cutoff ωmax = 0.5 eV (red); existing constraints from SENSEI [49], SuperCDMS HVeV [51], DAMIC [47], Xenon10 [14], DarkSide-50 [43], and Xenon1T [48] (shaded gray); and the theory target of a freeze-in model when the mediator is a kinetically-mixed dark photon [1,[62][63][64] (dashed blue). The corresponding plot for a heavy mediator is shown in the SM (Fig. S6).
the random phase approximation (RPA) at zero temperature.
Even within RPA, our formalism identifies two important corrections to the DM interaction rate from Ref. [6]. First, all interactions coupling to electron density are screened, including a light scalar mediator and a nonkinetically-mixed vector mediator. This unifies the reach for all models considered in Ref. [6]. Second, the analytic structure of the loss function imposed by causality implies a particular choice of branch cut in RPA differing from that used in Ref. [6] (see SM for details).
The latter correction improves the projected sensitivity of conventional superconductor detectors to DM scattering through a light mediator by several orders of magnitude at low masses. We can understand this by examining RPA in the kinematic regime q k F , ω qv F relevant for sub-MeV DM scattering near the Fermi surface, where k F is the Fermi momentum and v F is the Fermi velocity, respectively 3.5 keV and 6.8 × 10 −3 in Al.
The result is [65] RPA (q, ω) ≈ where λ TF 3.8 keV is the Thomas-Fermi screening length and ω p 15 eV is the plasma frequency. The imaginary part is typically smaller than the real part, so W(q, ω) scales as ω/q 3 1/q 4 ∼ ωq, a much softer screening than the q 4 implied from 1/| | 2 .
Moving beyond RPA, we use the results of Ref. [60], which fits to data a model containing both a 1-loop 'local field' correction to the electron vertex and a q-dependent plasmon width Γ p /ω p 0.1-0.3. The fit implies that the contribution from the ion polarizability in Al is small, justifying our approximation that only electrons contribute to the loss function. The projected reach for a 1 meV threshold is shown in Fig. 2 for a light mediator, with comparisons to previous results given in Fig. S5 of the SM. The orange band reflects theoretical uncertainty in the proper form of the loss function in the energy range of interest (see SM).
In most materials, the loss function features a plasmon with a Lorentzian lineshape peaked at ω p [53,55] and a low-energy tail (see Fig. 1 and SM). In the parametrization of Ref. [60], W(q = 0, ω) scales linearly with ω for ω ω p , and the plasmon tail dominates over the RPA contribution. Our results suggest that a kg-yr exposure of an Al target with a 1 meV threshold is sufficient to cover the entire freeze-in thermal relic target [1,[62][63][64] above 10 keV. However, this depends on the extrapolation of the plasmon tail to meV energies, and existing measurements only characterize the loss function at ω 100 meV [66]. Thus, additional measurements of W are crucial to accurately determine the sensitivity. There may also be contributions to W from coherent scattering with the Cooper pair condensate for energies ω 2∆, as well as finite-temperature effects. We leave investigation of these effects for future work [67].

SEMICONDUCTORS
In a typical semiconductor like silicon (Si) with a gap E g ∼ eV, an energy deposit ω E g requires a momentum deposit q ≥ E g /v χ ∼ keV for v χ ∼ 10 −3 , independent of the DM mass, as shown in Fig. 1. The size of the first Brillouin zone (BZ) in Si is 2π/a 2 keV, where a is the lattice constant. Thus, for ω 2 eV, DM is probing interatomic distances rather than delocalized electrons, and the electrons may be modeled as a FEG with an effective k F 2π/a set by the total valence electron density. This approximation is an excellent match to both density functional theory (DFT) calculations [68] and data [69] for q 5 keV and ω E g in Si [70]; for sufficiently large q (∼ 15 keV, see SM), the bound elec-tron orbitals give large-momentum tails not captured by the FEG.
Equation (1) and Fig. 1 show that at fixed ω, the rate receives contributions from W(q, ω) over many orders of magnitude in q for m χ 10 MeV, so the FEG approximation is best for a light mediator, where V (q) ∝ q −4 weights the integrand most toward small q. Our formalism thus suggests a generic explanation for the behavior of the DM-electron spectrum in the 5-15 eV range (2-4 electron-hole pairs in Si [3]) from light mediator exchange in any conventional semiconductor. The projected reach in Si under the FEG approximation with a 2e − threshold is shown in Fig. 2 for a light mediator.
The differences among various targets become most apparent when ω E g , where the band structure describing delocalized electrons with q 2π/a becomes important. In addition to band structure effects, there is also an irreducible contribution from the plasmon [71], where the tail extends into the kinematically allowed region for DM. This has important implications for rate predictions in currently-operating semiconductor detectors [49][50][51]. DFT calculations predict a rate which peaks in the 1-or 2-electron bin, corresponding to ω 8.3 eV, for all DM masses for which these energies are kinematically accessible [3]. Currently available measurements of W suggest the true rate in these few-electron bins may be somewhat larger. Near-gap effects are quite difficult to model [70], but in our formalism, they can be accounted for by making more precise measurements at ω E g and q E g /v χ .
On the other hand, for near-gap scattering in a narrowgap semiconductor (E g ∼ 10 meV), we have q min 10 eV 2π/a, so the delocalized electrons in the uppermost valence band dominate the behavior of the scattering rate as q → 0. We may understand the absence of screening in these systems through the Lindhard form of the dielectric function [65], which shows that (q, ω) has a finite limit as q → 0, with the imaginary part proportional to the interband transition matrix element. The lack of mobile charge carriers inhibits the screening present in metals. In the next section, we discuss an example of such a narrow-gap semiconductor: a Dirac material.

NOVEL MATERIALS
Our formalism suggests that optimal materials for sub-GeV DM detection will have a loss function with large support for ω < v χ q (Fig. 1). For an ordinary metal with an electron effective mass m * = m e , the loss function is maximized at large q when ω = qv F , where v F = k F /m * 10v χ . This is outside of the kinematicallyallowed region for DM scattering. For small q, collective modes such as the plasmon will dominate, but the plasmon is damped at momenta q > q c ω p /v F [65] due to decay into the particle-hole continuum. Therefore, DM can only excite the undamped plasmon if v χ > v F [32]. Here we explore two qualitatively different ways to achieve v F < v χ : Dirac materials, in which v F is not tied directly to free-electron properties, and heavy-fermion materials, where strongly-correlated electrons can create a Fermi surface with a large m * . Dirac materials, characterized by linear electronic dispersion ω(k) = v F k with widely-varying v F across materials [72], are promising targets for DM detection [16,20,27,28]. Consider a gapless isotropic Dirac material with a single Dirac cone and effective background dielectric constant κ ≡ Re[ (0, 0)]. In typical materials, Re( ) Im( ) over the relevant q and ω [27], and we may write the loss function as The loss function with a gap 2∆ is given in the SM; W Dirac (q, ω) is constant as q → 0 for all ω > 2∆, as anticipated. The loss function immediately displays two key features of scattering in Dirac materials [16]: small v F increases the rate, and scattering is forbidden if v χ < v F for ω = ω q . In Fig. 2, we show the sensitivity of an isotropic Dirac material for a light mediator. This analysis neglects many-body effects, including the plasmon contribution to the loss function. Dirac materials are expected to exhibit two tuneable plasmon modes distinct from the ordinary valence plasmon: a temperature-dependent mode which could lie in the O(meV) range [73][74][75], and a zero-temperature mode tuneable with chemical potential [76]. Therefore, measurements of the loss function in real materials are crucial to accurately estimate the scattering rate, since the plasmon contribution may dominate [77] as was the case for superconductors.
Another way to lower v F is to find materials with ordinary quadratic dispersion but large effective masses. As an example, a number of materials containing f -electrons are known as heavy-fermion systems because they display a Fermi surface with m * ∼ (10-100)m e [78][79][80]. These materials are expected to have a plasmon at energy ω * p T * , the Fermi temperature of the heavy electrons [81]. One such material is URu 2 Si 2 , a heavy-fermion superconductor with T * = 75 K = 6.5 meV and m * 6m e [82], from which one may estimate v F 6.5 × 10 −5 , ω * p T * = 6.5 meV, and q c ω * p /v F 100 eV. In reality, the measured loss function in URu 2 Si 2 [61] shows considerable anisotropy with Lorentzian peaks at either 4 meV or 6 meV depending on the direction of q, as well as a broad peak around 18 meV, which can also be interpreted as a heavy-fermion plasmon (see SM). Despite the extremely rich electron dynamics in this material, in our formalism we may compute the DM rate unambiguously once W is measured in the relevant kinematic regime.
The measured data (see SM, Fig. S1) show that W(ω) ∝ ω above the heavy-fermion plasmon peaks, consistent with the tail of the ordinary valence electron plasmon. However, in contrast to spectra from conventional superconductors or semiconductors, the measured loss function in URu 2 Si 2 shows rich structure which could be used to separate signals from backgrounds not due to fast-particle scattering. Integrating over ω from a threshold of 1 meV up to ω max = 74 meV, the maximum value where data exists, we obtain the projected reach in Fig. 2. The band spans measurements of W(q, ω) as q → 0 along two different crystal axes. We leave a full analysis of the anisotropic response to future work [67]. As expected, the reach in URu 2 Si 2 can surpass Al in the mass range 5-40 keV, where the DM kinetic energy is comparable to the heavy-fermion plasmon energies. Our reach estimates motivate further study of URu 2 Si 2 and similar materials as targets for light DM scattering.

IMPLICATIONS FOR EXPERIMENTS
The advantage of our formulation of the DM scattering rate is that no theoretical input from e.g. DFT is required to compute the scattering rate; the DM energy loss spectrum from spin-independent electron scattering may be precisely predicted from a measurement with an electromagnetic probe in the appropriate kinematic regime. For MeV-GeV DM, X-ray scattering covers the regime q ∼ keV and ω ∼ eV [56], while for keV-MeV DM, momentum-resolved electron energy loss spectroscopy (EELS) can cover q ∼ eV and ω ∼ meV [55,83]. These techniques are standard in condensed matter physics, and a rich literature on measurements of dielectric and loss functions already exists for a number of systems of interest.
The downside of this formalism is that it does not directly predict how many electron-hole pairs are created in the material per unit deposited energy, or how the energy is down-converted from plasmon excitations to charge and phonons. However, if individual quasiparticle contributions to W(q, ω) can be modeled, this information can be reconstructed. (For related work in the context of superconducting targets, see Ref. [84].) Moreover, the quasiparticle contributions may be determined empirically by correlating scattering events using an electromagnetic probe with the partition of excitations read out by the detector, as has been done for nuclear recoil calibrations at higher energy. We argue that these measurements should be considered the primary calibration mechanisms for DM-electron scattering, analogous to photoabsorption for bosonic DM absorption [57][58][59].
Finally, our work may be applied to unify the electronic and phonon descriptions of DM scattering with other sub-gap loss mechanisms that have not yet been explored, such as dielectric heating in insulators or coherent scattering off the superconducting condensate. Dielectric skin depth in the long-wavelength limit q → 0 is proportional to Re[ (ω)]/{ω Im[ (ω)]}, and thus materials with a small skin depth for THz photons and calorimetric readout should respond efficiently to DM-electron scattering, even for meV-scale energy deposits below the eVscale electronic band gaps. Many materials have THz absorption features, so high-resolution THz or infrared transmission spectra are likely fertile ground for exploring new materials for keV-scale DM scattering.
Note added. This work appeared simultaneously with Ref. [85], which also discusses the loss function as a tool for predicting DM scattering rates. Our main results with respect to the loss function are substantively similar, although Ref. [85] emphasizes comparisons with ab initio methods, whereas the present work emphasizes the utility of the loss function formalism in target selection for future experiments. Mariangela Lisanti, Andrea Mitridate, Lucas Wagner, and Kathryn Zurek for enlightening discussions. YK is indebted to Peter Abbamonte for relentlessly (and correctly!) emphasizing the importance of the loss function and plasmon excitations for dark matter scattering. The idea for this work was conceived via an email exchange during the "New Directions in Light Dark Matter" workshop at Fermilab, supported by the Gordon and Betty Moore Foundation and the American Physical Society, and via a Skype call taken during the workshop "Quantum Information and Systems for Fundamental Physics" at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. This project was supported in part by the Fermi National Accelerator Laboratory, managed and operated by Fermi In this Supplemental Material, we provide a number of derivations and further details to support the results in the main Letter. Section I derives our main result for the DM scattering rate in terms of the loss function. Section II outlines a number of simple analytic models for dielectric functions in various materials, and compares them to measured data for Al (representative of an ordinary superconductor), Si (a typical semiconductor), and URu 2 Si 2 (an example of a heavy-fermion superconductor with meV-scale plasmons). We compare the free-electron gas (FEG) model for Si to the spectrum computed using crystal form factors generated by the publicly-available QEdark code [3], and show good qualitative agreement in the range 5-15 eV. Section III is devoted to a detailed comparison of our results for superconductors with other results in the literature, justifying our claim of a stronger reach by several orders of magnitude compared to previous estimates, and gives the projected reach for heavy mediators.

I. SCATTERING RATE IN TERMS OF THE LOSS FUNCTION
Here we derive Eq. (1) and show how the scattering rate for all spin-independent DM-electron interactions is governed by the loss function. Suppose DM couples to electrons through a low-energy Hamiltonian of the form where the sum runs over all electrons in the target. Fourier transforming the potential, we can write the interaction Hamiltonian aŝ where the momentum-space electron density operator is defined aŝ By Fermi's Golden Rule (equivalently, the Born approximation), we can compute the transition rate Γ(v χ ) from the ground state |0 for a given incoming DM velocity v χ , treating the incoming and outgoing DM as plane waves with energy and momentum (E χ , p χ ) and (E χ , p χ ) respectively. We take the ground state to have zero energy without loss of generality. The transition rate is given by [25] Γ where |f is a final state with energy ω f and the sum runs over all possible final states of the system, and we recall that Note that the only assumption that was made here was thatĤ int is sufficiently weak compared to the unperturbed HamiltonianĤ 0 of the target system; this is the case in Ref. [54] for ordinary electron-electron scattering, so it must be the case for DM-electron scattering where the couplings are much weaker. Note this implies one cannot directly apply our result to regions of parameter space where DM and electrons are strongly coupled, as would be relevant for regions in parameter space where DM may not reach underground detectors due to multiple scattering [10-12, 18, 23].
The insight of Ref. [54] is to relate the density matrix element | f |ρ(q)|0 | 2 to an experimentally measurable quantity, the dielectric function (q, ω). The dielectric function is defined as the linear response of the target to the longitudinal electric field of a test charge. For simplicity and to elucidate the formalism, in this work we consider the case of an isotropic material where the dielectric function is a scalar rather than a tensor, and relegate the treatment of the anisotropic case to upcoming work [67]. Since a test charge will also perturb the electron density of the target, it can be shown that this is equivalent to defining the dielectric function as a density-density correlation function [86]. Of course, the electrons will also couple to ions, and strictly speaking the ion density operator should also appear in the measured loss function. In what follows, we will assume that these contributions are negligible, which is the approximation always made in the condensed matter literature. This assumption makes our formalism independent of the DM coupling to protons or neutrons, and even if ion contributions are significant, the measured loss function will give exactly the correct rate for a dark photon mediator.
Because the dielectric function is defined as the linear response of the system, the same assumptions are implicit in the setup of [54] (with an electromagnetic probe) as are present in the DM scattering setup: the test charge interactions are weak compared to the internal interactionsĤ 0 . Therefore the Coulomb potential of the test charge may be factored out in Fourier space, separating the (weak) perturbation due to the probe and the (possibly strong) response of the system to such a probe. The result is [54,86] Here we are using Heaviside-Lorenz conventions for the electron charge e as is common in high-energy physics, which differs from the Gaussian unit definition common in condensed matter physics by a factor of √ 4π. (Note also that Eq. (9) of [54] is missing a factor of π.) Plugging Eq. (S.7) into Eq. (S.5), we obtain our main result, Eq. (1). For reference, the sum over final states on the right-hand side of Eq. (S.7) is (up to factors of π) conventionally defined in the condensed matter literature as the dynamic structure factor.
Notice that we have made no assumptions whatsoever about the character of the final state |f . It is an exact eigenstate of the (in general very complicated) many-body condensed matter Hamiltonian, and the only requirement is that it represents some rearrangement of the electrons in the target so that it has a nonzero matrix element with the electron density operator with respect to the ground state. In this sense our treatment is distinct from Ref. [25], which defines a general dynamic structure factor (with slightly different normalization compared to the condensed matter conventions) very similar to the sum in Eq. (S.5), but one which is excitation-specific and requires quantization in terms of single-quasiparticle states in the case of electron scattering. When all many-body states are included, the structure factor defined in Ref. [25] for electron scattering is identical to the loss function defined through the complex dielectric function, is directly measurable without the need to compute single-particle wavefunctions, and automatically includes all in-medium effects. (Our formalism is philosophically similar to Ref. [87], which parameterizes non-relativistic potentials using only general principles such as the Källén-Lehmann spectral representation, without relying on the assumption of the perturbative exchange of a single mediator.) On the other hand, the formalism of Ref. [25] is useful when DM couples differently to electrons, protons, and neutrons than the photon, in an energy regime where density perturbations to both electrons and ions are relevant, as might be the case for sub-gap single-phonon excitations.
Finally, we note that other UV Lagrangians considered in Ref. [36] also generate non-relativistic potentials which couple to the electron density, but are often accompanied by other spin-or momentum-dependent operators which may complicate our arguments, so we focus on the case of spin-independent scattering. In particular, if DM is a Dirac fermion χ which couples to a scalar φ of mass m φ through the scalar current L ⊃ g χ φχχ, or to a vector V µ of mass m V through the vector current L ⊃ g χ V µχ γ µ χ, and if the mediator couples to electrons in an analogous fashion but with coupling g e , the resulting potential is the same in both cases [36]: (S.8) Similar formulas apply when DM is a complex scalar. Note that in contrast with Ref. [36], we leave the DM-electron coupling as its 'bare' value and place all in-medium corrections to this coupling entirely within the loss function.
In the case where g e ∝ e, as would be the case for a kinetically-mixed dark photon mediator or when the DM is millicharged, the factors of 1/e 2 cancel in Eq. (1) because the DM-induced perturbation to the electron density is exactly proportional to an ordinary electromagnetic probe.
For completeness, we give the expression for the energy spectrum from DM-electron scattering, where ρ T is the mass density of the target, η(v min ) is the mean inverse DM speed vmin d 3 v χ f (v χ )/v χ , and v min = ω q + q 2mχ is the minimum DM speed required to produce an excitation with momentum q and energy ω for DM of mass m χ . To compare with the literature, we take f (v χ ) to be the standard halo model with dispersion v 0 = 220 km/s, escape velocity v esc = 550 km/s, and Earth velocity v E = 232 km/s in the galactic frame. Integrating Eq. (S.9) over ω within the dynamic range of a given experiment gives the total scattering rate.

II. MODELS AND MEASUREMENTS OF THE LOSS FUNCTION
In our formalism, the detector response and its influence on the scattering rate are entirely captured by the complex dielectric function (q, ω) via the loss function W of Eq. (S.7) and Eq. (2). In principle, this quantity is directly measurable with electromagnetic probes in a given material. However, most measurements presently available in the literature are made at values of (q, ω) different than those of interest for the detection of light DM (see Fig. 1). Thus, for a first estimate of the scattering rate, we employ analytical approximations to the dielectric function. Important consistency checks can be implemented based on the fact that −1 is defined as a causal correlation function, and thus must have certain analytic properties. In particular, the following two 'sum rule' relations are satisfied exactly by W(q, ω) in the limit of an isotropic system [88]: (S.11) Equation (S.10) is effectively a manifestation of charge conservation, which explains the appearance of the plasma frequency 12) which is proportional to the total electron density n e in the FEG limit, while Eq. (S.11) follows from causality. Causality also implies that W(q, −ω) = −W(q, ω) [88], which has important consequences for the projected reach in superconductors, as we will see below.
A. RPA dielectric function for a homogeneous electron gas An analytic form for the dielectric function of a homogeneous electron gas can be derived from first principles under the random phase approximation (RPA). Here a word about terminology is in order: screening effects arise from Coulomb interactions between electrons, but in RPA these are embodied in the total scalar potential for the system which is solved for self-consistently [88]. Thus RPA captures only a certain subset of electron interactions without including electron-electron interactions directly in the Hamiltonian; in QFT language, it sums the series of ladder diagrams constructed from the 1-loop vacuum polarization to obtain the resummed photon propagator, but does not include higher-loop diagrams involving additional electron lines. This is the sense in which the electrons are treated as 'free' and RPA is sometimes referred to as the dielectric function for the free electron gas (FEG). Below we will consider further improvements to this approximation.
The resulting dielectric function at zero temperature is given by Eq. (5.4.21) of Ref. [65] as Here Log denotes the principal value of the natural logarithm, k F and v F are the Fermi momentum and Fermi velocity respectively, and Γ p is a free parameter controlling the width of the plasmon which can also be interpreted as a quasiparticle lifetime. The plasma frequency can also be written in the form where λ TF is the (inverse) Thomas-Fermi screening length. We expect the zero-temperature RPA result to be an excellent approximation for ω 2∆, where 2∆ is the superconducting gap. As mentioned in the main text, this approximation ignores possible enhancements to the loss function from scattering off of the condensate at energies near or below the gap, which will be considered in future work [67]. In the literature, Eq. (S.13) is known as the Lindhard dielectric function, though Lindhard's formalism may also be applied to semiconductors as well as metals; in what follows, we will use the terms 'Lindhard,' 'RPA,' and 'FEG' interchangeably to refer to Eq. (S. 13).
Observe that the arguments of the logarithms in Eq. (S.13) are in general complex. For some values of q and ω, these arguments lie along the negative real axis in the narrow-width limit Γ p → 0, and the imaginary part of then depends crucially on the choice of branch. The branch choice is fixed by the causality condition W(q, −ω) = −W(q, ω), which is automatic for positive real values of Γ p , but the Γ p → 0 limit is non-trivial. The causal result is given by Eq. (5.4.22b) of Ref. [65] as (S. 16) where Q ± = q 2k F ± ω qv F . The acausal branch prescription was employed in Ref. [6], which as we will see in Section III below, artificially suppresses the scattering rate for low DM masses.
The imaginary part of the Lindhard dielectric function naturally contains the plasmon as a Lorentzian peak at ω = ω p of width Γ p . For the purposes of light DM detection, kinematics favor energy deposits ω ω p . The plasmon has then typically been neglected in the literature in the computation of the scattering rate, i.e., the rate is computed in the limit Γ p → 0. However, for realistic values of Γ p , the tail of the plasmon peak may significantly contribute to or even dominate the loss function at the relevant values of ω.
For DM-electron scattering in semiconductors, if the deposited energy is O(5 eV) or greater, the minimum momentum transfer is q 5 keV independent of the DM mass (see Fig. 1 in the main text). Since k F 2π/a 5 keV for typical interatomic spacings a, this means that the behavior of this part of the spectrum will be determined by the loss function in the region q k F . For these values of q, the DM is probing length scales smaller than the distance between lattice sites, so we might expect that the inhomogeneities due to the lattice become unimportant and the response is similar to a FEG. For q > 2k F the loss function peaks when Q − ≈ 0, corresponding to ω = q 2 v F 2k F = q 2 2me , which is elastic scattering from free electrons at rest. For a given q, the loss function is nonzero over a range ∆ω 2qv F around the peak, reflecting the fact that electrons at the Fermi surface have a nonzero velocity. Note however that the loss function vanishes when |Q − | > 1, which can happen for sufficiently small ω at sufficiently large q. This is an artificial feature of the FEG which is not present in semiconductors, where the valence (and core) electron wavefunctions have a tight-binding character with a momentum-space tail that extends to arbitrarily large values. This regime corresponds to q Z eff /a 0 15 keV where a 0 is the Bohr radius and Z eff ≈ 4 is the effective nuclear charge felt by the valence electrons in Group 14 elements (carbon, silicon, and germanium). The large-q behavior is especially apparent in some materials like germanium, where the 3d shell may become energetically accessible for ω exceeding the binding energy. A corresponding feature is seen in the spectrum in models using tight-binding wavefunctions [4] as well as those using density functional theory (DFT) techniques [3].

B. Plasmon pole approximation and local field corrections
In the limit that the plasmon dominates, the dielectric function may be derived by modeling the atomic response as a damped harmonic oscillator. This is known as the Fröhlich model [89], and the result is Here c denotes the contribution from core electrons, which is assumed to be independent of q and ω, and ω g is an average band gap which can be set to zero for metals. The corresponding loss function features a Breit-Wigner-like peak, with the form This function satisfies the sum rules of Eqs. (S.10) and (S.11) with c = 1 and ω g = 0. Note that this form of the loss function is linear in ω for ω ω p . The low-energy loss function is also subject to effects which are not included in the Lindhard dielectric function. Ref. [60] (hereafter denoted 'GSRF') fits the plasmon in aluminum including a local-field correction and accounting for the polarizability of atomic cores χ core , resulting in a dielectric function of the form where G(q) is known as the exchange parameter and arises in the microscopic theory from 1-loop corrections to the electron-photon vertex [88]. Ref. [60] provides fits to G and Γ p as functions of q. Complex values of G produce damping, which influences the form of the loss function at small values of ω. However, Im G(q) = 0 can lead to unphysical negative values of the loss function at the smallest values of ω, thereby violating the positivity requirements imposed by the sum rules, and moreover G as computed in various microscopic theories tends to be real [88]. Following Ref. [60], we divide our treatment into two cases, one with complex-valued G ('damped') and one with real-valued G ('undamped').

C. Dielectric function for Dirac materials
Dirac materials are characterized by electrons with the approximately linear dispersion characteristic of relativistic Dirac fermions, rather than the usual quadratic dispersion expected at a band minimum. In real materials, there are typically two such bands, one below and one above the Fermi energy, with dispersions E ± (k) = ± v 2 F k 2 + ∆ 2 . Here, ∆ plays the role of the fermion mass and the Fermi velocity v F is the analogue of the speed of light; the gap at the Dirac point with k = 0 is 2∆. The band structure may be anisotropic, with different Fermi velocities along different lattice directions, but for pedagogical purposes we will focus here on isotropic materials; see Refs. [27,28] for a detailed investigation of anisotropic Dirac materials for DM detection.
In the approximation that only two nondegenerate bands contribute to the Dirac electron spectrum, the dielectric function may be computed using Lindhard's formalism in the Bloch wave basis [65]. At zero temperature, with the valence (−) band full and the conduction (+) band empty, this reads where |k; ± represents a Bloch wavefunction with crystal momentum k in the band − or +, V is the crystal volume, the factor of 2 is for spin degeneracy, and the integral is taken over the first Brillouin zone (BZ) in the continuum limit using the unit cell volume to regularize the momentum sum, k → V uc d 3 k/(2π) 3 . There are some complications with this procedure in the case of anisotropic materials [27], but it yields an accurate estimate for the imaginary part in isotropic materials, which is dominated by the smallest gaps and hence the bands other than the Dirac bands may be neglected. However, as noted in Ref. [27], Re[ (0, 0)] acts as a background dielectric constant receiving contributions from the entire BZ and thus cannot be reliably calculated analytically. We may therefore estimate the real part as simply Re( Dirac ) = κ 1 independent of q and ω over the relevant kinematic range. To obtain the imaginary part, we may use the identity Im(lim η→0 1 x−iη ) = πδ(x) and perform the integral using spinor wavefunctions with the matrix element given in Ref. [16]. Note that this is precisely analogous to performing the phase space integral over the valence and conduction bands in the single-particle formalism for determining the scattering rate; the dielectric function allows us to express the results of Ref. [16] in a more convenient and generalizable formalism. Equivalently, we may recognize that with the replacements ∆ → m e and v F → c, the imaginary part is identical to that of the 1-loop vacuum polarization in relativistic quantum electrodynamics (QED), which is proportional to the cross section for γ * → e + e − by the optical theorem [90]. The result is where the coefficient e 2 /(12π) is (up to a factor of π) the familiar 1-loop beta function coefficient of QED. Indeed, the physics of the dielectric function is the same in Dirac materials as it is in the true QED vacuum; the screening of bare charges due to Im( ) at q 2m e is known as the Uehling potential. As long as v F is not too small, Im( ) 1. (Otherwise perturbation theory would break down, as noted in Ref. [16].) Then if κ 1, we may approximate Im(−1/ ) ≈ Im( )/κ 2 and thus Setting ∆ = 0 gives Eq. (4) in the main text. The last factor may be explained as follows. In real materials, the Dirac band structure does not extend throughout the entire BZ, but deviates from linearity at some point. In Ref. [16] this was expressed as a momentum cutoff Λ, which is required to regularize the real part of Dirac . Here, since we are dealing with model functions rather than real materials, we instead impose a cutoff ω max on the depth of the Dirac band, which has typical values of ω max 0.5 eV in e.g. ZrTe 5 [16]. Finally, note that W Dirac violates the causality requirement W Dirac (q, −ω) = −W Dirac (q, ω). This indicates that W Dirac as computed here does not represent the entire loss function, and in particular (as noted in the main text) it is missing plasmon contributions.

D. Measurements of the loss function in various materials
Measurements of the loss function in the vicinity of the plasmon peak are available in the literature for certain materials, so it is already possible to fit the Fröhlich model directly to data and to assess the significance of the plasmon tail at ω ω p . Figure S1 (left) shows such a fit to measurements in Al. While the fit is excellent in the vicinity of the plasmon peak, the behavior at ω ω p should be viewed only as a benchmark: other physical effects contribute at these energies, notably those encapsulated by the Lindhard dielectric function which incorporates electron screening effects. See Fig. S4 and Section III below for further details.
High-precision measurements of the loss function at nonzero q have also been performed for Si using X-ray scattering [69]. The plasmon is clearly visible at small q, but here we focus on the behavior at large q. Figure S1 (center) shows the measured loss function along the [100] crystal direction (solid lines), compared to the RPA loss function Error bars indicate the accuracy with which the data points could be transcribed from Ref. [60]. Center: loss function for Si at large momenta q > 2π/a, measured from X-ray scattering in Ref. [69]. for the homogeneous electron gas taking ω p = 16.67 eV for the measured plasmon frequency [71]. While semiconductors and insulators do not, strictly speaking, have a Fermi velocity at zero temperature where there are no free carriers, we may regard v F as a tuneable parameter which governs the behavior of the loss function at small ω. With v F = 8.6 × 10 −3 , on the same scale as v F for typical metals, the fit is quite good, especially for ω < 25 eV. On the other hand, at q = 10 keV, the RPA loss function vanishes identically for ω < 12 eV, which is likely unphysical given that atomic tight-binding wavefunctions have support in this kinematic range. The purpose of this comparison is not to advocate for using this extremely simplified model-indeed, data should be used to compute DM rates whenever possible-but rather to demonstrate how in the absence of data a simple model may provide an accurate estimate for the light-mediator spectrum for ω ∈ [5 eV, 10 eV], where the rate integral is dominated by q ∈ [5 keV, 10 keV]. Indeed, the success of the RPA model suggests that this part of the spectrum from scattering in any semiconductor or insulator with eV-scale bandgaps is nearly universal, determined only by the valence electron density and an effective Fermi velocity. This model may be seen as an extension, accounting for screening, of earlier simplified models for scattering in semiconductors using atomic orbitals or tight-binding wavefunctions [2,4].
To complete our survey of sample loss functions, we show in Fig. S1 (right) the measured loss function at q = 0 for URu 2 Si 2 along the a and c crystal axes, measured with Fourier transform infrared spectrometry [61]. URu 2 Si 2 has been extensively studied for decades [91] due to its unusual 'hidden order' below 17.5 K, and thus has been synthesized as ultra-pure single crystals. Below T c = 1.5 K it behaves as a conventional superconductor [92]. A number of features are present below 20 meV which may be interpreted as heavy-fermion plasmons, as we discuss in the main text. Based on this interpretation, to perform our rate estimates in the main text, we extrapolate the loss function as independent of q out to q = q c 100 eV. Indeed, this is the standard approximation made in scattering experiments near the plasmon pole [71]. Then, we see from Eq. (S.9) that the spectrum is largely determined by the shape of the zero-momentum loss function W(ω), with the inverse mean speed η only serving to enforce the kinematic condition q > ω/v χ . All of the approximations we have made may easily be dropped once momentum-resolved data on W(q, ω) within the DM regions shown in Fig. 1 is available.
It is also interesting to note that at larger ω, the loss function is linear to an excellent approximation, in the c direction above 20 meV and in the a direction above 50 meV. In Fig. S1 we show a linear fit to both loss functions Error bars indicate the accuracy with which the data points could be transcribed from Ref. [69]. The shaded purple region represents the kinematically-allowed region for vχ = 10 −3 . The measured loss function agrees fairly well with both the loss function computed from the single-particle basis from QEdark [3] and the Lindhard FEG approximation in the range 5-10 eV for ω, but there are large differences at both small ω near the gap, and near the plasmon energy ωp 17 eV for small q.
with zero offset. In the Fröhlich model Eq. (S.18), the plasmon tail gives a loss function W F (q = 0, ω) ≈ ω × (Γ p /ω 2 p ) at small ω. The slope of the linear fit is consistent with Γ p /ω p 0.1 − 0.2 and ω p 15 eV, which would be reasonable parameters for the ordinary valence electron plasmon in a generic metal. This data therefore provides some preliminary indication that the linear tail of the plasmon in ordinary superconductors like Al may extend down to the meV scale. We emphasize again that dedicated measurements are needed to confirm this.

E. Semiconductor spectrum in the free-electron gas approximation
In order to relate the energy loss function to the crystal form factor [3,25], we compare to Eq. (1), which gives the relation between the dynamic structure factor S(q, ω) defined in Ref. [25] and the loss function via S(q, ω) = 2q 2 e 2 W(q, ω).

(S.24)
On the other hand, the dynamic structure factor in a semiconductor, computed in the basis of single-particle states, can be related to the crystal form factors |f ii kk G | 2 via [3,25] S(q, ω) = 2 where the momentum integral is taken over the first BZ, G runs over all reciprocal lattice vectors, and i and i run over all valence and conduction bands, respectively. Thus we can compute an equivalent loss function from QEdark [3] crystal form factors by Using Eq. (S.26), we can compare the measured loss function to the loss function computed in the single-particle basis by QEdark, as well as the Lindhard dielectric function for the FEG with the best-fit v F in Fig. S1. The results are shown in Fig. S2. Note that for a given q, the range of ω which is accessible is ω < qv χ , which only comprises a small piece of the total support of W(q, ω). Regardless, we see that QEdark tends to slightly underpredict the measured loss in the kinematically-allowed region. Furthermore, QEdark accurately reproduces the measured loss in the near-gap region ω ∈ [1 eV, 5 eV] where Lindhard fails to do so, as expected. On the other hand, QEdark fails to capture the plasmon which is seen in the measured loss function because the single-particle band structure states do not account for collective effects.
Overall, though, the nearly-linear shape of the measured loss function in the range ω ∈ [5 eV, 15 eV] is reproduced fairly well by the Lindhard model, and matches that of QEdark. We therefore expect that the spectral shape (though perhaps not the normalization) will be captured in this energy range by the simple Lindhard model for the loss function. Moreover, since the Lindhard model loss function goes to zero at sufficiently large q for small ω, and since the rate recieves contributions from all q > ω/v χ , we expect the Lindhard approximation to be best for a light mediator which weights the rate integrand by |V (q)| 2 ∝ 1/q 4 . The results are shown in Fig. S3. Indeed, the Lindhard FEG model matches the spectrum fairly well for the light mediator, roughly independent of the DM mass as long as the DM kinetic energy is well above the gap. The spectrum for a heavy mediator is a poorer match, especially at large ω where the kinematic mismatch between the FEG and the bound atomic wavefunctions becomes more important. We emphasize once again that these simple arguments are not meant to replace a measurement of W in the relevant kinematic range, which would predict the spectrum unambiguously. However, they do highlight a qualitative understanding of the spectrum in a limited energy range based on simple material properties like the effective v F , which may be useful for identifying other detector materials suitable for DM-electron scattering. Furthermore, the part of the spectrum where the FEG model performs best corresponds to the 2-electron bin in Si, which is of considerable practical importance to experiments: the 1-electron bin is typically dominated by backgrounds such as leakage current and Cherenkov radiation [37], while the rates in the bins with 3 or more electrons drop precipitously, at least based on estimates from the single-particle loss function. Integrating the FEG spectra from a threshold of ω = 4.7 eV, corresponding to a 2e − threshold in the model of Ref. [3], we obtain the reach curve shown in Fig. 2 in the main text.

III. UPDATED REACH PROJECTIONS FOR SUPERCONDUCTORS
In Ref. [5], the scattering rate in a superconductor is first computed treating the electrons as free particles, with screening included afterwards in Ref. [6] via a correction to the matrix element. We now show that the result of Ref. [6] at T = 0 is exactly reproduced by our Eq. (1) when (q, ω) is taken to be the Lindhard dielectric function in the limit of vanishing plasmon width.
In a relativistic formalism for single-particle scattering, the superconductor scattering rate is given by where q ≡ p χ −p χ denotes the 3-momentum transfer, p χ denotes the momentum of the scattered dark matter particle in the final state, and S(q, ω) (not to be confused with the dynamic structure factor defined in Eq. (S.24) above) characterizes the available phase space, to be defined shortly. The presence of | | 2 in the denominator of Eq. (S.27) accounts for screening and was treated in Ref. [6] as an in-medium modification to the dark photon propagator. In the non-relativistic limit, any interaction of the class considered in Eq. (S.8) gives rise to a matrix element of the form where q = |q|. Equation (S.27) is trivially transformed to an integral over q, and the rate becomes Thus, to agree with Eq. (1), it is sufficient to have S(q, ω) = 2q 2 e 2 Im (q, ω). (S.30) Equation (S.30) holds exactly in the low-temperature limit for the form of S used in Refs. [5,6], where the superconductor is treated as a free electron gas. In this case, S is given by Left: loss function for each of several models for Al, for q = 10 eV. The curve labeled 'Acausal' shows the loss function used in [6], which involves an unphysical choice of branch cut in the complex logarithm. The Fröhlich model fit is the same as that shown in Fig. S1, for which measured data are only available within the red band. The Lindhard model is the RPA dielectric function Eq. (S.13) with Γp = 0, and the GSRF models use fit parameters for Eq. (S.19) from Ref. [60], with 'undamped' corresponding to Im G = 0 and 'damped' corresponding to Im G = 0. The damped curve becomes negative at small ω, which is an unphysical consequence of the GSRF model. The curve labeled 'Data' shows the fit to q = 0 measurements provided by Ref. [66]. We use dashes to indicate the continuation of the fit beyond the range of measured data. The gray band shows the reference range of 1 meV-1 eV deposits. Right: recoil spectra corresponding to each of these loss functions, assuming (mχ, m φ,V ) = (10 keV, 1 µeV) and σe = 10 −39 cm 2 .
where f FD is the Fermi-Dirac distribution and the P i denote 4-momenta. We reserve p i for the magnitudes of 3momenta. The integration over p e is readily performed using the 3-momentum delta function. Writing the p e integral in spherical coordinates and performing the trivial integral over the azimuthal angle produces S(q, ω) = 2 p 2 e dp e d(cos θ) where θ denotes the angle between p e and q. The remaining delta function can be used to evaluate the integral over cos θ, but here care must be taken to enforce | cos θ| ≤ 1. With the appropriate Heaviside function, the final integral becomes S(q, ω) = dp e m e p e πq Now the zero-temperature Fermi-Dirac distribution can be inserted and the integral can be performed analytically. The result is where E q ≡ q 2 /2m e and E ± = E q ± qv F . The conditions in Eq. (S.34) are equivalent to those in Eq. (S.16), i.e., the imaginary part of the Lindhard dielectric function in the limit that the plasmon is infinitely long-lived. Equation (S.30) follows by direct comparison. Given this agreement between the single-particle and dielectric-function formalisms, Eq. (1) can reproduce prior calculations of the scattering rate in superconductors; essentially, the final-state phase space integral is pre-computed in Im( ). However, Eq. (1) is more flexible than the traditional calculation in that we are not limited to the narrowplasmon limit of the Lindhard dielectric function. Any model or measurement of the loss function can be inserted To evaluate the event rate in a superconducting detector, we take the velocity of the DM in the galactic frame to have a modified Maxwell-Boltzmann distribution, For our reach projections, we take v 0 = 220 km/s and v esc = 550 km/s, and we take Earth to have a velocity v E = 232 km/s in the galactic frame. This matches the conventions of Ref. [16]. In order to facilitate comparison with other results in the literature, we also show some results with v E = 0 and v esc = 500 km/s, matching the conventions of e.g. Refs. [5,6]. We refer to this as the 'simple halo' scenario. Finally, for illustrative purposes, we show selected results for a hypothetical halo with v 0 = v esc = 10 4 km/s and v E = 0. In this 'fast DM' scenario, the plasmon peak is kinematically accessible, and this is directly visible as a feature in the recoil spectrum. The various models for the loss functions in Al are shown in Fig. S4, together with the corresponding DM recoil spectra for a kg-yr exposure. The undamped GSRF model and the Lindhard model with Γ p = 0 correspond to the boundaries of the shaded region in Fig. 2. We also show the result obtained by choosing the acausal branch in the Lindhard dielectric function. It is clear from Fig. S4 that at low energies, a naive extrapolation of the plasmon tail dominates over the Lindhard loss function with its infinitely long-lived plasmon. Moreover, the energy range of interest for light DM detection is precisely where the effects of damping in the GSRF loss function become important. While the loss functions given here are valuable benchmarks, the true loss function likely falls somewhere between the Lindhard result and the plasmon tail. This is also suggested by fitting the measurements of Ref. [66], which go down to ω = 100 meV and lie somewhat below the plasmon tail. (See the purple line in Fig. S4.) To accurately predict the DM scattering rate, it is both essential and feasible to measure the loss function in the entire relevant regime of 1 meV < ω < 1 eV. Figure S5 shows updated reach curves for an aluminum superconductor target alongside the results of Refs. [5,6]. The reach curves are specified with respect to a reference cross section defined by σ e = 16πµ 2 eχ α e α χ (α EM m e ) 2 + m 2 φ 2 , (S. 36) where µ χe denotes the reduced mass of the electron-DM system, m φ is the mediator mass, and α e,χ = g 2 e,χ /(4π) in terms of the couplings which define the potential in Eq. (S.8). In Fig. S5, 'light mediator' means m φ α EM m e See text for details. The fast halo scenario is unrealistic and is shown for illustrative purposes only: in this case, the plasmon peak is kinematically accessible, and the recoil spectrum exhibits a corresponding kink at ω = ωp. The shaded areas indicate two fiducial experimental configurations, one sensitive to deposits 1 meV-1 eV (orange), and the other sensitive to deposits 1 eV-1 keV (green). Bottom left: projected reach in an Al superconductor for a 1 kg-yr exposure assuming a heavy mediator. Orange curves show the reach for the lowthreshold scenario, and green curves show the reach for the high-threshold scenario. Projections for the standard halo, simple halo, and fast halo scenarios are shown by the solid, dotted, and dot-dashed curves, respectively. Right: projected reach for a 1 kg-yr exposure of a Dirac material, Si, and the two Al superconductor configurations assuming a heavy scalar or vector mediator. The parameters of the Dirac material are taken as in Fig. 2, with gap 2∆ = 20 meV, Fermi velocity vF = 4 × 10 −4 , background dielectric constant κ = 40, and Dirac band cutoff ωmax = 0.5 eV. The projected reach for Si assumes a two-electron ionization threshold. The projected reach of URu2Si2 lies above the top edge of the plot. All curves assume the standard halo model. For the Al target, the shaded regions indicate the range of variation in different models of the loss function. The solid lines are computed using the GSRF loss function without damping, and the top of each shaded band is computed using the Lindhard loss function. An example of the target parameter space for thermal freeze-out through a heavy dark photon mediator [13] is shown in dashed blue.
(defined with respect to the ordinary electromagnetic fine-structure constant α EM 1/137) and 'heavy mediator' means m φ α EM m e . All reach projections are computed in the zero-temperature limit and assume that the detector is sensitive to deposits between 1 meV and 1 eV. We show reach curves for a high-threshold experiment sensitive to deposits 1 eV-1 keV for a heavy mediator in Fig. S6, along with recoil spectra for selected model points. We also illustrate the appearance of a feature in the recoil spectrum at ω p in the fast halo model, where the plasmon peak is kinematically accessible. To facilitate comparison with the literature, we show reach curves corresponding to an event rate of 3 kg −1 yr −1 , which corresponds roughly to a 95% C.L. constraint. Figure S5 in particular underscores the importance of properly treating the material response. For any interaction of the kind we consider in this work, screening is significant at low DM mass or for a light mediator. However, the implementation of screening in Ref. [6] overestimated the size of the effect for a vector mediator: at the lowest DM masses, the causal branch choice in the logarithms of Eq. (S.13) yields a rate as much as seven orders of magnitude greater than that produced by the acausal choice. Furthermore, accounting for the non-zero width of the plasmon peak further enhances the rate by an order of magnitude or more. The lingering uncertainty in analytical predictions of the loss function can be easily resolved by directly measuring the loss function in promising target materials.