First-principles predictions of Hall and drift mobilities in semiconductors

Carrier mobility is at the root of our understanding of electronic devices. We present a unified methodology for the parameter-free calculations of phonon-limited drift and Hall carrier mobilities in real materials within the framework of the Boltzmann transport equation. This approach enables accurate and parameter-free calculations of the intrinsic mobility and will find applications in the design of electronic devices under realistic conditions of strain and temperature. The methodology exploits a novel approach for incorporating the effect of long-range quadrupole fields in the electron-phonon scattering rates and capitalizes on a rigorous and efficient procedure for numerical convergence. The accuracy reached in this work allows us to assess the impact of common approximations employed in transport calculations, including the role of exchange and correlation functionals, spin-orbit coupling, pseudopotentials, Wannier interpolation, Brillouin-zone sampling, dipole and quadrupole corrections, and the relaxation-time approximation. We study diamond, silicon, GaAs, 3C-SiC, AlP, GaP, c-BN, AlAs, AlSb, and SrO, and find that our most accurate calculations predict Hall mobilities significantly larger than the experimental data in the case of SiC, AlAs, and GaP. We identify possible improvements to the theoretical and computational frameworks to reduce this discrepancy, and we argue that new experimental data are needed to perform a meaningful comparison, since almost all existing data are more than two decades old. By setting tight standards for reliability and reproducibility, the present work aims to facilitate validation and verification of data and software towards predictive calculations of transport phenomena in semiconductors.

Figure 1

(a) Relative deviation of the calculated carrier mobility with respect to our most accurate calculations when using one of the approximations described in Sec. 4b. Electron mobilities are in orange, and holes are in green. In each case, only one effect is considered, while all other settings are the same as for the reference calculation. We consider the following effects: (i) Long range part of the dynamical matrix (D) and electron-phonon matrix elements (G). The label “0” means that no long range part has been considered, and “DD,” “DQ,” and “QQ” mean that dipole-dipole, dipole-quadrupole and quadrupole-quadrupole contributions have been included. The labels “eD” and “eQ” mean that monopole-dipole and monopole-quadrupole interactions are included, respectively. (ii) Using the velocity computed in the local approximation via Eq. (26) instead of Eq. (27). (iii) Neglecting SOC. (iv) Employing the SERTA in Eq. (5) instead of Eq. (1). (v) Using the relaxed lattice parameter of the PBE exchange-correlation functional. (vi) Using the experimental lattice parameter with the LDA exchange-correlation functional. (vii) Using the PBE exchange-correlation functional and experimental lattice parameter, but with a different parametrization of the pseudopotentials. The gray boxes represent the first and third quartile, while the average values are shown with a small horizontal black line. The standard deviations are reported at the top of the figure. The electron mobility of GaAs is not reported in this figure since the values are off the chart. (b) Same as in (a), but the deviations are reported as absolute values in ${\mathrm{cm}}^2 {\mathrm{V}}^{-1} {\mathrm{s}}^{-1}$.

Figure 2

Comparison between the deformation potential of the ${\mathrm{\Gamma }}_1$ valence band of silicon from the direct DFPT (dots) and the Wannier interpolation with and without quadrupoles, for different k-point grids (the q-point grid is half of that). The components included in the dynamical matrix and in the electron-phonon matrix elements are labeled as “D” and “G,” respectively. The labels “DD,” “DQ,” and “QQ” indicate the dipole-dipole, dipole-quadrupole and quadrupole-quadrupole contributions, while “eD” and “eQ” denote the monopole-dipole and monopole-quadrupole interactions, respectively.

Figure 3

[(a)–(j)] Comparison between the deformation potential of the valence band manifold from the direct DFPT calculations (dots) and the Wannier interpolation including dipoles and quadrupoles (green lines) for diamond, silicon, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb, and SrO, respectively. The $\mathbf{k}$ point for the deformation potential is located at the valence band maximum.

Figure 4

Convergence of the SERTA (orange) and BTE (red) mobility and the SERTA (light green) and BTE (dark green) Hall factor as a function of the fine k-and q-point grid size, using the converged coarse grids from Fig. 19. [(a)–(e)] Electron and (f)–(j) hole Hall factor and mobility of diamond, silicon, GaAs, 3C-SiC, and AlP. [(k)–(o)] Electron and [(p)–(t)] hole Hall factor and mobility of GaP, cubic BN, AlAs, AlSb, and SrO. Both quantities converge linearly with the fine grid density, and can therefore be extrapolated to infinity. The shaded area denotes the interpolation range for the linear least squares fit of the Hall factor and mobility.

