Dark Matter Absorption via Electronic Excitations

We revisit the calculation of bosonic dark matter absorption via electronic excitations. Working in an effective field theory framework and consistently taking into account in-medium effects, we clarify the relation between dark matter and photon absorption. As is well-known, for vector (dark photon) and pseudoscalar (axion-like particle) dark matter, the absorption rates can be simply related to the target material's optical properties. However, this is not the case for scalar dark matter, where the dominant contribution comes from a different operator than the one contributing to photon absorption, which is formally next-to-leading-order and does not suffer from in-medium screening. It is therefore imperative to have reliable first-principles numerical calculations and/or semi-analytic modeling in order to predict the detection rate. We present updated sensitivity projections for semiconductor crystal and superconductor targets for ongoing and proposed direct detection experiments.

transitions between electronic states.
It has been widely appreciated that, for several well-motivated bosonic DM models, the absorption process is closely related to that of photon absorption, and the rate can be expressed in terms of the target material's optical properties, i.e. the (complex) conductivity or dielectric function. In fact, most studies on DM absorption so far have utilized this feature to make rate predictions by simply rescaling optical data. This approach is obviously attractive because it saves the labor of first-principles calculations, which can be technically challenging or resource-intensive, and because one can often make quick comparisons between target materials based on existing data.
Nevertheless, this data-driven approach has important limitations. First of all, conductivity/dielectric data are not always readily available, especially for newly proposed, more exotic materials, in which case one has to resort to first-principles calculations and/or semi-analytic modeling (this is the case, e.g. for Dirac materials studied in several recent works [25][26][27][28]). Meanwhile, and more importantly, the question of whether DM absorption for a particular model can be simply related to photon absorption is a nontrivial one, and explicit calculations are needed to establish the answer.
It is the purpose of this work to revisit the calculation of DM absorption via electronic excitations.
We critically examine the question above by carefully working out the matching between relativistic Lagrangians for DM-electron interactions and non-relativistic (NR) effective field theories (EFTs) (Sec. II), and computing in-medium self-energies to fully account for mixing and screening effects (Sec. III). This is a slightly different strategy than several previous calculations: by matching onto a NR EFT from the beginning instead of taking the NR limit of a relativistic calculation in the end, the power counting relevant for the absorption process becomes more transparent; also, the cryogenic nature of direct detection experiments allows us to perform the in-medium calculation in the zerotemperature limit and avoid the complications of thermal field theory. We will carry out the calculation for three widely-studied bosonic DM candidates: • Vector (e.g. dark photon) DM, which can be produced, for example, by inflationary fluctua-DM type Scalar (φψψ) Pseudoscalar (φψiγ 5 ψ) Related to dielectric? tions [41], by parent particle decays or coherent oscillations after reheating [42][43][44][45], or from a network of cosmic strings [46]. In this case, since the DM couples to electrons via the same vector currentψγ µ ψ as the photon does, its absorption rate is trivially a rescaling of the photon absorption rate.
• Pseudoscalar (e.g. axion-like particle) DM, which can be produced, for example, via the misalignment mechanism [47][48][49], from the decays of topological defects [50][51][52], or by a variety of other mechanisms (see e.g. Refs. [53][54][55][56][57]). While not immediately obvious (since the DM couples to a different current,ψiγ 5 ψ, than the photon does), it has been well-known that also in this case, there is a simple relation between DM and photon absorption [29]. We will recover this result in the NR EFT calculation. It is worth noting that the dominant contribution to NR pseudoscalar DM absorption actually comes from an operator generated at the next-toleading order (NLO) in the 1/m e expansion, because the leading order (LO) operator suffers a suppression by the DM's momentum q.
• Scalar DM, which can be produced via mechanisms similar to pseudoscalar DM mentioned above. It couples to the scalar currentψψ, which at LO coincides with the temporal component of the vector currentψγ 0 ψ. However, as we will see, the LO operator gives a q-suppressed contribution and, as in the pseudoscalar case, the rate is dominated by a NLO operator. Importantly, this NLO operator has a different structure than the photon coupling, and its contribution cannot be simply related to photon absorption, invalidating the data-driven approach.
We make the statements above on the scalar and pseudoscalar DM more concrete in Table I.
The fact that the DM absorption rate is not always relatable to the target material's optical properties highlights the necessity to go beyond the conventional data-driven approach. (The same can be said for DM scattering, for which the data-driven approach based on the dielectric function that has been advocated recently [58][59][60] covers only a limited set of DM interactions.) In this work, we consider two types of targets: • Semiconductor crystals with O(eV) gaps (Sec. IV), focusing on silicon (Si) and germanium (Ge) that are in use in current experiments (DAMIC [61][62][63], EDELWEISS [64][65][66], SENSEI [67][68][69], SuperCDMS [70][71][72][73][74][75][76]). We compute DM absorption rates using first-principles density functional theory (DFT) calculations of electronic band structures and wave functions, which are now publicly available [77]. The numerical calculation builds upon the EXCEED-DM framework [21] and we publish the "absorption" module of the program together with this work [78].
For all the materials under study, we find good agreement between our theoretical calculation and the data-driven approach for the DM models where both are valid, i.e. vector and pseudoscalar DM. This serves as an important validation of our calculations. In the case of scalar DM, we show explicitly how the data-driven approach fails to reproduce the leading contribution, and present our calculated sensitivity projections. In particular, for Al superconductor, our revised projected reach is much more optimistic than that found in Ref. [38], although somewhat weaker than the original estimate in Ref. [24].

