Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach

TL;DR

This paper presents a unified lattice-dynamical approach to model heat transport in both crystals and glasses, integrating quantum effects and bridging existing models, validated against silicon data.

Contribution

A novel theoretical framework that unifies heat transport modeling in crystals and glasses using lattice dynamics and Green-Kubo theory, including quantum effects.

Findings

Successfully bridges Boltzmann and Allen-Feldman models

Accurately predicts heat transport in silicon structures

Incorporates quantum mechanical effects naturally

Abstract

We introduce a novel approach to model heat transport in solids, based on the Green-Kubo theory of linear response. It naturally bridges the Boltzmann kinetic approach in crystals and the Allen-Feldman model in glasses, leveraging interatomic force constants and normal-mode linewidths computed at mechanical equilibrium. At variance with molecular dynamics, our approach naturally and easily accounts for quantum mechanical effects in energy transport. Our methodology is carefully validated against results for crystalline and amorphous silicon from equilibrium molecular dynamics and, in the former case, from the Boltzmann transport equation.

Figures

Equations

$${\kappa }_{\alpha \beta }=\frac{1}{V{k}_B{T}^2}\underset{0}{\overset{+\infty }{\int }}⟨J_{\alpha }(t)J_{\beta }(0)⟩dt,$$

$$J_{\alpha }=\frac{1}{2}\sum _{ij\beta \gamma }(R_{i\alpha }^{\circ }-R_{j\alpha }^{\circ }){\Phi }_{i\beta }^{j\gamma }{u}_{i\beta }{\dot{u}}_{j\gamma },$$

$${\kappa }_{\alpha \beta }=\frac{k_B}{V}\sum _{nm}{v}_{nm}^{\alpha }{v}_{nm}^{\beta }{\tau }_{nm}^{\circ },$$

$$v_{nm}^{\alpha }=\frac{1}{2\sqrt{{\omega }_n{\omega }_m}}\sum _{ij\beta \gamma }\frac{R_{i\alpha }^{\circ }-R_{j\alpha }^{\circ }}{\sqrt{M_i{M}_j}}{\Phi }_{i\beta }^{j\gamma }{e}_n^{i\beta }{e}_m^{j\gamma },$$

$${\tau }_{nm}^{\circ }=\frac{{\gamma }_n+{\gamma }_m}{({\gamma }_n+{\gamma }_m{)}^2+({\omega }_n-{\omega }_m{)}^2}+\mathcal{O}({ϵ}^2),$$

$$J_{\beta }=\frac{i}{2}\sum _{nm}{v}_{nm}^{\beta }{\omega }_m({\alpha }_n^{\ast }+{\alpha }_n)({\alpha }_m^{\ast }-{\alpha }_m).$$

$${\kappa }_{\alpha \beta }=\frac{1}{VT}\underset{0}{\overset{1/k_{B}T}{\int }}d\lambda \underset{0}{\overset{+\infty }{\int }}dt⟨{\hat{J}}_{\alpha }(t+i\hslash \lambda ){\hat{J}}_{\beta }(0)⟩,$$

$${\hat{J}}_{\beta }=\frac{i\hslash }{2}\sum _{nm}{v}_{nm}^{\beta }{\omega }_m({\hat{a}}_n^{†}+{\hat{a}}_n)({\hat{a}}_m^{†}-{\hat{a}}_m).$$

$${\kappa }_{\alpha \beta }=\frac{1}{V}\sum _{nm}{c}_{nm}{v}_{nm}^{\alpha }{v}_{nm}^{\beta }{\tau }_{nm}^{\circ },$$

$${\kappa }^{\prime }(\omega )=\frac{1}{3V}\sum _{\alpha }\sum _{nm}\Delta (\omega -{\omega }_n)c_{nm}(v_{nm}^{\alpha }{)}^2{\tau }_{nm}^{\circ },$$

$$\kappa (\omega )={\int }_0^{\omega }{\kappa }^{\prime }({\omega }^{\prime })d{\omega }^{\prime }.$$

$$\begin{array}{r}{\xi }_n=\sum _{i\alpha }\sqrt{M_i}{u}_{i\alpha }{e}_{ni}^{\alpha };{\pi }_n=\sum _{i\alpha }\frac{1}{\sqrt{M_i}}{\dot{u}}_{i\alpha }{e}_{ni}^{\alpha };{\alpha }_n=\sqrt{\frac{{\omega }_n}{2}}{\xi }_n+\frac{i}{\sqrt{2{\omega }_n}}{\pi }_n;\\ H=\frac{1}{2}\sum _n\left({\pi }_n^2+{\omega }_n^2{\xi }_n^2\right)=\sum _n{\omega }_n|{\alpha }_n{|}^2;J=\frac{i}{2}\sum _{nm}{v}_{nm}{\omega }_m({\alpha }_n^{\ast }+{\alpha }_n)({\alpha }_m^{\ast }-{\alpha }_m),\end{array}$$

