Heat conductivity from energy-density fluctuations

We present a method, based on the classical Green-Kubo theory of linear response, to compute the heat conductivity of extended systems, leveraging energy-density, rather than energy-current, fluctuations, thus avoiding the need to devise an analytical expression for the macroscopic energy flux. The implementation of this method requires the evaluation of the long-wavelength and low-frequency limits of a suitably defined correlation function, which we perform using a combination of recently-introduced cepstral-analysis and Bayesian extrapolation techniques. Our methodology is demonstrated against standard current-based Green-Kubo results for liquid argon and water, and compared with a recently proposed similar technique, which utilizes mass-density, instead of energy-density, fluctuations.


I. INTRODUCTION
Heat transport is one of the most fundamental offequilibrium processes occurring in nature, whose understanding sinks its roots into the classic works of Onsager in the thirties 1 and of Green and Kubo (GK) 2,3 in the fifties.In spite of the strong theoretical foundation that these works provide to the numerical simulation of transport phenomena in condensed matter, until recently widespread misconceptions had hindered the application of standard equilibrium simulation techniques to the numerical modeling of heat transport when atomic forces are derived from quantum-mechanical ab initio methods. 4While the crux of the difficulties was originally ascribed to the impossibility of assigning a well-defined value to the atomic energies that enter the definition of the macroscopic energy current-when atomic forces are derived from quantum mechanical methods-it was later realized that a unique partition of the total energy into local contributions is not possible even when atomic forces are derived from classical potentials.It was then demonstrated that, although the macroscopic current corresponding to different local partitions of the total energy depends on the specific partition being considered, the value of the resulting heat conductivity does not.This independence of the heat conductivity of the local representation of the total energies was dubbed gauge invariance of transport coefficients [5][6][7][8] and has far-reaching consequences in other branches of transport theory as well, such as, e.g., charge transport in ionic conductors. 9,10lthough gauge invariance provides a general framework to compute heat transport coefficients for general atomic forces derived from either electronic-structure theory or arbitrary-possibly machine-trained-interatomic potentials, the problem still remains that the expression for the macroscopic energy current for general classical force fields or quantum density functionals may be difficult to compute analytically and to implement numerically.For this reason, it would be desirable to devise a) Electronic mail: endrigo@sissa.it a method to compute the heat conductivity not relying on the macroscopic energy current, which is the standard ingredient of the GK formula for the transport coefficient, but rather on the energy density, which-albeit not uniquely defined because of gauge freedom-is easier to compute.A successful attempt in this direction has been recently made by Cheng and Frenkel (CF) 11 who have shown how to evaluate the heat conductivity of simple fluids, κ, by analyzing the wave-vector dependence of the widths of the Rayleigh (R) and Brillouin (B) peaks, as estimated from equilibrium molecular dynamics (EMD) simulations: [12][13][14] where D T = κ mncp and D ℓ are the heat and longitudinal diffusivities, respectively, R LP = cp cv − 1 is the (Landau-Planczek) ratio between the intensities of the Rayleigh and Brillouin peaks , which vanishes for incompressible fluids, c p and c v are the isobaric and isochoric specific heats, and m and n the molecular mass and average number density, respectively.The smallness of the Landau-Planczek ratio in the condensed phases (R LP < 0.01 in liquid water at ambient conditions), 15 however, makes the application of this method numerically demanding in liquids and practically prohibitive in solids.
In a pioneering, but not highly cited, 1994 paper, Palmer showed how the heat conductivity can be determined by a direct analysis of the energy-density fluctuations in EMD simulations. 16This method does not rely on the value of the Landau-Planczek ratio, and therefore can be equally applied to liquids and solids.However, its broad applicability is also plagued by the numerical and statistical problems related to the evaluation and long-wavelength extrapolation of the peak of the energy dynamical structure factor (see below for a precise definition) and, on a more conceptual level, by the intrinsic indeterminacy of the energy density at the molecular scale, discussed above.
Building on the work of Palmer, 16 in this paper we first reexamine the concept of gauge invariance in the specific context of the energy-density approach to heat transport, and then we show that a combination of recent advances in data analysis and signal processing, based on Bayesian inference 17 and cepstral analysis, 18,19 can be leveraged to overcome the numerical difficulties that have plagued the application of Palmer's approach.We argue that the proposed methodology naturally lends itself to a quantitative evaluation of the resulting statistical uncertainties and can be applied using simulation cells of moderate size and EMD runs of moderate length; we finally showcase it with three applications to the paradigmatic cases of liquid water and Argon, and solid amorphous Silica.

