# Fuelling the search for light dark matter–electron scattering
with spherical proportional counters

Louis Hamaide<EMAIL_ADDRESS>Christopher McCabe  Theoretical
Particle Physics and Cosmology Group, Department of Physics, King’s College
London, Strand, London, WC2R 2LS, United Kingdom

###### Abstract

Dark matter (DM) detectors employing a Spherical Proportional Counter (SPC)
have demonstrated a single-electron detection threshold and are projected to
have small background rates. We explore the sensitivity to DM–electron
scattering with SPC detectors in the context of DarkSphere, a proposal for a
300 cm diameter fully-electroformed SPC. SPCs can run with different gases,
so we investigate the sensitivity for five targets: helium, neon, xenon,
methane, and isobutane. We use tools from quantum chemistry to model the
atomic and molecular systems, and calculate the expected DM induced event
rates. We find that DarkSphere has the potential to improve current exclusion
limits on DM masses above 4 MeV by up to five orders of magnitude. Neon is the
best all-round gas target and provides good sensitivity to scenarios with both
light and heavy mediators. Gas mixtures, where methane or isobutane is added
to a noble gas, can extend the sensitivity at lower masses. Our study
highlights the currently untapped potential of SPCs to search for DM-electron
scattering in the MeV-to-GeV DM mass range.

††preprint: KCL-PH-TH-2021-31

##  I Introduction

Inferring more about the particle physics properties of dark matter (DM)
remains one of the most important open challenges in astro-particle physics
[1, 2]. It is in this context that the experimental program of DM direct
detection operates, since if a signal from DM is observed in a direct
detection experiment, we will learn about the value of the DM mass and DM
interaction with normal matter [3, 4, 5]. Traditionally, direct detection
experiments have focused on the search for nuclear recoils induced by
scattering from DM particles with the canonical WIMP mass range of around
5-1000 GeV (see, e.g., [6, 7, 8]). However, the lack of any conclusive DM
measurement has meant that during the past decade the mass range of interest
has broadened to encompass DM particles with a mass below and above this
narrow mass window [9].

The shift to searching for DM in the MeV to few-GeV mass range is challenging
experimentally, and a number of different avenues are being pursued. The
search for nuclear recoils induced by sub-5-GeV DM has continued by using
detectors with significantly smaller energy-detection thresholds, often in
combination with the utilisation of lighter target nuclei (see, e.g., [10, 11,
12, 13, 14, 15]). Alternatively, observable signals from sub-5-GeV DM can be
generated from rare processes that involve the emission of photons or ‘Migdal’
electrons from the recoiling atom (see, e.g., [16, 17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29]), or from the small flux of boosted DM that is non-
galactic in origin, or arises from interactions between DM and cosmic rays,
the Sun, or mesons (see, e.g., [30, 31, 32, 33, 34, 35, 36, 37]).

A complementary approach is to search for DM-electron scattering processes
that excite or ionise an electron from an atom or molecule, excite roto-
vibrational states or break chemical bonds in molecules, or excite an electron
across the band-gap in semi-conductors or other quantum materials (see, e.g.,
[38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 44, 48]). These
searches have the advantage that only a small amount of energy, of order the
ionisation energy or band-gap energy, needs to be transferred to the electron
to produce an observable signal. There already exists an impressive range of
existing constraints on the DM-electron cross section in the MeV-to-GeV mass
range, ranging from detectors designed specifically for this search channel to
detectors that have been designed first and foremost for nuclear recoil
searches but have been repurposed to search for DM-electron scattering [52,
53, 54, 55, 56, 57, 58, 59, 60].

It is in the latter context that we come to the main topic of this paper. The
New Experiments With Spheres-Gas (NEWS-G) detectors employ a spherical
proportional counter (SPC) [61, 62] and can be filled with different species
of atomic or molecular gases. Limits on nuclear-recoil signals from sub-GeV DM
have already been placed with a 60 cm diameter SPC [63], while a 140 cm
diameter detector, ‘S140’ (or sometimes SnoGlobe), which utilises ultra-
radiopure copper to reduce the background rate, is being commissioned at
SNOLAB [64, 65]. Funding has also been secured for ECuME, a planned follow-up
140 cm diameter detector that will benefit from a fully electroformed vessel
to lower the background rate, and a feasibility study of hosting a 300 cm
diameter fully electroformed SPC dubbed DarkSphere at the Boulby Underground
Laboratory is underway [66].

In its simplest form, the operation of the SPC is illustrated in fig. 1. The
detector consists of a grounded, spherical, metallic vessel filled with a gas
mixture and a spherical anode of radius $\sim 1$ mm is placed at the centre.
The anode produces an electric field that varies as $r^{-2}$ in an ideal SPC
and allows the ionisation electrons produced through particle interactions in
the gas volume to drift to the anode. When the electron reaches the anode, an
avalanche occurs, providing signal amplification, which is read out by a
grounded metallic rod. Larger SPCs will employ a multi-anode sensor (an
‘ACHINOS’ [67, 68]), which improves detector stability by decoupling the drift
field from the avalanche field, but the basic principle of operation is
similar to the simple SPC shown in fig. 1.111The electric field induced by a
60-anode ACHINOS, as envisaged for DarkSphere , varies from $\sim
0.1\,\mathrm{V/cm}$ at the outer sphere to $\sim 10\,\mathrm{V/cm}$ near the
anodes [69].

Owing to the low capacitance of the central sensor and the high amplification
gain near the anode, it has been demonstrated that this detector setup has a
single ionisation electron threshold [70, 71]. This threshold is key to the
$\mathcal{O}(15)$ eV threshold envisaged for nuclear recoil searches [71, 66],
where the recoiling nucleus produces electrons as it moves through the gas
volume. However, crucially, a single ionisation electron threshold also allows
for the possibility of searching for DM-electron interactions, where the
signal induced by the DM is a small number of electrons that are produced when
an atom or molecule in the gas mixture is ionised.

Figure 1: A schematic of the detection principle of DarkSphere, which can be
filled with a variety of atomic and molecular gases: we consider helium, neon,
xenon, methane, and isobutane. The outer sphere is a grounded, metallic shell,
and an anode is placed at the centre. DM scattering can ionise an electron
from an atom or molecule. After ionisation, the electron drifts to the anode
where an avalanche occurs, and the signal is read out. This setup has a
single-electron threshold, which is key to measuring the DM-induced signal.

In this paper, we quantify the range of DM masses and DM-electron scattering
cross sections that could be probed by DarkSphere, which we assume gives the
ultimate sensitivity from a large-SPC experiment. The DarkSphere SPC can
easily be filled with different gases; therefore, we study the response to
five different gases: helium, neon, and xenon, which are noble atoms commonly
used for DM experiments, and methane and isobutane (also known as
2-methylpropane), which may be more commonly known for their use as fuel. The
chemical formulae are He, Ne, Xe, CH4 and C4H10, respectively. Helium, neon,
and molecules with hydrogen and carbon atoms have been proposed because the
light nucleus of these atoms enhances the sensitivity to light DM nuclear-
recoil interactions. Meanwhile, xenon has been proposed as a target for
supernova-neutrino detection and neutrinoless double beta decay SPC searches
[72, 73]. For the noble gases, $\mathcal{O}(10\%)$ of methane or isobutane
will in general be added to the mixture as a ‘quench gas’. Quench gases are
used in proportional counters to improve detector stability by absorbing UV
photons that can be produced following the excitation of neutral atoms by
collisions between recoiling electrons or ions with the atoms in the gas.
Without the quench gas, the UV photons may lead to a loss of proportionality
and/or induce spurious signals (see e.g. [74] for a more extended discussion).
In contrast to the noble gases, the detectors could operate stably with pure
methane or pure isobutane as these gases can act as both the proportional and
quench gas. In addition to studying pure gas mixtures, which is advantageous
to clearly assess the advantages and disadvantages of each gas, we will also
study gas mixtures with $\mathcal{O}(10\%)$ of methane and isobutane as this
is the more likely gas-mixture to be employed in an operating detector.

Although the primary motivation for adding methane and isobutane comes from
the requirement of adding a quench gas, there are several reasons to believe
that these molecular targets may improve the sensitivity to DM-electron
interactions compared to atomic targets. Firstly, the ionisation energies in
molecules are often lower than the ionisation energies in the atomic targets
used in DM searches. For example, considering neon and methane which both have
10 bound-electrons, the ionisation energy of the outer six electrons in
methane is 13.6 eV compared to 23.1 eV for the outer six electrons in neon.
Isobutane has an even lower ionisation energy at approximately 11.1 eV,
compared with xenon’s 12.7 eV. Secondly, for a gas at the same pressure,
volume, and temperature, so that the number of moles is the same, a molecular
gas can benefit from a larger number of target electrons. For example,
counting the number of electrons with an ionisation energy less than 100 eV,
we find that neon has eight electrons, xenon has 18 electrons, while isobutane
has 26 electrons (although xenon has 54 atomic electrons in total, as we will
discuss, electrons with large ionisation energies cannot be ionised through DM
scattering). In this case, this means that the ratio of target electrons is
8:18:26 per mole of neon, xenon, and isobutane, respectively.

To the best of our knowledge, there has been no attempt to compare atomic and
molecular targets to see whether these apparent advantages are borne out.
Therefore, we will explore whether any improvements in the sensitivity to DM
masses and DM-electron scattering cross sections can be improved through the
use of methane and isobutane. While DM-electron scattering from atomic xenon
[41, 75, 76] has previously been considered, no results exist for helium,
neon, methane and isobutane.

The rest of this paper is structured as follows. In section II, we review the
formalism of electron-DM scattering and present our calculations of the bound-
and continuum-electron wave functions, as well as the ionisation form factors
for the atomic and molecular gases under consideration. In section III, we
present the DM-electron scattering rates and compare with the anticipated
DarkSphere background rate. Section IV presents our sensitivity projections
for DarkSphere and we discuss the possible advantages of gas mixtures over
single-species gases. We summarise our results and give our conclusions in
section V. Appendices provide more details on the numerical implementation of
our calculations, provide alternative calculations of the molecular ionisation
form factors, and compare our xenon calculations with others in the
literature.

Finally, many of our numerical results are available online. This includes
expressions for the bound-electron wave functions in position and momentum
space, interpolation tables for the ionisation form factors, and code to
calculate the scattering rates [77].

##  II Dark matter scattering with atomic and molecular electrons

We begin by presenting our calculation of the DM – electron scattering event
rate. In units of counts/keV/kg/year, the differential scattering event rate
for a DM particle of mass $m_{\mathrm{DM}}$ to ionise an electron from an atom
or molecule is

$\derivative{R}{E_{e}}=\frac{1}{m_{A}}\frac{\rho_{\mathrm{DM}}}{m_{\mathrm{DM}}}\sum_{i}w_{i}\,\derivative{\langle\sigma^{i\to
f}_{\rm{ion}}v_{\mathrm{DM}}\rangle}{E_{e}}\;,$ (1)

where $m_{A}$ is the mass of a target atom or molecule, $\rho_{\mathrm{DM}}$
is the local DM density, $w_{i}$ is the occupation number for electrons in the
atomic or molecular orbital labelled by the quantum number(s) $i$,
$\langle\sigma^{i\to f}_{\rm{ion}}v_{\rm{DM}}\rangle$ is the velocity-averaged
cross section for DM that scatters with an electron in the state labelled by
$i$ to the state $f$, and finally, $E_{\mathrm{e}}$ is the kinetic energy of
the continuum-electron. For an atom, $i$ would be the familiar $n,l,\cdots$
quantum numbers while for molecules, the quantum numbers depend on the
discrete symmetries of the system.

We follow the convention established in ref. [38] and parameterise the
differential cross section according to

$\derivative{\langle\sigma^{i\to f}_{\mathrm{ion}}v\rangle}{\ln
E_{e}}=\frac{\bar{\sigma}_{e}}{8\mu^{2}_{e}}\int_{q^{i}_{-}}^{q^{i}_{+}}q\,dq\,\absolutevalue{f^{i\to
f}_{\mathrm{ion}}}^{2}|F_{\rm{DM}}|^{2}g(v^{i}_{\mathrm{min}})\;.$ (2)

In this equation, $\bar{\sigma}_{e}$ is a model-independent reference cross
section, $\mu_{e}$ is the DM-electron reduced mass, $q$ is the amplitude of
the momentum transferred by the DM to the electron with upper and lower limits
$q^{i}_{+}$ and $q^{i}_{-}$, $F_{\rm{DM}}(q)$ is the dimensionless DM ‘form
factor’, $g(v^{i}_{\mathrm{min}}(q))$ is the DM velocity integral, and
$f^{i\to f}_{\mathrm{ion}}(E_{\mathrm{e}},q)$ is the dimensionless ionisation
form factor. For gas atoms or molecules in the DarkSphere SPC, the orientation
of the atoms and molecules will continuously change so we will always use the
rotationally averaged ionisation form factor in our calculations. Explicitly,
the general form of $f^{i\to f}_{\mathrm{ion}}(E_{\mathrm{e}},q)$ that we will
use is

$\begin{split}&\left|f^{i\to f}_{\rm{ion}}(E_{e},q)\right|^{2}\\\
&\quad=\left\langle\int d\Omega_{k_{e}}\frac{2k_{e}^{3}}{8\pi^{3}}\left|\int
d^{3}x\psi^{*}_{f}(\mathbf{x},\mathbf{k_{e}})e^{i\mathbf{q}\cdot\mathbf{x}}\psi_{i}(\mathbf{x})\right|^{2}\right\rangle,\end{split}$
(3)

where $k_{e}=\sqrt{2m_{e}E_{e}}$ is the momentum of the unbound electron
described by $\psi_{f}(\mathbf{x})$, the initial bound wave function is
$\psi_{i}(\mathbf{x})$, and the angled-brackets indicates that we are
averaging uniformly over all orientations of the atom or molecule. As
described further in appendix A, the rotational averaging procedure means that
$|f^{i\to f}_{\mathrm{ion}}|^{2}$ depends on the scalar $q=|\mathbf{q}|$
rather than the vector $\mathbf{q}$. We use the standard convention for the
bound-state wave functions that $\int
d^{3}x\,\psi^{*}_{i}(\mathbf{x})\psi_{i}(\mathbf{x})=1$ while the continuum
wave functions are normalised such that $\int
d^{3}x\,\psi^{*}_{f}(\mathbf{x},\mathbf{k})\psi_{f}(\mathbf{x},\mathbf{k^{\prime}})=(2\pi)^{3}\delta^{3}(\mathbf{k}-\mathbf{k^{\prime}})$.
Many authors work in the spherical-wave basis where the integral over the
solid angle of $k_{e}$ does not appear. However, we prefer to work in a basis
where the solid angle is explicit as it makes some of the manipulations in
appendix A more transparent. The final results are of course independent of
the choice of basis. It is also worth stating clearly that with our
convention, $f^{i\to f}_{\mathrm{ion}}(E_{\mathrm{e}},q)$ is the form factor
for a single electron (e.g. one electron in the $5p$-shell of xenon). This may
be different from other definitions in the literature that absorb the
occupation number $w_{i}$ into the definition of the form factor.

To provide the most straightforward comparison of the DarkSphere sensitivity
with other constraints, we use the Standard Halo Model (SHM) for the DM
velocity distribution with $v_{0}=235$ km/s [78, 79, 80, 81],
$v_{\rm{esc}}=544$ km/s (consistent with the latest estimates [82]), the Solar
peculiar velocity from ref. [83], and the average Earth velocity from ref.
[84]. The minimum speed parameter is

$v^{i}_{\rm{min}}(q)=\frac{E_{e}+I^{i}}{q}+\frac{q}{2m_{\rm{DM}}}\;,$ (4)

where $I^{i}$ is the ionisation energy of an electron in the initial-state
with quantum number(s) $i$, and $q^{i}_{\pm}$ are determined by solving eq.
(4) for $q$ with $v^{i}_{\rm{min}}=v^{\rm{max}}_{\rm{DM}}$, where
$v^{\rm{max}}_{\rm{DM}}$ is the maximum DM speed in the detector frame of
reference. This leads to

$q^{i}_{\pm}=m_{\rm{DM}}v^{\rm{max}}_{\rm{DM}}\left(1\pm\sqrt{1-\frac{E_{e}+I^{i}}{\tfrac{1}{2}m_{\rm{DM}}(v^{\rm{max}}_{\rm{DM}})^{2}}}\right).$
(5)