$${\alpha }_n(t)={\alpha }_n{\mathrm{e}}^{i{\omega }_{n}t},$$

$$J(t)J(0)=P_4(\alpha ,{\alpha }^{\ast };t)=-\frac{1}{4}\sum _{mnpq}{v}_{nm}{\omega }_m{v}_{pq}{\omega }_q({\alpha }_n^{\ast }(t)+{\alpha }_n(t))({\alpha }_m^{\ast }(t)-{\alpha }_m(t))({\alpha }_p^{\ast }+{\alpha }_p)({\alpha }_q^{\ast }-{\alpha }_q).$$

$$⟨ABCD⟩=⟨AB⟩⟨CD⟩+⟨AC⟩⟨BD⟩+⟨AD⟩⟨BC⟩,$$

$$⟨{\alpha }_n{\alpha }_m⟩=0;⟨{\alpha }_n^{\ast }(t){\alpha }_m⟩={\delta }_{mn}{g}_n(t),$$

$$⟨{\alpha }_n(t){\alpha }_m(t){\alpha }_p^{\ast }{\alpha }_q^{\ast }⟩=⟨{\alpha }_n(t){\alpha }_p^{\ast }⟩⟨{\alpha }_m(t){\alpha }_q^{\ast }⟩+⟨{\alpha }_n(t){\alpha }_q^{\ast }⟩⟨{\alpha }_m(t){\alpha }_p^{\ast }⟩=({\delta }_{np}{\delta }_{mq}+{\delta }_{nq}{\delta }_{mp})g_n(t)g_m(t).$$

$$⟨J(t)J(0)⟩=-\frac{1}{4}\sum _{nm}\left(\frac{{\omega }_n-{\omega }_m}{{\omega }_n}(g_n(t)g_m(t)+g_n^{\ast }(t)g_m^{\ast }(t))-\frac{{\omega }_n+{\omega }_m}{{\omega }_n}(g_n(t)g_m^{\ast }(t)+g_n^{\ast }(t)g_m(t))\right).$$

$$g_n(t)=\frac{k_{B}T}{{\omega }_n}{e}^{i({\omega }_n+i{\gamma }_n)t}.$$

$${\kappa }_{\alpha \beta }=\frac{k_B}{V}\sum _{nm}{v}_{nm}^{\alpha }{v}_{nm}^{\beta }{\tau }_{nm},v_{nm}^{\alpha }=\frac{1}{2\sqrt{{\omega }_n{\omega }_m}}\sum _{ij\beta \gamma }\frac{R_{i\alpha }^{\circ }-R_{j\alpha }^{\circ }}{\sqrt{M_i{M}_j}}{\Phi }_{i\beta }^{j\gamma }{e}_n^{i\beta }{e}_m^{j\gamma },$$

$$\begin{gathered} {\tau }_{nm}=\frac{({\omega }_n+{\omega }_m{)}^2}{4{\omega }_n{\omega }_m}\frac{{\gamma }_n+{\gamma }_m}{({\gamma }_n+{\gamma }_m{)}^2+({\omega }_n-{\omega }_m{)}^2} \\ +\frac{({\omega }_n-{\omega }_m{)}^2}{4{\omega }_n{\omega }_m}\frac{{\gamma }_n+{\gamma }_m}{({\gamma }_n+{\gamma }_m{)}^2+({\omega }_n+{\omega }_m{)}^2}. \end{gathered}$$

$$\begin{gathered} G_n(t)=\hslash ⟨a_n^{†}(t)a_n(0)⟩=\hslash n_k{e}^{i{\omega }_{n}t}, \\ {\tilde{G}}_n(t)=\hslash ⟨a_n(t)a_n^{†}(0)⟩=\hslash (n_k+1)e^{-i{\omega }_{n}t}. \end{gathered}$$

