Semiclassical Electron and Phonon Transport from First Principles: Application to Layered Thermoelectrics

Owing to the relative simplicity, reliability, reasonable computational effort, and extensive implementation in widely available software packages [5], DFT has rapidly emerged as the method of choice to find approximate solutions to the many-body problem and determine the properties of molecules and materials at the quantum-mechanical level [6]. DFT is conceptually rooted in the two theorems introduced by Hohenberg and Kohn [7], and its practical deployment follows the KS Hamiltonian [8],

$$\hat{H}^{KS}=-\frac{1}{2}\nabla^{2}+V(\mathbf{r})+V_{\textsubscript{H}}(\mathbf{r})+V^{xc}(\mathbf{r}),$$

(1)

where $V(\mathbf{r})$ is the external potential, $V_{\textsubscript{H}}(\mathbf{r})$ the Hartree potential, and $V^{xc}(\mathbf{r})$ the exchange-correlation potential. This Hamiltonian leads to a set of $n$ Schrödinger-like equations for the single-particle Kohn-Sham orbitals $\Psi(\mathbf{r})$ and accompanying energies $\epsilon$,

$$\hat{H}^{KS}\Psi_{i}(\mathbf{r})=\epsilon_{i}\Psi_{i}(\mathbf{r}).$$

(2)

The solution of the Kohn-Sham equations proceeds in self-consistent fashion by (i) starting with an educated guess of the electron density $n_{k}(\mathbf{r})$ to (ii) construct the Kohn-Sham Hamiltonian $\hat{H}^{KS}$, (iii) finding upon diagonalization the corresponding eigenvalues $\epsilon_{i}$ and eigenvectors, which are subsequently used to (iv) obtain the new density $n_{k+1}(\mathbf{r})$, until the difference between $n_{k+1}(\mathbf{r})$ and $n_{k}(\mathbf{r})$ does not exceed a given numerical tolerance.

A number of approximations to the exchange-correlation potential appearing in Equation 1 has been devised [9, 10], most notably (i) the Local Density Approximation (LDA) [11], in which the electron density is assumed to be locally the same as a spatially uniform electron gas with the same density, (ii) the Generalized Gradient Approximation (GGA) [12, 13, 14, 15] and Meta-Generalized Gradient Approximation (meta-GGA) [16, 17], and Laplacian of the electron density or the kinetic energy density, are taken into account to increase the accuracy over LDA, and (iii) hybrid functionals [18, 19, 20], in which a fraction of the orbital-dependent, exact Fock exchange is included in the otherwise bare GGA.

Although the exchange interactions between electrons can be determined exactly (cf. Hartree-Fock equations), correlation interactions are invariably approximated, a trait that is common to all the exchange-correlation functionals listed above. It is therefore not surprising that density-functional theory provides an inaccurate account of those physical situations where correlation effects dominate, for example, van der Waals interactions or localized electronic states. In addition, being DFT commonly considered a ground state theory (cf. Hoehnberg-Kohn theorems), predictions concerning excited states and band gaps are often quantitatively unreliable, see references [21, 22] for through discussions.

When periodic boundary conditions apply, like in the case of crystalline materials, one can take advantage of the Bloch theorem, which allows the single-particle wavefunctions to be labeled by a band index $n$ and the crystal momentum $\mathbf{k}$ and written as the product of a plane wave and a function that is periodic in the crystal lattice:

$$\phi_{n,\mathbf{k}}(\mathbf{r})=\psi_{n,\mathbf{k}}({\mathbf{r}})e^{i\mathbf{k}\cdot\mathbf{r}},\qquad\psi_{n,\mathbf{k}}(\mathbf{r}+\mathbf{R})=\psi_{n,\mathbf{k}}(\mathbf{r})\,.$$

(3)

Because of this periodicity, $\psi_{n,\mathbf{k}}(\mathbf{r})$ can be expanded in a Fourier series including only reciprocal lattice vectors, $\mathbf{G}$:

$$\psi_{n,\mathbf{k}}(\mathbf{r})=\sum_{\mathbf{G}}e^{i\mathbf{G}\cdot\mathbf{r}}\psi_{n,\mathbf{k}}(\mathbf{G})\,.$$

(4)

Thus, plane waves are the most natural basis set for periodic systems. According to Eq. (3), however, the evaluation of the solution at a single point requires the summation over an infinite number of $\mathbf{G}$ vectors. In practice, this expansion is truncated to an energy cutoff, $E\textsubscript{cut}$, implying that only solutions of kinetic energy

$$\frac{h^{2}}{2m_{e}}|\mathbf{k}+\mathbf{G}|^{2}<E\textsubscript{cut}$$

(5)

are evaluated. Delocalized plane waves further offer the advantage of completeness, and the convergence of the calculated properties can be systematically improved by increasing the energy cutoff ¹¹1Besides plane waves, localized basis sets consisting of atomic-like orbitals (e.g., Gaussian- or Slater-type functions) have found a widespread use, in particular in the computational chemistry community. Contrary to plane waves, fewer basis functions are often needed to achieve a reasonable accuracy, hence significantly decreasing the computational effort. However, localized basis sets are controlled by many parameters in addition to the energy cutoff, in a way that no systematic convergence can be attained.. A list of representative software packages implementing plane-wave density-functional theory for solid-state systems is given in Table 1.

Although extended Bloch orbitals are the most natural choice to describe the electronic states of a periodic system, an alternative approach based on the Wannier representation is convenient for an array of applications in electronic structure [32, 33]. Contrary to Bloch orbitals, Wannier functions provide a real-space representation of localized orbitals which are specified by a lattice vector $R$ and a band index $n$. The general transformation from Bloch orbitals into Wannier functions is

$$\ket{W_{\mathbf{R},n}}=\frac{V}{(2\pi)^{2}}\int_{\textnormal{BZ}}d\mathbf{k}\sum_{m=1}^{N}U_{m,n}^{\mathbf{k}}\psi_{m,\mathbf{k}}e^{{-i\mathbf{k}\cdot\mathbf{R}}},$$

(6)

where $V$ is the volume of the primitive cell and $U_{m,n}^{\mathbf{k}}$ is a gauge-fixing, $N$-dimensional unitary matrix. Additional freedom in choice of $U_{m,n}^{\mathbf{k}}$ means that the Wannier functions are not uniquely determined. A popular choice for the gauge is determined by minimizing the spread functional, $\Omega$:

$$\frac{\partial\Omega[U]}{\partial U}=0,$$

(7)

where

$$\Omega=\sum_{n}\left[\bra{W_{0,n}}r^{2}\ket{W_{0,N}}-\bra{W_{0,n}}r\ket{W_{0,n}}^{2}\right].$$

(8)