II. THEORY
The standard GK expression for the heat conductivity, κ, reads: [12][13][14] In Eq. (3) V and T are the system's volume and temperature, k B the Boltzmann's constant, J is any Cartesian component of the macroscopic average of the heat current density, J = 1 V V ȷq (r)dr, a hat indicates an implicit dependence on phase-space variables, Γ-as in J = J(Γ)and ⟨•⟩ indicates an equilibrium average over the initial conditions of a molecular trajectory.The heat current density is defined as: 12-14 ȷq = ȷe − (p + e)ȷ n /n, where ȷe and ȷn are the energy and molecular-number density currents, respectively, and p, e are the system's pressure and average energy, respectively.The explicit expressions of the energy and energy-current densities depend on the specific simulation framework (e.g.force-field vs. ab initio), as well as on gauge fixing (see below).A detailed derivation of the relevant formulas is presented, e.g., in Ref. 6.
In a one-component system, the macroscopic average of the molecular number current density is proportional to the total momentum, which, being conserved, can be assumed to vanish.As in this work we will only consider this class of systems, we will not make any difference between heat and energy densities, and we will therefore omit any suffixes aimed at distinguishing one from the other.In Eq. (4) χ ′′ is the imaginary part of the heat-temperature susceptivity, which, according to the fluctuation-dissipation theorem, is proportional to the energy dynamical structure factor (EDSF), C(k, ω): where e(k) is the space Fourier transform of the energy density, e(r).The equivalence between Eqs. ( 3) and ( 4) follows from the fluctuation-dissipation theorem, Eq. ( 5), and the continuity equation, which, in the frequency-wavenumber domain, reads ω e(k, ω) = k • ȷ(k, ω).Thanks to this relation, the density-density and current-current correlation functions can be expressed in terms of each other as: . By combining the Onsager's regression hypothesis 1 with a hydrodynamic description of the long-time evolution of the long-wavelength components of the energy-density fluctuations, the small-frequency/smallwavevector limit of the EDSF can be shown to read: 13,14 Inspection of Eq. ( 7) reveals that for large wavelengths the static limit of the EDSF reads: We conclude that the heat conductivity can be equivalently expressed as: Note the non-commutativity of the limits in Eq. ( 4), manifested by Eq. ( 7).This well-known singularity, which has been extensively discussed in the literature following Kadanoff and Martin's landmark paper, 12 should not come as a surprise, for the zero-frequency limit of a generic response function is a static susceptibility.It is an obvious consequence of time-reversal invariance that, strictly speaking, a static perturbation cannot induce a stationary current.The interpretation of this equation when the limits are performed in the correct order is as follows: While the long-time limit of the current induced by a temperature fluctuation vanishes at all finite wavelengths, as the system approaches equilibrium, at any finite time, the current tends to a finite, time-independent, value as the wavelength of the perturbation grows large.This behavior defines a stationary current in response to the temperature fluctuations.

A. Gauge invariance
The energy density, e(r), entering Eq. ( 6) is intrinsically ill-defined, as it is affected by a so-called gauge freedom, deriving from the insensitivity of the total energyand therefore of all the macroscopic thermal propertieson the addition of the divergence of a bounded vector field to the energy density: e(r) → e(r) + ∇ • p(r).The independence of the heat conductivity on the specific representation of the energy density is discussed, e.g., in Refs., 5,7,8,18 and demonstrated below in the context of the present approach.
In order to prove the gauge invariance of Eq. ( 9), we examine how a gauge transformation, e ′ (r) = e(r) + ∇ • p(r), affects the zero-frequency value of the EDSF, C(k, ω = 0).Under such a transformation, the integrand in Eq. ( 6) would go into: According to the Wiener-Kintchine theorem, 20,21 the Fourier transform of the time auto-correlation function of a zero-mean stationary stochastic process, X(t), in the large-time limit is proportional to the variance of the Fourier transform of the process, The ω = 0 limit of this relation reads: When X(t) is a conserved flux, the left-hand side of Eq. ( 13) is the Green-Kubo expression of a transport coefficient.This equation is a generalization of Einstein's celebrated relation between a particle's diffusivity and its velocity auto-correlation function. 22Einstein's formula was extended to generic transport coefficients by E. Helfand in 1960. 23ccording to these considerations, the ω = 0 value of the transformed EDSF can be computed as The second term of Eq. ( 14) is bounded from above by and therefore, since p(r) has a finite variance, it vanishes in the k → 0 limit.Since in k → 0 the limit both e(k, t) and p(−k, t ′ ) are real, the third term in Eq. ( 14) vanishes as well.While at finite k it is to be expected that κ(k) in Eq. ( 10) may be affected by the specific choices of p(r), its long-wavelength limit does not depend on it, thus providing an independent demonstration of the gauge invariance of the heat conductivity, particularly adapted to the present framework.

