key: cord-0059772-myc906dl
authors: Qiao, Chen-Kai; Lin, Shin-Ted; Chi, Hsin-Chang; Jia, Hai-Tao
title: Relativistic impulse approximation in the atomic ionization process induced by millicharged particles
date: 2021-03-19
journal: J High Energy Phys
DOI: 10.1007/jhep03(2021)184
sha: a38362a04c9f4bdaa2f4954b06bae9e8c9a98df6
doc_id: 59772
cord_uid: myc906dl

The millicharged particle has become an attractive topic to probe physics beyond the Standard Model. In direct detection experiments, the parameter space of millicharged particles can be constrained from the atomic ionization process. In this work, we develop the relativistic impulse approximation (RIA) approach, which can duel with atomic many-body effects effectively, in the atomic ionization process induced by millicharged particles. The formulation of RIA in the atomic ionization induced by millicharged particles is derived, and the numerical calculations are obtained and compared with those from free electron approximation and equivalent photon approximation. Concretely, the atomic ionizations induced by mllicharged dark matter particles and millicharged neutrinos in high-purity germanium (HPGe) and liquid xenon (LXe) detectors are carefully studied in this work. The differential cross sections, reaction event rates in HPGe and LXe detectors, and detecting sensitivities on dark matter particle and neutrino millicharge in next-generation HPGe and LXe based experiments are estimated and calculated to give a comprehensive study. Our results suggested that the next-generation experiments would improve 2-3 orders of magnitude on dark matter particle millicharge δ(χ) than the current best experimental bounds in direct detection experiments. Furthermore, the next-generation experiments would also improve 2-3 times on neutrino millicharge δ(ν) than the current experimental bounds.

Charge quantization is one of the most profound and fascinating open problems in physics. After exploring for several decades, the validity of charge quantization is still not revealed. All the observed elementary particles in Standard Model have quantized electric charge, but JHEP03(2021)184 the underling nature is still a mystery and can't be solved in the context of Standard Model. The electric charge quantization can be predicted in many theories beyond the Standard Model, i.e., grand unifications [1] , magnetic monopoles [2, 3] , and extra dimensions [4, 5] . However, no clear evidences have ever been provided as a confirmation for such theories.

Recently, a lot of studies propose a category of new particles, whose electric charge is tiny and non-quantized. These particles are named as "millicharged particles" [6] 1 and has stimulated a number of theoretical and experimental investigations [7, 8] . The experimental studies on millicharged particles can be carried out through reactor experiments [9] , positronium decays [10] , solar and celestial bodies observations [8] , cosmic microwave background (CMB) [11] , big bang nucleosynthesis (BBN), supernovas [12] , 21-cm line observations [13] , accelerator and collider experiments [14] [15] [16] . These experiments strictly constrain the mass of millicharged particles, as well as its electric charge over the past few years.

The millicharged particle has tiny electric charge, it could have electromagnetic interactions with target atoms or molecules. For instance, it can cause the following atomic ionization process

Therefore, in the direct detection experiments, the millicharged particles can be caught and detected from the signals produced by atomic ionization processes. In future superterranean or underground experiments, i.e., the China Dark matter EXperiment (CDEX) experiments located at China Jin-Ping underground Laboratory (CJPL) [17] , through searching the ionization, scintillation or heat signals in real detectors produced by atomic ionization processes, we can detect these unknown millicharged particles and give a constrain on their parameter space.

Previous studies on atomic ionization process induced by millicharged particles are often proceeded based on the free electron approximation (FEA) [9] . In the FEA, the calculation is the simplest and atomic many-body effects are neglected. The FEA approach works well in the high-energy transfer region, where atomic binding effects are negligible and the atomic electrons are approximately free. However, in the low-energy transfer region, the atomic many-body effects (including atomic binding effects, electron shielding, electron correlation effects and other many-body interactions) become dominant, thus the FEA breaks down and new approaches inspired by many-body physics are needed. It is worth noting that, the low-energy transfer region plays a crucial role in the experimental search for millicharged particles through atomic ionization processes, because the differential cross sections and reaction event rates in detectors have a dominant enhancement in this region. Particularly, detailed studies shown that, the FEA results may underestimate the differential cross section by more than 1 order in the low-energy transfer region, compared with the results from many-body physics [18] . Therefore, in this case, the FEA approach may produce large errors and it must be corrected through considering atomic many-body effects.

Apart from the FEA, there are other approaches used in the studies of atomic ionizations induced by millicharged particles. For instance, the equivalent photon approximation (EPA) is frequently used [9] . In the EPA approach, the atomic many-body effects can be partly considered. In this formulation, contributions coming from virtual photons in the electromagnetic interaction are equivalent to contributions from real photons. However, there is a fatal weakness: the EPA is valid only when energy transfer is extremely small (in the T → 0 limit). Therefore, a precise method which could deal with atomic many-body effects in the entire energy region is especially needed.

Recently, a new approach -the relativistic random-phase approximation -is applied to the studies of millicharged particles to treat the atomic many-body effects [18] [19] [20] [21] . However, this method is complicated in numerical calculations, especially when the incident particle energy is extremely high [22] . A relatively simple approach that could contain atomic many-body effects in the entire energy region would be helpful.

Inspired by the previous researches in many-body physics, in this work, we develop the relativistic impulse approximation (RIA) in the atomic ionization process induced by millicharged particles. The original ideas and framework of RIA approach are developed in previous years to handle a number of electromagnetic interactions in atomic physics, such as atomic Compton scattering [23] [24] [25] , electron impact [26] and other atomic processes [27, 28] . The atomic many-body effects can be treated effectively in the RIA approach. With the advantages of simplicity and flexibility, the RIA formulation has been widely applied to atomic and molecular physics [29] [30] [31] , condensed matter physics [32, 33] , nuclear and elementary particle physics [34] [35] [36] . In particular, in the conventional Monte Carlo simulation program Geant4 [37] , which are extensively used in nuclear and elementary particle physics experiments, many processes are treated using the formulation of RIA [38] [39] [40] .

In the present work, we develop the RIA approach for the atomic ionization process induced by millicharged particles. The formulation of RIA is derived, and the numerical calculations are obtained and compared with those from FEA and EPA approaches. These comparisons should give us information of the influences brought by atomic many-body effects, especially in the low-energy transfer regions. Our RIA approach developed in this work is quite general, irrelevant to the material composition of detectors and the underling nature of millicharged particles.

In this work, we consider two categories of millicharged particles: millicharged dark matter particles and millicharged neutrinos. On the one hand, dark matter problem is one of the most important topics in elementary particle physics, astrophysics, astronomy and cosmology. Currently, accumulating evidences in special rotation curves of spiral galaxies [50] [51] [52] , gravitational lensing [53] , large scale structure formation [54, 55] , cosmic microwave background and baryon acoustic oscillations [56, 57] have indicated that there are large amount of non-luminous dark matter in our universe [41, 58] . Therefore, direct detection of dark matter particle becomes an extremely significant and urgent work [41] . Traditionally, the most promising candidates for these unknown dark matter particles are weakly interacting massive particles (WIMPs), which only interact with ordinary matter through weak interaction beyond the Standard Model [41, 58, 59] . Other candidates, such as axion, sterile neutrino, mirror dark matter and so on [60] [61] [62] [63] , are also actively studied in JHEP03(2021)184 the dark matter detections. The millicharged particles, due to its ultra-tiny electromagnetic interaction, can successfully generate the dark matter relic abundance and give consistent results with astrophysical observations, which makes it become a candidate of dark matter particles [64] [65] [66] . On the other hand, neutrino physics also becomes a promising field in elementary particle physics, astrophysics, astronomy and cosmology, for it can reveal many aspects of physics beyond the Standard Model, e.g., baryon non-conservation [67] , matter-antimatter asymmetry [68, 69] , neutrino oscillation [70, 71] , seesaw mechanism [72] , and their Dirac or Majorana nature of fermions [42, 73] . Interestingly, there is an overlap, recent studies suggested that neutrino may become millicharged particles and they may have tiny electromagnetic interactions [74, 75] .