Figure 5

(a) Convergence of the SERTA and BTE electron mobility of c-BN as a function of the fine k-and q-point grid size. The coarse grid employed in these calculations is given in Fig. 19. We test various fixed smearing parameters as well as the adaptative smearing. (b) Electron mobility of c-BN at extrapolated infinite grid density, plotted as a function of smearing parameter. The dashed lines are the extrapolation to zero smearing. The use of adaptative smearing (stars) yields results close to the zero-smearing extrapolation.

Figure 6

Wannier-interpolated band structures [(a)–(e), (k)–(o)] and phonon dispersions [(f)–(j), (p)–(t)] of diamond, silicon, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb, and SrO, respectively. The gray dots are experimental neutron scattering data for diamond [111], silicon [112, 113], GaAs [114], 3C-SiC [115], AlP [116], GaP [117], c-BN [118], AlAs [114], AlSb [119], and SrO [120].

Figure 7

(a) Valence band deformation potential corresponding to c-BN from 0 to 2 eV/bohr. The label “D” refers to the dynamical matrix, the label “G” refers to electron-phonon matrix elements. The labels “DD,” “DQ,” and “QQ” indicate that dipole-dipole, dipole-quadrupole, or quadrupole-quadrupole contributions are included. The labels “eD” and “eQ” indicate monopole-dipole and monopole-quadrupole terms. The black circles are the reference DFPT calculations. (b) Same as in (a), but for the ${\mathrm{\Gamma }}_1$ valence band of 3C-SiC.

Figure 8

Absolute value of the mass-scaled eigendisplacement vectors $e_{\kappa \alpha ,\mathbf{q}\nu }$ in c-BN for the LA mode of the Boron atom in the $x, y$, and $z$ Cartesian direction along the $L-\mathrm{\Gamma }$ q-point line. The label “D” refers to the dynamical matrix while the labels “DD,” “DQ,” and “QQ” indicate that dipole-dipole, dipole-quadrupole, or quadrupole-quadrupole contributions are included. The black arrows are the reference DFPT calculations.

Figure 9

[(a)–(j)] Mode decomposition and total scattering rates at room temperature as a function of energy from the band edges for diamond, silicon, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb, and SrO, respectively. The energy range covers the energy window used in the calculations.

Figure 10

[(a)–(j)] Spectral decomposition of the electron (orange) and hole (green) scattering rates as a function of phonon energy at 300 K for diamond, silicon, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb, and SrO, respectively. The rates are calculated as angular averages for carriers at an energy of $3k_{\mathrm{B}}T/2=39$ meV away from the band edges. The dashed lines represent the cumulative integrals of the calculated rates, and add up to the carrier scattering rate ${\tau }_{3/2k_{\mathrm{B}}T}^{-1}$. The percentage indicates the contributions from acoustic modes.

Figure 11

Comparison between the room temperature first-principles BTE drift mobility (blue) and the mobility obtained with the Drude formula from Eq. (7) (green). The scattering lifetimes ${\tau }_{3/2k_{\mathrm{B}}T}$ are calculated as angular averages for carriers at an energy of $3k_{\mathrm{B}}T/2=39$ meV away from the band edge. The materials are arranged in decreasing order of their BTE drift mobility value.

Figure 12

Temperature dependence of the calculated electron [(a)–(e), (k)–(o)] and hole [(f)–(j), (p)–(t)] Hall factor for diamond, silicon, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb, and SrO, respectively. We show results within the SERTA (orange) and BTE (green). The calculations have been extrapolated to infinitely dense Brillouin zone grids. The gray dashed lines correspond to the isotropic approximation introduced in Eq. (15). Experimental data for diamond are taken from Konorova and Shevchenko [163]. Experimental data for GaAs are taken from Rode [164] and Stillman et al. [165], for a low magnetic field of 0.5 kG.

Figure 13

