Atomic responses to general dark matter-electron interactions

In the leading paradigm of modern cosmology, about 80% of our Universe's matter content is in the form of hypothetical, as yet undetected particles. These do not emit or absorb radiation at any observable wavelengths, and therefore constitute the so-called Dark Matter (DM) component of the Universe. Detecting the particles forming the Milky Way DM component is one of the main challenges for astroparticle physics and basic science in general. One promising way to achieve this goal is to search for rare DM-electron interactions in low-background deep underground detectors. Key to the interpretation of this search is the response of detectors' materials to elementary DM-electron interactions defined in terms of electron wave functions' overlap integrals. In this work, we compute the response of atomic argon and xenon targets used in operating DM search experiments to general, so far unexplored DM-electron interactions. We find that the rate at which atoms can be ionized via DM-electron scattering can in general be expressed in terms of four independent atomic responses, three of which we identify here for the first time. We find our new atomic responses to be numerically important in a variety of cases, which we identify and investigate thoroughly using effective theory methods. We then use our atomic responses to set 90% confidence level (C.L.) exclusion limits on the strength of a wide range of DM-electron interactions from the null result of DM search experiments using argon and xenon targets.


I. INTRODUCTION
One of the major unsolved mysteries in modern physics is the elusive nature of Dark Matter (DM) [1].Revealed via the gravitational pull it exerts on stars, galaxies, light and the Universe at large, DM is widely believed to be made of hypothetical particles which have so far escaped detection [2].Detecting the particles forming the Milky Way DM component and understanding their properties in terms of particle physics models is a top priority in astroparticle physics [3].The experimental technique known as direct detection will play a major role in elucidating the nature of DM in the coming years [4].It searches for signals of DM interactions in low background deep underground detectors [5,6].Such interactions could produce observable nuclear recoils via DM-nucleus scattering or induce electronic transitions via DM-electron scattering in the detector material [7,8].DM particles that are gravitationally bound to our galaxy do not have enough kinetic energy to produce an observable nuclear recoil if their mass is approximately below 1 GeV/c 2 .On the other end, they have enough kinetic energy to ionize an atom in a target material or induce the transition from the valence to the conduction band in a semiconductor crystal if their mass is above about 1 MeV/c 2 [9].Motivated by recent experimental efforts in the field of astroparticle physics [10], here we focus on DM particles of mass in the 1 -1000 MeV/c 2 range and on their interactions with the electrons in materials used in direct detection experiments.In the literature, this framework is referred to as the "sub-GeV" or "light DM" scenario [11].
In order to interpret the results of direct detection experiments searching for signals of DM-electron interactions, it is crucial to understand quantitatively how detector materials respond to general DM-electron interactions.Materials' responses to DM-electron interactions can be quantified in terms of the overlap between initial (before scattering) and final (after scattering) electron wave functions.The stronger the material response, the larger the expected rate of DM-electron interactions in the target.
In the pioneering works [8,9,18], detector materials' responses to DM-electron interactions have been computed under the assumption that the amplitude for DM scattering by free electrons only depends on the initial and final state particle momenta through the momentum transferred in the scattering.While this restriction significantly simplifies the calculation of the predicted DM-electron interaction rate in a given detector, it prevents computations of physical observables whenever the scattering amplitude exhibits a more general momentum dependence.A more general dependence on the particle momenta occurs in a variety of models for DM [33], including those where DM interacts with electrons via an anapole or magnetic dipole coupling [34].
In this article, we calculate the response of materials used in operating direct detection experiments to a general class of non-relativistic DM-electron interactions.We build this class of non-relativistic interactions by using effective theory methods developed in the context of DMnucleus scattering, e.g.[33,[35][36][37].In terms of dependence on the incoming and outgoing particle momenta, this class of interactions generates the most general amplitude for the non-relativistic scattering of DM particles by free electrons.The predicted amplitude is compatible with momentum and energy conservation, and is invariant under Galilean transformations (i.e.constant shifts of particle velocities) and three-dimensional rotations.At the same time, it depends on the incoming and outgoing particle momenta through more than only the momentum transferred.It can explicitly depend on a second, independent combination of particle momenta, as well as on the DM particle and electron spin operators.We present our results focusing on DM scattering from electrons bound in argon and xenon atoms, since these are the targets used in leading experiments such as XENON1T [15] and DarkSide-50 [12].However, most of the expressions derived in this work also apply to other target materials.Within this general theoretical framework, we find that the rate of DM-induced electronic transitions in a given target material can be expressed in terms of four, independent material response functions.They depend on scalar and vectorial combinations of electron momenta and wave functions.We numerically evaluate these response functions for argon and xenon targets, and use them to set 90% C.L. exclusion limits on the strength of DM-electron interactions from the current null result of the XENON10 [38], XENON1T [15] and DarkSide-50 [12] experiments.We also find that the contribution of each material response function to the rate of DMinduced electronic transitions is weighted by a corresponding combination of particle momenta which is independent of the initial and final electron wave functions.We refer to these weights as "DM response functions", in analogy to previous studies in the context of DM-nucleus scattering [36].
Previous works could only identify one of the four material response functions found here because they assumed that the amplitude for DM scattering by free electrons solely depends on the momentum transfer, e.g.[18].This previously found response function is not only generated by DM-electron scattering.It also arises when materials are probed with known particles, such as photons or electrons which interact via the electromagnetic force.In contrast, the three novel response functions found in this work can only arise from a general treatment of DM-electron interactions.We will show this using an effective theory approach to DM-electron interactions.Specific examples of scenarios where the new response functions appear include models where DM interacts via an anapole or magnetic dipole coupling to photons.
These general considerations lead us to the intriguing conclusion that there are "hidden" material properties that have not yet been probed and that can only be investigated if the external probe is DM.Consequently, the detection of DM particles at direct detection experiments might not only solve one of the most pressing problems in astroparticle physics, but also open up a new window on the exploration of materials' responses to external probes.From this perspective, the results of this work constitute the first step within a far-reaching program which has the potential to open an entirely new field of research.
This work is organized as follows.In Sec.II we derive a general expression for the rate of DM-induced electronic transitions in target materials which we then specialize to the case of general, non-relativistic DM-electron interactions, focusing on argon and xenon as target materials.In Sec.III, we express the rate of DM-induced electronic transitions in terms of DM and atomic response functions, while in Sec.IV we apply our results to compute physical observables for argon and xenon targets and set limits on the strength of DM-electron interactions from the current null results of the XENON10, XENON1T and DarkSide-50 experiments.We present our conclusions in Sec.V and provide the details of our calculations in the appendices.In addition, we provide the numerical code used to evaluate the atomic response functions [39] 1 .

II. DARK MATTER-INDUCED ELECTRONIC TRANSITIONS
In this section we derive a general expression for the rate of electronic transitions induced by the scattering of Milky Way DM particles in target materials (Sec.II A).This result applies to arbitrary DM-electron interactions (including relativistic interactions) and target materials.We then specialize our general expression to the case of non-relativistic DM-electron interactions (Sec.II B) and to argon and xenon as target materials (Sec.II C).We use these results to model the response of atoms to general DM-electron interactions in Sec.III and set limits on the strength of such interactions from the null results of operating direct detection experiments in Sec.IV.
1. Illustration of elastic (left) and inelastic (middle, right) DM-electron scatterings.In the inelastic case, the electron gets excited by DM, transitioning from an atomic bound state with negative energy E1 to a final state with higher energy E2, which is either still negative for excitation (b) or positive for ionization (c).

A. Electron transition rate
We start by calculating the rate of transitions from an initial state |i ≡ |e 1 , p = |e 1 ⊗|p to a final state |f ≡ |e 2 , p = |e 2 ⊗ |p , where initial and final DM particle states are labelled by the corresponding tridimensional momenta, p and p , respectively, while |e 1 (|e 2 ) denotes the initial (final) electron state.Here |e 1 and |e 2 are eigenstates of the electron's energy with eigenvalue E 1 and E 2 , but not of momentum.This is illustrated in the middle and right panel of Fig. 1.Following [18], single-particle states are normalized as and, similarly, e.g., e 1 |e 1 = V .The divergent volume factor, V , will not appear in the expressions for physical observables.
The first non-trivial term in the Dyson expansion of the S-matrix element corresponding to the |i → |f transition is where H I (x) (H S (x)) is the interaction Hamiltonian density in the interaction (Schrödinger) picture at the spacetime point x, and H 0 is the Hamiltonian of the DM-electron system with H I (x) set to zero which obeys We also used the identity where |k are eigenstates of the free-electron Hamiltonian.The asymptotic states |i and |f are eigenstates of H 0 because H I is assumed to be different from zero within a finite time window only (the so-called adiabatic hypothesis).In the non-relativistic limit, the corresponding eigenvalues, denoted here by E i and E f , respectively, are given by where q = p − p is the momentum transferred in the scattering, v (v) is the initial DM particle velocity (speed), m χ and m e are the DM and electron mass, respectively, and E 1 (E 2 ) is the energy of the initial (final) electron bound state.
where θ qv is the angle between the vectors q and v, while q = |q|.Notice that S f i has dimension [energy] −6 , because of our choice of single-particle state normalization, Eq. ( 1).
Let us now consider the virtual process |k, p → |k , p , where a DM particle of initial (final) momentum p (p ) scatters elastically off a free electron of initial (final) momentum k (k ) as shown in the left panel of Fig. 1.The first non-trivial term in the Dyson expansion of the S-matrix element corresponding to this process is where we defined Ẽf ≡ E k +E p and Ẽi ≡ E k +E p .The above S-matrix element can equivalently be expressed as follows where M is the amplitude for DM scattering by a free electron.Compared to Eq. (4.73) of [40], Eq. ( 8) differs by a factor of 1/ 2E p 2E k 2E k 2E p because of our different normalization of single-particle states (see Eq. ( 1)).By comparing Eq. ( 7) with Eq. ( 8), we obtain the following identity where in the non-relativistic limit By substituting the free scattering result of Eq. ( 9) back into the general S-matrix element in Eq. ( 2), we obtain where we used the non-relativistic expression for 2E p 2E k 2E k 2E p .We now introduce the electron wave functions 1) and Eq. ( 3), a simple calculation shows that the wave functions ψ 1 and ψ 2 are normalized to one.In terms of ψ 1 and ψ 2 , Eq. ( 10) can be written as follows, The probability for the transition |e 1 , p → |e 2 , p is then given by where the divergent factor T = dt arises when squaring the one-dimensional Dirac delta in Eq. ( 11).While T enters in Eq. ( 12) explicitly, physical observables will not depend on T .In order to simplify the notation, in Eq. ( 12) we replaced M(k, p, k + q, p − q) with M(k, p, q).The rate per unit DM number density for the transition |e 1 , p → |e 2 , p with p within (p , p + dp ) is therefore given by P(p)V /T .Indeed, there is just one DM particle with initial momentum p within (p , p + dp ) in a phase-space element of three-dimensional volume V [41].Finally, integrating P(p)V /T over the initial DM particle velocity distribution f χ and the transferred momentum q, and multiplying by the local DM number density, n χ2 , for the total rate of DM-induced |e 1 → |e 2 transitions we find3 Here, we defined the squared electron transition amplitude, |M 1→2 | 2 , in terms of ψ 1 , ψ 2 and the free scattering amplitude, where a bar denotes an average (sum) over initial (final) spin states.
In the following two subsections, we motivate and explore the implications of our assumptions for M, ψ 1 and ψ 2 .Here, we model M within a non-relativistic effective theory of DM-electron interactions (Sec.II B) and use initial and final electron wave functions derived for isolated atoms in argon and xenon targets [42] for ψ 1 and ψ 2 (Sec.II C).

