Transport and noise of hot electrons in GaAs using a semi-analytical model of two-phonon polar optical phonon scattering

We first review the ab-initio treatment of electron transport and electronic noise using the BTE Hamaguchi and Hamaguchi (2010). For a spatially homogeneous system with electric field $\vec{\mathcal{E}}$, the electron distribution function $f_{m\boldsymbol{k}}$ is governed by

where $f_{m\boldsymbol{k}}$ is the electron occupation in the state with band index $m$ and wave vector $\boldsymbol{k}$, $e$ is the fundamental charge, and $\mathcal{I}[f_{m\boldsymbol{k}}]$ is the collision term, which describes the scattering of electrons with phonons Bernardi (2016). In this work, index $m$ is dropped in all the following derivations for simplicity as only one band is relevant for electron transport in GaAs in the range of electric fields considered.

For non-degenerate electrons, the collision term can be linearized as Cheng et al. (2022):

where $\Theta_{\boldsymbol{k}^{\prime},\boldsymbol{k}}$ is e-ph collision matrix. The diagonal elements $\Theta_{\boldsymbol{k},\boldsymbol{k}}$ are equal to the total scattering rates as $\Theta_{\boldsymbol{k},\boldsymbol{k}}=\Gamma_{\boldsymbol{k}}=-\sum_{\boldsymbol{k}^{\prime}\neq\boldsymbol{k}}\Theta_{\boldsymbol{k}^{\prime},\boldsymbol{k}}$. With this linearization and a finite difference representation of the derivative operator $\nabla_{\boldsymbol{k}}$ Choi et al. (2021), Eq. 1 becomes a linear partial-differential equation which can be solved by numerical linear algebra. The equation for the steady distribution function $f_{\boldsymbol{k}}^{s}$ is given by Eq. 6 of Ref. Choi et al., 2021. Steady-state mean transport properties such as drift velocity can be calculated with the appropriate Brillouin zone sum using this distribution.

The current power spectral density (PSD) is used to characterize fluctuations in occupation about the mean distribution. The PSD is defined as the Fourier transform of the autocorrelation of the current density fluctuations (Eq. 12 of Ref. Choi et al., 2021). Following Ref. Choi et al., 2021, the current PSD at frequency $\omega$ can be calculated as

where $v_{\boldsymbol{k}\alpha}$ is the group velocity of the electron at state $\boldsymbol{k}$ along axis $\alpha$, $\mathcal{V}_{0}$ is the supercell volume, and $g_{\boldsymbol{k}\beta}$ is the effective distribution function Rustagi and Stanton (2014) which satisfies the following equation:

where $D$ is the finite difference representation of $\nabla_{\boldsymbol{k}}$ specified in Eq. 1, $f_{\boldsymbol{k}}^{s}$ is the steady-state occupation for the state at wave vector $\boldsymbol{k}$, and $V_{\alpha}$ is the drift velocity along axis $\alpha$.

The e-ph collision matrix is obtained from perturbation theory in orders of the e-ph interaction strength. For 1ph scattering of non-degenerate electrons, the non-diagonal scattering matrix elements are given by Cheng et al. (2022):

for $\boldsymbol{k}\neq\boldsymbol{k}^{\prime}$, where $g_{\nu}(\boldsymbol{k},\boldsymbol{q})$ is the e-ph scattering matrix element, $\epsilon_{\boldsymbol{k}}$ is the energy of the electronic state $\boldsymbol{k}$, $\omega_{\nu\boldsymbol{q}}$ is the frequency of phonon with mode $\nu$ and wave vector $\boldsymbol{q}$, and $N$ is the total number of $\boldsymbol{k}$ points in the Brillouin zone, $N_{\boldsymbol{q}}$ is the phonon occupation according to the Bose-Einstein statistics. The two delta functions are energy conservation conditions for the emission and absorption subprocesses, respectively.

The collision integral for 2ph scattering was derived in Ref. Lee et al., 2020, and the linearized form is given in Ref. Cheng et al., 2022. 2ph scattering can be divided into the two-emission (2e), one-emission-one-absorption (1e1a) and two-absorption (2a) subprocesses. Here, we rewrite the original formalism in Ref. Lee et al., 2020 to facilitate the derivations in the next section. Assuming non-degenerate electron statistics in Eq. 4 in Ref. Lee et al., 2020, we rewrite Eq. 12 in Ref. Cheng et al., 2022 by splitting $W^{(i)}$ to two parts:

where

where $\eta$ is a positive infinitesimal and $\Gamma_{\boldsymbol{k}+\boldsymbol{q}}=\Gamma^{(\mathrm{1ph})}_{\boldsymbol{k}+\boldsymbol{q}}+\Gamma^{(\mathrm{2ph})}_{\boldsymbol{k}+\boldsymbol{q}}$ is the total scattering rate of the intermediate state $\boldsymbol{k}+\boldsymbol{q}$. For the 2e and 2a subprocesses, the two terms give the same contribution after the summation. For the 1e1a subprocess, the two terms in Eq. 6 physically represent the emission-then-absorption (a-e) and absorption-then-emission (e-a) subprocesses, respectively. Finally, we exchange the summation order of $\boldsymbol{q}$ and $\boldsymbol{p}$ in Eq. 9 in Ref. Cheng et al., 2022 for the second term of Eq. 6 and arrange the equations to obtain:

for $\boldsymbol{k}\neq\boldsymbol{k}^{\prime}$, where $\alpha_{1}$ and $\alpha_{2}$ indicate whether the first and second phonon is emitted ($\alpha_{1,2}=1$) or absorbed ($\alpha_{1,2}=-1$), so that the four combinations of $\alpha_{1,2}=\pm 1$ describe the four subprocesses.

The term $\tilde{\Theta}^{(\alpha_{1},\alpha_{2})}_{\boldsymbol{k},\boldsymbol{q}\nu,\boldsymbol{p}\mu}$ in Eq. 8 is defined as

As the 2ph scattering rates depend on the intermediate state rates, the 2ph scattering must be calculated iteratively. In each iteration, the intermediate state rates of the previous iteration is used to update the 2ph scattering matrix $\Theta^{(\mathrm{2ph})}_{\boldsymbol{k}^{\prime},\boldsymbol{k}}$ and scattering rates $\Gamma^{(\mathrm{2ph})}_{\boldsymbol{k}}$.

Compared with 1ph scattering, 2ph scattering is much more computationally expensive, particularly for high-field transport which requires a larger energy window than for low-field transport. In Ref. Cheng et al., 2022, several approximations were made to make the computation feasible, including limiting the number of self-consistent iterations of the 2ph rates to three, restricting the off-shell extent $|\epsilon_{\boldsymbol{k}^{\prime}}-\epsilon_{\boldsymbol{k}+\boldsymbol{q}}+\alpha\hbar\omega_{\nu\boldsymbol{p}}|$ to 25 meV, and neglecting the interference term $\text{Re}(\tilde{W}_{\boldsymbol{k},\boldsymbol{q},\boldsymbol{p},\alpha_{2}}\tilde{W}_{\boldsymbol{k},\boldsymbol{p},\boldsymbol{q},\alpha_{1}}^{*})$ term in Eq. 9. Additionally, the maximum grid density that could be used was $200\times 200\times 200$. The effect of these approximations on the observable transport and noise properties was not assessed. In particular, the on-shell approximation neglects off-shell processes and thus underestimates the scattering rates. These approximations were mentioned as possible reasons for the PSD discrepancy in Ref. Cheng et al. (2022).

In this section, we introduce a semi-analytical model to treat 1ph and 2ph $\Gamma-\Gamma$ intravalley scattering by the Fröhlich interaction that retains the accuracy of the ab-initio formalism to within a few percent while allowing the approximations described above to be lifted. This model is based on the fact that over the range of wavevectors considered in this study, the $\Gamma$ valley in GaAs is nearly spherically symmetric, and $\Gamma\text{-}\Gamma$ scattering can be accurately described by using only the Fröhlich interaction Verdi and Giustino (2015). The model is valid only for $\Gamma$ intravalley scattering because $\Gamma$-L intervalley scattering lacks an analytic description of similar accuracy.