Wannier functions obtained following this procedure are referred to as Maximally Localized Wannier Functions (MLWF). Wannierization of DFT calculations is typically performed with the Wannier90 package [34, 35].

Within the Born-Oppenheimer (BO) approximation, the lattice-dynamical properties of a system can be determined by obtaining the eigenvalues, $E$, and the eigenvectors, $\Phi$, of the following Schrödinger equation [36]

$$\left(\sum_{I}\frac{\hbar^{2}}{2M_{I}}\frac{\partial^{2}}{\partial{\bf{R}}_{I}^{2}}+\hat{H}_{BO}({\bf{R}})\right)\Phi({\bf{R}})=E({\bf{R}})\Phi({\bf{R}})~{}.$$

(9)

Here, the index $I$ labels each ion in the system, while ${\bf{R}}_{I}$ and $M_{I}$ correspond to its coordinates and masses, respectively. For conciseness, $\mathbf{R}$ correspond to all ionic coordinates and $\hat{H}_{BO}$ is the BO electronic Hamiltonian that depends parametrically upon $\mathbf{R}$

$$\hat{H}_{BO}({\bf{R}})=\frac{-\hbar^{2}}{2m}\sum_{i}\frac{\partial^{2}}{\partial{\bf{r}}_{i}^{2}}+\frac{e^{2}}{2}\sum_{i\neq j}\frac{1}{|{\bf{r}}_{i}-{\bf{r}}_{j}|}-\sum_{iI}\frac{Z_{I}e^{2}}{|{\bf{r}}_{i}-{\bf{R}}_{I}|}+\frac{e^{2}}{2}\sum_{I\neq J}\frac{Z_{I}Z_{J}}{|{\bf{R}}_{I}-{\bf{R}}_{J}|}~{}.$$

(10)

Interatomic forces are computed from the derivative of the ground-state energy with respect to the ion positions, and can be calculated by using the Hellmann-Feynman theorem [37, 38]:

	$\displaystyle{\bf{F}}_{I}=-\frac{\partial E({\bf{R}})}{\partial{\bf{R}}_{I}}$	$\displaystyle=$	$\displaystyle-\braket{\Phi({\bf{R}})}{\frac{\partial(\hat{H}_{BO})}{\partial{\bf{R}}_{I}}}{\Phi({\bf{R}})}$		(11)
		$\displaystyle=$	$\displaystyle-\int n_{{\bf{R}}}({\bf{r}})\frac{\partial V_{{\bf{R}}}({\bf{r}})}{\partial{{\bf{R}}}_{I}}d{\bf{r}}-\frac{\partial E_{N}({{\bf{R}}})}{\partial{{\bf{R}}}_{I}}~{},$		(12)

where ${\bf{R}}_{I}={\bf{R}}_{p}+\xi_{\kappa\alpha p}$, with ${\bf{R}}_{p}$ one of the direct lattice vectors of a Born-von Kármán supercell composed of $N_{p}$ unit cells of the crystal, and $\xi_{\kappa\alpha p}$ is the position vector with Cartesian coordinate $\alpha$ of the nucleus $\kappa$ in the unit cell labeled by $p$. In Eq. (12) we have

$$V_{{\bf{R}}}({\bf{r}})=\sum_{iI}\frac{Z_{I}e^{2}}{|{\bf{r}}_{i}-{\bf{R}}_{I}|}~{},$$

(13)

which is the external potential, and $n_{{\bf{R}}}({\bf{r}})$ is the ground-state electron charge density for a given geometry. $E_{N}({{\bf{R}}})$ is the electrostatic interaction between nuclei

$$E_{N}({{\bf{R}}})=\frac{e^{2}}{2}\sum_{I\neq J}\frac{Z_{I}Z_{J}}{|{\bf{R}}_{I}-{\bf{R}}_{J}|}~{}.$$

(14)

Equation (12) assumes the external potential acting on the electrons is a differentiable function of the crystal coordinates. In particular, the equilibrium geometry is reached when forces acting on individual ion vanish, that is, ${\bf{F}}_{I}\equiv\frac{\partial E({\bf{R}})}{\partial{\bf{R}}_{I}}=0$. The vibrational properties of a system can be determined by the interatomic force constants, which correspond to the second derivatives of $E({\bf{R}})$ in relation to ion displacements and can be calculated via Hellmann-Feynman theorem as well

	$\displaystyle{\bf{C}}_{\kappa\alpha p,\kappa^{\prime}\alpha^{\prime}p^{\prime}}$	$\displaystyle=$	$\displaystyle\frac{\partial^{2}E}{\partial{\bf{\xi}}_{\kappa\alpha p}\partial{\bf{\xi}}_{\kappa^{\prime}\alpha^{\prime}p^{\prime}}}$		(15)
	$\displaystyle{\bf{C}}_{I,J}$	$\displaystyle=$	$\displaystyle\frac{\partial^{2}E({{\bf{R}}})}{\partial{\bf{R}}_{I}\partial{\bf{R}}_{J}}=-\frac{\partial{\bf{F}}_{I}}{\partial{\bf{R}}_{J}}$		(16)
		$\displaystyle=$	$\displaystyle\int\frac{\partial n_{{\bf{R}}}({\bf{r}})}{\partial{\bf{R}}_{J}}\frac{\partial V_{{\bf{R}}}({\bf{r}})}{\partial{\bf{R}}_{I}}d{\bf{r}}+\int n_{{\bf{R}}}({\bf{r}})\frac{\partial^{2}V_{{\bf{R}}}({\bf{r}})}{\partial{\bf{R}}_{I}\partial{\bf{R}}_{J}}d{\bf{r}}+\frac{\partial^{2}E_{N}({\bf{R}})}{\partial{\bf{R}}_{I}\partial{\bf{R}}_{J}}~{},$		(17)

where $\frac{\partial n_{{\bf{R}}}({\bf{r}})}{\partial{\bf{R}}}$ is the linear response of $n_{{\bf{R}}}({\bf{r}})$ to a distortion of the nuclear geometry, which is a fundamental result derived by De Cicco and Johnson [39] and by Pick, Cohen, and Martin [40]. This establishes that, within the adiabatic approximation, the electrons feel only a static phonon perturbation. It is important to note that, due to translational invariance, the interatomic force constants depend on $p$ and $p^{\prime}$ only through the difference ${\bf{R}}_{p}-{\bf{R}}_{p^{\prime}}$. The Fourier transform of the interatomic force constants

$$D_{\kappa\alpha,\kappa^{\prime}\alpha^{\prime}}({\bf{q}})=\frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}}\sum_{p}{\bf{C}}_{\kappa\alpha 0,\kappa^{\prime}\alpha^{\prime}p}e^{i{\bf{q}}\cdot{\bf{R}}_{p}}~{},$$

(18)