B. Non-relativistic Dark Matter-electron interactions
In this subsection, we investigate the non-relativistic limit of DM-electron interactions.Our goal is to show that the standard assumption in the context of DM direct detection, M = M(q), is highly restrictive.Rather than basing our argument on a specific model, we build a general non-relativistic effective theory where the active degrees of freedom are electrons and DM particles.The basic symmetries governing the DM-electron scattering, together with the assumption that both the DM particle and the electron are non-relativistic, will determine the allowed interactions between the active degrees of freedom in our effective theory, and allow us to show that M = M(q) is in general a too restrictive assumption.The symmetries that govern the nonrelativistic DM-electron scattering are the invariance under time and space translations (leading to energy and momentum conservation), the invariance under Galilean transformations (namely, constant shifts of particle velocities) and the invariance under three-dimensional rotations.Galilean invariance replaces the invariance under Lorentz boosts of relativistic theories.Let us now explore the constraints set by these symmetries on the amplitude for non-relativistic DM-electron scattering.
While the non-relativistic scattering of Milky Way DM particles by free electrons is characterized by the three-dimensional momenta k (k ) and p (p ) of the incoming (outgoing) electron and DM particle, respectively, momentum conservation and Galilean invariance imply that only two out of the four three-dimensional vectors are actually independent.A convenient choice of independent momenta is q = p − p , the momentum transferred in the scattering, and where µ e denotes the reduced mass of the DM particle and the electron, and v ≡ p/m χ is the incoming DM velocity.In the case of elastic DM-electron scattering, v ⊥ el • q = 0 because of energy conservation.In contrast, when DM scatters from electrons inelastically (and initial and final electrons populate different energy levels), v ⊥ el • q = 0.In this case, one could define such that v ⊥ inel • q = 0 by construction.While M can equivalently be expressed in terms of v ⊥ inel or v ⊥ el [33], we find that it is more convenient to use v ⊥ el as a basic variable when matching explicit models to the nonrelativistic effective theory built in this subsection (compare for example Eqs. ( 18) and ( 21) below).
To summarize, within our effective theory the free amplitude for non-relativistic DM-electron scattering reads where we emphasized the dependence of the amplitude M on the vectors q and v ⊥ el .In addition to q and v ⊥ el , M can also depend on matrix elements of the electron and DM particle spin operators, denoted here by S χ and S e , respectively.Within the non-relativistic effective theory of DM-electron interactions developed here, the amplitude M does not explicitly depend on the properties that characterize the particle mediating the interactions between DM and electrons, such as the mediator mass, m med .There are two scenarios in which our effective theory approach can be applied.In the first one, m med is much larger than the momentum transfer in DMelectron scattering (contact interaction).In the second one, the mediator is effectively mass-less (long-range interaction).These are the two limiting cases that we investigate here.For m med ∼ |q|, M is model dependent, and some of the considerations made here do not apply.
In Eq. (18) we are assuming that both DM particle and electron are non-relativistic.Indeed, Milky Way DM particles have a typical speed of the order of 10 −3 in natural units, and the typical speed of electrons bound in atoms is αZ eff , where α = 1/137 and Z eff is 1 for electrons in outer shells (the most interesting ones from a detection point of view) and larger for inner shells.Notice also that for m χ ≥ 1 MeV/c 2 the electron is the fastest and lightest particle in the non-relativistic DM-electron scattering, which implies that the typical momentum transferred in the scattering is of the order of Z eff αm e [18].Eq. ( 18) can then be expanded in powers of |q|/m e and |v ⊥ el |.Each term in this expansion of M must be both Galilean invariant and scalar under three-dimensional rotations.At linear order in v ⊥ el and at second order in q, there are fourteen combinations of q, v ⊥ el , S χ and S e fulfilling the above requirements when DM has spin 1/2 [35,36,43].Specifically, these three-dimensional scalar combinations are operators in the DM/electron spin space.They correspond to the first fourteen entries in Tab.I.The remaining entries were found to arise from the non-relativistic reduction of models for spin 1 DM coupling to nucleons via a vector mediator [37].We include them in Tab.I after replacing nucleon with electron operators.Based on the above considerations, the amplitude M can in general be . Scalar combinations of the three-dimensional vectors q, v ⊥ el , Sχ and Se contributing to the DM-electron scattering amplitude M in Eq. (19).They are operators in the DM/electron spin space and invariant under Galilean transformations.1χ and 1e are 2 × 2 identity matrices in the DM particle and electron spin space, respectively.The combinations O17, O18, O19 and O20 were found to arise from the non-relativistic reduction of models for spin 1 DM with a vector mediator.We follow the numbering used for nonrelativistic DM-nucleon interactions [37], and neglect O2, which is quadratic in v ⊥ el , and O16, which can be expressed as a linear combination of O12 and O15.For further details, see for example [33].
expressed as the linear combination4 where the reference momentum, q ref , is given by q ref ≡ αm e , the interaction operators O i are defined in Tab.I (with i running on subsets of Tab.I, depending on the DM particle spin), and the coefficients c s i (c i ) refer to contact (long-range) interactions.The case m 2 med |q| 2 corresponds to the limit of contact interaction, c i = 0.In contrast, for m 2 med |q| 2 the interaction between DM and electrons is long-range, i.e. c s i = 0. Angle brackets in Eq. (19) denote matrix elements between the twodimensional spinors ξ λ and ξ s (ξ λ and ξ s ) associated with the initial (final) state electron and DM particle, respectively.For example, O 1 = ξ s ξ s ξ λ ξ λ .An analogous expansion for M was found previously in studies of non-relativistic DM-nucleon interactions (see e.g.[33] for a review).
As an illustrative example, let us compute M for DMelectron interactions arising from the DM anapole mo-ment and match the result of this calculation to the nonrelativistic expansion in Eq. (19).The DM anapole moment is potentially important being the leading coupling between DM and photons if DM is made of Majorana fermions.It generates DM-electron interactions via the electromagnetic current.The Lagrangian of the model reads where χ is a Majorana fermion describing the DM particle, g is a dimensionless coupling constant, Λ a mass scale and F µν the electromagnetic field strength tensor.The corresponding non-relativistic scattering amplitude is given by where g e is the electron g-factor, while ξ λ and ξ s are two-component spinors for the electron and DM particle, respectively.The first line in Eq. ( 21) depends on v ⊥ el and is proportional to O 8 , whereas the second line is linear in q and proportional to O 9 .This example shows that the assumption commonly made in the interpretation of direct detection experiments searching for DM-electron interactions, M(q, v ⊥ el ) = M(q), is likely restrictive, and contributions to M depending on v ⊥ el are in general expected.
Let us now calculate the squared transition amplitude assuming that the expansion in Eq. ( 19) applies.This will lead us to introduce two atomic form factors depending on ψ 1 , ψ 2 and M. We start by noticing that the free electron scattering amplitude in Eq. ( 19) can also be written as where we used v ⊥ el = v − k/m e − q/(2µ χe ).This is not an approximation because Eq. ( 19) is at most linear in v ⊥ el .Substituting Eq. ( 22) into the squared transition amplitude in Eq. ( 14), we find While the first term in Eq. ( 23) has already been considered in previous studies, the second and third terms are new, and arise from our general modelling of the free electron scattering amplitude.Indeed, in the analysis of direct detection experiments searching for DM-electron interaction signals, it is standard to assume that the DMelectron free scattering amplitude M depends solely on the momentum transfer q and not on the velocity v ⊥ el and, therefore, on the electron's initial momentum k.In this case, the scattering amplitude M(q) can be taken out of the integral in Eq. ( 13), and only the first term in Eq. ( 23) contributes to the squared transition ampli- The advantage of this standard assumption is that M = M(q) implies a convenient factorization of atomic physics and DM physics.Following the notation of [18], for M = M(q) one can define a scalar atomic form factor, and obtain the transition amplitude, M 1→2 (q), by multiplying the free scattering amplitude by the form factor f 1→2 (q), Since the momentum space wave functions ψ 1 and ψ 2 have dimension [energy] −3/2 , f 1→2 (q) is a dimensionless quantity.
Besides the well-known scalar atomic form factor of Eq. ( 24), a new, vectorial atomic form factor appears in Eq. ( 23), The details of the evaluation of the scalar, f 1→2 , and vectorial, f 1→2 , atomic form factors are presented in App.B 3.