II. DARK MATTER COUPLINGS TO NON-RELATIVISTIC ELECTRONS
Since electrons in a detector are non-relativistic, it is convenient to perform the DM absorption calculation in the framework of NR EFT (see e.g. Refs. [79,80] for reviews). In this section, we work through the procedure of matching a relativistic theory of DM-electron interactions onto effective operators involving the NR electron field. The total Lagrangian of interest is Here L ψ is the Standard Model part that includes the electron ψ coupling to electromagnetism, L φ contains the standard kinetic and mass terms of the DM field φ, and we consider the following DM-electron interactions: where we have also indicated the relation between the coupling g and commonly adopted parameters d φee , g aee , κ in the literature. Note that there are two equivalent ways of writing the pseudoscalar coupling that are related by a field redefinition and integration by parts (IBP).
Let us first consider L ψ . Writing the electron field in the relativistic theory as We obtain We now define projection operators which satisfy P 2 ± = P ± , P + P − = P − P + = 0 and (P ± ) † = P ± . By using P ± γ 0 = γ 0 P ± = ±P ± and P ± γ i = γ i P ∓ , we can rewrite Eq. (5) as where ψ ± ≡ P ± ψ NR (thus ψ NR = ψ + + ψ − ). Integrating out the heavy field ψ − at tree level by solving its equation of motion (EOM), we arrive at the EFT for ψ + : where we have used We can readily identify the first four terms in Eq. (9), which come from LO in the 1/m e expansion, as the familiar electromagnetic interactions as in NR quantum mechanics. There are several operators at NLO in the 1/m e expansion, of which we have only written out the one involving ∂ t . This is the last term in Eq. (9), and is the only NLO term that will be relevant in what follows. Importantly, it gives a tree-level contribution to the wave function renormalization of the ψ + field. In NR EFT calculations, it is often convenient to adopt an operator basis where temporal derivatives in the quadratic part of the Lagrangian have been traded for spatial derivatives, so as to eliminate any non-trivial wave function renormalization factors at tree level. The field redefinition needed to go into this basis, at the order we are working here, is This field redefinition does not change the LO Lagrangian (the first four terms in Eq. (9)), but replaces the last term in Eq. (9) by NLO operators that do not contain ∂ t (and hence do not contribute to the wave function renormalization ofψ + ). We will not need the NLO operators for electron couplings to vector fields (photon and dark photon), 1 but the field redefinition in Eq. (11) that modifies the NLO Lagrangian will be important in the cases of scalar and pseudoscalar DM.
We are interested in the case where the photon field A µ consists of an electrostatic background Φ and quantum fluctuations A µ : The normalized NR fieldψ + can be expanded in energy eigenstates of the NR Schrödinger equation: whereĉ I,s are annihilation operators for NR electrons, and Note that the form of the background field in Eq. (12) assumes negligible spin-orbit coupling, in which case the two spin states s = ± for a given I are degenerate. From Eqs. (9) and (12), we can also deduce the electron's coupling to photon quanta A µ at LO in the NR EFT: whereψ † Let us now move on to the DM-electron interaction L int . For vector DM, we can simply replace e A µ → e A µ − g φ µ in the derivation above, and obtain: 1 As a side remark, in the special case of an electrostatic potential, A0 = Φ(x), A = 0, one can check that keeping all the NLO terms reproduces the familiar fine structure correction in NR quantum mechanics: L eff,NLO where the three terms are the relativistic kinetic energy correction, the Darwin term and spin-orbit coupling, respectively.
For the scalar and pseudoscalar cases, since L int contains an operator that has a different structure than all the operators in L ψ , there is no such simple replacement. In principle, we should have included L int when solving the EOM for ψ − in Eq. (8). However, if we are working at leading order in the DM-electron coupling g, it is sufficient to simply substitute Eq. (8) into L int . We therefore obtain, at LO in the NR expansion: We now show that these LO terms are not sufficient to capture the dominant contributions to DM absorption. The point is that our NR EFT is an expansion in ∇ me ∼ v e , and the power counting is such that momenta (and spatial derivatives) count as m e v e and energies (and time derivatives) count as m e v 2 e . For NR absorption, the energy deposition is ω m φ ∼ m e v 2 e . Meanwhile although the momentum transfer formally counts as m e v e , it is in fact much smaller: with v φ ∼ O(10 −3 ). Therefore, when the LO result contains factors of q, we need to work out the NLO terms and see if they may in fact dominate.
From Eq. (17) it is clear that such q suppression is indeed present in the pseudoscalar case. It is perhaps less obvious that the scalar case also suffers a q suppression, and its origin can be understood from charge conservation: the LO operator couples the scalar DM φ to the electron number densitŷ ψ † +ψ + = −ρ e /e (with ρ e the charge density carried by the electron), whose matrix elements vanish in the q → 0 limit because ρ e = q · J e /ω; technically this is manifest via the orthogonality of initial and final state electron wave functions, as we will see later in the paper. 2 Therefore, in both scalar and pseudoscalar cases, we need to expand L int up to NLO where there are several operators. Many of them will not be needed, though, because they are also q suppressed or involve too many fields to contribute to the in-medium self-energies to be computed in the next section. Including only the unsuppressed operators at NLO that contain up to four fields, we have These results were already summarized in Table I (for brevity we dropped the hat onψ + and omitted the last operator in the scalar case in that table -we will see that it gives vanishing contribution to DM absorption in an isotropic medium). The second term in the scalar case, where the DM φ couples replaced by its EOM solution Eq. (8)) and additional terms from the field redefinition in Eq. (11). We will see in the next section that this operator gives the dominant contribution to scalar DM absorption. Pseudoscalar DM absorption is likewise dominated by the NLO

