Absorption of Vector Dark Matter Beyond Kinetic Mixing

Massive vector particles are minimal dark matter candidates that motivate a wide range of laboratory searches, primarily exploiting a postulated kinetic mixing with the photon. However, depending on the high energy field content, the dominant vector dark matter (VDM) coupling to visible particles may arise at higher operator dimension, motivating efforts to predict direct detection rates for more general interactions. Here we present the first calculation of VDM absorption through its coupling to electron electric (EDM) or magnetic (MDM) dipole moments, which can be realized in minimal extensions to the Standard Model and yield the observed abundance through a variety of mechanisms across the eV\,-\,MeV mass range. We compute the absorption rate of the MDM and EDM models for a general target, and then derive direct detection constraints from targets currently in use: Si and Ge crystals and Xe and Ar atoms. We find that current experiments are already sensitive to VDM parameter space corresponding to a cosmological freeze-in scenario, and future experiments will be able to completely exclude MDM and EDM freeze-in models with reheat temperatures below the electroweak scale. Additionally, we find that while constraints on the MDM interaction can be related to constraints on axion-like particles, the same is not true for the EDM model, so the latter absorption rate must be computed from first principles. To achieve this, we update the publicly available program EXCEED-DM to perform these new calculations.


I. INTRODUCTION
Light vector particles are economical extensions to the Standard Model (SM) that require no stabilizing symmetries or mediator particles to account for the dark matter (DM) in our universe. The cosmological abundance of vector DM (VDM) can arise through a variety of mechanisms [1][2][3][4][5][6][7][8][9][10]. Minimally, VDM can be produced gravitationally through quantum fluctuations during inflation [11]. Alternatively, the VDM abundance can arise through its SM interactions via the "freeze-in" mechanism [12,13] which relates VDM production at early times to observable signatures in terrestrial laboratories.
The most commonly studied VDM interaction is kinetic mixing with the SM photon through the V µν F µν operator, where V µ is the VDM field, with mass m V , and V µν ≡ ∂ µ V ν − ∂ ν V µ is its field strength tensor. This interaction can populate the VDM through the infrared (IR) freeze-in mechanism, in which the DM is initially absent at reheating and builds up through sub-Hubble interactions as the universe expands. While this mechanism is elegant and predictive, it is excluded for nearly all m V by a combination of direct and indirect detection searches [14], so there is motivation to explore alternative possibilities.
In the absence of kinetic mixing, the leading, viable, V -SM interactions are the electric (EDM) and magnetic * krnjaicg@fnal.gov † ttrickle@fnal.gov dipole moment (MDM) operators, which have mass dimension five, σ µν = i 2 [γ µ , γ ν ], and Ψ is a charged SM fermion field represented as a Dirac spinor in the broken electroweak phase. Such operators can be the leading VDM interaction with SM particles if suitable new states are integrated out at energy scales above E > 1/d E,M (for a concrete model see Refs. [15,16]).
The phenomenology of the MDM model was fist studied in Ref. [17], where it was shown that VDM with keV < ∼ m V < ∼ MeV can viably freeze-in, while avoiding indirect detection, and warm DM constraints. Since the MDM operator has mass dimension five, the cosmological abundance depends on the reheat temperature, T RH . This UV sensitivity makes it possible to viably freeze-in VDM by leveraging the potentially large T RH to enhance cosmological production, while evading indirect detection constraints that exclude freeze-in through kinetic mixing.
While Ref. [17] mainly studied the indirect detection bounds on the MDM interaction, direct detection constraints were left for future work. In this paper, we extend this analysis to study VDM absorption onto atomic and crystal targets: • Atomic Targets: For keV < ∼ m V < ∼ MeV, large exposure liquid noble experiments, e.g., XENON [21][22][23] and DarkSide [26,27], are expected to be sensitive to absorption events when the VDM model couples to the electron. We will re- (1), in crystal Si (red) and Ge (blue) targets, and atomic Xe (green) and Ar (purple) targets. Crystal targets assume an exposure of 1 kg · yr, and atomic targets assume an exposure of 10 ton · yr. Dashed lines are constraints using rescaled photon absorption data. Photon absorption data for the crystal targets is a combination of Refs. [18,19], and data for the atomic targets is a combination of Refs. [19,20], shown in Fig. 3. Previous constraints on ALPs from XENONnT [21], XENON1T [22], XENON10/100 [23], and SuperCDMS [24] (shaded teal, orange, cyan, red, respectively) have been recast by converting those constraints to photon absorption cross sections, and using Eq. (26). Indirect detection bounds on V → 3γ from INTEGRAL [17,25] are shown in shaded blue. Gray dashed lines apply if the DM is produced via the freeze-in mechanism [13,17]. The warm DM limit (WDM) is taken from Ref. [17], and the lines labelled by T RH correspond to the necessary reheat temperature to generate the relic abundance via Eqs. (7) and (8).
fer to these targets as "atomic targets," since we approximate them as a collection of individual atoms, such that the total absorption rate is a simple sum of contributions from each atom. Atomic targets are especially interesting because they close the open window between the low mass (m V ∼ keV) warm DM constraints, and the indirect detection bounds which are dominant at higher (m V ∼ MeV) masses. As we will see, these targets also play a key role in testing the predictive freeze-in scenarios for a wide range of T RH .
• Crystal Targets: For lower masses, m V < ∼ keV, freeze-in is not a viable production mechanism since the DM would be too warm. However, there are a variety of production mechanisms which can populate the DM; see Sec. II A for more details. For these DM models, the energy levels of atomic targets are no longer suitable for efficient VDM absorption. Thus, for eV < ∼ m V < ∼ keV we also compute absorption rates and extract constraints for crystal Si and Ge targets with lower energy thresholds. These targets are utilized in several current and future experiments including, CDEX [28], DAMIC [29][30][31][32][33], EDELWEISS [34][35][36], SENSEI [37][38][39], and SuperCDMS [40][41][42].
In principle Ψ, in Eqs. (1) and (2), could be any charged SM fermion. However, since our main focus is VDM absorption onto direct detection targets, for the remainder of this work, we will only consider the electron coupling. Furthermore, throughout our analysis, we treat the EDM and MDM cases separately, though our results generalize easily to scenarios in which both d M and d E are nonzero.
Previous VDM absorption rate calculations [14,[43][44][45][46][47] have been focused on the kinetically mixed scenario, for which the absorption rate is simply related to the photon absorption rate. However, this relation does not hold in general; for generic interactions, the absorption Ge (kg · yr) Xe (10 ton · yr) Ar (10 ton · yr) FIG. 2: Projected 95% C.L. constraints (3 events, no background) on the d E parameter in the EDM model, Eq. (2), in crystal Si (red) and Ge (blue) targets, and atomic Xe (green) and Ar (purple) targets. Crystal targets assume an exposure of 1 kg · yr, and atomic targets assume an exposure of 10 ton · yr. Shaded regions with dashed outline are expected direct detection constraints derived by rescaling the constraints in analogy with Fig. 1. Indirect detection bounds from INTEGRAL [25] are shown in shaded blue. Gray dashed lines apply if the DM is produced via the freeze-in mechanism [13,17]. The warm DM limit is taken from Ref. [17], and the lines labelled by T RH correspond to the reheat temperature that yields the observed VDM abundance via Eqs. (7) and (8).
rate must be calculated from first principles [45,46]. For our operators of interest, Eqs. (1) and (2), we generalize the procedure outlined in Ref. [45] and find that, while the MDM absorption rate can be related to the photon absorption rate, the same is not true for the EDM interaction, which requires a dedicated calculation.
Moreover, we show that for 100 keV < ∼ m V even the absorption rate of the familiar kinetically mixed model cannot be related to the photon absorption rate due to the kinematic mismatch (see Sec. III). We perform the first principles absorption calculation by extending the publicly available code EXCEED-DM [47][48][49] with support for atomic targets, and make the modifications publicly available as well.
This paper is organized as follows. In Sec. II we discuss cosmological production mechanisms that can populate VDM in the early universe. In Sec. III we derive the electronic absorption rate of VDM in the MDM and EDM models. In Sec. IV we begin by comparing our first principles calculation of the VDM absorption rate to previously computed photoelectric cross section, and then we compute, and discuss, the direct detection constraints for both the MDM and EDM models. Lastly, in Sec. V we summarize our results and discuss how future work may extend experimental sensitivity to these scenarios.

II. COSMOLOGICAL PRODUCTION
Viable DM candidates that allow for absorption processes are generically out of equilibrium in the early universe, as the interaction strengths required to thermalize with the SM also induce particle lifetimes much shorter than the age of the universe. Thus, in this section, we briefly survey a variety non-thermal VDM production mechanisms, categorized according to whether or not the abundance arises from SM interactions.

A. Production From Additional BSM Fields
In a broad class of models, the non-thermal VDM abundance depends on the details of the very early universe. For example, if the V mass is nonzero during inflation and light compared to the Hubble rate, there is a cosmic abundance of longitudinally polarized vectors arising from inflationary fluctuations [11], where H I is the inflationary Hubble scale, Ω i ≡ ρ i /ρ c is the present day abundance fraction of species i, ρ c ≈ 4 × 10 −47 GeV 4 is the critical density, and Ω DM = 0.264 [50].
If V couples to additional fields that undergo non-trivial evolution in the early universe, the VDM production rate can be further be enhanced through parametric resonance, which can yield the observed DM abundance even if gravitational production through inflationary fluctuations is inefficient [7,10]. The VDM abundance can also arise from initial conditions via pre-inflationary misalignmnent [3,9,51]. However, as noted in Ref. [3], misalignment is inefficient at producing VDM abundance unless the vector is non-minimally coupled to gravity. Since these mechanisms populate VDM independently of its coupling to the SM, we remain agnostic to the UV details of such scenarios. Throughout this work, we assume that one of these mechanisms suffices to produce the abundance -particularly in the low mass (m V < ∼ keV) regime, where SM freeze-in (discussed below) production is excluded by structure formation bounds on warm DM.

B. Production From SM Freeze-In
The freeze-in mechanism postulates that DM is initially not populated when the SM radiation bath is created after inflation. Self-consistency requires the DM-SM interaction rate to be sub-Hubble at this time, so that DM does not thermalize with visible matter. In this class of models, the V abundance is, where ⟨Γ V ⟩ is the thermally averaged VDM production rate,n V is the number density the V particles would have if they were in chemical equilibrium at temperature T , H is the Hubble expansion rate during radiation domination, s is the entropy density, and a zero subscript represents a present day quantity. The integration range in Eq. (4) spans from the reheat temperature T RH to the temperature at which freeze-in production halts, which typically satisfies T IR = max(m V , m SM ), where m SM is the mass of the main SM species driving freeze-in production.

Excluding Renormalizable Freeze-In
Massive vectors that kinetically mix with the photon can be frozen in through this same interaction, while maintaining a cosmologically long lifetime, for couplings that yield the observed DM abundance. For m V < 2m e , the dominant decay channel is V → 3γ, which can be cosmologically metastable due to the sharp phase space suppression in the width for this process (∝ m 9 V /m 8 e ). However, the keV < ∼ m V < ∼ 2m e window is almost fully excluded by a combination of X-ray and direct detection limits [14]. For lighter (m V < ∼ keV) masses, these direct bounds can be evaded, but in this regime the VDM is too warm for viable structure formation.
If the vector particle is the gauge boson of an abelian SM extension -e.g., gauged B − L or L i − L j , where B and L are the baryon and lepton number, respectivelyit can couple directly to visible particles in the absence of kinetic mixing (see Ref. [52] for a review). However, in the absence of additional new field content at low energies, all anomaly free U (1) extensions require V couplings to neutrinos. Thus, for gauge couplings that would produce the observed DM abundance (e.g., g ∼ 10 −11 for the ALP and kinetically mixed dark photon models [14]), V →νν decays are prompt on cosmological timescales, so the vector is not a viable DM candidate.
Similar considerations apply to VDM whose population freezes in through other dimension four operators (e.g., V µΨ γ µ γ 5 Ψ). Since IR dominated freeze-in predicts a one-to-one correspondence between production and late time decay, the irreducible loop-level decay V → 3γ decay is comparably constrained by the same X-ray bounds that exclude kinetic mixing. Thus, any viable model of VDM freeze-in through SM interactions must involve operators beyond mass dimension four.

Freeze-In Through Dipole Operators
In light of the above considerations, we now consider freezing in VDM through MDM and EDM interactions. In the early universe, if electroweak symmetry is restored at high temperatures, the interactions in Eqs. (1) and (2) can be resolved as, where L is the first generation lepton doublet, e c is the right handed electron singlet, and H is the Higgs doublet. After electroweak symmetry breaking (EWSB), the Higgs doublet acquires a vacuum expectation value where v = 246 GeV, and the operators in Eqs. (5) and (6) reduce to the MDM and EDM dipole interactions in Eqs. (1) and (2), respectively.
Since these interactions are higher dimension operators, the VDM freeze-in abundance will depend on the reheat temperature of the universe after inflation, T RH . If T RH > 160 GeV then the universe is initially in the unbroken electroweak phase [53], and the freeze-in abundance accumulates through the interactions in Eqs. (5) and (6). As shown in Ref. [17], this yields, where VDM is produced via eh → eV and e + e − → hV reactions, and h is the Higgs field. Note that Eq. (7) holds for both the MDM and EDM models. If, on the other hand, T RH < 160 GeV, the radiation era begins in the broken electroweak phase, and VDM freeze-in proceeds through the interactions in Eqs. (1) and (2), for which the abundance satisfies [17], Here, the main freeze-in reactions are now eγ → eV and e + e − → γV . The key difference between Eq. (8) and Eq. (7) is the T RH dependence; the VDM abundance is more strongly dependent on T RH in the unbroken electroweak phase. Note that in both Eqs. (7) and (8), we only use one of the MDM/EDM operators to calculate the abundance. Additionally, similar to Eq. (7), this expression holds for both the MDM and EDM models. In Figs. 1 and 2 we show gray dashed contours corresponding to the d M/E necessary to generate the observed freeze-in abundance for various choices of T RH and m V > ∼ keV. For smaller masses, freeze-in produces warm VDM in conflict with the observed matter power spectrum on small scales [17], and therefore the abundance curves do not extend below m V ∼ keV.
As discussed earlier, for smaller masses the abundance must be set by other mechanisms. However, for this to be true, it must be the case that the VDM is not thermalized via pair annihilation and Compton-like scattering processes. Conservatively assuming T RH = 1 MeV leads to a cosmological consistency condition of d M/E < ∼ 10 −6 GeV −1 , so we do not plot above above this value in Figs. 1 and 2.

III. ABSORPTION RATE CALCULATION
In this section we calculate VDM absorption rates through the MDM and EDM interactions in Eqs. (1) and (2), respectively. Bosonic DM absorption rates in atomic targets have previously been calculated for the axion-like particles (ALP) and kinetically mixed dark photons [14,43,44]. However, since we are studying different DM interactions, the corresponding absorption rates need to be derived from first principles.
In our calculation, we follow the approach in Refs. [45,46], which extracted general absorption rates in terms of bosonic self-energies in the non-relativistic (NR) limit of the interaction Lagrangian. 1 The advantages of this approach are that it applies to any DM model or target electronic structure, and automatically incorporates any in-medium effects (although these are mainly important for crystal targets with O(eV) band gaps).
If the dark photon, V , does not mix with the photon, A, the optical theorem tells us that the absorption rate of the i th polarization of V , Γ i V , is given by, where Π i V V is the self-energy of the i th polarization. However, if V mixes with A (e.g., through a loop of electrons) then V and A are no longer the eigenstates of the theory, and the true eigenstates, V ′ and A ′ , are those that diagonalize the 2 × 2 self-energy matrix between V and A [45,46,54]. In this case, the VDM absorption rate is related to the imaginary part of the V ′ , V ′ self-energy, where we have assumed that V, A are perturbatively coupled. The absorption rate per unit exposure, averaged over the incoming DM polarizations, now becomes, where ρ V = 0.4 GeV/ cm 3 is the DM mass density, and ρ T is the target mass density. Assuming that the selfenergies are independent of polarization, shown explicitly in App. B for the isotropic targets of interest here, Eq. (12) becomes, and computing the absorption rate becomes a problem of evaluating the relevant self-energies, To calculate the self-energies we use an NR effective field theory (EFT) of electrons appropriate for energy and momentum transfers below the electron mass. This involves taking the NR limit of the QED Lagrangian, supplemented with the interaction terms from Eqs. (1) and (2). This procedure is a tedious, but straightforward, exercise performed in Ref. [45], and which we detail in App. A. The full expressions of the NR limit of the MDM and EDM Lagrangians, to O(m −2 e ), can be found in Eqs. (A22) and (A23), respectively.
While the full NR Lagrangians are relatively complicated, different approximations only leave a few important terms. First, we assume that the target has no spin ordering (i.e. there is no net electronic spin polarization) and that the electronic states are spin degenerate. This allows all of the self-energies to be written in terms of spin-independent matrix elements. Second, in typical targets the electron velocity, v e ∼ Zα > ∼ 10 −2 , is greater than the halo DM velocity of ∼ 10 −3 . This allows us to neglect many terms which are proportional to the DM momentum, q.
Explicitly, the terms which will give the dominant contribution to the absorption rates, via Π V V , are, for the MDM and EDM models, respectively, where σ are the Pauli matrices, and k = −i∇ is the electron momentum. These lead to the self-energies, where i, j are summed indices, ω is the energy flowing through the self-energy diagram, and v i are the components of v ≡ k/m e . TheΠ O1,O2 are then computed in terms of the target electronic structure [45][46][47], where V is the target volume, |I⟩, |F ⟩ are the initial and final electronic states, respectively, ω I , ω F are the initial and final state energies, respectively, T O ≡ ⟨F |O|I⟩ is the transition matrix element for Hermitian operator O, and δ is the width of the electron resonance. 2 This expression will take different forms in crystal and atomic targets, since the electronic states, |I⟩, |F ⟩, differ between them. Explicit forms for the transition matrix elements that defineΠ O1,O2 for crystal targets have been discussed in detail in Refs. [45][46][47]. In App. B we derive the results for atomic targets. Note that this definition ofΠ is a slight generalization from previous works to account for non-unit normalized final states, ⟨F |F ⟩ ̸ = 1. This is useful when working with a continuum of final states, as appropriate for atomic targets.
Additionally, one can show that starting from the complete Lagrangians in App. A, at leading order the V, A mixing self-energies are only non-zero in the MDM 2 Strictly speaking, there should be an additional phase factor in the definition of the transition matrix element: [45][46][47]. However, for absorption kinematics, q ≪ ω, the phase factor is generally negligible, except for O = 1 due to state orthonormality. Therefore when discussing T 1 we keep the phase factor. model, and are related to the photon self-energy, so the MDM model generates m V /m e suppressed mixing effects, while there are no mixing effects in the EDM model. With all of the self-energies computed, we can now compute the absorption rates for the MDM and EDM models by substituting Eqs. (16), (17) and (19) in to Eq. (13). While the expressions in terms of the selfenergies are identical between different targets, for the MDM model the rate can be written in terms of the photon self-energy since both the V, A mixing term and the imaginary part of the V, V self-energy are related to Π AA , Therefore, while the MDM model absorption rate can be written in terms of the photon self-energy, that does not necessarily imply that it is related to the photon absorption rate. The reason for this is kinematics: when a photon with energy ω is absorbed, the momentum absorbed by the target is q = ω. Therefore the photon absorption rate is determined by Π AA (q = ωq, ω), whereq is the direction of the incoming photon. However, when VDM with energy ω ≈ m V is absorbed, the momentum absorbed by the target is much smaller, where v V ∼ 10 −3 is the VDM velocity. Therefore, only when, can the VDM absorption rate be related to the photon absorption rate, at incoming photon energies of ω = m V . To understand when this is a good approximation, it is important to know that Π AA is a function of the matrix element T 1 . If the dipole approximation is valid, then Eq. (21) is also approximately true, since the Π AA (q → 0, m V ) depends on the approximated, right side of Eq. (22), and Π AA (q = m Vq , m V ) depends on the left side of Eq. (22) evaluated at q = m Vq . The dipole approximation is valid when qx ≪ 1, where x is a typical distance scale. For atomic targets, typical x values are a 0 /Z, where a 0 is the Bohr radius, and Z is the nuclear charge. Therefore for Eq. (21) to be valid, For the atomic targets of interest here, Xe and Ar, this implies that the absorption rate of VDM with mass m V , can only be related to the photon absorption rate, at ω = m V , for m V < ∼ 100 keV. Moreover, the most accurate VDM absorption rate calculation for m V > ∼ 100 keV would be a first principles calculation done in the dipole approximation, as opposed to rescaling the photon absorption rate.
While this is an important conceptual point, it is also at the boundary of interesting VDM parameter space, since indirect detection constraints from V → 3γ become important near m V ∼ MeV. Therefore, for the most interesting VDM masses, it is appropriate to relate the MDM model absorption rate to the photon absorption rate.
For crystal targets, the natural way to express this is to write Π AA is terms of the dielectric function, ϵ(ω) [45], while for atomic targets it is more natural to use the photoelectric cross section, σ pe . These are related by, where n T is the target number density. Using these relations the absorption rates are given by, where m T is the mass of the target atom. Therefore we see that while the MDM absorption rate can be related to photon absorption in a target, this is not true for the EDM model; a similar result was found for the scalar DM absorption model discussed in Ref. [45]. Therefore to make projections for the direct detection constraints on the EDM model, the absorption rate must be computed from first principles.

IV. DIRECT DETECTION CONSTRAINTS
From the discussion in Sec. III, we know that while the absorption rate of the MDM model can be related to photon absorption, via the dielectric in crystal targets or the photoelectric cross section in atomic targets, the EDM absorption rate must be computed from first principles. Since the first principles calculation relies on an assumption about the initial and final electronic states in the target, we begin by discussing the electronic configurations assumed for the crystal Si, Ge, and atomic Xe, and Ar targets used here.
For Si and Ge targets we use the publicly available electronic configuration from Ref. [55], which has been used previously [47] to compute the absorption rate for scalar, axion-like particle, and kinetically mixed dark photon models. Detailed information about the electronic configuration can be found here [56], and a longer discussion regarding modelling of the electronic configuration in this way can be found in Refs. [47,49]. The configurations use three different methods to approximate the electronic states. The deeply bound, "core" (all orbitals inclusively below 2p in Si and 3d in Ge) states are assumed to be solutions to the Hamiltonian of an isolated atom, which are computed using the RHF method [57]. These states are expanded in an STO basis, and tabulated values of the coefficients can be found in Ref. [58]. The states closer to the Fermi surface, including four valence bands below and 60 (82) conduction bands above in Si (Ge), are computed with density functional theory (DFT) methods, expanded in a Bloch basis with an E cut = 2 keV, and are all-electron reconstructed. Lastly, the highest energy states, with energies between 60 eV and 1 keV above the Fermi surface are treated as free plane waves. Lastly, following the treatment in Refs. [45,47], we model the electron width as δ = 0.2 + 0.1 ω.
Relative to the electronic states in crystal targets, those in atomic targets are much simpler. This is because all electrons in the target are tightly bound to an individual atom, and therefore the electronic states can be determined in isolation of the other atoms in the target. Similar to the deeply bound, "core" electron states in crystal targets, we use the results of an RHF calculation [58] for the initial electronic states. For the final states, we use the continuum solutions to the Hamiltonian with a V = −Z/r potential, where Z is the nuclear charge, sometimes known as "Coloumb wave functions" [59][60][61][62], with Z set by the binding energy of the initial electronic state. 3 This approximation for the initial and final electronic states has been used in previous studies of DM-electron interactions in Xe and Ar targets [59,63,64]. More details about the initial and final electronic states, and the conventions used here, can be found in App. C.
To compute the absorption rates in Si, Ge, Xe, and Ar targets we use EXCEED-DM [47][48][49]. While EXCEED-DM has been used extensively with the electronic configurations of Si and Ge, it did not previously support atomic targets. We added support for this class of targets, and have made these updates publicly available in a new version. Additionally, the MDM and EDM model absorption rate calculation in Si and Ge targets has also been added.
Before discussing the constraints on the MDM and EDM models, it is important, when possible, to verify the electronic configurations being used against measured photon absorption data. This has been done previously for Si and Ge [45], and therefore we focus on the Xe and Ar calculations. Using Eq. (26) we can compute the photoelectric cross section, σ pe , by rescaling the dark photon absorption rate in the MDM model. In Fig. 3 we compare the photoelectric cross sections computed with EXCEED-DM to a variety of other calculations and measurements [19,20,43,61]. Overall we find good agreement in both Xe and Ar targets, with O(1) discrepancies for ω < ∼ keV, and ω > ∼ 100 keV. The slight disagreement at low energies is somewhat expected due to our "isolated" atom approximation, which completely neglects  Fig. 1, is shown in solid red. Experimental measurements from Ref. [20] are shown in dashed green. Other experimental measurements, used in Ref. [43], to place constraints on other DM models with Xe targets, are shown in dashed purple. Photoelectric cross sections from the NIST database, computed with the XCOM program [19], are shown in dashed blue. Lastly, we show the photoelectric cross section for the K shell of Ar (dashed orange), computed under the dipole approximation in Ref. [61], to further illustrate the discrepancy between the photoelectric cross section computed with photon versus dark photon kinematics, as discussed in Sec. III.
the target environment when solving for the electronic wave functions. A more sophisticated approach (e.g., using the DFT formalism) for the electronic states involved in low ω absorption would likely reduce this discrepancy.
At high energies, the difference between the XCOM and EXCEED-DM calculations is due to the kinematic difference between photon and dark photon absorption discussed in Sec. III. This high energy discrepancy can also be understood in the context of the standard dipole approximation in Eq. (22). For photon absorption, when qx = ωx ≫ 1, or equivalently when ω > ∼ Z/a 0 , the dipole approximation breaks down, and one must use the exponential form of the operator in the transition matrix element. The XCOM calculation uses the exponential form, whereas the dashed orange curve from Ref. [61] uses the dipole approximation, which underestimates the photon absorption rate at high ω. Ref. [61] also computes σ pe in the exponential form and reaches same conclusion: the dipole approximation underestimates the photon absorption rate for large q = ω. However, while the dipole approximation is not appropriate for photon absorption at high ω, it is appropriate for dark photon absorption, since q ≪ ω, and therefore qx ≪ 1 is still valid even when ω ∼ MeV, contrary to photon absorption.
With verification that our electronic configurations reproduce other observables, we now discuss the constraints for the MDM model, shown in Fig. 1, as well as the EDM model, shown in Fig. 2. For both the MDM and EDM models, the lowest probeable DM mass is set by the minimum energy difference between initial and final states, since ω ≈ m V . For crystal targets this is the band gap, which is about 1.11 eV in Si, and 0.67 eV in Ge. In atomic targets this is the ionization energy (or negative of the binding energy) of the least bound electron, i.e., the 5p electron in Xe, with E I ≈ −12 eV and the 3p electron in Ar, with E I ≈ −16 eV. 4 The high DM mass cutoff in Si and Ge at m V ∼ 1 keV is due to the fact that the electronic configuration only includes final states with final energies up to a keV. That is, the cutoff is just an analysis cutoff, not a physical one. However, for m V > ∼ keV future iterations of the XENON experiments are expected to dominate the bounds due to their large exposure. Therefore it is unlikely that Si and Ge target projections will be important above for m V > ∼ keV. As discussed in Sec. II, for m V > ∼ keV, the cosmological abundance can be set by the freeze-in mechanism. This lower bound is set by constraints on structure for-mation [65][66][67], and is shown as a dashed gray vertical line in Figs. 1, 2. If the DM was lighter than this it would be too hot, and suppress structure formation on small scales. For masses larger than this, the abundance is then set by the reheat temperature, T RH , via Eqs. (7) and (8). As shown in Fig. 1, we find that current direct detection constraints, mainly XENONnT [21], set a lower bound on the viable reheat temperature for the MDM model, T RH > ∼ 100 GeV. The constraint increases to nearly T RH > ∼ 10 3 GeV for m V ∼ 1 keV. This nearly closes the previously open parameter space on MDM VDM production via sub-electroweak scale reheat temperatures. Roughly 100× greater exposure will be needed to make the same statement about T RH ∼ 10 3 GeV, due to the stronger scaling of the abundance with T RH above the electroweak scale, as shown in Eq. (7). Furthermore, low threshold analyses of atomic targets, along with even lower threshold constraints from crystal target experiments such as CDEX [28], DAMIC [29][30][31][32][33], EDEL-WEISS [34][35][36], SENSEI [37][38][39], and SuperCDMS [40][41][42] will constrain models producing low mass MDM coupled DM via the other mechanisms discussed in Sec. II.
Similar conclusions, to that of the MDM model, hold for the EDM model shown in Fig. 2; although the constraints are slightly weaker. This can be understood from Eqs. (26), (27): the EDM absorption rate is suppressed relative to the MDM scenario by a factor of v 2 e , and therefore constraints are weaker by, roughly, a factor of Zα. Since the EDM absorption rate is not related to the photon absorption rate, there are no official direct detection constraints. To get an estimate of where these would lie, we assume that the projections, solid lines in Fig. 2, can be rescaled by the same factor that would bring the solid lines in agreement with the edge of the shaded regions in Fig. 1. Essentially, we are assuming the same "detection efficiency" in both the EDM and MDM models. These rescaled, expected limits are shown as shaded regions with dashed edges in Fig. 2.

V. CONCLUSION
In this paper we have studied the absorption rate of vector dark matter particles that couple to electrons preferentially through electric and magnetic dipole moment operators. We use these results to place new limits on these scenarios from a variety of direct detection searches utilizing both crystal and atomic targets, and make projections for future searches, which are poised to improve experimental sensitivity to these interactions by several orders of magnitude across the eV-MeV mass range. This parameter space is particularly interesting because it covers masses and dipole couplings that can yield predictive cosmological freeze-in production through the same operator responsible for absorption reactions; freeze-in through kinetic mixing is nearly fully excluded for all choices of particle mass.
While it has been known for some time that the ab-sorption of the kinetically mixed dark photon, and ALPs, can be related to the photon absorption rate [14,43,44], this relation does not hold in general. Indeed, in Ref. [45] it was shown that for scalar DM this relationship does not exist. Following the ideas presented in Ref. [45], albeit with a different derivation discussed in App. A, in Sec. III we derive the NR limit of UV Lagrangians in Eqs. (1), (2) and use these to compute the self-energies and absorption rate. This procedure is general to any DM model, and therefore should be useful for future studies of different DM models. We find that while the absorption rate of the MDM model can be related to the photon absorption rate, this is not true for the EDM model, which therefore must be computed from first principles.
To compute the absorption rate of the EDM model, we modified EXCEED-DM [47][48][49], the program previously used in the first principles study of DM absorption on crystal targets [45][46][47]. We implemented two main improvements: first, we added support for absorption calculations involving atomic targets, i.e., transitions between bound and continuum states using the standard approximations for these electronic states [59,63]. Second, we added the EDM and MDM absorption rate calculation for crystal targets, e.g., Si and Ge (in addition to the atomic targets previously mentioned). Updates to the program are publicly available.
The results of these new calculations were discussed in Sec. IV, and shown in Figs. 1 and 2. We began by verifying the first principles calculation against other photoelectric cross section measurements and calculations, in Xe and Ar targets, in Fig. 3. We find good agreement, up to O(1) factors, with the largest discrepancies below ∼ 100 eV, and above ∼ 100 keV. At low energies, the discrepancy is likely due to a too simplified treatment of the electronic states. Future work using more advanced methods, e.g., density functional theory, should decrease the disagreement and also lead to more accurate scattering rate calculations [59,63]. At high energies, the difference is due to using the dipole approximation in the transition matrix element. However, as discussed in Sec. III, the dipole approximation is appropriate in the context of DM absorption, since kinematically the dark photon is depositing much less momentum, relative to a photon, for a given energy deposition. Therefore, strictly speaking, in this region even the kinetically mixed dark photon absorption rate cannot be related to the measured photon absorption rate. While this is a theoretically interesting point, the difference ends up being marginal; and, moreover, the DM masses for which this is important are ruled out by indirect detection, as seen in Figs. 1 and 2.
VDM coupling to electrons via MDM and EDM operators are simple, motivated extensions to the SM which can account for the DM abundance, and therefore must be searched for in every possible avenue. In this paper, we have shown how to compute the direct detection rate for these models, and aim to include them in future official analyses alongside ALP and kinetically mixed dark photon constraints. Even without considering MDM and EDM couplings to other SM fermions, there is interesting physics beyond the scope of this paper yet to be explored. Searching for VDM masses below m V ∼ O(eV) will require utilizing more novel excitations in low threshold experiments, such as phonons [68][69][70][71][72][73][74][75][76][77][78][79][80][81][82] and magnons [79,[82][83][84][85][86], or electronic excitations in more novel targets, such as small band gap crystals [46,[87][88][89][90] and superconductors [91][92][93]. Additionally, a detailed study of stellar cooling constraints is important, and may place stronger constraints for m V < ∼ 10 keV than the thermalization re-quirement discussed in Sec. II.

ACKNOWLEDGMENTS
We thank Bogdan Dobrescu, Simon Knapen, Duncan Rocha, and Anastasia Sokolenko for helpful conversations. This work is supported by the Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The computations presented here were conducted on the Wilson-Institutional Cluster at Fermilab.

Appendix A: Non-Relativistic Lagrangians
The initial step in any calculation of DM induced electronic excitation rates is to reduce the UV Lagrangian, written in terms of the four component Dirac field, Ψ, to the non-relativistic (NR) Lagrangian, written in terms of a light, two component field, ψ, which solves a charged particle Schrödinger equation. In general, this is a non-trivial problem for electrons in background electromagnetic potentials. To leading order in the DM-electron coupling, this can be done in two separate steps. The first step is to find the Ψ → ψ map which reduces to the NRQED Lagrangian starting from, where D µ = ∂ µ + ieA µ . The second step is to apply the same field transformation, Ψ → ψ, on a UV DM-electron interaction vertex, e.g.,ΨO UV Ψ, which creates a UV to NR map to one written in terms of the light field, ψ, where O UV is a 4 × 4 matrix, and O NR is a 2 × 2 matrix. In Sec. A 1 we derive the NRQED Lagrangian to O(m −2 e ). That is, we take the low energy, momentum limit of Eq. (A1) to find L NR QED , whose leading terms give the Schrödinger equation of a charged particle. While the NRQED Lagrangian has clearly been known for a long time, we rederive it starting from Eq. (A1) since we are using a different method to take the UV to NR limit than was done in previous calculations [45]. Therefore a complete derivation is useful in comparing to previous results, as well as a pedagogical introduction to NRQED. Then, in Sec. A 2, we apply the field transformation on the MDM and EDM interactions of interest in Eqs. (1) and (2), respectively. The final result of these calculations is summarized in Table I.

NRQED Lagrangian via Foldy-Wouthuysen Transformation
We begin with the NRQED step, which can be done in many ways [94]. One detailed recipe for doing this in the context of DM absorption on electrons is given in Ref. [45]. Here we use an alternative formulation, known in other contexts as a Foldy-Wouthuysen (FW) transformation [95][96][97][98][99]. For this problem, and for reasons that will become clear shortly, it is better to work in the Dirac basis, where, The main problem that we are trying to solve is that the QED Lagrangian in Eq. (A1) mixes the two, two-component fields inside Ψ, which makes solving the system for each two-component field difficult. To see this define, where ψ, ψ h are two component fields, and substitute Ψ back in to the QED Lagrangian, Only ψ h has a mass term, and therefore it is referred to as the "heavy" component, and ψ is the "light" component.
We also see that ψ and ψ h are coupled due to the presence of σ i D i . One approach to decouple the terms is to integrate out the heavy field, ψ h . This is the approach taken in Ref. [45]. The idea here, and that of the FW procedure, is to perform consecutive field redefinitions, at each order in 1/m e , where U i are some operators acting on Ψ, to remove the mixing between ψ and ψ h , defined (in the Dirac basis) as the upper and lower components, respectively, of Ψ on the right hand side of Eq. (A6). To do this efficiently it is important to identify the terms which mix ψ and ψ h , or equivalently, the upper and lower components of Ψ.
In the context of FW transformations, the operators which mix ψ and ψ h are known as odd operators, and diagonal operators are even. They are defined by their (anti) commutation relations with γ 0 . Specifically, an odd operator, O, satisfies, {O, γ 0 } = 0, and an even operator, E, satisfies [E, γ 0 ] = 0. The goal of the FW transformation is then to remove all odd operators from the Lagrangian at each order in 1/m e , thereby decoupling ψ and ψ h at any given order.
Specifically, the recipe is to find n Hermitian operators, {X 0 , . . . , X n−1 }, such that the field redefinition, removes all the odd operators to O(m −n e ). To expand the QED Lagrangian, Eq. (A1), to O(m −2 e ) we need to find X 0 , X 1 . One can show that the X 0 , X 1 which do this are, where [D µ , D ν ] = ieF µν . With X 0 , X 1 in hand the NRQED Lagrangian of ψ to O(m −2 e ) can be derived by substituting the X 0 , X 1 in Eq. (A8) to Ψ in Eq. (A7), the Ψ in to the QED Lagrangian in Eq. (A1). While seemingly inconspicuous, the phase factor e −imet plays an important role here. To see this explicitly, note that under the FW transformation, Furthermore, we can define P − ≡ 1 2 (1 − γ 0 ) and rewrite the QED Lagrangian without the e imet phase as, and, again, X 0 , X 1 are given by Eq. (A8). Therefore, as in the derivation of Eq. (A5), this phase factor introduces an operator projecting the mass on to the lower component of Ψ. Following the derivation in Eq. (A10) further, the ψ part of the QED Lagrangian becomes, where P + = (1 + γ 0 )/2 is only non-zero in the upper left diagonal component, the trace is in 2 × 2 block diagonal space, and we have ignored the ψ h terms. The Pauli matrices are given by σ, the momentum variables (k, p, ω ′ ) are shorthand for derivatives, p µ = i∂ µ , k acts on ψ, p and ω ′ act on A, p 0 acts on A 0 , and K ′ = 2k + p.

NR Limit of MDM and EDM Interactions
With the field redefinition which diagonalizes QED Lagrangian found in Eq. (A7), we can now use this to expand the interactions in Eqs. (1) and (2). However before this we can reach a more general result: the NR limit of a UV general operator, to any order in 1/m e (still only leading order in DM-electron coupling), where the outer most brackets, [· · · ], implicitly indicate taking the 2 × 2 upper diagonal component of the 4 × 4 matrix. All X act to the right, meaning that those on the left side of O will also act on O. To O(m −2 e ) this expression simplifies to,Ψ  1) and (2), in three component notation we can expand the operators as, and therefore we see that there are four O operators whose NR limit must be extracted: Each of the terms in Eq. (A15) can be further simplified, defining Γ ≡ γ 0 O, where (∂ i V ) indicates that the derivative operator only acts on V . The final NR expansion, to O(m −2 e ), is given in Table I.
We can now substitute the results in Table I to the terms in Eqs. (A16) and (A17) to find the NR limit of the Lagrangians given in Eqs. (1) and (2),  (1) and (2), respectively. We only keep terms involving two fields (excluding the electron field, ψ).
where, similar to Sec. A 1, the momentum are shorthand for derivatives acting on different fields and ω, q act on V , p, ω ′ act on A, and k acts on ψ.
Appendix B: Self-Energy Calculations With the NR limit of the QED, MDM, and EDM Lagrangians, given in Eqs. (A12), (A22) and (A23), respectively, we can now derive the self-energies needed to compute the rate in Eq. (13). To do this in full generality, one must first compute the diagonal self-energies in the component basis, i.e., Π µν V V and Π µν AA , then find the polarization vectors, ϵ µ λ , such that they are diagaonlized, i.e., where ϕ ∈ {V, A}. The off-diagonal self-energies, Π λλ ′ V A are then simply the off-diagonal self-energies in the component basis, Π µν V A , projected on to the polarization vectors. For the targets of interest (Xe, Ar, Si, and Ge) we can make a key simplifying assumption: isotropy. Assuming that the targets are isotropic the spatial components of the self-energies become, where the repeated k index on the right hand side is summed. Under this assumption it can be shown that the polarizations which diagonalize Π µν ϕϕ are given by the standard longitudinal and transverse polarization vectors, where q ± are two orthogonal vectors satisfying q · q ± = 0 and Q 2 = ω 2 − q 2 . The Ward Identity (WI) simplifies the calculation of projecting the component basis self-energies to the polarization basis. One can show that, and therefore the longitudinal self-energy can be computed from Π 00 alone. Additionally, the transverse components are only related to Π ij . Therefore these are the only self-energies we need to compute. In Secs. B 1, B 2 and B 3 we derive the leading order contributions to the self-energies in the component basis, and then project them on to the polarization basis, i.e., the T, L components, via the inverse of Eq. (B1). Finding the leading order contribution is a relatively straightforward exercise in power counting with respect to the dimensionless variables, v e = k/m e , v V = q/m V and m V /m e , where v e is the electron velocity, q is the dark photon momentum, and m V is the dark photon mass. The main subtlety, discussed in detail in Ref. [45], is due to state orthonormality which reduces the order of operators. For example, while an interaction of the form, A 0 ψ † ψ is naively O(1), the reduced self-energy,Π, Eq. (18), depends on and therefore this is a suppressed operator. At the end of each section we find that Π T ≈ Π L (these are the self-energies projected on to the transverse and longitudinal polarization vectors), and therefore Π λ is approximately independent of λ, justifying the simplification in Sec. III. Additionally, we find that the mixing between V, A in the EDM model is negligible, and therefore the section deriving this is absent.
The self-energies resulting from these calculations will be written in terms of some "reduced" self-energiesΠ O1,O2 , which are only dependent on the target electronic structure. Ref. [45,47] details how these are computed for crystal targets, and in Sec. B 4 we derive the formula for atomic targets.

QED
Starting from the NR limit of the QED Lagrangian in Eq. (A12), one can derive the photon self-energies, where the single O reduced self-energy,Π O , is given by, where I runs over filled electronic states. Note that this result was also derived in Ref. [45]. TheΠ 1 term can be reduced using a WI, ωΠ 0µ − q i Π iµ = 0, to, Under our isotropic target approximation we can replaceΠ k i ,k j → 1 3Π k i ,k i δ ij , thereby reducing the self-energies to, Lastly, projecting on to the polarizations gives, Note that explicit expressions forΠ O1,O2 are computed in Refs. [45][46][47] for crystal targets and in Appendix C below for atomic targets.

Magnetic Dipole Moment
For the MDM model two types of self-energies need to be computed: the dark photon self-energy with itself, Π V V , and the mixing of the dark photon with the photon via Π V A . For readability we split these two calculations in to the subsections below.

a. Dark Photon Self-Energy
Starting from the MDM interaction in Eq. (A22) the leading order contribution to the self-energies in the component basis are, The dominant term in Eq. (A22), when computing, e.g., Π ij V V , can be easily extracted. Simply take the q → 0 limit of terms involving Vψ † ψ and there is only a single remaining term. To further reduce these expressions we must make an additional approximation relative to the photon self-energy calculations. We assume that the target is not spin-ordered, and therefore the sums over initial and final states can be split in to I → i s , where s indexes the spin states. This allows the Pauli matrices to be traced over in the reduced self-energy expressions, e.g., in Eq. (B11), Note that our convention for non spin-ordered targets is to absorb the factor of two, from the Pauli spin matrix trace, in toΠ. With this substitution the equations in Eq. (B11) become, Using the isotropic approximation, analogous to Sec. B 1, these self-energies simplify to, Lastly, we project these into the polarization basis to reach our final results, b. Dark Photon-Photon Mixed Self-Energy The dark photon-photon mixed self-energy is computed from the MDM and QED Lagrangians in Eqs. (A22) and (A23), respectively. Note that while we will only compute Π V A , at this order Π V A = Π AV . The leading order mixed self-energies are given by, where, similar to the previous section, we have assumed that the target is not spin ordered, i.e.,Π σ = 0. This assumption also allows us to trace out the σ matrices, As done in Sec. B 1 we now replaceΠ 1 via a WI, and utilize our isotropic assumption to get to, Lastly, projecting on to components gives, which are related to the photon self-energies given in Sec. B 1.

Electric Dipole Moment
The dark photon self-energy in the EDM model is derived from the Lagrangian in Eq. (A23). In the component basis they are given by, Again assuming no spin ordering and target isotropy we can replace, leading to, Lastly, projecting on to components gives,

Atomic Target Self-Energy
In this subsection, we simplify the reduced self-energy, also given in Eq. (18), for the electronic states in an atomic target. More detailed specifics about the electronic states are discussed in Sec. C. While our focus will be on an atomic target with a single atomic species, the analysis here easily generalizes to multiple atomic species. In such a target there are N = n × V target atoms, where n is the number density and V is the target volume. Each atom hosts bound electronic states which are indexed by the quantum numbers n, ℓ, m. Assuming these states are spin degenerate, the sum over initial electron states in Eq. (B24) becomes, The discrete nature of this sum makes this substitution intuitive. On the other hand, the final states form a continuum indexed by k, ℓ ′ , m ′ where ℓ ′ , m ′ are angular quantum numbers, and k is a continuous momentum. The sum over k then becomes an integral, where δ(0) is a normalization coefficient which will drop out of the final formulas. The simplest way to take care of this is to define the initial (discrete, bound) and final (continuum) states with discrete and continuous normalizations, Substituting these sums and state normalizations leads to a reduced self-energy of, where δ is the electron width. Further simplifications can be made to isolate the imaginary part of this self-energy needed to compute the rate in Eqs. (26) and (27). In the limit of zero electron width, δ → 0, the Green's function, G, simplifies, and the imaginary part of the reduced self-energy becomes a simple sum. The transition matrix elements, T O , are derived in App. C.

Appendix C: Atomic Absorption Transition Matrix Elements
To compute the self-energies derived in Sec. B, the transition matrix elements, T O ≡ ⟨F |O|I⟩ must be computed. The details of this calculation for crystal targets, e.g., Si and Ge, has been discussed in Refs. [45,47], and similar calculations for transitions between atomic bound and continuum states have been performed in the scattering limit [59,63,100]. In this appendix we derive the transition matrix elements between the bound and continuum electronic states in an isolated atom (appropriate for Xe and Ar targets), in the absorption limit.
We begin by defining the initial, |I⟩, and final, |F ⟩, quantum states (mainly for comparison with other conventions). The initial, bound states are labeled by the standard quantum numbers, n, ℓ, m, and satisfy, Note that the states are dimensionless. The position space representation of these states is, where ⟨x|x⟩ = 1, and ψ nℓm has dimension eV 3/2 . We assume that the position space representation of these states can be further expanded as, where R nℓ has dimension eV 3/2 , and Y ℓm are the spherical harmonics, normalized to, In order to satisfy the orthonormality relationship in Eq. (C1), the R nℓ must satisfy, dr r 2 R * n ′ ℓ R nℓ = δ nn ′ .
A useful basis to expand R nℓ in is the "Slater type orbital" basis, where a 0 is the Bohr radius, C n,ℓ,j , n ℓ,j , Z ℓj are coefficients found by solving the isolated atom Hamiltonian. We use the coefficients tabulated in Ref. [58] to compute the results in the main text. The final states are taken to be the Coloumb wave function solutions to a −Z/r potential [59][60][61][62][63][64]. These states are labelled by k, ℓ, m, where k is a continuous index. Different conventions are reasonable for orthonormalizing these states; here we choose, The position space representation can be decomposed in a way analogous to the initial states, where R kℓ have dimension eV. Note that the R kℓ here differ from those defined in Ref. [59] by a factor of k/2π. Specifically, R kℓ = (k/2π) √ VR kℓ , whereR kℓ are defined in Ref. [59]. The R kℓ in Eq. (C9) satisfy, and are explicitly given by [60], R kℓ (r; Z F ) = 2 r C ℓ ρ ℓ+1 e −iρ 1 F 1 (ℓ + 1 − iη, 2ℓ + 2, 2iρ) (C11) where we have defined C l = 2 ℓ (2ℓ + 1)! e −πη/2 |Γ(ℓ + 1 + iη)| , ρ = kr , η = − Z a 0 k , and 1 F 1 (a, b, c) is the confluent hypergeometric function of the first kind. One common approximation for the Z F parameter [40,101,102] is to relate it to the binding energy of the state it was absorbed from, i.e., where n, ω I are properties of the initial states used when calculating the transition matrix elements, T . With these conventions, the transition matrix elements become In the next two subsections we derive explicit forms for the O needed to derive the atomic absorption rate for the magnetic dipole and electric dipole models studied in the main text.
To compute the second derivative of the initial state wave functions we simply need to use the identity in Eq. (C15) twice, where it is understood that the c, d superscripted with i are evaluated at (ℓ, m, k, q), and the those superscripted with j are evaluated at (ℓ + k, m + q, k ′ , q ′ ). Therefore the transition matrix element is given by, where I α,β IF is defined in Eq. (C36), and it is understood that the c, d superscripted with i are evaluated at (ℓ, m, k, q), and the those superscripted with j are evaluated at (ℓ+k, m+q, ∆ℓ−k, ∆m−q). The selection rules for this transition matrix element are |∆ℓ| ≤ 2 and |∆m| ≤ 2.