For our choice of parameters $v^{\rm{max}}_{\rm{DM}}=794~{}\mathrm{km/s}$. We
note that the SHM is of course an idealised model while the real Milky Way DM
halo contains much sub-structure not included within the SHM (e.g. see
discussion in refs. [85, 86, 87, 88, 89]). Finally, we set
$\rho_{\mathrm{DM}}=0.3~{}\mathrm{GeV}\mathrm{cm}^{-3}$ to facilitate
comparison with other experimental constraints.

###  II.1 Bound-state wave functions

As we have just discussed, determining the scattering rate requires the
ionisation form factor $f^{i\to f}_{\mathrm{ion}}(E_{\mathrm{e}},q)$, which in
turn depends on the bound and continuum wave functions. We will begin by
discussing our calculation of the bound-state wave functions. There are many
code packages capable of computing bound-state wave functions in atomic and
molecular systems using self-consistent field (SCF) methods. In this work, we
make use of PySCF [90, 91], a versatile open-source python-based quantum
chemistry package to calculate the bound-state wave functions for helium, neon
and xenon atoms and for methane and isobutane molecules. Although PySCF allows
for various SCF techniques to obtain the wave functions, we will stick to the
approach most often used in the DM literature: the non-relativistic Hartree-
Fock method.

It should be noted that the non-relativistic Hartree-Fock approach may have
limitations. Firstly, relativistic effects have been shown to become important
when the electron recoil energy is $\mathcal{O}(\mathrm{keV})$, which requires
that $q\langle r\rangle\gg 1$, where $\langle r\rangle\sim a_{0}/Z_{\rm{eff}}$
is the average position of the bound electron under consideration, $a_{0}$ is
the Bohr radius and $Z_{\rm{eff}}$ is the effective charge felt by the
electron. In this regime, non-relativistic calculations can underestimate the
rate by an order of magnitude or more [92, 75]. Using the comparison between
relativistic and non-relativistic calculations in ref. [93], we can estimate
that relativistic effects begin to become significant (meaning there is a
factor two or more difference) when $q\gtrsim 35Z_{\rm{eff}}$ keV. In this
work, we will consider electron recoil energies
$\mathcal{O}(10)~{}\mathrm{eV}$ and our calculations will be dominated by $q$
values $\sim 20~{}\mathrm{keV}$ (as we will see in fig. 3) so we do not expect
any significant error from not incorporating relativistic effects in our
calculations. Secondly, the Hartree-Fock method is limited in its accuracy
because it does not fully account for electron correlation. For instance, in
photoionisation cross section calculations with full shell atoms, the absence
of correlation effects limits the accuracy of Hartree-Fock calculations to
around 30% [94]. We therefore expect there to be a similar error in our
calculations.

PySCF expresses the wave functions in terms of a finite-dimensional Gaussian
basis set. In general, it is important to choose a large enough Gaussian basis
set to obtain accurate forms for the wave functions, particularly at small and
large values of $r$. For helium and neon, we therefore use the rather large
aug-cc-pV5Z basis set [95, 96], while for xenon, we use the Jorge-QZP basis
set [97] as aug-cc-pV5Z is not available for xenon. For methane and isobutane,
we use the 6-31G(d,p) set [98, 99, 100], which is sufficiently large to
capture the effects of the chemical bonds between the atoms. We have used
numerical coefficients from the Basis Set Exchange database [101, 102, 103].
For the methane and isobutane calculations, we also need to specify the
structure of the molecules. In our calculations, we use the geometric data
from the National Institute of Standards and Technology Computational
Chemistry Comparison and Benchmark Database (NIST CCCBDB) [104].

Helium (He) | Neon (Ne) | Methane (CH4)
---|---|---
Basis: aug-cc-pV5Z | Basis: aug-cc-pV5Z | Basis: 6-31G(d,p)
Tot. energy [au]†: $-2.8616$ | Tot. energy [au]†: $-128.5467$ | Tot. energy [au]†: $-40.2016$
Orbital |  $I_{\mathrm{HF}}$ [eV] |  $I_{\mathrm{exp}}$ [eV] | Orbital |  $I_{\mathrm{HF}}$ [eV] |  $I_{\mathrm{exp}}$ [eV] | Orbital |  $I_{\mathrm{HF}}$ [eV] |  $I_{\mathrm{exp}}$ [eV]
$1s^{2}$ | 24.98 | 24.6 | $2p^{6}$ | 23.14 | 21.7 | $1t_{2}^{6}$ | 14.80 | 13.6
|  |  | $2s^{2}$ | 52.53 | 48.5 | $2a_{1}^{2}$ | 25.66 | 22.9
|  |  | $1s^{2}$ | 891.79 | 870.2 | $1a_{1}^{2}$ | 304.96 | 290.8

|  | Isobutane (C4H10) | Xenon (Xe) |
---|---|---|---|---
|  | Basis: 6-31G(d,p) | Basis: Jorge-QZP |
|  | Tot. energy [au]†: $-157.3123$ | Tot. energy [au]†: $-7229.7195$ |
|  | Orbital |  $I_{\mathrm{HF}}$ [eV] |  $I_{\mathrm{exp}}$ [eV] | Orbital |  $I_{\mathrm{HF}}$ [eV] |  $I_{\mathrm{exp}}$ [eV] |
|  | $6a^{2}_{1}$ | 12.34 | 11.13 | $5p^{6}$ | 12.45 | 12.7 |
|  | $5e^{4}$ | 12.44 | 11.75 | $5s^{2}$ | 25.54 | 23.3 |
|  | $1a_{2}^{2}$ | 13.86 | 12.85 | $4d^{10}$ | 75.72 | 68.5 |
|  | $4e^{4}$ | 14.54 | 13.71 | $4p^{6}$ | 163.56 | 146.1 |
|  | $3e^{4}$ | 16.04 | 15.03 | $4s^{2}$ | 212.69 | 213.2 |
|  | $5a^{2}_{1}$ | 17.15 | 15.91 | $3d^{10}$ | 711.26 | 682.7 |
|  | $4a^{2}_{1}$ | 20.62 | 18.58 | $3p^{6}$ | 958.02 | 971.4 |
|  | $2e^{4}$ | 25.17 | 21.83 | $3s^{2}$ | 1087.7 | 1149 |
|  | $3a^{2}_{1}$ | 29.44 | 24.83 | $2p^{6}$ | 4839.8 | 4947 |
|  | $2a^{2}_{1}$ | 305.01 | — | $2s^{2}$ | 5132.0 | 5453 |
|  | $1e^{4}$ | 305.01 | — | $1s^{2}$ | 33321 | 34561 |
|  | $1a^{2}_{1}$ | 305.30 | — |  |  |  |

Table 1: A summary of the most pertinent details from our PySCF calculations.
In the second and third rows, we give the basis-set used and the total
converged energy in atomic units (au) for each atomic or molecule species. The
values $I_{\rm{HF}}$ and $I_{\rm{exp}}$ give the ionisation energy in eV from
the PySCF Hartree-Fock calculation and the experimental value from refs. [105,
106, 107, 108], respectively. The superscript in the orbital name gives the
number of electrons in that orbital. For isobutane, we do not have a reference
for the ionisation energies of the inner core orbitals (denoted by —).
Following ref. [107], we group the isobutane orbitals into three categories:
the outer valence orbitals ($6a_{1}$, $5e$, $1a_{2}$, $4e$, $3e$, $5a_{1}$),
the inner valence orbitals ($4a_{1}$, $2e$, $3a_{1}$), and the core orbitals
($2a_{1}$, $1e$, $1a_{1}$). The large number of orbitals (and therefore number
of electrons) in isobutane with ionisation energies less than
$30~{}\mathrm{eV}$ is much higher than in the other species.
† To provide for a more straightforward comparison of the total energy with
other computational approaches, we quote the total energy in atomic units,
where $1\,\mathrm{au}=27.211\,\mathrm{eV}$.

The main results of our PySCF calculations for helium, neon, methane,
isobutane and xenon are summarised in tab. 1.222An even more detailed summary
of the PySCF output is given on the GitHub repository [77], which also
includes the full analytic forms of the wave functions in position and
momentum space. The atomic orbitals are labelled by the principal and angular
momentum $n$ and $l$ quantum numbers. The labels for the methane and isobutane
orbitals follow from the discrete symmetry of the molecule. Methane orbitals
are eigenfunctions of irreducible representations of the discrete tetrahedral
group Td while the isobutane orbitals are eigenfunctions of irreducible
representations of the $C_{3v}$ point group.

For the atomic results, we can compare directly with the tabulated
computational results in ref. [109], which use Slater-type orbitals, and
against output from the ATSP atomic code [110]. We find that the total energy
and the individual orbital ionisation energies from our PySCF calculation,
$I_{\rm{HF}}$, agree very well, at the level of $0.001\%$ or better for helium
and neon, and at the level $0.03\%$ or better for xenon. We find excellent
agreement between the atomic wave functions when comparing them directly. We
also find very good agreement between the calculated methane total energy and
the computed ionisation energies in ref. [111], agreeing at the level of
$0.5\%$ or better. A similar level of agreement is found more generally
between our methane and isobutane total energy calculations and the suite of
calculations summarised in the NIST CCCBDB database [104]. Finally, we also
find good agreement between our spherically-averaged methane and isobutane
momentum-space wave functions and the various calculations presented in refs.
[112] and [107].

From tab. 1, we can also compare our calculated $I_{\rm{HF}}$ values with the
experimental values $I_{\exp}$. For atoms, we take the values from ref. [108]
and quote experimental values averaged over the fine-structure levels. The
experimental values for methane and isobutane are from refs. [106, 105, 107],
respectively, where we have quoted the vertical ionisation energies. We do not
have a reference for the experimental values of the inner-most core orbitals
of isobutane; however, these orbitals either do not contribute at all or
provide a sub-dominant contribution. The theoretical values are, however,
consistent with the inner-most methane values, which is what we would expect
given these are all carbon $1s$-like states.

We can take the analogy between methane and isobutane further by following
ref. [107], which groups the isobutane orbitals into three groups: the outer
valence orbitals ($6a_{1}$, $5e$, $1a_{2}$, $4e$, $3e$, $5a_{1}$), the inner
valence orbitals ($4a_{1}$, $2e$, $3a_{1}$), and the core orbitals ($2a_{1}$,
$1e$, $1a_{1}$). After this grouping, the average of the outer valence
ionisation energies is similar to $1t_{2}$ in methane, the inner valence
average is somewhat similar to $2a_{1}$ in methane, and as already noted the
core orbital average is similar to $1a_{1}$.

In general, we see that the theoretical values and experimental values do not
match exactly. This is the normal outcome of a non-relativistic Hartree-Fock
calculation since this approach misses many of the effects due to electron
correlation, which are required to obtain more accurate results. Electron
correlation can be accounted for with methods such as the multi-
configurational self-consistent field, which can also incorporate relativistic
effects (for a detailed discussion, see, e.g., ref. [113]). There are also
mature code packages such as the atomic code GRASP [114, 115] or the molecular
code BERTHA [116] that can calculate these effects. However, we leave the
exploration of multi-configurational self-consistent field techniques to
future work and in our calculations, we use the PySCF wave functions but,
where available, use the experimental values of the ionisation energy from
tab. 1.

Figure 2:  A visualisation of the outermost orbitals for helium, neon and
xenon, methane, and isobutane. The blue and orange regions indicate the
spatial extent of the electron orbital. The white and grey nodes represent the
positions of the atomic nuclei, and the connections between the nuclei in the
lower panels indicate the bonds. This visualisation was generated by
Mathematica [117] using the cubegen output from PySCF.

Finally, as the molecular orbitals are not usually encountered in the DM
literature, in fig. 2, we show a visualisation of the outermost orbitals for
helium, neon (and xenon), methane, and isobutane. The blue and orange shaded
regions indicate the typical spatial extent of the wave function, while the
smaller spheres indicate the position of the atomic nuclei. In the lower two
panels, the bonds between the hydrogen and carbon atoms are also shown.
Methane is a tetrahedron with four equivalent bonds between carbon at the
centre and the four outer hydrogen atoms. The tetrahedral geometry of the
highest occupied molecular orbital of methane is the result of
$sp^{3}$-hybridisation between electrons originally associated to the carbon
($2p$) and hydrogen ($1s$) states. The carbon atoms in isobutane form a
trigonal pyramid, where the carbon atom at the centre is bonded with three
carbon atoms at the base of the pyramid and one hydrogen atom at the top of
the pyramid. Each base carbon atom then has bonds with three hydrogen atoms,
forming tetrahedrons. The outermost orbital in isobutane is again formed as a
result of the sharing of electrons between the carbon $2p$ and hydrogen $1s$
electrons.

###  II.2 Continuum-electron wave functions

We now have the bound-state wave functions that enter eq. (3). Next, we
describe our calculations of the continuum-electron wave functions. For both
atomic and molecular systems, we look for outgoing states that have been
expanded into partial waves,

$\psi_{f}(\mathbf{x},\mathbf{k_{e}})=\sum_{l,m}i^{l}\frac{P_{k_{e}l}(r)}{r}Y^{*}_{lm}(\hat{\mathbf{x}})Y_{lm}(\hat{\mathbf{k}}_{e})\;,$
(6)

where $P_{k_{e}l}(r)$ is the radial part of the wave function and $Y_{lm}$ are
the spherical harmonics. The normalisation condition $\int
d^{3}x\,\psi^{*}_{f}(\mathbf{x},\mathbf{k})\psi_{f}(\mathbf{x},\mathbf{k^{\prime}})=(2\pi)^{3}\delta^{3}(\mathbf{k}-\mathbf{k^{\prime}})$
holds if the radial wave functions satisfy $\int
drP_{kl}(r)P_{k^{\prime}l^{\prime}}(r)=(2\pi)^{3}\delta_{ll^{\prime}}k^{-2}\delta(k-k^{\prime})$.
We impose this by following the standard procedure of ensuring that at large
values of $r$ our solution $P_{k_{e}l}(r)$ asymptotes to a sine function with
amplitude $4\pi/k_{e}$ [118].

We first discuss our numerical approach to finding atomic states, where, owing
to the spherical symmetry of atoms, the expansion in eq. (6) is exact. The
continuum-electron wave functions $P_{nl\to k_{e}l_{e}}(r)$ for an electron
initially in the state $n$, $l$ that has been promoted to an continuum state
with quantum numbers $k_{e}$, $l_{e}$ are obtained by numerically solving an
approximate form of the one-electron Hartree-Fock equation in the mean field
of the remaining electrons. Explicitly, the equation that we solve, written in
atomic units, is

$\begin{split}&\left[-\frac{1}{2}\frac{d^{2}}{dr^{2}}+\frac{l_{e}(l_{e}+1)}{2r^{2}}+V_{nl\to
k_{e}l_{e}}(r)\right]P_{nl\to k_{e}l_{e}}(r)\\\ &\qquad\qquad=E_{e}P_{nl\to
k_{e}l_{e}}(r)+\sum_{n^{\prime}l^{\prime}}\delta_{l_{e}l^{\prime}}\lambda_{k_{e}n^{\prime}}P_{n^{\prime}l^{\prime}}(r)\;,\end{split}$
(7)

where $\lambda_{k_{e}n^{\prime}}$ are Lagrange multipliers that are introduced
to ensure that radial wave functions with the same $l$ quantum numbers obey
the correct orthogonality relations.

The potential energy in the full Hartree-Fock equations contains both a local
and non-local term [119]. For the continuum electron, we employ Cowan’s
Hartree-plus-statistical-exchange (HX) method to approximate the non-local
exchange potential with a local potential, which dramatically reduces the
complexity of the problem while producing results in good agreement with the
full Hartree-Fock method [120, 121]. Cowan’s HX method is particularly suited
to ionisation problems since the potential energy is guaranteed to have the
expected $-1/r$ relation for an electron far from the atom.

In the HX method, the potential energy is approximated as

$V_{nl\to k_{e}l_{e}}(r)=-\frac{Z}{r}+V^{\rm{H}}_{nl\to
k_{e}l_{e}}(r)+V^{\rm{HX}}_{nl\to k_{e}l_{e}}(r)\;,$ (8)

where the first term is the familiar Coulomb interaction of the electron with
a nucleus of charge $Z$, $V^{\rm{H}}(r)$ is the classical potential energy of
an electron in the mean field of the other atomic electrons and
$V^{\rm{HX}}(r)$ is the local approximation to the exchange potential.