C. Ionization of isolated atoms in argon and xenon targets
Our treatment of DM-electron scattering in target materials has so far been general in terms of initial and final state electron wave functions, ψ 1 and ψ 2 , respectively.From now onward, however, we specialize our results to the case of DM-induced ionization of isolated atoms.In this case the initial state (formerly simply denoted by "1") is a bound state characterized by the principal, angular and magnetic quantum numbers (n, , m).The final state ("2") is a free electron state at large distance from the atom, but still affected by the remaining ion's presence at low distance.It is defined by the quantum numbers (k , , m ), where k is the asymptotic momentum of the electron and , m are its angular and magnetic quantum numbers.We obtain the total ionization rate R n ion of a full atomic orbital (n, ) by summing the transition rate R 1→2 over all occupied initial electron states and integrating over the allowed final electron states.
For the initial bound electrons, the rate has to be summed over all values of the magnetic quantum number m and multiplied by 2 to account for the spin degeneracy.In order to account for all the allowed electron final states, we have to act on R 1→2 with the integral operator [18] V 2 Here, E e = k 2 /(2m e ) is the ionized electron's final energy, and the number of final states with asymptotic momentum between k and k + dk is V d 3 k /(2π) 3 .Summarizing, the total ionization rate for the (n, ) orbital is given by and the corresponding final state electron ionization energy spectrum by The divergent factor V in Eq. ( 29) cancels with the 1/V factor arising from the normalization of ψ k m (see App. B 4). Substituting Eq. ( 13), we can rewrite the spectrum as 10 0.
FIG. 2. The ionization form factor f n ion (k , q) 2 defined in Eq. ( 33) for the outer electron orbital of argon (left) and xenon (right).
× dq q where, following [18], we integrated over cos θ qv while assuming that f χ (v) = f χ (v), and then replaced f χ (v) with f χ (v) in the final expression 5 .The squared ioniza-tion amplitude, M n ion 2 , appearing in Eq. ( 30) is defined as follows and can explicitly be expressed in terms of the amplitude Here, we defined ξ = ∆E 1→2 /(qv) + q/(2m χ v).In order to connect to the results and notation of previous works, in Eq. ( 32) we introduced the dimensionless ionization form factor, This ionization form factor is depicted in Fig. 2 for the outer atomic orbitals of argon and xenon.In Eq. ( 30), the minimum DM speed, v min , required to deposit an energy of ∆E 1→2 given a momentum transfer q is Finally, the wave functions of the initial and final state electrons can be expanded in spherical harmonics, Their radial parts, R n (r) and R k (r), are given in App.B 4 for the case of isolated argon and xenon atoms [42].While these wave functions have also been used in previous works on DM direct detection [12], we note that they are not fully applicable to dense liquid argon and xenon systems [44].Our neglect of the difference in electronic structure between the isolated atoms and the liquid state makes our results conservative since the electron binding energies in liquid nobles are smaller than those of the isolated atoms.Furthermore, this treatment of the wave functions can be improved by including relativistic corrections [23,45] and many-body effects [46].However, both improvements, in particular the relativistic corrections, are expected to have limited impact on theoretical predictions for sub-GeV DM.