III. COMPUTER SIMULATIONS
The method embodied in Eq. ( 9) has been benchmarked against standard GK calculations based on Eq. (3) in the cases of liquid water, liquid Argon and solid amorphous Silica.Water was modeled with a sample of 216 molecules interacting through the SPC/E force field 24 at T = 300 K and p = 0 bar.Argon was mod-eled with a sample of 864 atoms interacting through a Lennard-Jones potential 25 at T = 100 K and p = 0 bar.Amorphous Silica was modeled with a sample of 648 atoms interacting via the BKS potential 26 at 500 K.Long EMD simulations were run in the micro-canonical (NVE) ensemble after careful equilibration performed in a number of different ensembles.In all cases periodic boundary conditions (PBC) were used.All the technical details of the simulations are reported in Appendix A. In order to estimate the ω = 0 value of the ESDF in Eq. ( 10), the time series of the Fourier transform of the energy density, ẽ(k ℓ , t) = n ϵ n (t)e −ik ℓ •Rn(t) , was first collected along the EMD trajectory at regular time steps.In the preceding expression ϵ n (t) and R n (t) are the local energy and position of n-th atom, respectively, at time t.Leveraging the Wiener-Kintchine theorem, the ESDF is then evaluated as: where τ is the length of the EMD run and the integral is actually evaluated as a discrete Fourier transform, as explained, e.g., in Refs.18 and 19.In the signal-processing literature, a sample of a power spectral density (PSD) drawn from a finite and discrete time series is often referred to as a periodogram.In Eq. ( 16) and in the rest of this section a hat indicates that the quantity it is put on (the periodogram and the energy density, in the present case) is actually a sample drawn from a stochastic process.It can be demonstrated that the periodogram constitutes a set of stochastic variables that are uncorrelated for different frequencies and that, for a given frequency, ω j = 2πj τ , are distributed as: where the ξ's are independent χ 2 2 variates.From Eq. ( 17) it follows that the (energy-density) periodogram is an unbiased estimator of the PSD (the ESDF, in this case), C(k, ω j ) = C(k, ω j ).This estimator is not consistent, though, because the variance of any of its values is independent of the length of the EMD run.In order to get such a consistent estimator (at the price of introducing some bias), it is expedient to make use of the so called cepstral analysis method. 18,19,27The cepstrum is defined as the (inverse) discrete Fourier transform of the logarithm of the periodogram.The logarithm turns the multiplicative noise affecting the periodogram into an additive noise affecting the cepstrum; furthermore, the PSD, as well as its logarithm, are much smoother than the periodogram and the cepstrum itself.Therefore, a consistent estimate of the log-PSD can be obtained by applying a low-pass filter to the cepstrum.The optimal value of the low-pass quefrency 27 -i.e. the number of inverse Fourier coefficients of the cepstrum-is determined in such a way that the estimated bias introduced by the filter is of the order of the resulting statistical noise.This number can be obtained by any model-selection method, and in our applications we use the Akaike information criterion. 18,19,28ith the aim of evaluating Eq. ( 9), the EDSF was first sampled on a homogeneous grid of wave-vectors, {k ℓ }, with magnitude smaller than a preassigned cutoff, |k ℓ | ≤ k max , compatible with the chosen PBCs.Because of spherical symmetry, the EDSF only depends on the modulus of its wave-vector argument and can therefore be averaged accordingly.In practice, PBCs reduce the symmetry from spherical to cubic, thus making some PBC-allowed wave-vectors non-equivalent to others with the same magnitude.This artifact is eliminated by averaging the estimates of the EDSF corresponding to different wave-vectors of the same magnitude, thus further reducing the variance of any individual value.The set of wave-vectors with a same magnitude, k i , will be referred to as the k i -star of wave-vectors.In Fig. 1 we display the EDSF of liquid Argon at the the kstar of minimum magnitude compatible with the chosen PBCs, k min = 1.77 nm −1 , and the one of largest magnitude utilized in our Bayesian analysis (see below), k max = 11.3 nm −1 .The increasing sharpness of the peaks for increasing wavelength (decreasing wave-vector) is a signature of their hydrodynamic nature, as manifested by Eq. ( 7), and makes cepstral analysis increasingly difficult to perform, thus giving rise to a larger statistical incertitude.Because of isotropy, the wave-vector dependence of the conductivity, Eq. (10), is an even function of the squared modulus of k, which we assume to be well approximated by a polynomial of order M : In order to determine the coefficients of the polynomial fit in Eq. ( 18), w = {w 0 • • • w M }, we relied on a Bayesian inference method, by seeking to maximize their posterior probability, given the occurrence of the measured data, D: p (w|D) = p (D|w) p(w)/p(D), where p(D) and p(w) are the data and model-parameter priors, respectively.The data are the conductivities measured at each k-star, which are independent normal variates whose expectations, κ i , and variances, σ 2 i are estimated by the cepstral analysis procedure outlined above.The normality of the distribution of the measured values of κ(k) is showcased in Appendix B. The probability of the data conditional to the values of the coefficients of the polynomial is thus given by p (D|w) ∝ e −L(w) , where L(w) = i , is the likelihood function and Φ i is an array of monomials of even degree up to order 2M , evaluated at the k-stars, k i , for which the k-dependent conductivity has been estimated: We now make the further assumption that the prior distribution of the parameters is normal, p(w) = α 2π M 2 exp −α∥w∥ 2 , thus allowing the extrapolation to prevent over-fitting, which, in ordinary linear regression, would be equivalent to introducing a regularization term, α∥w∥ 2 .By leveraging again Bayes theorem, the optimal number of parameters, M , and regularization parameter, α, are determined as those that maximize the probability of their occurrence, conditionally to actual observation of the data set: The actual procedure that we followed to determine M , α, w, and, therefore, κ = w 0 is described in Refs.17 and 30.
In Fig. 2 we display the wave-vector dependence of Eq. ( 10), κ(k), computed for SPC/E water, liquid Argon and solid amorphous Silica via the cepstral analysis, and the k → 0 extrapolation of Eq. ( 9), κ, obtained from the Bayesian regression method.The results are in good accordance with the GK integral estimated via the cepstral analysis on the same trajectory, also reported in the figures.In Appendix C we study the finite size effects on the computation of Eq. (10).We studied the size convergence of Eq. ( 10) for SPC/E water modeled with samples of different linear box size, L: 18.5 Å, 31.1 Å and 37.4 Å; corresponding to 216, 1000 and 1728 molecules respectively.
In Fig. 3 we display the Bayesian extrapolation (k = 0) of the heat conductivity of liquid Argon as a function of the maximum magnitude of the wave-vectors considered in the fit, k c .The Bayesian extrapolation is robust and converges to the GK result computed on the same trajectory, also reported in the figures, when we enlarge the data set, thus demonstrating the numerical stability and statistical consistency of our method.

