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

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.

H eat 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) 1,2 . 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 2,3 . 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 4 , 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 3,5-8 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 2 and a generalization of the Allen-Feldman (AF) model for disordered system 9,10 that explicitly and naturally accounts for normal-mode lifetimes. The main ingredients of our approach are the matrix of the interatomic 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 11 . Our theory is thoroughly validated in crystalline and amorphous silicon by comparing its predictions with those of EMD simulations and BTE computations.

Results
Theory. The basis of our work is the GK theory of linear response 3,5-8 , which relates the heat conductivity to the ensemble average of the heat-flux auto-correlation function: where V and T are the system volume and temperature, k B is the Boltzmann constant, J α (t) the α-th Cartesian component of the macroscopic heat flux, which in solids coincides with the energy flux, and 〈⋅〉 indicates a canonical average over initial conditions 3 . 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, R i , and local energies, the harmonic approximation ϵ i can be defined as: jβ iα u jβ , M i being the mass of i-th atom, is the IFC matrix. By expressing the energy flux in terms of the u's, one obtains: 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 12,13 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: J ← J°. Note that it is essential to use equilibrium atomic positions in the definition of J°, 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 J − J°would not be a total time derivative of a stationary process and the value of the transport coefficient resulting from J°w ould be biased. By making use of Newton's equation of motion, the final expression for the harmonic heat flux reads 9,10 : 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 14 . By doing so, the resulting time integral would diverge, thus yielding an infinite conductivity, as expected for a harmonic crystal 15 . 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: jγ iβ e iβ n e jγ m ; ð4Þ where ω n and γ n are the harmonic frequency and decay rate of the n-th normal mode, and e α ni is the corresponding eigenvector of the dynamical matrix, P and ϵ indicates the ratio γ/ω, which vanishes in the harmonic limit. Equations (3)(4)(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: 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: 2ω n p π n , whose time evolution is α n ðtÞ ¼ α n ð0Þe Àiω n t . By doing so, the energy flux can be expressed in terms of the α amplitudes as By using this expression, the integrand in Eq. (1) becomes a Gaussian integral of a fourth-order polynomial in the α's and α*'s that, by means of the Wick's theorem 14 , can be cast into a sum of products of pairs of single-mode (classical) Green's functions, hα Ã n ðtÞα m ð0Þi ¼ δ nm g n ðtÞ and hα n ðtÞα m ð0Þi ¼ 0. In the purely harmonic approximation, one would have g n ðtÞ ¼ k B T ω n e iω n t . Anharmonic effects broaden the vibrational lines by a finite line-width, γ 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Þ ¼ k B T ω n e iðω n þiγ 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)(4)(5) is presented in the Supplementary Note 1.
To lowest order in the anharmonic interactions, vibrational linewidths can be computed from the classical limit of the Fermi golden rule, where n l is the Bose-Einstein occupation number of the l-th normal mode and V′′′ nlm ¼ ∂ 3 V ∂ξ n ∂ξ l ∂ξ m is the third derivative of the potential energy with respect to the amplitude of the lattice distortion along the lattice normal modes 11 .
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 α matrices of Eq. (4) can be expressed in terms of the matrix elements of the matrices , where the notations "(e, e′)" and "V ⋅ e" indicate scalar and matrix-vector products in the space of atomic displacements. In crystals, equilibrium atomic positions are characterized by a discrete lattice position, a i , and by an integer label, s i , indicating different atomic sites within a unit cell, Likewise, in the Bloch representation, normal modes can be labeled by a quasi-discrete wavevector, q, belonging to the first Brillouin zone (BZ), and by a band index, ν: n → (q n , ν n ). In particular, the IFC matrix and its eigenvectors can be expressed as tβ ðqÞ is the dynamical matrix of the system and η sα νq its eigenvectors: e α ν ¼ e iq n ÁR i η s i α ν n q n and P tβ D sα tβ ðqÞη tβ νq ¼ ω 2 νq η sα νq . When normalmode eigenvectors are chosen to be real, the v α matrices of Eq. (4) are real and anti-symmetric. In particular, v α nn ¼ 0 and a nonvanishing thermal conductivity results from the matrix elements of v α connecting (quasi-) degenerate normal modes, i.e. modes whose frequencies coincide within the sum of their line widths. In the Bloch representation, v α is anti-Hermitian and block-diagonal with respect to the wave-vector, q. Its diagonal elements are imaginary, though not necessarily vanishing. In this representation one has: v α νq;μp ¼ i At fixed q, the vibrational spectrum is strictly discrete i.e. it remains so even in the thermodynamic limit. The number of q points for which there exists a pair of distinct modes, (q, ν) and (q, μ) with ν ≠ μ, that are degenerate within the sum of their line-widths ðjω qν À ω qμ j ≲ γ qν þ γ qμ Þ 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 νq,μq ∝ ω νq − ω μq . Hence in the periodic case the τ°m atrix in Eq. (5) is strictly diagonal, τ qν;pμ ¼ δ qp δ νμ τ qν , where is the anharmonic lifetime of the (q, ν) 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: ∂q α is the group velocity of the ν-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 1 . 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 Supplementary Fig. 1).
Moving to the quantum regime is straightforward in our approach. To this end, we start from the quantum GK formula 3,5-8 , reading: whereĴ α are quantum heat-flux operators and 〈⋅〉 indicates quantum canonical averages. A quantum expression for the heat flux is obtained from its classical counterpart, Eq. (6), by making the substitutions: α ! ffiffi ffi h pâ and α Ã ! ffiffi ffi h pâ y ,â y andâ being normal-mode creation/annihilation operators, satisfying the standard commutation relations for bosons:â;â y Â Ã ¼ 1. Note that no ordering ambiguities arise when quantizing Eq. (6) because the v α nm matrices are antisymmetric, and they therefore vanish for n = m. The resulting expression for the quantum heat flux is:Ĵ 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: hâ y n ðtÞâ n ð0Þi ¼ n n e iðω n þiγ n Þt and hâ n ðtÞâ y n ð0Þi ¼ ðn n þ 1Þe Àiðω n Àiγ n Þt , n n ¼ 1= e hω n k B T À 1 being the Bose-Einstein distribution function. The final quantum-mechanical expression for the heat conductivity in the QHGK is:  4) and (5) for the classical case. τ nm , in particular, is only different from zero for jω n À ω m j ≲ γ n þ γ 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→∞ c nn = k B and the quantum formula, Eq. (9), reduces to Eq. (3). Further details are given in Supplementary Note 2.
Simulations. 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 16 , which describes well the thermal conductivity of bulk and nanostructured silicon, including a-Si 9,10,17-19 . 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 20,21 . The integral of the heat flux autocorrelation function in Eq. (1) can then be efficiently evaluated via cepstral analysis, as described in refs. 3,22 , which can be enhanced by averaging over multiple trajectories at low temperature (T ≤ 300 K). Details on the data analysis procedure followed here and on the estimate of the statistical errors is given in the Supplementary Note 3. The results of these calculations are reported in Fig. 1 and exhibit a weak non-monotonic dependence on T. Performing similar MD simulations on models of increasing size (4096 and 13,824 atoms) generated with the same protocol, we have verified that size effects on κ at 300 K amount to less than 15%, which is of the same order as the variation κ 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 ( Supplementary Fig. 2). The resulting strongly diagonally dominant form of the τ°matrices in Eq. (5) is also displayed in Supplementary Fig. 2.
The thermal conductivity obtained by QHGK is in excellent agreement with that computed by EMD for T ≤ 600 K (Fig. 1). At higher temperatures QHGK overestimates κ, as it neglects higherorder anharmonic effects. We point out that, in spite of the formal analogies with the AF model 9,10,23 and recent refinements thereof 24,25 , Eq. (3) entails no empirical parameters. It thus allows one to predict temperature trends dictated by anharmonic effects in good agreement with MD, without making any a priori distinction among propagating, diffusive, or localized vibrational modes. Conversely, in the AF model the temperature dependence lies only in the heat capacity term, therefore in the classical limit κ is temperature independent.
Similarly to the GK modal analysis approach 26 , based on classical MD, the transport character of the modes is dictated by the relative contribution from the diagonal and slightly offdiagonal terms of the v α nm matrix, weighted by τ nm (Supplementary Fig. 2). 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 Fig. 2 we report our results from quantum QHGK calculation for an amorphous Si model of 13824 atoms along with three sets of experimental data 27,28 . QHGK results are in excellent agreement with the measurements in Ref. 27 above 100 K. At lower temperature the estimate of κ 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). Specifically, we see a significant improvement in the estimate of κ at 50 K from the 1728-atom model (κ = 0.027 Wm −1 K −1 ) to the 13824 model (κ = 0.25 Wm −1 K −1 see Fig. 3). However, at 50 K and lower temperatures the latter is not converged yet. 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 Fig. 3, where we compare for different temperatures and model sizes the frequency-resolved differential conductivity, where Δ(ω) is a broadened approximation of the δ function and the other symbols have the same meaning as in Eq. (9), and the conductivity accumulation function defined as: The AF model can also reproduce κ for a-Si, but it is extremely sensitive to the empirical choice of the line broadening parameter (η). The impact of η on the resulting κ(T) is also shown in Fig. 2, which shows that the value of κ AF varies by a factor two by varying η between 0.01 and 0.5 meV in the temperature range considered. Whatever value is chosen for η, the AF model cannot reproduce the correct κ QM (T) of a-Si over the whole temperature range, in which we deem QHGK accurate (T ≤ 600 K), 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 Supplementary Fig. 3.  Fig. 2 Quantum thermal conductivity. Thermal conductivity computed for a 13824-atom supercell of a-Si using the quantum QHGK approach in the quantum regime (Eq. (9)), compared with the Allen-Feldman approach 9,10,23 and experimental data (yellow triangles and orange diamonds ref. 27 ), (green triangles ref. 28 ). The broadening η used in Allen-Feldman calculations is set equal for every normal mode Conclusions. In conclusion, we have introduced a unified approach to compute the lattice thermal conductivity 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. Although the numerical results of this work were obtained by evaluating Eqs. (3)-(5) using equilibrium positions R iα , second order force constants Φ jγ iβ , and line-widths γ n , computed at mechanical equilibrium (T = 0), it is possible to evaluate these same quantities through temperature-dependent statistical sampling approaches [29][30][31] , thus extending the reach of QHGK to systems with strong anharmonicity and high-temperature phases. 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. While this paper was being written we learnt that conclusions similar to ours were reached by Simoncelli et al., following a different approach based on a generalization of the BTE 32 .

Data availability
The data that support the plots within this paper are available from the corresponding author upon reasonable request.

Code availability
Computer codes are available from the corresponding author upon reasonable request.