Accurate electron-recoil ionization factors for dark matter direct detection in xenon, krypton and argon

While most scintillation-based dark matter experiments search for Weakly Interacting Massive Particles (WIMPs), a sub-GeV WIMP-like particle may also be detectable in these experiments. While dark matter of this type and scale would not leave appreciable nuclear recoil signals, it may instead induce ionization of atomic electrons. Accurate modelling of the atomic wavefunctions is key to investigating this possibility, with incorrect treatment leading to a large suppression in the atomic excitation factors. We have calculated these atomic factors for argon, krypton and xenon and present the tabulated results for use with a range of dark matter models. This is made possible by the separability of the atomic and dark matter form factor, allowing the atomic factors to be calculated for general couplings; we include tables for vector, scalar, pseudovector, and pseudoscalar electron couplings. Additionally, we calculate electron impact total ionization cross sections for xenon using the tabulated results as a test of accuracy. Lastly, we provide an example calculation of the event rate for dark matter scattering on electrons in XENON1T and show that these calculations depend heavily on how the low-energy response of the detector is modelled.


I. INTRODUCTION
As the astrophysical evidence for the existence of dark matter (DM) has strengthened, research into its identity has started to seep into many fields of physics.However, even after years of dedicated experiments, the nature of DM remains a mystery, with no confirmed detection to date [1].The widely held explanation is that DM is an undiscovered particle that interacts primarily via gravity, with a very weak coupling to ordinary matter that could be exploited as a detection route [2].Of the particle candidates, the Weakly Interacting Massive Particle (WIMP) is the most sought, with numerous experiments designed to search for WIMPs with GeV mass-scales and above (see, e.g., Refs.[3][4][5][6]).
Experiments utilising dual-phase time projection chambers (TPCs) are of particular interest.In these detectors, the bulk of the scintillating material is in a liquid phase with an applied electric field (sometimes called the 'drift' field), while the remaining material is in its gaseous phase above it, with a stronger electric field [21].Due to this set up, results from these types of detectors can come in the form of S1 and S2 signals.
The prompt scintillation signal, S1, occurs in the liquid phase when a collision between an incoming particle and an atom of the scintillator causes a release of photons.Due to the electric field in this section, if the collision instead results in atomic electrons being ionized, those electrons will drift upward, toward the gaseous phase (for more details, see, e.g., Refs.[5,6,21,22]).For the case of sub-GeV DM interacting with a scintillator, most research focuses on looking at the S2 signal (see, e.g., Refs.[23,24]).However, DM-electron interactions may have a higher chance of creating a detectable S1 signal than previous research suggests due to enhancements in the event rate [10].
To explore the possibility of atomic ionization we need to calculate atomic ionization factors, the details of which are discussed in Sec.III.However, the calculations present many difficulties.Depending on the details of the experiment, accurate atomic ionization factors are often required across many orders of magnitude of energy deposition (∼ eV to keV) and momentum transfer (∼ keV to MeV).Futhermore, as inaccurate description of the atomic wavefunctions can lead to errors of up to several orders of magnitude in the calculations, this prevents many common and convenient approximations from being used, as previously discussed in Ref. [10].
At high values of energy and momentum transfer, relativistic effects become crucial to the calculations [9,25].These effects can even dominate the calculations, as the parts of the the electron wavefunctions that are closest to the nucleus contribute the most to scattering, and so we need fully relativistic wavefunctions to accurately model the small-distance behaviour.
At moderate momentum transfer values, the smalldistance scaling of the atomic wavefunctions is again very important.This can lead to drastic errors in the calculations when the wavefunctions are approximated as hydrogen-like and scaled by the relevant factor (some-times referred to as "effective-Z" methods).
At low momentum transfer, but also arising at all scales to some degree, the form of the continuum wavefunctions is crucial to the calculations.Approximating these continuum wavefunctions as plane waves misses the significant Sommerfeld-like enhancement Refs.[8,10].The attractive potential of the nucleus means that plane waves do not have appropriate small-distance scaling in this area.For the continuum states to be unbound energy eigenstates, they must be found in a self-consistent atomic potential to ensure correct orthogonality to the bound electrons [17,26].
The need for accuracy also extends to the continuum state energy, as solving the Dirac equation in this region can be numerically unstable.Finally, at moderately low values of both energy deposition and momentum transfer, the atomic ionization factor can depend significantly on the atomic potential itself.
Our method applies the relativistic Hartree-Fock approximation, and accounts for the most important manybody effects.This approach addresses all of the above issues, allowing accurate calculation at the energy deposition and momentum transfer values relevant to DMelectron scattering.Additionally, we can test our method by calculation electron-impact ionization cross sections, which is a similar process to DM-impact ionization for low-mass WIMP-like particles.High accuracy experimental rates have been measured for xenon in the relevant impact energy regime, allowing a stringent test of the accuracy of our method.In the important of ∼ keV incident energies, our calculations agree with experiment substantially better than dedicated calculations that focused solely on electron impact rates.
We present tables of atomic excitation factors for argon, krypton, and xenon as supplementary materials [27] and on GitHub [28].These tables can be used in conjunction with DM form factors to calculate cross section and event rates, without the underestimates that arise from inaccurate atomic physics.An example of this process is provided as code, also on Github alongside the tables [28].Finally, we note that our code and technique, while developed for scattering, also applies for absorption (see, e.g., Refs.[26,29,30]), and may be beneficial in that case also.

