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.

Introduction

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)[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 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[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, kB 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, and local energies, , as  (ref. [3]), where in the harmonic approximation can be defined as: , M being the mass of i-th atom, its displacement from its equilibrium position, , α and β label Cartesian components, and is the IFC matrix. By expressing the energy flux in terms of the u’s, one obtains: , where . 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° would 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:where ω and γ are the harmonic frequency and decay rate of the n-th normal mode, and is the corresponding eigenvector of the dynamical matrix, , and indicates the ratio γ/ω, which vanishes in the harmonic limit. Equations (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: and , reading: . 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: , whose time evolution is . By doing so, the energy flux can be expressed in terms of the α amplitudes asBy 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, and . In the purely harmonic approximation, one would have . Anharmonic effects broaden the vibrational lines by a finite line-width, γ, which results in a finite imaginary part of the frequency and in a decay of the single-mode Green’s function, reading: . 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 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 is the Bose-Einstein occupation number of the l-th normal mode and 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 between normal-mode eigenvectors: , 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, and by an integer label, s, indicating different atomic sites within a unit cell, d: . 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, ν). In particular, the IFC matrix and its eigenvectors can be expressed as , where is the dynamical matrix of the system and its eigenvectors: and . When normal-mode eigenvectors are chosen to be real, the v matrices of Eq. (4) are real and anti-symmetric. In particular, and a non-vanishing 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: , where . 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 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 ∝ ω − ω. Hence in the periodic case the τ° matrix in Eq. (5) is strictly diagonal, , 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: , where 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: and , and being normal-mode creation/annihilation operators, satisfying the standard commutation relations for bosons: . Note that no ordering ambiguities arise when quantizing Eq. (6) because the 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: and , being the Bose-Einstein distribution function. The final quantum-mechanical expression for the heat conductivity in the QHGK is:with . For n = m this term reduces to the modal heat capacity . The other symbols are the same as in Eqs. (4) and (5) for the classical case. , in particular, is only different from zero for . 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 c = kB 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.

Fig. 1

Classical thermal conductivity. Comparison between the thermal conductivity of a-Si computed for a 1728-atom supercell by the classical Green–Kubo theory of linear response using either our QHGK approach (Eq. (3), green) or equilibrium molecular dynamics (purple). The vertical bars indicate statistical errors obtained by means of cepstral analysis, as explained in Ref. [22]. The k, k and k components of thermal conductivity tensor κ are averaged to obtain a value corresponding to an isotropic amorphous media

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 higher-order 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 off-diagonal terms of the matrix, weighted by (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:

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

Fig. 3

Thermal conductivity accumulation function. Conductivity accumulation function, κ(ω), and frequency-resolved differential conductivity, κ′(ω) (Eqs. (10) and (11)), computed for two different model sizes (N = 1,728 and N = 13,824 atoms) at temperatures T = 100 K and T = 600 K. Horizontal arrows in the upper panels indicate cumulative values of κ. κ′ is in units of WK−1 m−1 ps

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 κ 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 κ(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.

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 , second order force constants , and line-widths γ, computed at mechanical equilibrium (T = 0), it is possible to evaluate these same quantities through temperature-dependent statistical sampling approaches[29-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].