$$\langle\hat{J}(\tau)\hat{J}(0)\rangle=-\frac{1}{4}\sum\limits_{nm}\left(\frac{\omega_{n}-\omega_{m}}{\omega_{n}}\Bigl{(}G_{n}(\tau)G_{m}(\tau)+\tilde{G}_{n}(\tau)\tilde{G}_{m}(\tau)\Bigr{)}-\frac{\omega_{n}+\omega_{m}}{\omega_{n}}\Bigl{(}G_{n}(\tau)\tilde{G}_{m}(\tau)+\tilde{G}_{n}(\tau)G_{m}(\tau)\Bigr{)}\right),$$

$$\begin{gathered} {\kappa }_{\alpha \beta }=\frac{{\hslash }^2}{4VT}\sum _{nm}{v}_{nm}^{\alpha }{v}_{nm}^{\beta }(\frac{n_n-n_m}{\hslash ({\omega }_n-{\omega }_m)}\frac{{\gamma }_n+{\gamma }_m}{({\omega }_n-{\omega }_m{)}^2+({\gamma }_n+{\gamma }_m{)}^2}({\omega }_n+{\omega }_n{)}^2 \\ +\frac{e^{\beta \hslash ({\omega }_n+{\omega }_m)}-1}{\hslash ({\omega }_n+{\omega }_m)}\frac{2({\omega }_m-{\omega }_m{)}^2}{(e^{\beta \hslash {\omega }_n}-1)(e^{\beta \hslash {\omega }_m}-1)}\frac{{\gamma }_n+{\gamma }_m}{({\omega }_n+{\omega }_m{)}^2+({\gamma }_n+{\gamma }_m{)}^2}) \end{gathered}$$

$${\kappa }_{\alpha \beta }=\frac{1}{V}\sum _{nm}{c}_{nm}{v}_{nm}^{\alpha }{v}_{nm}^{\beta }{\tau }_{nm}^{\circ }.$$

Full text

Heat transport in solid insulators, either crystalline or disordered, is dominated by the dynamics of lattice vibrations. Far from melting, atomic displacements from equilibrium are much smaller than interatomic distances and they can thus be treated in the (quasi-) harmonic approximation. In crystals this observation enables a kinetic description of heat transport in terms of phonons that can be embodied in the Peierls–Boltzmann transport equation (BTE) *[SeeforinstanceChapterVIIof][]Ziman1960; Fugallo2018. In disordered systems the typical phonon mean free paths may be so short that the quasi-particle picture of heat carriers breaks down and BTE is no longer applicable, making it necessary to resort to molecular dynamics (MD), in either its nonequilibrium or equilibrium (EMD) flavors HMM; Fugallo2018. MD is of general applicability to solids, either periodic or disordered, and liquids. Yet, it may require long simulation times and it is subject to statistical errors, which are at times cumbersome to evaluate especially for systems at low temperatures, where lack of ergodicity may be an issue. Most importantly, MD cannot account for quantum-mechanical effects BedoyaMartinez:2014cq, which are instead easily treated in BTE, thus making the treatment of heat transport for glasses in the quantum regime, i.e. below the Debye temperature, a methodological challenge.

In this paper we present a novel approach to heat transport in insulating solids, which combines the Green-Kubo (GK) theory of linear response Green1952; *Green1954; *Kubo1957a; *Kubo1957b; HMM and a quasi-harmonic description of lattice vibrations, thus resulting in a compact expression for the thermal conductivity that unifies the BTE approach in the single-mode relaxation-time approximation (RTA) for crystals Fugallo2018 and a generalization of the Allen-Feldman (AF) model for disordered system Allen1989; *Allen1993 that explicitly and naturally accounts for normal-mode lifetimes. The main ingredients of our approach are the matrix of the inter-atomic force constants (IFC) computed at mechanical equilibrium and the anharmonic lifetimes of the vibrational modes, as computed from the cubic corrections to the harmonic IFCs using the Fermi’s golden rule Fabian1996. Our theory is thoroughly validated in crystalline and amorphous silicon by comparing its predictions with those of EMD simulations and BTE computations.

The basis of our work is the GK theory of linear response Kubo1957a; *Kubo1957b; HMM, which relates the heat conductivity to the ensemble average of the heat-flux auto-correlation function:

[TABLE]

where $V$ and $T$ are the system volume and temperature, $k_{B}$ is the Boltzmann constant, $J_{\alpha}(t)$ the $\alpha$-th Cartesian component of the macroscopic heat flux, which in solids coincides with the energy flux, and $\langle\cdot\rangle$ indicates a canonical average over initial conditions HMM. Far from melting, the energy flux and the lattice Hamiltonian of a solid, both crystalline and amorphous, can be expressed as power series in the atomic displacements, and Eq. (1) can be evaluated in terms of Gaussian integrals using standard field-theoretical techniques.