The classical potential energy in the frozen core approximation is

$V^{\rm{H}}_{nl\to
k_{e}l_{e}}(r)=\sum_{n^{\prime}l^{\prime}}(w_{n^{\prime}l^{\prime}}-\delta_{nl,n^{\prime}l^{\prime}})\int_{0}^{\infty}\frac{dr^{\prime}}{r_{>}}P^{2}_{n^{\prime}l^{\prime}}(r^{\prime})\,,$
(9)

where $r_{>}=\text{min}(r,r^{\prime})$ and $w_{nl}$ is the electron occupation
number. The sum is over all bound electrons, $P_{nl}(r)$ are the bound wave
functions from PySCF (i.e. in the ionisation calculation, we ‘freeze’ the wave
functions of the bound core electrons), and the Kronecker delta excludes the
bound electron that has been ionised from contributing to the potential.

The HX exchange potential energy in the frozen core approximation for
continuum states is

$V^{\rm{HX}}_{nl\to
k_{e}l_{e}}(r)=-\frac{k_{x}}{2}\left(\frac{\rho^{\prime}(r)}{\rho(r)}\right)\left(\frac{24\rho(r)}{\pi}\right)^{1/3},$
(10)

where $k_{x}$ is an $\mathcal{O}(1)$ numerical factor, and

$\displaystyle\rho(r)$
$\displaystyle=\sum_{n^{\prime}l^{\prime}}w_{n^{\prime}l^{\prime}}\frac{P_{n^{\prime}l^{\prime}}^{2}(r)}{4\pi
r^{2}}\;,$ (11) $\displaystyle\rho^{\prime}(r)$
$\displaystyle=\sum_{n^{\prime}l^{\prime}}\left[w_{n^{\prime}l^{\prime}}-\text{min}(2,w_{nl})\delta_{nl,n^{\prime}l^{\prime}}\right]\frac{P_{n^{\prime}l^{\prime}}^{2}(r)}{4\pi
r^{2}}\;.$ (12)

Again, the sum is over all bound electrons, while the $\text{min}(2,w_{nl})$
term ensures that the exchange potential vanishes for any two-electron
configuration, e.g., as found in helium, as exchange terms only arise between
electrons with parallel spin. As with the classical potential energy
calculation, $P_{nl}(r)$ are PySCF bound wave functions, and the Kronecker
delta excludes the self-interaction of the electron that has been ionised. We
fix $k_{x}=0.65$ as this leads to the most accurate results [120, 121, 122].

Next, we discuss our approach to calculating continuum states from molecules.
A rigorous calculation of the continuum-electron states in a molecule is a
much more challenging calculation. Firstly, the geometry of methane and
isobutane means that the potential energy is no longer exactly spherically
symmetric, and, secondly, the potential contains multiple poles at the
positions of the atomic nuclei. In this work, we take a simplified approach
and model the outgoing electron with a single-centred Coulomb wave function,
which is an analytic solution to the Schroedinger equation for a $-Z/r$
potential. This simplified approach has previously been used in the DM
literature for atomic states (see, e.g., [54, 123]). In atomic units, the
Coulomb function takes the form

$\begin{split}P_{kl}(r)&=\frac{4\pi}{2k}\frac{\left|\Gamma\left(\ell+1-\frac{iZ}{k}\right)\right|e^{\frac{\pi
Z}{2k}}}{(2\ell+1)!}(2kr)^{\ell+1}\\\ &\times
e^{-ikr}\,M\left(\ell+1+\frac{iZ}{k},2\ell+2;2ikr\right)\;,\end{split}$ (13)

where $M(a,b;x)$ is Kummer’s confluent hypergeometric function and $\Gamma(z)$
is the Gamma function.

The charge $Z$ remains a free parameter in the Coulomb function. We again opt
for a simple and transparent method using the ionisation energies to choose
$Z$ for the methane and isobutane orbitals. Using the hydrogen-like scaling
for the ionisation energy, $I=Z^{2}/n^{2}\times 13.6~{}\mathrm{eV}$ with $n=1$
for $1a_{1}$ and $1t_{2}$ and $n=2$ for $2a_{1}$, leads to $Z$ values of
$4.7$, $2.6$ and $1.0$ for the $1a_{1}$, $2a_{1}$ and $1t_{2}$ states in
methane, respectively. For isobutane, we follow ref. [107] and use the
grouping of orbitals into inner core states ($1a_{1}$, $1e$, $2a_{1}$), inner
valence states ($3a_{1}$, $2e$, $4a_{1}$), and outer valence states (all
remaining orbitals). From the ionisation energies and using $n=1$ for the
inner core and outer valence states and $n=2$ for the inner valence states, we
find that $Z=4.7$ for the inner core orbitals, $Z=2.7$ for the inner valence
orbitals, and $Z=1$ for the outer valence orbitals. We also have to specify
the centre of the Coulomb wave functions. In our calculations, the centre is
defined by the coordinate system from the NIST CCCBDB database [104], which
takes the centre-of-nuclear-charge as the origin. In methane, this corresponds
to functions centred on the carbon nucleus, while in isobutane, the centre is
in the orange region directly below the central carbon nucleus at the top of
trigonal pyramid in the lower right panel of fig. 2.

The Coulomb wave functions have some obvious limitations. Firstly, the $-Z/r$
Coulomb potential misses the $r$-dependent contributions from $V^{\mathrm{H}}$
and $V^{\mathrm{HX}}$ in eq. (8). Secondly, the assignment of the $Z$ value
for molecules is less robust than for atomic states. Thirdly, the Coulomb
functions are not guaranteed to have the correct orthogonality conditions for
the bound-state functions, as imposed by the Lagrange multipliers in eq. (7).
Fourthly, the Coulomb potential is spherically symmetric, while the methane
and isobutane potentials are non-spherical. And finally, the Coulomb waves may
not be centred at the centre-of-nuclear-charge but instead at one of the
nuclei.

Figure 3:  The dimensionless ionisation form factor for helium (top-left),
neon (middle-left), xenon (bottom-left), methane (middle-right) and isobutane
(bottom-middle). In all cases, we show results when the continuum-electron
kinetic energy $E_{e}$ is 20 eV. For clarity in the xenon panel, we only show
the results for the three outermost shells, $4d$, $5s$ and $5p$, as these
usually give the dominant contribution. In the isobutane panel, for clarity,
we show two results from the outer valence orbitals ($6a_{1}$ and $5e$) to
demonstrate that they are very similar, one inner valence ($2e$) and one core
orbital ($2a_{1}$). The diamonds and triangles on some lines indicate the
$q_{-}$ and $q_{+}$ values, respectively, for $m_{\mathrm{DM}}=30$ MeV. There
are similarities in terms of the peak position and overall shape between the
form factors for the outermost orbitals in helium ($1s$), neon ($2p$) and
xenon ($5p$), and between all three form factors in neon and methane. The
shape of the methane and isobutane orbitals are also somewhat similar. Note
however that the form factors for methane and isobutane drop-off more rapidly
than for neon as $q$ increases.

Employing a spherically averaged potential is a commonly used approximation
even in more advanced treatments of scattering problems involving molecules
(see, e.g., [124, 125, 126, 127, 128]) and typically leads to results with an
accuracy of approximately $50\%$ or better compared to more rigorous methods.
We explore the uncertainty induced by changing the $Z$ value and by moving the
origin of the Coulomb waves in appendix B.

###  II.3 Ionisation form factors

Having obtained the bound- and continuum-electron wave functions, we are now
able to calculate the rotationally averaged ionisation form factor $|f^{i\to
f}_{\rm{ion}}(E_{e},q)|^{2}$. The general expression was given in eq. (3);
however, to generate the results in this section, we use the equivalent forms
derived in appendix A as they are more amenable to numerical calculations. We
have made extensive use of the GNU Scientific Library-2.7 [129], Cubature
[130], and GNU Parallel [131] in our numerical calculations.

In fig. 3, we plot the results of our ionisation form factor calculation for
atoms (left column) and molecules (right column). For clarity, we show only a
subset of the orbitals for xenon and isobutane (bottom panels). In the xenon
panel, we show the outer three orbitals as these are the most important for DM
scattering. In the isobutane panel, we show two of the outer valence orbitals
($6a_{1}$ and $5e$) to show that they are almost identical, differing only
slightly at $q\lesssim 5$ keV. The other outer valence orbitals (not plotted)
are also very similar. Similarly, we plot one inner valence orbital ($2e$) and
one core orbital ($2a_{1}$) as the other inner valence and core orbitals are
similar.333Tabulated values of the ionisation form factors for all of the
helium, neon, xenon, methane and isobutane orbitals are available on the
GitHub repository [77].

In fig. 3, we show results when the continuum electron’s kinetic energy is
$E_{e}=20$ eV. At smaller values of $E_{e}$, the ionisation form factors have
a similar shape but with a smaller amplitude and a small shift by
$\sim\mathcal{O}(\mathrm{keV})$ to smaller values of $q$. This change in the
amplitude occurs because of the dependence of the electron kinetic energy in
the definition of the dimensionless ionisation form factor, eq. (3) (recall
that $E_{e}$ enters through $k_{e}$), while the shift to smaller values of $q$
occurs because a smaller continuum-electron kinetic energy requires a smaller
kick from the momentum transfer. Conversely, for larger values of $E_{e}$, the
amplitude of the form factors is larger and the peak shifts toward larger
values, e.g., when $E_{e}\sim 1$ keV, the peak of the outer orbitals shifts to
the range $q\sim 30$ to $50$ keV.

Comparing the atomic form factors, we see that the results of the outermost
orbitals ($1s$, $2p$ and $5p$ for helium, neon and xenon, respectively) share
similarities in the height and position of the peak values. This should not be
too surprising, given that outer-shell atomic electrons of different atoms
find themselves in a similar screened-Coulomb potential. Comparing neon to
methane and isobutane, we also see similarities between all three form
factors. Finding some similarities between neon and methane may not be too
surprising as both are systems with 10 electrons. Furthermore, after grouping
the isobutane orbitals into the outer valence, inner valence, and core
orbitals, we already observed similarities between the methane and isobutane
ionisation energies, and here we see similarities between the form factors.
However, there are two important differences particularly at small and large
values of $q$ that deserve further comment.

The first important difference, at $q\lesssim 10$ keV, is that the methane and
isobutane form factors do not show the steady increase with $q$ observed with
the atomic results. Instead, in the molecular form factors there is a gradual
flattening to a finite value as $q\to 0$. This is a result of using Coulomb
wave functions in the methane and isobutane calculations, which only
approximately obey the orthogonality conditions expected between the bound and
continuum wave functions; it should be the case that when $q=0$ the overlap
integral in eq. (3) is zero. Instead, when using Coulomb functions for the
continuum-electron, the overlap integral gives a constant non-zero value when
$q=0$.

However, for the DM scattering rates, this anomalous behaviour at small $q$ is
not important. We can understand why by considering the diamonds and triangles
marked on some of the lines in fig. 3. The diamonds show the $q_{-}$ value,
while the triangles show the $q_{+}$ values, which are the limits of
integration that appeared in eq. (2), for $m_{\mathrm{DM}}=30$ MeV. Orbitals
that do not have a diamond or triangle do not contribute to the differential
cross section at this value of $m_{\mathrm{DM}}$. For some orbitals (e.g.,
helium 1s), the triangle does not appear simply as a result of our plotting
scale for the form factor. For all targets, we see that the typical $q_{-}$
value is approximately 10 to 20 keV and therefore above the region where the
Coulomb functions could introduce an error in the molecular calculations from
the anomalous orthogonality conditions. Figure 3 shows this explicitly for
$m_{\mathrm{DM}}=30$ MeV and a single value of $E_{e}$, but we can use eq. (5)
to show that the smallest value of $q_{-}$ is approximately
$11~{}\mathrm{keV}\cdot(E_{e}+I^{i})/(30~{}\mathrm{eV})$, from which we
conclude that the anomalous region $q\lesssim 5$ keV does not enter the
integration over $q$ in eq. (2).

The second and more important difference between neon, methane and isobutane
is at larger values of $q$, where fig. 3 shows that the form factors drop-off
more rapidly as $q$ increases. This is particularly evident when comparing the
outer shell orbitals ($2p$, $1t_{2}$, and $6a_{1}$). This effect can be
understood as follows. Rewriting the form factor (eq. (3)) in momentum space,
we would find that $|f_{\rm{ion}}^{i\to f}|^{2}\sim$
$\absolutevalue{\tilde{\psi}_{i}(\mathbf{k_{e}}-\mathbf{q})}^{2}$, where
$\tilde{\psi}$ is the momentum-space wave function (see appendix A.1 for an
explicit realisation of this when the outgoing electron is modelled by a
plane-wave). This implies that large values of $q$ correspond to a probe of
the large-momentum tail of the momentum-space wave function. The outermost
momentum-space wave functions fall away more quickly for methane and isobutane
at large values of the momentum compared to neon, helium, and xenon, so the
molecular form factors fall away more rapidly.

##  III Comparison of scattering rates

Figure 4:  The differential event rates for DM-induced ionisation for DM with
mass 10 MeV (top), 30 MeV (middle), and 300 MeV (bottom). Rates are shown for
helium (dot-dashed orange), neon (light blue), xenon (dashed dark-blue),
methane (dotted green), and isobutane (dot-dot-dashed pink). In all three
cases, we show results for $\bar{\sigma}_{e}=10^{-41}$ cm2 and
$F_{\mathrm{DM}}=1$. Neon tends to have the largest rate, except at the lowest
mass (top panel), where both neon and helium suffer a suppression from the
relatively high ionisation energy of the outer shell. The features in the
xenon curve arise from the interplay of the $5p$, $5s$ and $4d$ orbitals. The
grey shaded region is an estimate of the background rate for DarkSphere
assuming that it is situated in the Large Experimental Cavern at the Boulby
Underground Laboratory.

With the ionisation form factors in hand, it is now straightforward to
calculate the differential cross section using eq. (2) and the differential
scattering rate using eq. (1). In fig. 4, we plot the scattering event rates
as a function of the ionised electron’s kinetic energy $E_{e}$ for the five
different gases assuming $\bar{\sigma}_{e}=10^{-41}~{}\mathrm{cm}^{2}$,
$F_{\rm{DM}}=1$ and for three values of the DM mass: $m_{\rm{DM}}=10$ MeV
(top), $m_{\rm{DM}}=30$ MeV (middle), and $m_{\rm{DM}}=300$ MeV (bottom). For
all targets except helium, more than one orbital contributes to the scattering
rate. This is most apparent in the xenon rate (dashed dark-blue line) in the
bottom panel, where the features at different energies arise from the
interplay of the $5p$, $5s$, and $4d$ orbitals, which have been summed
together.444Appendix C provides a comparison of our xenon scattering rate
calculation with other xenon calculations in the literature.

For $m_{\rm{DM}}=30$ and 300 MeV, we see that the largest rate occurs for neon
(solid blue line), except for the high-energy tail above $E_{e}\sim 60$ eV for
$m_{\mathrm{DM}}=300$ MeV where xenon takes over. Neon has a number of
advantages compared to the other gases. Compared to helium, there are three
times as many electrons available in the outer-shell. Compared to methane and
isobutane, its form factor falls off more slowly at large $q$, as discussed in
the previous section. And compared to xenon, the rate suffers less of a
suppression by the $1/m_{A}$ factor in the scattering rate (cf. eq. (1)).

Figure 5:  Similar to fig. 4 but with $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ and
$\bar{\sigma}_{e}=10^{-39}~{}\mathrm{cm}^{2}$. For this form factor, the rates
fall off much more rapidly as $E_{e}$ increases compared to $F_{\rm{DM}}=1$.
Methane and isobutane tend to dominate at low-energies (below $\sim 10$ eV)
while above this energy, the rate for neon or xenon is larger.

For DM masses between 300 MeV and 1 GeV, the scattering rates have the same
shape but scale with the $1/m_{\mathrm{DM}}$ factor in eq. (1). At
$m_{\rm{DM}}=10$ MeV (top panel), the situation is different. Both the neon
and helium rates (solid blue and dashed-dot orange lines, respectively) are
more suppressed relative to the other gases. This follows as a result of the
larger outer-shell ionisation energies in neon and helium compared to the
other targets (cf. tab. 1). Comparing methane and isobutane (dashed green and
dot-dot-dashed pink lines, respectively), the relative scattering rates again
depend on the competition between more electron targets in isobutane versus
the larger suppression from $1/m_{A}$.