The semi-analytical model uses the following approximations. First, the band structure is described using the Kane model Kane (1957) for a spherically symmetric, non-parabolic band. This description is accurately satisfied for the $\Gamma$ valley, with the Kane model bands deviating from the ab-initio band structure by at most 7% over the range of wave vectors considered ($\sim 0.1G$, where $G$ is the reciprocal lattice constant). Second, prior works have shown that $\Gamma\text{-}\Gamma$ e-ph scattering in GaAs is dominated by longitudinal optical (LO) phonons via the Fröhlich interaction Fröhlich (1954). We therefore neglect scattering processes involving other phonon branches and scattering mechanisms. The computed matrix elements $g_{\text{LO}}(\boldsymbol{k},\boldsymbol{q})$ for Fröhlich scattering are found to exhibit negligible anisotropy so that $g_{\text{LO}}(\boldsymbol{k},\boldsymbol{q})=g_{\text{LO}}(k)$, enabling an analytic form of $g_{\text{LO}}(k)$ to be fit to the ab initio values as described in Sec.III. In the range of wave vectors considered, this approximation is satisfied to within 3% Verdi and Giustino (2015). Third, we take the LO phonon frequency to be a constant $\omega_{\text{LO}}=35$ meV. In the range of phonon wave vectors $k\in(0,0.2G)$ considered here, this assumption is satisfied to within less than 0.3%.

We now discuss the treatment of 1ph and 2ph e-ph scattering based on these simplifications. The summation in Eq. 8 may be rewritten as an integral in the Brillouin zone over the intermediate wave vector $k_{m}$ by letting $\boldsymbol{q}\rightarrow k_{m}-\boldsymbol{k}$ and $\boldsymbol{p}\rightarrow\boldsymbol{k}^{\prime}-k_{m}$. Additionally, we exploit spherical symmetry to rewrite all the quantities in spherical coordinates as: $\Theta_{\boldsymbol{k}^{\prime},\boldsymbol{k}}=\Theta(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}})$ and $\Gamma_{\boldsymbol{k}}=\Gamma_{k}$. After some simplifications, we obtain the 1ph and 2ph collision matrices as:

and

where $\Omega_{\text{BZ}}$ is the Brillouin zone volume. $A_{\alpha}$ is the phonon occupation factor defined as

where $N_{\text{LO}}=(\exp(\hbar\omega_{\text{LO}}/k_{B}T)-1)^{-1}$ is the LO phonon occupation, and $I^{(\alpha_{1},\alpha_{2})}=I_{\text{1}}^{(\alpha_{1},\alpha_{2})}+I_{\text{2}}^{(\alpha_{1},\alpha_{2})}$ is decomposed to the non-interference part $I_{\text{1}}^{(\alpha_{1},\alpha_{2})}$ and the interference part $I_{\text{2}}^{(\alpha_{1},\alpha_{2})}$:

and

Equation 13 can be further simplified by writing the integration in spherical coordinates and separating the radius and angular part:

where $\tilde{I}_{\text{1}}^{(\alpha_{1},\alpha_{2})}$ is the angular part defined as:

where $\theta_{k_{m}}$, $\phi_{k_{m}}$ are the polar angle and azimuthal angle defining the intermediate wave vector $k_{m}$, respectively. Since $\tilde{I}_{\text{1}}^{(\alpha_{1},\alpha_{2})}$ is independent of the band structure and the self-energy, recomputation of this term in each 2ph iteration is not required. In practice, to significantly reduce the computational cost, $\tilde{I}_{\text{1}}^{(\alpha_{1},\alpha_{2})}(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}},k_{m})$ are precomputed on a grid of $k$, $k^{\prime}$, $\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$ and $k_{m}$ before the 2ph iteration. Here $k$, $k^{\prime}$, $\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$ are not independent of each other due to the energy conservation condition. Once $\tilde{I}_{\text{1}}^{(\alpha_{1},\alpha_{2})}(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}},k_{m})$ is computed on a predefined grid, $\Theta^{(\mathrm{2ph})}(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}})$ can be calculated according to Eqs. 11, 15 and 14. We note that such separation of the radius and spherical part is not valid for $\tilde{I}_{\text{2}}^{(\alpha_{1},\alpha_{2})}$, so an expensive iterative update is required. However, since $\tilde{I}_{\text{2}}^{(\alpha_{1},\alpha_{2})}$ is generally much smaller compared with $\tilde{I}_{\text{2}}^{(\alpha_{1},\alpha_{2})}$, we update $\tilde{I}_{\text{2}}^{(\alpha_{1},\alpha_{2})}$ every 10 iterations to decrease the cost of the self-consistent calculations.