is the (hermitian) dynamical matrix [41], where $M_{\kappa}$ is the mass of the ion $\kappa$. Because of the translational invariance, the lattice distortion is monochromatic, meaning that a phonon perturbation with wavevector $q$ does not induce a force response with wave vector ${\bf{q^{\prime}}}\neq{\bf{q}}$. The phonon eigenfrequencies of the phonon mode $\nu$ and wave vector ${\bf{q}}$, $\omega_{\nu{\bf{q}}}$, are obtained by diagonalizing the dynamical matrix

$$\sum_{\kappa^{\prime}\alpha^{\prime}}D_{\kappa\alpha,\kappa^{\prime}\alpha^{\prime}}({\bf{q}})e_{\kappa^{\prime}\alpha^{\prime},\nu}({\bf{q}})=\omega^{2}_{\nu{\bf{q}}}e_{\kappa\alpha,\nu}({\bf{q}})~{},$$

(19)

and the eigenvectors $e_{\kappa^{\prime}\alpha^{\prime},\nu}({\bf{q}})$ are the phonon polarizations.

We now briefly introduce density functional perturbation theory (DFPT), which can be used to calculate the interatomic force constants and vibrational spectra as well as the electron-phonon matrix elements. As presented in the last section, DFT allows for the calculation of the total energy and charge density, which can in turn be used to obtain a large number of experimental observables. DFPT allows one to obtain quantities that depend on derivatives of the total energy and charge density with respect to a small change in the potential. In particular, DFPT enables the calculation of different properties in condensed matter physics such as vibrational frequencies, elastic constants, dielectric tensors, Born effective charges, piezoelectric tensors and flexoelectricity [42, 43, 44, 45, 46, 47, 48, 49, 50, 51].

First-order perturbative approaches within DFT were proposed independently by different groups working in different contexts [52, 53, 54, 55, 56, 57]. Specifically, the investigation of perturbations in condensed matter physics, such as the response of fermionic systems to displacement of atoms within DFT, has been pionered by Zein [56], Baroni, Giannozzi and Testa [57], and Gonze [58]. A generalization to arbitrary order was introduced by Gonze and Vigneron [59], on the basis of the $2n+1$ theorem of perturbation theory.[60] Such perturbative approaches are based on different techniques comprising Green’s functions [57], the generalized Sternheimer equation [54, 61], or the Hylleraas variational scheme [62, 63]. Below we will follow the generalized Sternheimer approach. Connections to other techniques have been discussed extensively elsewhere [62].

Equation (17) demonstrates that the matrix of force constants are determined by the electron-density linear response, $\partial n$. Here we will write the KS eigenfunctions $\Psi_{n{\bf{k}}}$ in the Bloch form

$$\Psi_{n{\bf{k}}}({\bf{r}})=N_{p}^{-1/2}\psi_{n{\bf{k}}}({\bf{r}})e^{i{\bf{k}}\cdot{\bf{r}}}~{}.$$

(20)

with lattice periodic part $\psi_{n{\bf{k}}}$. The density $n$ is given by

$$n({\bf{r}})=\sum_{v{\bf{k}}}|\Psi_{v{\bf{k}}}({\bf{r}})|^{2}$$

(21)

where the band index $v$ indicates occupied states only.

Within DFPT, $\partial n$ is induced by the first-order variation of the KS potential, $\partial_{\kappa\alpha,{\bf{q}}}v^{KS}({\bf{r}})e^{i{\bf{q}}\cdot{\bf{r}}}$, and can be calculated through the first-order variation of the lattice periodic KS wave functions, $\partial\psi_{n{\bf{k}},\bf{q}}e^{i{\bf{q}}\cdot{\bf{r}}}$:

$$\partial n_{\kappa\alpha,{\bf{q}}}({\bf{r}})=\frac{2}{N_{p}}\sum_{v{\bf{k}}}\psi_{v{\bf{k}}}^{\ast}\partial\tilde{\psi}_{v{\bf{k}},{\bf{q}}}.$$

(22)

In Eq. (22), $\partial\tilde{\psi}_{v\mathbf{k},\mathbf{q}}=P^{C}\partial\psi_{v\mathbf{k},\mathbf{q}}$ is the projection of the first-order variation of the lattice periodic KS wave functions onto empty states, where $P^{C}=(1-P^{V})$ and $P^{V}=\sum_{v}\ket{\psi_{v{\bf{k}}+{\bf{q}}}}\bra{\psi_{v{\bf{k}}+{\bf{q}}}}$ is the projector onto filled valence states.

The first-order expansion of the KS equations gives the Sternheimer equation [62]

$$\left(\hat{H}^{KS}_{{\bf{k}}+{\bf{q}}}-\epsilon_{v{\bf{k}}}\right)\partial\psi_{v{\bf{k}},{\bf{q}}}=-\left(\partial_{\kappa\alpha{\bf{q}}}v^{KS}-\partial\epsilon_{v{\bf{k}}}\right)\psi_{v{\bf{k}}}~{},$$

(23)

in which the unperturbed wave functions were written by explicitly indicating the wavevector ${\bf{k}}$ and band index $v$, while the perturbed wavefunctions are projected onto the manifold of wave vectors ${\bf{k}}+{\bf{q}}$. $\hat{H}^{KS}_{{\bf{k}}+{\bf{q}}}=e^{-i({\bf{k}}+{\bf{q}})\cdot{\bf{r}}}\hat{H}^{KS}e^{i({\bf{k}}+{\bf{q}})\cdot{\bf{r}}}$ where $\hat{H}^{KS}$ is the unperturbed KS Hamiltonian, Eq. (1), and $\partial_{\kappa\alpha,{\bf{q}}}v^{KS}$ is the first-order variations of the KS potential

$$\partial_{\kappa\alpha{\bf{q}}}v^{KS}({\bf{r}})=\partial_{\kappa\alpha{\bf{q}}}V({\bf{r}})+e^{2}\int\frac{\partial_{\kappa\alpha{\bf{q}}}n({\bf{r^{\prime}}})}{|{\bf{r}}-{\bf{r^{\prime}}}|}e^{-i{\bf{q}}\cdot({\bf{r}}-{\bf{r^{\prime}}})}d{\bf{r^{\prime}}}+\frac{d_{\kappa\alpha{\bf{q}}}V^{xc}}{dn}\bigg{|}_{n=n({\bf{r}})}\partial_{\kappa\alpha{\bf{q}}}n({\bf{r}})~{}.$$

(24)