The energy flux can be expressed in terms of atomic positions, $\mathbf{R}_{i}$, and local energies, $\epsilon_{i}$, as $\mathbf{J}=\sum_{i}(\dot{\mathbf{R}}_{i}\epsilon_{i}+\mathbf{R}_{i}\dot{\epsilon}_{i})$ HMM, where in the harmonic approximation $\epsilon_{i}$ can be defined as: $\epsilon_{i}=\frac{M_{i}}{2}\sum_{\alpha}\left(\dot{u}_{i\alpha}\right)^{2}+\frac{1}{2}\sum_{j\alpha\beta}u_{i\alpha}\Phi_{i\alpha}^{j\beta}u_{j\beta}$, $M_{i}$ being the mass of $i$-th atom, $\mathbf{u}_{i}=\mathbf{R}_{i}-\mathbf{R}^{\circ}_{i}$ its displacement from its equilibrium position, $\mathbf{R}^{\circ}_{i}$, $\alpha$ and $\beta$ label Cartesian components, and $\Phi_{i\alpha}^{j\beta}=\left.\frac{\partial^{2}E}{\partial u_{i\alpha}\partial u_{j\beta}}\right|_{\mathbf{u}=0}$ is the IFC matrix. By expressing the energy flux in terms of the $\mathbf{u}$’s, one obtains: $\mathbf{J}=\mathbf{J}^{\circ}+\frac{d}{dt}\sum_{i}\mathbf{u}_{i}\epsilon_{i}$, where $\mathbf{J}^{\circ}=\sum_{i}\mathbf{R}^{\circ}_{i}\dot{\epsilon}_{i}$. The second term on the right-hand side of this expression is the total time derivative of a process that, in the absence of atomic diffusion, is stationary and of finite variance. A recently established gauge invariance principle for heat transport Marcolongo2016; *Ercole2016 states that such a total time derivative does not contribute to the thermal conductivity. We will therefore dispose of it and express the energy flux as: $\mathbf{J}\leftarrow\mathbf{J}^{\circ}$. Note that it is essential to use equilibrium atomic positions in the definition of $\mathbf{J}^{\circ}$, i.e. the positions describing the (metastable) mechanical equilibrium of any given model of an ordered or disordered system, rather than instantaneous ones. Otherwise, the difference $\mathbf{J}-\mathbf{J}^{\circ}$ would not be a total time derivative and the value of the transport coefficient resulting from $\mathbf{J}^{\circ}$ would be biased. By making use of Newton’s equation of motion, the final expression for the harmonic heat flux reads Allen1989:

[TABLE]

where the minimum-image convention is adopted for computing inter-atomic distances in periodic boundary conditions.

By inserting Eq. (2) into Eq. (1), the integrand is cast into the canonical average of a quartic polynomial in the atomic positions and velocities. In the harmonic approximation, this average reduces to a Gaussian integral, which can be evaluated by way of the Wick’s theorem *[SeeforinstanceSec.2.3in][]Negele1988. By doing so, the resulting time integral would diverge, thus yielding an infinite conductivity, as expected for a harmonic crystal Rieder:1967fi. In order to regularize this integral, we introduce anharmonic effects by renormalizing the single-mode correlation functions using the normal-mode lifetimes, as explained below. Our final result for the heat conductivity tensor reads:

[TABLE]

where $\omega_{n}$ and $\gamma_{n}$ are the harmonic frequency and decay rate of the $n$-th normal mode, and $e_{ni}^{\alpha}$ is the corresponding eigenvector of the dynamical matrix, $\sum_{j\beta}\frac{1}{\sqrt{M_{i}M_{j}}}\Phi^{i\alpha}_{j\beta}e_{nj}^{\beta}=\omega_{n}^{2}e_{ni}^{\alpha}$, and $\epsilon$ indicates the ratio $\gamma/\omega$, which vanishes in the harmonic limit. Eqs. (3-5) will be dubbed as the quasi-harmonic Green-Kubo (QHGK) approximation for the heat conductivity.