III. IN-MEDIUM SELF-ENERGIES AND ABSORPTION RATES
We now use the NR EFT derived in the previous section to compute DM absorption rates. Generally, the absorption rate of a state can be derived from the imaginary part of its self-energy. In a medium, care must be taken because of mixing effects. If the DM φ mixes with a SM state A in the medium (generalization to the case of mixing with multiple states is straightforward) then the self-energy matrix has to diagonalized to find the in-medium eigenstates: where Simple algebra shows that to O(g 2 ), The DM absorption rate is then derived from the imaginary part of the eigenvalue corresponding to the DM-like state,φ: where the wave function renormalization Zφ = 1 − d Re Πφφ d ω 2 −1 = 1 + O(g 2 ) has been approximated as unity. The total rate per unit target mass is given by where ρ T is the target's mass density, and ρ φ = 0.4 GeV/cm 3 is the local DM energy density. For non-relativistic DM, ω m φ , and ρ φ 1 2 m 2 φ φ 2 0 with the DM field amplitude defined by φ(x, t) = φ 0 cos(q · x − ωt). 3 Technically, the electron fields in the two equivalent expressions of the pseudoscalar coupling, gφψiγ 5 ψ and − g 2me (∂µφ)(ψγ µ γ 5 ψ), are not the same, but are related by a field redefinition. If one derives the NR EFT starting from − g 2me (∂µφ)(ψγ µ γ 5 ψ), this NLO operator is obtained directly from its µ = 0 component. On the other hand, if one derives the NR EFT from gφψiγ 5 ψ, a further field redefinition is needed to eliminate operators involving the background electrostatic potential Φ and arrive at the same operator coefficient shown in Eq. (18).
The calculation of self-energies generally involves two graph topologies: where a blob represents the sum of one-particle-irreducible (1PI) graphs. While we have drawn curly external lines for concreteness, they can each represent a scalar, pseudoscalar or vector. The operators O 1 , O 2 , O that the external fields couple to may carry Lorentz indices, in which case Π and Π inherit these indices. We discuss the calculation of these self-energy diagrams in App. A.
In the cases of interest here, A represents one of the polarizations of the SM photon, and Π AA is directly related to the target's complex conductivity/dielectric function, as discussed further below.
Since Π AA enters the absorption rate formula (Eq. (22)) as long as there is a nonzero mixing Π φA , let us examine this quantity in more detail before specializing to each DM model. The photon self-energy tensor Π µν is defined such that the effective action contains where A j (j = 1, 2, 3) represent the three components of A. As usual, we compute Π µν from the sum of 1PI graphs. From Eq. (25) it is clear that the sign convention here is such that i Π 00 , −i Π 0j and i Π ij are given by the sum of two-point 1PI graphs between A 0 A 0 , A 0 A j and A i A j , respectively.
From the photon-electron couplings in Eq. (15), we obtain Π µν in terms of Π O 1 ,O 2 and Π O defined in Eq. (23) and (24): where the velocity operator v j is defined by Here and in what follows, we suppress the arguments Q µ = (ω, q) of self-energy functions where there is no confusion. To arrive at the expression of Π ij in Eq. (26), we have simplified the spin trace assuming the electron loop does not involve non-trivial spin structures; for example, This assumption is obviously valid for one-loop self-energy diagrams. In the superconductor calculation in Sec. V, we will need two-loop self-energies with an internal phonon line; in that case the electronphonon coupling is spin-independent, so the same simplification applies.
We can write Π µν in terms of its polarization components as follows: where for Q µ = (ω, q) = (ω, qẑ). These are the three photon polarizations in Lorenz gauge Q µ e µ λ = 0, and coincide with the three physical polarizations of a massive vector with m 2 = Q 2 .
We will mostly focus on isotropic target materials in this work, and leave a discussion of the anisotropic case to App. B. For an isotropic medium, the 3 × 3 matrix Π λλ is diagonal: where Π T and Π L are the transverse and longitudinal photon self-energies, respectively. The photon self-energy tensor Π µν therefore has the following form: From the linear response relation J µ = −Π µν A ν 4 and Ohm's law J = σE = σ(iωA − iqA 0 ) we can relate Π T and Π L to the complex conductivity σ, which in turn is related to the complex dielectric ε via σ = iω(1 − ε) [17,23,26]: where Z L = ω 2 /Q 2 . The real part of the conductivity σ 1 ≡ Re σ (the imaginary part of the dielectric) gives the photon absorption rate in medium: We finally note that all the quantities introduced above -the complex conductivity σ and dielectric ε, and photon self-energies Π T , Π L can be simply computed from Π 1,1 : With the photon part of the self-energy calculation completed, we now move on to consider self-energies involving the DM and compute DM absorption rates.

A. Vector Absorption
Since a vector DM couples to electrons in the same way as the photon, albeit with a coupling rescaled by −g/e = −κ, we have Each of the three polarizations of φ mixes with the corresponding polarization of the photon. Therefore, for the transverse (longitudinal) polarization, we simply set . As a result, The total absorption rate for an unpolarized vector DM is obtained by averaging over the three The rate is semi-independent of the momentum transfer (and hence the DM velocity) since Π 1,1 generically scales as q 2 .
The result can also be written in terms of the material's complex conductivity/dielectric: with ε, σ 1 evaluated at ω = m φ , q = 0. One may think of as an in-medium screening factor, which suppresses the absorption rate compared to the obvious rescaling of photon absorption by κ 2 [14,24,31,81].