In the numerical calculations, we choose Ge and Xe elements as detector materials to study the atomic ionization processes induced by millicharged particles. The Ge and Xe elements are ideal materials for experimental detection of charged or neutral particles. With sufficient low threshold, large effective volume, high efficiency and ultra low background, high-purity germanium (HPGe) and liquid xenon (LXe) detectors are most generally used in particle physics experiments, especially for dark matter direct detections and neutrinoless double beta decay experiments [17, [41] [42] [43] [44] [45] [46] [47] [48] [49] . Concretely, in the present work, we study the ionization of Ge and Xe atoms by millicharged particles, calculating the differential cross section and the reaction event rate of these ionization processes in real HPGe and LXe detectors. The low-energy transfer and near threshold regions, where atomic many-body effects could have great impacts, are especially considered. The estimation of detecting sensitivities for dark matter particle and neutrino millicharge in next-generation HPGe and LXe based experiments is also provided according to the reaction event rate the and experimental background level.

Furthermore, in the actual ab initio calculations, the influence coming from relativistic effects of atomic electrons is also a critical point in dealing with atomic ionization process induced by millicharged particles. For the deep inner-shell atomic electrons, their motion around atomic nucleus could be very rapid with relativistic corrections at the order of Zα em , 2 where α em ≈ 1/137 is the fine structure constant for electromagnetic interactions [76, 77] . Particularly, for Xe atomic, this value could reach Zα em ≈ 0.4, suggesting that electron relativistic effects would play a significant role the same as electron many-body effects. Some previous studies have shown that the relativistic effects have non-negligible contributions to the scattering of atomic electrons with photon, ions, neutrinos and dark matter particles [18-20, 29-31, 78] . Recently, B. M. Roberts et al. concluded that the electron relativistic effect could give a large enhancement on cross sections as well as reaction event rates in the ionization process of atomic electrons induced by weakly interacting massive particles (WIMPs) [78] [79] [80] . These studies strongly imply that in the atomic ionization process induced by millicharged particles, which we are interested, the relativistic effects of atomic electrons may also become a non-negligible issue.

In atomic and molecular physics, the relativistic effect and atomic many-body effects

can be treated efficiently in the Dirac-Fock theory, which is the relativistic extension of the self-consistent Hartree-Fock method [81] [82] [83] [84] [85] . In this theory, the ground state wavefunctions are obtained by solving the fully relativistic many-body Dirac-Fock equation for atomic systems. The Dirac-Fock theory, since it was developed in the 1970s, has been widely applied to a number of atomic and molecular processes and has been confirmed by spectroscopic observations and scattering experiments in the past few decades. Therefore, in order to incorporate relativistic effects and many-body effects in the calculation of ground state wavefunctions and electron momentum distributions for atomic systems, we adopt the fully relativistic Dirac-Fock theory in this work. This paper is organised as follows: section 2 gives an introduction of the millicharged particle; section 3 briefly describes the general ideals for RIA approach; section 4 is devoted to theoretical derivation of RIA approach in the atomic ionization process induced by millicharged particles; numerical results and discussions are given in section 5 and section 6 for millicharged dark matter particle as well as millicharged neutrino; and conclusions and future perspectives are summarized in section 7. Furthermore, in the appendices, we give descriptions on free electron approximation (FEA), equivalent photon approximation (EPA) and the Dirac-Fock theory.

This section gives a brief introduction to the millicharged particles. The mechanism giving rise to the millichaged particle and the current experimental bounds for millicharged particles are mainly discussed.

The millicharged particle can be obtained from theories beyond the Standard Model [12, 14] . In particular, we will describe two mechanisms that could give rise to millicharged particles in this section.

First, millicharged particles can be generated in the extension of Standard Model by introducing an additional unbroken local U (1) h gauge group to the Standard Model gauge group [12, [86] [87] [88] . All Standard Model particles are singlets under the new gauge group U (1) h . We also add a massive hidden fermion χ charged under the new gauge group U (1) h only. Therefore, together with the Abelian gauge group U (1) Y in the Standard Model, there are two Abelian gauge groups: U (1) Y and U (1) h . The two gauge fields associated with gauge groups U (1) Y and U (1) h can couple to each other through the kinetic mixing. The Lagrangian for this model is: 3

Here, γ µ is the conventional Dirac-γ matrices, B µ is the gauge field of U (1) Y group, C µ is the gauge field of U (1) h group, f is Standard Model fermion and χ is the hidden fermion

are the field strength for B µ and C µ , respectively. Moreover, J B µ = g Sf γ µ f is the current associating with the U (1) Y gauge field B µ , and J C µ = g hχ γ µ χ is the current associating with the additional U (1) h gauge field C µ . The κ denotes the kinetic mixing parameter between two gauge fields B µ and C µ .

To make the physical picture clearer, we introduce the following two gauge fields A µ andà µ as the combination of gauge field B µ and C µ :

After the definition and re-coupling of A µ andà µ , the Lagrange density can be rewritten as [88] :

From this rearrangement, it can be clearly manifested that, the bosonic field A µ is coupled with currents J B µ and J C µ , while the bosonic fieldà µ is coupled with J C µ only. In this picture, our universe can be divided into two parts: the Standard Model sector and the "hidden sector". Accordingly, J B µ = g Sf γ µ f is the current in the Standard Model sector, while J C µ = g hχ γ µ χ can be viewed as "current" in the hidden sector, with A µ andà µ to be the ordinary photon and "dark photon" in the Standard Model sector and hidden sector. 4 The g h is the coupling between the "dark photon"à µ and the hidden sector fermion χ, and g S is the coupling between photon A µ and Standard Model fermion f . There is one important point should be noted, based on eq. (2.3), the "current" J C µ = g hχ γ µ χ in the hidden sector not only couples with the dark photonà µ , but also couples with photon A µ . Therefore, in this picture, a fermion χ living in the hidden sector not only acts as a charged particle in the hidden sector, but also behaves likes a charged particle in the Standard Model sector. Its electric charge can be determined through the coupling between "current" J C µ and photon A µ :

Assuming the kinetic mixing is extremely small, namely κ 1, then the electric charge of hidden sector fermion χ is tiny, which makes it to be a millicharged particle. In this case, the electric charge of the millicharged particle χ can be further simplified as q χ = δ χ e ≈ −κg h .

Figure 1. The current experimental bounds for millicharged particles. In this figure, the horizontal axis labels the mass of millicharged particles m χ , and the vertical axis labels the millicharge δ χ . Both the horizontal axis and vertical axis are plotted employing logarithmic coordinates. The upper and lower panels correspond to two different mechanisms: upper panel shows the exclusion region for models with an additional U (1) h gauge group; lower panel shows the exclusion region for models with right-handed massive fermions, which are SU (2) L singlets in the Standard Model gauge group. This figure present measurements and observations from astrophysics, cosmology, and particle physics experiments. The cosmological and astrophysical bounds from the sun (SUN) [8] , horizontal branch stars (HB) [12] , red-giant (RG) [12, 92] , white dwarf (WD) [12, 92] , supernova (SN1987A) [93] , cosmic microwave background (CMB) [11, 12, 94] and big bang nucleosynthesis (BBN) [12] are denoted as dashed lines. Direct detections from the underground experiments (XENON [95] and Super-Kamiokande [96] ), reactor experiments (TEXONO) [9, 97] , positronium decays (OPOS) [10] are displayed as solid lines. The experimental constrains from accelerators and colliders (COLL [7, 92] , SLAC [98] , LHC [14] , E613 [99] , MiniBooNE [100] , ArgoNeuT [101] ) are also presented in this figure as comparisons. Furthermore, our estimations of detecting sensitivity for millicharged dark matter particles in next-generation LXe based experiments calculated using our RIA approach are also shown in this figure.

In the aforementioned models, the "dark photon"à µ could acquire mass through the Higgs mechanism or Stückelberg mechanism [87, 88] . Therefore, the parameter space of "dark photon" naturally consists of its mass mà and kinetic mixing κ. In recent years, the studies of "dark photon" have attracted considerable attention [89, 90] . Furthermore, the millicharged particle can also be generated by other mechanisms. For instance, a class of models can be constructed by introducing right-handed massive fermion, which is a singlet under the SU (2) L Standard Model gauge group [8] . The Lagrangian for these models is:

where γ µ is the conventional Dirac-γ matrices, and B µ is the gauge boson in the U (1) Y Standard Model gauge group. In these models, the right-handed massive fermion χ is the millicharged particle with mass to be m χ , and its electric charge is related to the millicharge δ χ via q χ = δ χ e. Particularly, the neutrino millicharge, which will be discussed in section 6, can be obtained in this way by introducing right-handed Dirac neutrinos [75] . The parameter space of the millicharged particle is defined by (m χ , δ χ ). Many experimental investigations have strongly constrained the parameter space of millicharged particles. Figure 1 gives the current experimental bounds for millicharged particles. This figure presents measurements and observations from astrophysics, cosmology, and particle physics experiments. The cosmological and astrophysical bounds from the sun (SUN) [8] , horizontal branch stars (HB) [12] , red-giant (RG) [12, 92] , white dwarf (WD) [12, 92] , supernova (SN1987A) [93] , cosmic microwave background (CMB) [11, 12, 94] and big bang nucleosynthesis (BBN) [12] are denoted as dashed lines. Direct detections from the underground experiments (XENON [95] and Super-Kamiokande [96] ), reactor experiments (TEXONO) [9, 97] , positronium decays (OPOS) [10] are displayed as solid lines. The experimental constrains from accelerators and colliders (COLL [7, 92] , SLAC [98] , LHC [14] , E613 [99] , MiniBooNE [100] , ArgoNeuT [101] ) are also presented in this figure as comparisons. Furthermore, this figure also gives our estimations of detecting sensitivity for millicharged dark matter particles in next-generation LXe based experiments calculated using our RIA approach developed in this work.

In this section, we give an introduction of the RIA approach used in electromagnetic interactions in atomic physics. The general ideas, physical pictures, and theoretical formulation of the RIA approach are introduced in details. The development of the RIA approach in the atomic ionization process induced by millicharged particles is given in section 4.

In the formulation of RIA, due to atomic binding effects, the atomic bound electrons in an atom have a momentum distribution, which can be determined through its ground state wavefunctions. In the scattering process, electrons with different momentum scattered with incident particle independently, the interference term between different momentum JHEP03(2021)184 electrons is omitted for simplicity. With the advantages of simplicity and flexibility, the RIA formulation has been extensively used in many atomic physics processes, especially in the atomic Compton scattering [23] [24] [25] [102] [103] [104] [105] , electron impact [26] , and other atomic processes [27] [28] [29] .

In the following part, we will use the atomic Compton scattering

as an example to illustrate the general pictures and basic ideas for RIA formulation. In atomic Compton scattering, consider an incident photon with energy ω i and momentum k i scattering with an atomic bound electron with energy E i and momentum p i . After scattering, the energy and momentum of emitted photon are ω f and k f , and energy and momentum of final state electron are E f and p f , respectively. Then the doubly-differential cross section (DDCS) of Compton scattering in RIA formulation is given by [24, 25] :

where r 0 is the electron classical charge radius, E i = p 2 i c 2 + m 2 c 4 and E f = p 2 f c 2 + m 2 c 4 are the energies of initial and final state electrons, respectively. The functions K i and K f are defined as:

The function X(K i , K f ) is proportional to the reaction probability of the free electron Compton scattering γ + e → γ + e, which is the scattering between the incident photon γ and electron momentum eigenstate |p i . It is defined as:

Here, ρ(p i ) denotes the momentum distribution of atomic electrons, which is calculated through ground state wavefunctions. From eq. (3.2), it is easy to see that in the RIA formulation electrons with different momentum eigenstate |p i scattered with photon γ independently, and the interference terms between different momentum eigenstates are omitted. In this approach, atomic many-body effects are mainly reflected in the momentum distribution of atomic bound electrons.

In the previous studies, Roland Ribberfors et al. pointed out that the reaction probability function X(K i , K f ) in eq. (3.2) is a slow-varying function with respect to the integration variable p i . Therefore, it can be pulled out of the integration [24, 25, 34] . Successively, Ribberfors et al. made an approximation for function X(K i , K f ): Coordinate system XY Z and xyz. Coordinate system XY Z is chosen such that the Z axis is along the direction of initial photon γ, and X axis can be chosen as arbitrary direction perpendicular to the Z axis. While the coordinate system xyz is chosen such that the z axis represents the momentum transfer direction. After Compton scattering, the momentum of the scattered photon is denoted as k f , and the momentum transfer vector q is defined as q ≡ k f − k i . with K i (p z ) and K f (p z ) defined as:

In the above expressions, q is the modulus of the momentum transfer vector q ≡ k f − k i , and p z is the projection of the electron's initial momentum on the momentum transfer direction

with energy E(p z ) defined by

In the above calculations, the coordinate system xyz is chosen such that the z axis represents the momentum transfer direction in the Compton scattering process. The coordinate systems xyz and XY Z are defined and illustrated in figure 2 . In many literatures [25, 34, 104] , a convenient approximation for p z component is proposed as follows:

This approximation works well for small p z values, however, it can cause notable discrepancies for large p z values.

Using the above assumptions, the DDCS of atomic Compton scattering process in the RIA formulation is given by:

In this expression, the correction factor J(p z ) in the DDCS is called as the atomic Compton profile [106] J(p z ) ≡ ρ(p i )dp x dp y (3.11) with ρ(p i ) to be the ground state electron momentum density of the atomic system. For most of the atomic systems, the momentum distribution is spherical symmetric, then the atomic Compton profile reduces to

In this work, we only consider the spherical symmetric cases, and we use a fully relativistic Dirac-Fock theory to calculate the ground states of atomic systems and obtain their atomic Compton profiles. The DDCS of atomic Compton scattering in the RIA formulation in eq. (3.10) can be further simplified. In previous studies, an alternative and simpler approximation of the reaction probability function X(K i , K f ) was made by taking the p z → 0 limit of X(p z ), which finally gives its FEA value (also called as the Klein-Nishina value) [25, 104] 

Therefore, the simplified results of DDCS for atomic Compton scattering in RIA formulation can be expressed as:

From eq. (3.10) and eq. (3.14), it is obvious that the DDCS of atomic Compton scattering in the RIA approach factorizes into two parts

The factor Y RIA is dependent on the kinematical and dynamical properties of atomic Compton scattering, and it is irrelevant to the electronic structure of target materials. The correction factor J(p z ), known as the Compton profile, is related to the momentum distributions of electrons in the atomic or molecular ground state. In the RIA approach, all the atomic many-body effects can be incorporated into atomic Compton profiles.

Given the DDCS in atomic Compton scattering, the differential cross section with respect to the energy transfer can be calculated through the integration

with T = ω i − ω f to be the energy transfer in Compton scattering. With the atomic many-body effects incorporated into atomic Compton profiles, the RIA formulation could overcome the shortcomings in the FEA formulation. Therefore, it is a practical approach to calculate the Compton scattering in the low-energy transfer region or near photoionization threshold region. With the advantages of simplicity and flexibility, the RIA formulation has been widely applied to atomic [25] , condensed matter [32, 33] , nuclear and elementary particle physics [34, 36] . In particular, in the Monte Carlo simulation program Geant4 [37] , which are extensively used in nuclear and elementary particle physics, several algorithms employ the RIA approach to treat the Compton scattering process [38] [39] [40] . Furthermore, atomic Compton profile J(p z ) can also reflect some important information in condensed matter physics and material science, i.e. the electronic structure [107, 108] , electron momentum distribution [108, 109] , electron correlation [33, 110] , band structure, and Fermi surface [111, 112] .

In this section, we develop the RIA approach to the atomic ionization process induced by millicharged particles

The derivation of the RIA formulation for the atomic ionization process is presented in detail. The differential cross section of the atomic ionization process is focused and discussed. The general results of the doubly-differential cross section (DDCS) are given in subsection 4.1, and simplified results of DDCS are given in subsection 4.2. In subsection 4.3, we give comments on our newly developed RIA approach for atomic ionization process induced by millicharged particles. Finally, the explicit expressions for differential cross section with respect to energy transfer are presented in 4.4.

For the atomic ionizations induced by millicharged particles χ + A → χ + A + + e − , consider the millicharged particle with energy E χ and momentum k χ . After the ionization, the energy and momentum of millicharged particle become E χ and k χ , respectively. Similar to the cases discussed in section 2, in RIA formulation, the electron in an atom has a momentum distribution ρ(p i ), which is calculated through the electron momentum wavefunction of atomic ground state. Further, atomic electrons with different momentum p i JHEP03(2021)184 scatter independently with millicharged particle χ, and interactions between electrons with different momentum are omitted. After the ionization process, final state electron gets momentum p f . Therefore, in the RIA approach, the DDCS of the atomic ionization process χ + A → χ + A + + e − is calculated by summing over contributions from all possible momentum p i :

where Ω χ is the solid angle for scattered millicharged particles,

are energies of initial and final state electrons, respectively. The function X is proportional to the reaction probability of the scattering between millicharged particle and electron momentum eigenstate, namely the scattering process χ + e − → χ+ e − with electron momentum eigenstate |p i . In the atomic Compton scattering, the function X is given by eq. (3.4) in section 2. In the atomic ionization process induced by millicharged particles, the probability function X should be calculated through the scattering amplitude of χ + e − → χ + e − . This process is very similar to the Rutherford scattering process p + + e − → p + + e − . In analogy with the Rutherford scattering, the function X can be written as [113] :

The function X in eq. (4.2) is directly obtained from the probability of Rutherford scattering with the replacement: p → χ, m p → m χ and e → q χ = δ χ e. Here, s, t, u are Mandelstam variables defined as:

According to the property of Mandelstam variables s + t + u = 2m 2 e c 4 + 2m 2 χ c 4 , the variable u can be simplified as:

Choosing an appropriate coordinate system will benefit the numerical calculation. In this work, the coordinate system xyz is chosen similar to the cases in atomic Compton scattering (in figure 2 ). In the xyz system displayed in figure 3 , the z axis represents the momentum transfer direction. After introducing such coordinate system, the momentum component p z is determined by energy and momentum conservations p µ i + k µ χ = p µ f + (k χ ) µ . The explicit expression for p z is given by: figure 2 ). The coordinate system XY Z is chosen such that the Z axis is along the direction of initial momentum k χ for millicharged particle, and X axis can be chosen as arbitrary direction perpendicular to the Z axis. The momentum for the scattered millicharged particle is denoted as k χ , and the momentum transfer vector q χ in atomic ionization process is defined as q χ ≡ k χ − k χ . In this coordinate, the axis z represents the momentum transfer direction.

with q χ to be the modulus of the momentum transfer q χ = k χ −k χ in the scattering process

and E(p z ) = m 2 e c 4 + p 2 z c 2 . Assuming that millicharged particles are massive particles, then the initial and final state momentum k χ and k χ can be calculated as:

Similar to the cases in Compton scattering, from the energy and momentum conservations, it can be revealed that p z and E(p z ) are exactly the minimal energy and momentum of the initial state electrons activated in the ionization process, namely

Furthermore, the momentum component p z can be approximated as:

However, it should be noted that eq. (4.9) only valid when p z is sufficiently small.

Similar to the cases in Compton scattering, the probability function X in the integrand of eq. (4.1) is can be averaged and pulled out of the integration as Roland Ribberfors et al. did in Compton scattering in reference [24, 25] (see eq. (3.5) in section 2). Concretely, we can take the following approximation:

And the corresponding values of Mandelstam variables can be expressed by:

Obviously, the above approximation of probability function X in eq. (4.10) made from eq. (4.2) indicates that the electron initial momentum p i is specified only in the momentum transfer direction z, while momentum components in other directions p x and p y are omitted for simplicity. Based on the above assumptions, we substitute the approximation (4.10) into eq. (4.1) and simplify the energy of electron as

Finally, the DDCS of atomic ionization process induced by millicharged particles can be expressed as:

In this subsection, we provide a simpler version in the calculation of DDCS in atomic ionization process induced by millicharged particles. Simpler results of DDCS can be achieved by making more simplified approximation for probability function X in the calculation of eq. (4.1). For instance, similar to the cases in atomic Compton scattering, an alternative and simpler approximation of function X in eq. (4.10) can be made by taking the p z → 0 limit of X(s(p z ), t(p z ), u(p z )), which finally gives:

Correspondingly, the 3 Mandelstam variables s, t, u can be further simplified as:

Using the approximation (4.13), the DDCS of atomic ionization process induced by millicharged particles can be further simplified as:

From the above results in eqs. (4.12) and (4.15), it is evident that the DDSC of the atomic ionization process induced by millicharged particles can be summarized as:

This results of DDCS for atomic ionization process induced by millicharged particle in eq. (4.16) is similar to the cases in atomic Compton scattering introduced in section 3 (the eq. (3.15)). From eq. (4.12), it is clearly that the DDCS of the atomic ionization process also factorizes into two parts: the factor Y RIA and the atomic Compton profile J(p z ). The kinematical and dynamical property of the atomic ionization process is incorporated in factor Y RIA , irrespective of the material elements and electron structures in atomic systems. The correction from the atomic effects and electronic structures is mainly incorporated into the atomic Compton profile J(p z ) as in atomic Compton scattering.

It should be noted that the above results in eqs. (4.1) -(4.15) only correspond to the single electron systems. However, the detector atom is usually a multi-electron system and consists of electrons form different subshells. Summing over contributions from all subshell electrons, the DDCS of the multi-electron atomic system can be calculated as:

In the above expressions, J njl (p z ) is the atomic Compton profile for subshell (njl)

The E B njl is the atomic binding energy of subshell (njl), Z njl is the number of electron in

When the energy transfer T = E χ − E χ is less than the subshell binding energy E B njl , electron in subshell (njl) is inactive in the atomic ionization process induced by millicharged particles χ + A → χ + A + + e − , and this subshell gives zero contributions in the DDCS.

Similarly, given the DDCS of the atomic ionization process induced by millicharged particles, the differential cross section with respect to the energy transfer T in this process can be obtained through the integration

with T = E χ − E χ to be the energy transfer for atomic ionization process induced by millicharged particles. Put the DDCS in eq. (4.17) and eq. (4.18) into the integration, we finally get the explicit expressions for differential cross section with respect to energy transfer

To summarise, this section gives the theoretical derivation of our RIA approach in the atomic ionization process induced by millicharged particles. A promising feature is that our approach is quite general, depend neither on the underling nature or mechanism of millicharged particles, nor on the composition of detector materials. Therefore, it can be extensively applied to the studies of millicharged particles. In this work, we also develop a numerical program based on the above approach. The numerical calculations are presented in the next two sections for millicharged dark matter particles and millicharged neutrinos.

This section is devoted to the numerical results of the atomic ionization process induced by millicharged dark matter particles. Based on our RIA approach derived in section 4, a numerical program is developed utilizing the basic Fortran language. In subsection 5.1, the differential cross sections with respect to energy transfer are obtained for Ge and Xe atom, and results from our RIA approach are compared with those from FEA and EPA approaches. In subsection 5.2, the differential reaction event rates in HPGe and LXe detectors are given for typical experimental environments. Furthermore, in subsection 5.3, we give an estimation of the detecting sensitivities on dark matter particle millicharge δ χ in next-generation HPGe and LXe based experiments. These numerical results presented in this section can shed light on theoretical investigations as well as experimental explorations for millicharged dark matter particles.

In this subsection, we provide numerical calculations on differential cross section with respect to energy transfer for the atomic ionization process induced by millicharged dark matter particles. Figure 4 shows the differential cross section dσ/dT for atomic ionization process induced by high-energy millicharged dark matter particles. In this figure, the mass and initial energy of millicharged particle are chosen as m χ c 2 = 1 KeV and E χ = 1 MeV, respectively. The millicharge of dark matter particle is chosen to be δ χ = 10 −12 as a typical example. The numerical results come from FEA, EPA and RIA approaches are given in this figure for comparisons. Among these approaches, the FEA results calculated through eq. (A.3) neglect all the atomic many-body effects, and they could provide a approximate result only in the high-energy transfer region, in which the atomic electron is nearly free and atomic binding effects become very weak. However, FEA results fail to give a precise prediction in the low-energy transfer region because atomic many-body effects have a strong effect on the atomic ionization process. The simplified FEA results, which are calculated through eq. (A.4), can be reduced from FEA results when m χ m e , T m e c 2 , T E χ are satisfied. The EPA results calculated using eq. (B.7) could provide a more precise results than the FEA results in the ultra-low-energy transfer region by including atomic manybody effects partly. The EPA approach can be derived from quantum field theory when energy and momentum transfer are extremely small, namely in the T → 0 limit, and it becomes invalid in high-energy region. The introduction of FEA and EPA approaches is give in the appendices. Our RIA approaches developed in this work could deal with atomic many-body effects in the entire region, regardless of the underlining nature of millicharged particles and the composition of detector materials.

From figure 4 , it can be clearly manifested that the differential cross sections dσ/dT from FEA, EPA and RIA calculations all diminish as energy transfer T increases. In the low-energy transfer region, both EPA and RIA results acquire larger cross sections than FEA results, indicating that the atomic many-body effects, including the atomic binding, JHEP03(2021)184 . Differential cross sections of atomic ionization process for Ge and Xe atoms induced by millicharged dark matter particles. The mass and initial energy of millicharged particle are chosen as m χ c 2 = 1 keV and E χ = 1 MeV, and the dark matter particle millicharge is chosen to be δ χ = 10 −12 . In this figure, we compare the numerical results on differential cross section dσ/dT calculated in the FEA, EPA and RIA approaches. The red solid lines correspond to the FEA results calculated through eq. (A.3); red dashed lines represent the simplified FEA results calculated through eq. (A.4); blue lines stand for the EPA results calculated from eq. (B.7); black solid lines show the RIA results calculated using eq. (4.22); and black dashed lines present the simplified RIA results calculated using eq. (4.23).

electron shielding and electron correlation, could greatly enhance the atomic ionization process induced by millicharged particles and enlarge their differential cross sections. Particularly, the EPA results for Xe atom present large peak when energy transfer T ∼ 100 eV. This is because in the EPA approach, the differential cross section dσ/dT for atomic ionization process induced by millicharged particles is proportional to the photonabsorption cross section, as shown in appendix B (see eq. (B.7)). For photoabsorption cross section, there is a giant resonance for 4d electrons of Xe atom in the 100 eV region [114] [115] [116] [117] [118] . 6 For our RIA results, in the low-energy transfer region, our RIA results get larger cross section than those from FEA results; while in the high-energy transfer region, the RIA results do not exhibit notable differences with respective to FEA results. The physical reason can be explained naturally: when energy transfer T is much larger than the atomic binding energy E B 1s for 1s electron (which is 11.1 keV for Ge atom and 34.5 keV for Xe atom), the atomic effects can be neglected and the atomic electron is approximately free. However, when energy transfer T is sufficient low and is comparable to the atomic binding energy E B 1s , atomic binding, electron shielding as well as electron correlation effects become dominant. In these cases, atomic electrons can no longer be treated as free electrons, which lead to large deviations between RIA and FEA results in the low-energy transfer region. Furthermore, figure 4 also indicates that our RIA results are approaching to the EPA results when energy transfer T is extremely small, especially in the T → 0 limit. It can be viewed as a demonstration for the validity and availability of our RIA approach developed in the present work. Figure 4 also manifested that, for large incident particle energy E χ = 1 MeV, the simplified RIA results calculated using eq. (4.23) converge to the RIA results calculated using eq. (4.22) in the entire region of energy transfer T . Therefore, for high-energy millicharged dark matter particles, among the approximations of probability function X in the integrand in eq. (4.1), the more simplified approximation X ≈ X sim is good enough, and it does not lead to notable deviations compared with the more accurate approximation X ≈ X(s(p z ), t(p z ), u(p z )).

In figure 5 , we given the differential cross sections for low-energy millicharged dark matter particles. The initial energy of millicharged dark matter particle is chosen as E χ = 10 keV, while the dark matter particle millicharge is chosen to be δ χ = 10 −12 the same as in figure 4 . The upper and lower panels correspond to m χ c 2 = 1 keV and m χ c 2 = 10 eV, respectively. In these cases, the incoming particle energy E χ is not large enough to make the atomic many-body effects negligible. Therefore, our RIA results does not converge to the FEA results in the entire region of T , but our RIA results still come close to the EPA results in the ultra-low-energy T → 0 limit. From this figure, we can also observe that the EPA results acquire larger cross sections for smaller millicharged dark matter particle mass m χ . There is another notable point should be noted: large differences between the simplified FEA results and the FEA results emerge in energy range T > 0.1 keV. In the appendix A, it would be clarified that the simplified FEA results 6 For Xe atom, there is also a peak in T ∼ 700 eV region due to the resonance for 3d electrons of Xe atom in the photoabsorption cross section [119] . The peak of 3d electrons for Xe atom in T ∼ 700 eV is relatively smaller than that of 4d electrons in T ∼ 100 eV JHEP03(2021)184 For low-energy millicharged dark matter particles, e.g. for incident particle energy E χ = 10 keV in figure 5 , the condition T E χ is not satisfied when T > 0.1 keV. Therefore, only FEA results calculated through eq. (A.3) are reasonable in such cases. Compare with the FEA results, the simplified FEA results calculated through eq. (A.4) overestimate the differential cross sections of the atomic ionization process for Ge and Xe atoms. Furthermore, for small JHEP03(2021)184 incident particle energy E χ = 10 keV, there are discrepancies between simplified RIA results and RIA results when energy transfer T > 1 keV. For low-energy millicharged dark matter particles, the more simplified approximation X ≈ X sim of probability function X could bring about some deviations, and it is better to use the more accurate approximation X ≈ X(s(p z ), t(p z ), u(p z )) to evaluate the differential cross section dσ/dT .

Particularly, results in figure 4 and figure 5 could reflect some sort of generality. In various approaches, i.e. FEA, EPA and RIA, the differential cross section dσ/dT for the atomic ionization process induced by millicharged dark matter particles is proportional to δ 2 χ . For the same incoming energy, results correspond to other millicharge δ χ can be obtained from figure 4 and figure 5 by proportional magnifying or shrinking the results by δ 2 χ times.

In this subsection, we shall give the numerical calculations of differential reaction event rate for atomic ionization process induced by millicharged particles in typical super-terranean or underground experiments.

In a typical experimental environment, such as CDEX experiment located in CJPL as well as other super-terranean or underground experiments, the differential reaction event rate in detectors for atomic ionization process induced by millicharged particles can be expressed as:

where ρ A is the number density of detector atoms, φ χ is the total flux of millicharged dark matter particles, and dφ χ /dE χ is the flux spectrum at a given incoming energy E χ . In eq. (5.1), E min χ and E max χ are the maximal and minimal energy of the millicharged dark matter particles that could enter into the detectors. From the eq. (5.1), it can be clearly manifested that the energy spectrum of the ionization process dσ/dT and the flux spectrum of the millicharged particles dφ χ /dE χ totally determine the differential reaction event rate dR/dT in a typical experiment environment.

In this work, for simplicity, we assume that incoming millicharged dark matter particles all come from the cosmic rays. 7 Although the behaviour of dark matter particles in the cosmic rays is still an open question, several recent studies suggested that the millicharged dark matter particles could be accelerated analogous to Standard Model charged particles in cosmic rays [9, 96] through the Fermi acceleration mechanism [120] [121] [122] . Therefore, as a result, the flux spectrum of millicharged dark matter particle dφ χ /dE χ obeys a simple power low [9, 96] :

There are other sources could give rise to millicharged dark matter particle flux, e.g., using theoretical calculation and experimental measurements, reference [9] also consider the millicharged particles come from nuclear reactors as well as earth atmosphere.

where α is the power index with a fixed value α = 2.7 [96] , and sr represents the steradian.

A key point should be mentioned is that: eq. (5.2) is satisfied under certain conditions, demanding that the millicharged dark matter particle should be ultra relativistic [96] . In this work, the minimal incoming energy of millicharged dark matter particle is chosen to be E min χ = 10 m χ c 2 in the numerical calculations, the same as in reference [9] .

In typical experimental environments, the differential reaction event rates for atomic ionization process induced by millicharged particles are given in the figure 6 for HPGe and LXe detectors. The numerical results obtained from FEA, EPA, and RIA approaches are displayed for comparisons. The FEA, EPA and RIA results on differential event rates dR/dT are obtained by integrating the differential cross sections dσ/dT through eq. (5.1). The differential cross sections dσ/dT are calculated using FEA, EPA and RIA approaches as in subsection 5.1. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. 8 The mass of millicharged dark matter particle is set as m χ c 2 = 1 keV, and its millicharge is chosen to be δ χ = 8 × 10 −9 for HPGe detector and δ χ = 4 × 10 −9 for LXe detector, respectively.

From figure 6 , it is indicated that the reaction event rates in HPGe and LXe detectors decrease rapidly as energy transfer T becomes higher. Therefore, the low-energy transfer region is dominant in the atomic ionization process induced by millicharged dark matter particles, and this region should be pay close attention to in the direct detection experiments for millicharged dark matter particles. For this reason, figure 6 only presents the event rates in T < 10 keV region, and the T > 10 keV region is omitted. Particularly, figure 6 shows that, for RIA results with millicharged dark matter particle mass m χ c 2 = 1 keV, the differential event rate in HPGe detector at energy range T ∼ 0.1 keV is about dR/dT ∼ 0.1 cpkkd, and the differential event rate in LXe detector at energy range T ∼ 0.5 keV is roughly dR/dT ∼ 10 −4 cpkkd.

The same as figure 4 in subsection 5.1, results in figure 6 can also reveal the influences come from atomic many-body effects that act on the atomic ionization process induced by millicharged particles. In low-energy transfer region, the EPA and RIA approaches obtain more reaction event rates in HPGe and LXe detectors, compared with those from FEA results. This exhibits a similar tendency with the differential cross sections presented in figure 4 in subsection 5.1. In the low-energy transfer region, the atomic binding, electron shielding as well as electron correlation effects can greatly enhance the atomic ionization processes induced by millicharged particles. This would bring desirable news to the next-generation direct detection experiments for millicharged dark matter particles. Furthermore, in the T → 0 limit, in which range the EPA approach is derived, our RIA results do not appear large deviations from the EPA results, both for HPGe and LXe detectors. It can also be provided as an indication for the validity of our RIA approach developed in this work in the low-energy transfer region. The differential reaction event rates for atomic ionization process induced by millicharged particles. The mass of millicharged dark matter particle is set as m χ c 2 = 1 keV, and its millicharge is chosen to be δ χ = 8 × 10 −9 for HPGe detector and δ χ = 4 × 10 −9 for LXe detector, respectively. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. The red solid lines correspond to the FEA results; red dashed lines represent the simplified FEA results; blue lines stand for the EPA results; black squares display the RIA results; gray triangles show the simplified RIA results.

when energy transfer T > 0.2 keV. This is caused by the differential cross section in FEA approach. The simplified FEA results on differential cross section calculated using eq. (A.4) converge to the FEA results calculated using eq. (A.3) only when m χ m e , T m e c 2 , T E χ are satisfied. However, for the millicharged dark matter particle coming from cosmic rays, the flux spectrum represent a power low as in eq. (5.2). When energy of incident particle is lower, the flux spectrum becomes larger. In typical superterranean or underground experiments, there are large amount of low-energy millicharged dark matter particles entering into the HPGe and LXe detectors, which would destroy the condition T E χ and make the simplified FEA results inappropriate. 9 When energy transfer T becomes higher, the condition T E χ is harder to satisfy, and there are more differences between the simplified FEA results and FEA results. This is similar to the cases of differential cross sections discussed in subsection 5. 

According to the calculations of reaction event rates in HPGe and LXe detectors, we can give an estimation of the detecting sensitivity for dark matter particle millicharge δ χ in next-generation HPGe and LXe based direct detection experiments. The estimation is carried out according to the following assumptions:

• Consider the dark matter particle with millicharge δ χ = δ 0 , if the calculated reaction event rates in energy range above the experimental threshold surpass the experimental background, then the signals from atomic ionization process induced by millicharged particles can be catched and identified effectively. In this case, next-generation experiments have the ability to detect dark matter particles with millicharge δ χ = δ 0 . 10

• On the other hand, if the calculated reaction event rates in energy range above the experimental threshold are less than the experimental background, then the atomic ionization signals would not be effectively identified and next-generation experiments could't set a constrain on dark matter particles with millicharge δ χ = δ 0 .

In the numerical calculations, for any dark matter particle mass m χ , we calculate the differential event rates dR/dT in HPGe and LXe detectors using FEA, EPA and RIA methods at a given dark matter particle millicharge δ χ , then adjust the value of millicharge 9 The minimal energy of incident millicharged dark matter particle is chosen as E min χ = 10 mχc 2 in the numerical calculations. In figure 6 , the minimal energy corresponds to E min χ = 10 mχc 2 = 10 keV, and the condition T Eχ is not satisfied with great accuracy when T > 0.2 keV. 10 Only when reaction event rate overwhelm the experimental background in energy range above the experimental threshold, signals produced from the atomic ionization process induced by millicharged particles could be catched and identified effectively. Otherwise, the energy transfer is too small to track the atomic ionization signals in detectors, or these atomic ionization signals may be overwhelmed by background signals and couldn't be identified and analyzed effectively.

HPGe Based Experiments m χ c 2 detecting sensitivity on dark matter particle millicharge δ χ keV FEA Results EPA Results RIA Results 0.01 3.2 × 10 −9 3.5 × 10 −10 1 × 10 −9 0.1 6.0 × 10 −9 1.9 × 10 −9 2 × 10 −9 1

1 δ χ such that reaction event rates in HPGe and LXe detectors in energy region above the experimental threshold surpass the experimental backgrounds. Finally, we can obtain the lower limit of milliarge δ χ satisfying the above conditions. This is the estimation of detecting sensitivity on dark matter particle millicharge δ χ in the next-generation HPGe and LXe based direct detection experiments.

The estimation of detecting sensitivity on dark matter particle millicharge δ χ in the next-generation HPGe and LXe based experiments is shown in table 1 for several dark matter particle mass m χ . The results from the FEA, EPA and RIA calculations are given in this table for comparisons. For the HPGe based next-generation experiments, the energy threshold and background level have been assumed as 100 eV and 0.1 cpkkd, respectively. For the LXe based next-generation experiments, the energy threshold and background level are assumed as 500 eV and 10 −4 cpkkd, respectively. From this table, it can be clearly shown that, for several dark matter particle mass, the detecting sensitivities of millicharge δ χ calculated from RIA and EPA approaches are much larger than those calculated from the FEA approach. In subsection 5.1 and subsection 5.2, we have learned that the atomic any-body effects can greatly enhance the atomic ionization process induced by millicharged dark matter particles in the low-energy transfer region, leading to the increase of differential cross sections dσ/dT as well as differential reaction event rates dR/dT in this region. Therefore, for the same experimental background, atomic many-body effects make it more JHEP03(2021)184 easy to let reaction event rates surpass the experimental background, which eventually leads to a more strong constrain on dark matter particle millicharge δ χ . This would be beneficial for direct detection of millicharged dark matter particles in next-generation experiments. These results shown that atomic many-body effects would play a significant role in the electromagnetic interactions of millicharged particles, and it may open an new window for the explorations of millicharge particles. Furthermore, with relatively lower experimental background, the next-generation LXe based experiments could set a lower bound on dark matter particle millicharge δ χ , no matter which approach is employed in the numerical calculations. For HPGe based experiments, the EPA results get smaller dark matter millicharge δ χ than RIA results in the low-mass cases (m χ c 2 ≤ 0.1 keV). While for LXe based experiments, the EPA results get smaller millicharge δ χ than RIA results in all cases (10 eV < m χ c 2 < 100 keV) because of the giant resonance for d electrons in the photoabsorption cross section of Xe atom [114] [115] [116] [117] [118] [119] . Detailed numerical results giving rise to the detecting sensitivities in figure 1 are presented in the appendix D.

The estimations of detecting sensitivity on dark matter particle millicharge δ χ in this work can contribute to the parameter space of millicharged dark matter particles. In the figure 1, we also present our estimations of detecting sensitivity on millicharge δ χ in RIA calculations for next-generation LXe based experiments. From figure 1 , the indirect searches from astronomy and cosmology set stronger constrains on dark matter particle millicharge δ χ . While the direct detection experiments and accelerator/collider experiments, i.e. XENON10, TEX, OPOS, COLL, SLAC, LHC in figure 1, set looser bounds for millicharge δ χ . However, it is remarkable that the next-generation LXe based direct detection experiments would greatly increase the detecting sensitivity on dark matter particle millicharge δ χ . In the range 10 eV < m χ c 2 < 100 keV, The current best experimental bound in direct detection experiments and accelerator/collider experiments is roughly δ χ ∼ 10 −5 to δ χ ∼ 10 −7 , which is 2-3 order of magnitude larger than our estimation for next-generation LXe based experiments.

There is one point need to be mentioned: the calculations of reaction event rates in HPGe and LXe detectors as well as the calculations of detection sensitivities on dark matter particle millicharge δ χ in next-generation HPGe and LXe based experiments are just a leading order estimation. In our numerical calculations, we have made some simplified assumptions. The electromagnetic interactions between millicharged dark matter particles in cosmic rays and the charged particles in earth atmosphere, as well as the electromagnetic interactions between millicharged dark matter particles and atoms and molecules in environmental rocks, are not taken into considerations. These interactions may lead to an upper bound in the parameter space of millicharge dark matter particles, as inidicated in reference [9] . If the millicharge of dark matter particle is much too large, then the electromagnetic interactions between millicharged dark matters in cosmic rays and the charged particles in earth atmosphere would be too strong, which lead to tremendous attenuation of dark matter particle flux in the atmosphere. As a result, it will prevent millicharged dark matter particles entering into HPGe and LXe detectors in super-terranean or underground experiments.

In section 1, it is revealed that neutrino physics is becoming a rising field in many branches of science. Recently, many studies suggested that neutrinos may have tiny electromagnetic interactions [74, 75, 123] , and they may have millicharge as well as magnetic moment. Theoretical and experimental explorations on neutrino millicharge and magnetic moment is becoming more and more attractive, and a number of researches on this area emerge in recent years [18, 124, 125] .

As discussed in section 4, the RIA approach we developed in this work is irrelevant to the underling nature and mechanism of millicharged particles. In principle, our approach can also be applied to the study of millicharged neutrinos. In this section, we use our RIA approach to study atomic ionization process induced by millicharged neutrinos. The numerical results on differential cross section dσ/dT , differential reaction event rate dR/dT , and detecting sensitivity on neutrino millicharge δ ν in next-generation direct detection experiments are presented similar to the cases of millicharged dark matter particles discussed in section 5.

There are several kinds of sources which may contribute to millicharged neutrinos: reactor neutrinos, cosmological neutrinos, solar neutrinos, atmospheric neutrinos, supernova neutrinos, cosmogenic neutrinos, and active galactic nucleus (AGN) produced neutrinos [121, 126, 128] . The reactor neutrinos become dominant only when laboratory is near the nuclear reactors, and supernova neutrinos become notable when supernova is activated, i.e. supernova 1987A burst. For other sources, the cosmological neutrinos mainly appear in ultra-low energy range (below 1 eV), while atmospheric neutrinos, cosmogenic neutrinos and AGN produced neutrinos all centered in ultra-high energy range (above GeV). More details of neutrino sources and their flux can be found in reference [126, 127] . Therefore, in energy range 100 eV ≤ E ν ≤ GeV, which is sensitive to HPGe and LXe detectors and is of great interests in direct detection experiments, solar neutrino have the maximal flux and can be viewed as the main source of millicharged neutrinos. In this work, for simplicity, we only consider solar neutrinos as the source of millicharged neutrinos. Contributions from other sources are leaving for future studies.

There are several channels which can produce solar neutrinos [128, 129] : displayed in figure 7 . Furthermore, it should be mentioned that the pp neutrinos and 7 Be neutrinos give predominant contributions to solar neutrino flux, and this two channels contribute to 98% of solar neutrinos [21, 121, 128] .

In this subsection, we take 7 Be solar neutrinos as an example to study the differential cross section with respect to energy transfer in the atomic ionization process induced by millicharged neutrinos. The differential cross section dσ/dT of atomic ionization process for Ge and Xe atoms is presented in figure 8 . The 7 Be solar neutrinos have discrete spectrum, and the incoming neutrino energy is located in E ν = 384 keV and E ν = 862 keV with branch ratios 89.5% and 10.5% [121] . Form figure 8, it can be clearly shown that the results for millicharged neutrinos are very similar to the results for millicharged dark matter particles given in subsection 5.1. The differential cross section of atomic ionization process induced by millicharged neutrinos diminish as energy transfer T increases, and FEA, EPA, RIA results show the same tendency. When the energy transfer T is smaller than the atomic binding energy for 1s electron (which is 11.1 keV for Ge atom and 34.5 keV for Xe atom), the differential cross sections calculated using RIA and EPA methods are larger than those in FEA results, indicating that the atomic many-body effects tend to intensify the electromagnetic interaction for millicharged neutrinos in low-energy transfer region. In the ultra-low-energy transfer JHEP03(2021)184 region, namely the T → 0 limit, our RIA results are near the EPA results, which shows the validity of our methods in this region. When the energy transfer T is sufficiently large, the EPA approach breaks down and it underestimate the differential cross sections, while the RIA results converge to the FEA results because the atomic many-body effects are negligible in high-energy transfer region. Furthermore, from figure 8, it is also indicated that the simplified RIA results converge to the RIA results in the entire region of energy transfer T . Therefore, for 7 Be solar neutrinos with energies E ν = 384 keV and E ν = 862 keV, among the approximations of function X in the integrand of eq. (4.1), the more simplified approximation X ≈ X sim is adequate and does not lead to large deviations compared with the more accurate approximation X ≈ X(s(p z ), t(p z ), u(p z )). This is similar to the cases of high-energy millicharged dark matter particles presented in figure 4 .

Similarly, the results in figure 8 can also reflect some sort of generality. In various approaches, i.e. FEA, EPA and RIA, the differential cross section dσ/dT of the atomic ionization process induced by millicharged neutrinos is proportional to δ 2 ν . For the same incoming energy E ν , results correspond to other neutrino millicharge δ ν can be obtained by proportional magnifying or shrinking the results in figure 8 by δ 2 ν times. From the numerical calculations in subsection 5.1 and subsection 6.1, we can draw a conclusion that the atomic ionization processes, whether induced by millicharged dark matter particles or millicharged neutrinos, exhibit similar tendency. The differential cross section dσ/dT of atomic ionization process induced by millicarged particles drops rapidly as energy transfer becomes higher. In the low-energy transfer region, atomic binding, electron shielding and electron correlation effects could greatly enhance the atomic ionization process induced by millicharged particles. Our RIA approach developed in this work is appropriate in the entire region of energy transfer. For high energy millicharged particles, our RIA results on differential cross section show small discrepancies with EPA results in the ultra-low energy region, and our RIA results successfully converge to the FEA results in the high-energy transfer region, where the atomic effects are weak and atomic electron can be treated as free electron approximately. 11

Similar to the calculations in subsection 5.2 for millicharged dark matter particles, in a typical experiment environment, the differential reaction event rate in HPGe and LXe detectors for atomic ionization process induced by millicharged neutrinos can be expressed similar to eq. (5.1):

where dφ ν /dE ν is the neutrino flux spectrum. As we have discussed in the beginning of this section, in the energy range relevant to direct detection experiments, which is from keV to GeV, solar neutrino is the main source for millicharged neutrinos. Since the pp channel and 7 Be channel contribute to 98% of solar neutrinos, we can omit contributions from other channels. Therefore, the solar neutrino flux spectrum can be simplified as:

For 7 Be reaction neutrino, the flux spectrum dφ Be ν /dE ν is discrete with energy located at 384 keV and 862 keV. For pp reaction neutrino, the flex spectrum dφ pp ν /dE ν is continuous. The flux spectra dφ pp ν /dE ν and dφ Be ν /dE ν and can be obtained either by fitting the corresponding curves in figure 7 , or from the solar neutrino databases [131] [132] [133] [134] [135] [136] . Figure 9 shows the differential reaction event rates dR/dT for atomic ionization process induced by millicharged neutrinos for HPGe and LXe detectors in typical super-terranean or underground experimental environments. The neutrino mass is chosen to be m ν c 2 = 0.1 eV, and the neutrino millicharge is set as δ ν = 10 −12 for HPGe detectors and δ ν = 2.5 × 10 −13 for LXe detectors. 12 The numerical results from FEA, EPA and RIA approaches are given in this figure for comparison. Similar to the cases for millicharged dark matter particles, figure 9 also indicates that the differential event rates for atomic ionization process induced by millicharged neutrinos reduce significantly as energy transfer T increases, both in HPGe and LXe detectors. Therefore, to search the millicharged neutrino in direct detection experiments, we should focus on the low-energy transfer region. Furthermore, figure 9 shows that, in the low-energy transfer region, the differential event rates dR/dT calculated using RIA and EPA approaches are larger than those from FEA results, indicating the atomic many-body effects could greatly enhance the atomic ionization process induced by millicharged neutrino in the low-energy transfer region. This totally agree with the conclusions for millicharged dark matter particles discussed in subsection 5.2. For our RIA results, they converge to the FEA results as energy transfer T increases. Meanwhile, our RIA results slowly approach to the EPA results when energy transfer T becomes very small, but the convergence between RIA results and EPA results in the ultra-low-energy region (in T → 0 limit) is not as good as those of millicharged dark matter particles presented in the figure 6 in subsection 5.2, as well as in figures 10-14 in appendix D.

Similarly, according to the calculated reaction event rates in HPGe and LXe detectors, we can give an estimation of the detecting sensitivity on neutrino millicharge δ ν in next-generation HPGe and LXe based experiments. In the next-generation direct detection experiments, we assume that the energy threshold and background level for HPGe based experiments would be 100 eV and 0.1 cpkkd. Meanwhile, the energy threshold and background level for LXe based experiments would reach 500 eV and 10 −4 cpkkd. From figure 9 , it is clearly manifested that, for neutrino mass m ν c 2 = 0.1 eV, the reaction event rates in HPGe and LXe detectors in energy range above the experimental thresholds successfully 12 The neutrino millicharge δν is adjusted such that the reaction event rates in HPGe and LXe detectors are comparable to the experimental background levels in next-generation HPGe and LXe based experiments. Figure 9 . The differential reaction event rates for atomic ionization process induced by millicharged neutrinos for HPGe and LXe detectors. The neutrino mass is chosen to be m ν c 2 = 0.1 eV, and the neutrino millicharge is set as δ ν = 10 −12 for HPGe detectors and δ ν = 2.5×10 −13 for LXe detectors. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. The same as in figure 6 , the red solid lines correspond to the FEA result; red dashed lines represent the simplified FEA results; blue lines stand for the EPA results; black squares display the RIA results; gray triangles show the simplified RIA results.

suppress the experimental background levels. Therefore, the next-generation direct detection experiments have the ability to push the detecting sensitivity of neutrino millicharge to δ ν ∼ 10 −12 for HPGe based experiments and δ ν ∼ 2.5 × 10 −13 for LXe based experiments.

There is one important point should be noted. Since the flux spectrum of solar neutrino is irrelevant to the neutrino mass m ν , both the reaction event rates in HPGe and LXe detectors and the estimated detecting sensitivities on neutrino millicharge δ ν in nextgeneration experiments do not have an obvious dependency on neutrino mass m ν . 13 This is different with the cases of millicharged dark matter particles. Firstly, the dark matter particle flux spectrum in eq. (5.2) is manifestly mass and energy dependent. Secondly, if the mass of millicharge dark matter particle is smaller, then the minimal energy E min χ of imcoming dark matter particles, which has been set as E min χ = 10m χ c 2 to ensure the ultrarelativistic property of millicharge dark matter particles, becomes lower and give rise to large numbers of low-energy millicharged dark matter particles according to the power law in the flux spectrum in eq. (5.2). These points make the reaction event rates in HPGe and LXe detectors and the estimated detecting sensitivity on dark matter particle millicharge δ χ in next-generation HPGe and LXe based experiments highly depend on dark matter particle mass m χ .

Based on the numerical calculations presented in this subsection, in table 2, we give the comparison between our estimated detecting sensitivities in next-generation experiments and current experimental bounds on neutrino millicharge δ ν in the direct detection experiments. From this table, the estimated detecting sensitivity on neutrino millicharge δ ν in next-generation experiments is roughly 2-3 times smaller than the current best experiment bound. 14 Therefore, the next-generation HPGe and LXe based experiments have the potential to make a great progress on the detecting ability of millicharged neutrinos.

Similar to the cases of millicharged dark mater particles presented in subsection 5.3, the numerical calculations of reaction event rates in HPGe and LXe detectors as well as the calculations of detection sensitivities on neutrino millicharge δ ν in next-generation HPGe and LXe based experiments are just a leading order estimation. In our calculations, the interactions between millicharged neutrinos and charged particles in earth atmosphere as well as the interactions between millicharged neutrinos and atoms and molecules in environmental rocks are not taken into considerations. These interactions may lead to an upper bound in the parameter space of millicharged neutrinos, as discussed in subsection 5.3.

13 Furthermore, the differential cross section of atomic ionization process induced by millicharged neutrinos does not show obvious differences when neutrino mass mν varies. The mass of neutrino is much too small compared with its energy Eν in solar neutrino spectrum, and it is approximately massless in these cases.

14 Very recently, an excess of electron recoil events was reported in the XENON1T experiment [47] . Amir N. Khan interpreted these signals to be nonstandard neutrino interactions, and a constrain was given on neutrino millicharge δν based on the electron recoil excess in XENON1T experiment [141] . Khan's results indicated that neutrino millicharge would be δν = (1.7 − 2.3) × 10 −12 . However, other studies also suggested that the electron recoil excess maybe caused by some experimental backgrounds, which were ignored in the experimental analysis [142, 143] . Therefore, more experimental data are needed to confirm this excess.

HPGe Detector Group detecting sensitivity on neutrino millicharge δ ν Our Results 2.5 × 10 −12 (FEA Results) 7 × 10 −13 (EPA Results) 1 × 10 −12 (RIA Results) TEXONO [18] 2.1 × 10 −12 GEMMA [137] 1.5 × 10 −12 / 2.7 × 10 −12 (based on different methods) LXe Detector Group detecting sensitivity on neutrino millicharge δ ν Our Results 4 × 10 −13 (FEA Results) 9 × 10 −14 (EPA Results) 2.5 × 10 −13 (RIA Results) XENON1T [46] 6.4 × 10 −13 PandaX [138] 2.06 × 10 −12 Projected DARWIN [139] 2.4 × 10 −13 Projected LZ [140] 2.8 × 10 −13 Table 2 . Detecting sensitivity on neutrino millicharge δ ν in the next-generation direct detection experiments. The results for HPGe and LXe based experiments in the FEA, EPA and RIA calculations are given in this table. In the HPGe based next-generation experiments, the energy threshold and background level would be 100 eV and 0.1 cpkkd, respectively. The energy threshold and background level for LXe based next-generation experiments could reach 500 eV and 10 −4 cpkkd. The current experimental bounds on neutrino millicharge δ ν obtained from direct detection experiments are also given for comparisons. In this table, for the XENON1T, PandaX, projected DARWIN and projected LX experimental bounds, we present the results in reference [21] analyzed using the experimental data.

In this work, we develop the RIA approach in the atomic ionization process induced by millicharged particles. Our approach is inspired and benefited from many-body physics, because the RIA approach was originally invented in atomic and molecular physics to study the various electromagnetic interactions in atoms and molecules. In the experimental detections for millicharged particles, atomic many-body effects play a crucial role and they cannot be arbitrarily neglected in theoretical calculations, especially when energy transfer in the scattering process is not very high (comparable to atomic binding energies). Our new developed RIA approach could effectively handle with atomic many-body effects in the entire energy region, therefore it can give a more precise results than traditional FEA approach in the study of millicharged particles. With this superiority, our new developed RIA approach would be helpful and have impacts on the theoretical predictions and experimental investigations for millicharged particles. In the present work, the formulation of RIA is derived for atomic ionization process induced by millicharged particles, and a numerical program is developed based on our RIA approach. The numerical results obtained using our RIA approach are compared with those from FEA and EPA approaches.

Concretely, we study the atomic ionization process induced by millicharged dark matter particles as well as millicharged neutrinos in HPGe and LXe detectors. The differential cross section with respect to energy transfer, the differential reaction event rate in HPGe and LXe detectors, and the estimated detecting sensitivity for next-generation HPGe and LXe based experiments are presented in this work. In addition, to incorporate relativistic effects of atomic electrons from the beginning, the fully relativistic Dirac-Fock theory is used to obtain the ground state wavefunctions, electron momentum distributions and atomic Compton profiles for atomic systems. For the differential cross section dσ/dT , when energy transfer is smaller than the binding energies E B 1s for 1s electron (which is 11.1 keV for Ge atom and 34.5 keV for Xe atom), our RIA results present large discrepancies with respect to the FEA results. In this region, the atomic many-body effects play a significant role in the atomic ionization process induced by millicharged particles, which lead to the breaking down of FEA approach. On the other hand, in high-energy transfer region, where electrons can be treated as free electrons approximately, there are no notable differences between our RIA results and FEA results. In the ultra-low-energy transfer region where EPA approach is derived (namely in T → 0 limit), our RIA results does not exhibit large differences with respect to EPA results. The above conclusions show that our RIA approach developed in this work is valid and could treat the atomic many-body effects efficiently in the entire region of energy transfer T . Another important point is that both RIA and EPA results give larger cross sections than FEA results in the low-energy transfer region, indicating that atomic manybody effects could greatly enhance the atomic ionization process induced by millicharged particles. Furthermore, for the reaction event rates dR/dT in HPGe and LXe detectors, numerical calculations present the similar phenomenon and tendency with those obtained from differential cross sections.

According to the calculated reaction event rates in HPGe and LXe detectors, we give an estimation of the detecting sensitivity on dark matter particle and neutrino millicharge δ χ and δ ν in next-generation HPGe and LXe based experiments. The energy threshold and background level for next-generation HPGe based experiments are postulated as 100 eV and 0.1 cpkkd, and the energy threshold and background level for next-generation LXe based experiments are postulated as 500 eV and 10 −4 cpkkd, respectively. Our numerical results show that, with relatively lower backgrounds, the next-generation LXe based experiments probably could give a better constrain on neutrino and dark matter particle millicharge. For millicharged dark matter particles, the next-generation LXe based experiments would give a tremendous improvement on dark matter particle millicharge δ χ in mass range 0.01 ≤ m χ c 2 ≤ 100 keV for direct detection experiments, with estimated detecting sensitivity 2-3 orders of magnitude smaller than the current best experimental bound in direct detection experiments. For millicharged neutrino, the estimated detection sensitivity of neutrino millicharge δ ν for next-generation LXe based experiments would reach δ ν ∼ 2.5 × 10 −13 , which is improved by roughly 2-3 times than the current best experimental bound.

In particular, our RIA approach developed in the present work is quite general, neither depends on the underling nature and mechanisms for millicharged particles, nor on the atomic and molecular composition of detector materials. In this work, we choose HPGe

and LXe detectors to given a comprehensive study on the atomic ionization process induced by millicharged particles. Other detector materials, such as liquid argon (LAr), sodium iodide (NaI) and cesium iodide (CsI), still deserve to study in the following works. Furthermore, the physical ideas and formulation for RIA approach can also be applied to other electromagnetic interactions for millicharged particles, e.g. the atomic Compton scattering between millicharged particles and dark photons, the magnetic moment interactions for millicharged neutrinos, as well as other processes relevant to the millicharged particles and detector atoms or molecules. We wish our work could enlarge the understanding of millicharged particles and push forward studies in the related fields.

In this appendix, we give a brief description of the free electron approximation (FEA) in the atomic ionization process induced by millicharged particles. In the FEA formulation, atomic electrons are treated as free electrons, and atomic bindings, electron shielding, electron correlation as well as other many-body effects are neglected.

In the atomic ionization process induced by millicharged particles χ+A → χ+A + +e − , assuming the electric charge of millicharged particle is q χ = δ χ e, the differential cross section for this process in the FEA formulation can be expressed as [9] :

where m e and m χ are the mass of electron and millicharged particle, E χ is the incoming energy of millicharged particle, and T = E χ − E χ is the energy transfer in the atomic ionization process. When the millicharged particle mass m χ and energy transfer T are both sufficiently small, namely the conditions m χ m e , T m e c 2 , T E χ are satisfied, the above differential cross section in eq. (A.1) for ultra-relativistic millicharged particles (with m χ c 2 E χ ) can be simplified as [18] :

The above results in eqs. (A.1) and (A.2) only correspond to the single-electron system. To calculate the differential cross section for multi-electron atomic systems, contributions from all subshell electrons should be summed over to give the following results

with E B njl and Z njl to be the atomic binding energy and number of electron in (njl) subshell. Similarly, when m χ m e , T m e c 2 , T E χ are satisfied, the differential cross section

can be simplified as:

The FEA formulation works well when the energy transfer T is much larger than the atomic binding energy, in which cases the atomic many-body effects can be neglected and electrons are approximately free. Previous studied have confirmed that when energy transfer T is comparable to the atomic binding energy, FEA results underestimate the differential cross section in the atomic ionization process induced by millicharged particles [9, 19, 124] , which implies atomic effects could intensify the atomic ionization process χ + A → χ + A + + e − . This conclusion is consistent with our numerical results presented in section 5 and section 6.

Equivalent Photon Approximation (EPA) is an approaches widely used in nuclear and elementary particle physics, especially in Quantum Electrodynamics (QED) and Quantum Chromodynamics (QCD) calculations [144, 145] .

In the EPA formulation, considering the process e − + X → e − + Y . The electron scatters with particle X by exchanging virtual photons. In the low-momentum transfer limit q → 0, the contribution from the longitudinal polarized virtual photons vanishes, and the contribution coming from virtual photons can be equivalent to those from real photons. Therefore, the cross section of the whole process e − + X → e − + Y can be calculated from the cross section of its subprocess γ + X → Y [113, 144] :

where α em is the conventional fine-structure constant, s = (p µ e + p µ X ) 2 is the total energy square in the center-of-mass frame, z = q/p e is the ratio between and virtual photon momentum q (the momentum transfer in the whole process e − + X → e − + Y ) and electron momentum p e , and σ(γ + X → Y ) is the cross section for subprocess γ + X → Y . In the eq. (B.1), f γ (z) is the Weizsacker-Williams distribution function defined as:

It can be viewed as the probability of finding a photon with momentum q = p z from the incident electron beam [144] . Similarly, in the atomic ionization process induced by millicharged particles χ + A → χ + A + + e − , the EPA approach connects its cross section with the cross section of photoionization process γ + A → A + + e − . In analogy with eq. (B.1), the cross section of JHEP03(2021)184 atomic ionization process χ + A → χ + A + + e − can be expressed as:

Here, α χ is the "fine structure constant" in the electromagnetic interactions induced by millicharged particles. It is defined as:

where q χ = δ χ e is the electric charge of the millicharged particle. To extract differential cross section in the EPA formulation, we consider the case that both energy transfer T and momentum transfer q is sufficiently small, namely in the limit

In this case, the total cross section for photoionization process can be simplified as [9, 19, 124] :

with σ γ abs (T ) to be the total cross section for photoabsorption process at incident photon energy T , and R 0 T (q 2 = 0) to be the atomic transverse response function for on-shell real photons in zero-momentum transfer cases. In this work, we only duel with the case that the mass of millicharged particle m χ is tiny and much smaller than the mass of detector atom m A . In this case, laboratory frame can be viewed as the center-of-mass frame approximately, and the total energy square in the center-of-mass frame can be simplified as s = (p µ e + p µ X ) 2 ≈ E 2 χ . 15 In eq. (B.3), a divergent part arise in the z → 0 limit. This divergent part can be cancelled by the higher order corrections with the help of renormalization. Finally, after tidies calculations, the differential cross section of the atomic ionization process induced by millicharged particles in the EPA approach can be expressed as [19] :

From the discussions above, it can be clearly shown that the derivation of eq. (B.7) requires that the energy transfer T and the momentum transfer q to be sufficiently small (in the q → 0 and T → 0 limits). Therefore, the EPA formulation only works well in the ultra-low-energy transfer region (T → 0). When energy transfer T is large, EPA results would bring about large discrepancies, and this point has been confirmed by recent researches [19] . This conclusion is also consistent with our numerical results presented in section 5 and section 6.

In this appendix, we give an introduction of the relativistic Dirac-Fock theory. We will focus on the construction of Dirac-Fock Hamiltonian for atomic systems, and how to calculate ground state wavefunctions, electron momentum distributions and atomic Compton profiles using the Dirac-Fock theory.

The Dirac-Fock theory [81] [82] [83] 146] , which is a relativistic extension of the nonrelativistic Hartree-Fock self-consistent method, is commonly used in ab initio calculations in atomic and molecular physics. In the last few decades, it has been confirmed by a number of experiments and has become a milestone in atomic and molecular physics [84, 85, 147, 148] . In this work, the Dirac-Fock theory is used to obtain the ground state wavefunctions, electron momentum distributions and Compton profiles for atomic systems.