By expressing $\partial n$ in terms of a sum over the whole spectrum of the unperturbed Hamiltonian, both occupied and empty states, it becomes evident that contributions to $\partial n$ coming from only occupied states will vanish [44]. Only perturbations that couple the occupied-state manifold with the empty-state manifold will contribute. At the same time, the right-hand side of Eq. (23) depends only on occupied states, while its left-hand side is ill-conditioned because the possibility of null eigenvalues of the linear operator. In order to make Eq. (23) nonsingular, both sides are projected onto the empty conduction states by the projector $P^{C}$. Additionally, by adding $P^{V}P^{C}\partial\psi_{v}=0$ to the left-hand side of Eq. (23) we remove any null eigenvalues without changing the equation:

$$\left(\hat{H}^{KS}_{{\bf{k}}+{\bf{q}}}+\alpha P^{V}_{{\bf{k}}+{\bf{q}}}-\epsilon_{v{\bf{k}}}\right)\partial\tilde{\psi}_{v{\bf{k}},{\bf{q}}}=-(1-P^{V}_{{\bf{k}}+{\bf{q}}})\partial_{\kappa\alpha,{\bf{q}}}v^{KS}\psi_{v{\bf{k}}}~{}.$$

(25)

In practice, Eqs. (22), (24) and (25) are solved self-consistently, just like the KS equations within DFT. They form a generalized linear problem, since the KS perturbation, $\partial v^{KS}$, depends linearly on the electron-density linear response, $\partial n$, which in turn is a linear functional of the variation of the KS wavefunctions, $\partial\psi$. In practice, the initial KS perturbation is set to be equal to the external potential, Eq. (13). By solving the Sternheimer equation one obtains the induced electron-density linear response, Eq. (22), which is used to get first-order perturbations of the Hartree and exchange-correlation potentials. Equation (24) then defines a new first-order perturbation of KS potential in the Sternheimer equation. This procedure is repeated until the convergence of $\partial n$ is reached.

Within DFPT the responses to perturbations of different wavelengths are decoupled, which allows one to calculate vibrational properties at any wavevector without resorting to large supercells.[61] In some cases this can be a strength when compared to other methods such as frozen-phonon (FP)[64, 65] or molecular dynamics (MD) methods [66, 67]. However, in the case of surfaces, interfaces, or defects, both FP and MD are suitable because the systems naturally require large supercells. Another strength of FP is that it is much simpler to implement. On the other hand, ab initio MD explicitly includes temperature dependencies beyond the harmonic approximation, which is particularly important for systems with soft phonons that are dynamically unstable at 0 K.[68] Recent developments of machine learning force fields have enabled accurate MD simulations with reduced computational cost.[69] In what follows we focus on calculations within the framework of DFPT.

The justification for the applicability of the Bloch-Boltzmann formalism to the dynamics of electrons and phonons in materials rests on the concept of quasiparticles within Landau Fermi liquid theory [78, 79]. Electrons can be considered to be wavepackets that obey Newtonian laws of motion [80]. For a wavepacket with well defined momentum, $p=\hbar k$, and hence small uncertainty on the scale of the Fermi momentum, $\Delta k\ll k_{F}$, the uncertainty principle guarantees that $\Delta x\Delta(\hbar k)\sim\hbar\rightarrow\Delta x\sim\frac{1}{\Delta k}\gg\frac{1}{k_{F}}\sim a$, where $a$ is the lattice constant. Thus, the typical size of the electron wavepacket is much larger than the lattice constant, $\Delta x\gg a$. When this picture is valid, the electron distribution function can be defined over phase space cells larger than $\hbar^{3}$, the uncertainty principle is not violated, and the BTE can be justified. It assumes that classical external fields are varying sufficiently slowly in space and time to justify the use of electron wavepackets. This approach enables us to relate the transport properties to the Bloch band structure.

Purely quantum limitations have to be taken into account and the validity condition of the BTE can be related to the Peierls inequality [81] for non-degenerate semiconductors, $\frac{\hbar}{\tau}\ll k_{B}T$. This means that the time interval in which the distribution function evolves should be much smaller than the relaxation time. For metals, this condition can be smoothed by considering $\frac{\hbar}{\tau(\epsilon_{F})}\ll\epsilon_{F}$ [81], which can be related to the Mott-Ioffe-Regel criterion [82]

$$k_{F}l\gg 1~{},$$

(26)

where $k_{F}$ is the Fermi wave vector and $l$ is the mean free path of the carriers. When this criterion is not fulfilled, such as in strongly correlated bad metals [83], the BTE is no longer applicable since the system of interest may not support the existence of long-lived quasiparticles as described by Landau Fermi liquid theory [84]. Additionally, band transport based on the BTE can fail in describing the case of strong electron-phonon coupling leading to the small polaron limit, in which the charge transport is dominated by thermally activated polaron hopping.[85] Also, another subtle assumption that justifies the BTE is the random phase approximation. Within this approximation, the distribution function is only given by diagonal terms of the density matrix, while off-diagonal elements are neglected. Kohn and Luttinger provided a rigorous mathematical formalism to justify that off-diagonal elements are indeed suppressed by considering an ensemble average over a random distribution of impurities.[86, 87]

The Boltzmann Transport Equation describes the propagation of electron or phonon distribution functions, $f_{n{\bf{k}}}({\bf{r}},t)$ or $N_{\nu{\bf{q}}}({\bf{r}},t)$, respectively. The distribution functions give the probability that an electron (phonon) occupies a state with momentum ${\bf{k}}$ and band $n$ (momentum ${\bf{q}}$ and branch $\nu$) at position ${\bf{r}}$ and time $t$. Within this picture we assume that when the electron (phonon) BTE is solved, the phonon (electron) system remains in equilibrium, which relies on the Bloch assumption. Coupled electron-phonon transport that takes into account the drag effect of non-equilibrium electrons (phonons) on the phonons (electrons) involves the solution of a coupled electron–phonon BTE, which can be accomplished by a recently released solver called elphbolt.[88] In the examples that follow we assume that the drag effect can be ignored.

In the diffusive transport limit and the presence of a temperature gradient ($\nabla T$) and applied electric (${\bf{E}}$) and magnetic fields (${\bf{B}}$), carrier transport properties can be obtained by solving the semiclassical BTE for the nonequilibrium carrier distribution function $f_{n,{\bf{k}}}=f(\epsilon_{n,{\bf{k}}})$

$$\frac{\partial f_{n,{\bf{k}}}}{\partial t}+{\bf{v}}_{n,{\bf{k}}}\cdot\nabla_{{\bf{r}}}f_{n,{\bf{k}}}-\frac{{\bf{F}}}{\hbar}\cdot\nabla_{{\bf{k}}}f_{n,{\bf{k}}}=\left(\frac{\partial f_{n,{\bf{k}}}}{\partial t}\right)_{coll}~{}.$$

(27)

