Solar Active-Sterile Neutrino Conversion with Atomic Effects at Dark Matter Direct Detection Experiments

The recent XENON1T excess can be explained by the solar active-sterile neutrino conversion with bound electrons via light mediator. Nevertheless, the atomic effects are usually omitted in the solar neutrino explanations. We systematically establish a second quantization formalism for both bound and ionized electrons to account for the atomic effects. This formalism is of great generality to incorporate various interactions for both neutrino and dark matter scatterings. Our calculation shows that the change in the cross section due to atomic effects can have important impact on the differential cross section. It is necessary to include atomic effects in the low-energy electron recoil signal at dark matter direct detection experiments even for energetic solar neutrinos. With the best-fit values to the XENON1T data, we also project the event rate at PandaX-4T, XENONnT, and LZ experiments.

In addition, light neutrinos can also provide an explanation. Especially, the Sun is a natural source of energetic neutrino fluxes around Earth. While the SM neutrino interactions can only produce a flat event spectrum, a low energy peak arises if there is a BSM interaction with light mediator. Possible realizations include magnetic dipole moment, charge radius and anapole [3,[46][47][48][49][50][51][52][53][54], all with the massless photon as mediator. In addition, light scalar [55][56][57] and vector [55][56][57][58][59][60][61][62][63] mediators can also achieve the same purpose. It is interesting to observe that only scalar and vector mediators have been used to explain the low energy peak while the pseudo-scalar one was claimed to be incapable of achieving the goal [56]. However, if neutrino scattering with electron into a massive sterile neutrino in the final state, a light pseudo-scalar mediator can also produce a peak in the low recoil energy [64]. The O(100) keV sterile neutrino can also be applied with dipole magnetic moment interactions to explain the excess [65].
For an O(keV) electron recoil energy, the momentum transfer is comparable with the atomic energy of a heavy element such as Xe. The electrons inside atom can no longer be treated as a free particle. Instead, one should use quantum wave functions to describe the electron distribution. It is the electron cloud rather than a single electron particle that participates in the scattering process. This has been extensively explored for the DM scattering with electron in recent years [66][67][68][69][70][71]. The effect of the initial-and final-state electron wave functions can be summarized into a K-factor [71] or sometimes f -factor [66] as function of the momentum transfer q. The size of atomic effect is typically O(1) and hence not negligible.
Being usually unnoticed, the atomic effect with electron cloud has also been considered in the study of neutrino electromagnetic properties. Since the massless photon is exactly a light mediator, the momentum transfer in neutrino scattering with photon mediation is intrinsically suppressed in the same way as the low energy electron recoil at DM detectors. So the atomic effect cannot be neglected in the study of neutrino electromagnetic properties. For example, the ionization effect due to the neutrino scattering with bound electrons is obtained by considering the binding energy and the initial wave function in light atoms [72,73]. Later, the effect of atomic potential on the final-state electron is also taken into consideration when calculating the cross section [74][75][76]. Although the atomic effect was claimed to be small [76], recent studies realize that it is actually not negligible [77,78]. Even spin effects have been studied very lately [79].
In this work, we study the neutrino scattering with bound electrons into a massive sterile neutrino to explain the observed electron recoil peak at (2 ∼ 3) keV and make projections for the future DM experiments. In Sec. II, we first summarize the formalism of calculating the neutrino scattering with both free and bound electrons in the language of second quantization. Especially, the bound electron is directly quantized using annihilation and creation operators without involving the inappropriate plane waves or equivalently momentum eigenstates. Sec. III shows the cross section of neutrino electron scattering including the atomic effect and compare it with the free electron scattering case. The results are further used in Sec. IV to fit the observed data with χ 2 minimization to constrain the parameter space and the prospect of detecting such signal at future experiments is projected in Sec. V. Finally we summarize and conclude in Sec. VI.

II. NEUTRINO-ELECTRON SCATTERING WITH ATOMIC EFFECTS
As mentioned above, the atomic effects need to be considered in order to make a precise study of the electron recoil signal from DM or neutrino scattering. The initial-state electron is bound inside the atom instead of being free. In contrast, the final-state electron is knocked out of the atom and becomes ionized leaving a recoil signal in the detector. With a recoil energy around (2 ∼ 3) keV, the ionized electron is not completely free but is subject to the atomic Coulomb potential. One needs to consider the atomic effect for both the initial-and final-state electrons. We first try to establish a unified second quantization description of the bound and ionized states in Sec. II A and then use it to calculate the atomic K-factor in Sec. II B. This formalism of second quantization can accommodate general interactions and apply not just for neutrino but also DM scatterings.

A. Second Quantization of Bound and Ionized Electron States
An electron trapped in the Coulomb potential is no longer a free particle and hence cannot be described by plane wave with fixed momentum. Instead, the conserved variable is the energy eigenvalue including both kinetic energy and Coulomb potential. To be concrete, the bound electron field e B (x) is in general a function of spatial coordinates [80], where a nlm (a † nlm ) is the creation (annihilation) operator for the bound state with principal (n), angular (l), and magnetic (m) quantum numbers. Since positron never enters our discussion, we can safely omit b † nlm in (1) from the beginning. When acting on the vaccum state |0 , the creation operator a † nlm gives the bound state |nlm ≡ a † nlm |0 . There is no need to involve the √ 2E prefactor that is necessary for a relativistic particle to keep Lorentz covariance but reduces to a constant √ 2m and hence can be combined into normalization for a non-relativistic particle. Similarly, we follow the convention of second quantization to define the anti-commutation relations, {a nlm , a † n l m } = δ nn δ ll δ mm and {a nlm , a n l m } = {a † nlm , a † n l m } = 0. With dimensionless discrete δ functions, the creation and annihilation operators are also dimensionless, [a nlm ] = [a † nlm ] = 0, and the bound electron field has the same dimension as its wave function, [e B (x)] = [ψ nlm (x)]. The normalization condition of the wave function, ψ † nlm (x)ψ n l m (x)d 3 x = δ nn δ ll δ mm , fixes the field dimension to be [e B (x)] = 3/2 which is consistent with quantum field theory (QFT). Although the momentum integration is replaced by a summation nlm , the second quantized field e B (x) is still a linear combination of energy and angular momentum eigenstates. However, the second quantized field e B (x) in (1) is different from the usual formalism of a free particle in quantum field theory [81] with the evolution phase e −iE nl t containing only the energy eigenvalue and time dependence while the spatial dependence in the electron wave function ψ nlm (x) cannot be factorized out as a simple complex phase.
For an ionized electron that is still under the influence of the atomic Coulomb potential, its energy is continuously distributed. The corresponding ionized electron field contains an integration over the asymptotic momentum, without involving a principal quantum number. The asymptotic kinetic energy T r ≡ |p| 2 /2m e is the one that we can experimentally measure. Comparing with the bound case in (1), the only difference is that the energy eigenvalue changes from a discrete principal quantum number n to a continuous T r or equivalently the asymptotic momentum |p|. Similar to the summation over discrete variables, the asymptotic momentum is integrated. In other words, the second quantized field e I (x) is also a linear combination of energy and angular momentum eigenstates. With oneto-one correspondence between the asymptotic energy T r and the asymptotic momentum |p|, the latter is also a well defined physical quantity. The ionized electron behaves essentially as a free particle when |x| → ∞ and its wave function reduces to a plane wave [82]. The direction of the asymptotic momentum can also be used to label an ionized electron in the similar way as its magnitude. The Fourier transformation for the asymptotic state at infinity distance (|x| → ∞) is, However, the DM direct detection experiments cannot tell the directional information and are sensitive to only the magnitude |p| or the recoil energy T r . Namely, the relevant wave function is ψ Tr (x) ≡ dΩ p ψ(p)e −ip·x and the remaining phase space integration is exactly the |p| 2 d|p| used in (2), ψ(x) = |p| 2 d|p| (2π) 3 ψ Tr (x). The original anti-commutation relation {a plm , a † p l m } = (2π) 3 δ ll δ mm δ (3) (p − p ) reduces to {a Trlm , a † T r l m } = (2π) 3 δ ll δ mm δ(|p| − |p |)/|p| 2 after integrating away the angular information of the asymptotic momentum p. Correspondingly, the operator dimension is, [a Trlm ] = −3/2 and the wave function ψ Trlm (x) for ionized state is dimensionless. This is consistent with the wave function normalization, ψ † Trlm (x)ψ T r l m (x)d 3 x = (2π) 3 δ ll δ mm δ(|p| − |p |)/|p| 2 [68]. Using the non-relativistic dispersion relation, T r = |p| 2 /2m e , one can prove that, [δ(|p| − |p |)/|p| 2 ] × |p| 2 d|p| = δ(T r − T r )dT r and the normalization condition be- [82] to make everything consistent. Similar to the bound state, an ionized state |T r lm ≡ a † Trlm |0 is created from vacuum by a † Trlm without involving the energy prefactor √ 2E. The wave function can be constructed from the field, 0|e B (x)|nlm = ψ nlm (x)e −iE nl t for the bound state and 0|e I (x)|T r lm = ψ Trlm (x)e −iTrt for the ionized one. In addition to the spatial distribution, the wave function should also contain the spin information. The spinor wave function ψ N (x) is a solution of the covariant Dirac equation in the presence of electromagnetic field, i / ∂ − m − e / A ψ N (x) = 0, where N ≡ (nlm) and N ≡ (T r lm) for the bound and ionized states, respectively. For a stationary energy eigenstate, ψ N (x) = e −iE N t ψ N (x), the time dependence can be removed from the Dirac equation, For clarity, we have only kept the electric potential V (|x|) in the atom as a central force. In the non-relativistic limit, the spatial and spin parts separate into [80], for spin s and other quantum numbers N . The information of spatial distribution and spin is represented by the single-valued wave function φ N (x) and the two-component spinor ξ s , respectively. The normalization condition can be rewritten in terms of For a heavy element such as Xenon, the binding energy is typically O(1 ∼ 10) keV. According to the Virial theorem, the electron kinetic energy is of the same size, which is roughly 0.1% ∼ 1% of the electron mass, m e ≈ 511 keV. In other words, the effect of special relativity or spinor structure represented by the gradient expansion is roughly 0.1% ∼ 1%.

B. Atomic Effects in Neutrino-Electron Scattering
With the second quantization formalism of the bound and ionized states established, we can follow the usual QFT calculation of transition rate and differential distributions. For a general four-fermion coupling (νΓ ν ν s )(ēΓ e e) between neutrino and electron, the transition matrix element T is, with φ generally denoting a mediator with mass M and p ν (p s ) is the solar (sterile) neutrino momentum. For generality, Γ ν , Γ e ≡ 1, γ 5 , γ µ , γ 5 γ µ , σ µν denote all the possible Lorentz structures. Correspondingly, the mediator φ can be a scalar (S), pseudo-scalar (P), vector (V), axial-vector (A), or even tensor particle. The T matrix element should also contain coupling constants which are omitted for simplicity. Note that the tensor case is included just for completeness when illustrating the atomic effects but not discussed when being applied to the realistic case of sterile-active neutrino conversion. This is because the tensor current is typically mediated by a spin-2 particle such as the graviton that is beyond the scope of this work. The other types with scalar, pseudo-scalar, or vector mediator will be elaborated in the coming Sec. III and Sec. IV.
When acting field operators on the initial |nlm, p ν ≡ a † nlm a † pν |0 and final |T r lm, p s ≡ a † Trlm a † ps |0 states, the a nlm operator from e B (y) annihilates a bound state while the a † Trlm operator fromē I (y) creates an ionized state at position y. The T matrix becomes, after integrating away first the coordinate x and the momentum transfer q. The x integration produces a δ-function of four momentum, δ (4) (p ν −p s −q), to impose energy momentum conservation on the neutrino vertex, q = p ν − p s . From the bound and ionized electron fields, one can extract an explicit energy dependence, ∆E nl ≡ T r − E nl , as the exponential factor e i∆E nl t . This allows imposing energy conservation on the electron vertex by integating away the time component y 0 , Together with the δ (4) (p ν − p s − q) function on the neutrino side, energy is conserved both locally and globally while the momentum conservation only applies to the neutrino vertex. This is the key difference between the calculations of the scattering process with free or bound electron. For free electron, the spatial dependence of their wave functions can also factorize out as exponential factors and the integration d 4 y can impose both energy and momentum conservations. The physical reason behind this is that a bound electron does not have a definite momentum, especially in the coordinate representation. Consequently, we can only get a product of the initial and final wave functions.
Note that the δ E ≡ δ(∆E nl − E ν + m 2 s + |p ν | 2 + |q| 2 − 2 |p ν q| cos θ qv ) function in (7) for energy conservation can be moved outside of the d 3 y integration since it only depends on energy and momentum. The spatial integration with bound and ionized wave functions, A e ≡ d 3 ye iq·y ψ Trl m (y)Γ e ψ nlm (y), is essentially the source of the atomic K-factor. Nevertheless, the wave functions ψ nlm and ψ Trl m have spinor ξ that needs to be singled out in order to obtain the M matrix element. For a scalar interaction (Γ e = 1), the electron amplitude A e divides into two parts, when expanded to the linear order of i∇f N,s . The two-component spinor ξ that is momentum independent has been reexpressed,ū(m e )u(m e ) = 2m e ξ † ξ, in terms of the electron spinor u(m e ) at rest. The spinor and spatial wave function then factorize into two parts, is the so-called atomic form factor [66,68,71], The factorization of the scattering amplitude into two parts, one for the spinor and the other for the atomic form factor, is a generic feature when expanding to the leading order. For all the possible electron bilinears, The key feature here is that the initial-and final-state electrons are not momentum eigenstates or free-particle solution of the Dirac equations. Instead, the bound and ionized spinors are originally a function of the spatial coordinates (4). There is no need to involve initial and final momentum for electrons at all. The only momentum that can appear is the momentum transfer q. Especially for the pseudo-scalar case, the electron spin or equivalently γ i with a spatial index is relevant. However, only the inner product q·γ can appear as a combination to keep the scalar nature of this vertex. This factorization is derived with second quantization for the bound and ionized electron wave functions in a systematic way.
With A e factorization, the T matrix element (7) is composed of several contributions, T ≡ Mf Trl m nlm (q)δ E , where M is the scattering matrix element, Note that M is a function of the incoming neutrino momentum p ν and the momentum transfer q while the sterile neutrino momentum p s = p ν − q is a dependent variable. Following the derivations in the Section 4.5 of [81], the scattering cross section with a bound electron can be also expressed in terms of M, for a single electron target. The prefactor 1/(2l + 1) accounts for the degenerate electrons with m = −l, −l + 1, · · · , l − 1, l while |M| 2 ≡ |M| 2 /2 is averaged over the electron spin. For the high energy solar neutrino that contains mainly the left-handed neutrinos, there is no need to average over the neutrino spin. Since the asymptotic kinetic energy T r = |p| 2 /2m e is the one that can be experimentally measured, it is much more convenient to express the phase space integration The original phase space integration for the sterile neutrino momentum p s has also been replaced by the momentum transfer q due to momentum conservation of the neutrino vertex. In addition, the combined neutrino and atom system has rotational invariance around the incoming neutrino momentum p ν . So the azimuthal angle of q can be integrated away to give 2π. The zenith angle θ qν integration is reduced by the δ E to δ E d cos The allowed range of the scattering angle, cos θ qν ≤ 1, gives the |q| integration range, Putting things together, the differential cross section of the recoil energy T r becomes, The atomic effect can now be parameterized as the so-called K-factor, The 1/(2l + 1) factor accounts for the average over the initial magnetic quantum number m. Usually, the direct detection experiments can only distinguish the energy deposit but not the angular quantum numbers m for the initial-and {l , m } for the final-state electrons. So these quantum numbers are summed over when defining the K-factor in (14). For completeness, we rewrite the K-factor more explicitly [67], 1. The left panel shows |q|K nl (Tr, |q|) (dashed) and the phase space ratio |q|K nl (Tr, |q|)/meTr (solid) for the 2p, 3p, 4p and 5p orbitals from top to bottom. Three typical recoil energies Tr = 5 keV (black), 10 keV (blue), and 20 keV (green) are shown for illustration. The peak locations, |q| = √ 2meTr, estimated from the scattering with a free electron are shown as vertical gray lines. As a function of the electron recoil energy Tr, the right panel shows the phase space ratio |q|K nl (Tr, |q|)d|q|/meTr for the 2p (black solid), 3p (red dashed), 4p (blue dotted) and 5p (green dash-dotted) orbitals in the Xe atom with solar neutrino energy Eν = 300 keV and the sterile neutrino mass ms = 100 keV.
with the wave functions divided into radial and angular parts, φ nlm (r) = R nl (r)Y lm (r) and φ Trl m (r) = R Trl (r)Y l m (r). For convenience, the spatial integration variable y has been replaced by r and r is radius. The Bessel function j L (|q|r) originates from the exponential e iq·y of (9). After integrating over the solid angles, only the radial functions [82] R nl for the bound and R Trl for the ionized electrons survive. Most importantly, the K-factor becomes no longer a function of q but its magnitude |q|.
Note that the normalization of K-factor varies in literature [66,68,71,83]. In our definition (15), the summation over m and m gives a factor of 2l + 1 and 2l + 1, respectively. That 2l + 1 factor has been canceled with the one in (14). With this convention, (15) defines the average K-factor for a single electron in the nl state and (13) is the cross section per electron. In principle, the final-state angular quantum number l should sum up to infinity. However, in practice the summation is cut off at sufficiently large l ∼ O(100) and higher contributions are neglected [71,83] To see the basic features of the K-factor more clearly, let us compare the scattering with free and bound electrons. If the electron is not tightly bound in the atom, or equivalently n → ∞, the whole process should reduce to the scattering with a free electron. The differential cross section of the scattering with a free electron is, with |p| 2 = 2m e T r . For easy comparison, we have kept the momentum transfer integration which can be removed by the δ(|p| − |q|) function. Comparing (16) with (13), we can see that the Kfactor reduces to |q|δ(|p| − |q|)/2 in the free electron scattering limit. For the scattering with a free electron, the two-body phase space reduces to a single integration over either the recoil energy T r or equivalently the momentum transfer |q|. But for the bound electron case, the phase space has double integration, over T r and |q| which are no longer correlated. This is because the initial bound electron does not have definite momentum but a distribution. Consequently, the contribution from  different momentum transfer should be integrated. Besides |M| 2 /32πm 2 e |p ν | 2 T r , the phase space integrations are m e T r and |q|K nl (T r , |q|)d|q| for the free and bound electron cases, respectively. Assuming constant |M| 2 , the atomic enhancement can be roughly measured by the ratio between the phase space integrations, |q|K nl (T r , |q|)d|q|/m e T r .
The left panel of Fig. 1 shows the atomic factor |q|K nl (T r , |q|) (dashed) and the atomic factor ratio |q|K nl (T r , |q|)/m e T r (solid) as functions of the momentum transfer |q|. The curves are narrower for outer shell electrons (larger n) and the typical width keeps growing when n becomes smaller. The half width at the half height σ |q| = (50, 20, 10, 4) keV for the (2p, 3p, 4p, 5p) electron can be directly read off from the solid lines. The solid and dashed lines share the same width since the only difference between them is a constant factor. For given electron shell, the width is independent of T r . The width is exactly a manifestation of the electron motion inside atom. Pointing in all directions, the initial electron momentum smears the momentum transfer |q| and hence the recoil energy T r . The larger initial momentum, the larger smearing effect. With binding energy |E b | = (4.53, 0.95, 0.16, 0.012) keV for (2p, 3p, 4p, 5p) orbitals [71], the initial electron momentum can be estimated by non-relativistic dispersion relation, |p i | ≈ 2m e |E b | =(68, 31, 13, 3.5) keV correspondingly. These numbers roughly match the width read off from the curves. Although the bound electron motion broadens the curves, the central values are largely unaffected, especially for the outer shells. The peak position can be estimated by the momentum transfer |q peak | ≈ √ 2m e T r of the free electron scattering. For the three curves of T r = (5, 10, 20) keV, the estimated peak positions are |q peak | ≈ (71, 101, 143) keV as indicated by three vertical grey lines from left to right. With larger n, the peak position becomes closer to the free electron momentum transfer.
The right panel of Fig. 1 shows the integrated phase space ratio as functions of the recoil energy T r . The ratio starts from a relatively small value for vanishing T r and converges to 1 with increasing T r . This happens because for small T r the allowed momentum transfer is suppressed to have smaller phase space [66]. In contrast, the electron approaches a free particle to have larger phase space with increasing T r . For larger n, the integrated phase space ratio approaches 1 faster, since the electron is less tightly confined to the atom.
The left panel of Fig. 2 shows the non-zero terms inside the double summation of (15) as a function of the angular indices (l , L) for the 2p, 3p, 4p, and 5p orbitals. Those terms violating |l ± L| = l , and 5p. The first row is purely for the initial electron radial wave function r 2 |R nl | 2 which appears in the integration for normalization. The next three rows are for the convolution between the initial and final electron radial functions r 2 R * nl R Tr l (blue solid) and r 2 R * nl R Tr l jL(|q|r) with L = 0 (red dashed) for electron recoil energies Tr = (5, 10, 20) keV. For illustration, we take the peak value |q| = √ 2meTr of the momentum transfer. The horizontal axises are in the unit of Bohr radius a0 = 1/meα. are zero and omitted in the plot for simplicity. With increasing n, the peak shifts to the right. The 2p curve peaks around l = 5 and drops to zero at l = 15. For comparison, the 5p curve peaks around l = 95 and is non-zero even for l > 190. The exact evaluation of the K nl factor in (15) requires double summation over the principal quantum numbers, l and L, to infinity which is highly time-consuming. To increase efficiency, we cut the summation at large enough l max and L max to guarantee precision.
The relative precision can be estimated from the right panel of Fig. 2. Each curve shows the ratio between the partial sum up to (l max , L max ) divided by the one up to (200, 201) for the 2p, 3p, 4p, and 5p orbitals. With increasing (l max , L max ), the ratio converges to 1. As n increases, larger l max and L max are needed. For the 2p, 3p, 4p, and 5p orbitals to achieve sub-percentage precision, (l max , L max ) needs to be at least (10,11), (30,31), (70,71), and (190, 191), respectively.
The integrand in (15) is highly oscillatory and illustrated in Fig. 3. The top panels show purely the initial bound electron wave functions as a convenient combination r 2 |R nl | 2 (black solid) that should appear in normalization integration. While the initial wave functions are quite regular with only one major peak, the final ionized electron wave functions oscillate a lot. This oscillatory behavior clearly appears in the product combination r 2 R * nl R Trl (blue solid). Since the bound and ionized states are orthogonal to each other, the product oscillates symmetrically around the vanishing value so that the integration should also vanish. The third component, the Bessel function j L (|q|r), enters the atomic form factor to render a nonzero K-factor (15). The full combination r 2 R * nl R Trl j L (|q|r) (red dashed) has asymmetric oscillation so that its integration can be nonzero. With larger n, l, and T r , the oscillation becomes more frequent. To accurately perform the integration in (15), a large number of grid points over r is necessary. A sub-percentage precision requires 3000 grid points in the range r ∈ [0, 10a 0 ] where a 0 ≡ 1/m e α is the Bohr radius.

III. NEUTRINO SCATTERING INTO A MASSIVE STERILE NEUTRINO WITH LIGHT MEDIATORS
With the second quantization formalism for the atomic effect established in the previous section, we can proceed to calculate the recoil spectrum (13) by implementing the concrete scattering matrix element |M| 2 as input from the particle physics side. For the solar neutrino scattering into a massive sterile neutrino in the final state, we first derive the differential cross section in Sec. III A and then show the modification due to atomic effects in Sec. III B.

A. Neutrino Scattering with Free or Bound Electrons
The XENON1T excess in the low energy region can be explained by the solar neutrino scattering with electron via a scalar or pseudo-scalar mediator [64].
In addition to the recoil electron, a massive sterile neutrino appears in the final state. For the scattering with a free electron, the cross section is a function of the recoil energy T r , Pseudo-Scalar : In addition to the neutrino energy E ν , m s and m φ are the sterile neutrino and scalar/pseudo-scalar mediator masses, respectively. Contrary to the common expectations [56] with massless neutrinos, a massive sterile and a light mediator (m 2 φ 2m e T r m 2 s ) can introduce 1/T r enhancement even for the pseudo-scalar mediator [64] by a factor of m 2 s /2m e T r . For a low energy electron recoil signal, such as the XENON1T excess at T r ≈ (2 ∼ 3) keV, the mediator mass should satisfy an upper bound m φ (45 ∼ 55) keV to receive enhancement. For the scattering with a bound electron, the electron part resembles the free case for the scalartype vertex and receives a momentum insertion for the pseudo-scalar one as shown in (11), When calculating the trace, we should remember that the spinor for electron u(m e ) is the one without momentum. But it does not mean this approximation has no information of the electron momentum. Actually, the electron momentum appears as the gradient operator in (4), which when acting on the integration will generate the momentum transfer q. This is also the origin of q · γ in the pseudoscalar matrix element. Putting the scalar and pseudo-scalar cases together, the spin-averaged matrix element is, |M S,P | 2 = 2 4m 2 e |y ν S y e S | 2 + |q| 2 |y ν P y e P | 2 |q| 2 + m 2 Although derived for the bound electron case, (20) can reproduce the leading terms of the free electron case (18) with the following replacements. The initial electron is at rest and the energy difference is exactly the recoil kinetic energy, ∆E nl = T r , of the final electron. In addition, the momentum transfer magnitude |q| is uniquely related to the recoil energy, |q| 2 = T 2 r + 2m e T r , due to the on-shell condition of the final-state electron.
On the other hand, a bound electron in the initial state has no definite momentum. There is no way to use the on-shell condition for the final-state ionized electron to correlate the energy change ∆E nl ≡ T r − E nl and the momentum transfer magnitude |q|. The differential cross section (13) is then an integration over all the possible momentum transfers, Pseudo-Scalar : The similar scenario with a light vector boson mediator can also explain the XENON1T excess [64]. For simplicity, we consider only the situation where Z couples to the left-handed neutrino, while the electron coupling can have either vector (g e V ) or axial-vector (g e A ) current coupling with Z . The differential cross section of neutrino scattering with a free electron is, where the + (−) signs are for the vector and axial-vector interactions, respectively. The sterile neutrino mass term m 2 s and hence the second term in the numerator can be important only when m s becomes comparable with the electron mass m e and the neutrino energy E ν . In this paper, we focus on the mass region m s ∼ O(100) keV. Then with 4m e E 2 ν dominating the numerator, a 1/T 2 r enhancement is possible for a light enough mediator, m 2 Z 2m e T r . For comparison, the differential cross section of neutrino scattering with a bound electron is, with ∓ for the vector and axial-vector couplings, respectively.

B. The Cross Section Enhancement from Atomic Effects
By implementing the K-factor elaborated in Sec. II B, we obtain the differential cross sections (21) and (24) for scalar and vector mediators. To see the features clearly, we first evaluate dσ nl /dT r for the transition of a single initial electron in a given nl bound state to the ionized electron with recoil energy T r , shown as colorful lines in Fig. 4. For comparison, we also show the free electron case as black solid line. We take the scalar and pseudo-scalar mediators to illustrate the physics picture while the vector case is similar. The difference is quite significant. Especially, the smearing effect makes the sharp peak much fatter and lower mainly due to the electron motion in the atom. For smaller principal quantum number n, the reduction is much stronger. One reason is the phase space ratio due to the atomic K-factor as shown in Fig. 1. Another factor comes from the scattering matrix element. The electron in the atomic Coulomb potential has a negative binding energy E b which is exactly E nl . Intuitively, It is more difficult for an electron with larger binding energy to be recoiled off. If the energy change ∆E nl = T r − E nl in (20) dominates, especially for the inner shells, the scattering matrix element can be greatly suppressed. This is because the momentum transfer peaks at |q| = √ 2m e T r with large spread above the peak position as shown in Fig. 1. So the modification of the differential cross section is a combined result of the K-factor phase space integration reduction and the scattering matrix element suppression. On the other hand, the electron in outer shells is loosely trapped with marginal suppression and is very similar to the free case.
In addition, the smearing effect extends the spectrum beyond the cut-off. For free electron scattering, the recoil energy is bounded from below, T r ≥ T − r as defined in our earlier paper [64]. A nonzero lower limit T − r ≈ m 4 s /8m e E 2 ν arises when expanding the sterile neutrino mass to the fourth power. For the E ν = 300 keV used to plot the black solid line, the differential cross section vanishes around 1.7 keV. However, bound electrons have distributed momentum in all directions and hence can smear the momentum transfer to extend the recoil energy T r even down to 0 keV. Furthermore, the 2p or 3p electrons even have a non-zero differential cross section at vanishing T r . As shown in (13), the phase space integration is |q|d|q|. Since the lower limit E ν − (E ν − ∆E nl ) 2 − m 2 s of |q| with ∆E nl > 0 is always nonzero for the scattering of bound electron into an ionized one, the phase space would not disappear even for T r = 0 keV. This behavior can also explain the sudden drop of the 5p curve when approaching vanishing T r in Fig. 1.
Solar neutrino can scatter with electrons of different quantum numbers inside the Xenon atom. For the free electron scattering, the total dσ/dT r is simply N e = 54 times of the differential cross section for a single electron where N e is the electron number in the Xenon atom. The total cross section for the bound electron scattering is a sum over all the differential cross sections dσ nl /dT r , Free : dσ dT r = 54 × dσ 0 dT r , vs Bound : , and the shifted one including recombination energies (solid). The enhancement from the free electron scattering cross section is shown as ratio R in the lower panels for comparison. Since the couplings only appear as overall factor, we adopt |y ν y e | = |g ν g e | = 1 for simplicity.
Of the weight 2(2l + 1), the factor 2 for the bound case accounts for the two electron spin degrees of freedom with the same nlm quantum numbers. Note that the summation over the magnetic quantum number m results in the 2l + 1 factor. As illustrated by the blue lines in Fig. 4, the 4s, 4p, and 4d curves with the same principal quantum number n = 4 but different angular quantum numbers have similar shapes and amplitudes. More generally, orbitals with the same principal quantum number have roughly the same differential cross sections. The n = 5 orbitals have higher and sharper peaks than the n = 4 ones. However, the electron number of the former, 2(5s) + 6(5p) = 8, is less than half of the later, 2(4s) + 6(4p) + 10(4d) = 18. So the total differential cross section of atomic scattering (thick dashed and solid) in Fig. 5 is mainly contributed by the 4s, 4p, and 4d curves (blue) in Fig. 4. Since the differential cross section dσ nl /dT r for a single electron is reduced for all nl quantum numbers, the total result dσ A /dT r is also reduced from the free electron case. As shown in the lower panel, the reduction is roughly a factor of 0.5 in the low-energy region (T r 5 keV) and gradually recovers to 1 in the high energy region (T r 10 keV) for the scalar or pseudo-scalar mediators. For vector and axial-vector interactions, atomic effects bring the same features.
A key difference between the scalar/pseudo-scalar and vector/axial-vector interactions is that the latter has much larger cross section by almost a factor of 10 ∼ 15, as shown in the upper panels of Fig. 5. As discussed below (23), the matrix element part is approximately 4m e E 2 ν /(2m e T r +m 2 Z ) 2 for the free electron scattering with vector mediator. For comparison, the scalar case has m e (2m e T r + m 2 s )/(2m e T r + m 2 φ ) 2 . With tiny mediator mass, m 2 φ , m 2 Z 2m e T r , the difference mainly comes from the numerator part. The ratio between vector and scalar matrix elements is 4m e E 2 ν /m e (2m e T r +m 2 s ). We can see that the vector case has a major contribution 4m e E 2 ν . With E ν ≈ 300 keV, m s ≈ 150 keV, and |q| 2 ∼ 2m e T r ≈ 3000 keV, the vector cross section is naturally enhanced by a factor around 10.
In addition to the free (dotted) and bound (dashed) electron curves, Fig. 5 also shows the shifted results. This is because the energy deposit ∆E nl is not just the electron recoil energy T r but also the binding energy E b = E nl . The later is released when an ambient electron is attracted by the positively charged Xenon atom to fill the hole left by the ionized electron [84]. The differential cross section is then shifted to the right as the solid lines. The binding energy is typically O(0.01) keV for outer shells [84] and can be ignored in comparison to the recoil energy. But the binding energy can be as large as O(10) keV for most inner shells and hence seems able to induce significant horizontal shift, the differential cross section dσ nl /dT r is suppressed by the same large binding energy in the first place as discussed above. So the energy shift does not affect the total differential cross section dσ/dE vis much as shown in the plot.

IV. EXPERIMENTAL CONSTRAINTS
The XENON1T Collaboration has recently observed an excess in the electron recoil spectrum around (2 ∼ 3) keV [2]. The experimental data is shown as black points in Fig. 6 while the total background from radioactive materials present in the detector and solar neutrino scattering is shown as the red curve for comparison. Note that the background is roughly flat and hence can not explain the excess at low energy. It is possible for this excess to arise from some new physics and the data points can be used to constrain the possible new physics model parameters.
In this paper, we give a more detailed evaluation of the constraint on the solar neutrino scattering with electrons into a massive sterile neutrino [64]. The atomic effects are also included to give a realistic result. We use the following χ 2 function, to evaluate the best-fit value and sensitivities. The χ 2 function contains two major contributions: the first term from the data bins and the second is nuisance parameter c for the background normalization with uncertainty σ c = 2.6% [2]. In each bin, N data i is the data point, N bg i the expected background, and N NP i the expected events due to new physics (NP) contribution. For the light mediator case considered in this paper, N NP i is a function of the mediator mass (m φ or m Z ), the sterile neutrino mass (m s ), and couplings (|y ν S,P y e S,P | or |g ν V,A g e V,A |). On the theoretical side, the event rates are calculated by the convolution of the differential crosssection defined in (25) and the solar neutrino flux. We use the pp-neutrinos since they have the largest flux at low energy [85]. The detector's energy resolution is also taken into account by a Gaussian parametrization with variance given by σ Tr = a/ T r /keV + b with a = 31.71 ± 0.65 keV and b = 0.15±0.02 keV [86]. In addition, the detector efficiency at a given recoil energy T r is extracted  FIG. 7. The XENON1T constraints for scalar (up), pseudo-scalar (bottom) interaction couplings |y ν S,P y e S,P | versus ms (left) and m φ (right) at 1σ, 2σ, and 3σ C.L. The black curves are for the free electron approximation and red ones for the bound electron scenario whose best-fit points are indicated with black (free) and red (bound) stars, respectively. The mediator (sterile neutrino) mass m φ (ms) is marginalized over 0 keV < m φ < 175 keV (0 keV < ms < 250 keV).
from [2]. Fig. 7 shows the best-fit values and exclusion curves obtain by fitting the predicted event numbers N NP i with data using (26). The left and right panels are obtained by fixing two parameters, |y ν S,P y e S,P | and m φ (m s ), and minimizing the χ 2 function over the remaining m s (m φ ). For comparison, each panel shows both free (black curves) and bound (red areas) electron scenarios. The bound-to-free ratio R = 0.5 for the scalar mediator in Fig. 5 seems be quite different from 1 and can have significant consequence. However, it is mainly due to the broadening effect of the electron motion inside atom. With detector smearing effect, the sharp peak for the free scattering case would also become broader and lower to more closely resemble the bound case. Consequently, the coupling best-fit value almost does not change. Take the scalar interaction as an example, the coupling best fit marginally shifts from the free electron case (|y e S y ν S | = 5.03 × 10 −13 ) to the bound one, |y e S y ν S | = 5.13 × 10 −13 . Similarly, the coupling strength changes from 5.63 × 10 −12 to 5.33 × 10 −13 for pseudo-scalar mediator. The background-only hypothesis (|y e S,P y ν S,P | = 0) is excluded by more than 3σ C.L. for the scalar interaction but is still inside the 2σ region for the pseudo-scalar case. This is because the pseudo-scalar interaction has a smooth total differential cross section (blue curves in the left panel of Fig. 5) and is not as different from the background as the scalar one.
In addition to the DM direct detection, the coupling with electron y e S,P also receives constraints from other experiments and astrophysical observations. The torsion balance experiments impose a constraint |y e S,P | < 10 −23 at 95% C.L. for the mediator mass below 4 × 10 −6 eV [87], which is 10 orders smaller than the current DM direct detection constraints of O(10 −13 ). In other words, an ultra-light mediator cannot explain the XENON1T excess. The mediator mass up to O(0.01) eV is also excluded by the Red Giant (RG) and Horizontal Branch (HB) stellar cooling constraint   |y e S,P | < 10 −16 [88] at 95% C.L. that is 3 orders smaller than the DM direct detection ones. The Big Bang nucleosynthesis (BBN) constrains the mediator mass above 10 −2 eV. The φ production in the early Universe increases the relativistic degrees of freedom and accelerates the universe expansion. A faster expansion reduces the deuterium abundance. The φ production can only be evaded if |y e S,P | < 5 × 10 −10 at 95% C.L. or m φ > 1 MeV [89]. For the coupling with neutrino, the meson decay experiments require |y ν S,P | 10 −3 if m s < m π [90][91][92][93]. Altogether, the combined constraint is |y e S,P y ν S,P | 5 × 10 −13 in the mass region m φ > 10 −2 eV (the gray lines in Fig. 7). The best-fit value of the scalar case is allowed while the pseudo-scalar one is in tension.
The best-fit value of sterile neutrino mass almost does not change for the scalar interaction, the best fit shifts slightly from m s = 147.4 keV to 128.2 keV after implementing atomic effects. For the pseudo-scalar interaction, the best fit remains almost the same at m s = 102.6 keV. In other words, the dependence on the sterile neutrino mass is not sensitive to the atomic effects. Similar feature also applies for the light mediator mass whose best-fit value is 0 keV in all situations. This is because the propagator contribute a factor ∝ 1/(|q| 2 − ∆E 2 nl + m 2 φ ) 2 in the differential cross section which typically grows with decreasing m φ for given |q| 2 . In other words, a smaller mediator mass leads to a higher peak at lower recoil energy which is preferred by the XENON1T data. However, the light mediator mass cannot really be zero due to various constraints [87,89]. To be on the safe side, we adopt m φ = 10 keV which is still inside the 1σ range when predicting the event rates in Sec. V. A mediator mass of 10 keV avoids the constraints mentioned before. This modification does not make big difference in the direct detection measurement. To be more concrete, the χ 2 only increases by 1 and 1.5 for the scalar and pseudo-scalar cases, respectively. Fig. 8 shows the results for vector and axial-vector interactions. We can see that the basic features are similar to the scalar and pseudo-scalar counterparts. The coupling |g ν V g e V | changes from 1.7×10 −13 to 1.8 × 10 −13 for vector while |g ν A g e A | changes from 1.6 × 10 −13 to 1.5 × 10 −13 for axial vector interactions. This change can also be understood by the combination of bound-to-free cross section ratio and the detector smearing effect in a similar way as the scalar/pseudo-scalar cases above. Both vector and axial-vector interactions can exclude the background-only hypothesis by more than 3σ. The best-fit value of sterile neutrino mass is around 150 keV and the mediator mass still prefers a tiny value. Similarly to the scalar and pseudo-scalar cases, we take m Z = 10 keV below when making projections for the future experiments.
The most stringent model-independent constraint on the coupling with electron, |g e V,A | < 5×10 −10 , also comes from BBN [59]. In the sterile neutrino mass region m s < m π , the leptonic pion decay imposes a bound |g ν V,A | < 0.014 on the coupling with neutrino [92,94]. The combined result is |g e V,A g ν V,A | < 7 × 10 −12 . For both the vector and axial-vector interactions, the best-fit points are well below this constraint. Both cases can escape the constraints. In fact, the bound is one order of magnitude higher than the best-fit points and is not visible in Fig. 8

V. PREDICTIONS FOR FUTURE DM EXPERIMENTS
Although the XENON1T excess can be explained by new physics, the significance is not large enough and there is still no definite conclusion for the DM-electron interaction [2]. Especially, the tritium background is also possible explanation [2,3]. More data is necessary to obtain a conclusive result. In this section, we use the best-fit values from the current XENON1T data to make prediction for future experiments.
We focus on three major DM direct detection experiments with liquid Xenon target. First, the PandaX-4T experiment [95] that has already started running in 2021 is an upgrade of PandaX-II with a fiducial mass of 2.8 t and a factor of 10 improvement in the sensitivity. Next is the XENONnT experiment [96] upgraded from XENON1T. XENONnT has a fiducial mass of 4 t and can reduce its background by a factor of 6. Finally, the LUX-Zeplin (LZ) experiment is a combination of two existing experiments LUX [97] and Zeplin-III [98]. Its fiducial mass can reach 5.6 t and also has very low background. In addition, we make predictions with the same efficiency of PandaX-II [99] for PandaX-4T and XENON1T [2] for XENONnT while for LZ the Fig. 3 9 shows the predicted event rates with atomic effects for future experiments PandaX-4T (red), XENONnT (blue), and LZ (green). The left panel shows the results for scalar (solid) and pseudoscalar (dashed) interactions while the right panel shows the vector (solid) and axial-vector (dashed) cases. In these predictions, the couplings and the sterile neutrino mass m s are assigned the best-fit values with the XENON1T data. But the light mediator masses (m φ and m Z ) take a universal value of 10 keV. With atomic effects, the event rates for scalar, vector and axial-vector interactions have a more conspicuous low energy excess around (2 ∼ 3) keV. The event rate for the pseudo-scalar interaction gradually increases towards low energy. Fig. 10 shows the expected sensitivity of the scalar, pseudo-scalar (left) and vector, axial-vector (right) couplings as function of the sterile neutrino mass m s , for the three future DM experiments PandaX-4T, Lux-Zeplin, and XENONnT. Their nominal exposure are 5.6, 15.3, and 20 ton-years, respectively. The results represent a benchmark where the radioactivity background is reduced to a negligible level. Only the irreducible solar neutrino background is considered in the analysis. There is no new physics presenting in the pseudo data. For each m s , we first fix m φ , m Z = 10 keV to obtain a ∆χ 2 as function of the coupling constants. From this one-dimensional ∆χ 2 , we obtain the 95% C.L. limit with ∆χ 2 = 3.84 on couplings as shown in Fig. 10.
The sensitivity varies among PandaX-4T, XENONnT, and Lux-Zeplin due to different exposure time of 2 years, 5 years, and 1000 days, respectively. All experiments can probe the best-fit values in Fig. 7 and Fig. 8. They can also improve the sensitivity by roughly one order of magnitude. For example, the best-fit value of pseudo-scalar interaction is 5.6 × 10 −12 at XENON1T and the constraint becomes O(10 −13 ) in Fig. 10. With this improvement, the future DM experiments can confirm or falsify the pseudo-scalar explanation to the XENON1T excess.

VI. CONCLUSIONS
The XENON1T electron recoil excess stimulated various explanations with new physics including solar neutrino scattering with light mediators. Although the atomic effects are commonly considered for those dark matter scenarios, they are omitted in the solar neutrino explanations. In the first part of this work, we establish a systematic second quantization formalism for both the initial-state bound and final-state ionized electrons. This approach introduces the atomic K-factor in a natural way with general interactions for both neutrino and dark matter scatterings. The K-factor calculation includes a summation over the final-state angular quantum number l to infinity. In practice, one needs to truncate the summation at some maximum value l max and this introduces some calculation error. A sub-percent precision in the atomic K-factor requires at least l max = 200 for the electron recoil energy of O(10) keV.
This formalism is then applied to the solar neutrino scattering with bound electrons into a sterile neutrino through scalar, pseudo-scalar, vector, and axial-vector interactions. The atomic effects decrease the cross section by 1 ∼ 5 times to modify the recoil energy spectrum. Besides, the electron momentum distribution inside atom smears the momentum transfer of scattering. Consequently, the differential cross sections become smoother and broader than the free case. We then use the bound electron scattering scenario to fit the XENON1T excess. The scalar, vector, and axialvector interactions are preferred over the background-only hypothesis by more than 3σ with best-fit coupling constants of O(10 −13 ). Meanwhile, the pseudo-scalar case is favored over the backgroundonly hypothesis by less than 2σ with a coupling constant of O(10 −12 ). In all cases, the best-fit value of sterile neutrino mass is around 100 keV. Our results demonstrate for the first time that the atomic effects can not be ignored when using solar neutrino scattering to explain the XENON1T electron recoil excess.
Finally, we project the signal event rate and sensitivity at PandaX-4T, XENONnT, and LZ. These future experiments improve the sensitivity by roughly one order. The current best-fit values and the 1σ C.L. regions of the scalar, vector, and axial-vector mediators are marginally below the constraints from astrophysics and pion decay. For the pseudo-scalar case, most regions of the parameter space are in tension with the current constraint. With a factor of 10 improvement in sensitivity, the solar active-sterile neutrino conversion with bound electrons via light mediators as an explanation to the XENON1T excess can be tested.