In the Dirac-Fock theory, the total Hamiltonian for atomic systems is given by [83, 84, 148] :

Here, h (a) is the single-particle Hamiltonian for the a-th election, which includes the kinetic energy for a-th electron and the Coulomb potential between the atomic nuclei and this electron. The h (ab) is the two-particle Hamiltonian, which is the interaction between the a-th and b-th electrons. In the relativistic case, the single-particle Hamiltonian is given by: where α (a) is the conventional Dirac-α matrices for a-th electron, β is the Dirac-β matrix, the symbol ∇ (a) represent the gradient operator for a-th electron, and r (a) is the radius between atomic nuclei and the a-th electron. In this work, we only consider the leading order Coulomb interactions between electrons. Therefore, the two-particle Hamiltonian can be simply expressed as: 16

where r (ab) is the distance between the a-th and b-th electron.

With the the single-particle Hamiltonian h (a) and two-particle Hamiltonian h (ab) given in eqs. can be expressed as [83, 84, 148] :

In the following part, we give a description on how to calculate the ground state energies and wavefunctions for atomic systems in the Dirac-Fock theory. To satisfy the Pauli exclusion principle, in Dirac-Fock theory, the total ground state wavefunctions for atomic systems can be constructed through the Slater determinant of single-electron wavefunctions Ψ(r (1) , r (2) , · · · , r (Z) ) = 1