In order to establish Eq. (3), we first express the energy flux in Eq. (2) in terms of normal-mode coordinates and momenta, defined as: $\xi_{n}=\sum_{i\alpha}\sqrt{M_{i}}u^{i}_{\alpha}e_{ni}^{\alpha}$ and $\pi_{n}=\sum_{i\alpha}\dot{u}_{i}^{\alpha}e_{ni}^{\alpha}/\sqrt{M_{i}}$, reading: $J^{\alpha}=\sum_{nm}v_{nm}^{\alpha}\sqrt{\omega_{n}\omega_{m}}\xi_{n}\pi_{m}$. It is then convenient to express these normal-mode coordinates and momenta in terms of classical complex amplitudes, reminiscent of the quantum boson ladder operators and defined as: $\alpha_{n}=\sqrt{\frac{\omega_{n}}{2}}\xi_{n}+\frac{i}{\sqrt{2\omega_{n}}}\pi_{n}$, whose time evolution is $\alpha_{n}(t)=\alpha_{n}(0)\mathrm{e}^{-i\omega_{n}t}$. By doing so, the energy flux can be expressed in terms of the $\alpha$ amplitudes as

[TABLE]

By using this expression, the integrand in Eq. (1) becomes a Gaussian integral of a fourth-order polynomial in the $\alpha$’s and $\alpha^{*}$’s that, by means of the Wick’s theorem *[SeeforinstanceSec.2.3in][]Negele1988, can be cast into a sum of products of pairs of single-mode (classical) Green’s functions, $\langle\alpha^{*}_{n}(t)\alpha_{m}(0)\rangle=\delta_{nm}g_{n}(t)$ and $\langle\alpha_{n}(t)\alpha_{m}(0)\rangle=0$. In the purely harmonic approximation, one would have $g^{\circ}_{n}(t)=\frac{k_{B}T}{\omega_{n}}\mathrm{e}^{i\omega_{n}t}$. Anharmonic effects broaden the vibrational lines by a finite line-width, $\gamma_{n}$, which results in a finite imaginary part of the frequency and in a decay of the single-mode Green’s function, reading: $g_{n}(t)=\frac{k_{B}T}{\omega_{n}}\mathrm{e}^{i(\omega_{n}+i\gamma_{n})t}$. By plugging this expressions into the lengthy formula that results from applying Wick’s theorem to the integrand of Eq. (1) and performing the time integral, after some cumbersome but straightforward algebra we get Eq. (3). A full derivation of Eqs. (3-5) is presented in the Supplemental Material (SM), Sec. S1.

To lowest order in the anharmonic interactions, vibrational linewidths can be computed from the classical limit of the Fermi golden rule, \gamma_{n}=\frac{\pi\hbar^{2}}{8\omega_{n}}\sum_{ml}\frac{|V^{\prime\prime\prime}_{nml}|^{2}}{\omega_{m}\omega_{l}}\bigl{[}\frac{1}{2}(1+n_{m}+n_{l})\delta(\omega_{n}-\omega_{m}-\omega_{l})+(n_{m}-n_{l})\delta(\omega_{n}+\omega_{m}-\omega_{l})\bigr{]}, where $n_{l}$ is the Bose-Einstein occupation number of the $l$-th normal mode and $V^{\prime\prime\prime}_{nlm}=\frac{\partial^{3}V}{\partial\xi^{n}\partial\xi^{l}\partial\xi^{m}}$ is the third derivative of the potential energy with respect to the amplitude of the lattice distortion along the lattice normal modes Fabian1996.