IV. CONCLUSIONS
In this work we have integrated a method originally proposed by Palmer in 1994, 16 to compute the heat conductivity in extended systems from energy-density fluctuations, with modern data and signal analysis techniques that considerably enhance its applicability.As a downside, as already anticipated by Palmer, the numerical efficiency of the methodology presented here appears to be poorer than that of the standard GK procedure based on the fluctuations of the energy flux.In the specific case of Argon, we estimated that, in order to achieve the same statistical accuracy of the standard GK procedure, one would need EMD trajectories one-order of magnitude as long.The enhanced statistical performance of currentbased methods for estimating transport coefficients, as  10), of SPC/E water at T = 300 K and 0 bar, as evaluated from the EMD simulation, kmin = 3.40 nm −1 .Central panel: same for liquid Argon at T = 100 K and 0 bar, kmin = 1.77 nm −1 .Lower panel: same for solid amorphous Silica at T = 500 K, kmin = 3.04 nm −1 .The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use.Cubic average refers to an average performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry.The red line indicates the Bayesian fit with the associated estimated uncertainty.At k = 0 we report the Bayesian extrapolation and the GK estimate, on the same EMD trajectory, of the heat conductivity.opposed to density-based methods, can be ascribed to the necessity of performing long-wavelength extrapolations when utilizing the latter.Avoiding any such extrapolation would significantly enhance the performance of any statistical analysis, no matter how smart this may be.Of course, this numerical superiority may, or may not, be tarnished by the need to devise and implement an explicit expression for the macroscopic energy flux.
In particular, it may be challenging, if not impossible, to devise an analytical form of the energy flux when ad- vanced energy functionals, such as hybrid or meta-GGA, are used for ab initio MD.We believe that a densitybased approach, such as the one examined in the present paper or in CF's, 11 would be mandatory in these cases.The relevance of these cases is however increasingly limited by the emergence of artificial-intelligence methods to train classical force fields with quantum mechanical accuracy, 31 which allow for a fast and accurate evaluation of various transport coefficients using the standard GK approach, including the heat conductivity. 32,33Furthermore, all the density-based methods proposed so far are at present limited to one-component systems, whereas a systematic extension of the standard GK one has been recently established. 34In summary, we conclude that standard current-based methods are preferable in all those cases where an explicit expression for the energy flux is available.In the rare cases where it is not (e.g. for ab initio MD with advanced functionals), density-based methods may be the way to go.The method discussed here and the one proposed by CF 11 are conceptually similar, but the latter is in practice only applicable when the Landau-Planczek ratio is sizeable, i.e. for compressible fluids, whereas the former can be applied to solids and incompressible fluids as well.We also hint at the superiority of the data analysis technique introduced in Ref. 30 and utilized here, which allows for an easy control and optimal estimate of the statistical accuracy of the transport coefficients.This technique would readily integrate into CF's approach, 11 as well as in any methodology relying on the long-wavelength extrapolation of correlation functions estimated for finite, periodic models.through an iter-atomic Lennard-Jones potential whose parameters are fitted to experiments. 25After a preliminary equilibration in the NPT ensemble, we sampled the system from a 2-ns NVE trajectory using a time step of 1 fs.In this case, the energy density was computed and dumped at intervals of 5 fs.Solid amorphous Silica was simulated at 500 K using a sample of 648 atoms modeled via the BKS potential. 26The initial configuration was taken from Ref. 36 and optimized so as to make atomic forces smaller than a preassigned threshold of 10 −8 eV/ Å.After careful equilibration for 125 ps in the NVT ensemble, we simulated the system for 1 ns in the NVE ensemble with a timestep of 0.25 fs.The energy density was computed and dumped at intervals of 5 fs.All simulations were performed with the parallel code LAMMPS 37,38 using the velocity-Verlet algorithm.10), of SPC/E water at 300 K and 0 bar, as evaluated from the EMD simulations for three different linear box sizes, L: 18.5 Å, 31.1 Å and 37.4 Å; corresponding to 216, 1000 and 1728 water molecules respectively.The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use.Cubic averages, represented as crosses, refer to the averages performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry.At k = 0 we report, as dots, the Bayesian extrapolations of the heat conductivity at each size.