In the expression, u (a) (r (a) ) is the single-electron wavefunction for a-th electron, and r (a) is the position of a-th electron (with the center of atomic nucleus set as the coordinate origin).

In this work, we only consider the spherical symmetric atomic systems. Therefore, the single-electron wavefunction for atomic ground state with definite quantum number (nκm) = (njlm), which is also called as the Dirac orbital, has the following form [83, 84, 151] :

where G nκ (r) and F nκ (r) are the large and small components respectively, Ω κm (θ, φ) is normalized spherical spinor defined as:

where Y lm (θ, φ) is the spherical harmonics, lm − µ; 1 2 µ|jm is the Clebsch-Gordan coefficient, and χ µ is a spinor with s = 1/2 and s z = µ.

In many cases, only the radial part of Dirac orbital need to be focused, and the angular part can be separated and neglected for simplicity. Therefore, we can introduce the following two-component radial Dirac orbital:

After the introduction of Dirac orbital as well as its large and small components F nκ = F njl , G nκ = G njl , the Dirac-Fock equations for atomic systems can be expressed and solved

routinely. The total energies and ground state wavefunctions for atomic systems as well as the energy eigenvalues and single-electron wavefunctions for individual electrons can be obtained.

In the Dirac-Fock theory, the total energy for atomic system is calculated by solving the eigen-equation (1) , r (2) , · · · , r (Z) ) = E Dirac-Fock atom Ψ(r (1) , r (2) , · · · , r (Z) ) (C.9)