B. Pseudoscalar Absorption
A pseudoscalar does not mix with the photon due to parity mismatch, 5 and we simply have The pseudoscalar self-energy Π φφ is defined such that the effective action contains Therefore, −i Π φφ is given by the sum of two-point 1PI graphs. From the pseudoscalar coupling in Eq. (18), we find, again after simplifying the spin trace as in Eq. (28): Comparing with Eq. (26), we see that Im Π φφ for a pseudoscalar is closely related to the photon polarization Π µν : Note that the Π 1 term in Π jj is purely real and thus does not appear in the equation above. Also, since ω m e , we can drop the last term. Writing Π µν in terms of Π T and Π L as in Eq. (33) and setting g = g aee , we find For NR absorption, (36)), and therefore, As in the vector DM case, the absorption rate can be written solely in terms of Π  Table I).
We can further recast the pseudoscalar DM absorption rate in terms of the photon absorption rate σ 1 = Re σ = ω Im ε and reproduce the standard result [14,15,24,29]: We remark in passing that pseudoscalar absorption has also been studied in the context of solar axion detection; in that case, the relativistic kinematics ω m φ means that the Im Π L term in Eq. (45)

C. Scalar Absorption
For scalar DM, we need to compute explicitly both Im Π φφ and its mixing with the photon Therefore, −i Π φφ , −i Π 0 φA and i Π j φA are given by the sum of two-point 1PI graphs between φφ, φA 0 and φA j , respectively. From the scalar coupling in Eq. (18) and photon coupling in Eq. (15), we find: As in the photon case, the self-energies are related by the Ward identity Q µ Π µ φA = 0: where we have used the first relation in Eq. (29). One can explicitly check that Eq. (53) holds between the one-loop-level expressions for the self-energies in Eqs. (A7) and (A8).
For an isotropic medium, we must have Π j φA ∝ q j because there is no special direction other than q. 6 So the mixing only involves the photon's longitudinal component. Therefore, Π AA in the rate formula Eq. (22) should be set to Π L = m 2 φ e 2 q 2 Π 1,1 (see Eq. (36)), and Π φA should be set to where we have used the Ward identity to trade q j Π j φA for ω Π 0 φA . Substituting the expressions for Im Π φφ , Π φL and Π L above into Eq. (22), and applying the NR absorption kinematics ω 2 Q 2 = m 2 φ , we find We see that the result for scalar absorption, Eq. (55), depends on Πv2 ,v 2 , Πv2 ,1 , Π 1,v 2 in addition to Π 1,1 . If we had kept only the LO operator φψ † +ψ + in the calculation above, we would obtain Eq. (55) with Πv2 ,v 2 , Πv2 ,1 , Π 1,v 2 set to zero, which coincides with q 2 m 2 φ times the vector DM absorption rate in Eq. (39). Just as in the vector DM case, the contribution of the LO operator φψ † +ψ + to scalar DM absorption is screened due to in-medium mixing [38]. However, the formally NLO operator As we will see in the next two sections, generically Π 1,1 , Πv2 ,1 ∼ q 2 while Πv2 ,v 2 ∼ q 0 . It is thus clear from Eq. (55) that the absorption rate of a NR scalar DM is in fact dominated by the Πv2 ,v 2 term: Importantly, this term (overlooked in several previous calculations of scalar DM absorption [38][39][40]) is not directly proportional to the photon absorption rate and is unscreened. We emphasize that the suppression of LO operator's contribution is specific to the case of non-relativistic DM absorption, where q ω; for absorption of a relativistic scalar (q ω) or scalar-mediated scattering (q ω), the LO operator φψ † +ψ + indeed gives the dominant contribution.
To summarize, in this section we have derived DM absorption rates in terms of in-medium selfenergies of the form Π O 1 ,O 2 , as defined in Eq. (23). (Contributions from the other graph topology, Eq. (24), have been eliminated using the Ward identity.) Both vector and pseudoscalar absorption involve a single self-energy function Π 1,1 ∝ Π L (see Eqs. (39) and (46)), and the rates can be simply 6 We note in passing that the Π v j term in Π j φA is q independent and must therefore vanish in an isotropic medium. This is why we have omitted the φ A · ψ † + ← → ∇ψ+ operator in Eq. (18), which only contributes to this term, from Table I. related to the (complex) conductivity/dielectric (see Eqs. (40) and (47)). In these cases, the datadriven approach based on the measured conductivity/dielectric is viable, and we can also use optical data to calibrate our theoretical calculations based on DFT or analytic modeling. On the other hand, for scalar DM absorption, additional self-energy functions Πv2 ,v 2 , Πv2 ,1 , Π 1,v 2 enter (see Eq. (55)), and the rate is not directly related to photon absorption. In this case, the data-driven approach fails and theoretical calculations are needed.
In the next two sections, we compute the self-energies Π 1,1 , Πv2 ,v 2 , Πv2 ,1 , Π 1,v 2 in crystal and superconductor targets, respectively, which then allow us to derive the absorption rates of vector, pseudoscalar and scalar DM in these targets. Our main results for Si, Ge and Al-superconductor and comparing the sizes of the terms. Here , while the remaining terms define Rv2 ,1 . Next, Fig. 2 shows the projected reach for the pseudoscalar and vector DM models, where we see good agreement between our theoretical calculations (solid curves) and rescaled optical data (dashed curves). Lastly, Fig. 3 shows our calculated reach for scalar DM and compares the Al-SC results with previous work [14,38]. These results will be discussed in detail in the following sections.