To complete the 2ph iteration, the last quantities to be computed are the total 1ph and 2ph scattering rates:

where $\text{type}=\mathrm{1ph},\ \mathrm{2ph}$ indicates the type of scattering. We also perform the radius integration over $k^{\prime}$ analytically to integrate the delta functions in Eq. 10 and  11. The angular integrations in Eq. 17 are performed numerically; details are provided in Sec. III.

The computational flow of the semi-analytical model is as follows. First, we generate a grid of $k$, $\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$, $k_{m}$ and calculate the corresponding $k^{\prime}$ from the energy conservation conditions for each subprocess. Second, we calculate $\Theta^{\mathrm{(1ph)}}(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}})$ by Eq. 10 and $\Gamma_{k}^{(\mathrm{1ph})}$ by Eq. 17. Then, we calculate $\tilde{I}_{\text{1}}^{(\alpha_{1},\alpha_{2})}(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}},k_{m})$ by Eq. 16. Finally, we perform the self-consistent 2ph iterations through Eq. 15, 14, 11 and 17 until convergence, where Eq. 14 is calculated every 10 iterations.

From the perspective of computational cost, the semi-analytical model reduces the number of integration variables in the 2ph scattering rate calculation from 9 (integration over $\boldsymbol{k}$, $\boldsymbol{k}^{\prime}$, $k_{m}$) to 5 ($k$, $\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$, $k_{m}$) due to the spherical symmetry, and avoids the recomputation of $I_{\text{1}}^{(\alpha_{1},\alpha_{2})}$ in the 2ph self-consistent iterations due to the separation of radius and angular integration in Eq. 15. This reduction allows for the use of a denser grid for the intermediate state integration and thereby reduces the discretization error. Therefore, the total scattering rates can be calculated with negligible discretization error compared with the ab-initio calculation. However, since the semi-analytical model is only for $\Gamma$-$\Gamma$ scattering, a discretized scattering matrix ($\Theta_{\boldsymbol{k}^{\prime},\boldsymbol{k}}$) is still needed to compute the drift velocity and current PSD which are affected by $\Gamma$-L intervalley scattering. Therefore, the discretization error in the final state integration cannot be avoided for the present calculations. Nevertheless, the semi-analytical model still decreases the discretization error of the intermediate state integration and treats the full 2ph scattering term. The differences between the ab-initio calculation and the semi-analytical model for $\Gamma$ intravalley scattering are summarized in Table 1.

The ab-initio calculation parameters are identical to those in our previous work Cheng et al. (2022). In brief, electronic structure and e-ph matrix elements are computed using Density Functional Theory (DFT) and Density Functional Perturbation Theory (DFPT) with Quantum Espresso Giannozzi et al. (2009, 2017) with an $8\times 8\times 8$ coarse grid, 72 Ryd plane wave energy cutoff, a relaxed lattice parameter of 5.556 Å, and a non-degenerate carrier concentration of $10^{15}$ cm^-3. Following our previous work Cheng et al. (2022), we apply a band structure correction for both the $\Gamma$ valley and the L valley. For the $\Gamma$ valley, we use a spherically symmetric Kane model band structure Conwell and Vassell (1968) with an experimental effective mass of 0.067$m_{e}$ and a non-parabolicity of 0.64 Ryd^-1 Lundstrom (2000). For the L valley, we shift the energy by 50 meV to achieve the experimental $\Gamma$-L valley separation of 300 meV. Wannier interpolation in Perturbo Zhou et al. (2021) is then applied to interpolate the e-ph matrix elements to a finer grid of $200\times 200\times 200$. After the e-ph matrix elements are obtained, the 1ph and 2ph scattering matrices are computed according to Eqs. 5 and 8. The delta functions in Eqs. 5 and 9 are approximated by a Gaussian function with a standard deviation of 5 meV.