Put the Dirac-Fock Hamiltonian in eq. (C.4), atomic total wavefunction in eq. (C.5) and single-electron wavefunction in eq. (C.6) into eq. (C.9). After separation of variables for the angular part θ and φ, the total energy for atomic system can be expressed as [81, 82] :

(C.10) where p is the abbreviation for subshell (n p κ p ) = (n p j p l p ), Z p = 2j p + 1 is the number of electrons in subshell p = (n p κ p ) = (n p j p l p ). The f k (pq) and g k (pq) are angular coefficients, which are calculated through Wigner-3j coefficients, and k 0 , k 1 , k 2 are defined in reference [84] . In the eq. (C.10), integral I(pq) gives rise to the one-body interaction, and Slater integral R k (pqrs) represents the two-body interaction. The explicit expressions for I(pq) and R k (pqrs) can be found in reference [84] .

The total energy value for atomic system E Dirac-Fock atom relies on the large and small components of single-electron wavefunctions through the integrals I(p) and R k (pqrs). The ground state wavefunctions for atomic system in eq. (C.5) as well as the ground state wavefunction for individual electron in eq. (C.6) should minimize the total energy E Dirac-Fock atom in eq. (C.10). After the variational method, the differential equations for large and small components of single-electron wavefunctions become [81] [82] [83] :