In order to show that our QHGK expression for the thermal conductivity, Eq. (3), reduces to the predictions of the BTE-RTA in crystals, we first notice that the $v^{\alpha}$ matrices of Eq. (4) can be expressed in terms of the matrix elements of the matrices \bigl{(}V^{\alpha}\bigr{)}_{i\gamma}^{j\delta}=\frac{R^{\circ}_{i\alpha}-R^{\circ}_{j\alpha}}{2\sqrt{M_{i}M_{j}}}\Phi_{i\gamma}^{j\delta} between normal-mode eigenvectors: v^{\alpha}_{nm}=\bigl{(}e_{n},V^{\alpha}\cdot e_{m}\bigr{)}/\sqrt{\omega_{n}\omega_{m}}, where the notations “$(e,e^{\prime})$” and “$V\cdot e$” indicate scalar and matrix-vector products in the space of atomic displacements. In crystals equilibrium atomic positions are characterised by a discrete lattice position, $\mathbf{a}_{i}$, and by an integer label, $s_{i}$, indicating different atomic sites within a unit cell, $\mathbf{d}_{s}$: $\mathbf{R}^{\circ}_{i}=\mathbf{a}_{i}+\mathbf{d}_{s_{i}}$. Likewise, in the Bloch representation, normal modes can be labelled by a quasi-discrete wavevector, $\mathbf{q}$, belonging to the first Brillouin zone (BZ), and by a band index, $\nu$: $n\to(\mathbf{q}_{n},\nu_{n})$. In particular, the IFC matrix and its eigenvectors can be expressed as $\frac{1}{\sqrt{M_{i}M_{j}}}\Phi_{j\beta}^{i\alpha}=\sum_{\mathbf{q}}\mathrm{e}^{i\mathbf{q}\cdot(\mathbf{R}^{\circ}_{i}-\mathbf{R}^{\circ}_{j})}D_{t\beta}^{s\alpha}(\mathbf{q})$, where $D_{t\beta}^{s\alpha}(\mathbf{q})$ is the dynamical matrix of the system and $\eta^{s\alpha}_{\nu\mathbf{q}}$ its eigenvectors: $e^{\alpha}_{nu}=\mathrm{e}^{i\mathbf{q}_{n}\cdot\mathbf{R}^{\circ}_{i}}\eta^{s_{i}\alpha}_{\nu_{n}\mathbf{q}_{n}}$ and $\sum_{t\beta}D^{s\alpha}_{t\beta}(\mathbf{q})\eta_{\nu\mathbf{q}}^{t\beta}=\omega_{\nu\mathbf{q}}^{2}\eta^{s\alpha}_{\nu\mathbf{q}}$. When normal-mode eigenvectors are chosen to be real, the $v^{\alpha}$ matrices of Eq. (4) are real and anti-symmetric. In particular, $v^{\alpha}_{nn}=0$ and a non-vanishing thermal conductivity results from the matrix elements of $v^{\alpha}$ connecting (quasi-) degenerate normal modes, i.e. modes with frequencies that coincide within the sum of their line widths. In the Bloch representation, $v^{\alpha}$ is anti-Hermitian and block-diagonal with respect to the wave-vector, $\mathbf{q}$. Its diagonal elements are imaginary, though not necessarily vanishing. In this representation one has: v^{\alpha}_{\nu\mathbf{q},\mu\mathbf{p}}=i\frac{\delta_{\mathbf{q}\mathbf{p}}}{\sqrt{\omega_{\nu\mathbf{q}}\omega_{\mu\mathbf{p}}}}\bigl{(}\eta_{\nu\mathbf{q}},D^{\alpha}(\mathbf{q})\cdot\eta_{\mu\mathbf{q}}\bigr{)}, where $D^{\alpha}(\mathbf{q})=\frac{\partial D(\mathbf{q})}{\partial q^{\alpha}}$. At fixed $\mathbf{q}$, the vibrational spectrum is strictly discrete i.e. it remains so even in the thermodynamic limit. The number of $\mathbf{q}$ points for which there exists a pair of distinct modes, $(\mathbf{q},\nu)$ and $(\mathbf{q},\mu)$ with $\nu\neq\mu$, that are degenerate within the sum of their line-widths $(|\omega_{\mathbf{q}\nu}-\omega_{\mathbf{q}\mu}|\lesssim\gamma_{\mathbf{q}\nu}+\gamma_{\mathbf{q}\mu})$ is vanishingly small, because, in practice, this quasi-degeneracy can only occur close to high-symmetry lines. Furthemore, for such few pairs, one can prove that $v_{\nu\mathbf{q},\mu\mathbf{q}}\propto\omega_{\nu\mathbf{q}}-\omega_{\mu\mathbf{q}}$. Hence in the periodic case the $\tau^{\circ}$ matrix in Eq. 5 is strictly diagonal, $\tau^{\circ}_{\mathbf{q}\nu,\mathbf{p}\mu}=\delta_{\mathbf{q}\mathbf{p}}\delta_{\nu\mu}\tau^{\circ}_{\mathbf{q}\nu}$, where $\tau^{\circ}_{\mathbf{q}\nu}=\frac{1}{2\gamma_{\mathbf{q}\nu}}$ is the anharmonic lifetime of the $(\mathbf{q},\nu)$ normal mode. We conclude that, for periodic systems in the Bloch representation, the double sum in Eq. (3) can be cast into a single sum over diagonal terms, reading: $\kappa_{\alpha\beta}=\sum_{\mathbf{q}\nu}v_{\nu\mathbf{q}}^{\alpha}v_{\nu\mathbf{q}}^{\beta}\tau_{\nu\mathbf{q}}$, where v_{\nu\mathbf{q}}^{\alpha}=\frac{1}{2\omega_{\nu\mathbf{q}}}\bigl{(}\eta_{\nu\mathbf{q}},D^{\alpha}(\mathbf{q})\cdot\eta_{\nu\mathbf{q}}\bigr{)}=\frac{\partial\omega_{\nu\mathbf{q}}}{\partial q^{\alpha}} is the group velocity of the $\nu$-th phonon branch. This is the final expression for the thermal conductivity of a crystal in QHGK, which remarkably coincides with the solution of BTE-RTA.Ziman1960 We tested the QHGK approach against BTE-RTA for a crystalline silicon supercell of 1728 atoms, with a lattice parameter of 5.431 Å. The two calculations, performed with different codes, give the same results, as expected by the proven equivalence of the two methods for crystalline systems (see Figure S1 in SM).