II. THEORY
We consider inelastic scattering between a nonrelativistic DM particle and an atom.We may model the DM-electron interaction with an effective Yukawa coupling, where α χ is the DM-electron coupling strength, and µ is the inverse interaction length scale.This form of potential will result from the non-relativistic limit of a vector or scalar coupling between electrons and DM particles, in which case µ may be recognised as µ = m v c/ , where m v is the mass of the mediator particle.This potential will reduce to a Coulomb-like potential if the mediator is massless, or to a contact interaction if it is heavy.The differential cross section for an atomic electron in initial state nκm to be ionized into final state εκ m , can be written as where E is the energy deposition, v is a fixed DM velocity, q is the momentum transfer, E H ≡ m e c 2 α 2 ≈ 27.2 eV is the Hartree energy unit, which we introduce to to make K dimensionless, K is the dimensionless atomic excitation factor [10], n is the principal quantum number of the bound electron, and ε is the energy of the ionized electron, which we can write as ε = E − I nκ , with I nκ denoting the ionization energy of state nκ.The integration limits, q ± , are the allowed range of momentum transfer1 Lastly, κ is the Dirac quantum number, defined in terms of the quantum numbers for orbital angular momentum, l, and total angular momentum, j: κ = (l − j)(2j + 1).The atomic excitation factor is proportional to the chance of the transition from state nκm to εκ m occurring due to the interaction with a DM particle, and is defined in Ref. [10] for a vector electron coupling as where m is the magnetic quantum number.We stress that the final states are energy eigenstates, not momentum eigenstates.It is therefore natural to normalise on the energy scale, such that The continuum wavefunctions then have dimension: The factor E H is introduced to make K dimensionless.Care should be taken comparing ionization factors calculated in different works, which may choose different normalisation.
We then sum over all electrons to get the total atomic excitation factor, K, which we can use to reach the velocity-averaged differential cross section, expressed as where σe is the free electron cross section at a fixed momentum transfer of q = a −1 0 ≈ 3.6 keV, which we introduce following Ref.[31,32] to simplify comparisons between results (a 0 is the Bohr radius; see Ref. [10] for full expression linking σe back to α χ parameter of Eq. ( 1)).The DM speed distribution is denoted as f , which we assume to be a Maxwell-Boltzmann distribution with the Standard Halo Model assumptions (see, e.g.,Refs.[33,34]).We note that electron recoil spectrum may be particularly sensitive to details of the velocity disctributions, see, e.g., Refs.[9,10,35].Lastly, F χ is the DM form factor.
As F χ is able to be separated from the atomic excitation factor, the atomic excitation factors are largely independent of the DM model, and we only need to consider the electron coupling.For the example calculations in Sec.V, we use a DM form factor that is applicable to both vector and scalar couplings, as defined in Ref. [10] as While the effective potential that we began with can be applied to vector and scalar interactions, Eq. ( 4) is only applicable to the vector case.We can reach scalar, pseudoscalar, or pseudovector electron couplings if we replace e iq•r in the matrix element with γ 0 e iq•r , γ 0 γ 5 e iq•r , or γ 5 e iq•r , respectively.

