Electron BTE

Introduction to the BTE

Let \(f_{\lambda}\) be the out-of-equilibrium electron occupation number, where \(\lambda = (\boldsymbol{k},b)\) labels both electronic wavevectors and band index. To set up the electron BTE, we expand the electron occupations to first order,

\[f_{\lambda} = \bar{f}_{\lambda} + \delta f_{\lambda}\]

where \(\bar{f}_{\lambda}\) is the Fermi–Dirac distribution function and we introduced \(\delta f_{\lambda}\) is a small out-of-equilibrium population of electrons due to the presence of a field.

The linearized electronic BTE for a system exposed to an applied electric field and thermal gradient can then be written,

\[e \boldsymbol{v}_{\lambda} \cdot \boldsymbol{E} \frac{\partial \bar{f}_{\lambda}}{\partial \epsilon} + \boldsymbol{v}_{\lambda} \cdot \boldsymbol{\nabla} T \frac{\partial \bar{f}_{\lambda}}{\partial T} = - \sum_{\lambda'} \Omega_{\lambda\lambda'} \delta f_{\lambda'}\]

where the first term describes the diffusion due to an externally applied electric field \(\boldsymbol{E}\), the second the diffusion due to a temperature gradient, and the third term is the linearized scattering operator.

Electron scattering rates

The electron scattering matrix \(\Omega_{\lambda,\lambda'}\) can be computed using electron-phonon scattering as,

\[\begin{split}\Omega_{\boldsymbol{k}b,\boldsymbol{k}'b'} =& \frac{1}{\tau_{kb}} \delta_{kb,k'b'} + (1-\delta_{kb,k'b'}) \frac{2\pi}{N_k\hbar} \sum_{\boldsymbol{q}\nu} |g_{bb'\nu}(\boldsymbol{k},\boldsymbol{k}')|^2 \\ &\times \bigg[ ( 1 - \bar{f}_{\boldsymbol{k}b} + \bar{n}_{\boldsymbol{q}\nu}) \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} + \hbar \omega_{\boldsymbol{q}\nu}) \\ &+ (\bar{f}_{\boldsymbol{k}b} + \bar{n}_{\boldsymbol{q}\nu}) \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} - \hbar \omega_{\boldsymbol{q}\nu}) \bigg] \delta(\boldsymbol{k}-\boldsymbol{k}'+\boldsymbol{q}),\end{split}\]

with

\[\begin{split}\frac{1}{\tau_{kb}} =& \frac{2\pi}{N_k\hbar} \sum_{b'\boldsymbol{k}',\nu \boldsymbol{q}} |g_{bb'\nu}(\boldsymbol{k},\boldsymbol{k}')|^2 \times \bigg[ (1-\bar{f}_{\boldsymbol{k}'b'} + \bar{n}_{\boldsymbol{q}\nu}) \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} - \hbar \omega_{\boldsymbol{q}\nu}) \\ &+ (\bar{f}_{\boldsymbol{k}'b'} + \bar{n}_{\boldsymbol{q}\nu}) \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} + \hbar \omega_{\boldsymbol{q}\nu}) \bigg] \delta(\boldsymbol{k}-\boldsymbol{k}'+\boldsymbol{q}).\end{split}\]

where \(\bar{n}\) is the Bose–Einstein distribution function for phonons with wavevector, mode index \(\boldsymbol{q}\nu\) and energy \(\omega_{\boldsymbol{q}\nu}\). This scattering matrix requires us to know the phonon and electron energies, as well as the electron-phonon coupling, \(g(\boldsymbol{k},\boldsymbol{k+q})\), on a fine (interpolated) mesh, as discussed in Wannier interpolation of electron-phonon coupling.

The Kronecker delta function on momentum is enforces exactly, while the Dirac delta for conservation of energy is instead approximated with methods as described in the section Dirac delta approximations.

Symmetrization of the scattering matrix

Since the scattering matrix \(\Omega_{\lambda,\lambda'}\) is not inherently symmetric (which we will need if we want to have a real symmetric matrix, as in the realxons solution), we perform the transformation,

\[\tilde{\Omega}_{\lambda \lambda'} = \Omega_{\lambda \lambda'} \sqrt{ \frac{\bar{f}_{\lambda'} (1-\bar{f}_{\lambda'})}{\bar{f}_{\lambda} (1-\bar{f}_{\lambda})} }\]

which results in the matrix with diagonal matrix elements:

\[\tilde{\Omega}_{\lambda \lambda} = \Omega_{\lambda \lambda} = \frac{1}{\tau_{\boldsymbol{k}b}}\]

and for the off-diagonal terms:

\[\begin{split}\tilde{\Omega}_{\boldsymbol{k}b,\boldsymbol{k}'b'} =& - \frac{2\pi}{V N_k} \sum_{\boldsymbol{q}\nu} |g_{bb'\nu}(\boldsymbol{k},\boldsymbol{k}')|^2 \delta(\boldsymbol{k}-\boldsymbol{k}'+\boldsymbol{q}) \\ &\times \bigg[ \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} + \hbar \omega_{\boldsymbol{q}\nu}) + \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} - \hbar \omega_{\boldsymbol{q}\nu}) \bigg] \frac{1}{2 \sinh{\big( \frac{\hbar \omega_{\boldsymbol{q}\nu}}{2 k_B T} \big) } } .\end{split}\]

This matrix is symmetric and has a number of interesting physical properties (e.g. the eigenvalues corresponding to the exact relaxation times of the electronic bath). Computationally, the symmetric matrix can be used in a conjugate gradient method that maximises the electrical and thermal conductivity, and guarantees the existence of eigenvalues.

In Phoebe, instead of solving the original BTE problem in the form \(\sum_{\lambda'} \Omega_{\lambda,\lambda'} \delta f_{\lambda'} = b_{\lambda}\), we solve the symmetrized problem:

\[\sum_{\lambda'} \tilde{\Omega}_{\lambda,\lambda'} \delta \tilde{f}_{\lambda'} = \tilde{b}_{\lambda}\]

with

\[\delta f_{\lambda} = ( \bar{f}_{\lambda} (1-\bar{f}_{\lambda}) )^{-\frac{1}{2}} \delta f_{\lambda}\]

and

\[\tilde{b}_{\lambda} = ( \bar{f}_{\lambda} (1-\bar{f}_{\lambda}) )^{-\frac{1}{2}} b_{\lambda}\]

Onsager coefficients

In the electronic case, there are a handful of transport coefficients we’d like to solve the BTE to find, such as the electrical conductivity, \(\sigma\), the electronic part of the thermal conductivity, \(\kappa_e\), and the Seebeck coeffieint, \(S\). These quantities are defined in terms of the Onsager coefficients.

We assume that the response to the applied electric field and thermal gradient is linear in these external fields:

\[\delta f_{\lambda} = \sum_{i} \delta^i f^E_{\lambda} E_i + \delta^i f^T_{\lambda} \nabla_i T\]

After computing the out-of-equilibrium population, the charge and heat flux density can be computed as:

\[\boldsymbol{J} = \frac{e g_s}{V N_k} \sum_{\lambda} \boldsymbol{v}_{\lambda} f_{\lambda}\]

and

\[\boldsymbol{Q} = \frac{g_s}{V N_k} \sum_{\lambda} (\epsilon_{\lambda}-\mu) \boldsymbol{v}_{\lambda} f_{\lambda}\]

where \(g_s\) is the spin degeneracy.

We can decompose these to write,

\[\boldsymbol{J} = L_{EE} \boldsymbol{E} + L_{ET} \boldsymbol{\nabla} T\]

\[\boldsymbol{Q} = L_{TE} \boldsymbol{E} + L_{TT} \boldsymbol{\nabla} T\]

From this, the electrical conductivity \(\sigma\), the thermal conductivity \(k\), the Seebeck coefficient \(S\) and the mobility \(\mu\) are:

\[\sigma = L_{EE}\]

\[k = L_{TT} - L_{TE} L_{EE}^{-1} L_{ET}\]

\[S = - L_{EE}^{-1} L_{ET}\]

\[\mu = \frac{\sigma}{n}\]

where \(n\) is the carrier concentration.

Solutions of the electron BTE

Largely, these solvers follow the equivalent section in the phonon BTE section, where they may be described in more detail. For further details and references on any specific solver, we suggest you visit the equivalent phonon sections, as well. Here, we again establish methods of finding the solution vector to the BTE, \(f\), but in this case, we have two: \(f^T\) and \(f^E\), for each field.

RTA Solution

The relaxation time approximation (RTA) approximates the scattering matrix to be diagonal,

\[A_{ \boldsymbol{k}b,\boldsymbol{k}b } \approx \frac{1}{ \tau_{\boldsymbol{k}b} }\]

which results in a trivial solution to the Boltzmann transport equation,

\[e \boldsymbol{v}_{\lambda} \cdot \boldsymbol{E} \frac{\partial \bar{f}_{\lambda}}{\partial \epsilon} + \boldsymbol{v}_{\lambda} \cdot \boldsymbol{\nabla} T \frac{\partial \bar{f}_{\lambda}}{\partial T} = - \frac{1}{ \tau_{\lambda} } \delta f_{\lambda}.\]

Solving separately for the response to the electric field and the thermal gradient, we find,

\[\delta^i f^E_{\lambda} = - e v^i_{\lambda} \frac{\bar{f}_{\lambda}(1-\bar{f}_{\lambda})}{k_B T} \tau_{\lambda}\]

\[\delta^i f^T_{\lambda} = - v^i_{\lambda} \frac{(\epsilon_{\lambda}-\mu)\bar{f}_{\lambda}(1-\bar{f}_{\lambda})}{k_B T^2} \tau_{\lambda}\]

MRTA Solution

Like the RTA, the momentum relaxation time approximation (MRTA) considers the scattering matrix to be diagonal. However, in the case of the MRTA, we calculate the transport coefficients using only processes which degrade the electrical current. This is accounted for using an additional term in the calculation of the relaxation times, for example, in the electron-phonon case,

\[\begin{split}\frac{1}{\tau^{\mathrm{MRTA}}_{kb}} =& \frac{2\pi}{N_k\hbar} \sum_{b'\boldsymbol{k}',\nu \boldsymbol{q}} |g_{bb'\nu}(\boldsymbol{k},\boldsymbol{k}')|^2 \times \bigg[ (1-\bar{f}_{\boldsymbol{k}'b'} + \bar{n}_{\boldsymbol{q}\nu}) \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} - \hbar \omega_{\boldsymbol{q}\nu}) \\ &+ (\bar{f}_{\boldsymbol{k}'b'} + \bar{n}_{\boldsymbol{q}\nu}) \delta(\epsilon_{\boldsymbol{k}b} - \epsilon_{\boldsymbol{k}'b'} + \hbar \omega_{\boldsymbol{q}\nu}) \bigg] \bigg(1 - \frac{v_{\boldsymbol{k}b}\cdot v_{\boldsymbol{k}'b'}}{|v_{\boldsymbol{k}b}||v_{\boldsymbol{k}'b'}|} \bigg) \delta(\boldsymbol{k}-\boldsymbol{k}'+\boldsymbol{q}).\end{split}\]

where the term \(1 - \frac{v_{\boldsymbol{k}b}\cdot v_{\boldsymbol{k}'b'}}{|v_{\boldsymbol{k}b}||v_{\boldsymbol{k}'b'}|}= 1 - cos\theta\), where theta is the scattering angle. This solution therefore effectively only counts the contributions of electrons which undergo backscattering processes.

Iterative solution

Note

Generally, we recommend the variational method over this for numerical stability. However, the iterative solution allows us to use symmetries as in Chaput, PRL (2013)., which can decrease cost.

To better understand this method, please have a look first at the counterpart iterative solution in the phonon documentation. In short, the electron BTE consists in two linear algebra problems:

\[m^{i}_{\lambda} = - \sum_{\lambda'} A_{\lambda\lambda'} \delta f_{\lambda}^E\]

\[n^{i}_{\lambda} = - \sum_{\lambda'} A_{\lambda\lambda'} \delta f_{\lambda}^T\]

where

\[m^{i}_{\lambda} = e v_{\lambda}^i \frac{\partial \bar{f}_{\lambda}}{\partial \epsilon}\]

\[n^{i}_{\lambda} = v_{\lambda}^i \frac{\partial \bar{f}_{\lambda}}{\partial T}\]

The iterative scheme solves these two independent linear algebra problems with a geometric series,

\[\delta^i f^E_{K} = \sum_{K} \left(-\frac{1}{\boldsymbol{A}^{\mathrm{out}}} \boldsymbol{A}^{\mathrm{in}}\right)^{K} \frac{1}{\boldsymbol{A}^{\mathrm{out}}} \: m^i\]

and

\[\delta^i f^T_K = \sum_{K} \left(-\frac{1}{\boldsymbol{A}^{\mathrm{out}}} \boldsymbol{A}^{\mathrm{in}}\right)^{K} \frac{1}{\boldsymbol{A}^{\mathrm{out}}} \: n^i\]

where \(K\) is an iteration index, \(A^{in}\) is the off-diagonal part of the scattering matrix, and \(A^{out}\) is the diagonal part of the scattering matrix. In the code, the two problems are solved together, as we compute the action on the two different vectors at the same time.

Note that, like any geometric series, this algorithm may not converge.

Iterative solution: Variational/Conjugate Gradient Method

Again, this solver is very similar to the phonon case (and we recommend you read more there as well). The only difference for electronic systems is that we need to solve two problems simultaneously, one for the electric field response and one for the response to the thermal gradient.

For the variational method, we can define the variational thermal conductivity, in closed-circuit conditions, as:

\[k^\mathrm{V}(\delta f^T) = - 2 \mathcal{T}({\delta f^T})\]

where

\[\mathcal{T}(\delta f^T) = \frac{1}{2} \sum_{\lambda \lambda'} {\delta f^T_{\lambda}} \cdot{\boldsymbol A_{\lambda\lambda'}} {\delta f^T_{\lambda'}} - \sum_{\lambda} {\boldsymbol n_{\lambda}} \cdot {\delta f^T_{\lambda}}\]

The variational electrical conductivity is defined similarly as:

\[\sigma^\mathrm{V}(\delta f^E) = 2 \mathcal{E}({\delta f^E})\]

where

\[\mathcal{E}(\delta f^E) = \frac{1}{2} \sum_{\lambda \lambda'} {\delta f^E_{\lambda}} \cdot{\boldsymbol A_{\lambda\lambda'}} {\delta f^E_{\lambda'}} - \sum_{\lambda} {\boldsymbol m_{\lambda}} \cdot {\delta f^E_{\lambda}}\]

These two functionals are the minimization targets of a conjugate gradient method. Knowing this, the variational method is exactly the same as the phonon case, with the proper substitution of the vector b with either \(m\) or \(n\).

As in the case of the simple iterative method, we solve the two equations (response to electric field and thermal gradient) at the same time, as it allows us to minimize the number of times the scattering matrix is evaluated (the most expensive step).

Relaxons Method (Direct Diagonalization)

We first diagonalize the scattering matrix,

\[\frac{1}{N_k} \sum_{\lambda'} A_{\lambda\lambda'} \theta_{\lambda'\alpha} = \frac{1}{\tau_{\alpha}} \theta_{\lambda\alpha}\]

where \(\theta\) are eigenvectors, \(\alpha\) are eigenvalue indices, and \(\frac{1}{\tau_{\alpha}}\) are eigenvalues. We first build the auxiliary quantities:

\[\delta^i f^E_{\alpha} = \sum_{\lambda} \frac{\partial \bar{f}_{\lambda}}{\partial \epsilon} v_{\lambda}^i \theta_{\lambda \alpha} \tau_{\alpha}\]

\[\delta^i f^T_{\alpha} = \sum_{\lambda} \frac{\partial \bar{f}_{\lambda}}{\partial T} v_{\lambda}^i \theta_{\lambda \alpha} \tau_{\alpha}\]

From these, we can compute the solutions of the BTE as:

\[\delta f^E_{\lambda} = \frac{1}{N_k V} \sum_{\alpha} f^E_{\alpha} \theta_{\lambda \alpha}\]

\[\delta f^T_{\lambda} = \frac{1}{N_k V} \sum_{\alpha} f^T_{\alpha} \theta_{\lambda \alpha}\]

Wigner correction to the electron BTE

The theory developments for the Wigner corrections to the electron BTE are described in Materials Today Physics 19, 100412 (2021). The Wigner transport equation is,

\[\begin{split}\frac{\partial f_{bb'}(\boldsymbol{x},\boldsymbol{k},t)}{\partial t} &+ \frac{i}{\hbar} \Big[ \mathcal{E}(\boldsymbol{k}) + \boldsymbol{D}(\boldsymbol{k})\cdot\boldsymbol{E} , f(\boldsymbol{x},\boldsymbol{k},t) \Big]_{bb'} + \frac{1}{2} \Big\{ \boldsymbol{v}(\boldsymbol{k}) , \cdot \frac{\partial f(\boldsymbol{x},\boldsymbol{k},t)}{\partial \boldsymbol{x}} \Big \}_{bb'} \\\\ &+ e \boldsymbol{E} \cdot \frac{\partial f_{bb'}(\boldsymbol{x},\boldsymbol{k},t)}{\partial \boldsymbol{k}} = -\frac{\partial f_{bb'}(\boldsymbol{x},\boldsymbol{k},t)}{\partial t} \bigg|_{coll}\end{split}\]

where \(f_{bb'}(\boldsymbol{x},\boldsymbol{k},t)\) is the Wigner distribution function, \({ \cdot,\cdot }\) indicates an anticommutator, \([ \cdot,\cdot ]\) indicates a commutator, \(v_{bb'}(\boldsymbol{k})\) is the velocity operator, and we defined the matrix \(\mathcal{E}(\boldsymbol{k})_{bb'} = \delta_{bb'} \epsilon_{\boldsymbol{k}b}\) and \(\mathcal{D}(\boldsymbol{k})_{bb'} = (1-\delta_{bb'}) d_{\boldsymbol{k}bb'}\) is a matrix of electronic dipoles. The electronic dipole can be computed as:

\[\boldsymbol{d}_{\boldsymbol{k},bb'} = - i e \frac{\boldsymbol{v}_{bb'}(\boldsymbol{k})}{\epsilon_{b}(\boldsymbol{k})-\epsilon_{b'}(\boldsymbol{k})} , \quad \text{for }b \neq b'\]

The scattering operator acts on the diagonal Wigner distribution as the BTE scattering operator, instead it acts on the off-diagonal components with a decay term:

\[\frac{\partial f_{bb'}(\boldsymbol{x},\boldsymbol{k},t)}{\partial t} \bigg|_{coll} = (1-\delta_{bb'}) \frac{\Gamma_{b}(\boldsymbol{k}) + \Gamma_{b'}(\boldsymbol{k})}{2} f_{bb'}(\boldsymbol{x},\boldsymbol{k},t) + \delta_{bb'} \frac{1}{V} \sum_{\boldsymbol{k}'b'} A_{\boldsymbol{k}b,\boldsymbol{k}'b'} f_{b'b'}(\boldsymbol{x},\boldsymbol{k}',t)\]

where \(\Gamma_b(\boldsymbol{k}) = \frac{2\pi}{\tau_{\boldsymbol{k}b}}\) are the electronic linewidths.

To solve the Wigner transport equation, just like we did for the BTE, we assume linear response and separate the response to electric field and thermal gradient \(f = f^E E + f^T \nabla T\). The diagonal part of the Wigner transport equation is exactly equal to the BTE, and can be solved using one of solvers described above. The off-diagonal part of the Wigner distribution function can be solved easily with a little algebraic manipulation.

The related transport coefficients are defined as:

\[L_{EE}^{ij} = \frac{e g_s}{V N_k} \sum_{\boldsymbol{k}b} \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{E_j}(\boldsymbol{k}) \Big\}_{bb}\]

\[L_{ET}^{ij} = \frac{e g_s}{V N_k} \sum_{\boldsymbol{k}b} \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{T_j}(\boldsymbol{k}) \Big\}_{bb}\]

\[L_{TE}^{ij} = \frac{g_s}{V N_k} \sum_{\boldsymbol{k}b} \big( \epsilon_{b}(\boldsymbol{k})-\mu \big) \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{E_j}(\boldsymbol{k}) \Big\}_{bb}\]

\[L_{TT}^{ij} = \frac{g_s}{V N_k} \sum_{\boldsymbol{k}b} \big( \epsilon_{b}(\boldsymbol{k})-\mu \big) \frac{1}{2} \Big\{ v^i(\boldsymbol{k}) , f^{T_j}(\boldsymbol{k}) \Big\} _{bb}\]