During the computation of scattering matrices, a phonon frequency cutoff of 2 meV is applied to neglect phonons with low frequencies. An energy cutoff of 360 meV above the conduction band minimum is used to reduce the number of k points in the Brillouin zone integration. The 2ph calculation applies the on-shell approximation by restricting the off-shell extent to 25 meV. Following Ref. Choi et al. (2021), the linear systems of equations representing the Boltzmann equation are solved by the generalized minimal residual method (GMRES).

The band structure used in the semi-analytical model is the same as that in the ab-initio calculation. The LO phonon energy is taken to be $\omega_{\text{LO}}=35$ meV. The function $g_{\text{LO}}(k)$ is obtained by a weighted averaged of $g_{\text{LO}}(\boldsymbol{k}_{i},\boldsymbol{p}_{i})$:

where the standard deviation $\sigma=5\times 10^{-4}$ Ryd, the summation is over all the on-shell processes, and the $g_{\text{LO}}(\boldsymbol{k}_{i},\boldsymbol{p}_{i})$ are calculated by the Wannier interpolation. The spherical coordinates integration in Eq. 16 is defined such that the the $\theta_{\boldsymbol{k}_{m}}=0$ direction is orthogonal with both $\boldsymbol{k}$ and $\boldsymbol{k}^{\prime}$. 200 grid points are used for each $\theta_{\boldsymbol{k}_{m}}$ and $\phi_{\boldsymbol{k}_{m}}$ integration. The radius integration in Eq. 15 is transformed into the integration of $\epsilon_{k_{m}}$ and performed using an adaptive integration range with 120 grid points. The integration range is $(\epsilon_{k^{\prime}}+\alpha_{2}\hbar\omega_{0}-6\hbar\Gamma_{k_{m}},\epsilon_{k^{\prime}}+\alpha_{2}\hbar\omega_{0}+6\hbar\Gamma_{k_{m}})$, corresponding to a width of $12\hbar\Gamma_{k_{m}}$ around the center of the Lorentzian function in the denominator of Eq. 15. The relative residual error from this choice of integration limits is estimated as $1/12^{2}\approx 0.7\%$. In the final state integration of both 1ph and 2ph, the angular integrations in Eq. 17 are performed with 200 grid points. All the above numerical integrations are performed on uniform grids using the midpoint rule. The 2ph calculation is performed with 20 iterations, with the interference term Eq. 14 updated every 10 iterations. The relative difference between the 10th and 20th iterations is less than 1%, indicating convergence of the iterative process. The discretization of $\Theta(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}})$ to $\Theta_{\boldsymbol{k}^{\prime},\boldsymbol{k}}$ is performed by the regular grid interpolation provided in scipy Virtanen et al. (2020).

We first present the 1ph and 2ph $\Gamma$ valley scattering rates versus energy obtained by ab-initio calculation and the semi-empirical model for GaAs at 300K in Fig. 1. For both 1ph and 2ph rates, the ab-initio calculations and the semi-analytical model are in quantitative agreement. Specifically, we observe a rapid increase of the 1ph and 2ph rates at $\hbar\omega_{\text{LO}}\approx 35$ meV followed by a nearly constant trend. The degree of agreement between the semi-analytical model and the mean values of the ab-initio calculation is notable considering the semi-analytical model includes off-shell processes and the interference term, both of which are neglected in the ab-initio calculations. The agreement can be attributed in part to the cancellation of errors between the limitation on the iteration number and the on-shell approximation in the ab-initio calculation. The third iteration of the 2ph rates yields values that are overestimated from the converged value by about 9%, while the on-shell approximation and the non-interference approximation are found to underestimate the rates by around 3% and 5%, respectively. These approximations offset each other to yield good agreement between the two approaches. Overall, the contribution of off-shell processes is found to make only a minor contribution to the $\Gamma$ intravalley 2ph scattering rates even up to energies of 300 meV.

The major difference between the ab-initio and the semi-analytical rates is the variation of the individual rates in the ab-initio calculation in a given energy range, which is due to the relatively low grid density used in the ab-initio calculation. As explained in Sec. II.3, the semi-analytical model uses a significantly finer grid, leading to negligible variations in individual scattering rates in the same energy range. Although anisotropy could in principle lead to similar variations of the ab-initio rates, this contribution is negligible owing to the high spherical symmetry of the band structure and e-ph matrix elements (about 3% as mentioned in Sec. II.3).