III. CALCULATIONS A. Hartree-Fock approximation
As finding the atomic excitation factor includes a matrix element calculation with both the bound and continuum wavefunctions, we need to accurately form the electron orbitals.We use the relativistic Hartree-Fock (HF) method, which is both self-consistent and includes the electron exchange interaction.The final HF potential is then used to solve the Dirac equation at each step in a range of energy deposition values for the continuum wavefunctions.We write the HF potential as the sum of the direct and exchange potentials, where the first term corresponds to the direct potential, the second term to the exchange potential, i denotes the bound electron state with quantum numbers {n i , κ i , m i }, N c is the total number of electrons, and finally First, the Hartree-Fock equations are solved selfconsistently for the N c bound electrons.Then, the wavefunctions for the continuum electrons are found in the frozen Hartree-Fock potential by directly solving the Dirac equation, including the exchange term.(A small deviation from the frozen potential, known as the holeparticle interaction, will be discussed below.)We use the energy normalisation condition for the continuum states (5).In practice, this is achieved by solving the Dirac equation by integrating outwards to very large distances from the nucleus, where the wavefunctions take a spherical wave form analogous to the hydrogen-like case, and comparing to analytic solutions [36].

B. Calculation of atomic ionization factors
In the Dirac basis, we write the bound state orbitals as and similarly for the the continuum state orbitals where f and g are the large and small components of the Dirac wavefunction, and Ω is a two-component spherical spinor.We can express Ω as where j 1 m 1 j 2 m 2 |JM denotes a Clebsch-Gordon coefficient, Y is a spherical harmonic, χ is a spin eigenstate, and s z is the electron spin, such that the sum runs over To calculate the matrix element in Eq. ( 4), we use irreducible spherical tensors to expand the exponential term: where where j L is a spherical bessel function of the first kind and L is the multipolarity.With this, and angular reduction rules [37], we can then write the atomic excitation factor as where R is the radial integral and C is an angular coefficient.In the case of a vector electron coupling, R can be expressed as (16) If we instead consider the scalar case, where we replace e iq•r with γ 0 e iq•r , this will result in a radial integral of which can be re-written in a form that is more numerically stable at low q due to the orthogonality condition: as discussed below.For the pseudoscalar case, replacing e iq•r with γ 0 γ 5 e iq•r gives and finally, for the pseudovector case, replacing e iq•r with γ 5 e iq•r gives (20) In the cases of both vector and scalar electron couplings, C can be expressed as where [J] ≡ 2J + 1, the term in parenthesis is a Wigner 3-j symbol, and Π L ll the parity selection rule (it is unity if l + l + L is even and zero otherwise).The angular coefficient for the pseudovector and pseudoscalar cases are similar, but we replace κ with κ = −κ and l with l = |κ + 1/2| − 1/2.
Example calculations for the total atomic factors for each of these couplings can be seen in Figs. ( 1) and ( 2) as functions of energy deposition and momentum transfer, respectively.Tables of these factors are included as supplementary material.Alternatively, the tables can be found on GitHub [28], alongside example code that uses the tables to calculate cross sections as seen in Sec.IV and Sec.V.
As a stability test for the calculations, we also computed K using an approximate local potential instead of the HF potential.The local potential used was the nuclear potential plus a parametric potential, making the calculation far simpler and less prone to numerical instabilities.The local potential leads to very similar atomic excitation factors, giving a point of comparison that highlights possible numerical issues or errors in the calculation.

