The wavefunction reconstruction effects in calculation of DM-induced electronic transition in semiconductor targets

The physics of the electronic excitation in semiconductors induced by sub-GeV dark matter (DM) have been extensively discussed in literature, under the framework of the standard plane wave (PW) and pseudopotential calculation scheme. In this paper, we investigate the implication of the all-electron (AE) reconstruction on estimation of the DM-induced electronic transition event rates. As a benchmark study, we first calculate the wavefunctions in silicon and germanium bulk crystals based on both the AE and pseudo (PS) schemes within the projector augmented wave (PAW) framework, and then make comparisons between the calculated excitation event rates obtained from these two approaches. It turns out that in process where large momentum transfer is kinetically allowed, the two calculated event rates can differ by a factor of a few. Such discrepancies are found to stem from the high-momentum components neglected in the PS scheme. It is thus implied that the correction from the AE wavefunction in the core region is necessary for an accurate estimate of the DM-induced transition event rate in semiconductors.


Introduction
The existence of the dark matter (DM) has been well established on both theoretical and phenomenological grounds in astrophysics and cosmology. However, its nature still remains a mystery from a particle physical perspective. Among various theoretical hypotheses, the generic weakly interacting massive particles (WIMPs) have long been one of the prevailing particle DM candidates. Their typical masses ranging from GeV to TeV, and coupling strengths with the standard model (SM) particles at the weak scale, not only naturally account for the observed relic abundance in the context of thermal freeze-out, but also make the WIMP direct detection possible for the present-day technologies [1]. In the past three decades tremendous progress has been made in relevant field and the sensitivity of the underground detectors will eventually reach the solar neutrino background in the coming years, leaving only a limited unexplored parameter region for the direct detection of the WIMPs [2][3][4].
Facing this situation, both theorists and experimentalists begin to shift their interest to other new directions. The sub-GeV DM as an alternative consideration, with DM masses in the MeV to GeV range, has attracted increasing attention in recent years, for its theoretical motivations, and more importantly, its experimental feasibility. In the sub-GeV DM paradigm, the DM particle is expected to reveal itself by triggering ionization signals via the DM-electron scattering in semiconductor detectors (e.g., the germanium-based detector EDELWEISS [5] and SuperCDMS [6] consisting of both silicon and germanium crystals) with energy thresholds as low as a few eV, which translates to a DM mass detection threshold down to the MeV scale. To realistically quantify the electronic excitation rate in the crystal poses a challenge for theoretical estimation of relevant experimental sensitivity, since the outer-shell electrons participating in the excitation are delocalized in the crystal, and hence the isolated atomic description turns invalid.