IV. DARK MATTER ABSORPTION IN CRYSTALS
In this section, we specialize to the case of crystal targets that are described by band theory. It which in real and momentum space read, respectively: where the sum runs over all reciprocal lattice vectors G. These are related by upon applying the standard dictionary between discrete and continuum expressions: We now examine the matrix element i , k | O 1 e iq·x |i, k involved in Eq. (A7) for thev 2 and 1 operators; i, k| O 2 e −iq·x |i , k is completely analogous. For thev 2 operator, we simply obtain For NR absorption in the mass range of interest here, m φ 100 eV, the momentum transfer q ∼ 10 −3 m φ ∼ meV m φ eV is well within the 1BZ (O(keV)). This implies that Umklapp processes where G = G do not contribute, so (lattice) momentum conservation simply dictates k = k + q. At leading order in q we can set k = k, and Eq. (61) simplifies to superconductor (Al-SC) targets for the vector and pseudoscalar DM models defined in Eq. (3), assuming 1 kg-yr exposure. We compare our theoretically calculated reach (solid) against the data-driven approach utilizing the target material's measured conductivity/dielectric [83,84] (dashed). For Si and Ge, the data-driven approach was taken in previous works [14,15], with which we find good agreement. For Al-SC, our theoretical calculation reproduces the results in Ref. [24] (dotted) up to the choice of overall normalization factor. Also shown are existing direct detection limits from XENON10/100 [15], stellar cooling constraints from the Sun (assuming Stückelberg mass for vector DM) [85] and white dwarfs (WD) [86], and pseudoscalar couplings corresponding to the QCD axion in KSVZ and DFSZ (for 0.28 ≤ tan β ≤ 140) models [87].
For the 1 operator, additional care is needed since i , k | e iq·x |i, k vanishes in the q → 0 limit: |i , k and |i, k are distinct energy eigenstates and therefore orthogonal. At O(q), we have i , k | e iq·x |i, k iq · i , k | x |i, k . A numerically efficient way to compute this matrix element is to trade the position operator for the momentum operator via its commutator with the Hamiltonian H = p 2 2me + V (x): Substituting in the wave functions, we find: where It is convenient to define the following crystal form factors, via which the Bloch wave functions  Fig. 2, the projections here cannot be derived from the target's optical properties. Differences compared to Hochberg et al. [24] and Gelmini et al. [38] in the Al-SC case are discussed in detail in Sec. V. Also shown are existing constraints from fifth force [88] and red giant (RG) cooling [89]. enter DM absorption rates (at leading order in q): Note that they differ from the crystal form factor used in spin-independent DM scattering [12,17,21]: The absorption kinematics simply set the k and G vectors of the initial and final states to be the same; also, powers of (k + G) appear as follows from the effective operators.
The crystal form factors defined above allow us to write the self-energies in a concise form. For the operators 1 andv 2 , the spin trace is trivial and simply yields a factor of two. Each pair of valence/conduction states between which a transition can happen contributes to two terms in the sum over electronic states, because either i, k or i , k can be a valence or conduction state. Combining the two terms for each pair, we obtain where δ → 0 + . We see explicitly that Π 1,1 ∼ q 2 and Πv2 ,v 2 ∼ q 0 , as already alluded to in Sec. III. The other two self-energies, Πv2 ,1 and Π 1,v 2 , take the form of q · F + O(q 2 ), where F is a target-dependent function that involves f i i, k and f i i, k . In the absence of a special direction, we must have F = 0 and therefore, Πv2 ,1 , Π 1,v 2 ∼ O(q 2 ). Working out the leading O(q 2 ) contribution to these self-energies To calculate the DM absorption rates and make sensitivity projections, we use DFT-computed electronic band structures and wave functions for Si and Ge [77], including all-electron reconstruction up to a cutoff of 2 keV; see Ref. [21] for details. We adopt the same numerical setup as the "valence to conduction" calculation in Ref. [21], and include also the 3d states in Ge as valence (treating them as core states gives similar results). The finite resolution of the k-grid means we need to apply some kind of smearing to the delta functions coming from the imaginary part of Eqs. (67) and (68). This is done in practice by setting δ in Eqs. (67) and (68) to a finite constant 0.2 eV, which we find appropriate for a 10 × 10 × 10 k-grid for the majority of the DM mass range. We implement our numerical calculation as a new module "absorption" of the EXCEED-DM program [78].
We present the projected reach for the three DM models in Figs. 2 and 3, assuming 3 events (corresponding to 95% CL) for 1 kg-yr exposure without including background, together with existing constraints on these models for reference. The solid curves are our theoretical predictions; they are obtained using the rate formulae Eqs. (39), (46) and (56)  For vector and pseudoscalar DM, we can alternatively take the data-driven approach, using Eqs. (40) and (47), respectively, to derive the rate from the measured conductivity/dielectric. As in Ref. [14,15], we use the measured optical data from Ref. [83]. Results from this data-driven approach are shown by the dashed curves; they are the same as in Ref. [14,15] upon inclusion of backgrounds.
For Si, the solid and dashed curves are very close to each other for m φ 3 eV; the theoretical calculation (solid curves) systematically overestimates the rate as m φ approaches the band gap (1.2 eV) because of the smearing procedure discussed above. For Ge, we see the same systematic discrepancy close to the band gap (0.67 eV); also, the theoretical calculation predicts a sharper plasmon peak (corresponding to a smaller Im Π 1,1 near the plasmon frequency) compared to data. Aside from these issues, we view the overall good agreement between the solid and dashed curves in the vector and pseudoscalar cases as a validation of our DFT-based theoretical calculation in the majority of DM mass range. Importantly, this gives credence to the reach curves we have calculated in the scalar DM case, where the data-driven approach does not apply, though one has to keep in mind that our calculation systematically overestimates the rate for DM masses below about 3 eV because of the smearing issue.