C. Many-body effects
The most important many-body effect (that is, deviation from the frozen-core Hartree Fock approximation) is the hole-particle interaction.Physically, this effect arises due to the deviation of the (direct) Hartree-Fock potential for the ionized electron compared to those in the core, as in Fig. 4. In practical Hartree-Fock calculations for occupied core states, the self-interaction term is included in the direct potential; this is then exactly compensated by the corresponding term in the exchange potential, e.g., by setting i = a in Eq. ( 8).However, this cancellation does not apply to an electron that has been excited out of FIG. 3. Velocity-averaged differential cross section for xenon when accounting for the hole-particle interaction (solid line), and when excluding it ("frozen" Hartre-Fock, dashed line).In this case, we have set the DM mass to be mχ = 1 GeV.The DM form factor is set to Fχ = 1, corresponding to a heavy mediator as per Eq. 7. We use a vector electron coupling with free electron cross-section σe = 10 −37 cm 2 .FIG. 4. Goldstone diagrams for the direct (left), exchange (middle), and hole-particle (right) contributions to the atomic potential.The wavy line is the Coulomb interaction, external lines are either bound atomic states or unbound continuum states, internal lines are core states (holes).The hole-particle effect arises to the deviation of the direct potential for the ejected continuum electron from that of the bound electrons.
the core.Therefore, the self-interaction term should be removed for the excited states.While this hole-particle interaction term makes a very small difference when looking at the ionization of bound electrons with high ionization energies (close to the nucleus), the impact becomes more obvious for electrons further from the nucleus, and for lower-energy scatterings.From Fig. 3, we can see that this contribution to the atomic factors carries through to the cross section calculation, resulting in a more significant discrepancy as we move to lower energies.
Beyond Hartree-Fock and hole-particle effects, the next most important many-body correction to matrix elements of the external field is the core polarisation, often referred to as the relativistic random phase approximation with exchange (RPA) [38].The lowest-order core polarisation diagrams are shown in Fig. 5. Physically, this effect arises due to the combined action of the external field and inter-electron Coulomb interaction.In the present case, it manifests in the possibility that when the dark matter particle interacts with an electron in one state, an electron in a different state may become ionized FIG. 5. Goldstone diagrams for core polarisation correction to matrix elements, showing the direct (left) and exchange (right) contributions.The diagram notation is the same as in Fig. 4, with the dotted line representing interaction with the external field (DM particle), and the internal forward lines representing virtual excited atomic states.In the all-order RPA method, each of the external field vertices are replaced with these four diagrams; this process is repeated until the matrix elements converge.
due to its Coulomb interaction with the first electron.This may be particularly important in cases where, for example, only s-states have appreciable e iqr matrix elements (due to their non-zero wavefunctions at the nucleus), but they are energetically inaccessible.In that case, an energetically accessible p, d-state may be ejected via Coulomb interaction with an s-electron that interacts with the dark matter.
In the many-body perturbation theory diagrams, there is an implied sum over intermediate states, with the backward lines denoting the core (bound, occupied) atomic states, and the forward lines representing the full spectrum of excited (unoccupied) bound and continuum states.To approximate this spectrum, we form an approximately complete, though finite, basis by diagonalising a set of B-splines over the atomic Hamiltonian (see, e.g., [39,40]).We use the dual-kinetic-balance basis as introduced in Ref. [41], which offers a high level of convergence and stability.Our technique is well tested, and produces high-accuracy results across a range of atomic systems [42].
The lowest-order calculation typically overestimates the core-polarisation contribution, and often significantly.For accurate calculations, all-orders (in the Coulomb interaction) calculations are required.In the RPA approximation, this is achieved iteratively, by replacing each external field vertex in the four RPA diagrams (Fig. 5) with the four RPA diagrams.This process is continued iteratively until convergence is reached.This is equivalent (up to first-order in external field) to including the action of the external field into the potential when solving the Hartree-Fock equations (known as time-dependent Hartree-Fock method [43]).We note that the RPA equations must be iterated separately for each value of momentum transfer q, and each multipolarity, L (see Eq. ( 13)).We thus include the all-orders RPA effects only for a subset of the parameter range to check its contribution, which we find to be small.
The lowest-order core polarisation effects has the largest impact at small energy deposition values, and gives at most a correction of a factor of two.The allorders RPA corrections significantly reduce this.After integrations, the uncertainty in Hartree-Fock calculations from excluding RPA is of order 20%, which is more than sufficient for the current purpose.
As a final consideration, we check for errors due to imperfect orthogonality.Due to numerical uncertainties, the exact orthogonality is not guaranteed between the core and continuum states in practical calculations, particularly when the hole-particle interaction is included.This can be source of significant errors, so must be checked.We can check for this by explicitly enforcing orthogonality between the continuum states and the states in the core using the standard procedure: |a → |a − |b b|a .In our calculations, this makes very little difference, due to the already good orthogonality achieved in the self-consistent Hartree-Fock procedure.However, in methods where the self-consistency of the potential cannot be guaranteed, this check is essential.
Another option for ensuring orthogonality between core and continuum states is to make a substitution in the matrix element of e iq•r → e iq•r − 1 in Eq. 4 (or, equivalently, j L → j L − 1 in Eq. 15).While this corrects errors caused by orthogonality in regions of low momentum transfer, it can shift the error to high momentum transfer.As we need accurate atomic factors across many orders of magnitude of momentum transfer for DM-electron scattering, this substitution alone cannot ensure orthogonality for this case.This is especially true when accounting for the hole-particle interaction.As a result, we opt to use the previous procedure for ensuring orthogonality, as this addresses issues across the required range.