The electronic band velocity of the carrier in the state $\{n,{\bf{k}}\}$ with energy $\epsilon_{n,{\bf{k}}}$ is given by ${\bf{v}}_{n,{\bf{k}}}=\frac{1}{\hbar}\nabla_{{\bf{k}}}\epsilon_{n,{\bf{k}}}$, considering the diagonal matrix elements of the velocity operator and ignoring the Berry curvature.[89] The external force is given by ${\bf{F}}=e({\bf{E}}+{\bf{v}}_{n,{\bf{k}}}\times{\bf{B}})$, where $e$ is the absolute value of the charge of the carriers. In this equation, the temporal evolution of the electron distribution function results from a balance between drift and collision terms. The drift term on the left-hand side of Eq. (27) is an external field-driven flow of the space and momentum variables, while the collision term on the right-hand side accounts for any relevant scattering mechanisms within the diffusive or hydrodynamic regime. In Sect. IV.3 we will elaborate on the analogous BTE that can be written for the phonon distribution function, $N_{\nu{\bf{q}}}({\bf{r}},t)$, where the phonons follow the Bose-Einstein distribution at equilibrium and are not subject to the external force, $\bf{F}$. The electronic and phononic scattering mechanisms will be discussed in Sect. IV.

The knowledge of $f_{n,{\bf{k}}}$ allows for the evaluation of the charge current density,

$${\bf{j}}=-\frac{2e}{V}\sum_{n}\sum_{k}{\bf{v}}_{n,{\bf{k}}}f_{n,{\bf{k}}}=-\frac{2e}{(2\pi)^{3}}\sum_{n}\int{{\bf{v}}_{n,{\bf{k}}}f_{n,{\bf{k}}}d{\bf{k}}}~{},$$

(28)

and the heat energy flux density,

$${\bf{j}}_{Q}=\frac{2}{V}\sum_{n}\sum_{{\bf{k}}}\left(\epsilon_{n,{\bf{k}}}-\mu\right){\bf{v}}_{n,{\bf{k}}}f_{n,{\bf{k}}}\\ =\frac{2}{(2\pi)^{3}}\sum_{n}\int{\left(\epsilon_{n,{\bf{k}}}-\mu\right){\bf{v}}_{n,{\bf{k}}}f_{n,{\bf{k}}}d{\bf{k}}}~{},$$

(29)

where the factor of 2 appears due to the electron spin, $V$ is the crystal’s volume, and $\mu$ is the chemical potential.

In the next section we discuss the solution of the BTE in the relaxation time approximation for the case of charge carriers; the solution for phonons is analogous. The iterative approach to solving the BTE will be presented in Sect. III.4.

The collision term on the right-hand side of Eq. (27) drives the system towards a steady state. This scattering term can be expressed by introducing the per-unit-time probability, $W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})$, of the transition of the charge carrier from the state $\{n,{\bf{k}}\}$ to state $\{j,{\bf{{k}^{\prime}}}\}$, as a result of a particular scattering mechanism. From the principle of detailed balance, the number of charge carriers coming into the state $\{j,{\bf{{k}^{\prime}}}\}$ from $\{n,{\bf{k}}\}$ is the same as the number coming out from $\{j,{\bf{{k}^{\prime}}}\}$ into $\{n,{\bf{k}}\}$. That is, the scattering processes are equally likely going forward or in reverse. This allows us to write

$\displaystyle\left(\frac{\partial f_{n,{\bf{k}}}}{\partial t}\right)_{coll}=\sum_{j,{\bf{{k}^{\prime}}}}[W(j,{\bf{{k}^{\prime}}}|n,{\bf{k}})f_{j,{\bf{{k}^{\prime}}}}\left(1-f_{n,{\bf{k}}}\right)-W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})f_{n,{\bf{k}}}\left(1-f_{j,{\bf{{k}^{\prime}}}}\right)]~{}.$

(30)

This condition is a consequence of time-reversal symmetry of the microscopic equations of motion and is sufficient to ensure a positive local entropy production rate for systems out of equilibrium. It is described via the H-theorem[90] by writing the results as a function of H, the negative of the entropy.[91, 92] The intricate dependency of $f_{n,{\bf{k}}}$ on the linear response distribution function of all other states, $f_{j,{\bf{{k}^{\prime}}}}$, complicates the solution of the BTE and thus various forms of the RTA are usually applied. Poncé et.al. [76] reviewed different levels of approximations such as the momentum relaxation time approximation (MRTA), the self-energy relaxation time approximation (SERTA)[93] and the lowest-order variational approximation (LOVA)[94] or the Ziman resistivity formula for metals.[95]. In general, both the computational cost and the accuracy decreases from the former to the latter. Here we focus on SERTA, which includes forward scattering from $\{n,{\bf{k}}\}$ into $\{j,{\bf{{k}^{\prime}}}\}$ and neglects backward scattering rates into the state $\{n,{\bf{k}}\}$. MRTA partially includes backward scattering rates by considering a geometrical factor that depends on the scattering angle.[76] On the other hand, LOVA includes only an average of the state- and momentum-resolved total decay rates combined with Drude’s formula.[96]

Assuming that the system is close enough to local equilibrium, the non equilibrium distribution function, $f_{n,{\bf{k}}}$, differs only slightly from that of the equilibrium state, $f_{n,{\bf{k}}}^{(0)}$; namely $\Delta f(n,{\bf{k}})=\lvert f_{n,{\bf{k}}}-f_{n,{\bf{k}}}^{(0)}\rvert\ll f_{n,{\bf{k}}}^{(0)}$. Consequently, $f_{n,{\bf{k}}}$ can be expanded to first order as

$$f_{n,{\bf{k}}}=f_{n,{\bf{k}}}^{(0)}-{\tau}_{n,{\bf{k}}}{\bf{v}}_{n,{\bf{k}}}\cdot{\bf{\Phi_{0}}}(\epsilon)\left(\frac{\partial f^{(0)}}{\partial\epsilon}\right)~{},$$

(31)

where ${\bf{\Phi_{0}}}(\epsilon)=-e{{\varepsilon}}-\frac{\epsilon-\mu}{T}\nabla T$ is the generalized disturbing force (dynamic and static) causing the deviation from the equilibrium distribution, and ${{\varepsilon}}={\bf{E}}+(1/e)\nabla\mu=-\nabla(\phi_{0}-(\mu/e))$ is the gradient of the electrochemical potential. Using this approximation, Eq. (30) can be written in the RTA as

$$\left(\frac{\partial f_{n,{\bf{k}}}}{\partial t}\right)_{coll}=-\frac{\Delta f(n,{\bf{k}})}{\tau_{n,{\bf{k}}}}~{},$$

(32)

where