where ε nκ is the energy eigenvalue for subshell (nκ) = (njl). The Y nκ (r) is the direct potential acts on large and small components of Dirac orbital for (nκ) subshell, while X • Take the solutions G (1) nκ , F (1) nκ , ε (1) nκ in the first step as new trial solutions, then calculate the direct and exchange potentials Y nκ (r), X 

nκ , ε (2) nκ , the same as in the first step.

• . . .

• Repeat the above procedures routinely. When the energy eigenvalue in the i step ε (i) nκ converges to the energy eigenvalue in the i+1 step ε (i+1) nκ , the correct energy eigenvalue ε nκ and the ground state wavefunction u nκ = (G nκ , F nκ ) for each subshell electron are solved from Dirac-Fock equation, and this iterative algorithm is self-consistent.

Once Dirac-Fock equations are solved, the large component G nκ = G njl , small components F nκ = F njl as well as the energy eigenvalue ε nκ = ε njl for different subshell electrons are obtained. Therefore, the corresponding electron momentum wavefunctions are given by the following Fourier transformation [152] :

where φ G njl , φ F njl are the large and small components of electron momentum wavefunctions of (njl) subshell, and j l (pr) is the spherical Bessel function. Based on electron momentum wavefunctions, the momentum distribution of electrons in atomic system is calculated as follows [83, 105] : ρ njl (p) = |φ njl (p)| 2 = (φ G njl (p)) 2 + (φ F njl (p)) 2 (C.13a)