The high grid density in the semi-analytical model enables features in the scattering rates to be observed that cannot be discerned in the ab-initio calculations, including the previously mentioned rapid increase of 1ph and 2ph rates at $\hbar\omega_{\text{LO}}$, and also a small but evident kink at $2\hbar\omega_{\text{LO}}$ (about 70 meV).

We now analyze each of these observations. For the 1ph rates, the increase at $\hbar\omega_{\text{LO}}$ is because LO phonon emission from an electron may only occur above an energy of $\hbar\omega_{\text{LO}}$. For the 2ph rates, the situation is more complicated due to the existence of four subprocesses (2e, e-a, a-e, 2a) in 2ph scattering. To better understand the features in the 2ph scattering rates, the scattering rates of the four subprocesses are plotted separately in Fig. 1. We observe that the increase of the total 2ph rates at $\hbar\omega_{\text{LO}}$ can be attributed to the e-a subprocess due to a similar reason with the emission subprocess in 1ph process, namely that the e-a subprocess requires the electron to have energy exceeding $\hbar\omega_{\text{LO}}$. The kink at $2\hbar\omega_{\text{LO}}$ comes from the cancellation between the increase of the 2e rates and the decrease of e-a rates. The increase of the 2e rates is due to the emission of 2 LO phonons at energies higher than $2\hbar\omega_{\text{LO}}$. The decrease of the e-a rates is due to the increase of the intermediate state rates in the denominator of Eq. 7 (or Eq. 15). Specifically, for an e-a subprocess with $\epsilon_{\boldsymbol{k}}=2\hbar\omega_{\text{LO}}$, the corresponding intermediate state has the band energy $\epsilon_{\boldsymbol{k}+\boldsymbol{q}}=\hbar\omega_{\text{LO}}$, where an increase of both $\Gamma^{(\mathrm{1ph})}$ and $\Gamma^{(\mathrm{2ph})}$ occurs as explained above. In fact, the effect of the intermediate state self-energy leads to change at any integer multiple of $\hbar\omega_{\text{LO}}$, but the magnitudes are decaying with increasing energy such that they cannot be observed at high energy.

Another observation in Fig. 1 is the difference between the e-a rates and a-e rates asymptotically decreases to zero with increasing energy. This can be understood by analysis of Eq.  11. Specifically, in the process $\boldsymbol{k}_{i}\rightarrow\boldsymbol{k}_{m}\rightarrow\boldsymbol{k}_{f}$, the relative differences between norms of the state vectors $k_{i}$, $k_{m}$, $k_{f}$ become small at high energies (above a few multiples the LO phonon energy), so that the factor $I^{(\alpha_{1},\alpha_{2})}(k,k^{\prime},\theta_{\boldsymbol{k},\boldsymbol{k}^{\prime}})$ becomes insensitive to the subprocess type. Thus, their magnitudes are fully determined by the phonon occupation factor $A_{\alpha_{1}}A_{\alpha_{2}}$ defined in Eq. 12 which satisfy $\Gamma^{(2a)}/A_{-1}^{2}=\Gamma^{(e-a)}/A_{-1}A_{+1}=\Gamma^{(a-e)}/A_{-1}A_{+1}=\Gamma^{(2e)}/A_{+1}^{2}$. From this relationship, we find $\Gamma^{(e-a)}=\Gamma^{(a-e)}$. Furthermore, a common ratio $\Gamma^{(2a)}/\Gamma^{(e-a)}=\Gamma^{(e-a)}/\Gamma^{(2e)}=A_{-1}/A_{+1}$ can also be obtained for the subprocess rates at the high energy region. This relationship is observed in Fig. 1.