$$\frac{1}{{\tau_{n,{\bf{k}}}}}=\sum_{{\bf{{k}^{\prime}}}}\sum_{j}W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})\\ \left(\frac{1-f^{(0)}_{j,{\bf{{k}^{\prime}}}}}{1-f^{(0)}_{n,{\bf{k}}}}-\frac{f^{(0)}_{n,{\bf{k}}}}{f^{(0)}_{j,{\bf{{k}^{\prime}}}}}\frac{\Delta f(j,{\bf{{k}^{\prime}}})}{\Delta f(n,{\bf{k}})}\right)~{},$$

(33)

considering both the absence of quantization effects and that $W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})$ does not depend on ${\bf{E}}$, ${\bf{B}}$, or $\nabla T$. Substituting Eq. (31) into Eq. (33) we obtain

$$\frac{1}{\tau_{n,{\bf{k}}}}=\sum_{{\bf{{k}^{\prime}}}}\sum_{j}W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})\\ \frac{1-f^{(0)}(\epsilon^{\prime})}{1-f^{(0)}(\epsilon)}\left(1-\frac{{\tau}_{j,{\bf{{k}^{\prime}}}}}{\tau_{n,{\bf{k}}}}\frac{{\bf{v}}_{j,{\bf{{k}^{\prime}}}}\cdot\bf{\Phi_{0}}(\epsilon^{\prime})}{{\bf{v}}_{n,{\bf{k}}}\cdot\bf{\Phi_{0}}(\epsilon)}\right)~{}.$$

(34)

Thus, in the steady-state limit of a homogeneous system with no magnetic field, Eq. (27) simplifies to

$${\bf{v}}_{n,{\bf{k}}}\cdot\nabla_{{\bf{r}}}f_{n,{\bf{k}}}-\frac{e{\bf{E}}}{\hbar}\cdot\nabla_{{\bf{k}}}f_{n,{\bf{k}}}=-\frac{\Delta f(n,{\bf{k}})}{\tau_{n,{\bf{k}}}}~{},$$

(35)

from which the non equilibrium distribution function is obtained provided that $\tau_{n,{\bf{k}}}$ does not depend on ${\bf{E}}$ or $\nabla T$.

A common additional approximation is the assumption that $W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})$ does not depend on ${\bf{k}}$ and ${\bf{{k}^{\prime}}}$ separately, but only on their magnitudes and the angle between them, $W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})=W_{n,j}(\lvert{\bf{k}}\rvert,\lvert{\bf{{k}^{\prime}}}\rvert,{\bf{k}}\cdot{\bf{{k}^{\prime}}})$. Then, Eq. (34) may be rewritten as [97]

$$\frac{1}{{\tau_{n,{k}}}}=\sum_{{\bf{{k}^{\prime}}}}\sum_{j}W(n,{\bf{k}}|j,{\bf{{k}^{\prime}}})\left(1-\frac{{\bf{k}}\cdot{\bf{{k}^{\prime}}}}{k^{2}}\right)~{},$$

(36)

assuming that charge carrier scattering is purely elastic and the dispersion relation is an arbitrary spherically symmetric function of the magnitude of the wavevector (not necessarily parabolic), so that $\epsilon(\lvert{\bf{k}}\rvert)=\epsilon(\lvert{\bf{{k}^{\prime}}}\rvert)$. Although Eq. (36) was derived for isotropic bands, results obtained from it have been used to study transport properties of chalcogenides, which are anisotropic [98, 99, 100]. This is possible because the transport properties along the different directions in these materials are mutually independent, except for the magnetoresistance, which critically depends on the anisotropy.

In some cases it is possible to go beyond the RTA approach to the BTE by using an iterative solution method. From perturbation theory, if the interaction potential is weak enough, the scattering can be treated in the Born approximation and Fermi’s golden rule can be used to determine the transition probability of the electronic scattering process $\{\ket{n{\bf{k}}}\}\rightarrow\{\ket{m{\bf{k}}+{\bf{q}}}\}$, describing absorption of a phonon with mode $\{\nu\}$ and wave vector ${\bf{q}}$:

$$W_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}=\frac{2\pi}{\hbar}|g_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}|^{2}\delta(\epsilon_{n{\bf{k}}}+\hbar\omega_{\nu{\bf{q}}}-\epsilon_{m{\bf{k}}+{\bf{q}}})~{}.$$

(37)

The transition probability exhibits the microreversibility property which stems from the time-reversal invariance of the microscopic equations of motion. Thus the last equation is equal to the reverse transition $W^{n{\bf{k}},\nu{\bf{q}}}_{m{\bf{k}}+{\bf{q}}}$ from $\{\ket{m{\bf{k}}+{\bf{q}}}\}$ to $\{\ket{n{\bf{k}}}\}$ by emitting a phonon. The transition rate at equilibrium, which is the transition per unit time is given by

$$\Pi_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}=\Pi^{n{\bf{k}},\nu{\bf{q}}}_{m{\bf{k}}+{\bf{q}}}=f_{n{\bf{k}}}^{0}(1-f_{m{\bf{k}}+{\bf{q}}}^{0})N_{\nu{\bf{q}}}^{0}W_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}~{}.$$

(38)

By writing the canonical form of the scattering term [95] of the BTE for the case of el-ph coupling

$$\frac{\partial f_{n{\bf{k}}}}{\partial t}\bigg{|}_{scatt}=-\sum_{\nu{\bf{q}}}\left(\Pi_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}+\Pi_{n{\bf{k}}}^{m{\bf{k}}+{\bf{q}},-\nu{\bf{q}}}\right)\left(\chi_{n{\bf{k}}}-\chi_{m{\bf{k}}+{\bf{q}}}\right)$$

(39)

where $\chi_{n{\bf{k}}}=\frac{f_{n{\bf{k}}}-f_{n{\bf{k}}}^{0}}{f_{n{\bf{k}}}^{0}(1-f_{n{\bf{k}}}^{0})}=\frac{q{\bf{E}}}{k_{B}T}\cdot{\bf{G}}_{n{\bf{k}}}$ and

$$\Pi_{n{\bf{k}}}^{m{\bf{k}}+{\bf{q}},-\nu{\bf{q}}}=\frac{2\pi}{\hbar}|g_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}|^{2}f_{n{\bf{k}}}^{0}(1-f_{m{\bf{k}}+{\bf{q}}}^{0})(1+N^{0}_{-\nu{\bf{q}}})\delta(\epsilon_{n{\bf{k}}}-\hbar\omega_{-\nu{\bf{q}}}-\epsilon_{m{\bf{k}}+{\bf{q}}})~{}.$$

(40)

In the steady-state and approximating the drift term by keeping only the linear terms in ${\bf{E}}$, the BTE can be linearized as [101]