Moving to the quantum regime is straightforward in our approach. To this end, we start from the quantum GK formula Green1952; *Green1954; Kubo1957a; *Kubo1957b; HMM, reading:

[TABLE]

where $\hat{J}_{\alpha}$ are quantum heat-flux operators and $\langle\cdot\rangle$ indicates quantum canonical averages. A quantum expression for the heat flux is obtained from its classical counterpart, Eq. (6), by making the substitutions: $\alpha\to\sqrt{\hbar}\hat{a}$ and $\alpha^{*}\to\sqrt{\hbar}\hat{a}^{\dagger}$, $\hat{a}^{\dagger}$ and $\hat{a}$ being normal-mode creation/annihilation operators, satisfying the standard commutation relations for bosons: $\left[{\hat{a}},{\hat{a}}^{\dagger}\right]=1$. Note that no ordering ambiguities arise when quantizing Eq. (6) because the $v^{\alpha}_{nm}$ matrices are antisymmetric, and they therefore vanish for $n=m$. The resulting expression for the quantum heat flux is:

[TABLE]

The computation of the heat conductivity proceeds exactly as in the classical case, except for the expressions for the single-mode Green’s functions. In the quantum case they read: $\langle\hat{a}_{n}^{\dagger}(t)\hat{a}_{n}(0)\rangle=n_{n}\mathrm{e}^{i(\omega_{n}+i\gamma_{n})t}$ and $\langle\hat{a}_{n}(t)\hat{a}^{\dagger}_{n}(0)\rangle=(n_{n}+1)\mathrm{e}^{-i(\omega_{n}-i\gamma_{n})t}$, $n_{n}=1\left/\left(e^{\frac{\hbar\omega_{n}}{k_{B}T}}-1\right)\right.$ being the Bose-Einstein distribution function. The final quantum-mechanical expression for the heat conductivity in the QHGK is:

[TABLE]

with $c_{nm}=\frac{\hbar\omega_{n}\omega_{m}}{T}\frac{n_{m}-n_{n}}{\omega_{m}-\omega_{n}}$. For $n=m$ this term reduces to the modal heat capacity $c_{n}=k_{B}\left(\frac{\hbar\omega_{n}}{k_{B}T}\right)^{2}\frac{e^{\hbar\omega_{n}/k_{B}T}}{(\mathrm{e}^{\hbar\omega_{n}/k_{B}T}-1)^{2}}$. The other symbols are the same as in Eqs. (4-5) for the classical case. $\tau^{\circ}_{nm}$, in particular, is only different from zero for $|\omega_{n}-\omega_{m}|\lesssim\gamma_{n}+\gamma_{m}$. Following the same derivation as for the classical case, one can prove that for periodic crystals Eq. (9) reduces to BTE-RTA. Furthermore, in the classical limit, one has $\lim_{T\to\infty}c_{nn}=k_{B}$ and the quantum formula, Eq. (9), reduces to Eq. (3). Further details are given in Sec. S2 of SM.

We validate our QHGK approach by testing the results of Eqs. (3) and (9) against MD simulations for amorphous silicon. Interatomic interactions are modeled using the empirical bond-order Tersoff potential Tersoff:1989tr, which describes well the thermal conductivity of bulk and nanostructured silicon, including a-Si Allen1993; He:2011wq; He:2012tq; Larkin2014. We consider a 1728-atom model of a-Si, generated by MD by quenching from the melt. Several EMD simulations where then run at different temperatures, as described in SM Fan:2015ba; Fan2017. The integral of the heat flux autocorrelation function in Eq. (1) can then be efficiently evaluated via cepstral analysis, as described in Refs. HMM and Ercole2017, which can be enhanced by averaging over multiple trajectories at low temperature ($T\leq 300$ K). Details on the data analysis procedure followed here and on the estimate of the statistical errors is given in the S3 section of the SM. The results of these calculations are reported in Figure 1 and exhibit a weak non-monotonic dependence on $T$. Performing similar MD simulations on models of increasing size (4,096 and 13,824 atoms) generated with the same protocol, we have verified that size effects on $\kappa$ at 300 K amount to less than 15$\%$, which is of the same order as the variation $\kappa$ among different models with the same size. The computation of the IFC matrix, normal-mode frequencies and lifetimes is described in detail in SM, where we also display the resulting dependence of lifetimes on temperature (Figure S2 of SM). The resulting strongly diagonally dominant form of the $\tau^{\circ}$ matrices in Eq. (5) is also displayed in Figure S2 of SM.