JHEP01(2019)149
Step by step, several seminal studies help make concrete progress on this issue [7][8][9][10]. In refs. [8,10], where large momentum transfers are assumed, the initial (valence) states are simply approximated as corresponding free atomic orbitals, while the scattered ones are described by plane waves. The calculations are performed with relevant atomic Roothaan-Hartree-Fock (RHF) wavefunctions in ref. [10]. Another strategy is to calculate relevant electronic states under the framework of the density functional theory (DFT) [11,12], as adopted in ref. [9], in which electronic states are described with the band structures in the context of solid state physics. The authors of ref. [9] adopted software package Quantum Espresso [13] to perform the ab initio pseudopotenial calculations on silicon and germanium semiconductors, hence putting the estimation of electronic excitation rates onto a realistic ground, and encouraging a recent wave of studies on a wider range of novel materials for the detection of even lighter DM candidates [14][15][16][17][18][19][20][21][22].
However, in this study our interest is focused on the effectiveness of the pseudopotenial approaches adopted in the DFT implementation. While the periodic electronic wavefunctions are conveniently expressed with the plane wave (PW) basis set in a typical DFT computation of the electrons in crystalline solid, the rapidly varying wavefunction near the nucleus requires a huge number of PWs to be presented, even though they are irrelevant to most of the low-energy physical processes (e.g., the chemical bonding properties). The purpose of introducing pseudopotenial is thus to soften the oscillatory core wavefunction and reduce the number of PWs accordingly, without disturbing the electrons outside the core region. Although the pseudization of the valence electrons has proved an effective and efficient method in calculations of a wide range of electronic properties, the effectiveness of such well-adopted procedure needs to be scrutinized in the study of electronic excitation process, where transition energies are often deposited through deep atomic inelastic scattering processes. In particular, in order to excite a valence electron far beyond the band gap, the incident DM particle needs to penetrate into the core region where the pseudo wavefunction description may no longer suffice for an accurate analysis. Consequently, a realistic description of inner electrons becomes necessary. In fact, similar investigations have existed for long in literature, for instance, exploring discrepancy between the pseudo (PS) and the all-electron (AE) descriptions in optical properties of semiconductors [23], and in X-ray Compton profiles for semiconductors and metals [24,25]. Following these literatures in which the AE wavefunctions are restored from various reconstruction schemes, we use the term reconstruction effects referring to discrepancies in calculation originating from the PS and AE approaches.
For DM particles in the MeV-mass range, such investigations are even more imperative, because compared to the case of optical excitation, due to different dispersion relations, a much larger momentum transfer will be cost onto the massive particle given the same amount of deposited energy in solids. However, high-momentum components of the PW expansion that would response to such large transferred momentum are systematically abandoned in the pseudopotenial calculations, and this may result in a noticeable discrepancy between the two approaches. Common pseudopotenial-like methods include norm-conserving (NC) [26], ultrasoft (US) [27], and projector augmented wave (PAW) [28], etc. For example, in calculation of the DM-induced transition event rates, the US and PAW JHEP01(2019)149 pseudopotenials were adopted in ref. [9] and ref. [17], respectively. As a benchmark study, in this work we choose the PAW method as a concrete example, since the PAW scheme not only well describes electrons in the interstitial region, but also are capable of restoring the AE wavefunctions in the core region. Besides, it is convenient to make comparison within the same framework, and qualitative results can apply to other pseudopotenial methods.
This paper is outlined as follows. We begin section 2 by giving a summary of concepts and formulae associated with the calculation of the electronic transition rates induced by DM particles. In section 3 we first take a brief review of the PAW method, and then generate appropriate PAW datasets for the computation of the band structures. By inputting these PAW files into the ABINIT package, in section 4 we calculate the DM-induced excitation events with both the AE and PW approaches for silicon and germanium semiconductors, respectively. The ratios of the AE event rates to the PW ones are calculated and the reconstruction effects are explored. We conclude in section 5.
We use natural units = c = 1 in all formulae in this paper, while velocities are expressed with units of km · s −1 in text for convenience.

DM-induced electronic transition rates
In this section, we give a summary of the formulae concerning the DM-induced electronic transition event rates in the semiconductor targets. For simplicity, in this work we only focus on the scenario where the amplitude for DM particle scattering with a free electron is constant, i.e., M free (q) ∝ 1, and hence the DM-electron cross section σ χe = |M free | 2 µ 2 χe /π is also a constant (with µ χe = m e m χ / (m e + m χ ) being the reduced mass of DM-electron pair), independent of DM transferred momentum q and incident velocity w. As a result, the transition event rate for a DM particle with incident velocity w to excite a bound electron from level 1 to level 2 can be expressed in the following schematic form [9]: where the form factor encodes the excitation probability for the bound states in the electron system, with ψ 1 (x) and ψ 2 (x) being the wave functions of the initial and final electron states, respectively, and ∆E 1→2 denotes the relevant energy difference. Then after taking into account the DM halo profile, the excitation event rate from state 1 to 2 can be expressed as 3)

JHEP01(2019)149
where the bracket · · · denotes the average over the DM velocity distribution, ρ χ represents the DM local density, f χ (w,q) is the DM velocity distribution with unit vectorq as its zenith direction, g χ (w;q) ≡ w 2 f χ (w;q), and Θ is the Heaviside step function. While φq w is the azimuthal angle of the spherical coordinate system (q; w), the polar angle θq w has integrated out the delta function in above derivation. For the given q and ∆E 1→2 , function w min determines the minimum kinetically accessible velocity for the transition: In practice, the velocity distribution can be approximated as a truncated Maxwellian form in the galactic rest frame, i.e., f χ (w,q) ∝ exp − |w , with the earth's velocity v e = 230 km · s −1 , the dispersion velocity v 0 = 220 km · s −1 and the galactic escape velocity v esc = 544 km · s −1 .
When applied to the periodic lattice of a crystalline solid, the energy level can be labeled with discrete band index i, and crystal momentum k defined within the first Brillouin Zone (BZ), corresponding to the eigen wave functions of the Kohn-Sham (KS) orbitals. By taking advantage of the Bloch's theorem, these wave functions can be written as a product of the phase factor e ik·x and lattice periodic function u ik (x) as the following, where Ω and V denote the volume of the unit cell (UC) and of the whole crystal, respectively, and u ik (x) is represented as a linear combination of PWs denoted by reciprocal lattice vectors G's. For convenience, the periodic wave functions {u ik } are normalized within the UC, while {ψ ik } are normalized on the whole crystal. Thus, the square of the form factor form factor eq. (2.3) can be explicitly written as where the integral Ω d 3 x (. . .) is performed over the UC. Then, after inserting eq. (2.6) into eq. (2.3) and taking into account the two degenerate spin states, the total excitation event rate can be expressed as where the summations are over both the initial valence band states and the final conduction states. On occasions where only the bonding properties are of concern, the KS wave functions are usually calculated with the pseudopotential method to avoid invoking too many PWs for complete representation of the oscillatory wave functions near the nucleus, while keeping intact the description of the delocalized electrons outside the core regions. At variance with those cases, however, the excitation rate in eq. (2.7) explicitly depends on the true wave functions within the core region. Thus, it is tempting to give a quantitative comparison between the excitation event rates computed respectively from the PS and AE approaches. For this purpose, in this work we choose the PAW method for the benchmark study. For one thing, the PAW formalism is closely resemblant to the widely used US pseudopotential method [27], also generating soft pseudopotential and giving an even more efficient and reliable description of material electronic structure [29]. For another, within the same PAW framework, the true AE wave functions can be restored consistently from the built-in transformation. These two reasons make the PAW framework a natural host for the comparison between the PS and AE approaches concerning the calculations of the excitation event rates.

PAW method
In this section we will take a brief review on the PAW method proposed originally in ref. [28]. The purpose is to calculate relevant physical quantities through transforming the AE wave function problem into a computationally convenient PS one. To avoid confusion, it should be noted that in the context the term "all electron" (AE) is only referred to the true wave functions, in contrast to the pseudopotential method that smoothens the valence wave functions near the nucleus, rather than the scheme where all the core electrons are put on the equal footing in solving the Kohn-Sham (KS) equations. In fact, the PAW scheme is also based on the frozen-core approximation. If necessary, unfreezing of lower-lying core states is straightforward and convenient in the PAW method [29]. Imagine one expands the true AE KS single particle wave function |Ψ in terms of a set of true atomic wave functions {|φ a i } inside an atom-specific augmentation sphere Ω a centered at the atom site a, and a smooth wave function |ψ in the interstitial region, respectively, as the following, where the AE partial waves {|φ a i } are the eigen wavefunctions of the KS Schrödinger equation for an isolated atom, and expansion coefficients {|c i } along with smooth wave function |ψ are left to be determined. The index i stands for the main and the angular momentum quantum numbers (n, ℓ, m). From {|φ a i }, one can also generate a complete set of corresponding PS partial waves {|φ a i } and an accompanying set of projector functions {|p a i } inside the augmentation region Ω a . By construction they are required to fulfill the condition over the augmentation region. Outside the augmentation sphere Ω a , the auxiliary smooth partial waves {|φ a i } are constructed so that they are identical to the their AE counterparts {|φ a i }. Since the AE and PS partial waves coincide on the augmentation sphere, it is natural to make a smooth continuation of wave function |ψ inwards with smooth PS partial waves in the following manner, As a result, the initial AE wave function can be recast in terms of the PS one through the transformation From this particular form, the effective Hamiltonian for the PS wave function can be expressed accordingly asˆ which in turn yields the transformed eigenvalue equation for |Ψ , with the same eigenvalues as the original KS equation. Now it is the PS wave function |Ψ that appears in the modified KS equation eq. (3.6), which is expected to significantly facilitate the computation because the strong oscillating components have been separated out. Consequently, once the PS wave function |Ψ is obtained, its AE counterpart can be easily reconstructed via eq. (3.4). Practical applications show that two partial waves per orbital are usually sufficient to obtain accurate results [28,29]. Further details of the PAW method can be found in refs. [28,29]. Now we delve into the computational details of the generation of the PAW datasets for two most common semiconductor targets, silicon and germanium. In practice, we use the atompaw code [30] to obtain the sets of atomic AE and PS partial waves, along with the projectors in a procedure similar to the one originally suggested by Blöchl [28]. Two partial waves are generated for each angular momentum quantum number ℓ, and the GGA-PBE exchange-correlation functional [31] is chosen for the calculation. The atomic AE partial wave is expressed in the spherical coordinates as a product of a radial function times a spherical harmonic function: where v at KS [n] (r) is self-consistent atomic potential and ǫ nℓ is the energy eigenvalue. Based on these AE wave functions, corresponding radial PS partial wavesφ nℓ (r) and projector functions p nℓ (r) can then be constructed. Relevant results for silicon and germanium are summarized below.

Si
In generation of the PAW datasets of the silicon, the 3s and 3p electrons are taken as the valence electrons, while others in the closed shell are approximated as the frozen core. In figure 1 we present the partial waves and energy derivatives for the silicon atom, following the atompaw recipe. Specifically, we adopt a radius of the augmentation sphere r c = 2. Bohr, and two partial waves (so as the projectors) for each angular momentum quantum number ℓ = 0 and ℓ = 1 are generated, based on two reference energies of a bound and a scattering states. For ℓ = 0 (1), the reference energies are respectively the eigen energies of the bound state ǫ 3s = −0.795 Ry (ǫ 3p = −0.300 Ry) and the scattering state of the silicon atom ǫ s = 5. Ry (ǫ p = 4. Ry). In third row shown are the dimensionless logarithmic derivatives D ℓ for corresponding AE and PS atomic partial waves, which is defined as where ψ ℓ (ǫ, r) is either the AE partial wave function of eq. (3.8) at an arbitrary reference energy ǫ, or the PS counterpart from the transformed version eq. (3.6). For large enough augmentation radius r c , the quantity D ℓ (ǫ, r c ) uniquely determines the phase shifts of the partial waves in scattering theory. Thus the coincidence of the logarithmic derivatives implies that the constructed PAW Hamiltonian well recovers the scattering properties of a single atom, indicating a good performance of the PAW implementation in the energy range of interest.

Ge
In a parallel fashion, the corresponding partial waves and energy derivatives for the germanium atom are presented in figure 2. Similarly to case of silicon, the 4s and 4p orbitals are treated as the valence states, while the closed shell is approximated as the frozen core.
With radius of the augmentation sphere r c = 1.9 Bohr, two reference energies of the bound orbital ǫ 4s = −0.862 Ry (ǫ 4p = −0.286 Ry) and unbound state ǫ s = 1. Ry (ǫ p = 1.2 Ry) are chosen to generate the ℓ = 0 (1) partial waves and projectors. For germanium, it is also shown that the logarithmic derivatives of the AE and PS atomic partial waves coincide with each other quite closely.

Reconstruction effects on transition event rate
In this section, we first utilize the PAW datasets generated from the atompaw code to calculate the band structures of silicon and germanium, as well as both the corresponding AE and PS eigen wave functions with the ABINIT package [32], from which we then analyze and compare the calculations of relevant DM-induced transition event rates. Fifteen adjacent bands from both the valence and conduction states across the band gaps are presented in figure 3 for silicon and germanium, respectively, along a specific Brillouin zone path L-Γ-X-Γ. for convenience, and the band-gap energies E g , namely, the energies of the conduction band minimums (CBM) are calibrated to the experimental values for silicon (1.12 eV) and germanium (0.67 eV), respectively, using the scissor operation. In calculation, we adopt an energy cutoff E cut = 96 Ry, a value sufficiently large to obtain converged band structures in the PAW framework. This energy cutoff translates to a cap on the radius of the reciprocal vector G at G max = √ 2m e E cut . Moreover, the fast Fourier transformation (FFT) of the electronic density, or square of the wave functions, which makes its appearance in eq. (2.7), requires PW components up to a reciprocal vector radius twice as large as G max JHEP01(2019)149 in order to avoid the wrap-around errors. So in practice the charge density is presented in a box consisting of (2G max a/π) 3 ≈ 64 3 grids. In this work we use the experimental lattice constants for silicon a Si = 10.263 Bohr and germanium a Ge = 10.619 Bohr, respectively.
In addition, to illustrate the source of possible error from the pseudopotential method for the calculation of transition event rate eq. (2.7), in figure 4 we make an elaborate comparison of the valence electron densities between the AE and PS approaches inside the UC. For concreteness, we take the VBM (at the Γ point k = 0) for instance, and choose three specific paths within the UC of the two crystals, which are x = r 3 a 3 (in the upper row), x = 1 4 a 1 + 1 4 a 2 + r 3 a 3 (in the middle row), and x = 1 2 a 1 + 1 2 a 2 + r 3 a 3 (in the bottom row) with 0 ≤ r 3 < 1, respectively, where the primitive vectors {a 1 , a 2 , a 3 } = a 0, 1 2 , 1 2 , a 1 2 , 0, 1 2 , a 1 2 , 1 2 , 0 are expressed in Cartesian coordinates and in terms of the lattice constant a. Corresponding electronic densities are normalized within the UC (see eq. (2.5)) with condition |u VBM (x)| 2 Ω d 3 r = 1. As shown in the first and second rows in figure 4, the AE electronic density (red) features rapid oscillations when the path passes through the nuclei, while outside the augmentation radius the PS density (blue) well coincides with the AE one. When the path does not intersect the augmentation sphere, as shown in the third row, the two wave functions are completely equal.
Another observation from figure 4 is that the oscillatory behavior of the wave function near the nucleus is so remarkable that the present real grids fall short of resolution to continuously depict the variation details. If the kinetics allows for a large transferred momentum in transition, these highly oscillating wave functions may couple with the highfrequency phase factor terms {e −iG·x } in the expression of transition rate eq. (2.7), and hence more G-vector terms can contribute to the total transition rate. The maximum transferred momentum for the given energy difference E i ′ k ′ − E ik and DM incident velocity JHEP01(2019)149 Figure 4. Comparison of the valence electron densities of the highest valence band state at Γ, for silicon (left) and germanium (right), determined from the AE (red ) and PS (blue) approaches within the PAW framework, respectively. In the upper row, the densities are along the path x = r 3 a 3 (0 ≤ r 3 < 1) inside the UC, while in middle and the bottom rows presented are those along the paths x = 1 4 a 1 + 1 4 a 2 + r 3 a 3 and x = 1 2 a 1 + 1 2 a 2 + r 3 a 3 , respectively. See text for details.

JHEP01(2019)149
w can be drawn from the Heaviside function Θ in eq. (2.7) as the following, In contrast to the case of optical excitation, this massive particle kinetic relation translates to an average momentum transfer q of O (c/w 0 ) ≈ O 10 3 larger given the same transition energy, where w 0 = 300 km · s −1 represents the typical speed of the galactic DM particles. As a consequence, a larger modulus of G-vectors G is accordingly required to resolve the short-range details in this deep inelastic scattering process. To get a sense, we take the excitation channel VBM→CBM in silicon crystal for instance, where a 50 MeV-DM particle with a typical velocity w 0 = 300 km · s −1 can produce a transferred momentum as large as 0.099 MeV, which corresponds to an energy cutoff E cut ≈ 176 Ry, far exceeding the typical value of 30 ∼ 60 Ry in pseudopotential PW calculations. While these large-G phase factors {e −iG·x } decouple from the smooth PS wave functions, they may still resonate with the rapidly oscillatory AE ones. Therefore, in order to investigate the AE reconstruction effects on the calculation of DM-induced transition rate, we need to choose an appropriate set of parameters in a sense that the maximum momentum transfer should be compatible with the resolution of the real grid adopted in computation. For illustration, we present the largest modulus of integer number vectors for the diamond structure, i.e., n G = q max (E g ) a/(2π) 1 versus two parameters: DM mass m χ and incident velocity w, for silicon and germanium crystals in the upper row in figure 5. The kinetic energy cutoff E cut = 96 Ry adopted in this study corresponds to the constraint n G ≤ 32 (in black dashed line) in the ABINIT implementation of PAW. That is to say, the resolution associated with the cutoff E cut = 96 Ry can only describe the electronic excitation processes in colored parameter regions below the black dashed lines in figure 5, while in the white areas corresponding excitations are kinetically forbidden.
After determining the parameter region compatible with the given energy cutoff, now we can compute the DM-induced electronic excitation rates for the silicon and germanium semiconductor detectors. In calculation we truncate the energy transfer E e at 10 eV, which corresponds to including only the transition events from the 4 valence bands to the 7 lowest conduction bands. Such constraint is presented in white dotted lines in the upper row in figure 5 for silicon and germanium, respectively. Parameters (m χ , w) above the white line can induce excitation events with transition energies above 10 eV, and hence more unoccupied energy bands are required for exact calculation. In addition, although the transition rate eq. (2.7) depends on the orientation of the crystal target with respect to galaxy, for simplicity however, in calculation the velocity distribution is approximated as an isotropic one for parameter (m χ , w), i.e., g χ (v,q) = (4π) −1 δ (v − w). Besides, following the strategy in ref. [9], we equivalently express eq. (2.7) in terms of transferred momentum q and energy transfer E e in order to make the scan of the parameter space computationally

JHEP01(2019)149
more efficient. In short, after inserting the monochromatic velocity distribution the total transition rate for (m χ , w) can be recast as In realistic computation, several numerical modifications are made to eq. (4.2) as follows.
• The integrals over q and E e are replaced by the summations over uniform discrete bins of q and E e . We use the bin widths ∆E = 0.058 eV and ∆q = 78.6 eV, and the integrand is evaluated at the central value in each relevant bin. The delta functions are discretized to a normalized step function. For example, if a specific parameter E e falls into the n-th energy bin, then the function δ ( with E n being the central value of this bin. • As a standard DFT numerical procedure, the integrals of the continuous k-points in the first BZ are replaced by the summations over a uniform discrete mesh of representative k-points in the following manner: 4) where N k is the number of k-points sampled in the first BZ. In this study, we use a homogeneous set of 6 × 6 × 6 k-points.
can be understood as the Fourier transformation of the square term u * i ′ k ′ (x) u ik (x) within the UC, which is practically performed using the discrete FFT technique, so the reciprocal G-vectors have a natural truncation consistent with the grid number in the real space. As mentioned above, the modulus number of the G-vectors n G can extend to a value no more than 32 in our analysis.
In the middle row of figure 5, we present the calculated transition event rates per kilogramday from the PS wavefunctions for silicon (left) and germanium (right), respectively. A DM-electron cross section σ χe = 10 −40 cm 2 is assumed, and only the parameter regions below the transition energy constraint (white dotted line in the upper row) and the G-  are the comparisons between the theoretical excitation event rates derived from the AE and PS approaches for the two crystals, in terms of the ratio of the AE event rate to its PS counterpart for each parameter pair m χ and w, i.e., η ≡ R AE (m χ , w) /R PS (m χ , w). In large part of the scanned region, the two calculated even rates are found to be quite consistent. Interestingly, however, in areas where large momentum transfer (or equivalently, large n G ) can be induced, the estimates based on the AE and PS approaches begin to differ with each other, and this trend seems to continue into larger G-vector parameter region. Therefore, the actual electronic distribution near the nucleus is of remarkable relevance for the estimation of the DM-induced excitation rate, especially in the area where a large momentum transfer is favored in the excitation process.
In order to demonstrate the origin of such deviation, in figure 6 we present the comparison of the form factor F (q, E e ) in eq. (4.2), between the AE and PS approaches for silicon. This non-dimensional quantity is convenient for direct contrast because all the integrand function values in pixels are accounted with equal weight when summed up to obtain the total event rate. The momentum transfer q is expressed in terms of the silicon reciprocal lattice 2π/a Si . We choose two benchmark points representative of two distinct scenarios in the bottom row of figure 5, i.e., when the AE-based estimate of event rate differs with its PS counterpart, and when the AE-based estimate of event rate agrees with its PS counterpart.
For the first case, we take the parameter pair m χ = 70 MeV, w = 150 km · s −1 as example. At this point one has η = 2.01, and the kinetics confines the integration in eq. (4.2) within the area in the red contour in figure 6, only where the Heaviside step function yields non-zero value. In addition, an important observation from figure 6 is that as we expected, the AE wavefunctions near the nucleus, reconstructed from the PAW transformation eq. (3.4), indeed play a non-negligible role in transitions where large transferred momentum participates. While the PS approach gives an accurate description of the exci-

JHEP01(2019)149
tation process in the parameter area below q ≈ 5 (2π/a Si ), (which is referred to as small-q region, and the rest area is referred to as large-q region, for convenience) the relevant form factor falls off so rapidly towards the large-q region that the PS scheme may cause remarkable error in the calculation of event rates. Comparing the AE case (left) to the PS one (right), one finds that the large ratio η can be attributed to the fact that the contour on one hand covers only a small portion of the parameter area where the PS-based description remains effective, and on the other it contains a sizable part of the large-q parameter area. These two factors effectively enchance the contrast.
As for the other case, we pick the parameter pair m χ = 10 MeV, w = 400 km · s −1 (η = 0.85) for illustration. Since half of the kinetically allowed region (contoured in blue line) overlaps the small-q area, based on the above reasoning, it is understandable that the corrections from the AE approach are not so remarkable under this circumstance. Moreover, the galactic DM velocity distribution also imposes a kinetic constraint on the E eq plane. Actually, given the truncated Maxwellian distribution introduced in section 2, and the fact that the w min (q, E e ) reaches its minimum as m χ → ∞ (see eq. (2.4)), the minimal kinetically accessible transferred momentum q for recoil energy E e equals E e / (v esc+ v e ), which corresponds to the yellow dashed line presented in figure 6. The region below the line is kinetically forbidden. The implication is that large recoil energy requires large momentum transfer accordingly, and thus an accurate description of the AE electronic wave functions becomes necessary.

Conclusions and discussion
In this work, we have investigated the implications of the PAW reconstruction procedure on the estimation of electronic transition event rate in two most common semiconductor targets used in the DM direct detection, namely, silicon and germanium semiconductors. By use of the PAW method, one can restore the AE electronic wavefunctions from the PAW PS wavefunctions without much extra computational loads. To this end, we first utilize the atompaw code to generate appropriate PAW datasets including the PAW pseudopotential, and input these files into the ABINIT package to compute the electronic structure and relevant AE and PS wavefunctions. Based on the wavefunctions, we then calculate relevant DM-induced excitation event rates and make comparisons between the AE and PS approaches. It is found that the PS approach turns inapplicable in describing excitation processes involving large transferred momentum, because the oscillating components that would couple with large momentum transfer are deliberately spared in calculation. Even in the limited parameter region scanned in our study, the two calculated event rates are found to differ with a factor of two. So this work serves as a caveat in relevant calculations that an accurate estimate of the DM-induced transition event rate should better be implemented with the AE approaches.
Similar experimental strategies can be extended to the detection of the sub-MeV DM particles. For instance, in ref. [17] the authors discussed the theoretical prospects for detecting sub-MeV DM particles by use of the three-dimensional Dirac materials with typical band gap E g ∼ meV. Discussion of the AE reconstruction effects can also be applied JHEP01(2019)149 parallelly to this case. However, considering the dependence of maximum momentum transfer on the DM mass at the sub-MeV scale (see eq. (4.1)), the high-momentum PWs from reconstruction do not contribute to the transition event rates, so in such case the reconstruction effects are expected to be insubstantial.
Finally, a thorough estimate of the total event rate should involve convoluting with the realistic DM velocity distribution, and thus requires a broader recoil energy range and a larger energy cutoff (see figure 5), which unfortunately, is beyond the capabilities of our cottage-industry work style. As a more efficient solution, computation of the crystal form factor should be processed within the well-established software packages, e.g., ABINIT in this study. We leave it for future work.