In fig. 5, we plot the scattering event rates for the five gases but this time
for $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$, where $\alpha$ is the fine-structure
constant. This form for $F_{\rm{DM}}$ is appropriate for models that have a
light-particle mediating the interaction between DM and electrons. In fig. 5,
we also choose a larger value of the reference cross section,
$\bar{\sigma}_{e}=10^{-39}~{}\mathrm{cm}^{2}$, but this is still below current
bounds for this form of $F_{\rm{DM}}$ (cf. fig. 8). The top, middle, and
bottom panels again show the scattering rates for $m_{\rm{DM}}=10$, 30, and
300 MeV, respectively.

Comparing figs. 4 and 5, there are several obvious differences. Firstly, the
rates in fig. 5 all drop off more quickly as $E_{e}$ increases. This is
especially noticeable at $m_{\rm{DM}}=300$ MeV. The rapid fall-off occurs
because an increase in $E_{e}$ leads to an increase in $q_{-}$, and when the
ionisation form factors are multiplied by the additional $1/q^{4}$ factor, the
cross section integrands drop off even more quickly with $q$, so the
scattering rates are even more sensitive to $q_{-}$. Secondly, in fig. 5 the
largest rates at small values of $E_{e}$ occur for methane and isobutane at
all three values of the mass shown, while for $F_{\rm{DM}}=1$, methane and
isobutane were the largest only for $m_{\rm{DM}}=10$ GeV. This is again
directly related to the sensitive dependence on $q_{-}$. As neon and helium
have larger values of the ionisation energy, this in turn means that the
$q_{-}$ values are larger (cf. eq. (5)). Although xenon also has a low
ionisation energy, the larger $1/m_{A}$ suppression means that its rate is
below methane and isobutane.

The sensitivity to DM-electron scattering will ultimately be limited by the
background rate. The S140 and ECuME SPC experiments are projected to achieve
background rates of $\mathcal{O}(1)~{}\\!\\!\text{cts/keV/kg/day}$ and
$\mathcal{O}\left(10^{-1}\right)~{}\\!\\!\text{cts/keV/kg/day}$, respectively,
in the energy range below 1 keV [65]. DarkSphere aims to achieve a background
rate of $\mathcal{O}\left(10^{-2}\right)~{}\\!\text{cts/keV/kg/day}$ below 1
keV. The grey line and shaded region in figs. 4 and 5 show an energy-
independent background rate of $10^{-2}~{}\\!\text{cts/keV/kg/day}$ below 1
keV, which we will use in our sensitivity calculations for all gas mixtures.

A preliminary study [66] shows that a background rate of
$10^{-2}~{}\\!\text{cts/keV/kg/day}$ could be achieved with a fully
electroformed-underground copper shell and by shielding the $3$ m diameter of
DarkSphere with a $2.5$ m thick water tank in the Large Experimental Cavern
(LEC) at the Boulby Underground Laboratory. After fully electroforming the
copper shell underground, it is expected that contaminants from 210Pb and
progeny decays would contribute a background rate $\sim 2\times
10^{-5}~{}\\!\text{cts/keV/kg/day}$ below 1 keV [66], which is significantly
smaller than the targeted rate. Instead, the dominant background is expected
to arise from environmental backgrounds, specifically from neutrons, photons
and muons that scatter in the SPC. Muons are present in cosmic rays while the
neutrons and photons are generated in the cavern rock, either from
radioactivity or from cosmic-ray interactions (see refs. [132, 133, 134, 135]
for related studies). A Geant4 [136] simulation assuming a $2.5$ m thick water
tank in the LEC at Boulby surrounding the DarkSphere SPC, and assuming a
90%–10% He–C4H10 gas mixture, predicts a photon, neutron and muon-induced
environmental background rate below 1 keV of $\sim 5\times
10^{-3}~{}\\!\text{cts/keV/kg/day}$ [66]. The dominant contribution arises
from 2-to-3 MeV photons, which produce an almost-flat spectrum below 1 keV.

Although both the contaminants in the copper shell and the environmental
backgrounds are below the targeted rate, there are two further sources of
background that could increase the background rate above
$10^{-2}~{}\\!\text{cts/keV/kg/day}$ below 1 keV. The first is from potential
radioactive contaminants in the gas mixture, with ${}^{14}\mathrm{C}$ and
${}^{3}\mathrm{H}$ being of particular concern (see, e.g., ref. [137]).
Tritium can be produced in neon when exposed to cosmic rays and estimates in
ref.[135] show that the gas should remain on the surface for less than one
month to ensure that the background from ${}^{3}\mathrm{H}$ remains sub-
dominant. Methane or isobutane could provide more of a challenge since they
are composed solely of C and H isotopes and either gas will generally be
present in their role as a quench gas. Both methane and isobutane are produced
from underground natural gas deposits, which will have been underground for
geological timescales, so the intrinsic contamination from ${}^{14}\mathrm{C}$
and ${}^{3}\mathrm{H}$ should only be introduced when extracted, manipulated
on the surface, or transported. The NEWS-G Collaboration will need to work
with gas-manufacturing companies to ensure that sufficiently pure gas is
procured and that time on the surface is minimised.

The second background source that could increase the rate beyond
$10^{-2}~{}\\!\text{cts/keV/kg/day}$ is the excess low-energy background below
$\mathcal{O}(100)~{}\\!\mathrm{eV}$ that has been observed in many DM
experiments [138]. Indeed, a test run of the S140 detector with methane gas
found a large excess of events in the single-electron channel (corresponding
to an energy below $\sim 28$ eV), with an event rate orders of magnitude
larger than in the two-electron channel. The single-electron events cannot
currently be removed by cuts, so physics searches may be limited to the two-
electron channel (corresponding to an energy $\gtrsim 28$ eV) [138].

With these caveats in mind, we will proceed assuming the energy-independent
background rate of $10^{-2}~{}\\!\text{cts/keV/kg/day}$ below 1 keV as shown
in figs. 4 and 5. The smallest value of the cross section that can be
constrained will be approximately determined by the point at which the signal
rate is the same magnitude as the background rate (we will use a more rigorous
likelihood-based method in the next section to more carefully quantify the
sensitivity). The reference cross section
$\bar{\sigma}_{e}=10^{-41}~{}\mathrm{cm}^{2}$ used in fig. 4 is below current
constraints for the three values of the mass shown (cf. fig. 7), so since all
of the rates in fig. 4 are above the background rate for some values of
$E_{e}$, this indicates that DarkSphere has the potential to probe new regions
of parameter space. Similarly, we see that the signal rates in fig. 5 are
above the background rate for a cross section below current bounds, so we
anticipate that DarkSphere has the potential to probe new regions of parameter
space for the $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ form factor.

Figure 6:  The total (integrated) number of events above a given threshold
energy as a function of the threshold energy for DM with mass 10 (top), 30
(middle), and 300 MeV (bottom). The number of events is shown for
$F_{\mathrm{DM}}=1$ in the left panel and $\bar{\sigma}_{e}=10^{-40}$ cm2, and
for $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ and
$\bar{\sigma}_{e}=10^{-38}~{}\mathrm{cm}^{2}$ in the right panel. We have
assumed the baseline DarkSphere parameters: an SPC with a 3 m diameter filled
with the target gas at a pressure of 5 bar, temperature of $20^{\circ}$C, and
taking data for 300 days. Over most of the parameter space xenon gives the
largest number of events as it has the largest target mass (380.8 kg), while
helium gives the smallest number of events as it has the smallest target mass
(11.6 kg). The lines in the right panel drop off more rapidly as
$E_{{\rm{th}}}$ increases compared to the left panel, indicating that a
smaller threshold energy is required to observe events when
$F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$.

##  IV Sensitivity projections

In this section, we will present projections for DarkSphere in the DM mass –
cross section parameter space. Before showing the results, we should clarify
the expected behaviour inside the DarkSphere detector after ionisation from
the atom or molecule has occurred (see ref. [74] for a more in-depth
discussion). After ionisation, the electric field will cause the ionised
electron to drift to the anode, where a Townsend avalanche is triggered when
the electron is within a few millimetres of the anode. This avalanche produces
a large number of secondary electron-ion pairs that are detected as a charge
pulse. However, there are two additional processes that may induce additional
primary electrons, which if produced, will also drift toward the anode and
increase the size of the pulse. Firstly, if the kinetic energy of the ionised
electron is larger than the energy required to create a new electron-ion pair
in the gas, additional primary electrons can be produced. Secondly, if an
electron is ionised from an inner-orbital of the atom or molecule, a higher-
energy bound electron will transition to the hole, and additional energy will
be released, e.g., in the form of a photon or Auger electron. An Auger
electron is obviously an additional electron, while more primary electrons can
be produced if the photon or Auger-electron energy is larger than the
electron-ion pair threshold energy.

Ideally, we would show scattering rates in terms of the experimentally
measurable parameters such as the number of primary electrons created or the
number of photoelectrons detected. However, this requires a precise
characterisation of the detector response at very low energies, which is not
yet available for all of the gases and gas mixtures that we will consider
(although important steps in this direction have been taken in ref. [71,
135]). Therefore, we will instead adopt a simpler approach and discuss the
experimental sensitivity only in terms of the kinetic energy of the ionised
electron ($E_{e}$). This means that our treatment is somewhat conservative as
it ignores the contribution from any additional primary electrons that could
be produced, e.g., in the case that an inner-orbital is ionised. This is
likely to lead to an under-estimate of the sensitivity, particularly for xenon
where the difference between the orbital ionisation energies can be relatively
large ($\gtrsim 100$ eV).

Characterising the sensitivity estimates in terms of $E_{e}$ also means that
we should define the analysis threshold in terms of $E_{e}$. To investigate
the dependence of the number of signal events on the threshold energy, in fig.
6 we have plotted the number of events above a given threshold energy as a
function of the threshold energy ($E_{\rm{th}}$). The left panels are for
$F_{\rm{DM}}=1$ and $\bar{\sigma}_{e}=10^{-40}~{}\mathrm{cm}^{2}$, while the
right panels are for $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ and
$\bar{\sigma}_{e}=10^{-38}~{}\mathrm{cm}^{2}$. The top, middle, and lower
panels are for $m_{\mathrm{DM}}=10,~{}30$, and $300~{}\mathrm{MeV}$,
respectively. To calculate the total number of events, we have assumed the
baseline detector and operating parameters for DarkSphere from [66]; this
assumes an SPC with a 3 m diameter (⌀) filled with the target gas at a
pressure of 5 bar, temperature of $20^{\circ}$C and taking data for 300 days.
The same pressure for the different gases implies that the total mass inside
the detector is different for each gas. We find that there should be
approximately 2900 mol for each gas, corresponding to 11.6, 58.5, 380.8, 46.5
and 168.6 kg for helium, neon, xenon, methane, and isobutane, respectively.

From fig. 6 we observe the general pattern that xenon gives the largest number
of events and helium gives the smallest number. This follows primarily from
the large difference in the target mass of the two targets. We also see that
for the chosen cross sections, which lie below current exclusion bounds, the
total number of events can be rather large ($\gtrsim 1000$ events). When
$m_{\rm{DM}}=10$ MeV (top panels), we see for both choices of $F_{\rm{DM}}$
that the number of events falls off rapidly as the threshold energy increases,
and falls to zero completely for $E_{\rm{th}}\gtrsim 20$ eV. For higher masses
and when $F_{\rm{DM}}=1$, the left-middle and left-lower panels show that the
drop off in the number of events as the threshold increases is much shallower,
such that the number if events is still above one (and still significantly
above one for neon and xenon gases) at $E_{\mathrm{th}}\simeq 30$ eV.
Conversely, the right panels show that when $F_{\rm{DM}}=(\alpha
m_{e}/q)^{2}$, the drop off in the number of events remains significant, even
at larger masses. This is a reflection of the behaviour observed previously in
the discussion of fig. 5.

###  IV.1 Single-species projections

We show in fig. 7 our projections for the 90% confidence limit (C.L.)
exclusion limits for pure gas mixtures assuming $F_{\rm{DM}}=1$ for the
baseline detector and operating parameters for DarkSphere from [66]. As a
reminder, this assumes an SPC with a 3 m diameter filled with the target gas
at a pressure of 5 bar, temperature of $20^{\circ}$C and taking data for 300
days with the background rate shown in figs. 4 and 5.

The two panels in fig. 7 show projections under two assumptions for the
detection threshold. The upper panel in fig. 7 assumes a threshold energy
$E_{\rm{th}}=1$ eV. This means that essentially all electrons that have been
ionised will be pulled by the electric field in the DarkSphere drift region
towards the anode, where the signal will be read out. This threshold is chosen
to mimic a single-electron search threshold since electrons with a kinetic
energy around $10$ eV or smaller will not have sufficient energy to generate
additional electrons in the gas. This is an idealised scenario as in a real
detector, there will be an imperfect efficiency to detect single electrons.
For instance, the test run of the NEWS-G S140 detector achieved a 50%
detection efficiency for single electrons [138]. This calculation also assumes
that the Geant4 background model is sufficiently accurate all the way to 1 eV
and includes all sources that contribute to the single-electron background. As
discussed in section III, this is also an idealised scenario and therefore the
results in the upper panel of fig. 7 should be understood as an indication of
the absolute best-case sensitivity given the DarkSphere detector and operating
parameters.

Figure 7:  Projected 90% C.L. exclusion limits for $F_{\rm{DM}}=1$ when pure
mixtures of only helium, neon, xenon, methane, or isobutane are used. We have
assumed the background rate shown in fig. 4 together with the baseline
DarkSphere parameters: an SPC with a 3 m diameter filled with the target gas
at a pressure of 5 bar, temperature of $20^{\circ}$C, and taking data for 300
days. The upper panel assumes a threshold energy $E_{\rm{th}}=1$ eV, chosen to
mimic a single-electron threshold, while the bottom panel assumes
$E_{\rm{th}}=28$ eV, to mimic a two-electron threshold. Neon performs well at
higher DM masses but is worse at lower DM masses where xenon, methane and
isobutane benefit from lower ionisation energies. Grey bands show current
exclusion limits, so we see that DarkSphere has the potential to exclude
unexplored parameter space with any of the five gas targets shown.

Arguably a more realistic scenario is shown in the lower panel in fig. 7,
which assumes a threshold energy $E_{\rm{th}}=28$ eV. This energy is chosen to
mimic a detection threshold of at least two electrons. As discussed in section
III, the test run of the S140 experiment found a large excess of events in the
single-electron channel, with an event rate orders of magnitude larger than in
the two-electron channel, so the two-electron channel may be preferable for
physics searches [138]. The justification for $28$ eV as a proxy for a two-
electron threshold follows from a consideration of the $W$-values. Recall that
the $W$-value is the mean energy to create a new electron-ion pair and its
value depends on the gas under consideration.555The $W$-value is the mean
energy and not the minimum energy to create a new electron-ion pair. The
minimum energy is lower than $W$ [139], and is around the value of the
ionisation energies in tab. 1. The asymptotic value for xenon is $W\sim 22$ eV
[139], for methane and isobutane is $W\sim 28$ eV [139], and for neon is
$W\sim 37$ eV [140], while for helium, it is slightly higher at $W\sim 48$ eV
[141].666The difference between the $W$-value at approximately $30$ eV and the
asymptotic value is $\lesssim 10\%$ for a wide range of gases (including
methane) [142].  Although in fig. 7 we show results for a single-species gas,
in reality, the noble gases would always be mixed with a small component of
methane or isobutane to act as quench gases. For gas mixtures, the $W$-value
for the mixture is much closer to the lowest $W$-value of the gases in the
mixture [140, 143]. For example, $W\approx 28$ eV for a
$\mathrm{Ne}+\mathrm{CH}_{4}\,(2\%)$ mixture [71]. Therefore, we expect that
at or above approximately $28$ eV for all gases the ionised electron has
enough energy to be above the mean energy to create a new electron-ion pair.
In the lower panel, we again assume the idealised scenario where the
efficiency factor is 100% above the threshold energy.