The thermal conductivity obtained by QHGK is in excellent agreement with that computed by EMD for $T\leq 600$ K (Figure 1). At higher temperatures QHGK overestimates $\kappa$, as it neglects higher-order anharmonic effects. We point out that, in spite of the formal analogies with the AF model Allen1989; Allen1993; Feldman:1993tn and recent refinements thereof Donadio:2010kp; Zhu:2016ks, Eq. (3) entails no empirical parameters. It thus allows one to predict temperature trends dictated by anharmonic effects in good agreement with MD  footnote2 , without making any a priori distinction among propagating, diffusive, or localized vibrational modes. Similarly to the GK modal analysis approach Lv:2016dr, based on classical MD, the transport character of the modes is dictated by the relative contribution from the diagonal and slightly off-diagonal terms of the $v^{\alpha}_{nm}$ matrix, weighted by $\tau^{\circ}_{nm}$ (Figure S2). The generality of QHGK is expected to have a major impact for the study of weakly disordered systems, which are beyond the scope of applicability of approaches based on the BTE and the AF model.

QHGK is a general theory that allows one to accurately calculate thermal transport in both crystals and glasses at a full quantum mechanical level. In Figure 2 we report our results from quantum QHGK calculation for an amorphous Si model of 13824 atoms along with three sets of experimental data Cahill1994; Zink:2006wt. QHGK results are in excellent agreement with the measurements in Zink:2006wt above 100 K. At lower temperature the estimate of $\kappa$ is affected by finite size effects, related to insufficient sampling of low-frequency acoustic modes: at 25 K these effects are so important, as to make the estimated conductivity almost vanish (see below) footnote. In fact, in order to eliminate finite-size effects, in our approach it would be necessary that in any relevant frequency range the density of vibrational states is larger than the normal-mode lifetimes, so that as many quasi-discrete normal modes as possible fall withing a line-width. In the low-frequency region, which is the most populated one in the quantum regime, this condition is hindered by the vanishing of both the density of states per unit volume and normal-mode line-widths. This effect is showcased in Figure 3, where we compare for different temperatures and model sizes the frequency-resolved differential conductivity,

[TABLE]

where $\Delta(\omega)$ is a broadened approximation of the $\delta$ function and the other symbols have the same meaning as in Eq. (9), and the conductivity accumulation function defined as:

[TABLE]

The AF model can also reproduce $\kappa$ for a-Si, but it is extremely sensitive to the empirical choice of the line broadening parameter ($\eta$). The impact of $\eta$ on the resulting $\kappa(T)$ is also shown in Figure 2, which shows that the value of $\kappa_{AF}$ varies by a factor two by varying $\eta$ between 0.01 and 0.5 meV in the temperature range considered. Whatever value is chosen for $\eta$, the AF model cannot reproduce the correct $\kappa_{QM}(T)$ of a-Si over the whole temperature range, in which we deem QHGK accurate ($T\leq 600K$), and it cannot give the correct decreasing trend at high temperature by construction. The predictions of the QHGK for the thermal conductivity of a-Si in the classical and fully quantum-mechanical regimes are compared in Figure S3 of SM.

In conclusion, we have introduced a unified approach to compute the lattice thermal con ductivity of both amorphous and crystalline systems. This quasi harmonic approach connects in a seamless fashion the AF model for disordered systems and the BTE-RTA for crystals. QHGK provides a significant improvement in generality over the Allen-Feldman model for disordered systems and is analytically proven to be equivalent to BTE for periodic systems. Classical QHGK calculations were validated against MD simulations for a-Si, and yield satisfactory agreement over a wide temperature range. Quantum QHGK can be deemed predictive at low temperature, not only for glasses and crystals but also for partially disordered systems, for which parameter-free models were up to now unavailable. The technique proposed in this work paves the way to robust calculations of heat transport in systems with any kind of structural order, including materials with point defects, extended defects and nanostructuring, without relying on any implicit knowledge of either their underlying symmetry, or the character of the vibrational modes, and without empirical parameters.