Calculated temperature dependence of electron [(a)–(e), (k)–(o)] and hole [(f)–(j), (p)–(t)] mobilities, and comparison to experimental data for diamond, silicon, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb, and SrO, respectively. The calculated drift mobilities are in red, the Hall mobilities are in green, and the experimental datapoints are in gray. The dashed lines indicate the power law fits with exponents reported in Table S2 from Ref. [108] and obtained on the 150 to 500 K range. The references for the experimental data are as follows: diamond [139, 166, 167, 168, 169, 170], silicon [171, 172, 173, 174, 175], GaAs [164, 176, 177, 178, 179, 180], 3C-SiC [181, 182, 183, 184], AlP [185, 186], GaP [187, 188, 189], AlAs [190, 191, 192, 193], and AlSb [194, 195, 196].

Figure 14

Comparison between theoretical and experimental (a) temperature exponents and (b) mobilities. The values are reported in Table S2 of Ref. [108]. The mean absolute error (MAE) and mean absolute relative error (MARE) are indicated as inset in the bottom figure and do not include the electron mobility of GaAs.

Figure 15

Visual comparison of the room temperature theoretical drift mobilities using the self-energy relaxation time approximation (SERTA), the Boltzmann transport equation (BTE), and the Hall mobilities using the BTE. Materials are sorted by mobility decreasing towards the right hand side using the Hall BTE as reference. The shaded region indicates the range of measured electron (-e) and hole (-h) mobilities using the values reported in Table S2 of Ref. [108].

Figure 16

Comparison between theoretical and experimental mobility with experiment. Black triangles indicate the impact of rescaling the matrix elements while blue crosses indicate additional effective masses rescaling. The mean absolute error (MAE) and mean absolute relative error (MARE) are indicated as inset in the bottom figure and do not include the electron mobility of GaAs.

Figure 17

Conduction band edge of silicon—including spin-orbit coupling (SOC)—where the zero energy has been set at the $X$ high-symmetry point and computed using DFT (black line). (a) Electronic bands obtained using Wannier interpolation with $\mathbf{k}$-point grids ranging from ${10}^3$ to ${16}^3$ and constructed from an initial set of four $s{p}^3$ projections on each Si atom. Both the valence and conduction bands have been wannierized, totalling 16 Wannier functions due to SOC. (b) Electronic bands obtained using Wannier interpolation with $\mathbf{k}$-point grids ranging from ${10}^3$ to ${16}^3$ and constructed from an initial set of one $s$ and five $d$ projections on one Si atom only. Only the conduction bands have been wannierized, totalling 12 Wannier functions due to SOC.

Figure 18

The first two rows present the Wannier functions in real space for the conduction and valence bands, respectively; generated using the xcrysden visualization program [219]. In the case of conduction bands, we show the $\pm 1.0{\AA }^{-3/2}$ isosurface everywhere, except in the case of SrO for which we show the isosurface at $\pm 2.0{\AA }^{-3/2}$. In the case of valence bands, we show the isosurfaces at $\pm 0.2, \pm 0.5{\AA }^{-3/2}, \pm 1.0$, and $\pm 1.0{\AA }^{-3/2}$ for Si, diamond, 3C-SiC and GaAs, respectively; and the isosurface at $\pm 2.0{\AA }^{-3/2}$ for all the other materials. By denoting the spinor components as $[\varphi ,\psi ]$, in the figure we show $\left|\varphi \right|\times \mathrm{sign}\left(\Re \right\{\varphi \left\}\right)$. The next five rows show the spatial decays of the electronic Hamiltonian, velocity matrix, dynamical matrix and electron-phonon vertex in the Wannier representation, in the limiting case of $R_p=0$ or $R_e=0$, as a function of $R=|R_e-R_p|$. The data points correspond to the largest value taken over the Wannier functions in each of the two unit cells at a distance $R$, and are normalized such that the value is 1 at $R=0$. The pink and purple colors indicate the wannierization for the conduction and valence manifold, respectively. The size of the k- and q-point grids used for each manifold are also indicated in pink and purple colors, respectively.

Figure 19