In the calculation of the 90% C.L. limits, we use a binned-likelihood approach
and employ the exclusion-limit test-statistic and ‘Asimov’ data approach from
ref. [144]. We bin the data into equally sized bins of width 14 eV from the
threshold energy up to 99 (98) eV for $E_{\rm{th}}=1\,(28)$ eV. We use an
energy resolution of the form
$\sigma(E^{\rm{eV}}_{e})/E^{\rm{eV}}_{e}=\alpha+\beta/\sqrt{E^{\rm{eV}}_{e}}$,
where $\sigma(E^{\rm{eV}}_{e})$ and $E^{\rm{eV}}_{e}$ have units of eV, and
$\alpha=4.4\times 10^{-3}$ and $\beta=5.28$ are estimated by fitting to the
270 and 2822 eV lines in [71]. In the case of $E_{\rm{th}}=28$ eV, since
$\sigma(E^{\rm{eV}}_{e}=28\,\mathrm{eV})\simeq 28\,\mathrm{eV}$, to prevent
upward fluctuations biasing the limits, we multiply the rate by the rapidly
falling error function,
$\tfrac{1}{2}(1+\erf{[(E^{\rm{eV}}_{e}-28\,\mathrm{eV})/(\sqrt{2}\cdot
3\,\mathrm{eV})}])$, which only allows upward fluctuations from a small energy
range below the threshold. Finally, we also include a 10% uncertainty in the
overall normalisation of the background rate.

At DM masses above 100 MeV in both panels of fig. 7, we see that neon and
xenon constrain smaller values of the cross section compared to helium,
methane, and isobutane. At smaller masses, gases that have the lowest
ionisation energies perform best, so around 10 MeV, isobutane and methane
perform better than neon. Comparing the upper and lower panels, we see that
the effect of the higher threshold is an overall shift to higher cross
sections by a factor $\sim 3$ for xenon, $\sim 7$ for neon, $\sim 18$ for
helium, $\sim 27$ for methane, and $\sim 49$ for isobutane, accompanied by a
more severe weakening in the limits for DM masses below around 20 MeV. The
grey regions in fig. 7 show the current constraints from PandaX-II [145],
SENSEI [60], XENON10 [53], and XENON1T [56], so we see that, even for the
higher threshold scenario, the projected limits from all gases exclude regions
of unexplored parameter space.

Figure 8:  Similar to fig. 7 but with $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$. When
$E_{\rm{th}}=1$ eV, we see that methane and isobutane can outperform the noble
atom targets but suffer a larger suppression factor when moving to
$E_{\rm{th}}=28$ eV. Grey bands show current exclusion limits, so we again see
that DarkSphere has the potential to exclude regions of unexplored parameter
space with any of the five gas-targets shown.

In fig. 8, we show the analogous 90% C.L. limits for $F_{\rm{DM}}=(\alpha
m_{e}/q)^{2}$. Again, we show projections for $E_{\rm{th}}=1$ eV (top panel)
and $E_{\rm{th}}=28$ eV (bottom panel) to mimic a single-electron and multi-
electron search. The same assumptions about the background and resolution and
the same binned-likelihood method as used in the $F_{\rm{DM}}=1$ calculation
are used to generate these limits. In the upper panel of fig. 8, we see that
neon, xenon, methane, and isobutane all have a similar performance for DM
masses above 100 MeV, while helium is higher by a factor $\sim 10$. Comparing
the upper and lower panels, we find a more dramatic reduction in the
sensitivity relative to the $F_{\rm{DM}}=1$ scenario. Here, the shift to
higher cross sections is a factor $\sim 100$ for xenon, neon, and helium and a
factor $\sim 1000$ for methane and isobutane. The larger reduction in
sensitivity follows from the more rapid drop-off in the scattering rate with
energy when $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$, as was highlighted previously
in the discussion surrounding figs. 5 and 6. The grey-shaded regions show the
current constraints from SENSEI [60] and XENON10 [53] on this parameter space,
and we again observe that the projected limits from all gases exclude regions
of unexplored parameter space.

###  IV.2 Mixed-species projections

Figure 9:  Projected 90% C.L. exclusion limits for gases consisting of pure-
helium (solid orange) and pure-neon (solid blue), as well as the more
realistic scenario of gas mixtures: 90% helium and 10% isobutane in dashed
orange and 90% neon and 10% methane in dashed blue. The coloured shading
indicates the exclusions limits for a fraction of the molecular gases between
0% and 10%. In both panels, we assume the baseline DarkSphere parameters and
$E_{\rm{th}}=1$ eV. The upper panel has $F_{\rm{DM}}=1$, while the lower panel
has $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$. Adding a 10% component of isobutane to
helium significantly improves the sensitivity over helium alone, especially in
the lower panel. Adding a 10% component of methane to neon has a smaller
impact, although still leads to a small improvement at low DM masses.

In this final section, we show projections involving mixtures of gases. As we
have stated previously, this is the most likely mode of operation in a real
experiment as the noble gases need to run with some amount of methane or
isobutane, since they act as quenchers within the proportional counter [74].
Although it is in principle possible to run with a pure mixture of methane or
isobutane as they can act as both the ‘fill’ and quench gas, as both methane
and isobutane are flammable, in reality a pure mixture of either is unlikely
to pass the safety constraints in an underground laboratory.

In fig. 9, we show the projected 90% C.L. exclusion limits for helium and neon
(solid orange and blue lines) as well as for two gas mixtures: 90% helium and
10% isobutane (‘$\mathrm{He}+10\%\mathrm{C}_{4}\mathrm{H}_{10}$’) as dashed
orange; and 90% neon and 10% methane
(‘$\mathrm{Ne}+10\%\mathrm{C}\mathrm{H}_{4}$’) as dashed blue. A smaller
percentage of methane or isobutane in the mixture will lie in the shaded
regions between the dashed and solid lines. We have chosen these mixtures as
it is anticipated that this neon–methane mixture will be used in S140 and
ECuME, while the initial DarkSphere studies have assumed this helium–isobutane
mixture. These mixtures imply a total mass of 57.3 kg for neon–methane and
27.3 kg for helium–isobutane in the DarkSphere spherical chamber.

The upper panel in fig. 9 shows the projected limits for $F_{\rm{DM}}=1$,
while the lower panel shows $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$. For
$F_{\rm{DM}}=1$, the benefit of adding isobutane to helium starts to become
apparent for DM masses below about 30 MeV. Above this mass, the pure-helium
and helium–isobutane mixture lines are almost indistinguishable. The benefit
of adding methane to neon is much smaller, with only a noticeable change in a
relatively narrow window of DM masses between approximately 4 and 8 MeV. At
higher DM masses, the difference between neon and a neon–methane mixture is
almost indistinguishable. The situation changes for $F_{\rm{DM}}=(\alpha
m_{e}/q)^{2}$. In this case, the sizeable increase in sensitivity of the
helium–isobutane mixture compared to pure-helium is apparent over the full
range of DM masses. Even for neon, we see in the lower panel that there is
more of an advantage of running with a neon–methane mixture for DM masses
below approximately 20 MeV.

The limits in fig. 9 both assume $E_{\rm{th}}=1$ eV, however, a similar
pattern is also observed for gas mixtures with $E_{\rm{th}}=28$ eV (not
shown): we again find there is an advantage to the helium–isobutane mixture
over a pure-helium gas, while for neon, the benefit of adding methane is
smaller.

For clarity, we have only shown two gas mixtures. However, a helium–methane
mixture will be similar to a helium–isobutane mixture, while a neon–methane
mixture will be similar to a neon–isobutane mixture. This is for two reasons.
Firstly, the benefits of methane and isobutane arise largely because of their
lower ionisation energies relative to helium and neon, and importantly, the
smallest ionisation energies for both methane and isobutane are relatively
similar (cf. tab. 1). Secondly, we have demonstrated explicitly that methane
and isobutane have a similar level of sensitivity in the analysis of pure gas
targets (cf. figs. 7 and 8).

With these considerations in mind, we can also determine that the benefit of
adding methane or isobutane to xenon will be even smaller than for neon. At
higher DM masses, neon and xenon behave similarly (cf. figs. 7 and 8), so, as
was the case with neon, we find little increase in sensitivity to a xenon
mixture compared to pure xenon at higher DM masses. Meanwhile, at lower DM
masses, xenon already benefits from a relatively low ionisation energy that is
comparable to the methane and isobutane ionisation energies. These
considerations are borne out by explicit calculations (not shown) in which the
90% xenon and 10% methane or isobutane projected limits are almost
indistinguishable from the pure-xenon projected limit. Therefore, for
constraining DM-electron scattering, there is no meaningful increase in
sensitivity when adding methane or isobutane to xenon (although they are still
useful experimentally to act as quenchers in the proportional counter).

##  V Summary and outlook

The NEWS-G Collaboration employs a large spherical proportional counter (SPC)
filled with atomic or molecular gases to search for dark matter (DM). It has
been demonstrated that this detector technology has a single-ionisation
electron threshold, which allows it to set competitive constraints on DM-
nucleon interactions in the few-GeV mass range. However, a single ionisation
electron threshold should allow for the possibility of constraining DM-
electron interactions (cf. fig. 1). In this paper, we investigated this
possibility in the context of DarkSphere, a proposal for a 300 cm diameter SPC
with a fully electroformed-underground copper shell located in the Boulby
Underground Laboratory. Proportional counters are most commonly operated with
a mixture of gases. Therefore, we considered several noble gases – helium,
neon and xenon – plus methane and isobutane, all of which have been proposed
as the fill and quench gases for operation in the SPC.

While several independent DM-electron event rate calculations involving xenon
atoms have been presented in the literature, new calculations for helium,
neon, methane, and isobutane were required. The crucial quantity for
calculating the event rate is the dimensionless ionisation form factor, which
is calculated from the bound- and continuum-electron wave functions. We made
use of the PySCF quantum chemistry package to calculate the bound-state wave
functions, and the output was summarised in tab. 1 and fig. 2. To find the
continuum-electron wave functions, we used the HX-method for the noble atoms,
while for methane and isobutane, we employed the simpler and less accurate
(analytic) approximation of a Coulomb wave with effective charges determined
using a hydrogenic scaling of the ionisation energies. The resulting
dimensionless ionisation form factors were shown in fig. 3, and the DM-
electron scattering rates were given in figs. 4 and 5.

Projected exclusion limits for the DarkSphere proposal, presented in figs. 7,
8, and 9, show that this detector technology has the potential to constrain
large swathes of new parameter space. Existing constraints on the DM-electron
cross section for DM masses above 4 MeV can be improved by up to five orders
of magnitude in the case of a low-background, single electron search.

We found that neon is the best all-round gas target although it suffers
somewhat at low DM mass relative to xenon because of the comparatively large
outer-shell ionisation energy. In this context, adding a 10% component of
methane or isobutane to neon helps somewhat, as the lower ionisation energies
of methane and isobutane allow the gas-mixture to match xenon’s sensitivity at
low DM mass. In contrast, we find that adding methane or isobutane to helium
can result in an increase in sensitivity. Helium has the largest ionisation
energy of any atom, and with only two electrons, it has the lowest sensitivity
of the five gases that we considered. Adding a 10% component of methane or
isobutane to helium can potentially lead to a significant improvements in the
sensitivity, especially for the $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ scenario
(cf. fig. 9).

Given their ability to aid in the search for DM-electron interactions, further
refinements of the methane and isobutane calculations are warranted. Our
approach of treating the continuum-electron as a single-centred Coulomb wave
could be refined, for instance, by solving the Schroedinger equation in the
spherically averaged potential of the molecule, along the lines of the method
discussed in ref. [128] or by rewriting the single particle states and
electron density with an auxiliary single-centre basis set that could also
incorporate non-spherical contributions.

Our discussion has been centred around DarkSphere, the largest of the proposed
SPC experiments. We end by speculating on the sensitivity to DM-electron
scattering that may be achievable with S140, an experiment that uses a smaller
diameter SPC and where the background rate is higher. The S140 detector is
currently being commissioned at SNOLAB, and the background rate is expected to
be $\mathcal{O}(1)$ cts/kg/keV/day [135]. Using figs. 4 and 5, we can estimate
the sensitivity by finding the cross section when the background and signal
rates are comparable. For $F_{\rm{DM}}=1$, we therefore anticipate that S140
could achieve sensitivity to $\bar{\sigma}_{e}\sim 10^{-40}$ or $10^{-41}$
cm2, while for $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$, $\bar{\sigma}_{e}\sim
10^{-36}$ or $10^{-37}$ cm2 may be possible, depending on the precise value of
the detection efficiency that can be achieved. For both forms of
$F_{\rm{DM}}$, this would be an improvement on current exclusion limits for DM
heavier than approximately 10 MeV. This underlines the currently untapped
potential of SPCs to search for DM-electron scattering by sub-GeV DM, and we
advocate that this potential should be explored in greater depth.

Data Access Statement: The data supporting the findings reported in this paper
are openly available from the GitHub repository in [77].

###### Acknowledgements.

We are particularly indebted to Kostas Nikolopoulos, Ioannis Katsioulas,
Patrick Knights, and Jack Matthews for discussions about the DarkSphere
detector and for providing the simulated background rates in the Boulby
Underground Laboratory; George Booth and Robert Anderson for help with running
PySCF and for lending us books on quantum chemistry; Peter Cox, Matthew Dolan,
and Harry Quiney for discussions on atomic and molecular physics; and to
Thomas Edwards for comments on the manuscript. L.H. is supported by the
Cromwell Scholarship at King’s College London. C.M. is supported by the
Science and Technology Facilities Council (STFC) Grants No. ST/N004663/1 and
No. ST/T000759/1. For the purpose of open access, the authors have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising from this submission.

##  Appendix A Expressions for the ionisation form factor

In section II, we gave the general form for the rotationally averaged
ionisation form factor for a single electron. To recap, the general expression
for the form factor is

$\begin{split}\left|f^{i\to f}_{\rm{ion}}(E_{e},\mathbf{q})\right|^{2}&=\int
d\Omega_{k_{e}}\frac{2k_{e}^{3}}{8\pi^{3}}\\\ &\times\left|\int
d^{3}x\,\psi^{*}_{f}(\mathbf{x},\mathbf{k_{e}})\,e^{i\mathbf{q}\cdot\mathbf{x}}\,\psi_{i}(\mathbf{x})\right|^{2},\end{split}$
(14)

where $k_{e}=|\mathbf{k_{e}}|=\sqrt{2m_{e}E_{e}}$ is the momentum of the
continuum-electron and $d\Omega_{k_{e}}$ is the integration over the direction
of the outgoing electron. Note that this expression depends on $\mathbf{q}$.
However, it does not take into account that the atoms or molecules in a gas
will be found at an arbitrary and constantly changing orientation. In our
results, we therefore average over all possible orientations of the atom of
molecule. By averaging uniformly over all orientations encoded by the Euler
angles $\alpha$, $\beta$ and $\gamma$,

$\displaystyle\left|f^{i\to f}_{\rm{ion}}(E_{e},q)\right|^{2}$
$\displaystyle=\left\langle\left|f^{i\to
f}_{\rm{ion}}(E_{e},\mathbf{q})\right|^{2}\right\rangle$ (15)
$\displaystyle\begin{split}&=\frac{1}{8\pi^{2}}\int_{-1}^{1}\\!\\!d\cos\alpha\int_{0}^{2\pi}\\!\\!d\beta\int_{0}^{2\pi}\\!\\!d\gamma\\\
&\times\left|f^{i\to
f}_{\rm{ion}}(E_{e},\mathbf{q},\alpha,\beta,\gamma)\right|^{2}\;.\end{split}$
(16)

As a result of this averaging procedure, the ionisation form factor now only
depends on $q$, the magnitude of $\mathbf{q}$.

Although eqs. (14) and (15) are the most general expressions, in practice
there are more convenient expressions that can be used in calculations. In
this appendix, we provide derivations of the more convenient expressions.

###  A.1 Ionisation form factor with a plane wave

Although we only consider scattering rates where the outgoing electron is a
plane wave in appendix B, it is a useful starting point as the expressions are
the most straightforward to derive. We begin with the function

$\psi_{f}(\mathbf{x},\mathbf{k_{e}})=e^{i\mathbf{k}_{e}\cdot\mathbf{x}}\;,$
(17)

which trivially satisfies our normalisation convention for continuum wave
functions: $\int
d^{3}x\,\psi^{*}_{f}(\mathbf{x},\mathbf{k})\psi_{f}(\mathbf{x},\mathbf{k^{\prime}})=(2\pi)^{3}\delta^{3}(\mathbf{k}-\mathbf{k^{\prime}})$.
Substituting this into eq. (14), we recognise that we can express the form
factor in terms of the momentum-space wave function of the initial state,