We now examine the transport and noise properties from each model. For the semi-analytical model results, the $\Gamma\text{-}\Gamma$ block of the scattering matrix is calculated by the semi-analytical model in Eq. 10 and 11, while the $\Gamma$-L and L-L blocks are those of the ab-initio calculation. Figures 2 and 2 display the drift velocity and normalized current PSD, respectively, from the ab-initio calculations and the semi-analytical model. The experimental measurements are also plotted for comparison. In Fig. 2, the ab-initio calculation and the semi-analytical model give a similar prediction of the drift velocity versus electric field up to 5 kVcm^-1. Both give the low-field mobility of around 7000 cm²V^-1s^-1, which agrees with the experimental value at about 8000 cm²V^-1s^-1 to within around 15%. The similarity between the ab-initio and semi-analytical results is expected due to the agreement of their scattering rates as discussed in Sec. IV.1.

In Fig. 2, the PSD obtained from different experimental measurements reveal a non-monotonic pattern characterized by an initial decrease, followed by a marked rise around the commencement of negative differential mobility, and a subsequent decrease. The origin of this trend was explained in Ref. Cheng et al. (2022). However, both the ab-initio calculation and the semi-analytical model predict the PSD to be nearly independent of electric field and thus fail to predict the characteristic PSD peak at about 3 kVcm^-1.

The lifetimes of photoexcited hot electrons in GaAs have been experimentally measured at 10 K from an analysis of the linewidths of peaks from continuous-wave luminescence spectroscopy. Fasol et al. (1990) In this section, we compare the ab-initio, semi-analytical and experimental scattering rates at cryogenic temperatures. In Fig. 3 and 3, we show the $\Gamma$ and $\Gamma$-L scattering rates obtained by ab-initio calculation, the semi-analytical model, and experiment. The calculations were performed at cryogenic temperatures to enable comparison with experiment. The experimental scattering rates and error bars are converted from the corresponding lifetimes and error bars directly reported in Ref. Fasol et al. (1990). Since the semi-analytical model is only valid for $\Gamma\text{-}\Gamma$ scattering, only the experimental and ab-initio results are shown in Fig. 3.

Figure 3 shows that the experimental and theoretical values for the $\Gamma$-valley scattering all yield a nearly constant value between 100 meV and 325 meV. The experimental rates are about 8 ps^-1 in this energy range, while the semi-analytical 1ph and (1+2)ph calculations give about 3.5 ps^-1 and 6 ps^-1, respectively. This result affirms that the 2ph scattering makes a non-negligible contribution to electron scattering in GaAs. Similarly with Sec. IV.1, the ab-initio calculations give mean values of the scattering rates that agree quantitatively with the semi-analytical model but with substantial scatter about the mean. Such observation again suggests that the approximations in the ab-initio calculations do not result in qualitative deviations.

Figure 3 shows the $\Gamma$-L intervalley scattering rates of the ab-initio calculations and the photoluminescence experiments. According to Ref. Fasol et al., 1990, the experimental $\Gamma$-L rates are obtained by

where $\Gamma_{\Gamma-\Gamma}$ is taken as a constant estimated by fitting the data in Fig. 3. To make a consistent comparison with the experimental estimation, the ab-initio $\Gamma$-L intervalley scattering rates in Fig. 3 are also calculated by Eq. 19 instead of being directly calculated from the scattering matrix. It is found that the experimental intervalley rates is about 6 ps^-1 in the energy range from 340 meV and 400 meV, while the ab-initio 1ph and (1+2)ph give around 2 ps^-1 and 3 ps^-1, respectively. Although the additional 2ph calculation decreases the deviation from experiment results, an underestimation of a factor of two is still observed. This discrepancy could be attributed to the need for an off-shell extent larger than 25 meV in the 2ph calculation, owing to the larger intermediate state scattering rates (appearing in the denominator of Eq. 7 and Eq. 15) at energies above 300 meV. However, at present a larger off-shell extent is computationally infeasible.

A comparison of our computed cryogenic rates with those measured from photoluminescence experiments indicates that the $\Gamma$ rates agree to within 25%, but the $\Gamma$-L intervalley rates are underestimated by around a factor of two. Despite this underestimate, prior work suggests that this effect is unlikely to reconcile the PSD discrepancy. Specifically, Monte Carlo simulations of electron transport in GaAs with a three-valley $\Gamma$-L-X model have found that increased intervalley scattering suppresses the PSD feature (see Fig. 7 of Ref. Požela and Reklaitis (1980)). Therefore, although the possibility cannot be definitively excluded at present, including intervalley scattering processes beyond those treated already is not expected to yield improved agreement with the PSD.

