Heat transport in silicon from first-principles calculations

Using harmonic and anharmonic force constants extracted from density functional calculations within a supercell, we have developed a relatively simple but general method to compute thermodynamic and thermal properties of any crystal. First, from the harmonic, cubic, and quartic force constants, we construct a force field for molecular dynamics. It is exact in the limit of small atomic displacements and thus does not suffer from inaccuracies inherent in semiempirical potentials such as Stillinger-Weber's. By using the Green-Kubo formula and molecular dynamics simulations, we extract the bulk thermal conductivity. This method is accurate at high temperatures where three-phonon processes need to be included to higher orders, but may suffer from size scaling issues. Next, we use perturbation theory (Fermi golden rule) to extract the phonon lifetimes and compute the thermal conductivity $\kappa$ from the relaxation-time approximation. This method is valid at most temperatures, but will overestimate $\kappa$ at very high temperatures, where higher-order processes neglected in our calculations also contribute. As a test, these methods are applied to bulk crystalline silicon, and the results are compared and differences are discussed in more detail. The presented methodology paves the way for a systematic approach to model heat transport in solids using multiscale modeling, in which the relaxation time due to anharmonic three-phonon processes is calculated quantitatively, in addition to the usual harmonic properties such as phonon frequencies and group velocities. It also allows the construction of an accurate bulk interatomic potentials database.

Figure 1

Total energy as an atom is moved in the (a) [100], (b) [110], and (c) [111] directions. DFT results are compared with the force field and the Stillinger-Weber potential. MD234 refers to the force field in which all harmonic, cubic, and quartic terms are included, while MD2 refers only to the harmonic force field, etc.

Figure 2

Plot of the ensemble-averaged (over 27 initial conditions) heat-current autocorrelation as a function of time, and its integral for the $10\times 10\times 10$ supercell. Vertical units for the integrated autocorrelation are in W/mK, and the autocorrelation (blue dots) has been multiplied by a constant to be on scale.

Figure 3

Phonon band structure of Si using force constants up to the fifth neighbor shell. Plus signs represent experimental data of Nelin and Nilsson.[30] Left portion of figure shows the DOS, and right portion of figure shows the rescaled Gruneisen parameters [$100\times (\gamma -3)$].

Figure 4

Top line in blue is the DOS associated with two-phonon creation or annihilation (${\mathrm{DOS}}_2^{-}$), and bottom line in green is the DOS associated with one-phonon emission or absorption (${\mathrm{DOS}}_2^{+}$). The red line is the contribution of the matrix elements defined in Eq. (19). The peak at 500 cm${}^{-1}$ is the main reason for smaller lifetimes of optical modes.

Figure 5

Lifetimes of the six branches in Si at 277 K vs. frequency on a logarithmic scale. Top plot is for umklapp and bottom plot is for normal processes. The quadratic dependence of the acoustic modes can be noticed for normal processes, while umklapp processes seems to scale as $1/{\omega }^3$.

Figure 6

Cumulative contributions of phonons to the thermal conductivity at 277 K from the $18\times 18\times 18$ $k$-point mesh data. Left plot is according to the wavelengths, and right plot is according to the MFPs. Both differential and cumulative thermal conductivities are shown in blue and red, respectively. For comparison, the extrapolated (to infinite $k$-point mesh) and experimental $\kappa$ are also shown with horizontal lines at 166 and 174 W/mK, respectively.

Figure 7

Thermal conductivity of pure Si crystal from Eq. (16) vs. temperature. Inset: Extrapolation to infinitely dense $k$-points.

Figure 8

Quantum and classical thermal conductivities of Si vs. temperature [from Eq. (16) with $k$-point mesh $=$ 18]. The contribution of each mode is also being plotted.