$\left|f^{i\to f}_{\rm{ion}}(E_{e},\mathbf{q})\right|^{2}=\int
d\Omega_{k_{e}}\frac{2k_{e}^{3}}{8\pi^{3}}\
\left|\tilde{\psi}_{i}(\mathbf{k}_{e}-\mathbf{q})\right|^{2},$ (18)

where, to be explicit about our convention,

$\tilde{\psi}(\mathbf{p})=\int
d^{3}x\,e^{-i\mathbf{p}\cdot\mathbf{x}}\psi_{i}(\mathbf{x})\;.$ (19)

In practice, it more convenient to rewrite this by first noting that

$\left|\tilde{\psi}_{i}(\mathbf{k}_{e}-\mathbf{q})\right|^{2}=\int
d^{3}p\;\delta^{3}(\mathbf{p}-(\mathbf{k}_{e}-\mathbf{q}))\left|\tilde{\psi}_{i}(\mathbf{p})\right|^{2}\;.$
(20)

and then integrating the delta-function over $d\Omega_{k_{e}}$ and
$d\cos\theta_{p}$ to obtain

$\begin{split}\left|f^{i\to
f}_{\rm{ion}}(E_{e},\mathbf{q})\right|^{2}&=\frac{1}{4\pi^{3}}\frac{k_{e}^{2}}{q}\int_{0}^{2\pi}\\!\\!d\phi_{p}\int_{p^{-}}^{p^{+}}\\!\\!p\,dp\left|\tilde{\psi}_{i}(p,\phi_{p})\right|^{2}\end{split}.$
(21)

In this expression, we have aligned $\mathbf{q}$ in the $z$-direction and
$\cos\theta_{p}$ is the angle between $\mathbf{p}$ and $\mathbf{q}$. The
integration over the delta-function enforces the relation
$\cos\theta_{p}=(k_{e}^{2}-p^{2}-q^{2})/(2pq)$ and sets the limits of
integration, $p^{\pm}=|k_{e}\pm q|$. For the rotationally averaged molecular
plane wave form factors shown in appendix B, we calculate

$\begin{split}\left|f^{i\to
f}_{\rm{ion}}(E_{e},q)\right|^{2}&=\frac{1}{8\pi^{2}}\int_{-1}^{1}\\!\\!d\cos\alpha\int_{0}^{2\pi}\\!\\!d\beta\int_{0}^{2\pi}\\!\\!d\gamma\\\
&\times\frac{1}{4\pi^{3}}\frac{k_{e}^{2}}{q}\int_{0}^{2\pi}\\!\\!d\phi_{p}\int_{p^{-}}^{p^{+}}\\!\\!p\,dp\\\
&\times\left|\tilde{\psi}_{i}(p,\phi_{p},\alpha,\beta,\gamma)\right|^{2}\;,\end{split}$
(22)

where the integration over the Euler angles is carried out using Monte Carlo
techniques. As a result of this averaging procedure, the ionisation form
factor now only depends on $q$, the magnitude of $\mathbf{q}$.

For full shell atoms, averaging over the $m$ quantum number is equivalent to
rotational averaging. Therefore before proceeding, we pause to check that we
recover the result in ref. [38] for a spherically symmetric atom with full
shells. Writing
$\tilde{\psi}_{nlm}(\mathbf{p})=\chi_{nl}(p)Y_{lm}(\theta_{p},\phi_{p})$,
where $\chi_{nl}$ is the radial part of the momentum wave function and
$Y_{lm}(\theta_{p},\phi_{p})$ is the spherical harmonic, we find, after
summing over the spin and magnetic quantum numbers,

$\displaystyle\begin{split}\sum_{ms}\left|f^{i\to
f}_{\rm{ion}}(E_{e},\mathbf{q})\right|^{2}&=\frac{2}{4\pi^{3}}\frac{k_{e}^{2}}{q}\int_{p^{-}}^{p^{+}}\\!\\!p\,dp\left|\chi_{nl}(p)\right|^{2}\\\
&\times\int_{0}^{2\pi}\\!\\!d\phi_{p}\sum_{m}\left|Y_{lm}(\theta_{p},\phi_{p})\right|^{2}\end{split}$
(23)
$\displaystyle=\frac{2l+1}{4\pi^{3}}\frac{k_{e}^{2}}{q}\int_{p^{-}}^{p^{+}}\\!\\!p\,dp\left|\chi_{nl}(p)\right|^{2},$
(24)

which is the formula quoted in ref. [38].

###  A.2 Ionisation form factor with a Coulomb wave

For the Coulomb wave result, we start by writing the partial wave expansion in
the form

$\psi_{f}(\mathbf{x},\mathbf{k_{e}})=\sum_{l}\sum_{m=-l}^{l}\psi_{k_{e}lm}(\mathbf{x})Y_{lm}(\hat{\mathbf{k}}_{e})\;.$
(25)

Substituting this form into eq. (14) and using the key result

$\int
d\Omega_{k_{e}}Y^{*}_{lm}(\hat{\mathbf{k}}_{e})Y_{l^{\prime}m^{\prime}}(\hat{\mathbf{k}}_{e})=\delta_{ll^{\prime}}\delta_{mm^{\prime}}$
(26)

allows us to move the sum over $l$ and $m$ outside the modulus squared, and we
obtain

$\left|f^{i\to
f}_{\rm{ion}}(E_{e},\mathbf{q})\right|^{2}=\frac{2k_{e}^{3}}{8\pi^{3}}\sum_{lm}\left|\int
d^{3}x\,\psi^{*}_{k_{e}lm}(\mathbf{x})e^{i\mathbf{q}\cdot\mathbf{x}}\psi_{i}(\mathbf{x})\right|^{2}.$
(27)

Concretely then, for the Coulomb wave, we use the expression

$\psi_{k_{e}lm}(\mathbf{x})=\frac{P_{kl}(r)}{r}Y^{*}_{lm}(\hat{\mathbf{x}})$
(28)

where $P_{kl}(r)$ is given in eq. (13). Any phase factors (including the
$i^{l}$ term) are removed by the absolute magnitude squared in eq. (27). In
our calculations, we sum over $l$ and $m$ until the result converges to a few
percent accuracy. Typically, this occurs for $l\lesssim 5$; however, in some
cases, we need to sum up to $l\sim 15$. The rotationally averaged result is
again obtained by averaging over the Euler angles, which we do using Monte
Carlo integration.

In passing, we note that eq. (27) can also be used to calculate the form
factor for plane waves and this approach bypasses the need to calculate the
momentum-space wave function. From the plane wave expansion, we have that

$\psi_{k_{e}lm}(\mathbf{x})=4\pi j_{l}(k_{e}r)Y^{*}_{lm}(\hat{\mathbf{x}})\;,$
(29)

where $j_{l}(kr)$ is the spherical Bessel function. In our calculations
however, we found it easier to use the expressions in appendix A.1, especially
at large values of $q$.

###  A.3 Ionisation form factor for atoms

The expressions in eqs. (21) and (27) make no assumption about the symmetry of
the bound-state wave function. For atomic systems, we can exploit the
spherical symmetry of the atom to arrive at another compact expression that is
in use in the literature. We start by substituting the expressions

$\displaystyle\psi_{i}(\mathbf{x})$
$\displaystyle=\frac{P_{nl}(r)}{r}Y_{lm}(\hat{\mathbf{x}})$ (30)
$\displaystyle\psi^{*}_{k_{e}l^{\prime}m^{\prime}}(\mathbf{x})$
$\displaystyle=\frac{P_{k_{e}l^{\prime}}(r)}{r}Y_{l^{\prime}m^{\prime}}(\hat{\mathbf{x}})$
(31) $\displaystyle e^{i\mathbf{q}\cdot\mathbf{x}}$
$\displaystyle=\sum_{LM}4\pi
i^{L}j_{L}(qr)Y_{LM}(\hat{\mathbf{x}})Y^{*}_{LM}(\hat{\mathbf{q}})$ (32)

into eq. (27). Note that here we have assumed that $P_{k_{e}l^{\prime}}(r)$ is
real. We can simplify the resulting expression by using various identities of
the spherical harmonics,

$\displaystyle\begin{split}&\int
d\Omega_{x}\,Y_{lm}(\hat{\mathbf{x}})Y_{l^{\prime}m^{\prime}}(\hat{\mathbf{x}})Y_{LM}(\hat{\mathbf{x}})\\\
&=\sqrt{\frac{(2l+1)(2l^{\prime}+1)(2L+1)}{4\pi}}\\\
&\qquad\times\begin{pmatrix}l&l^{\prime}&L\\\
0&0&0\end{pmatrix}\begin{pmatrix}l&l^{\prime}&L\\\
m&m^{\prime}&M\end{pmatrix}\;,\end{split}$ (33)
$\displaystyle\sum_{M}Y_{LM}(\hat{\mathbf{q}})\,Y^{*}_{LM}(\hat{\mathbf{q}})=\frac{2L+1}{4\pi}\;,$
(34)

where the term in brackets is the Wigner 3-j symbol, which satisfies

$\sum_{mm^{\prime}}\begin{pmatrix}l&l^{\prime}&L\\\
m&m^{\prime}&M\end{pmatrix}\begin{pmatrix}l&l^{\prime}&L^{\prime}\\\
m&m^{\prime}&M^{\prime}\end{pmatrix}=\frac{\delta_{LL^{\prime}}\delta_{MM^{\prime}}}{2L+1}\;.$
(35)

Figure 10:  Differential event rates for DM-induced ionisation from methane
(left panels) and isobutane (right panels) for various continuum-state
approximations: Coulomb wave (solid blue), plane wave (dot-dashed teal), and
plane wave with the Fermi function using the same $Z_{\rm{eff}}$ values that
enter the Coulomb calculation (pink dashed). The DM mass is fixed to 30 MeV.
The top panels show $F_{\rm{DM}}=1$, while the bottom panels show
$F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$. Note the change in the scattering rate
scale between the top and bottom panels.

Pulling all of these results together, we find that the form factor averaged
over the initial-state $m$ quantum number is

$\displaystyle\left|f^{nl\\!\to\\!f}_{\mathrm{ion}}(E_{e},q)\right|^{2}\\!$
$\displaystyle=\frac{1}{2l+1}\sum_{m=-l}^{l}\left|f^{nl\to
f}_{\rm{ion}}(E_{e},\mathbf{q})\right|^{2}$ (36)
$\displaystyle\begin{split}&=\frac{k_{e}^{3}}{4\pi^{3}}\sum_{l^{\prime}L}\left[\\!\begin{pmatrix}l&l^{\prime}&L\\\
0&0&0\end{pmatrix}\\!\right]^{2}(2l^{\prime}+1)(2L+1)\\\ &\times\left|\int
drP_{k_{e}l^{\prime}}(r)P_{nl}(r)j_{L}(qr)\right|^{2}.\end{split}$ (37)

This matches the result in ref. [52] after multiplying by the occupation
number $w_{nl}=2(2l+1)$ since ref. [52] includes $w_{nl}$ in their definition
of the form factor. In our atomic physics calculations, we use this result and
sum both $l^{\prime}$ and $L$ in eq. (LABEL:eq:app:fionatom) up to 20, which
gives a result that converges to better than 2% accuracy for helium and neon,
and better than 3% for xenon. As the $P_{k_{e}l^{\prime}}(r)$ and $j_{L}(qr)$
terms can oscillate very rapidly for large values of $k_{e}$ and $q$, we are
careful to determine all of the real roots and supply these to our integration
algorithm to ensure accurate numerical results.

##  Appendix B Methane and isobutane rates: alternative approaches

In section II we calculated the methane and isobutane rates with the
assumption that the continuum-electron wave functions were Coulomb waves
centred on the centre-of-nuclear-charge of the molecule. We also made a
simplifying assumption for the charge $Z$, which is a free parameter in the
Coulomb function. In this appendix, we investigate three alternative
approximations used in the continuum-electron wave functions and consider
their impact on the scattering rate. Firstly, we consider a more basic
approximation where we treat the continuum-electron as a plane wave. Secondly,
as isobutane is a large molecule, we consider the impact of centring the
Coulomb waves at one of the carbon atoms at the base of the trigonal pyramid
rather than on the centre-of-nuclear-charge of the molecule. Finally, we
consider the impact of changing the values of the charge $Z$ in the Coulomb
function.

###  B.1 Plane wave approximation

In atomic calculations at low electron energies, the plane wave approximation
gives results that are less accurate than calculations with a Coulomb wave or
that solve the Schroedinger equation in the potential generated by the bound-
state electrons (see e.g., ref. [93]). Nonetheless, a comparison of the plane
wave results with the Coulomb wave results can give a crude estimate of the
accuracy of the calculations.

In atomic calculations, the plane wave calculations are improved somewhat by
including the Fermi function as a prefactor in eq. (2). The Fermi function is
defined as

$\displaystyle F(k_{e})$
$\displaystyle=\left(\frac{\psi^{\rm{Coulomb}}_{f}(\mathbf{x}\to
0,\mathbf{k_{e}})}{\psi^{\rm{Plane}}_{f}(\mathbf{x}\to
0,\mathbf{k_{e}})}\right)^{2}$ (38)
$\displaystyle=\frac{2\pi\eta(k_{e})}{1-\exp[-2\pi\eta(k_{e})]}\;,$ (39)

which is the square of the ratio of the Coulomb wave to the plane wave at the
origin [146], and where

$\eta(k_{e})=Z_{\rm{eff}}\frac{\alpha\,m_{e}}{k_{e}}$ (40)

is the Sommerfeld parameter (here given in natural units, while we gave the
Coulomb function in eq. (13) in atomic units).

As with the Coulomb wave calculation presented in section II, we centre the
plane wave on the centre-of-nuclear-charge of the molecule, and we must supply
a value of $Z_{\rm{eff}}$ that enters eq. (40). We follow the same procedure
used in section II, which leads to $Z_{\rm{eff}}$ values of $4.7$, $2.6$, and
$1.0$ for the $1a_{1}$, $2a_{1}$, and $1t_{2}$ states in methane,
respectively, and $Z_{\rm{eff}}=4.7$ for $1a_{1}$, $1e$, $2a_{1}$,
$Z_{\rm{eff}}=2.7$ for $3a_{1}$, $2e$, $4a_{1}$, and $Z_{\rm{eff}}=1$ for the
remaining outer valence orbitals in isobutane.

In fig. 10, we plot the methane (left panels) and isobutane (right panels)
scattering event rates in the Coulomb wave and plane wave approximations (with
and without the Fermi function). The results are shown for $F_{\rm{DM}}=1$
(top panels) and $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ (bottom panels) with
$m_{\rm{DM}}=30$ MeV. Similar relative differences are found between the rates
for DM masses in the range 5 to 1000 MeV.

In all cases, the Fermi function improves the agreement between the plane wave
and the Coulomb wave. For methane, the agreement over most of the energy range
between the Coulomb rate and plane wave rate with the Fermi function is at
about the level of 30%. Without the Fermi function, the Coulomb and plane wave
rates differ by more than an order of magnitude. The differences are larger
for isobutane, where there is a factor of approximately 2-to-3 difference
between the Coulomb and plane wave rate with the Fermi function. As the plane
wave with Fermi function rates are all larger than the Coulomb rates, the
corresponding projected exclusion limits that would be obtained using the
plane wave with Fermi function would be a factor of approximately two-to-three
lower than the limits shown in section IV.

Figure 11:  Comparison of the differential event rates for DM-induced
ionisation from isobutane when the continuum-electron Coulomb waves are
centred at the centre-of-nuclear-charge (solid dark-blue lines), and at a
carbon atom at the base of the trigonal pyramid (dashed pink lines). The
comparison is made for a DM mass of 30 and 1000 MeV. The top panel shows
$F_{\rm{DM}}=1$, while the bottom panel shows $F_{\rm{DM}}=(\alpha
m_{e}/q)^{2}$. Note the change in the scattering rate scale between the top
and bottom panels.

###  B.2 Alternative centre for the Coulomb waves

In section II the Coulomb waves were centred at the centre-of-nuclear-charge.
While in methane the centre-of-nuclear-charge corresponds to the position of
the carbon atom, in isobutane the centre-of-nuclear-charge is a point in space
between the carbon atoms. In the case of isobutane, it is therefore possible
that after ionisation the positive charge could be localised on a specific
carbon atom rather than the centre-of-nuclear-charge. To investigate the
impact of this possibility on the scattering rate prediction, we have
recalculated the ionisation form factors when the continuum-electron Coulomb
waves are centred on a carbon atom at the base of the trigonal pyramid.
Specifically, we centre the Coulomb waves on the carbon nucleus at
$(2.38,-1.37,-0.19)$ Bohr in our Cartesian coordinate system, as defined by
the NIST CCCBDB database [104].