III. GENERAL DARK MATTER AND ATOMIC RESPONSE FUNCTIONS
In this subsection, we reformulate the expressions found in Sec.II C in order to obtain compact equations which are suitable for numerical applications.By doing so, we also investigate and characterize the atomic response to DM-electron interactions of the general form considered here.We start by observing that the squared ionization amplitude, |M n ion | 2 , can in general be rewritten as follows Here, we define a set of DM response functions R n i v ⊥ el , q me , which are functions of the couplings in Eq. ( 19), and atomic response functions W n i (k , q), which encapsulate the information on the electron's initial and final state wave functions 6 .By means of Eq. ( 37), we are again able to disentangle the particle physics input from that of the atomic physics.However, our general treatment of DM-electron interactions predicts not one but four atomic response functions.At this point, we summarize the final expressions for the DM and atomic response functions.Their detailed derivations can be found in App.A and B.
The four DM response functions, R n i v ⊥ el , q me , are given by They indirectly depend on the (n, ) quantum numbers through the variable ξ = ∆E 1→2 /(qv) + q/(2m χ v), see also App. C. In the above expressions, v ⊥ el is evaluated at k = 0 due to the non-relativistic expansion of M. The four atomic response functions W n i (k , q), which we define next in Eq. (39), are one of the main results of this work: The first atomic response function W n 1 (k , q) can be identified with the ionization form factor f n ion (k , q) 2 commonly used in the sub-GeV DM detection literature [18].The atomic response functions W n j (k , q), j = 2, 3, 4 were not considered in previous studies and can only arise from a general treatment of DM-electron interactions.While W n 4 /W n 1 roughly scales as particle momenta to the power of 4, and W n 2,3 /W n 1 as particle momenta to the power of 2, numerically we find that the products R n j W n j are often of comparable magnitude for several DM-interaction types (see below for further details).
In order to clarify the physical meaning of the new atomic responses, let us evaluate them in the planewave limit, where V are eigenstates of the electron momentum operator with eigenvalues k and k , respectively.In this limit, initial and final state electron are not bound to an atom and propagate as free particles.Substituting ψ 1 (x) and ψ 2 (x) in Eqs.(39a), (39b), (39c) and (39d), we find In the plane-wave limit, we can define a laboratory frame in which the initial electron is at rest and, therefore, k = 0.In this frame, W n j (k , q) = 0, j = 2, 3, 4 and Eq.(37) reduces to the expression for the modulus squared of the amplitude for DM scattering by a point-like proton found in [47] (with the proton mass replaced by the electron mass).
Consequently, the new atomic responses describe distortions in the ionization spectrum induced by the fact that the initial state electron obeys a momentum distribution with a finite dispersion, being the electron bound to an atom.
Figs. 3 and 4 show the four atomic response functions, W n j , j = 1, . . ., 4, for the 3p and 5p atomic orbitals of argon and xenon, respectively.

IV. FIRST APPLICATION TO DIRECT DETECTION
The non-relativistic effective theory of DM-electron interactions presented in the previous sections culminated in a general expression for the electron ionization energy spectrum of isolated atoms due to generic DM-electron interactions given by the Eqs.( 30) and (37).This main result allows us to make predictions for direct searches of sub-GeV DM particles by using an almost modelindependent framework and a generic expression for the scattering amplitude in Eq. (19) in terms of effective operators O i .
In the following subsections, we apply this effective theory to evaluate ionization spectra.
There are a number of large-scale dual-phase time-projectionchamber (TPC) detectors with xenon and argon targets, which is why we focus on these two elements.By reinterpreting the observational data by XENON10 [13,38], XENON1T [15], and DarkSide-50 [12], we can set exclusion limits on either the effective coupling constants c s i and c i of Eq. ( 19) in front of individual operators, or on the parameters of specific DM models generating combinations of interaction operators, as in the case of DM with anapole and magnetic dipole couplings.In the latter case, one has to compute the elastic scattering amplitude between a DM particle and an electron, take the nonrelativistic limit, and identify the different contributions with one or several of the effective operators in Tab.I. We have already shown the example of anapole interactions, which corresponded to a combination of O 8 and O 9 , in Sec.II B.
It is also worth pointing out how the most-used interaction model in the sub-GeV DM literature, i.e. the "dark photon" model, emerges in this framework.The dark photon model extends the Standard Model of particle physics by an additional U(1) gauge group under which the DM particle is charged.Interactions between the DM particle and electrically charged particles in the Standard Model arise from a "kinetic mixing" term in the interaction Lagrangian between the photon and the gauge boson A associated with the new U(1) gauge group [9,48], i.e. εF µν F µν .Here, ε is a coupling constant and F µν the A field strength tensor.This model can arguably be considered as the "standard model" of sub-GeV DM detection.The new massive gauge boson A is called the dark photon and acts as the interaction's mediator.Its mass is usually assumed to be either much larger than the typical momentum transfer q ref in the scattering process (contact interaction), or much lower (longrange interaction).By considering the non-relativistic DM-electron scattering amplitude, this interaction can be identified with O 1 , and the connection between the effective couplings in Eq. ( 19) and fundamental parameters is given by for contact and long-range interactions, respectively.
Here, g D is the gauge coupling of the additional "dark" U(1) gauge group, while e denotes the charge of the electron.
As anticipated, the scattering amplitude arising from a given DM model can correspond to either one or several of the effective operators in Tab.I.This is why we start by focusing on single operators, compute and study the corresponding ionization spectra, and set exclusion limits on the constants c s i and c i for individual O i .To illustrate the case of amplitudes being linear combinations of two or more of the effective operators, we study three relevant exemplary cases of DM-electron interactions, the anapole, magnetic dipole, and electric dipole interaction.These arise as leading and next-to-leading terms in an effective theory expansion of the DM coupling to the ordinary photon [34].We briefly review these models in App.D.