$${\bf{v}}_{n{\bf{k}}}f_{n{\bf{k}}}^{0}(1-f_{n{\bf{k}}}^{0})+\sum_{m,\nu{\bf{q}}}\left(\Pi_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}+\Pi_{n{\bf{k}}}^{m{\bf{k}}+{\bf{q}},-\nu{\bf{q}}}\right){\bf{G}}_{m{\bf{k}}+{\bf{q}}}={\bf{G}}_{n{\bf{k}}}\sum_{m,\nu{\bf{q}}}\left(\Pi_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}+\Pi_{n{\bf{k}}}^{m{\bf{k}}+{\bf{q}},-\nu{\bf{q}}}\right)~{}.$$

(41)

The relaxation time approximation is reached when the left hand side summation is neglected and then ${\bf{G}}_{n{\bf{k}}}^{\tau}={\bf{v}}_{n{\bf{k}}}\cdot\tau_{n{\bf{k}}}$, where $\tau_{n{\bf{k}}}$ is given by Eq. (53). Alternatively, Eq. (41) can be solved by starting with the RTA solution as an initial guess, ${\bf{G}}_{n{\bf{k}}}^{0}={\bf{G}}_{n{\bf{k}}}^{\tau}$, and iterating using the following relation:

$${\bf{G}}_{n{\bf{k}}}^{i+1}={\bf{G}}_{n{\bf{k}}}^{0}+\frac{\tau_{n{\bf{k}}}}{f_{n{\bf{k}}}^{0}(1-f_{n{\bf{k}}}^{0})}\left(\Pi_{n{\bf{k}},\nu{\bf{q}}}^{m{\bf{k}}+{\bf{q}}}+\Pi_{n{\bf{k}}}^{m{\bf{k}}+{\bf{q}},-\nu{\bf{q}}}\right){\bf{G}}^{i}_{m{\bf{k}}+{\bf{q}}}~{}.$$

(42)

This iterative approach has been implemented in EPW [102] for electron transport, ShengBTE [103] for thermal transport, and for both charge carriers and phonons in Perturbo[104] and Phoebe[105].

The tensorial formalism based on the Onsager-de Groot-Callen model[106, 107, 108, 109, 110] is appropriate to discuss thermoelectric (TE) effects for anisotropic materials using the thermodynamics of irreversible processes and linear response theory. In this model, the flux of charge carriers and thermal energy is described in terms of a kinetic matrix and generalized forces. (For reviews see Goupil [111] and Feldhoff.[112]) In summary, the off-diagonal coupling between the electronic current density, ${\bf{j}}$, and heat energy flux density, ${\bf{j}}_{Q}$, can be written:

$$\begin{bmatrix}{\bf{j}}\\ {\bf{j}}_{Q}\end{bmatrix}=\begin{bmatrix}{\bf{L^{11}}}&{\bf{L^{12}}}\\ {\bf{L^{21}}}&{\bf{L^{22}}}\end{bmatrix}\cdot\begin{bmatrix}{{\varepsilon}}\\ -\frac{\nabla T}{T}\end{bmatrix}$$

(43)

in which, ${\bf{L^{11}}}$, ${\bf{L^{12}}}$, ${\bf{L^{21}}}$, ${\bf{L^{22}}}$ are the moments of the generalized transport coefficients. The Onsager reciprocity relations guarantee ${\bf{L^{12}}}={\bf{L^{21}}}$  [106, 107, 110]. These kinetic coefficients are defined by

$$\Lambda^{(\alpha)}(\mu;T)=e{{}^{2}}\int\Xi(\epsilon,\mu,T)(\epsilon-\mu)^{\alpha}\left(-\frac{\partial f^{(0)}(\mu;\epsilon,T)}{\partial\epsilon}\right)d\epsilon~{},$$

(44)

with ${\bf{L^{11}}}=\Lambda^{(0)}$, ${\bf{L^{21}}}={\bf{L^{12}}}=-(1/e)\Lambda^{(1)}$, and ${\bf{L^{22}}}=(1/e^{2})\Lambda^{(2)}$, in which $\Xi(\epsilon,\mu,T)$ is the transport distribution kernel (TDK) given by

$$\Xi(\epsilon,\mu,T)=\int\sum_{n}{{\bf{v}}_{n,{\bf{k}}}\otimes{\bf{v}}_{n,{\bf{k}}}{\tau}_{n,{k}}}(\mu,T)\delta(\epsilon-\epsilon_{n,{\bf{k}}})\frac{d{\bf{k}}}{8\pi^{3}}~{}.$$

(45)

In the dynamic steady state, thermoelectric properties can be obtained by considering specific experimental conditions. For an isothermal situation ($\nabla T=0$), the charge current, ${\bf{j}}$, obeys Ohm’s law and the kinetic coefficient tensor can be identified with the electrical conductivity tensor, $\sigma=\Lambda^{(0)}$. As pointed out by Feldhoff [112], even with the isothermal condition, ${\bf{j}}$ is accompanied by an entropy current whose magnitude and direction depend on the Seebeck coefficient, $S$, which can be viewed as a quantity that measures entropy flow per unit charge. $S$ is obtained by requiring the electric current to vanish, so the electrochemical potential gradient and the thermal gradient are balanced. Then $S$ is given by the ratio between them, $S=(eT)^{-1}\Lambda^{(1)}/\Lambda^{(0)}$.

As one of the most sensitive probes of the carriers in a material, the Seebeck measurement is related to the heat per carrier over temperature or the entropy per carrier, as pointed out earlier. It suggests that $S$ can provide information about $i)$ the sign of the charge of the carriers and $ii)$ the characteristic energy associated with the carriers. However, the former result does not hold in a few cases, such as noble metals, where $S$ and the Hall coefficient have diverging signs.[113] This has been attributed to the complex energy dependence of the mean free path due to the electron-phonon scattering.[114]. The second result generally holds and can be used, for example, to distinguish metals and semiconductors by comparing the magnitude and temperature behavior of $S$. While metals exhibit values of $S$ that decrease with temperature and have magnitudes much smaller than $k_{B}/e\approx 87\mu V/K$, semiconductors present much larger magnitudes that increase with temperature.[113]

The zero electric current condition yields Fourier’s law, in which the charge carrier contribution to the thermal conductivity tensor is given by

$$\kappa_{el}=(e^{2}T)^{-1}\left({\Lambda^{(1)}\cdot{\Lambda^{(0)}}^{-1}}\cdot{\Lambda^{(1)}}-\Lambda^{(2)}\right)~{}.$$

(46)

The second term within brackets corresponds to the thermal conductivity due to the carrier transport under isoelectrochemical conditions, while the first term is the power factor ($PF=\sigma S^{2}$) related to the thermal conductivity that couples to the charge current. This is responsible for thermoelectric conversion, that is, the transfer of energy from an entropy current to an electric current or vice versa.[115, 112]