Convergence of BTE electron [(a)–(e), (k)–(o)] and hole [(f)–(j), (p)–(t)] mobility as a function of the coarse k-point grid size (magenta dots) for diamond, silicon, GaAs, 3C-SiC, AlP, GaP, cubic BN, AlAs, AlSb, and SrO, respectively. The coarse q-point grid and the fine k and q-point grids are given in the inset. The vertical gray line indicates that convergence has been achieved, and represents the coarse grid employed in the remainder of this work. The orange dots represent the BTE mobility after the rescaling defined by Eq. (8).

Figure 20

Convergence of the CARTA electron [(a)–(e), (k)–(o)] and hole [(f)–(j), (p)–(t)] mobility as a function of the coarse k-point grid size (magenta dots). We only show the relative convergence since the CARTA mobility contains the arbitrary parameter $g^2$, see Eq. (8). The coarse q-point grid and the fine k and q-point grids are given in the inset. The vertical gray line indicates that convergence has been achieved, and represents the coarse grid that is used in the remainder of this work. The orange dots are the mobility effective masses, and are referred to the vertical axis on the right. The horizontal orange and magenta dashed lines are the DFT reference values.

Figure 21

Comparison between the coarse-grid convergence in the present work and in Ref. [16], for the electron mobility of (a) Si and (b) GaP as a function of the coarse q-point grid. All calculations are performed within SERTA to be consistent with Ref. [16]. The coarse $\mathbf{k}$-point grid used in Ref. [16] is fixed to ${18}^3$, while we used a coarse $\mathbf{k}$-point grid that is twice the one reported in this figure for the q-point grid. We used equal fine k- and q-point grids of ${150}^3$ for Si and ${120}^3$ for GaP, while Ref. [16] used a fine ${72}^3\mathbf{k}$-point and a ${144}^3\mathbf{q}$-point grid for Si and a ${78}^3\mathbf{k}$-point and a ${156}^3\mathbf{q}$-point grid for GaP.

Figure 22

Simplified schematic (two-dimensional $2\times 2$ k/q grid) of the electron-phonon matrix elements where ${\mathbf{r}}_n$ and ${\mathbf{r}}_m$ indicate the position of the electronic Wannier centres, and ${\mathit{\tau }}_{\kappa }$ is the position of the atom $\kappa$ within a primitive cell. The square lattice indicates the unit cells of the crystal and ${\mathbf{R}}_p$ and ${\mathbf{R}}_{p^{\prime \prime }}$ are the lattice vectors for the electronic and phononic grids, respectively. When performing the electron or phonon Fourier transform, we use the direct lattice point $p$ or $p^{\prime \prime }$, respectively, if the trial point $\mathbf{t}={\mathbf{R}}_p+{\mathbf{r}}_m-r_n$ or ${\mathbf{t}}^{\prime \prime }={\mathbf{R}}_{p^{\prime \prime }}+{\mathit{\tau }}_{\kappa }-{\mathbf{r}}_n$ is closer to $\mathbf{0}$ than any other supercell center $\mathbf{S}$. The red square denotes the central supercell.

Figure 23

Acoustic phonon dispersions of 3C-SiC near the zone center, interpolated from a coarse $4\times 4\times 4\mathbf{k}$-point and $\mathbf{q}$-point grid using a $\mathrm{\Gamma }$-centered fixed Wigner-Seitz cell with dipole-dipole correction (dash gray line), an optimal Wigner-Seitz cell with dipole-dipole correction (red line), and an optimal Wigner-Seitz cell with dipole-dipole, dipole-quadrupole and quadrupole-quadrupole corrections (blue line). The black dots are reference datapoints obtained via direct DFPT calculations.

Figure 24

Deformation potential of 3C-SiC without spin-orbit coupling of the first electronic band ($m=n=1$) on a coarse $3\times 3\times 3\mathbf{k}$-point and $\mathbf{q}$-point grid including dipole and quadrupole corrections. The Wigner-Seitz cell is either $\mathrm{\Gamma }$-centered or constructed using ${\mathbf{R}}_p+{\mathbf{r}}_m-{\mathbf{r}}_n\in (i_p\times j_p\times o_p)$ for the electrons and ${\mathbf{R}}_{p^{\prime \prime }}+{\mathit{\tau }}_{\kappa }-{\mathbf{r}}_n\in (i_{p^{\prime \prime }}\times j_{p^{\prime \prime }}\times o_{p^{\prime \prime }})$ for the phonons, respectively. The black dots are reference datapoints obtained via direct DFPT calculations.