A. Individual effective operators
Let us consider the case where the free electron scattering amplitude, M, consists of only one of the O i operators, or, in other words, only one of the couplings c s, i in Eq. ( 19) is different from zero.On the one hand, this is instructive and allows us to study the impact of single operators in isolation.On the other hand, this situation might also arise from specific models for DM interactions, with the dark-photon model being the most obvious example.
From Tab.I, we can identify the operators that generate the three new atomic responses, i.e.W n j (k , q), j = 2, 3, 4, by their dependence on v ⊥ el .They are O 3 , O 5 , O 7 , O 8 , O 12 , O 13 , O 14 , and O 15 .Fig. 5 shows the total ionization spectrum, dR ion /d ln E e , for two example operators, O 6 (top panels) and O 13 (bottom panels).In the case of the O 6 operator, we compute the ionization spectrum for argon.For the O 13 operator, we focus on xenon.Furthermore, the left panels in Fig. 5 refer to contact interactions, whereas the right panels refer to long-range interactions.Finally, in all panels solid, longdashed and short-dashed lines refer to different DM particle masses.We compute the total ionization spectrum by summing the ionization spectrum for each atomic orbital (n, ), Eq. ( 30), over the five outermost occupied orbitals for xenon (4s,4p,4d,5s,5p), and all occupied orbitals for argon (1s, 2s, 2p, 3s and 3p).From the overall scale of the spectra, we can draw two conclusions.Firstly, we find that ionization spectra are smaller in the case of long-range interactions than for contact interactions due to the additional q 2 ref /q 2 factor in the scattering amplitude.This suppression is more severe for large DM masses than for m χ around 1 MeV/c 2 .Secondly, we find that the magnitude of the ionization spectrum associated with O 6 (and scattering on argon) is comparable with that of O 13 (and scattering on xenon).This is expected, because O 6 is quadratic in q while O 13 is linear in q and v ⊥ el , and at the same time scattering on argon and xenon occurs with comparable probabilities.
To make the last observation more precise and complete, we now extend the above comparison to the 14 operators in Tab.I describing spin 1/2 DM.For each of these operators, we compute the associated ionization spectrum considering scattering on both argon and xenon targets.Fig. 6 shows the integrated ionization spectra in argon and xenon that we find by integrating dR ion /d ln E e from zero up to the maximum kinematically allowed energy.We present our results for three different DM particle masses, namely 10, 100, and 1000 MeV/c 2 .By comparing the magnitude of the ionization spectra for a given operator, we find that the integrated ionization spectra in argon and xenon are always comparable, while the xenon ionization rate slightly exceeds the rate in argon in most but not all cases.
By comparing the ionization spectra of different operators for a given target material, we identify O 1 and O 4 as the "leading-order" (LO) operators.This is expected, since the corresponding terms in Eq. ( 38) are not suppressed by powers of q or v ⊥ el , and, in addition, they generate the first atomic response function W n 1 , which is numerically found to be the one of largest magnitude.The operators O 7 to O 12 can be regarded as "next-toleading-order" (NLO) operators, in that they are linear in either q or v ⊥ el .Similarly, the "next-to-next-to-leadingorder" (NNLO) operators (O 5 , O 6 , O 13 , and O 14 ) are quadratic in q or in a combination of q and v ⊥ el .Compared to the LO and NLO cases, the total ionization spectrum of NNLO operators is further suppressed.Lastly, the only N 3 LO operator we include in our analysis is O 15 , which is quadratic in q and linear in v ⊥ el .The ionization spectrum due to this operator is the lowest and suppressed by up to ten orders of magnitude compared to O 1 for identical couplings.
This hierarchy of operators is essentially identical to that found when classifying DM-nucleus interactions within the effective theory of DM-nucleon interactions [35,36].In that context, the equivalent of the O 1 and O 4 operators defined here correspond to the "familiar" spin-independent (SI) and spin-dependent (SD) interactions -a standard benchmark in the analysis of DMnucleus scattering data.Similarly to the SI/SD paradigm dominating the literature on direct DM searches via nuclear recoils for a long time, the dark photon model underlies most studies of sub-GeV DM searches nowadays.
Another interesting question in the context of isolated operators regards the relative contribution of the four different terms in Eq. (37).Each of these terms is the product of a DM response function given in Eq. ( 38) and an atomic response function given in Eq. (39).In order  long-range interac�on FIG. 7. The four R n j W n j , j = 1, . . ., 4 contributions to the ionization energy spectrum (here generically denoted as "responses") for the two example couplings c s, 7 (left) and c s, 15 (right).We present our results for argon (above) and xenon (below), as well as for contact interactions with mχ = 10 MeV (left) and long-range interactions with mχ = 100 MeV (right).A dagger † indicates negative contributions.
to answer the above question, in Fig. 7 we show the ionization energy spectrum for two example operators, each of which generates multiple atomic responses.Specifically, we show the total argon and xenon ionization spectra for O 7 , contact interaction and m χ = 10 MeV/c 2 (left panels), and for O 15 , long-range interaction and m χ = 100 MeV/c 2 (right panels).As one can see from Fig. 7, we find that in the case of contact interactions of type O 7 the four terms in Eq. (38) give comparable contributions to the total ionization spectrum.This is the result of a compensation between DM and atomic response functions.Indeed, while the first atomic response function W n 1 generally dominates over the second and third, which are in turn larger than W n 4 , in the case of contact interactions of type O 7 we find an inverse hierarchy for the DM response functions R n j , j = 1, 2, 3, which scale as (v ⊥ el ) 2 , (v ⊥ el • q)m e /q 2 and q 0 , respectively.In this respect, the operator O 7 is not unique (see for example O 8 , just to name one).From this example, we conclude that the four terms in Eq. ( 39) can in general contribute in a comparable manner and must therefore be taken into account in the interpretation of experimental data.
The contribution of an individual response may also be negative, an example of which is shown in the right panels of Fig. 7 for the long-range ionization spectrum of O 15 7 .The effective operator O 15 generates the DM responses 1, 3, and 4. In this example, the contribution to the total ionization spectrum from the first atomic response function, W n 1 , is suppressed by the first DM response function, R n 1 , which renders the R n 1 W n 1 contribution to Eq. (39) completely negligible.In addition, the contribution of R n 4 W n 4 is as anticipated negative, which originates from the sign of the c 15 term in R n 4 in Eq. (38d).
We conclude this subsection with an analysis of present exclusion limits on the c s, i coupling constants from direct detection data.Despite tremendous experimental efforts by a multitude of direct detection experiments (see [4] for a review), no conclusive DM signal has been established, and DM continues to evade direct searches to this day.The largest of these experiments are dualphase time-projection-chambers (TPC) with xenon and argon targets, which is why we are focussing on these elements in this paper.Here, we re-interpret the null results by the xenon target experiments XENON10 [38] and XENON1T [15], as well as the argon experiment ) from XENON10 [13,38](blue), XENON1T [15](dark blue) and DarkSide-50 [12](purple) on a selection of individual effective couplings for both contact and long-range interactions.The dashed gray lines indicate constant reference cross sections σe, see Eq. (43).Note that the constraints on long-range interactions depend on a choice of the reference momentum transfer q ref .
DarkSide-50 [12] and set 90% C.L. exclusion limits on the individual effective coupling constants.The details underlying the computation of the limits for these three experiments are summarized in App.E. Fig. 8 shows our 90% C.L. exclusion limits for three selected interaction operators, one LO (O 1 ), one NLO (O 8 ), and one NNLO operator (O 13 ), both for contact (left) and long-range interactions (right) 8 .The dashed gray lines indicate constant values of a reference DM-electron scattering cross section that we define as follows, With this definition the constraints on c 1 can be directly compared with previous works which described the DMelectron interaction using the dark photon model, see e.g.[12][13][14].

B. Anapole, magnetic dipole, and electric dipole interactions
As an illustration of how different combinations of effective operators, O i , can arise from specific models 8 Due to the chosen parametrization of the scattering amplitude in Eq. ( 19), it is important to remember that the limits on the dimensionless long-range interaction coupling constants c i depend on the choice for the reference momentum transfer q ref = αme.
for DM-electron interactions, we consider three examples: the anapole, magnetic dipole, and electric dipole interaction.Electric and magnetic dipole interactions have previously been studied in the context of direct DM searches via electron scatterings [13,16,49].These studies used different approximations and did not exploit our new atomic response functions.In App.D, we provide explicit Lagrangians for these types of interactions.Therein we also present the associated non-relativistic DM-electron scattering amplitudes, which we then match on to the effective theory expansion in Eq. (19).In this way we can read off the connection between the Lagrangian parameters and our effective couplings.
Let us start with the anapole interaction.From the amplitude in Eq. ( 21), we find that this model generates a linear combination of the O 8 and O 9 operators in the nonrelativistic limit.Only contact interactions are generated and the corresponding effective coupling constants are given by where g is a dimensionless coupling constant, and Λ is the energy scale at which the anapole interaction term is generated.
For magnetic dipole interactions, we find that the scattering amplitude in Eq. (D5) can be expressed as a linear combination of four interaction operators, namely O 1 , O 4 , O 5 , and O 6 .Here, the effective coupling constants FIG. 9. Ionization spectra for anapole, magnetic dipole, and electric dipole interactions shown in black.We set the DM mass to 100 MeV and focus on xenon.The different blue lines show the contributions of the individual effective operator.Note that in the case of magnetic dipole interactions, the sum of the four individual contributions exceeds the total spectrum due to a negative interference between O4 and O6 in Eq. (38a).
are given by The situation is particularly simple for electric dipole interactions, with the scattering amplitude given in Eq. (D6).It is linear in O 11 with the single effective coupling, Fig. 9 shows the xenon ionization spectra for anapole, magnetic dipole and electric dipole interactions.In the three cases we assume a DM particle mass of 100 MeVc 2 .Different blue lines correspond to contributions from the individual interaction operators to the total ionization spectrum (shown as a black solid line).Most interesting is the spectrum of the magnetic dipole interaction.While one might expect the interaction operator O 1 to dominate the spectrum (since it is one of the two LO operators), its contribution is in fact negligible due to the small relative size of the four effective coupling constants.This can be seen e.g. by the ratio c s 1 /c s 4 ≈ 10 −3 for m χ = 100 MeV/c 2 .
Another interesting aspect of the magnetic dipole interaction spectrum is the fact that the contribution of the single operator O 4 exceeds the total result.This is explained by the interference term between O 4 and O 6 in the first DM response function in Eq. (38a), whose overall contribution is negative due to the sign of c 6 in Eq. (45).
We conclude this chapter by presenting our 90% C.L. exclusion limits on the ratio between dimensionless coupling constant g and energy scale Λ (or Λ 2 ) for the anapole, magnetic dipole, and electric dipole interaction.Our exclusion limits are shown in Fig. 10 and based on a re-interpretation of the null results reported by the XENON10, XENON1T, and DarkSide-50 experiments 9 .For these specific DM-electron interaction models, we find that exclusion limits from DarkSide-50 are generically less stringent, as compared to those arising from XENON10 and XENON1T data.While XENON10 always sets the most stringent limits for DM masses below about 20-50 GeV/c 2 (depending on the assumed interaction), the minimum coupling constant value that XENON10 can exclude is comparable to the minimum coupling probed by XENON1T only for electric dipole interactions.