V. DARK MATTER ABSORPTION IN SUPERCONDUCTORS
We now turn to the case of conventional superconductors described by BCS theory. Within this simple free-electron model, the self-energies Π O 1 ,O 2 are real at one-loop level; it is easy to see that two electronic states differing by energy ω and momentum q cannot be both on-shell when ω q. Therefore, the leading contribution to the imaginary part arises at two loops, and we have While these are derived for normal conductors, we expect them to carry over to the superconductor case; proportionality to n e (the total number of electronic states within the Fermi sphere) implies insensitivity to deformations within O(∆) of the Fermi surface. We also note in passing that, via Eq. (36), we obtain the familiar result for the photon self-energies [90,91]: p ≡ e 2 ne m * is the plasma frequency squared. For the imaginary part Im we expect the dominant contribution to come from two-loop diagrams with an internal phonon line for a high-purity sample (otherwise impurity scattering may also contribute). These are associated with φ (or γ) + e − → e − + phonon processes by the optical theorem, and can be computed by the standard cutting rules, as we detail in App. A 3. We model the (acoustic) phonons with a linear dispersion, ω q = c s q where c s is the sound speed, and neglect Umklapp processes which amounts to imposing a cutoff on the phonon momentum, q max = q D ≡ ω D /c s with ω D the Debye frequency. The electron-phonon coupling, in our normalization convention, is given by C e-ph q √ 2ω q ρ T [24,92,93], with C e-ph ∼ O(ε F ) a constant with mass dimension one. Accounting for the superconducting gap, we obtain, for ω q: where E(z) = 1 0 dt 1−z 2 t 2 1−t 2 is the complete elliptic integral of the second kind. For energy depositions much higher than the gap, ω 2∆, the elliptic integral E(1) = 1 drops out and we reproduce the results for a normal conductor; see App. A 3 for details.
With the expressions of self-energies above, we can use Eqs. (39), (46) and (55) to calculate the absorption rates for vector, pseudoscalar and scalar DM. We consider an aluminum superconductor (Al-SC) target, for which the relevant material parameters are listed in Table II. We use the same numerical values as in Ref. [24] for ε F , ω p , ∆, ω D , c s , and determine the electron-phonon coupling C e-ph from resistivity measurements [94,95] as explained in App. A 3. For scalar DM, we again confirm the dominance of the Rv2v2 term in Eq. (57), as seen in Fig. 1  to Si and Ge to lower m φ . The solid red curves are obtained from the self-energy calculations discussed above; the underlying model has a UV cutoff ω max ∼ 0.5 eV where we truncate the curves.
Low-temperature conductivity data are available between 0.2 eV and 3 eV [84]. For the vector and pseudoscalar DM models, we also present the reach following the data-driven approach in this mass range (dashed curves), obtained by using Eqs. (40) and (47) with σ 1 (= Re σ = ω Im ε) taken from Ref. [84] and Re ε set to 1 − ω 2 p ω 2 . Between 0.2 eV and 0.5 eV where both theoretical (solid) and datadriven (dashed) predictions are shown, they are in reasonable agreement, with the latter stronger by about 40% for both κ and g aee at 0.2 eV. The difference is presumably a result of approaching the UV cutoff of the theoretical calculation, and possibly also the neglect of Umklapp contributions. For scalar DM, the data-driven approach is not viable, and we present our theoretical prediction up to 0.5 eV. We also show the reach curves obtained in the previous literature [24,38] for comparison, and discuss the differences in what follows.
Comparison with previous calculations. The calculation of DM absorption in superconductors was first carried out in Ref. [24], where the 2→2 matrix element for φ + e − → e − + phonon was evaluated at leading order in q. For vector and pseudoscalar DM, our results agree with Ref. [24] as seen in Fig. 2, up to a minor numerical prefactor understood as follows. Ref. [24] chose the value of the electron-phonon coupling C e-ph such that the photon absorption rate (i.e. conductivity σ 1 ) matches the experimentally measured value at ω = 0.2 eV. In this work, we instead determine C e-ph via the λ tr parameter following Refs. [94,95], which results in a slightly lower value and hence the slight mismatch observed in Fig. 2.
The more significant numerical difference in the scalar case between our results and Ref. [24], as seen in Fig. 3, can be traced to two sources. First, the numerically dominant effect is that Ref. [24] did not distinguish m * and m e , while we have kept the vacuum mass m e in the operator coefficients and used the effective mass m * for the electron's dispersion and phase space; the two masses differ by about a factor of three in Al-SC. Note that the difference between m e and m * does not affect the vector and pseudoscalar absorption rates as they only depend on Π 1,1 , which is independent of m * /m e . Second, Ref. [24] dropped a factor of (1 − x) 2 in the scalar absorption matrix element when taking the soft phonon limit; this results in an O(1) difference on the projected reach that is numerically subdominant. One can easily verify these two points by evaluating the integral in Eq. (A38) using 3 in the last line; this would reproduce the analytic relation presented in Ref. [24] between scalar and photon absorption rates in the limit ω 2∆.
More recently, Ref. [38] revisited scalar DM absorption and claimed that in-medium effects lead to a significantly weaker reach. We reiterate that while in-medium mixing with the photon screens the contribution from the LO operator 1, the leading contribution to scalar absorption comes instead from the NLO operatorv 2 that is not screened. In fact, the screening factor in Ref. [38] was (correctly) derived for the 1 operator but inconsistently applied to the dominant contribution coming from thē v 2 operator as obtained in Ref. [24]. As a result, Ref. [38] significantly underestimated the reach as we can see from Fig. 3.