In fig. 11, we show the differential scattering rate for two DM masses, 30 and
1000 MeV, and for $F_{\rm{DM}}=1$ (top) and $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$
(bottom). The solid lines show the result of the calculation in section II
while the dashed lines show the rates when the centre is the carbon atom at
the base of the trigonal pyramid. The difference between the two calculations
is 40% or smaller for both masses and both choices of $F_{\rm{DM}}$. Similar
results are found for all other DM masses considered in this work.

Figure 12:  Comparison of the differential event rates for DM-induced
ionisation from isobutane for different values of $Z$, which is a free
parameter in the Coulomb function. Solid blue shows the rates when common $Z$
values are assigned to the groups of inner core, inner valence and outer
valence orbitals; dashed pink shows the rates when
$Z=n\sqrt{I/13.6~{}\mathrm{eV}}$, where $I$ is the ionisation energy; and dot-
dashed teal is when $Z=2\times n\sqrt{I/13.6~{}\mathrm{eV}}$. The comparison
is made for a DM mass of 30 MeV and 1000 MeV. The top panel shows
$F_{\rm{DM}}=1$ while the bottom panel shows $F_{\rm{DM}}=(\alpha
m_{e}/q)^{2}$. Note the change in the scattering rate scale between the top
and bottom panels.
Figure 13:  Left panels: Comparison of the expected number of events from this
work, QEDark [147] and Pandey et al [76] for DM-electron scattering in xenon,
for $m_{\rm{DM}}=500$ MeV, a 1000 kg-years exposure, and for $F_{\rm{DM}}=1$
(top) and $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ (bottom). Right panel: An
orbital-by-orbital comparison of the expected number of events from this work
and QEDark for $m_{\rm{DM}}=500$ MeV and $F_{\rm{DM}}=1$.

###  B.3 Alternative values of the effective charge

The charge $Z$ remains a free parameter in the Coulomb function. We used the
values $Z=4.7$, $Z=2.6$, and $Z=1.0$ for the $1a_{1}$, $2a_{1}$, and $1t_{2}$
states in methane, which were obtained using the formula for the hydrogen-like
scaling for the ionisation energy. For isobutane we made a further simplifying
approximation, which was to assign $Z=4.7$ for the inner core orbitals
($1a_{1}$, $1e$, $2a_{1}$), $Z=2.7$ for the inner valence orbitals ($3a_{1}$,
$2e$, $4a_{1}$), and $Z=1$ for the outer valence states (all remaining
orbitals).

Here, we investigate two alternative choices for the $Z$ values that enter the
isobutane calculation. Firstly, we assign different $Z$ values to each
orbital, calculated using $I=Z^{2}/n^{2}\times 13.6~{}\mathrm{eV}$, where $I$
is the experimental ionisation energy (when available or theoretical value
otherwise) and $n=1$ for the inner core and outer valence states, and $n=2$
for the inner valence states. The resulting scattering rates for two DM
masses, 30 and 1000 MeV, are shown as the dashed pink lines in fig. 12. These
are compared with the rates from section II shown in solid blue, where common
$Z$ values were assigned to the groups of inner core, inner valence, and outer
valence orbitals. The rates differ by 25% or less for both DM masses shown and
for both choices of $F_{\rm{DM}}$. Secondly, we consider the (possibly
unphysical) scenario where the value of $Z$ is arbitrarily doubled, such that
$Z=2\times n\sqrt{I/13.6~{}\mathrm{eV}}$, in order to obtain intuition for how
the rates scale with $Z$. The resulting rates are shown as the dot-dashed teal
lines in fig. 12. Comparing with the dashed pink lines, we find that the
scattering rates where $Z$ has been doubled are between a factor of $1.5$ and
$2$ larger

##  Appendix C Comparison with other calculations

The literature contains several calculations of the DM-electron scattering
rate for xenon. In this appendix, we provide a comparison with the results
from QEDark (specifically, Xe_v2.0) [147], a commonly used package in limit-
setting calculations (see, e.g., [56]), and with Pandey et al [76], which
provides both a relativistic and non-relativistic calculation and reports some
deviations from the QEDark results.

In making the comparisons, we use the same astrophysical parameters employed
in those papers: $\rho_{\rm{DM}}=0.4\,\mathrm{GeV}/\mathrm{cm}^{3}$,
$v_{0}=220\,\mathrm{km/s}$, and $v_{\mathrm{esc}}=544\,\mathrm{km/s}$. We also
follow the same prescription used in those works to calculate the expected
number of electrons $n_{e}$ from a given $dR/dE_{e}$ recoil spectrum, namely,
we use the xenon binding energies and the input parameters
$W=13.8\,\mathrm{eV}$, $f_{e}=0.83$, and $f_{R}=0$ from ref. [53].

The solid blue lines in fig. 13 show the spectrum of the expected number of
events from our non-relativistic calculations for $F_{\rm{DM}}=1$ (top-left
panel) and $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$ (bottom-left panel). We have
assumed a DM mass of 500 MeV and a total exposure of 1000 kg-years. In the
same panels, we show the non-relativistic result from QEDark (Xe_v2.0),
together with the non-relativistic and relativistic calculations from Pandey
et al.

Comparing firstly with QEDark (pink dashed lines), we find good agreement for
a small number of electrons ($n\lesssim 4$) but a deviation at larger values.
This is most obvious when $F_{\rm{DM}}=1$ where the difference reaches a
factor of approximately three for $n=14$. The difference at larger $n_{e}$ is
less pronounced when $F_{\rm{DM}}=(\alpha m_{e}/q)^{2}$. To better understand
the difference in the $F_{\rm{DM}}=1$ case, in the right panel of fig. 13, we
have plotted the contributions from the individual orbitals for our results
and for QEDark. This clearly shows that the difference arises from the $4p$
shell. While our calculation shows that $4d$ and $4p$ both contribute at a
similar level at $n_{e}\sim 14$, the QEDark result for $4p$ is approximately
an order of magnitude smaller. There is also a substantial difference with
$4s$, but this orbital gives a sub-dominant contribution to the total rate.

Turning next to a comparison with Pandey et al (green dashed and orange dot-
dashed lines), we find agreement with the non-relativistic calculation at the
level of 30% or better for $n_{e}\gtrsim 3$. Both our result and the Pandey et
al non-relativistic calculation deviate from QEDark at larger values of
$n_{e}$. The relativistic calculation shows an enhancement over all non-
relativistic calculations at larger $n_{e}$, which is understood as arising
from the more accurate characterisation of the electron wave function at very
small distances in the relativistic calculation (see, e.g., ref. [92]).
Comparing the results at small electron number ($n_{e}<3$), we find that our
result and that of QEDark are about a factor of two higher than Pandey et al.
The $n_{e}=1$ and $n_{e}=2$ results (equivalent to low electron kinetic
energy) are most sensitive to the effective potential used in the calculation.
The Pandey et al calculation only includes the classical potential energy in
their calculation of the continuum wave function. In contrast, our result
includes the potential energy and an effective exchange potential. Rerunning
our calculations with $k_{x}=0$ in eq. (10) (so that only the classical
potential energy contributes), we find that the scattering rate falls just
below the Pandey et al result for $n_{e}<3$. This confirms that the $n_{e}=1$
and $n_{e}=2$ rates are sensitive to the effective potential, so differences
between our result and Pandey et al should be expected since we use different
models of the effective potential.