According to Ref. Cheng et al., 2022, another possibility is that the contribution of electron-two-phonon (e-2ph) interaction Kocevar (1980) is not considered. Here, we make a qualitative estimation of the magnitude of this effect based on the Fröhlich mechanism for electron scattering. According to Ref. Kocevar, 1980, the e-ph Hamiltonian up to the second order can be written as

where $\kappa$ is the index of atom in a unit cell, $\boldsymbol{R}$ is the unit cell position, $\boldsymbol{u}_{\boldsymbol{R}\kappa}$ is the corresponding phonon-induced displacement, and $V$ is the electron potential. In the long wavelength limit, the electric potential for the Fröhlich interaction can be obtained by assigning a point dipole to each atom Verdi and Giustino (2015). For acoustic phonons with the same displacements for atoms in the same unit cell, the net dipole moment will be zero and no scattering will occur. Similarly, the electric potential for e-2ph interaction can be obtained by assigning a point quadrupole to each atom. Following the same logic, the net quadrupole moment will be zero if the quadrupoles are induced by two acoustic phonons or two optical phonons, which means that the simultaneous e-2ph interaction based on the Fröhlich interaction can only be induced by an acoustic and optical phonon.

We estimate the order of magnitude of such simultaneous e-2ph interaction involving an acoustic and optical phonon. A full derivation can be found in the Appendix. The final estimated scattering rates in the $\Gamma$ valley is

where $\omega_{\text{O}}$ and $\omega_{\text{A}}$ are frequencies of optical and acoustic phonons at the edge of Brillouin zone, $N_{\text{O}}$ and $N_{\text{A}}$ are the corresponding phonon occupations, $\Omega_{\text{BZ}}$ is the Brillouin zone volume, $M$ is the average atom mass in a unit cell, $\epsilon_{\infty}$ is the high-frequency permittivity, and $Z$ is the Born effective charge of a single atom. For a typical $k$ such that $\epsilon_{k}=200$ meV, the e-2ph scattering rates at 300 K can be estimated as $10^{-2.5}$ ps^-1, which is about 3.5 orders of magnitude smaller than the ab-initio or semi-analytical 2ph scattering rates obtained in this work. Thus, we conclude that the effect of the simultaneous e-2ph interaction based on the Fröhlich interaction can be neglected.

Finally, we consider an alternate mechanism for the PSD peak which does not rely solely on intervalley scattering. The earliest studies of negative differential resistance in GaAs arose from the observation of current instabilities at electric fields approaching a threshold value of around 3 kVcm^-1. Gunn (1964) These instabilities were attributed to the formation of space charge domains associated with the negative differential resistance. The typical Boltzmann formalism used to describe charge transport from first principles would not include the contribution of such effects because it neglects real-space gradients and space charge effects which are essential to the instability.

Space charge instabilities manifest as current fluctuations, and therefore the nucleation of space-charge domains could explain the PSD peak around the threshold field. However, an inconsistency with this explanation is that the increase in PSD begins at a field below the threshold value for NDR as in Fig 2. This inconsistency can be accounted for by considering the possibility that the local electric field exceeds the threshold even though the average field does not. Such a possibility was investigated theoretically by McCumber and Chynoweth McCumber and Chynoweth (1966), who found that the dipole layer generation process was sensitive to inhomogeneities such as doping fluctuations that would arise from purely random Poisson statistics. The dipole layer was found to propagate even if the average uniform field was less than the nominal threshold field based on the velocity-field characteristics.

These considerations support the explanation of the PSD peak in terms of instabilities associated with the local electric field exceeding the threshold field for negative differential resistance space charge domain formation. This finding has implications for the use of transport studies to determine intervalley scattering strength. In particular, Monte Carlo methods have been used for decades for this purpose in GaAs by interpreting transport and noise measurements, and in those simulations, noise was assumed to arise from solely intervalley scattering. Our result indicates that this approach to studying intervalley scattering is not valid because the physical origin of noise differs from that assumed in the model. Instead, photoluminescence methods which directly provide an electronic lifetime as in Ref. Fasol et al. (1990) should be employed.