Finally, the atomic Compton profile defined in eq. (3.11) can be calculated through the integration of electron momentum distributions ρ(p) and ρ njl (p).

J njl (p z ) = 2π ∞ |pz| pρ njl (p)dp (C.14a) J(p z ) = 2π ∞ |pz| pρ(p)dp = njl N njl J njl (p z ) (C.14b)

The atomic Compton profile, when plugged into eqs. (4.17)-(4.18) and eqs. (4.22)-(4.23), can give the differential cross section for the atomic ionization process induced by millicharged particles in RIA approach.

In this appendix, the detailed numerical results on reaction event rates in HPGe and LXe detectors calculated using FEA, EPA and RIA approaches are given for different dark matter particle mass m χ and their millicharge δ χ . The estimations of detecting sensitivity on dark matter millicharge δ χ for next-generation HPGe and LXe based experiments, which has been summarized in subsection 5.3, are obtained from these results. From the numerical results in figures 10-14, it is indicated that our RIA results are approaching to the FEA results when energy transfer T is large. For large dark matter particle mass m χ , i.e. m χ c 2 =100 keV, this tendency becomes more apparent. On the other hand, the difference between our RIA results and EPA results becomes tiny in the ultra-low-energy transfer region (namely in the T → 0 limit), especially for LXe detectors. This is consistent with the conclusions in subsection 5.2, indicating the validity of our RIA approach in the entire region of energy transfer T .

Furthermore, figures 10-14 also show that, for smaller dark matter particle mass m χ , the differences between FEA and EPA results become larger in the low-energy transfer region. When the dark matter particle mass m χ reduces, the discrepancies between the FEA results and simplified FEA results become more notable. In the appendix A, it is shown that the simplified RIA results reduced to the FEA results only when conditions m χ m e , T m e c 2 , T E χ are satisfied. For smaller dark matter particle mass m χ , the minimal energy for incoming millicharged dark matter particle, which is chosen as E min χ = 10 m χ c 2 in our numerical calculations, becomes lower. Therefore, there are large amount of millicharged dark matter particles entering into HPGe and LXe detectors, and the condition T E χ is harder to satisfy in this case.

JHEP03(2021)184 Figure 10 . The differential reaction event rates for the atomic ionization process induced by millicharged particles in HPGe and LXe detectors. The mass of millicharged dark matter particle is set as m χ c 2 = 10 eV, and its millicharge δ χ is chosen such that the reaction event rates in HPGe and LXe detectors match the experimental background levels for next-generation experiments, respectively. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. The same as in figure 6 and figure 9 , the red solid lines correspond to the FEA results; red dashed lines represent the simplified FEA results; blue lines stand for the EPA results; black squares display the RIA results; gray triangles show the simplified RIA results.