## References

  * [1] European Astroparticle Physics Strategy 2017-2026, https://www.appec.org/roadmap.
  * Bertone and Hooper [2018] G. Bertone and D. Hooper, Rev. Mod. Phys. 90, 045002 (2018), arXiv:1605.04909 [astro-ph.CO] .
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Edwards  _et al._ [2018] T. D. P. Edwards, B. J. Kavanagh, and C. Weniger, Phys. Rev. Lett. 121, 181101 (2018), arXiv:1805.04117 [hep-ph] .
  * Billard  _et al._ [2022] J. Billard _et al._ , Rept. Prog. Phys. 85, 056201 (2022), arXiv:2104.07634 [hep-ex] .
  * Baudis [2016] L. Baudis, Annalen Phys. 528, 74 (2016), arXiv:1509.00869 [astro-ph.CO] .
  * Marrodan Undagoitia and Rauch [2016] T. Marrodan Undagoitia and L. Rauch, J. Phys. G 43, 013001 (2016), arXiv:1509.08767 [physics.ins-det] .
  * Schumann [2019] M. Schumann, J. Phys. G 46, 103003 (2019), arXiv:1903.03026 [astro-ph.CO] .
  * Battaglieri  _et al._ [2017] M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter_ (2017) arXiv:1707.04591 [hep-ph] .
  * Arnaud  _et al._ [2018a] Q. Arnaud _et al._ (EDELWEISS), Phys. Rev. D 97, 022003 (2018a), arXiv:1707.04308 [physics.ins-det] .
  * Agnese  _et al._ [2019] R. Agnese _et al._ (SuperCDMS), Phys. Rev. D 99, 062001 (2019), arXiv:1808.09098 [astro-ph.CO] .
  * Abdelhameed  _et al._ [2019a] A. H. Abdelhameed _et al._ (CRESST), Eur. Phys. J. C 79, 630 (2019a), arXiv:1902.07587 [astro-ph.IM] .
  * Abdelhameed  _et al._ [2019b] A. H. Abdelhameed _et al._ (CRESST), Phys. Rev. D 100, 102002 (2019b), arXiv:1904.00498 [astro-ph.CO] .
  * Alkhatib  _et al._ [2021] I. Alkhatib _et al._ (SuperCDMS), Phys. Rev. Lett. 127, 061801 (2021), arXiv:2007.14289 [hep-ex] .
  * Aguilar-Arevalo  _et al._ [2020] A. Aguilar-Arevalo _et al._ (DAMIC), Phys. Rev. Lett. 125, 241803 (2020), arXiv:2007.15622 [astro-ph.CO] .
  * Kouvaris and Pradler [2017] C. Kouvaris and J. Pradler, Phys. Rev. Lett. 118, 031803 (2017), arXiv:1607.01789 [hep-ph] .
  * McCabe [2017] C. McCabe, Phys. Rev. D 96, 043010 (2017), arXiv:1702.04730 [hep-ph] .
  * Ibe  _et al._ [2018] M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, JHEP 03, 194, arXiv:1707.07258 [hep-ph] .
  * Dolan  _et al._ [2018] M. J. Dolan, F. Kahlhoefer, and C. McCabe, Phys. Rev. Lett. 121, 101801 (2018), arXiv:1711.09906 [hep-ph] .
  * Akerib  _et al._ [2019] D. S. Akerib _et al._ (LUX), Phys. Rev. Lett. 122, 131301 (2019), arXiv:1811.11241 [astro-ph.CO] .
  * Armengaud  _et al._ [2019] E. Armengaud _et al._ (EDELWEISS), Phys. Rev. D 99, 082003 (2019), arXiv:1901.03588 [astro-ph.GA] .
  * Bell  _et al._ [2020] N. F. Bell, J. B. Dent, J. L. Newstead, S. Sabharwal, and T. J. Weiler, Phys. Rev. D 101, 015012 (2020), arXiv:1905.00046 [hep-ph] .
  * Liu  _et al._ [2019] Z. Z. Liu _et al._ (CDEX), Phys. Rev. Lett. 123, 161301 (2019), arXiv:1905.00354 [hep-ex] .
  * Baxter  _et al._ [2020] D. Baxter, Y. Kahn, and G. Krnjaic, Phys. Rev. D 101, 076014 (2020), arXiv:1908.00012 [hep-ph] .
  * Essig  _et al._ [2020] R. Essig, J. Pradler, M. Sholapurkar, and T.-T. Yu, Phys. Rev. Lett. 124, 021801 (2020), arXiv:1908.10881 [hep-ph] .
  * Liang  _et al._ [2020] Z.-L. Liang, L. Zhang, F. Zheng, and P. Zhang, Phys. Rev. D 102, 043007 (2020), arXiv:1912.13484 [cond-mat.mes-hall] .
  * Knapen  _et al._ [2021] S. Knapen, J. Kozaczuk, and T. Lin, Phys. Rev. Lett. 127, 081805 (2021), arXiv:2011.09496 [hep-ph] .
  * Araujo  _et al._ [2022] H. M. Araujo _et al._ , arXiv:2207.08284 [hep-ex] (2022).
  * Cox  _et al._ [2022] P. Cox, M. J. Dolan, C. McCabe, and H. M. Quiney, arXiv:2208.12222 [hep-ph] (2022).
  * Kouvaris [2015] C. Kouvaris, Phys. Rev. D 92, 075001 (2015), arXiv:1506.04316 [hep-ph] .
  * An  _et al._ [2018] H. An, M. Pospelov, J. Pradler, and A. Ritz, Phys. Rev. Lett. 120, 141801 (2018), [Erratum: Phys.Rev.Lett. 121, 259903 (2018)], arXiv:1708.03642 [hep-ph] .
  * Emken  _et al._ [2018] T. Emken, C. Kouvaris, and N. G. Nielsen, Phys. Rev. D 97, 063007 (2018), arXiv:1709.06573 [hep-ph] .
  * Bringmann and Pospelov [2019] T. Bringmann and M. Pospelov, Phys. Rev. Lett. 122, 171801 (2019), arXiv:1810.10543 [hep-ph] .
  * Ema  _et al._ [2019] Y. Ema, F. Sala, and R. Sato, Phys. Rev. Lett. 122, 181802 (2019), arXiv:1811.00520 [hep-ph] .
  * Alvey  _et al._ [2019] J. Alvey, M. Campos, M. Fairbairn, and T. You, Phys. Rev. Lett. 123, 261802 (2019), arXiv:1905.05776 [hep-ph] .
  * Cappiello and Beacom [2019] C. Cappiello and J. F. Beacom, Phys. Rev. D 100, 103011 (2019), arXiv:1906.11283 [hep-ph] .
  * Herrera and Ibarra [2021] G. Herrera and A. Ibarra, Phys. Lett. B 820, 136551 (2021), arXiv:2104.04445 [hep-ph] .
  * Essig  _et al._ [2012a] R. Essig, J. Mardon, and T. Volansky, Phys. Rev. D 85, 076007 (2012a), arXiv:1108.5383 [hep-ph] .
  * Graham  _et al._ [2012] P. W. Graham, D. E. Kaplan, S. Rajendran, and M. T. Walters, Phys. Dark Univ. 1, 32 (2012), arXiv:1203.2531 [hep-ph] .
  * Hochberg  _et al._ [2016a] Y. Hochberg, Y. Zhao, and K. M. Zurek, Phys. Rev. Lett. 116, 011301 (2016a), arXiv:1504.07237 [hep-ph] .
  * Essig  _et al._ [2016] R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto, T. Volansky, and T.-T. Yu, JHEP 05, 046, arXiv:1509.01598 [hep-ph] .
  * Hochberg  _et al._ [2016b] Y. Hochberg, M. Pyle, Y. Zhao, and K. M. Zurek, JHEP 08, 057, arXiv:1512.04533 [hep-ph] .
  * Derenzo  _et al._ [2017] S. Derenzo, R. Essig, A. Massari, A. Soto, and T.-T. Yu, Phys. Rev. D 96, 016026 (2017), arXiv:1607.01009 [hep-ph] .
  * Essig  _et al._ [2017a] R. Essig, J. Mardon, O. Slone, and T. Volansky, Phys. Rev. D 95, 056011 (2017a), arXiv:1608.02940 [hep-ph] .
  * Cavoto  _et al._ [2018] G. Cavoto, F. Luchetta, and A. D. Polosa, Phys. Lett. B 776, 338 (2018), arXiv:1706.02487 [hep-ph] .
  * Hochberg  _et al._ [2018] Y. Hochberg, Y. Kahn, M. Lisanti, K. M. Zurek, A. G. Grushin, R. Ilan, S. M. Griffin, Z.-F. Liu, S. F. Weber, and J. B. Neaton, Phys. Rev. D 97, 015004 (2018), arXiv:1708.08929 [hep-ph] .
  * Kurinsky  _et al._ [2019] N. A. Kurinsky, T. C. Yu, Y. Hochberg, and B. Cabrera, Phys. Rev. D 99, 123005 (2019), arXiv:1901.07569 [hep-ex] .
  * Essig  _et al._ [2019] R. Essig, J. Pérez-Ríos, H. Ramani, and O. Slone, Phys. Rev. Research. 1, 033105 (2019), arXiv:1907.07682 [hep-ph] .
  * Geilhufe  _et al._ [2020] R. M. Geilhufe, F. Kahlhoefer, and M. W. Winkler, Phys. Rev. D 101, 055005 (2020), arXiv:1910.02091 [hep-ph] .
  * Blanco  _et al._ [2021] C. Blanco, Y. Kahn, B. Lillard, and S. D. McDermott, Phys. Rev. D 104, 036011 (2021), arXiv:2103.08601 [hep-ph] .
  * Kahn and Lin [2022] Y. Kahn and T. Lin, Rept. Prog. Phys. 85, 066901 (2022), arXiv:2108.03239 [hep-ph] .
  * Essig  _et al._ [2012b] R. Essig, A. Manalaysay, J. Mardon, P. Sorensen, and T. Volansky, Phys. Rev. Lett. 109, 021301 (2012b), arXiv:1206.2644 [astro-ph.CO] .
  * Essig  _et al._ [2017b] R. Essig, T. Volansky, and T.-T. Yu, Phys. Rev. D 96, 043017 (2017b), arXiv:1703.00910 [hep-ph] .
  * Agnes  _et al._ [2018] P. Agnes _et al._ (DarkSide), Phys. Rev. Lett. 121, 111303 (2018), arXiv:1802.06998 [astro-ph.CO] .
  * Agnese  _et al._ [2018] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 121, 051301 (2018), [Erratum: Phys.Rev.Lett. 122, 069901 (2019)], arXiv:1804.10697 [hep-ex] .
  * Aprile  _et al._ [2019] E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 251801 (2019), arXiv:1907.11485 [hep-ex] .
  * Aguilar-Arevalo  _et al._ [2019] A. Aguilar-Arevalo _et al._ (DAMIC), Phys. Rev. Lett. 123, 181802 (2019), arXiv:1907.12628 [astro-ph.CO] .
  * Blanco  _et al._ [2020] C. Blanco, J. I. Collar, Y. Kahn, and B. Lillard, Phys. Rev. D 101, 056001 (2020), arXiv:1912.02822 [hep-ph] .
  * Arnaud  _et al._ [2020] Q. Arnaud _et al._ (EDELWEISS), Phys. Rev. Lett. 125, 141301 (2020), arXiv:2003.01046 [astro-ph.GA] .
  * Barak  _et al._ [2020] L. Barak _et al._ (SENSEI), Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO] .
  * Giomataris and Vergados [2004] Y. Giomataris and J. D. Vergados, Nucl. Instrum. Meth. A 530, 330 (2004), arXiv:hep-ex/0303045 .
  * Giomataris  _et al._ [2008] I. Giomataris _et al._ , Journal of Instrumentation 3, P09007 (2008).
  * Arnaud  _et al._ [2018b] Q. Arnaud _et al._ (NEWS-G), Astropart. Phys. 97, 54 (2018b), arXiv:1706.04934 [astro-ph.IM] .
  * Balogh  _et al._ [2021] L. Balogh _et al._ (NEWS-G), Nucl. Instrum. Meth. A 988, 164844 (2021), arXiv:2008.03153 [physics.ins-det] .
  * Balogh  _et al._ [2022] L. Balogh _et al._ (NEWS-G), arXiv:2205.15433 [physics.ins-det] (2022).
  * Balogh  _et al._ [2023] L. Balogh _et al._ , arXiv:2301.05183 [hep-ex] (2023).
  * Giganon  _et al._ [2017] A. Giganon _et al._ , JINST 12 (12), P12031, arXiv:1707.09254 [physics.ins-det] .
  * Giomataris  _et al._ [2020] I. Giomataris _et al._ , JINST 15 (11), 11, arXiv:2003.01068 [physics.ins-det] .
  * Katsioulas  _et al._ [2022a] I. Katsioulas, P. Knights, I. Manthos, J. Matthews, T. Neep, K. Nikolopoulos, and R. Ward, JINST 17 (08), C08025.
  * Bougamont  _et al._ [2012] E. Bougamont _et al._ , JMP 3, 57 (2012), arXiv:1010.4132 [physics.ins-det] .
  * Arnaud  _et al._ [2019] Q. Arnaud _et al._ (NEWS-G), Phys. Rev. D 99, 102003 (2019), arXiv:1902.08960 [physics.ins-det] .
  * Vergados and Giomataris [2007] J. D. Vergados and Y. Giomataris, Phys. Atom. Nucl. 70, 140 (2007), arXiv:astro-ph/0511470 .
  * Meregaglia  _et al._ [2018] A. Meregaglia _et al._ , JINST 13 (01), P01009, arXiv:1710.04536 [physics.ins-det] .
  * Knoll [2000] G. F. Knoll,  _Radiation detection and measurement, Third edition_ (John Wiley & Sons, Inc., 2000).
  * Roberts  _et al._ [2016a] B. Roberts, V. Dzuba, V. Flambaum, M. Pospelov, and Y. Stadnik, Phys. Rev. D 93, 115037 (2016a), arXiv:1604.04559 [hep-ph] .
  * Pandey  _et al._ [2020] M. K. Pandey _et al._ , Phys. Rev. D 102, 123025 (2020), arXiv:1812.11759 [hep-ph] .
  * [77] Atomic and molecular orbitals and form factors on GitHub, https://github.com/mccabech/AtMolDM.
  * Reid and Brunthaler [2004] M. J. Reid and A. Brunthaler, Astrophys. J. 616, 872 (2004), arXiv:astro-ph/0408107 .
  * Gravity Collaboration  _et al._ [2019] Gravity Collaboration, R. Abuter, _et al._ , Astronomy & Astrophysics 625, L10 (2019), arXiv:1904.05721 [astro-ph.GA] .
  * Abuter  _et al._ [2021] R. Abuter _et al._ (Gravity Collaboration), Astronomy & Astrophysics 647, A59 (2021), arXiv:2101.12098 [astro-ph.GA] .
  * Baxter  _et al._ [2021] D. Baxter _et al._ , Eur. Phys. J. C 81, 907 (2021), arXiv:2105.00599 [hep-ex] .
  * Deason  _et al._ [2019] A. J. Deason _et al._ , MNRAS 485, 3514 (2019), arXiv:1901.02016 [astro-ph.GA] .
  * Schoenrich  _et al._ [2010] R. Schoenrich, J. Binney, and W. Dehnen, Mon. Not. Roy. Astron. Soc. 403, 1829 (2010), arXiv:0912.3693 [astro-ph.GA] .
  * McCabe [2014] C. McCabe, JCAP 02, 027, arXiv:1312.1355 [astro-ph.CO] .
  * O’Hare  _et al._ [2018] C. A. O’Hare, C. McCabe, N. W. Evans, G. Myeong, and V. Belokurov, Phys. Rev. D 98, 103006 (2018), arXiv:1807.09004 [astro-ph.CO] .
  * Evans  _et al._ [2019] N. W. Evans, C. A. O’Hare, and C. McCabe, Phys. Rev. D 99, 023012 (2019), arXiv:1810.11468 [astro-ph.GA] .
  * O’Hare  _et al._ [2020] C. A. O’Hare, N. W. Evans, C. McCabe, G. Myeong, and V. Belokurov, Phys. Rev. D 101, 023006 (2020), arXiv:1909.04684 [astro-ph.GA] .
  * Buch  _et al._ [2020] J. Buch, M. A. Buen-Abad, J. Fan, and J. S. C. Leung, Phys. Rev. D 102, 083010 (2020), arXiv:2007.13750 [hep-ph] .
  * Radick  _et al._ [2021] A. Radick, A.-M. Taki, and T.-T. Yu, JCAP 02, 004, arXiv:2011.02493 [hep-ph] .
  * Sun  _et al._ [2017] Q. Sun _et al._ , WIREs Comput Mol Sci 8, e1340 (2017).
  * Sun  _et al._ [2020] Q. Sun _et al._ , The Journal of Chemical Physics 153, 024109 (2020).
  * Roberts  _et al._ [2016b] B. Roberts, V. Flambaum, and G. Gribakin, Phys. Rev. Lett. 116, 023201 (2016b), arXiv:1509.09044 [physics.atom-ph] .
  * Roberts and Flambaum [2019] B. Roberts and V. Flambaum, Phys. Rev. D 100, 063017 (2019), arXiv:1904.07127 [hep-ph] .
  * Kennedy and Manson [1972] D. J. Kennedy and S. T. Manson, Phys. Rev. A 5, 227 (1972).
  * Woon and Dunning [1994] D. E. Woon and T. H. Dunning, The Journal of Chemical Physics 100, 2975 (1994).
  * Dunning [1989] T. H. Dunning, The Journal of Chemical Physics 90, 1007 (1989).
  * Ceolin  _et al._ [2013] G. A. Ceolin, R. C. de Berrêdo, and F. E. Jorge, Theoretical Chemistry Accounts 132, 1339 (2013).
  * Ditchfield  _et al._ [1971] R. Ditchfield, W. J. Hehre, and J. A. Pople, J. Chem. Phys. 54, 724 (1971).
  * Hehre  _et al._ [1972] W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys. 56, 2257 (1972).
  * Hariharan and Pople [1973] P. C. Hariharan and J. A. Pople, Theor. Chim. Acta 28, 213 (1973).
  * Feller [1996] D. Feller, J. Comput. Chem. 17, 1571 (1996).
  * Schuchardt  _et al._ [2007] K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li, and T. L. Windus, J. Chem. Inf. Model. 47, 1045 (2007).
  * Pritchard  _et al._ [2019] B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibsom, and T. L. Windus, J. Chem. Inf. Model. 59, 4814 (2019).
  * NIST Computational Chemistry Comparison and Benchmark Database [ III] NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101 http://cccbdb.nist.gov/, 10.18434/T47C7Z (Release 21, August 2020, Editor: Russell D. Johnson III).
  * Thomas [1970] T. D. Thomas, The Journal of Chemical Physics 52, 1373 (1970).
  * Bieri and Asbrink [1980] G. Bieri and L. Asbrink, Journal of Electron Spectroscopy and Related Phenomena 20, 149 (1980).
  * Deng  _et al._ [2001] J. K. Deng _et al._ , The Journal of Chemical Physics 114, 882 (2001).
  * Williams [2009] G. Williams,  _X-Ray Data Booklet, Sec.(1.1)_ (Lawrence Berkeley National Laboratory, Lawrence Berkeley National Laboratory, Berkeley CA 94720, U.S.A., 2009).
  * Bunge  _et al._ [1993] C. Bunge, J. Barrientos, and A. Bunge, Atomic Data and Nuclear Data Tables 53, 113 (1993).
  * Fischer [2000] C. F. Fischer, Computer Physics Communications 128, 635 (2000).
  * Brundle  _et al._ [1970] C. R. Brundle, M. B. Robin, and H. Basch, The Journal of Chemical Physics 53, 2196 (1970).
  * Hamel  _et al._ [2002] S. Hamel, P. Duffy, M. E. Casida, and D. R. Salahub, Journal of Electron Spectroscopy and Related Phenomena 123, 345 (2002).
  * Grant [2007] I. P. Grant, ed.,  _Relativistic Quantum Theory of Atoms and Molecules_ (Springer New York, 2007).
  * Jonsson  _et al._ [2007] P. Jonsson, X. He, C. Froese Fischer, and I. P. Grant, Comput. Phys. Commun. 177, 597 (2007).
  * Froese Fischer  _et al._ [2019] C. Froese Fischer, G. Gaigalas, P. Jönsson, and J. Bieroń, Comput. Phys. Commun. 237, 184 (2019).
  * Quiney  _et al._ [1998] H. M. Quiney, H. Skaane, and I. P. Grant, Adv. Quantum Chem. 32, 1 (1998).
  * Wolfram Research, Inc. [2021] Wolfram Research, Inc., Mathematica, Version 12.3.1 (2021).
  * Bates and Seaton [1949] D. R. Bates and M. J. Seaton, Monthly Notices of the Royal Astronomical Society 109, 698 (1949).
  * Amusia and Chernysheva [1997] M. Y. Amusia and L. V. Chernysheva,  _Computation of Atomic Processes_ (IoP Publishing Ltd, 1997).
  * Cowan [1967] R. D. Cowan, Phys. Rev. 163, 54 (1967).
  * Cowan [1981] R. D. Cowan,  _The theory of atomic structure and spectra_ (University of California Press, 1981).
  * Condon and Odabasi [1980] E. Condon and H. Odabasi,  _Atomic Structure_ (Cambridge University Press, 1980).
  * Catena  _et al._ [2020] R. Catena, T. Emken, N. A. Spaldin, and W. Tarantino, Phys. Rev. Res. 2, 033195 (2020), arXiv:1912.08204 [hep-ph] .
  * Jain [1986] A. Jain, Phys. Rev. A 34, 3707 (1986).
  * Baluja  _et al._ [1992] K. L. Baluja, A. Jain, V. D. Martino, and F. A. Gianturco, Europhysics Letters (EPL) 17, 139 (1992).
  * Das  _et al._ [2014] T. Das, A. D. Stauffer, and R. Srivastava, The European Physical Journal D 68, 102 (2014).
  * Er-Jun  _et al._ [2007] M. Er-Jun, M. Yu-Gang, C. Xiang-Zhou, F. De-Qing, S. Wen-Qing, and T. Wen-Dong, Chinese Physics 16, 3339 (2007).
  * Mahato  _et al._ [2019] D. Mahato, L. Sharma, A. D. Stauffer, and R. Srivastava, The European Physical Journal D 73, 189 (2019).
  * Galassi  _et al._ [2021] M. Galassi _et al._ , GNU Scientific Library Reference Manual (3rd Ed.) (2021).
  * Johnson [2020] S. Johnson, Cubature v1.0.4 (2020).
  * Tange [2022] O. Tange, GNU Parallel 20220722 (‘Roe vs Wade’) (2022).
  * Kudryavtsev  _et al._ [2003] V. A. Kudryavtsev, N. J. C. Spooner, and J. E. McMillan, Nucl. Instrum. Meth. A 505, 688 (2003), arXiv:hep-ex/0303007 .
  * Smith  _et al._ [2005] P. F. Smith, D. Snowden-Ifft, N. J. T. Smith, R. Luscher, and J. D. Lewin, Astropart. Phys. 22, 409 (2005).
  * Malczewski  _et al._ [2013] D. Malczewski, J. Kisiel, and J. Dorda, J. Radioanal. Nucl. Chem. 298, 1483 (2013).
  * Brossard [2020] A. Brossard,  _Optimization of spherical proportional counter backgrounds and response for low mass dark matter search_ , Theses, Université Paris-Saclay; Queen’s University at Kingston (2020).
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Amare  _et al._ [2018] J. Amare _et al._ , Astropart. Phys. 97, 96 (2018), arXiv:1706.05818 [physics.ins-det] .
  * Adari  _et al._ [2022] P. Adari _et al._ , SciPost Phys. Proc. 9, 001 (2022), arXiv:2202.05097 [astro-ph.IM] .
  * Combecher [1980] D. Combecher, Radiation Research 84, 189 (1980).
  * Parks  _et al._ [1972] J. E. Parks, G. S. Hurst, T. E. Stewart, and H. L. Weidner, The Journal of Chemical Physics 57, 5467 (1972).
  * Dalgarno  _et al._ [1999] A. Dalgarno, M. Yan, and W. Liu, The Astrophysical Journal Supplement Series 125, 237 (1999).
  * Katsioulas  _et al._ [2022b] I. Katsioulas, P. Knights, and K. Nikolopoulos, Astropart. Phys. 141, 102707 (2022b), arXiv:2105.01414 [hep-ex] .
  * Vinagre and Conde [2000] F. L. R. Vinagre and C. A. N. Conde, Journal of Applied Physics 88, 5426 (2000).
  * Cowan  _et al._ [2011] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011), [Erratum: Eur.Phys.J.C 73, 2501 (2013)], arXiv:1007.1727 [physics.data-an] .
  * Cheng  _et al._ [2021] C. Cheng _et al._ (PandaX-II), Phys. Rev. Lett. 126, 211803 (2021), arXiv:2101.07479 [hep-ex] .
  * Abramowitz and Stegun [1972] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ (United States Department of Commerce, National Bureau of Standards (NBS), 1972).
  * [147] QEDark package, https://github.com/tientienyu/QEdark.