V. SUMMARY
We computed the response of isolated xenon and argon atoms to general, non-relativistic DM-electron interactions.This is a key input to the interpretation of DMelectron scattering data.For example, dual-phase argon and xenon targets are used in direct detection experiments searching for signals of non-relativistic interactions between Milky Way DM particles and electrons in the target materials.We modelled the DM-electron scattering by formulating a non-relativistic effective theory of DM-electron interactions which significantly extends the currently favoured framework, namely, the dark photon model.Within our effective theory description of DMelectron interactions, the free amplitude for DM-electron scattering not only depends on the momentum transferred in the scattering (as in the dark photon model), it can also depend on a second, independent combination of particle momenta as well as on the DM particle and electron spin operators.Quantitatively, we defined the atomic response to an external probe in terms of the overlap between the initial and final state electron wave functions.
As a first application of the atomic responses found in this work, we computed the rate at which Milky Way DM particles can ionize isolated argon and xenon atoms in target materials used in operating direct detection experiments.We found that the final state electron ionization energy spectrum can in general be expressed as the sum of four independent terms.Each of these terms is given by the product between an atomic response function, which depends on the initial and final state electron wave functions, and a DM response function, which only depends on kinematic variables, coupling constants, and the energy gap between initial and final electron states.We investigated the structure and relative strength of the four terms contributing to the total ionization spectrum, describing under what circumstances they can be generated in non-relativistic DM-electron scattering processes in argon and xenon targets.The effective theory approach developed in this work has proven very useful in the classification and characterization of the individual atomic responses and contributions to the total ionization energy spectrum.At the same time, the generality of the formalism developed here allows for straightforward applications to most specific DM interaction models.
We then used the argon and xenon atomic responses found here to compute 90% C.L. exclusion limits on the strength of DM-electron interactions from the null result reported by the XENON10, XENON1T and DarkSide-50 direct detection experiments.DarkSide-50 is the only direct detection experiment considered here that uses an argon target.For contact DM-electron interactions, we found that XENON10 generically sets the most stringent constraints for DM particle masses below O(10) MeV and that XENON1T is the leading experiment at larger masses.The DM particle mass above which XENON1T dominates depends on the specific DM-electron interaction model.For long-range interactions, we found cases in which XENON10 sets the strongest exclusion limits on the entire mass window explored here.Both for contact and long-range interactions, the DarkSide-50 experiment currently sets 90% C.L. exclusion limits which are generically less stringent than those from XENON10 and XENON1T.
Considering the plethora of proposed DM-electron scattering targets, the results presented in this work should be thought of as a first step within a far-reaching program which aims at exploring the response to general DM-electron interactions of condensed matter systems used in (or proposed for) DM direct detection experiments.Indeed, many of the expressions derived in this work would also apply to target materials different from isolated argon and xenon atoms, such as semiconductor crystals and quantum materials.Finally, we would like to stress that only one of the four atomic response functions computed here can arise from standard electromagnetic interactions between external electrons or photons and electrons bound to target materials.The remaining three atomic responses can only be generated when DM is the external probe used to investigate the properties of condensed matter systems.From this perspective, our findings show that the detection of Milky Way DM particles has the potential to open up a new window on the exploration of materials' responses to external probes, revealing as of yet hidden properties of matter.This appendix serves as a bridge between Eq. ( 32) and Eq.(37).In particular, we identify and derive the atomic and DM response functions given in Eqs.(38) and (39).This involves the computation of the three terms in Eq. ( 32), given more explicitly by where we average (sum) over initial (final) spins of the DM particle and the electron.For the amplitude we subsitute the general expression of Eq. ( 19), where we consider the first 14 operators listed in Tab.I.
When averaging (summing) the terms in Eqs.(A1) to (A3) over the initial (final) spin of the DM particle and the electron, the non-vanishing terms are of the form [47], Here, s(s ) is the initial (final) spin's index of either the electron or DM particle, and A and B are two general vectors.All cross terms which are linear in S i vanish.It is important to keep in mind that the scalar product q • v ⊥ el , which occurs frequently, does not vanish since we are not considering an elastic scattering process.

Response 1
We evaluate the spin sums in Eq. (A1) for the scattering amplitude including the first 14 effective operators.For the evaluation of the electron and DM spin sums, we make use of the identities of Eq. (A4).We find For real couplings (c i = c * i for all i) the squared amplitude simplifies to The square of the amplitude can be identified as the first DM response function, which was given in Eq. (38a) since the first atomic response function is of the form as we saw in Eq. (39a).

Response 2
By using the spin-sum identities of Eqs.(A4), the second term given by Eq. (A2) results in To keep the expressions more manageable, we introduced the vector A ≡ f 1→2 (q)f * 1→2 (q).For real couplings and a real vector A, this simplifies to The first line corresponds to the second response in Eq. ( 37), since we recognize the second atomic response function, Indeed, the vector A ≡ mm A is real which justifies our assumption above.
Furthermore, it is not obvious that the second term in Eq. (A10), which is proportional to (v ⊥ el • A), becomes part of the second DM response function as well.The reason for this is the fact that the vector A is (anti-)parallel to q, which allows us to write Hence, we obtain the final expression for the second DM response function, This is the expression we presented in Eq. (38b) of this paper.Note that the c 3 , c 5 , and c 15 terms as well as the c 12 c 15 interference term cancel due to the use of Eq. (A12).

Response 3 and 4
Lastly, we evaluate the third term in a similar fashion.Substituting the 14 first operators into Eq.(A3), we find For real couplings, this becomes From the definition of the third and fourth atomic response functions in Eqs.(39c) and (39d), we can directly assign the first and second line of Eq. (A15) to the third and fourth DM response functions, These two response functions were presented in the paper as Eqs.(38c) and (38d).The third line of Eq. (A15) does not contribute to ionization rates, since the vector mm (f 1→2 (q) × f * 1→2 (q)) vanishes.
This completes the derivation of the atomic and DM response functions.
Appendix B: Evaluation of the atomic response functions

Useful identities
This chapter contains a number of mathematical identities necessary for the evaluation of the atomic response functions.
a. Spherical harmonic addition theorem The spherical harmonics, Y m (θ, φ) are given by where P m (cos θ) are the associated Legendre polynomials, and N m is a normalization constant ensuring that dΩ (Y m ) * Y m = δ δ mm .The normalization constant is given by The spherical harmonic addition theorem states that for two given unit vectors defined by the spherical angle coordinates (θ 1 , φ 1 ) and (θ 2 , φ 2 ), the sum over all values of m yields [56] m=− Here, P (x) is a Legendre polynomial, and ω is the angle between the two vectors given by cos ω = cos θ 1 cos θ 2 + sin θ 1 sin θ 2 cos(φ 1 − φ 2 ) .(B3b) b. Plane-wave expansion A plane-wave e ik•x can be expressed as a linear combination of spherical waves [57], where ω is the angle between k and x.Using the spherical harmonic addition theorem in Eq. (B3), we can re-write this equation as = 4π c. Identities with the Wigner 3j-symbol The angular integral of the product of three spherical harmonics can be expressed in terms of the Wigner 3j-symbol [56], The 3j-symbols satisfy the following orthogonality relation [57], where ∆(j 1 , j 2 , j) is equal to 1 if (j 1 , j 2 , j) satisfy the triangular condition, i.e. |j 1 − j 2 | ≤ j 3 ≤ j 1 + j 2 , and zero otherwise.