JHEP03(2021)184 Figure 11 . The differential reaction event rates for the atomic ionization process induced by millicharged particles in HPGe and LXe detectors. The mass of millicharged dark matter particle is set as m χ c 2 = 100 eV, and its millicharge δ χ is chosen such that the reaction event rates in HPGe and LXe detectors match the experimental background levels for next-generation experiments, respectively. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. The same as in figure 6 and figure 9 , the red solid lines correspond to the FEA results; red dashed lines represent the simplified FEA results; blue lines stand for the EPA results; black squares display the RIA results; gray triangles show the simplified RIA results.

Figure 12. The differential reaction event rates for the atomic ionization process induced by millicharged particles in HPGe and LXe detectors. The mass of millicharged dark matter particle is set as m χ c 2 = 1 keV, and its millicharge δ χ is chosen such that the reaction event rates in HPGe and LXe detectors match the experimental background levels for next-generation experiments, respectively. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. The same as in figure 6 and figure 9 , the red solid lines correspond to the FEA results; red dashed lines represent the simplified FEA results; blue lines stand for the EPA results; black squares display the RIA results; gray triangles show the simplified RIA results.

Figure 13. The differential reaction event rates for the atomic ionization process induced by millicharged particles in HPGe and LXe detectors. The mass of millicharged dark matter particle is set as m χ c 2 = 10 keV, and its millicharge δ χ is chosen such that the reaction event rates in HPGe and LXe detectors match the experimental background levels for next-generation experiments, respectively. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. The same as in figure 6 and figure 9 , the red solid lines correspond to the FEA results; red dashed lines represent the simplified FEA results; blue lines stand for the EPA results; black squares display the RIA results; gray triangles show the simplified RIA results. Figure 14 . The differential reaction event rates for the atomic ionization process induced by millicharged particles in HPGe and LXe detectors. The mass of millicharged dark matter particle is set as m χ c 2 = 100 keV, and its millicharge δ χ is chosen such that the reaction event rates in HPGe and LXe detectors match the experimental background levels for next-generation experiments, respectively. In this figure, the horizontal axis represents the energy transfer T , and the vertical axis represents the differential event rate dR/dT in unit of cpkkd. The same as in figure 6 and figure 9 , the red solid lines correspond to the FEA results; red dashed lines represent the simplified FEA results; blue lines stand for the EPA results; black squares display the RIA results; gray triangles show the simplified RIA results.

Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

Unity of All Elementary Particle Forces

Quantised singularities in the electromagnetic field

Quantum Field Theory of Dirac Monopoles and the Charge Quantization Condition

The Atomicity of Electricity as a Quantum Theory Law

Magnetic Monopoles in Kaluza-Klein Theories

Millicharged particles

Limits on particles of small electric charge

Minicharged Particles from the Sun: A Cutting-Edge Bound

Constraints on millicharged particles with low threshold germanium detectors at Kuo-Sheng Reactor Neutrino Laboratory

An Improved Limit on Invisible Decays of Positronium

Constraints on millicharged particles from Planck data

Dark Radiation constraints on minicharged particles in models with a hidden photon

Reviving Millicharged Dark Matter for 21-cm Cosmology

LHC probes the hidden sector

Probing dark matter particles at CEPC

Millicharged particles at electron colliders

Limits on Light Weakly Interacting Massive Particles from the First 102.8 kg × day Data of the CDEX-10 Experiment

Constraints on millicharged neutrinos via analysis of data from atomic ionizations with germanium detectors at sub-keV sensitivities

Constraining neutrino electromagnetic properties by germanium detectors

Low-energy electronic recoil in xenon detectors by solar neutrinos

Discovery potential of multiton xenon detectors in neutrino electromagnetic properties

Private conversation with L

Compton Scattering of X Rays from Bound Electrons

Relationship of the relativistic Compton cross section to the momentum distribution of bound electron states

Incoherent-x-ray-scattering functions and cross sections (dσ/dΩ) incoh by means of a pocket calculator

Double and triple differential cross sections for K-shell ionisation by relativistic electron impact

Resonant transfer and excitation in ion-atom collisions

Binary-encounter electrons observed at 0 o in collisions of 1-2-MeV/amu H + , C 6+ , N 7+ , O 8+ , and F 9+ ions with H 2 and He targets

Inelastic scattering of quasifree electrons on O 7+ projectiles

An overview of the theories used in Compton Scattering Calculations

Compton scattering revisited

Compton scattering and the study of electron momentum density distributions

Electron correlation effects on Compton profiles of copper in the GW approximation

Fast sampling algorithm for the simulation of photon Compton scattering

Fernández-Varea, Overview of physical interaction models for photon and electron transport used in Monte Carlo codes

Measurement of Low Energy Ionization Signals from Compton Scattering in a Charge-Coupled Device Dark Matter Detector

Recent developments in Geant4

GEANT4 Physics Reference Manual

A simple model of photon transport

A low energy bound atomic electron Compton scattering model for Geant4

Dark matter direct-detection experiments

Neutrinoless double beta decay and neutrino physics

Results from the Super Cryogenic Dark Matter Search Experiment at Soudan

Improved Limit on Neutrinoless Double-β Decay of 76 Ge from GERDA Phase II

Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment

Dark Matter Search Results from a One Ton-Year Exposure of XENON1T

Excess electronic recoil events in XENON1T

Search for Neutrinoless Double-β Decay with the Complete EXO-200 Dataset

Search for Majorana Neutrinos near the Inverted Mass Hierarchy Region with KamLAND-Zen

Die Rotverschiebung von extragalaktischen Nebeln

Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions

The Extended Rotation Curve and the Dark Matter Halo of M33

A direct empirical proof of the existence of dark matter

Formation of Galaxies and Large Scale Structure with Cold Dark Matter

The Evolution of Large Scale Structure in a Universe Dominated by Cold Dark Matter

Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation

Planck 2013 results. XVI. Cosmological parameters

Particle dark matter: Evidence, candidates and constraints

Dark Matter Candidates from Particle Physics and Methods of Detection

Axion Dark Matter, Proton Decay and Unification

Sterile neutrino Dark Matter

Mirror dark matter: Cosmology, galaxy structure and direct detection

Millicharged Atomic Dark Matter

Self-interacting asymmetric dark matter coupled to a light massive dark photon

Solving the small-scale structure puzzles with dissipative dark matter

The Neutrino Mass and Family, Lepton and Baryon Nonconservation in Gauge Theories

The Origin of the matter -antimatter asymmetry

Cosmic microwave background, matterantimatter asymmetry and neutrino masses

Evidence for oscillation of atmospheric neutrinos

Neutrino Mass and Mixing with Discrete Symmetry

Neutrino Mass and New Physics

Neutrinoless Double-Beta Decay: a Probe of Physics Beyond the Standard Model

Neutrino electromagnetic properties

Neutrino electromagnetic interactions: a window to new physics

Theoretical Atomic Physics

Handbook of Theoretical Atomic Physics

Ionization of atoms by slow heavy particles, including dark matter

Dark matter scattering on electrons: Accurate calculations of atomic excitations and implications for the DAMA signal

Electron-interacting dark matter: Implications from DAMA/LIBRA-phase2 and prospects for liquid xenon detectors and NaI detectors

Relativistic self-consistent fields

Relativistic atomic wave functions

A multiconfiguration relativistic Dirac-Fock program

GRASP: A general-purpose relativistic atomic structure program

Multiconfiguration Dirac-Hartree-Fock Calculations with Spectroscopic Accuracy: Applications to Astrophysics

Two U(1)'s and Epsilon Charge Shifts

A Stueckelberg extension of the standard model

The Stueckelberg Z-prime Extension with Kinetic Mixing and Milli-Charged Dark Matter From the Hidden Sector

Direct Detection Constraints on Dark Photon Dark Matter

Results and perspectives in dark photon physics

MeV Dark Matter Complementarity and the Dark Photon Portal

Updated bounds on millicharged particles

Astrophysical constraints on minicharged particles

Narrowing the window for millicharged particles by CMB anisotropy

Dark Matter Detectors as Dark Photon Helioscopes

Dark Cosmic Rays

Search for millicharged particles in reactor neutrino experiments: A Probe of the PVLAS anomaly

Search for millicharged particles at SLAC

Scattering of Dark Particles with Light Mediators

Millicharged particles in neutrino experiments

Improved Limits on Millicharged Particles Using the ArgoNeuT Experiment at Fermilab

Relationship of the relativistic Compton cross section to the electron's velocity distribution

Relationship of the relativistic Compton cross section to the momentum distribution of bound electron states. II. Effects of anisotropy and polarization

X-ray incoherent scattering total cross sections and energy-absorption cross sections by means of simple calculation routines

Relativistic Impulse Approximation in Compton Scattering

Hartree-Fock Compton Profiles for the Elements, At. Data and Nucl

Compton profiles and electronic properties of Nd

Analytical reconstruction of momentum density from directional Compton profiles

Electron momentum density and Compton profile by a semi-empirical approach

Evidence of instantaneous electron correlation from Compton profiles of crystalline silicon

Compton profiles and band structure calculations of IV-VI layered compounds GeS and GeSe

Proposal to determine the Fermi-surface topology of a doped iron-based superconductor using bulk-sensitive Fourier-transform Compton scattering

Quantum Field Theory and the Standard Model

Relaxed relativistic random-phase-approximation calculations of photoionization amplitudes and phases for the 4d subshell of xenon

Absolute cross sections for the photoionization of 4d electrons in Xe + and Xe 2+ ions

Application of the relativistic time-dependent density functional theory to the photoionization of xenon

Spin-orbit-interaction activated interchannel coupling in dipole and quadrupole photoionization

Photoionization of Xe and Rn from the relativistic random-phase theory

Strong Electron Correlation in Photoionization of Spin-Orbit Doublets

Particle acceleration at astrophysical shocks: A theory of cosmic ray origin

Particle Astrophysics

Cosmic Rays and Particle Physics

Electromagnetic neutrinos in laboratory experiments and astrophysics

Ionization of hydrogen by neutrino magnetic moment, relativistic muon, and WIMP

Atomic ionization of germanium by neutrinos from an ab initio approach

High-Energy Neutrino Astrophysics: Status and Perspectives

Grand Unified Neutrino Spectrum at Earth: Sources and Spectral Components

Solar Neutrinos: Status and Prospects

Solar models: Current epoch and time dependences, neutrinos, and helioseismological properties

Solar models with accretion. I. Application to the solar abundance problem

Gallium solar neutrino experiments: Absorption cross-sections, neutrino spectra, and predicted event rates

What do we (not) know theoretically about solar neutrino fluxes?

New solar opacities, abundances, helioseismology, and neutrino fluxes

New bounds on neutrino electric millicharge from limits on neutrino magnetic moment

Limits on Axion Couplings from the First 80 Days of Data of the PandaX-II Experiment

Neutrino physics with multi-ton scale liquid xenon detectors

Projected WIMP sensitivity of the LUX-ZEPLIN dark matter experiment

Can Nonstandard Neutrino Interactions explain the XENON1T spectral excess?

XENON1T Excess: Some Possible Backgrounds

Carbon 14 C and Tritium as possible background sources in XENON1T

An Introduction to Quantum Field Theory

Quantum Electrodynamics

Dirac-Fock atomic electronic structure calculations using different nuclear charge distributions, At. Data Nucl

Relativistic Quantum Theory of Atoms and Molecules

New Version: grasp2K relativistic atomic structure package

Breit Interaction in Multielectron Atoms

Convergence of the Breit interaction in self-consistent and configuration-interaction approaches

Theory of angular distribution and spin polarization of photoelectrons

Relativistic Dirac-Hartree-Fock Photon Incoherent Scattering Functions, At. Data Nucl. Data Tables

We acknowledge the helpful discussions with L. Singh and Henry T. Wong. This work was supported by the National Key Basic Research and Develop Program