VI. CONCLUSIONS
In this paper we revisited the calculation of electronic excitations induced by absorption of bosonic DM. Specifically, we focused on O(1 -100) eV mass DM for Si and Ge targets that are in use in current experiments, and sub-eV mass DM that a proposed Al superconductor detector will be sensitive to.
We utilized an NR EFT framework, where couplings between the DM and electron in a relativistic theory are matched onto NR effective operators in a 1/m e expansion. We then computed absorption rates from in-medium self-energies, carefully accounting for mixing between the DM and the photon.
For crystal targets like Si and Ge, we used first-principles calculations of electronic band structures and wave functions based on density functional theory, and implemented the numerical rate calculation as a new module "absorption" of the EXCEED-DM program [21,78]. For BCS superconductors, we adopted an analytic model as in Refs. [22][23][24]  For scalar DM, however, we showed that the dominant contribution is not directly related to photon absorption. One therefore cannot simply rescale optical data to derive the DM absorption rate.
Importantly, the familiar coincidence between scalar and vector couplings,ψψ ψ γ 0 ψ, holds only at leading order in the NR EFT. For non-relativistic scalar DM φ, matrix elements of the leading order operator are severely suppressed by the momentum transfer q ∼ 10 −3 m φ . The dominant contribution comes instead from a different operator that is formally next-to-leading-order in the NR EFT expansion, and does not suffer from in-medium screening. We presented reach projections for scalar DM based on our theoretical calculations. Notably, for Al superconductor, the reach we found is much more optimistic than the recent estimate in Ref. [38].
It is straightforward to extend the calculation presented here to anisotropic targets and materials with spin-dependent electronic wave functions (as can arise from spin-orbit coupling); we will investigate this subject in detail in an upcoming publication. Another future direction is to calculate phonon and magnon excitations from DM absorption via in-medium self-energies in a similar EFT framework, refining and extending the calculation in Ref. [96]. Finally, in-medium self-energies are also relevant for DM detection via scattering; one can carry out a calculation similar to what we have done here, but in a different kinematic regime, to include in-medium screening corrections in the study of DM-electron scattering via general EFT interactions [97]. Foundation.
Appendix A: Self-energy Calculations

General Result for the One-loop Self-energy
At one-loop level, the self-energies defined in Eqs. (23) and (24) are given by where the external states (drawn with curly lines for concreteness) can be either spin-0 or spin-1, and the internal electronic states I, I are summed over. Using the in-medium Feynman rules (see e.g. Ref. [93]) we obtain, for the first diagram: where V is the total volume, "tr" represents the spin trace, and δ I ( ) ≡ δ sgn(ε I ( ) − ε F ) with δ → 0 + .
Note that the iδ prescription for electron propagators is different from the vacuum theory, and depends on whether the state is above or below the Fermi energy ε F ; using the correct iδ prescription is crucial for ensuring causality. Meanwhile, the matrix elements coming from the vertices are and likewise for I| O 2 e −iq·x |I . Here O 1,2 are matrices in spin space, and may involve spatial derivatives acting on the electronic wave functions. For example, for the velocity operator defined in Eq. (27) (which is proportional to the identity matrix in spin space), we have We can evaluate the energy integral in Eq. (A2) by examining the pole structure of the integrand in the complex plane. If δ I and δ I have the same sign (i.e. if both I and I are above or below the Fermi energy), the two poles are on the same side of the real axis and they have opposite residues; the integral therefore vanishes upon closing the contour via either +i∞ or −i∞. So we must have one state above the Fermi energy and one below it, in which case there is one pole on each side of the real axis; closing the contour via either +i∞ or −i∞ to pick up the residue at one of the poles, we obtain Here ω I I ≡ ε I − ε I , and δ → 0 + . All cases discussed above can be concisely summarized as: where f I , f I are the occupation numbers (equal to 1 for states below the Fermi energy, 0 for states above it), and δ I I ≡ δ sgn(ω I I ). We therefore obtain As the simplest example, setting O 1 = O 2 = 1 in Eq. (A7), we obtain Π 1,1 , and hence the dielectric via Eq. (36), which reproduces the familiar Lindhard formula (see e.g. Ref. [98] and recent discussions in Refs. [58,59]). Now move on to the second diagram in Eq. (A1). While we have shown in Sec. III that contributions to absorption rates from this diagram can be eliminated using the Ward identity, we present its result here for completeness and also to allow for an explicit check of the Ward identity. In this diagram, the electron propagator starts and ends at the same time point and time-ordering becomes ambiguous.
The correct prescription is to take the normal-ordered product of creation and annihilation operators, and the loop is simply proportional to the electron number operator [93]. Again writing the result in terms of occupation number f I , we find Note that Π O is purely real at all orders. With Eqs. (A7) and (A8) one can readily verify the relations implied by the Ward identity, Eqs. (29) and (53).

Real Part of the One-loop Self-energy in a Metal
We now apply Eq. (A7) to the case of a metal. As discussed in Sec. V, we model the electrons near the Fermi surface of a metal as free quasiparticles with an effective mass m * and energy eigenstates labeled by momentum. The sum over I, I becomes integrals over k, k , and we have Note that the iδ in the denominator is irrelevant since the intermediate states cannot go on-shell and Let us first consider Π 1,1 . For the matrix element part, we have Therefore, Expanding in small q and integrating by parts, we find where the gradients are in k space, and n e = 2 d 3 k (2π) 3 f k is the free electron density. We can calculate Πv2 ,1 in a similar way. The matrix element part again yields a momentumconserving delta function, and the integrand can then be expanded in small q. We find where we have used f k = Θ(k F − k), and n e = 2 (2π) 3