D. Approximation of Atomic Excitation Factors
We present an informative approximate approach to presenting the ionization form factors that is valid for q 1/a 0 .When the energy deposition goes above the ionization of a bound electron, this electron becomes 'accessible'; the atomic excitation factor for that electron is also zero below this point.For an increasing energy deposition past this point, the atomic excitation factor is relatively constant while the continuum energy remains small, so long as q is reasonably large.This is because for large q, only the low-r part of the electron wavefunction may contribute.In this region, the energy of the bound or continuum electron is insignificant compared to the nuclear potential |ε| |Z/r|, and the Dirac equation is independent of energy.In this case, we may approximate K nκ (E, q) as a step function of energy, allowing it to be expressed as where Θ is the Heaviside step function, and Knκ is the atomic excitation factor that is dependent on momentum transfer at a fixed energy deposition.Thus, above the ionization energy, K nκ is equal to Knκ .This approximate method loses accuracy for very small values of momentum transfer (q 0.1 MeV).For the typical momentum transfer values for DM-electron scattering, inaccuracies found in the low-q region do not have a significant impact on cross-section calculations due to integrating over q, as seen in Eq. ( 6), but care should be taken if high accuracy is needed in this region.

IV. ELECTRON IMPACT IONIZATION
We have also calculated the total cross section for the case of an atomic electron being ionized by an incoming free electron (referred to as 'electron impact' (EI) ionization in the literature).With existing experimental data, the calculations being applied to this interaction type serves as a test of the accuracy of the code.In the case that m χ → m e , m v → 0, and α χ → α, it can be seen that the electron-impact and DM-induced ionizations are very similar atomic processes.
To compare to experimental data, we can calculate the total cross section, where E i is the incident energy of the projectile electron, and E is again the energy deposition.
As found in past literature, calculations of the total cross sections for this type of scattering are, on average, overestimated when compared to experimental val-ues (see, e.g., Ref. [47] and references therein).For testing DM-induced ionization, the relevant incident energy scale is m χ v 2 ∼ (GeV)(10 −3 ) 2 ∼ keV.Our calculations are presented on Fig. 6, along with experimental data, and existing calculations for comparison.
When the incident electron energy is low (E i 100 eV), we can see that our calculations (and those of other groups) overesimate the total ionization cross section.This divergence from experiment largely stems from our use of the Born approximation [44], which is only valid for very weak interactions.As the DM-electron interactions that we consider in this work falls into this category [25], the approximation holds, and does not carry the same importance as it does when looking at electron impact ionization.Further, the electron-impact ionization case is complicated by the presence of an exchange term in the interaction (which we did not account for), since the projectile is an electron which may be exchanged with a bound atomic electron.Since no such term exists in the DM case, the accuracy for the DM scattering case is expected to be higher.
Importantly, in the region where energies are closer to that of an DM-electron scattering event (E i 1 keV), we can see in Fig. 6 that the case where we have subtracted the hole-particle interaction gives results that are closest to experiment.This shows an improvement on the accuracy of existing results in this high-energy region and highlights the importance of subtracting the hole-particle interaction term.This is both experimental verification of our atomic wavefunctions, and also of the approximations used to derive the scattering formula.

V. EXAMPLE EVENT RATE CALCULATION
For DM-electron scattering, the differential event rate is, where n T is the number of target atoms in the detector per unit mass (reciprocal of atomic mass), m χ is the DM mass, and ρ χ is the local energy density of DM, which we take to be ∼ 0.4 GeV cm −3 [48].However, we need to consider that not all events are feasibly detectable and account for capability of the detector itself.This allows us to find observable event rates for specific energy regions and then compare to the event rates seen in experiments.
As the detection capabilities vary for each detector depending on the equipment, the set-up of the detector and the scintillating material, the model will also change per experiment.As an example of this calculation, we follow Ref. [49] and start by modelling the detector resolution of XENON1T as a Gaussian with standard deviation, where a = (0.310 ± 0.004) √ keV and b = 0.0037±0.0003.