The electron–phonon interaction is a key factor in determining functional properties of materials such as thermoelectric transport properties. The electron-phonon coupling hamiltonian can be derived within DFT by expanding the KS Hamiltonian to first-order in $\Delta\xi_{\kappa p}$, the displacements of the nuclei from their equilibrium positions, $\xi^{0}_{\kappa p}$:

$$V^{KS}(\{\xi_{\kappa p}\})=V^{KS}(\{\xi^{0}_{\kappa p}\})+\sum_{\kappa\alpha p}\frac{\partial V^{KS}}{\partial\xi_{\kappa\alpha p}}\Delta\xi_{\kappa\alpha p}~{}.$$

(47)

Extending the expansion to second-order results in Debye-Waller terms [116] that will not be discussed in the present review, since these terms are purely real and do not affect el-ph relaxation times. [117] The operator for the phonon perturbation potential is given in terms of quantized normal mode coordinates using phonon annihilation and creation operators ($b_{\nu\mathbf{q}},b^{\dagger}_{\nu\mathbf{q}}$):

$$V^{KS}=V^{KS}(\{\xi^{0}_{\kappa p}\})+\frac{1}{N_{p}^{1/2}}\sum_{\nu{\bf{q}}}\Delta_{\nu{\bf{q}}}V^{KS}(\hat{b}_{\nu{\bf{q}}}+\hat{b}^{\dagger}_{-\nu{\bf{q}}})~{},$$

(48)

where

$$\Delta_{\nu{\bf{q}}}V^{KS}=e^{i{\bf{q}}\cdot{\bf{r}}}\Delta_{\nu{\bf{q}}}v^{KS}~{},$$

(49)

and

$$\Delta_{\nu{\bf{q}}}v^{KS}=c_{\nu{\bf{q}}}\sum_{\kappa\alpha p}\left(\frac{M_{0}}{M_{\kappa}}\right)^{\frac{1}{2}}e_{\kappa\alpha,\nu}({\bf{q}})e^{-i{\bf{q}}\cdot({\bf{r}}-{\bf{R}}_{p})}\frac{\partial V^{KS}}{\partial\xi_{\kappa\alpha}}\bigg{|}_{{\bf{r}}-{\bf{R}}_{p}}~{},$$

(50)

with $c_{\nu{\bf{q}}}$ being the zero-point displacement amplitude. Using the electron annihilation and creation operators, ($a_{n\mathrm{k}},a^{\dagger}_{n\mathrm{k}}$), the first-principles el-ph perturbation Hamiltonian is

	$\displaystyle\hat{H}^{el-ph}$	$\displaystyle=$	$\displaystyle\sum_{n{\bf{k}},n^{\prime}{\bf{k^{\prime}}}}\braket{\phi_{n^{\prime}{\bf{k^{\prime}}}}}{V^{KS}(\{\xi_{\kappa p}\})-V^{KS}(\{\xi^{0}_{\kappa p}\})}{\phi_{n{\bf{k}}}}\hat{b}_{\nu{\bf{q}}}\hat{a}^{\dagger}_{n{\bf{k}}}\hat{a}_{n^{\prime}{\bf{k^{\prime}}}}$		(51)
		$\displaystyle=$	$\displaystyle\frac{1}{N_{p}^{1/2}}\sum_{{\bf{k}},{\bf{q}}mn\nu}g_{mn\nu}({\bf{k}},{\bf{q}})\hat{a}^{\dagger}_{n{\bf{k}}}\hat{a}_{n^{\prime}{\bf{k^{\prime}}}}(\hat{b}_{\nu{\bf{q}}}+\hat{b}^{\dagger}_{-\nu{\bf{q}}})~{},$

where the electron-phonon matrix element has been defined as

$$g_{mn\nu}({\bf{k}},{\bf{q}})=\braket{\psi_{m\mathbf{k}+\mathbf{q}}}{\Delta_{\nu{\bf{q}}}v^{KS}}{\psi_{n{\bf{k}}}}~{},$$

(52)

which can be computed by DFPT (see Sect. II.2), using Eq. (22) [118]. From the electron-phonon Hamiltonian, the el-ph relaxation times (RT) can be derived in different ways, such as perturbatively via Feynman-Dyson diagram techniques [119, 120] or nonperturbatively via Hedin-Baym equations within quantum field theory [121, 122, 118]. Also, within the Born approximation, the RT can be derived by considering all scatterings in which an electron emits or absorbs one phonon via Fermi’s golden rule, along with the rate equation for the time-dependent electron distribution functions [96]. The latter approach was derived in Section III.4. In all cases the RT is given by the imaginary part of the electron self energy, $\Sigma$:

$$\begin{split}\operatorname{Im}\left[\Sigma_{n,{\bf{k}}}(\epsilon,T)\right]=\pi\sum_{m,\nu}\int_{BZ}\frac{d{\bf{q}}}{\Omega_{BZ}}|g_{mn\nu}({\bf{k},{\bf{q}}})|^{2}\\ \times\Bigg{[}\left[N_{\nu{\bf{q}}}(T)+f_{m{\bf{k}}+{\bf{q}}}\right]\delta(\epsilon-(\epsilon_{m{\bf{k}}+{\bf{q}}}-\epsilon_{F})+\hbar\omega_{{\bf{q}}\nu})\\ +[N_{\nu{\bf{q}}}(T)+1-f_{m{\bf{k}}+{\bf{q}}}]\delta(\epsilon-(\epsilon_{m{\bf{k}}+{\bf{q}}}-\epsilon_{F})-\hbar\omega_{{\bf{q}}\nu})\Bigg{]}~{},\end{split}$$

(53)

where $g_{mn\nu}({\bf{k},{\bf{q}}})$ are the el-ph matrix elements, $\Omega_{BZ}$ is the volume of the Brillouin zone (BZ) and the Dirac $\delta$-functions enforce energy conservation for emission or absorption of a phonon with wavevector ${\bf{q}}$ and mode $\nu$ and energy $\hbar\omega_{\nu{\bf{q}}}$. Equation (53) contains the dynamical structure of the electron-phonon interaction on the scale of the phonon energy. The temperature dependence comes from the electron and phonon distribution functions. Eq. (53) corresponds to the imaginary part of the lowest-order Feynman diagram for the electron self-energy shown in Fig. 2(c). Migdal theory demonstrated that non-adiabatic terms are less important in higher order diagrams, justifying the truncation at the lowest order diagram [123, 124].

Since the RTs are resolved for different bands and ${\bf{k}}$-points, after integrating over all phonon modes and ${\bf{q}}$-points, first-principles calculations can provide rich microscopic information. In order to converge the el-ph RTs when calculating transport properties, very dense ${\bf{k}}$ and ${\bf{q}}$ meshes are needed [101, 125], so interpolation schemes have to be used instead brute-force DFPT calculations. In Sect. V.1 we discuss interpolation schemes developed for this aim.