FIG. 1 .
FIG. 1. Upper panel: EDSF of liquid Argon computed at PBC-compatible k-vectors of minimum magnitude, kmin = 1.77 nm −1 .Lower panel: same, but for the k-vectors, kmax = 11.3 nm −1 .The window-filtered line refers to a moving average 29 performed width a window with of ∆ = 0.2 THz.The shaded red area indicates the statistical incertitude.

2 FIG. 2 .
FIG.2.Upper panel: Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar, as evaluated from the EMD simulation, kmin = 3.40 nm −1 .Central panel: same for liquid Argon at T = 100 K and 0 bar, kmin = 1.77 nm −1 .Lower panel: same for solid amorphous Silica at T = 500 K, kmin = 3.04 nm −1 .The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use.Cubic average refers to an average performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry.The red line indicates the Bayesian fit with the associated estimated uncertainty.At k = 0 we report the Bayesian extrapolation and the GK estimate, on the same EMD trajectory, of the heat conductivity.

FIG. 3 .
FIG.3.Upper panel: Bayesian prediction of liquid Argon at 100K and 0 bar, at k = 0 as a function of the maximum kvector included in the fit, kc, confronted with the GK integral computed on the same trajectory with the cepstral analysis.Lower panel: Number of parameters of the Bayesian regression fit as a function of the maximum k-vector included in the fit.

1 FIG. 4 .
FIG.4.Distribution of the wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300K and 0 bar as evaluated from multiple 200ps-long EMD segments extracted from a 20ns long trajectory and the pvalue of the Shapiro-Wilk normality test.The solid lines represent a Gaussian curve fitted on the distribution of the heat conductivities at the wavevector presented.

FIG. 5 .
FIG.5.Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at 300 K and 0 bar, as evaluated from the EMD simulations for three different linear box sizes, L: 18.5 Å, 31.1 Å and 37.4 Å; corresponding to 216, 1000 and 1728 water molecules respectively.The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use.Cubic averages, represented as crosses, refer to the averages performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry.At k = 0 we report, as dots, the Bayesian extrapolations of the heat conductivity at each size.