Imaginary Part of the Two-loop Self-energy in a Metal
The one-loop self-energies calculated above are purely real: both electrons cannot go on-shell if their energies and momenta differ by Q µ = (ω, q) with ω q. The leading contribution to Im Π O 1 ,O 2 (Q) comes from two-loop diagrams with an internal phonon line. In this section, we compute them first in the case of a normal conductor, and then discuss the corrections needed in the superconductor case when ω approaches the gap 2∆.
Cut diagrams. There are three contributing diagrams: Here each propagator is labeled by a four-momentum that consists of the energy it carries and the momentum label of the electron or phonon state. In each diagram, we denote four-momentum flowing into the O 1 vertex from the electron propagator as K µ = (ε, k), and the phonon four-momentum (with direction indicated by the arrow) as Q µ = (ω , q ). The electron and phonon propagators are denoted by iG and iG ph , respectively, with where δ ε = δ sgn(ε−ε F ), δ → 0 + , and ω q = c s q . The electron-phonon vertex y q = For ω > 0, the on-shell condition requires ω > 0, ε < ε F and ε + ω − ω > ε F ; this corresponds to a process where an electron jumps out of the Fermi sphere by absorbing Q µ = (ω, q) while emitting a phonon to conserve momentum. We therefore obtain where it is understood that ε (the energy components of K) is set to k 2 2m * in the final expression. The second diagram, Eq. (A16), is completely analogous. Cutting the propagators G K+Q , G ph Q and G K+Q , we obtain where we have shifted the integration variable k → k − q to arrive at the last line.
For the last diagram, Eq. (A17), there are two possible cuts: through G K , G ph Q , G K+Q−Q and through G K+Q , G ph Q , G K−Q . Carrying out the same procedure as above, we obtain where we have shifted the integration variable k → k + q and then changed q → −q (assuming the phonon energies ω q and electron-phonon couplings y q depend only on the magnitude but not the direction of q ) in the second term.
Adding up Eqs. (A22), (A23) and (A24), we obtain Small q expansion. As in the previous section, we expand the integrand in small q. The electron propagators become: where we have used the energy-conserving delta function to eliminate ω q in G K−Q . Therefore, at leading order in q , where an identity operator in spin space is understood, and we have again used energy conservation to simplify the expression in the O 1 =v 2 case. Note in particular how the O(q 0 ) terms cancel in the case of O 1 = 1, such that this LO operator gives a q-suppressed contribution. The other is completely analogous, so we obtain, after taking the spin trace (which simply yields a factor of two) and substituting in , ω q = c s q : Including the gap. We have presented the calculation of cut diagrams assuming a normal metal for simplicity. Accounting for pairing of electrons in the BCS theory introduces a slight modification in the final result in the form of a coherence factor [94]. Concretely, for the imaginary part of two-loop self-energies computed above, this amounts to replacing in Eq. (A29), where we have abbreviated k + q − q ≡ k and defined k ≡ k 2 2m * − ε F , E k ≡ 2 k + ∆ 2 (and similarly for k , E k ). The energy of the electron-hole pair is therefore constrained to be The k integral. We now perform the k integral: The integrand depends only on the magnitude of k and the angle θ between k and q − q. So the azimuthal angle integral simply yields a factor of 2π and we can use the δ function to perform the integral over cos θ. The argument of the δ function has two roots in cos θ (corresponding to k = ±| k |), both of which are within the range [−1, 1] in most of viable phase space. Noting that Changing the integration variable from k to E k , we find where a factor of two comes from combining contributions from the two values of k above and below k F that correspond to the same E k , and we have abbreviated E k , E k , k , k to E, E , , , with E = ω−ω q −E. The integral over E can be reduced to elliptic integrals via E = 1 2 ω−ω q +t (ω−ω q −2∆) : to simplify notation, and are the complete elliptic integrals of the first and second kind, respectively. In the ∆ → 0 limit, corresponding to a normal conductor, we have α → 1, β → 0, E(1) = 1, and I → m 2 * (ω−ω q ) 2π|q −q| .
The q integral. The remaining integral over the phonon momentum is Expanding q |q −q| = 1 + q ·q q 2 + . . . and keeping the leading nonvanishing term, we can easily carry out the angular integration. Finally, changing the radial integration variable to x = csq ω , we obtain where the upper limit is set by the requirements ω − ω q ≥ 2∆ and ω q = c s q ≤ ω D (Debye frequency).
Determination of C e-ph . We use resistivity measurements [95] to determine C e-ph . In Refs. [94,95], a parameter λ tr is introduced for the electron-phonon coupling, which is defined by The function α 2 tr F (ω ) is in turn defined from the conductivity of a normal conductor, Note that the normalization convention in Ref. [94] is such that 4πσ 1 there equals σ 1 in our notation.
However, it is straightforward to extend the calculation to anisotropic targets. In this appendix, we discuss the modifications needed to go beyond the isotropic limit.
First, the in-medium photon self-energy matrix Π λλ may have nonzero off-diagonal entries, and its eigenvalues can be found by diagonalization [26]: A DM state φ may mix with all three photon polarizations, and Eq. (22) generalizes to where Π φi is obtained from Π µ φA by first projecting onto e µ ±,L and then rotating into the diagonal basis. For vector DM φ, the same rotation in Eq. (B1) diagonalizes also the φφ and φA self-energy matrices. So each of the three polarizations of φ mixes only with the one corresponding photon polarization, and we simply replace Π T,L in Eq. (22) by Π 1,2,3 and average over the three polarizations to obtain the rate in the anisotropic case: For pseudoscalar DM, still assuming spin-degenerate electronic states, we obtain from Eq. (44): which generalizes Eq. (45). Note that while anisotropy allows for a nonzero mixing between the DM φ and the photon (via its coupling to the electron's magnetic dipole), its contribution to absorption rate is at O(q 2 ) and negligible. On the other hand, if the electronic states are not spin-degenerate (e.g. due to spin-orbit coupling), one would need to explicitly compute additional matrix elements of the spin operator Σ between the spin part of the wave functions, and the absorption rate cannot be written in terms of components of the photon self-energy matrix. Also, mixing between the DM and the photon becomes relevant in this case.
For scalar DM, anisotropy may introduce mixing with all three photon polarizations, and Eq. (B2) applies. The final result, however, is still expected to be dominated by the Πv2 ,v 2 term from Π φφ , and we therefore have the same formula, Eq. (56), as in the isotropic case: