Re-understanding of the deformation potential constant in the single crystal silicon

Since the deformation potential (DP) theory was proposed by Bardeen and Shockley in 1950 [1], it has been applied to describe the interaction between electron and phonon [2, 3, 4] and evaluate transport properties of semiconductors successfully [5, 6, 7]. In particular, its combination with density functional theory (DFT) provides a favorable mean for explaining experiments and predicting new materials in the thermoelectric community. However, there are still some disagreements about the calculation for the DP constant of bulk materials. The focus of disagreements is that no precise reference energy level is found for obtaining the absolute band energy of the different deformation structures.

As is well known, the effective one-electron Schrödinger (or Dirac) equation with respect to the eigenstates and eigenvalues ($\Psi_{i}$ and $\varepsilon_{i}$) is expressed as [8]

$$[K_{op}+V_{eff}(\textbf{r})]\psi_{i}(\textbf{r})=\epsilon_{i}\psi_{i}(\textbf{r}),$$

(1)

where the effective Coulomb potential operator is the sum of the exchange correction potential ($\mu_{xc}$) and Coulomb potential ($V_{c}$)

$$V_{eff}(\textbf{r})=\mu_{xc}(\textbf{r})+V_{c}(\textbf{r}).$$

(2)

The $V_{c}$ is given by

$$V_{c}(\textbf{r})=\int\frac{\rho(\textbf{r}^{\prime})}{|\textbf{r}-\textbf{r}^{\prime}|}d\textbf{r}^{\prime}-\sum_{\alpha}\frac{Z_{\alpha}}{|\textbf{r}-\textbf{R}_{\alpha}|}.$$

(3)

Obviously, there are the Coulomb singularities at the electron and nucleus sites ($\textbf{r}^{\prime}$ and $\textbf{R}_{\alpha}$). In order to solve the Eq. (1), there are many methods proposed to deal with the singularities. However, it is impossible to avoid introducing an uncertain constant when solving the Poisson’s equation in any case. In fact, introducing the uncertain constant ignores the effects of both volume deformation and atomic species on the periodic potential. As a result, it makes the band energies of different materials and of the same material with different volumes, unable to be compared quantitatively, although the influence on the total energy can be canceled by special mathematical processing [9, 10, 11].

For the calculation of the DP constant of one material, it is not necessary to consider the influence of atomic species. Therefore, finding a physical quantity that is not affected by the volume deformation as the reference energy level is the key to calculating the DP constant. At present, there are average potential [12, 13, 14, 15] and core state level [16, 17] alignments to define the reference energy level. The former is also known as the supercell method, where the band structures are aligned according to the average potential from the core electrons in the supercell composed of the compressed and expanded parts. The latter is based on the assumption that the energy level of the core electron is hardly affected by the volume deformation.

Usually, the core state level alignment method is adopted to explore thermoelectric properties due to the requirement of less computational resources. In fact, the core state levels are sensitive to the volume deformation, and thus have non-negligible DP values [16]. Consequently, the DP constants corrected by the core state level alignment are small and thus lead to much high mobility deviating from the experimental results. The disagreement between theory and experiment cannot be completely attributed to the theoretical calculation, because the recognized accurate experimental method for determining the DP constant has not been proposed.

Experimentally, the present determinations of the DP constant are only performed indirectly with a significant amount of assumptions. Therefore, the experimental DP constant for the same material has great differences. Taking GaAs as an example, Nolte et al. used the impurity level as the reference level to obtain the DP value of -9.3 eV for the conduction band minimum (CBM) [18], however, its value is 7.0 eV from the mobility measurements at low temperature [19]. Interestingly, a value of $16.0\pm 0.5$ eV for CBM is obtained from high-temperature mobility [20]. Moreover, Pfeffer et al. [21] inferred from the optical absorption in Se-doped GaAs bulk [22] that the DP value is as high as 15.7 eV. A much higher value of 20 eV is accepted and adopted in the Supporting Information of the recent literature [23]. In a word, the experiment has not so far given satisfactory results about the DP constant for GaAs.

Theoretically, the DP constant of a material also varies greatly. In the original literature of DP theory [1], Bardeen and Shockley first deduced 11.3 eV for the DP value of the valence band maximum (VBM) of cubic silicon from measured mobility data in the polycrystalline silicon. Latter, the values -10.2 eV [24], -7.9 eV [25], and 2.05 eV [17] are proposed by different researchers. It was claimed in some literatures [14, 16, 26, 27] that the obtained $\sim$1.5 eV along the [100] direction is very consistent with the experimental value $1.8\pm 0.7$ eV [28] for the p-type silicon. However, it is puzzled that we have not found relevant experimental data and only see the value of DP constants ($3.3\pm 0.7$ eV) for the n-type silicon in Ref. [28]. Furthermore, it is worth mentioning that the experimental data is obtained in the heavily doped n-type silicon with the high concentration of $5\times 10^{21}cm^{-3}$ and its Fermi level is located above the CBM. Therefore, the measured DP constant is for the Fermi level above the CBM, however, the theoretical DP constant is for the CBM [14]. This implies that the theory and experiment are not consistent, even may deviate.

In the present work we first redefined the relevant physical parameters in the mobility formula, based on the analysis and summary of previous theoretical and experimental data on the DP constant. Then we present a strategy to evaluate the correction of the DP constants of materials with the aid of all-electron first-principles approach. Finally, we employ the modified mobility formula to calculate the mobility of the single crystal silicon for verifying the effectiveness of the strategy.

According to the DP theory, the carrier mobility $\mu$ of non-polar crystals with the properties of isotropy and non-degeneracy is expressed as [1]

$$\mu_{e(h)}=\frac{2\sqrt{2\pi}e\hbar^{4}\rho v^{2}}{3(k_{B}T)^{3/2}m^{5/2}_{e(h)}(\lambda_{e(h)})^{2}},$$

(4)

where $k_{B}$, $e$ and $\hbar$ represent, respectively, the Boltzmann constant, the electronic charge and the reduced Planck constant, and $\rho$,$v$, $m_{e(h)}$ and $\lambda_{e(h)}$ are the mass density, the longitude acoustic velocity, the effective mass and DP constant for the n- or p-type (electron/hole) carrier (namely, for the CBM or valence band maximum (VBM)).

The $v$ of the isotropic polycrystal is govern by [29]

$${{\textit{v}}}={{\left(\frac{3{{B}}+4G}{3\rho}\right)}^{1/2}},$$

(5)

where $B$ and $G$ are the bulk modulus and the shear modulus, respectively. Sometimes the following expression, neglecting the shear strain contribution to the longitude acoustic velocity, is adopted.

$${{\textit{v}}}={{\left(\frac{{{B}}}{\rho}\right)}^{1/2}}.$$

(6)

This is reasonable because small shear strain does not induce the volume change and therefore cannot generate effective DP value.

The longitude acoustic velocity $v^{\beta}$ along the $\beta$ direction in the single crystal is approximatively expressed by the elastic constant $c_{\beta}$ and the mass density,

$$v^{\beta}=\sqrt{\frac{c^{\beta}}{\rho}}.$$

(7)

For the anisotropic single crystal, the effective mass of electron/hole along the $\beta$ direction ($m^{\beta}_{e(h)}$) can be obtained by calculating the band effective mass at the CBM/VBM

$$\frac{1}{m^{\beta}_{e(h)}}=\frac{\partial^{2}\varepsilon}{\hbar^{2}\partial(k^{\beta})^{2}}|_{\varepsilon=\varepsilon_{CBM(VBM)}},$$

(8)

where $\epsilon$ and $k^{\beta}$ stand for the energy of band edge and the wave vector along the $\beta$ direction.

We define three types of DP constants

$$\lambda_{e(h)}^{X}=\frac{\mathrm{d}\varepsilon}{\mathrm{d}\delta^{X}}|_{\varepsilon=\varepsilon_{CBM(VBM)}},$$

(9)

where $\delta^{X}$ represents different types of strains ($X=L,S\ or\ V$ for the uniaxial, biaxial or volumetric strain), $\varepsilon_{CBM(VBM)}$ denotes the energy at the CBM/VBM

	$\displaystyle\delta^{L}=\frac{L-L_{0}}{L_{0}},$		(10)
	$\displaystyle\delta^{S}=\frac{S-S_{0}}{S_{0}},$		(11)
	$\displaystyle\delta^{V}=\frac{V-V_{0}}{V_{0}},$		(12)

where $L_{0}$, $S_{0}$ and $V_{0}$ ($L$, $S$ and $V$) are length, area and volume of crystals without (with) the strains. Note that the DP constant $\lambda^{\beta}_{e(h)}$ can be obtained by applying the uniaxial strain, and thus is equivalent to the $\lambda^{L}_{e(h)}$. Also, we call $\lambda^{X}_{e(h)}$ ($L$, $S$ and $V$) as the uniaxial, biaxial and volumetric DP constant, respectively.

Usually, the CBM/VBM of the crystal with high symmetry has many equivalent points that possess the same energy and commonly determine the transport properties. The points can be divided into two categories (copoint or non-copoint degeneracy) based on whether they are at the same position. For instance, the n-type single crystal silicon has six non-copoint degenerate CBMs which symmetrically locate at the a*, b* and c* axes of the reciprocal space. Each two CBMs on the same axis have the same electronic structural properties along the same direction.

Each CBM has an ellipsoidal isoenergetic surface near it, and independently participates in electrical transport. Moreover, each CBM exhibits anisotropic transport behavior because the effective mass varies along different directions of the ellipsoidal isoenergetic surface. However, the total transport of the six CBMs is isotropic. For ith CBM, we assume that its contribution to electrical conductivity along one special $\beta$ direction ($\sigma_{i}^{\beta}$) is equivalent to the contribution of a spherical isoenergetic surface with the effective mass of the ellipsoidal isoenergetic surface along the $\beta$ direction, rather than the assumption of constant relaxation time. Thus, the total electronic conductivity along the $\beta$ direction can be expressed as

$$\sigma^{\beta}=\sum_{i=1}^{N_{v}}{\sigma^{\beta}_{i}}=\sum_{i=1}^{N_{v}}{n^{\beta}_{i}qu^{\beta}_{i}}={nqu^{\beta}},$$

(13)

where $N_{V}$, $u^{\beta}$ and n are the number of band degeneracies at band edge, the total real mobility and the total real carrier density that is described by

$$n=\sum_{i=1}^{N_{v}}{n_{i}}=\frac{V_{c}}{2\pi^{2}}(\frac{2m_{i}}{\hbar^{2}})^{\frac{3}{2}}\varepsilon^{\frac{1}{2}}f,$$

(14)

Here, $n_{i}$ is the real carrier density supported by the ith band degeneracy, $\varepsilon$, $V_{c}$ and $f$ are the energy of band edge, the volume of unit cell and the Fermi-Dirac distribution function, and $m_{i}$ is the local DOS effective mass of ith band degeneracy. The $m_{i}$ is determined by the shape of the isoenergetic surface near the band edge. The $m_{i}$ of ith spherical isoenergetic surface is band effective mass at CBM/VBM, which is equal to the effective mass along arbitrary $\beta$ direction ($m_{i}^{\beta}$). For elliptical isoenergetic surface, it can be written as

$$m_{i}=\sqrt[3]{m_{i}^{x}m_{i}^{y}m_{i}^{z}},$$

(15)

where $m_{i}^{x},m_{i}^{y}$ and $m_{i}^{z}$ are the effective masses along the three principal axes of the ellipsoidal isoenergetic surface. The $m_{i}$ is the result of equating the volume of an ellipsoid to a sphere. Mathematically speaking, it is the geometric mean of the effective masses along the three direction. The $m_{i}$ expression for the case of anisotropy has significant limitations. For example, when the effective mass along the one direction is infinitely close to zero, the $m_{i}$ is close to zero. Obviously, this is not in line with the facts. Therefore, it seems more reasonable to use geometric mean of the effective masses to describe the DOS near band edge.

$$m_{i}=\frac{m_{i}^{x}+m_{i}^{y}+m_{i}^{z}}{3},$$

(16)

In Eqs.(13), $n^{\beta}_{i}$ and $u^{\beta}_{i}$ are the equivalent carrier density and the equivalent mobility of ith band degeneracy given by

$\displaystyle n^{\beta}_{i}$

$\displaystyle=\frac{V_{c}}{2\pi^{2}}(\frac{2m^{\beta}_{i}}{\hbar^{2}})^{\frac{3}{2}}\varepsilon^{\frac{1}{2}}f,$

(17)

$\displaystyle\mu_{i}^{\beta}$

$\displaystyle=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3(k_{B}T)^{3/2}(m_{i}^{\beta})^{5/2}(\lambda^{\beta})^{2}}.$

(18)

Here, $m^{\beta}_{i}$ represents the effective masses of ith band CBM/VBM along the same $\beta$ direction. Substituting Eqs. (14), (17) and (18) into (13), the total real mobility is obtained

$\displaystyle\mu^{\beta}$

$\displaystyle=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3(k_{B}T)^{3/2}(m_{[m]}^{\beta})^{5/2}(\lambda^{\beta})^{2}}.$

(19)

We define the $m_{[m]}^{\beta}$ as the mobility effective mass along the $\beta$ direction,

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{N_{v}}\frac{1}{m_{i}^{\beta}}}{\sum\limits_{i=1}^{N_{v}}(m_{i})^{\frac{3}{2}}}.$

(20)

When the $N_{v}$ band degeneracies have the same isoenergetic surfaces, one can get the more common variation of Eq. (20),

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}=\frac{\frac{1}{N_{v}}\sum\limits_{i=1}^{N_{v}}\frac{1}{m_{i}^{\beta}}}{\frac{1}{N_{v}}\sum\limits_{i=1}^{N_{v}}(m_{i})^{\frac{3}{2}}}=\frac{\frac{1}{m_{c}}}{m_{i}^{3/2}}=\frac{1}{m_{c}m_{i}^{3/2}}.$

(21)

The $m_{c}$ can be called the conductivity effective mass. In fact, it is the mathematical harmonic mean. For the weak anisotropy, the ratio of $m_{c}$ and $m_{i}$ is close to 1. Thus, the other vibration[30] of Eq. (20) is

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}\approx\frac{1}{m_{i}^{5/2}}.$

(22)

Moreover, when the same isoenergetic surfaces locate at the same position, intervalley scattering can proceed smoothly like intravalley scattering. Therefore, the Eq. (18) should be writen as

$\displaystyle\mu_{i}^{\beta}$

$\displaystyle=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3N_{V}(k_{B}T)^{3/2}(m_{i}^{\beta})^{5/2}(\lambda^{\beta})^{2}}.$

(23)

and the Eq. (20) has a variant formula

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{N_{v}}\frac{1}{m_{i}^{\beta}}}{N_{V}\sum\limits_{i=1}^{N_{v}}(m_{i})^{\frac{3}{2}}}=\frac{1}{N_{V}m_{c}m_{i}^{3/2}}$

(24)

For n-type silicon, when the $\beta$ direction is parallel to the $a^{*}$ axis, the $m^{\beta}_{1}$=$m^{\beta}_{4}$, $m^{\beta}_{2}$=$m^{\beta}_{5}$ and $m^{\beta}_{3}$=$m^{\beta}_{6}$ correspond to $m_{l}$, $m_{t}$ and $m_{t}$ calling the longitude and transverse effective masses. Six $m_{i}$ values is equal and it is given by

$\displaystyle m_{i}=\sqrt[3]{m_{t}m_{t}m_{l}}.$

(25)

or

$\displaystyle m_{i}=\frac{m_{t}+m_{t}+m_{l}}{3}.$

(26)

Note that the effect of the intervalley scattering on the mobility is not taken account into in n-type silicon, it is because the CBMs possess vastly different momentum (namely locate at different sites), which hinders the intervalley scattering.

Different from the n-type silicon, the p-type silicon has three VBMs converging at $\Gamma$ point. The three degenerate bands contain two heavy bands and one light bands, so $m_{i}^{\beta}$ (i=1, 2, 3) is labled by $m_{hh}$, $m_{hh}$ and $m_{lh}$ respectively calling heavy and light band effective masses. Moreover, the carrier transition occurs between the same momentums and between the same energies, namely, the scattering near the VBMs is approximatively elastic. We assumed that the intervalley scattering is as well as the intravalley scattering. Therefore, the scattering rate is three times the case without intervalley scattering, and the hole mobility of the single crystal is given by

$$\mu_{h}^{\beta}=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3(k_{B}T)^{3/2}(m_{[m]h}^{\beta})^{5/2}(\lambda^{\beta}_{h})^{2}}.$$

(27)

With the help of Eq. (28), the $m_{[m]h}^{\beta}$ is roughly expressed as

$\displaystyle\frac{1}{(m_{[m]h}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{3}\frac{1}{m_{i}^{\beta}}}{3\sum\limits_{i=1}^{3}(m_{i})^{\frac{3}{2}}}=\frac{1}{3m_{c}m_{i}^{3/2}}.$

(28)

In fact, its accurate expression is

$\displaystyle\frac{1}{(m_{[m]h}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{3}(m_{i}^{\beta})^{1/2}}{(\sum\limits_{i=1}^{3}(m_{i})^{3/2})^{2}}.$

(29)

It is worth mentioning that the $\lambda^{\beta}$ is replaced by the $\lambda^{V}$ in the fact calculations because it is difficult to accurately calculate the deformation potentials of different degenerate band edges. Therefore, the anisotropy of the mobility comes from the contribution of the anisotropy of effective mass and elastic constant that respectively characterize the properties of electrons and phonons.

The relative total energies $\Delta E$ of the strained silicon with reference to the unstrained structure are all positive values shown in Figs. 1(a)–1(c). This indicates that the relative error for calculated lattice constant should be less than 1$\%$. It is expected that the calculated lattice constant 5.46 Å for Si agrees well with the experimental value 5.42 Å [31].

The improvements of total energies for the strained structures are related to the breaking of symmetry and the interaction between atoms. Compared to volumetric strain, the uniaxial and biaxial strains have greater impacts on total energy of the system, because the both lead to the reduction of structural symmetry. The uniaxially and biaxially strained structures belong to space group ($I4_{1}$, No.141), nevertheless, the volumetrically strained structure maintains the same structural symmetry ($Fd\overline{3}m$, No.227) with the unstrained structure. At the same strain strength, the compression more than the tensility can affect the total energy of silicon, because the repulsive interactions between atoms are more sensitive to the relative positions than the attractive interactions. For instance, the $\Delta E$ at the $-1\%$ uniaxial strain is 0.0025 eV obviously larger than 0.0012 eV at the $1\%$.

The calculated bandgap of unstrained silicon is 1.229 eV, close to the experimental value of 1.17 eV at 0 K [32]. The quit small bandgap deviating is mainly due to the calculated lattice constant slightly larger than the experimental data. The good agreement is because the band structure in our work is obtained within the modified Becke-Johnson exchange potential and the generalized gradient approximation for the correlation (mBJ-GGA) [33], using the full potential method as implemented in the WIEN2k code [34]. However, the mBJ-GGA strategy has a shortcoming that no exchange and correlation energy functional is defined [35], so that it cannot be used for all calculations in connection with the total energy of the system, such as total energy and geometric optimization calculations. As a consequence, the total energies in this work are obtained by the standard GGA strategy.

It is found in Figs. 1(d) and 1(e) that under uniaxial/biaxial strain the bandgap is positively related to the compressive volume, and negatively to the expansive volume. It is because both the compressive and tensile strains have different effects on the conduction and valence bands. However, the compressive and tensile strains have almost similar effects on the band edge so that the bandgap and volumetric strain are positively correlated throughout the overall region (see Fig. 1(f)). Without any energy corrections in Figs. 1(g)–1(i), we found the energy of the band edge goes down as the volume increases, and the standard linear relationship between the band edge energy and the volumetric tensile/compressive strain. However, the CBM energy is more sensitive to the uniaxial tensile strain, on the contrary, the VBM energy is more sensitive to the uniaxial compressive strain, which is also the case in the biaxial strain. The corresponding DP constants are summarized in Table 1. Under compressive strain, the $|\lambda_{VBM}^{L}|$ is as high as 14.90 eV that is about seven times of $|\lambda_{CBM}^{L}|$, and is twice $|\lambda_{VBM}^{L}|$ under tensile strain. Moreover, the overall $|\lambda_{VBM}^{X}|$ is larger than the overall $|\lambda_{CBM}^{X}|$. Especially, the difference between the overall $|\lambda_{VBM}^{L}|$ and $|\lambda_{CBM}^{L}|$ reaches as large as 4.59 eV.

Previous studies typically used the core state energy as the reference energy level to correct the band edge energy. In pseudopotential method, the use of pseudopotential to handle core electrons makes it impossible to obtain accurate core state energy. Therefore, the deep valence (DV) state energy is adopted as the reference energy level [36], and the DV state energy and the band edge energy are obtained by the analyses of the density of states (DOS) and the band structure along the high-symmetry points containing VBM and CBM, respectively. Due to differences in calculation methods and parameters, the DOS and band energies usually do not correspond perfectly. This increases the uncertainty of the calculation results on the DP constants. In this work both the DV state energy and the band edge energy are from the self-consistent calculation with very dense k points (30000 points in total). We take the average energy of the lowest energy band as the DV state energy. Figs. 1(j)–1(l) show the DV state energy and the strain have a negative linear relationship. The influence of the strain types on the DV state energy can be almost ignored owing to small differences between the $\lambda_{dv}^{L}$, $\lambda_{dv}^{S}$ and $\lambda_{dv}^{V}$ (-11.66 eV, -11.68 eV and -11.69 eV). More interestingly, the $\lambda_{dv}^{X}$ values under the compressive and expanded strains are almost the same, and the difference between them is less than 0.5 eV.

It is worth mentioning that we settle on the full potential rather than pseudopotential method for electronic structure calculations and thus easily obtain the DP constants of core states (1S and 2S) summarized in Table 2. The core state energy like the DV state energy is insensitive to strain types (see Fig.S1 of the Supporting Information), and the differences between them are less than 0.02 eV.

After correcting the DP constants at the VBM and CBM using the core state energy, we have found some abnormal conclusions as follows: (i) the DP constants change from negative to positive values; (ii) the ratios of $\lambda^{X}_{VBM}$ to $\lambda^{X}_{CBM}$ increase significantly; (iii) the $|\lambda^{X}_{CBM}|$ values are much smaller than those of the $|\lambda^{X}_{VBM}|$, which is the opposite situation before the correction, because of the reference energy levels. In addition, different strain types with the same strength lead to the difference that is comparable to, even larger than the corrected DP constants. For instance, the difference between the $\lambda^{S}_{VBM}$ and $\lambda^{L}_{VBM}$ is as large as 2.07 eV significantly larger than the $\lambda^{L[1S]}_{VBM}$ and $\lambda^{L[2S]}_{VBM}$. More importantly, choosing different core state energies as the reference levels makes a difference to the DP constants, leading to a large dissimilarity about 1 eV. Evidently, this means that the strain has great influence on the deep core states.

Surprisedly, the corrected DP constants of the CBM, as well as the other theoretical values, are in agreement well with the experimental data [28] that was determined in the heavily As-doped Si with the n-type carrier density as high as 5$\times 10^{21}cm^{-3}$. Unfortunately, the theoretical works neglect an important experimental parameter of the carrier density related with the Fermi energy level. Based on Boltzmann theory, the Fermi energy level corresponds to the carrier density $n$ that is given by

$$n=n_{0}-\int{f(\varepsilon,\varepsilon_{F},T)D(\varepsilon)d\varepsilon},$$

(30)

where D, $\varepsilon_{F}$ and $\emph{n}_{0}$ are the total density of states, the Fermi energy level and the number of the valence electrons in the unit cell. Fig. 2 shows that the Fermi energy level of silicon at the fixed carrier density is insensitive to temperature, however, is highly sensitive to the carrier density at the fixed temperature. When the carrier density changes from $2.7\times 10^{19}cm^{-3}$ to $5\times 10^{21}cm^{-3}$, the Fermi energy level rises from 5.94 eV (corresponding to CBM) to 6.60 eV. This implies that the experimental DP constants in Ref. [28] is for the Fermi energy level rather than for the band edge where it is theoretically predicted in Refs. [26, 27, 16]. Therefore, the high degree of agreement between theoretical and experimental data strongly indicates that the theoretical predictions are problematic. This is also confirmed by the fact that the calculated mobility with the corrected DP constant severely deviates from the experimental data.

In this work we consider that the strain affects not only core electrons but also valence electrons because the strain affects the entire periodic potential especially near the silicon nucleus. Fig. 3(a) shows that both the strength and shape of the periodic potential have undergone significant changes in the entire region when the volume $V_{0}$ is expanded to 301 $V_{0}$. For instance, the periodic potential strength at the $V_{0}$ is lower than -6.95 eV, however, it rises to roughly 0 eV at the 301 $V_{0}$, and the potential shape change from the curve like the spring to an approximate straight line. In fact, the two periodic potential values are not absolute due to the different zero-potential points. However, the potential difference can actually reflect the change in the periodic potential. The potential difference in Fig. 3(b) shows the shape and strength of the periodic potential have changed even if a small volume change of 1% occurs. Especially, the closer to the silicon nucleus, the greater the periodic potential changes.

The electron density $\rho$(r) is determined by the probability density related to the wave function that is independent of the zero-potential point. The electron density difference truly reflects that the electron density around the Si nucleus has more obvious change (see Fig. 3(c)). This change originates from the atomic interaction change caused by volume variation, and inevitably leads to the change for the core state energy. Generally, the energy of the core electron is determined by the zero-potential point and the atomic interactions composed of the exchange correction interaction and the Coulomb interaction.

In this work we expand the cell volume large enough so that the effect of the atomic interactions on the core electrons can be neglected. Obviously, the change in the core state energy is mainly due to the zero-potential rather than atomic interactions. Therefore, the correction for DP constant is obtained easily.

Considering limited computing resources, we used the standard GGA rather than mBJ-GGA strategy to perform the self-consistent calculation with rigorous charge and energy convergence limits ($10^{-6}$ e and $10^{-6}$ eV). It is worth noting to ensure the $R_{mt}\cdot K_{max}$ value in all calculations is 7.0 by setting suitable dimension parameters of WIEN2k code. Fig. 4 indicates the relationship between the core state energy and the volumetric strains seems to be linear. The slope values for 1S and 2S states are -0.024601 eV and -0.024605 eV at 451 $V_{0}$, and decrease down to -0.021944 eV and -0.021987 eV at 501 $V_{0}$, respectively.

We consider the slight slope results from the zero-potential change induced by the volumetric strain, and assume that the zero potential $E^{0}$ and the 1S core state energy $E^{1S}$ are approximately equal, and have the property of Coulomb potential

$$E^{0}\propto\frac{k}{r}.$$

(31)

Therefore, when the volume change from $V_{N}=NV_{0}$ to (1+$\delta^{V}$) $V_{N}$, the zero potential difference $\Delta E^{0}_{N}$ can be described by

$$\Delta E^{0}_{N}\propto\frac{k}{\sqrt[3]{(1+\delta^{V})V_{N}}}-\frac{k}{\sqrt[3]{V_{N}}},$$

(32)

So the correction for the DP constants at $V_{N}$ is about

$$\lambda_{N}^{V}=\frac{\Delta E^{0}_{N}}{\delta^{V}}\approx\frac{\Delta E^{1S}_{N}}{\delta^{V}}.$$

(33)

And then one can obtain the correction for the DP constant at $V_{0}$ according the equation

$$\lambda_{1}^{V}=\frac{\Delta E^{0}_{1}}{\delta^{V}}\approx\frac{k}{\sqrt[3]{(1+\delta^{V})V_{0}}}-\frac{k}{\sqrt[3]{V_{0}}}\approx\sqrt[3]{N}\lambda^{V}_{N}$$

(34)

Herein, taking N=501, the $\lambda_{1}^{V}$ is only about -0.1743 eV that is insignificant compared to the $\lambda^{V}_{e(h)}$ without the correction. It is worth mentioning that we also calculated the slope values at other large volumes, and found that most of them follow this trend of the slope value decreasing with the volume. However, a few slopes do not follow this trend, implying that the uncertainty of zero-potential, but it is certain that the corresponding $|\lambda_{1}^{V}|$ is not beyond 0.3 eV.

In addition, the effective mass has an important effect on the mobility, but it is difficult to calculate because it is sensitive to the step size during the numerical derivative process. The dependence of the effective mass with step size for is summarised in Table 3. We average the last five values as the effective mass values. The electron effective mass as well as the light-hole mass is in good agreement well with measured and other calculated values. The heavy-hole mass 0.34 $m_{0}$ is obviously larger than the other theoretical values, however, it is much closer to the experimental value 0.43 $m_{0}$ [38]. We adopted the calculated elastic constant $c_{11}$ = 153 GPa that is also in good agreement with experimental and theoretical values [42, 43, 44, 45]. Substituting the relevant parameters into the Eqs. (19) and (27) and using the local DOS effective masses calculated by Eqs. (16) and (15) respectively, the hole mobility at 300 K is calculated to be 704 cm²/Vs and 770 cm²/Vs that reasonably agrees with experimental value of 510 cm²/Vs [40] (See Table 3). Adopting the mobility effective mass calculated by Eq. (29), we obtain 1771 $cm^{2}$/Vs that agrees well with experimental data.

In summary, we improve the mobility formula based on the DP theory, and combine with the first-principle method to calculate the electron and hole mobilities for single crystal silicon. In the improved mobility formula, the mobility effective mass is redefined, more importantly, a method is proposed and applied to calculate the corrections for the DP constant. It is found that the corrections for the DP are slight and thus can be neglected. Our data on the mobility of the single crystal silicon is in good agreement with experimental results, and thus the method can be effectively applied to the prediction of the mobility in bulk materials.

This work is supported by the Scientific and Technological project of Nanchang Institute of Science and Technology (Grant No. NGKJ$-$22$-$09).