---
abstract: 'We herein present a [*first-principles*]{} formulation of the [*Green-Kubo*]{} method that allows the accurate assessment of the phonon thermal conductivity of solid semiconductors and insulators in equilibrium [*ab initio molecular dynamics*]{} calculations. Using the virial for the nuclei, we propose a unique [*ab initio*]{} definition of the heat flux. Accurate size- and time convergence are achieved within moderate computational effort by a robust, asymptotically exact extrapolation scheme. We demonstrate the capabilities of the technique by investigating the thermal conductivity of extreme high and low heat conducting materials, namely diamond Si and tetragonal ZrO$_2$.'
author:
- Christian Carbogno
- Rampi Ramprasad
- Matthias Scheffler
bibliography:
- 'bib2.bib'
title: |
    *Ab initio* Green-Kubo Approach\
    for the Thermal Conductivity of Solids
---

Macroscopic heat transport is an ubiquitous phenomenon in condensed matter that plays a crucial role in a multitude of applications, e.g., energy conversion, catalysis, and turbine technology. Whenever a temperature gradient $\nabla T({\bm{R}})$ is present, a heat flux ${\bm{J}}({\bm{R}})$ spontaneously develops to move the system back toward thermodynamic equilibrium. The temperature- and pressure-dependent thermal conductivity ${\bm{\kappa}}(T,p)$ of the material describes the proportionality between heat flux and temperature gradient (Fourier’s law) $${\bm{J}}({\bm{R}}) = - {\bm{\kappa}}(T,p)\cdot\nabla T({\bm{R}}) \;.
\label{Fourier}$$ In insulators and semiconductors, the dominant contribution to ${\bm{\kappa}}(T,p)$ stems from the vibrational motion of the atoms (phonons) [@AshcroftMermin]. In spite of significant efforts, a parameter-free, accurate theoretical approach that allows to assess the thermal conductivity tensor both in the case of weak and strong anharmonicity is still lacking: Studies of model systems via classical molecular dynamics (MD) based on force fields (FF) can unveil general rules and concepts [@Donadio:2009fo]. However, the needed accuracy for describing anharmonic effects is often not correctly captured by FFs [@Broido:2005kf] and trustworthy FFs are generally not available for “real” materials used in scientific and industrial applications.

Naturally, first-principles electronic-structure theory lends itself to overcome this deficiency by allowing a reliable computation of the inter-atomic interactions. However, severe limitations affect the approaches that have hitherto been employed in [*ab initio*]{} frameworks for studying the thermal conductivity of solids: (a) Approaches based on the [ *Boltzmann Transport Equation*]{} [@Broido:2007ua; @deKoker:2009cr; @Tang:2010vu; @Garg:2011hi] account for the leading, lowest order contributions to the anharmonicity. Accordingly, these approaches are justified at low temperatures, at which they also correctly describe relevant nuclear quantum effects. At elevated temperatures and/or in the case of strong anharmonicity, this approximation is however known to break down [@Ladd:1986tv; @Turney:2009ur]. (b) [*Non-Equilibrium*]{} approaches [@Gibbons:2009hd; @Gibbons:2011je; @Stackhouse:2010vg] require to impose an artificial temperature gradient, which becomes unreasonably large ($\gg 10^9$ K/m) in the limited system sizes accessible in [*first-principles*]{} calculations. Especially at high temperatures, this can lead to non-linear artifacts [@Tenenbaum:1982wu; @Schelling:2002jl; @He:2012hg] that prevent the assessment of the linear response regime described by Fourier’s law.

In this letter, we present an [*ab initio*]{} implementation of the [*Green-Kubo*]{} (GK) method [@Kubo:1957tf], which does not suffer from the aforementioned limitations [@Schelling:2002jl; @Broido:2007ua], since ${\bm{\kappa}}(T,p)$ is determined from [*ab initio*]{} molecular dynamics simulations (aiMD) in thermodynamic equilibrium that account for anharmonicity to all orders. Hitherto, fundamental challenges have prevented an application of this technique in a first-principles framework: Conceptually, a definition of the heat flux associated with vibrations in the solid is required; numerically, the necessary time- and length scales need to be reached. First, we succinctly describe how we overcome the conceptual hurdles, i.e., the unique [*ab initio*]{} definition of the microscopic heat flux (and its fluctuations) for solids. Second, we discuss how this allows to overcome the numerical hurdles by introducing a robust extrapolation scheme, so that time- and size convergence is achieved within moderate computational effort. Third, we validate our formalism and demonstrate its wide applicability by investigating the thermal conductivity of diamond Si and tetragonal ZrO$_2$ ($\mathit{P4_2/nmc}$), two materials that feature especially large/low thermal conductivities due to being particularly harmonic/anharmonic.