Gradient of a wave function
In order to compute the vectorial atomic form factor in Eq. (B18), we need to evaluate the gradient of a wave function ψ nlm (x) as defined in spherical coordinates (r, θ, φ) by Eq. ( 35) which is given by This gradient can be expressed in terms of vector spherical harmonics (VSH), where the VSH are defined following the conventions of Barrera et al. [58], The Cartesian components of the VSHs can themselves be expanded in spherical harmonics with coefficients closely related to Clebsch-Gordan coefficients [57], with the non-vanishing coefficients c (1) .
Here, we used the notation of [57], A Similarly, the second VSH can be expanded as with coefficients, d −1 m+1 = id (2) +1 m+1 = id (2) Eqs. (B11) and (B13) allow us to the express the gradient of the wave function ψ lmn in a compact expansion in spherical harmonics, where e i , (i = 1, 2, 3) are the three Cartesian unit vectors.

The atomic form factors
In Eqs. ( 24) and ( 26), we introduced the scalar and vectorial atomic form factors.In this section, we describe the details of evaluting these form factors.
First, we denote the initial ('1') and final ('2') state electron wave functions by their respective quantum numbers (n, , m) and (k , , m ).Consequently, the atomic form factors take the form Moving to position space via a Fourier transformation, we find In Eqs.(35) and (36) we expressed the initial and final state electron wave functions in terms of spherical coordinates x(r, θ, φ), where the angular dependence is given by the spherical harmonics, We postpone the discussion of the radial components of the wave functions to App.B 4. and start with the evaluation of the scalar atomic form factor f 1→2 .Afterwards, we discuss the vectorial atomic form factor f 1→2 .a. Scalar atomic form factor Focussing on the scalar atomic form factor for now, we replace the exponential e ix•q with the plane-wave expansion of Eq. (B4).
where we absorbed the radial integral into This integral is the spherical Bessel transform of R * k (r)R n (r), i.e. the product of the radial components of the final and initial wave functions.For its evaluation, we refer to App.B 5.
The integral over three spherical harmonics in the last expression can be re-written in terms of the Wigner 3jsymbols via Eq.(B5), At this point, we can already evaluate the first atomic response W n 1 given by Eq. (39a), which essentially given by the squared modulus of the scalar atomic form factor summed over m, m , The orthogonality of the Wigner 3j-symbols, Eq. (B6), allows us to sum over the L and M indexes, Lastly we apply the spherical harmonic addition theorem given in Eq. (B3), Note that the dependence on the momentum transfer's direction has disappeared.The explicit form of the first atomic response W n 1 (q, k ) for an atomic orbital (n, ) given by Eq. (39a) (or equivalently the 'ionization form factor |f n, ion | 2 defined in Eq. ( 33)) is obtained by adding the remainder of the final state phase space10 , see Eq. ( 27), In order to evaluate the three new atomic responses we need to follow similar steps for the vectorial atomic form factor as for the scalar one.
b. Vectorial atomic form factor We start from Eq. (B18) and substitute the gradient of the initial wave function using Eq.(B14), As for the scalar atomic form factor, the exponential in the above equation can be expressed in terms of spherical harmonics via Eq.(B4), where we defined two new radial integrals, which will require extra attention in App.B 5. We again apply Eq. (B5) to explicitly perform the angular integral, It can be useful to evaluate the above expression in a particular coordinate frame, namely the one in which the z-axis points towards q, such that θ q = 0. Then we can use Y M L (θ q = 0, φ q ) = (2L + 1)/(4π)δ M 0 and sum over M , We explicitly checked that this choice of coordinate frame has no impact on the final values of the atomic response functions.
With the scalar and vectorial atomic form factors given in the forms of Eqs.(B22) and (B32), it is possible to evaluate the three new atomic response functions given by the Eqs.(39b) to (39d).However, before doing so, we have to specify the radial components of the initial and final state electron wave functions without which we cannot evaluate the radial integrals defined in Eqs.(B21), (B29), and (B30).

Initial and final state wave functions
For the radial part of the initial state wave function ψ n m (x), we assume Roothaan-Hartree-Fock (RHF) ground state wave functions expressed as linear combinations of Slater-type orbitals, The coefficients C j n , Z j , and n j are tabulated for the atomic orbitals of argon and xenon in [42] together with their respective binding energies E n B .Furthermore, a 0 denotes the Bohr radius.
The final electron state is described by a positive energy continuum solution of the Schrödinger equation with a hydrogenic potential −Z eff /r [12,59].
Here, 1 F 1 (a, b, z) is the so-called Kummer's function, or confluent hypergeometric function of the first kind.Note that both radial wave functions are purely real.It should also be noted that the factor of 1/ √ V involving the volume always cancels due to the application of the integral operator in Eq. ( 27), or equivalently with the factor of V in Eqs.(39a) to (39d).
The factor Z eff is determined by matching the binding energy E n B of the ionized orbital [13,14,20], With the definition of the radial parts of the wave functions in place, we can continue with the evaluation of the radial integrals.

Numerical evaluation of the radial integrals
It is fair to say that the largest obstacle to the evaluation of the scalar and vectorial atomic form factors is the computation of the radial integrals 11 The highly oscillatory behavior of the spherical Bessel function can hinder the success of standard numerical integration methods, especially for large values of q.
Below, we outline two of the methods we used to solve this integral.The first one is exact, but requires the numerical integration of an oscillating function, the second one is approximate, but fully analytical.
The python code we used for the evaluation of the radial integrals, the atomic form factors, and the atomic response functions is publicly available [39].
a. Numerical solution One of the most straightforward ways to evaluate the radial integrals is the use of specific numerical integration methods of numerical libraries like e.g.NumPy [52], or computer algebra systems such as Wolfram Mathematica [55].The numerical integration method "DoubleExponential" of the Wolfram language yielded reliable results.However, for large values of the momentum transfer, the method becomes critically slow.
It is usually benefical for the performance of the numerical integration to split up the integration domain, with ∆r = a 0 and truncate the sum at a finite i after the result has converged to a desired tolerance.Since the integrand in our case is highly oscillatory and this method does not ensure that we integrate from one root to the next, it is possible that an accidental cancellation of positive and negative contributions over one sub-interval of the domain can be misinterpreted as convergence.To avoid this subtle problem, we always verify the series' convergence on two consecutive terms.b.Analytic solution The second method effectively replaces the integral with an infinite sum, which, if it converges sufficiently fast, yields accurate and prompt results.This method worked best when applied to large values of q, i.e. exactly where the numerical integration is problematic. 11These types of integrals are known in the mathematical literature as spherical Bessel and Hankel transforms.
The Kummer function (or confluent hypergeometric function of the first kind), which appears in Eq. (B34), can be written as a power series, with x (s) ≡ x(x + 1) . . .(x + s − 1) being the rising factorial.In combination with the analytic form of the initial wave function in Eq. (B33), we find that the above radial integrals are infinite sums of simpler integrals of the form where the values of the parameters α and β depend on whether we evaluate I 1 (q), I 2 (q) or I 3 (q).Importantly, this integral can be solved analytically (provided that (L + α) > −1), where the hypergeometric function 2 F 1 (a, b, c, z) appears which is defined via Provided that the sum in Eq. (B40) converges after a finite number of terms, we can compute the radial integrals analytically.For large q the convergence occurs quickly.For lower values, we run into numerical precision problems, since we have to include very large numbers of terms of (B40).In this case, the numerical approach described in the previous section performs better.
As an example, by using the analytical method described here, the first radial integral, I 1 (q), can now be written as follows where we defined The sum can be truncated when the series converges to a given tolerance, provided that it converges sufficiently quickly.