FIG. 7.
Example event rate calculations before (dotted lines) and after (solid lines) smearing with the Gaussian and correcting for detection efficiency [49], corresponding to Eqs. ( 24) and ( 26), respectively.The spikes seen in the calculated event rate are due to deeper shells becoming accessible as the energy deposition increases, which we can also see in the cross section.A cut-off has been applied to the observable event rate at 0.5 keV to indicate the minimum energy threshold.Again, σe = 10 −37 cm 2 .
We use this Gaussian, g σ , to smear the theoretical event rate, given in Eq. ( 24).The event rate also needs to be corrected for the efficiency of the detector, which we accounted for by fitting the total efficiency as a function of energy, ε(E), as given in Fig. 2 of Ref. [49].Combining these steps together, we can express the observable differential event rate as From Fig. 7, we can see that correcting for the specific detector dampens the event rate as expected.For the lowest mass case, 0.1 GeV, the calculated event rate drops to zero at a lower energy than the observable event rate peaks at.This is due to the detector response being modelled as a Gaussian, which allows very low energy signals to 'leak' through to higher energy regions and be considered detectable.While this is plausible at higher energies [50], this introduces large uncertainties in the observable event rate at low energies, where it's not appropriate to model the response using a simple Gaussian.
From Fig. 7, we can see that the theoretical event rate peaks below 1 keV for all DM mass cases.Although the tail ends of each Gaussian will have very low magnitude in this region, the magnitude of the event rate causes the observable event rate to be high as well.This implies that modelling the detector in this way overestimates the observable event rate when looking at any scattering type that has large contributions in regions of low energy deposition.In these cases, care should be taken to determine the strength of the effect that the use of the Gaussian has on the observable event rate, and whether FIG. 8. Example rate calculations for mχ = 0.1 GeV, with a vector electron coupling and a heavy mediator.The scattering rate (dashed line) represents the case of a perfect detector, matching the theoretical dR/dE, whereas the observable event rate (solid lines) represent the rate that would be seen by a detector with varying energy resolution modelling.The case of a = 0.310 matches that of XENON1T, while a = 0.2 and a = 0.1 are used as comparison to highlight the sensitivity to any change in the Gaussian.Again, σe = 10 −37 cm 2 .
an overestimation is being introduced by events happening far below threshold.If the effect is large, a more accurate option is to directly simulate the specific detector, such as the use of NEST (Noble Element Simulation Technique) [51][52][53][54], or Obscura [55].This is particularly important to consider for scattering cases involving DM, where the nature of the particle is unknown, as an overestimation in the predicted event rate could lead to DM models being excluded preemptively.

VI. CONCLUSION
We have presented tables of atomic excitation factors as functions of momentum transfer and energy deposition for DM-electron interactions with vector, scalar, pseudovector, and pseudoscalar couplings for argon, krypton, and xenon [27,28].These take relativistic effects into account and provides an accurate depiction of the atomic physics involved.As such, they can be combined with a DM model of choice to calculate cross sections and event rates, without risking the underestimates that are common when neglecting the atomic physics.
We have tested the code for any numerical instabilities and errors by calculating the atomic factor with an approximate potential.We have also tested the code for accuracy by applying the calculations to electron impact ionization and found good agreement with experimental results.
As these atomic factors may be used to calculate cross sections and event rates, we have presented an example of an event rate calculation specific to the XENON1T experiment.While accurate atomic physics is necessary to reach reliable event rates, we have also highlighted the importance of the modelling of the detector, particularly the low-energy response.

FIG. 1 .
FIG.1.Comparison of the total atomic excitation factor as a function of energy deposition for xenon for vector (solid line), scalar (dotted line), pseudovector (dashed line), and pseudoscalar (dash-dotted line) electron couplings at a fixed momentum transfer of q = 4 MeV.

FIG. 6 .
FIG.6.Calculations of total electron impact ionization cross sections for xenon with (HF) and without (Frozen HF) accounting for the hole-particle interaction, and comparison to calculations by Bartlett et al.[44] and experimental results from Sorokin et al.[45] and Wetzel et al.[46].