For a given pressure $p$, volume $V$, and temperature $T$, the fluctuation-dissipation theorem, which is the only assumption entering the GK formula, relates the cartesian components $\alpha\beta$ of the thermal conductivity tensor $${\kappa}_{\alpha\beta}(T,p) = \frac{V}{{k_{\mbox{\tiny B}}}T^2}\lim_{\tau\rightarrow\infty}  \int\limits_0^{\tau} \left\langle {G}[{\bm{J}}]_{\alpha\beta}(\tau') \right\rangle_{(T,p)} \;  d\tau' 
\label{GreenKubo}$$ to the time-(auto)correlation functions $${G}[{\bm{J}}]_{\alpha\beta}(\tau) =  \lim_{t_0\rightarrow\infty} \frac{1}{t_0}\int\limits_0^{t_0-\tau} {J}_\alpha(t)\; {J}_\beta(t+\tau)\; dt 
\label{AC}$$ of the heat flux ${\bm{J}}(t)$. In Eq. (\[GreenKubo\]), ${k_{\mbox{\tiny B}}}$ is the Boltzmann constant and $\left\langle\cdot\right\rangle_{(T,p)}$ denotes an ensemble average that is performed by averaging over multiple correlation functions ${G}[{\bm{J}}]_{\alpha\beta}(\tau)$, which are individually computed from different MD trajectories (time span $t_0$, microcanonical ensemble [@MACDOWELL:1999fe]) using Eq. (\[AC\]).

First, the GK method requires a consistent definition of the heat flux ${\bm{J}}(t)$. Common FF-based formalisms [@MACDOWELL:1999fe; @Hardy:1963gh; @Helfand:1960uy] achieve such a definition by partitioning the total, i.e., kinetic and potential, energy of the system $E      = \sum_{I} E_{I}$ into contributions $E_I$ associated with the individual atoms $I$. Using their positions ${\bm{R}}_I$, the energy density associated with the nuclei is $e({\bm{R}}) = \sum_I E_I \delta({\bm{R}}-{\bm{R}}_I)$ with the Delta distribution $\delta({\bm{R}})$. With this subdivision, the integration of the continuity equation $\partial e({\bm{R}},t)/\partial t + \nabla \cdot {\bm{j}}({\bm{R}},t) = 0$ for the heat flux density ${\bm{j}}({\bm{R}},t)$ reveals that the total heat flux $${\bm{J}}(t)=\frac{1}{V} \frac{d}{dt}\sum_{I} {\bm{R}}_I E_I
\label{Bary}$$ is related to the motion of the energy barycenter. Conceptually, the required partitioning is straightforward for FFs and challenging in a first-principles framework [@Yu:2011kv; @Marcolongo:2016dn], but in neither of the cases unique [@Howell:2012fh; @Marcolongo:2016dn]. Using a combined nuclear and electronic energy density, Marcolongo, Umari, and Baroni recently proposed a non-unique formulation of the heat flux and used it to study the convective heat flux in liquids from first principles [@Marcolongo:2016dn]. However, their approach is numerically unsuitable for the conductive thermal transport in solids, which features much longer lifetimes and mean free paths. We overcome this limitation by disentangling the different contributions to the heat flux and thus finding a unique definition for the conductive heat flux in solids. In turn, this allows to establish a link to the quasi-particle (phonon) picture of heat transport and thus to overcome finite time and size effects, as described in the second part of this letter.

For this purpose, we perform the time derivative in Eq. (\[Bary\]) analytically [@MACDOWELL:1999fe; @Hardy:1963gh] $$\begin{aligned}
{\bm{J}}(t) & = & \underbrace{\frac{1}{V}\sum_{I} \dot{{\bm{R}}}_I E_I}_{{\bm{J}}_c(t)} + \underbrace{\frac{1}{V}\sum_{I} {\bm{R}}_I \dot{E_I}}_{{\bm{J}}_v(t)} \;.
\label{disentangle}\end{aligned}$$ The first term ${\bm{J}}_c(t)$, which describes convective contributions to the heat flux, requires an energy partitioning scheme, but gives no contributions to the conductivity tensor in solids [@Howell:2012fh], as mass transport is negligible. Conversely, the dominant virial or conductive term ${\bm{J}}_v(t)$ does not require an [*ad hoc*]{} partitioning of the energy $E$. Only derivatives of the energy $E_I$ enter ${\bm{J}}_v(t)$, so that the forces ${\bm{F}}_I=-{\partial U}/{\partial {\bm{R}}_I}$, i.e., the gradients of the potential energy surface $U$, naturally disentangle the individual atomic contributions in a unique fashion – both in FF and *ab initio* frameworks. By rewriting the individual contributions in terms of relative distances ${\bm{R}}_{IJ} = {\bm{R}}_I - {\bm{R}}_J$ we get a definition of the virial flux that is compatible with periodic boundary conditions $$\begin{aligned}
{\bm{J}}_v(t) & = &  -\frac{1}{V}   \sum_{I,J} ({\bm{R}}_I-{\bm{R}}_J) (\nabla_{{\bm{R}}_{IJ}}  U)      \cdot \dot{{\bm{{R}}}}_I \nonumber\\
             & = &   \sum_{I} {\bm{\sigma}}_I \cdot \dot{{\bm{R}}}_I \;.
\label{HFV}\end{aligned}$$ As discussed in the detailed derivation provided in the Supp. Mat., the latter notation highlights that the terms in ${\bm{J}}_v(t)$ are the individual atomic contributions ${\bm{\sigma}}_I$ to the stress strensor ${\bm{\sigma}} = \sum_{I} {\bm{\sigma_I}}$.

In [*density-functional theory*]{} (DFT), the potential-energy surface $$U = E^{\text{DFT}} + \frac{1}{2}\sum_{I,J\neq I} \frac{Z_IZ_J}{|{\bm{R}}_I-{\bm{R}}_J|}$$ is given by the total energy $E^{\text{DFT}}$ of the electrons with their ground-state density $n({\bm{r}})$ plus the electrostatic repulsion between the nuclei with charges $Z_I$. Use of the [Hellman-Feynman]{} theorem leads to a definition for the cartesian components $\alpha\beta$ of the virials for the individual nuclei $$\begin{aligned}
{\sigma}_I^{\alpha\beta} & = &  - \frac{Z_I}{V}\left(\int d{\bm{r}}\;n({\bm{r}}) \frac{({r}^\alpha-{R}_I^\alpha)({r}^\beta-{R}_I^\beta)}{|{\bm{r}}-{\bm{R}}_I|^3} \right. \nonumber\\
 && \left. - \frac{1}{2}\sum_{J\neq I} Z_J \frac{({R}_J^\alpha-{R}_I^\alpha)({R}_J^\beta-{R}_I^\beta)}{|{\bm{R}}_J-{\bm{R}}_I|^3} \right)\label{sigma}\;,\end{aligned}$$ whereby all electronic contributions stem from the interaction with the ground state electron density $n({\bm{r}})$. In turn, this enables a straightforward and unique evaluation of ${\bm{J}}_v(t)$ using Eq. (\[HFV\]), since neglecting the convective term ${\bm{J}}_c(t)$ from the beginning allows to integrate out the internal electronic contributions to the heat flux (see Supp. Mat.). This holds true also in our practical implementation of the virial and the analytical stress tensor [@Knuth:2015kc], since Pulay terms and alike that can arise can again be associated to individual atoms. Since both Eq. (\[HFV\]) and Eq. (\[sigma\]) are exact and non-perturbative, evaluating ${\bm{J}}_v(t)$ along the [*ab initio*]{} trajectory accounts for the full anharmonicity.

![Early (a) and late (b) decay of the heat flux autocorrelation function (HF-ACF) of Silicon computed in a 64 atom supercell with DFT-LDA at a temperature of 960 K (trajectory length $\sim207$ ps). The green line (${G}[{\bm{J}}_v]$) employs the virial [*ab initio*]{} heat flux ${\bm{J}}_v(t)$ that incorporates all anharmonic effects, whereas the blue and orange lines show the HF-ACFs ${G}[{\bm{J}}_v^{ha}]$ and ${G}[{\bm{J}}_v^{qp}]$ for approximate heat fluxes computed for the exact same trajectory, but imposing the harmonic approximation, i.e., using ${\bm{J}}_v^{ha}(t)$ and ${\bm{J}}_v^{qp}(t)$ defined via Eq. (\[VirialHarmonic\]) and (\[Peierls\]), respectively.[]{data-label="ACFSIraw"}](Fig1.pdf){width="0.85\linewidth"}

To validate our implementation of the proposed approach in the all-electron, numeric atomic orbital electronic structure code [*FHI-aims*]{} [@Blum:2009fe] we compare the heat flux autocorrelation function (HF-ACF) computed from first principles ${G}[{\bm{J}}_v]$ with the respective harmonic HF-ACF ${G}[{\bm{J}}_v^{ha}]$ by evaluating the approximate virial heat flux ${\bm{J}}_v^{ha}(t)$ using the harmonic force constants $\Phi_{IJ}^{\alpha\beta}=\partial^2 U / \partial{R}_I^\alpha \partial{R}_J^\beta$. In the [*harmonic approximation*]{}, the virials $$\left({\sigma_I}^{\alpha\beta}\right)^{ha-{\bm{R}}} = \frac{1}{2V} \sum_{J\neq I}  \Phi_{IJ}^{\alpha\beta}(\Delta {R}_I^\alpha - \Delta {R}_J^\alpha) ({R}_{I}^\beta - {R}_J^\beta) 
\label{VirialHarmonic}$$ depend only on the positions and displacements from equilibrium $\Delta {\bm{R}}_I={\bm{R}}_I-{\bm{R}}_I^{eq}$ [@Ladd:1986tv], so that ${\bm{J}}_v^{ha}(t)$ can be evaluated using Eq. (\[HFV\]) and (\[VirialHarmonic\]) along the exact same first-principles trajectory used to compute ${\bm{J}}_v(t)$. As an example, Fig. \[ACFSIraw\] shows such a comparison: ${G}[{\bm{J}}_v]$ and ${G}[{\bm{J}}_v^{ha}]$ closely resemble each other and become equal for large time-lags $\tau$, which demonstrates the validity of the introduced first-principles definition of the heat flux and its applicability in [*ab initio*]{} GK calculations.

However, Fig. \[ACFSIraw\] also neatly exemplifies the severe computational challenges of such first-principle GK simulations: Due to the limited time scales accessible in aiMD runs, thermodynamic fluctuations dominate the HF-ACF, which in turn prevents a reliable and numerically stable assessment of the thermal conductivity via Eq. (\[GreenKubo\]). Furthermore, achieving convergence with respect to system size is numerically even more challenging, as classical MD studies based on FFs [@He:2012hg; @Howell:2012fh] have shown, so that [*ab initio*]{} GK simulations of solids appear to be computationally prohibitively costly. However, as we will show below, the computational effort can be reduced by several orders of magnitude by a correct extrapolation technique employing a proper interpolation in reciprocal space.

![Thermal conductivity ${\kappa}_{\alpha\alpha}$ of Si at 300 and 1000 K computed with the Tersoff-FF. Values extrapolated (left: in time; right graph: in time and size) with our aiGK method are denoted by circles, whereas the dashed lines show the values resulting from a brute force evaluation of ${G}[{\bm{J}}_v]$ using only ${\bm{J}}_v$, i.e., without any extrapolation. 6912 ${\bm{q}}$-points corresponding to a $12\times12\times12$ cubic supercell were used for the size extrapolation in the right plot.[]{data-label="ConvergenceSi"}](Fig2.pdf){width="0.85\linewidth"}

For this purpose, we first note that in the harmonic approximation the HF-ACF can be equivalently [@Ladd:1986tv] evaluated in reciprocal space using the heat flux definition in the phonon picture [@Peierls:1929jv]: $${\bm{J}}_v^{ha-{\bm{q}}}(t) = \frac{1}{2V}\sum_{s{\bm{q}}} n_{s}({\bm{q}},t) \; \omega_{s}^2({\bm{q}})\; {\bm{v}}_{s}({\bm{q}}) \;.
\label{Peierls}$$ Here, the sum goes over all reciprocal space points ${\bm{q}}$ commensurate with the chosen supercell; $\omega_{s}({\bm{q}})$ are the eigenfrequencies and ${\bm{v}}_{s}({\bm{q}})$ are the group velocities of the phonon mode $s$, which are obtained by Fourier transforming and diagonalizing the mass-scaled force constant matrix $\Phi_{IJ}^{\alpha\beta}$ introduced in Eq. (\[VirialHarmonic\]). The time-dependent phonon amplitudes and thus the occupation numbers $n_{s}({\bm{q}},t)$ can be extracted from the MD trajectory using the techniques described in Ref. [@McGaughey:2004iz]. Accordingly, we can reformulate the HF-ACF as $\hat{{G}} = {G}[{\bm{J}}_v] - {G}[{\bm{J}}_v^{ha-{\bm{R}}}] + {G}[{\bm{J}}_v^{ha-{\bm{q}}}]$. For fully time and size converged calculations, $\hat{{G}}$ equals ${G}[{\bm{J}}_v]$; for underconverged calculations (cf. Fig. \[ACFSIraw\]), $\hat{{G}}$ exhibits significantly less thermodynamic fluctuations, since the phases of the individual modes do not enter Eq. (\[Peierls\]).

Even more importantly, this formalism enables a straightforward size-extrapolation by extending the sum over (the finite number of commensurate) reciprocal space points ${\bm{q}}$ in Eq. (\[Peierls\]) to a denser grid. The required frequencies $\omega_{s}({\bm{q}}')$ and group velocities ${\bm{v}}_{s}({\bm{q}}')$ can be determined on arbitrary ${\bm{q}}'$-points that are not commensurate with the supercell by Fourier interpolating the force constants $\Phi_{IJ}^{\alpha\beta}$ [@Parlinski:1997kr]. In the same spirit, we introduce the dimensionless quantity $$\Delta{n}_{s}({\bm{q}},\tilde{t}) = \frac{n_{s}({\bm{q}},t = \tilde{t}/\omega_{s}({\bm{q}}))- \left\langle n_{s}({\bm{q}}) \right\rangle}{\left\langle n_{s}({\bm{q}})\right\rangle}\;,
\label{RenormN}$$ which accounts for the fact that the equilibrium fluctuations of the occupation numbers are proportional to their equilibrium value $\left\langle n_{s}({\bm{q}}) \right\rangle = {2{k_{\mbox{\tiny B}}}T}/\omega_s^2({\bm{q}})$ and that the natural time scale $\tilde{t}$ of each mode is determined by its frequency $\omega_{s}({\bm{q}})$. As a consequence, we inherently account for the typical $1/\omega^2$ dependence of the phonon lifetimes [@Glassbrenner:1964fb; @Garg:2011hi] so that the respective ACFs ${G}[\Delta{n}_{s}({\bm{q}},\tilde{t})]$ become comparable and can be accurately interpolated in ${\bm{q}}$-space (see Suppl. Mat.).

We validate this approach by applying it to FF simulation of Si based on the Tersoff potential, which are known to be particularly challenging to convergence [@He:2012hg; @Howell:2012fh]. As shown in Fig. \[ConvergenceSi\] for 300 and 1000 K, the proposed interpolation yields remarkable improvements with respect to size- and time-convergence: Panel (a) shows the dependence of ${\bm{\kappa}}$ on the trajectory length, while panel (b) shows the dependence of ${\bm{\kappa}}$ on the supercell size. Compared to traditional brute force GK simulations, we achieve reliable values for ${\bm{\kappa}}$ with sizes as small as 64 atoms and with trajectory lengths as short as 200 ps. This translates into a computational speed-up of more than three orders of magnitude, which in turn enables [*ab initio*]{} Green-Kubo calculations (aiGK) with reasonable numerical effort.

![Thermal conductivity ${\kappa}_{\alpha\alpha}$ computed using our aiGK method for Si (LDA green and PBEsol orange triangles) and tetragonal ZrO$_2$ (LDA: red squares; PBEsol: blue diamonds). Lattice expansion using the quasi-harmonic approximation and the tetragonal-cubic phase transition in ZrO$_2$ [@Carbogno:2014wa] are taken into account. Tabulated values of $\kappa$ are listed in the Suppl. Information. Black circles denote experimental results for single crystals (or extrapolated to this limit) compiled from Ref. [@Glassbrenner:1964fb; @Bisson:2000tu; @Youngblood:1988ws; @Raghavan:1998va]. The lines are generated by fitting the data with ${\kappa}(T)=a+b/T^c$. []{data-label="SiZrO2"}](Fig3.pdf){width="0.85\linewidth"}

![Accumulated thermal conductivity versus mean free path for Si (Tersoff/DFT-LDA) and ZrO$_2$ (DFT-PBEsol) at different temperatures.[]{data-label="VACF"}](Fig4.pdf){width="0.85\linewidth"}

Next, we apply our aiGK technique to compute the temperature-dependent thermal conductivity of Si (64-atom cell) and tetragonal ZrO$_2$ (96-atom cell) both for the LDA and PBEsol functional. Six trajectories with a duration of 200 ps were used for each data point. In both cases (cf. Fig. \[SiZrO2\]), we obtain good agreement with experimental data [@Glassbrenner:1964fb; @Bisson:2000tu; @Youngblood:1988ws; @Raghavan:1998va]. For Si, we note that our aiGK extrapolation technique is responsible for up to 50% of ${\bm{\kappa}}$, especially at low temperatures. Conversely, the thermal conductivities for ZrO$_2$ are mostly size- and time-converged within the aiMD regime, so that corrections stemming from the extrapolation are always smaller than $10$%. As the accumulated thermal conductivities in Fig. \[VACF\] show, the reason for this behavior are the exceptionally low phonon lifetimes or short mean free paths in ZrO$_2$. Interestingly, we can trace back this behavior to the peculiar, very anharmonic dynamics of zirconia at elevated temperatures. Under such thermodynamic conditions, the oxygen atoms and the lattice tetragonality spontaneously re-orient along a different cartesian direction in a switching mechanism [@Carbogno:2014wa]. These switches between different local minima of the potential-energy surface constitute a severe violation of the harmonic approximation, given that the dynamics does no longer evolve around one minimum. This is also reflected in the respective dependence on the functional: At the LDA level, the respective barriers are underestimated [@Carbogno:2014wa], so that ${\bm{\kappa}}$ decreases more drastically and at lower temperatures.

In summary, we presented an [*ab initio*]{} implementation of the Green-Kubo method that is applicable for the computation of thermal conductivities in solids, since the proposed unique definition of the heat flux is derived from the virial theorem and based on the local stress tensor. The developed extrapolation technique, which significantly lowers the required computational effort, enables to perform such computations within moderate time- and length scales for the aiMD trajectories. By this means, we are able to accurately compute thermal conductivities for both extremely harmonic and anharmonic materials on the same footing and at arbitrarily high temperatures. In particular, we are able to investigate materials with very low thermal conductivities and high anharmonicities (thermal barriers), for which perturbative treatments relying on the approximate validity of the harmonic approximation would fail. Accordingly, the proposed technique enables for the first time to perform accurate first-principles studies of such materials, which play a pivotal role in a multitude of scientific and technological applications, e.g., as thermal barrier coatings [@Levi:2004vh] and thermoelectric elements [@Snyder:2008dc].

The authors thank J. Meyer, S. K. Estreicher, and C. G. Van de Walle for fruitful discussions. We are grateful to S. Lampenscherf and C. Levi for pointing out to us that a general, non-perturbative theory of heat transport was missing, so far. R. R. acknowledges the Alexander von Humboldt Foundation. The project received funding from the Einstein foundation (project ETERNAL), the DFG cluster of excellence UNICAT, and the European Union’s Horizon 2020 research and innovation program under grant agreement no. 676580 with The Novel Materials Discovery (NOMAD) Laboratory, a European Center of Excellence.