Appendix C: Scattering kinematics
In this appendix, we summarize a few basic kinematic relations for DM-electron scatterings.
The initial and final energies of the inelastic scattering process are given in Eqs. ( 4) and (5).By the conservation of energy one can show that From this expression, it is clear that for a given momentum transfer q, the minimum DM speed necessary for an electron ionization with ∆E 1→2 = E B + k 2 /(2m e ) can be obtained by setting cos θ qv = 1, Since the local DM velocity distribution of the galactic halo is assumed to have a maximum speed v esc , DM particles are only able to ionize a bound electron with binding energy E B if v min < v max = v esc + v ⊕ .Thereby, we find the range of momentum transfers which can contribute to the ionization, From the requirement that v min < v max = v esc + v ⊕ , we can also derive the maximum final momentum of the electron for a given momentum transfer q, k < 2m e v max q − q 2 2m χ + E B , (C6) which were shown as white lines in Figs. 3 and 4.
Finally, there are a few identities involving v ⊥ el , which are necessary for the evaluation of the DM response functions in the Eqs.(38).We can write v ⊥ el as Hence the square of the v ⊥ el (with k = 0) is given by We will also need the scalar product with the momentum transfer q, For the last two relations, we used energy conservation via Eq.(C1).The dependence on the binding energies via ∆E 1→2 also explains why the DM response functions R n i carry the atomic quantum numbers (n ) as indices.
Appendix D: Anapole, magnetic and electric dipole dark matter The anapole, magnetic dipole and electric dipole DM models are described by the following interaction Lagrangians where χ (ψ) is a Majorana (Dirac) spinor describing the DM particle, g a dimensionless coupling constant and Λ a mass scale.For simplicity, we use the same notation for the coupling constants of the three interactions, but they can in principle be different.The Lagrangians in Eq. (D3) induce a DM-electron coupling via the electromagnetic current, J ν , and Maxwell equations, ∂ µ F µν = eJ ν , where ψ e is the four component electron spinor, F µν = ∂ µ A ν − ∂ ν A µ is the photon field strength tensor, and A ν the photon field.Taking the non-relativistic limit of the free field solution of the equations of motion for ψ (or χ), and using Gordon's identity to express the matrix element of the current J ν between single particle electron states in terms of electromagnetic form factors, at leading order we find the following amplitudes for DM scattering by free electrons + g e ξ †s S χ ξ s • i q m e × ξ †λ S e ξ λ , (D4) M magnetic = eg Λ 4m e δ s s δ λ λ + 16m χ m e q 2 iq • v ⊥ el × ξ †s S χ ξ s δ λ λ − 8g e m χ q 2 q • ξ †s S χ ξ s q • ξ †λ S χ ξ λ − q 2 ξ †s S χ ξ s • ξ †λ S χ ξ λ , (D5) Here, the notation is the same one used in the main body of the article.Matching these expressions on to the effective theory expansion in Eq. ( 19), one can express the effective coupling constants given in Sec.IV B in terms of g and Λ.In the numerical calculations, we set g e = 2.
Appendix E: Details of direct detection and exclusion limits The exclusion limits presented in this paper are based on a re-interpretation of published data and results by DarkSide-50 [12], XENON10 [38], and XENON1T [15].In this section, we briefly summarize the necessary details to compute the exclusion limits.
An essential input for the prediction of ionization rates is the local density and velocity distribution of the DM particles of the galactic halo, which first enters Eq. ( 13).The local DM number density is simply given by where we set ρ χ = 0.4 GeV/cm 3 [60].For the velocity distribution, we use the standard halo model (SHM) which approximates the distribution by a truncated Maxwell-Boltzmann distribution that has been boosted into the Earth's rest frame, Here, N esc ≡ erf(v esc /v 0 ) − 2(v esc /v 0 ) exp(−v 2 esc /v 2 0 )/ √ π is a normalization constant, and the standard choices for the other parameters are v 0 = 220km sec −1 for the Sun's circular velocity [61], and v esc = 544km sec −1 for the galactic escape velocity [62].For the speed of the Earth/the observer in the galactic rest frame, we choose v ⊕ ≈ 244km sec −1 .
In all cases, we have to translate the energy spectrum given in Eq. (30)  In modelling the probability P(n e |E e ) for n e electrons to reach the gas phase of the TPC given an initial electron energy of E e , we follow [13,14].
For the three experiments, we obtain the 90% C.L. exclusion limits by applying Poisson statistics to each event bin independently.The used parameters for the distribution's mean and width are the secondary-scintillation gain factor g 2 = 27 (33) and the associated width factor σ S2 = 6.7 (7) for XENON10 (XENON1T) [14,15,65,66].For XENON10, the binned observed numbers of events are given on the left hand side of Tab.II, which correspond to an exposure of 15 kg days.The efficiency ε(S2) is given by the product of a flat cut efficiency of 92% [38] multiplied by the trigger efficiency given in Fig. 1 of [13].
In the case of XENON1T, we extracted the number of observed events passing all the cuts from Fig. 3 of [15].They are listed on the right side of Tab.II.The efficiency ε(S2) is the product of a flat efficiency of 93% and the different efficiencies given in Fig. 2 of [15] 12 .The exposure E is given by where the target's radius is R = 47.9cm, the height of the search's volume is ∆z = 20cm (corresponding to z ∈ [−30cm, −10cm]), the density of liquid xenon is taken to be ρ Xe = 3.1g cm −3 , and the time of the search data ∆t = 180.7 days [15].This yields an exposure of ∼ 80755 kg days.Compared to the official XENON1T limits on the standard DM-electron scattering cross section (corresponding to our effective coupling c 1 ), our constraints are weaker by a factor of a few and therefore conservative, which can potentially be explained by the additional background subtraction the XENON collaboration performs.

DarkSide-50
Just like XENON10 and XENON1T, DarkSide is a dual-phase TPC.However, instead of xenon, it uses ar-gon as target [20].The DarkSide-50 experiment had an exposure of 6786 kg days, and a threshold of n e = 3 electrons.The number of observed events we base our limits on have been extracted from Fig. 3 of [20].The official constraints on the standard DM-electron interaction (corresponding to the effective coupling c 1 ) turn out to be stronger than our own limits by a factor of a few.

FIG. 3 .
FIG. 3. The four atomic responses for the outer atomic orbital of argon.Above the white dotted/dashed/solid line on the left plots, the minimum speed vmin in Eq. (34) exceeds the maximum speed vmax = v⊕ +vesc for a DM mass of 10 MeV/100 MeV/→ ∞ respectively (see App. C and E).The top left panel shows the same function as the left panel of Fig. 2, but now with the final state electron asymptotic momentum k on the y-axis.Note the different color bar scales in the four panels.

FIG. 4 .
FIG.4.The four atomic responses for the outer atomic orbital of xenon.The white lines have the same meaning as in Fig.3.Again, the function in the top left panel was shown before in the right panel of Fig.2.However, here the final state electron asymptotic momentum k is on the y-axis.

13 ℓ =10 - 3 FIG. 5 .
FIG. 5. Selected ionization spectra for argon for contact and long-range interactions with coupling c6 (top), and for xenon for contact and long-range interactions with coupling c13 (bottom), for three different DM masses.The faint vertical lines indicate the kinematic cut-off of the spectra.

FIG. 6 .
FIG.6.Integrated ionization spectra for individual DM-electron interaction operators divided by the corresponding coupling constant squared.We present results for argon (purple) and xenon (blue), as well as for contact interactions (top) and longrange interactions (bottom).A comparison of the different histograms shows the operators' relative strength when coupling constants are set to the same value.As indicated for the case of the O1 operator, we assume a DM mass of 10/100/1000 MeV for the left/middle/right histograms.

For
XENON10 and XENON1T, the number of observed signals given in bins of the ionization signal (S2) are available and listed in Tab.II.It is necessary to have the spectrum in terms of S2, i.e. the number of photoelectrons (PEs) in the photomultiplier tubes, number of electrons n e , we assume that the resulting number of PEs follows a Gaussian distribution with mean n e g 2 and width√ n e σ S2[14,63,64],P(S2|n e ) = Gauss(S2|n e g 2 , √ n e σ S2 ) .(E5) DM particle with initial momentum p and a bound electron in the |e 1 state scatter to a group of final states with DM momenta in the infinitesimal interval (p , p + dp ) and the electron in the |e 2 state is given by |S f i | 2 /V 4 times the number of states in the (p + dp ) interval, namely 4, where we divide by V 4 in order to obtain unit-normalized single-particle states from Eq. (1) (remember that S f i has dimension[energy] −6 ).Consistently, the probability, P(p), that a

TABLE II .
Events at XENON10 (left) and XENON1T (right) observed in the given S2 bins.