Thermal transport of glasses via machine learning driven simulations

Accessing the thermal transport properties of glasses is a major issue for the design of production strategies of glass industry, as well as for the plethora of applications and devices where glasses are employed. From the computational standpoint, the chemical and morphological complexity of glasses calls for atomistic simulations where the interatomic potentials are able to capture the variety of local environments, composition, and (dis)order that typically characterize glassy phases. Machine-learning potentials (MLPs) are emerging as a valid alternative to computationally expensive ab initio simulations, inevitably run on very small samples which cannot account for disorder at different scales, as well as to empirical force fields, fast but often reliable only in a narrow portion of the thermodynamic and composition phase diagrams. In this article, we make the point on the use of MLPs to compute the thermal conductivity of glasses, through a review of recent theoretical and computational tools and a series of numerical applications on vitreous silica and vitreous silicon, both pure and intercalated with lithium.

The pursuit of improving the thermal conductivity properties of amorphous solids is central to contemporary materials science and engineering [1].Glasses, characterized by their lack of crystalline order, possess unique attributes that make them invaluable across a wide range of applications.One of the prominent features of this class of materials is a inherently low thermal conductivity, a result of their disordered structure.This property is useful in various fields, such as aerospace engineering [2,3], electronics [4], and pharmaceutical industries [5].In contrast, specific industrial applications demand a nuanced balance in the thermal properties of glasses.For example, nuclear reactors and nuclear weapon decommissioning generate radioactive waste [6] that must be safely stored for exceptionally long time.This waste can be solidified through vitrification, preventing accidental radionuclide release thanks to the amorphous structure of glasses, which provides radiation protection and outstanding chemical durability, thus enabling thousands of years of safe storage [7].Here, effective heat management is crucial, as high thermal conductivity enhances vitrification efficiency, influencing melt rate and glass homogeneity [8,9].Moreover, in long-term storage, elevated heat conductivity rapidly dissipates decaygenerated heat, avoiding issues like high-temperatureassisted crystallization, porosity, and cracks [10].During the last decades, a huge effort has been put forward to access the structural and thermodynamical properties of the glasses employed for nuclear waste vitrification, mainly (boro)silicates, from both the experimental [8][9][10][11] and the computational sides [12][13][14][15][16][17][18].Nevertheless, a microscopic description of thermal conduction in these materials beyond the celebrated Cahill-Pohl model [19], as a function of temperature and composition, is still missing.The lack of computational studies on thermal conduction is shared by another important application, namely solid-state batteries, where several designs leverage amorphous solid electrolytes [20][21][22][23][24]. Here, the glassy phase is also characterized by a diffusive species (like Li + or Na + ions) which poses further challenges for the microscopic simulation of heat transport, since lattice methods cannot be formally applied due to the lack of well defined positions of mechanical equilibrium.
Machine-learning (ML) is an increasingly popular tool in materials modeling due to its ability to train on extensive datasets precise models able to match a wide array of target properties.One route is to leverage datasets of mechanical and functional glass properties to discover new materials [25,26], or to predict end-properties such as solubilities [27], dissolution rates [28], or transition temperatures [29], to name a few [30].By taking a microscopic approach, one can build surrogate models able of ab initio interatomic potentials to drive simulations to sample all kinds of properties of amorphous materials that can be accessed by atomistic simulations [31][32][33][34][35][36][37][38][39][40][41][42][43].
In this article, we provide an overview of the microscopic theory of heat transport in glasses.Our focus is on methods that effectively utilize ML potentials (MLPs) for computing thermal conductivity.We explore two primary approaches: lattice dynamics, suitable for solids at temperatures significantly below their melting point, enabling the incorporation of quantum-mechanical effects in heat transport; and equilibrium molecular dynamics (EMD) simulations, a versatile tool for sampling material properties, which however is limited in its ability to account for the aforementioned quantum effects.

A. Thermal transport in glasses
Heat transport is characterized by the value of the thermal conductivity, κ, whose linear-response value is given by the Green-Kubo (GK) formula [44][45][46][47] arXiv:2402.06479v1[cond-mat.dis-nn]9 Feb 2024 where Ω is the system's volume, k B the Boltzmann constant, T the temperature, and J q the heat flux.The factor 1/3 comes from the assumption of isotropy.The temperature-dependent behavior of κ in glasses exhibits three distinctive, universally recognized patterns [48].At extremely low temperatures, specifically when T ≲ 2 K, the predominant scattering mechanism involves quantum tunneling between various local minima within the glass energy landscape, resulting in κ ∼ T 2 [49][50][51].As the temperature reaches a few tens of kelvins, thermal conductivity increases and eventually reaches a plateau value.Despite the absence of a firmly established theoretical consensus in the literature, this phenomenon appears to be linked to the transition from a regime dominated by quantum processes to one where propagating waves are scattered by random disorder [48,[50][51][52].Beyond the plateau, the behavior of the thermal conductivity is governed by the anharmonic decay of normal modes [53,54], and κ starts increasing again until it saturates to its high-temperature value.We focus here on the latter range of temperatures, which is relevant for applications and can be investigated by means of atomistic simulations.
There are different techniques for computing κ.Equilibrium methods require the calculation of the heat flux appearing in Eq. ( 1) [47], or the energy density [55].Equivalently, one can compute the energy flux along with other relevant mass fluxes, as discussed below.A versatile approach for obtaining these fluxes is through EMD simulations, which allow to sample the energy flux at all orders within anharmonic perturbation theory.A downside of EMD simulations is that it is currently impossible to sample the energy flux including nuclear quantum effects [56], despite extensive efforts are being made in this direction [see, e.g.[57][58][59] and citations therein].In solids, an alternative approach lies in exploiting lattice dynamics, at the cost of neglecting high-order anharmonic interactions.At temperatures well below the melting point, atomic nuclei undergo small oscillations around well-defined equilibrium positions.This behavior allows us to represent the dynamics with harmonic normal modes, which, in crystals, are described as phonon quasiparticles.In the harmonic approximation, the normal modes feature infinite lifetimes, resulting in infinite thermal conductivity, regardless of the presence of (harmonic) perturbations, such as disorder [60].However, when anharmonic interactions are appropriately accounted for, e.g., through perturbation theory [53,54], they induce temperature-dependent frequency shifts and broadenings, impacting phonon lifetimes and ultimately determining the thermal conductivity of materials.
Using lattice dynamics to express the heat flux in Eq. (1) in terms of phonon creation and annihilation operators, one can compute the GK formula for κ in the quasi-harmonic Green-Kubo (QHGK) approximation [53,61] as where ω µ and γ µ are the angular frequency and anharmonic linewidth of the µth normal mode, respectively; v µν is a generalized velocity matrix, and is a generalized two-mode isochoric heat capacity, n(ω) = [e ℏω/(kBT ) − 1] −1 being the Bose-Einstein (BE) occupation function.
The QHGK approach employs two interconnected approximations [61,62].The first is the dressed bubble approximation, where four-point correlation functions among phonon creation and annihilation operators are factorized into products of two-point correlation functions, neglecting vertex corrections [61].This implies that normal modes decay independently, each interacting with a common heat bath.
The second approximation, termed Markovian, disregards memory effects in the heat bath-normal mode interaction [61].The combination of these approximations leads to four-point correlation functions being expressed using single-body greater Green's functions, denoted by g > µ (t) = −i(n µ + 1)e −iωµt−γµ|t| .The quasi-harmonic hypothesis requires γ 2 µ /ω 2 µ ≪ 1, implying that only nearlydegenerate modes with |ω µ − ω ν | ≲ γ µ + γ ν significantly contribute to heat conductivity.The QHGK approximation works on crystals and glasses alike, reducing to the result of the Boltzmann Transport Equation in the former case, and-at the expense of neglecting anharmonic effects-to the Allen-Feldman (AF) model of harmonic disordered solids in the latter [53,61,63].Neglecting anharmonic effects is no trivial matter-it remarkably transforms the finite bulk thermal conductivity of a glass into an infinite quantity, thereby emphasizing the AF calculations' qualitative accuracy driven solely by size effects [60,64].
If compared with crystals, glasses present an additional complexity in numerical computations due to their inherent aperiodic nature [64].The customary application of periodic boundary conditions (PBC) to finite simulation cells, from which quantities in Eq. ( 2) are derived, requires cautious consideration of size effects.While the physical symmetry of crystals allows for calculations on a fine mesh in reciprocal space, in glasses the same approach is not feasible without introducing spurious contributions to the thermal conductivity [64].Nonetheless, size extrapolation remains possible by leveraging the asymptotic Debye expression of the thermal conductivity of propagons-the low-frequency normal modes in glasses characterized by wave-like properties, including well-defined dispersion with small broadening [64,65].
Using lattice dynamics requires to be able to compute second and third-order interatomic force constants, i.e., the second and third partial derivatives of the potential energy with respect to nuclear coordinates at equilibrium [63].This can be achieved through finite differences [66,67] or perturbation theory [68], with the former method being more common.For large systems (i.e., containing tens of thousands of atoms) calculating interatomic force constants is computationally feasible using empirical force fields, but becomes impractical when pursued ab initio due to the N 2 atoms scaling of second and N 3 atoms scaling of third-order force constants.In principle, MLPs offer a partial remedy, significantly expediting interatomic force computations.Regrettably, practical implementation is often spoiled by pronounced numerical noise in interatomic force constants, impeding the computation of anharmonic linewidths through Fermi's golden rule [63], as it will be showcased below by a toy example involving vitreous silica.Inaccurate anharmonic linewidths severly hinder the computation of heat conductivities, thereby rendering the entire lattice dynamical workflow unfeasible.This concern is particularly relevant at low frequencies [69], even for MLPs with good performances within the medium to high-frequency range.A potential strategy to mitigate this issue entails substituting finite-difference third-order derivatives with analytical ones, exploiting the differentiabilty of MLP descriptors through, e.g., automatic differentiation [43].
Therefore, while the challenge of deriving interatomic force constants from MLPs is under scrutiny, their accuracy in reproducing ab initio results can be leveraged in MD simulations through the GK equation, Eq. (1).

B. Thermal conductivity from molecular dynamics simulations
The GK formula requires an expression for the heat flux.Equivalently, one can combine the energy flux with all the independent mass fluxes of the chemical species in the material [47,70].The latter option is often preferable, since all the needed fluxes are readily available from EMD simulations, while the heat flux also necessitates computing partial enthalpies [71] through postprocessing [72].
The energy flux is the first spatial moment of the time derivative of the energy density [47]: where R ℓ and P ℓ are the position and linear momentum of the ℓth atom, respectively, and f ℓ is force acting on it.There is no a priori correct way to define the energy density of a condensed matter system.However, this ambiguity does not pose an issue, as this quantity is not directly measurable in experiments.Moreover, quantities dependent on it, such as thermal conductivity, remain unaffected by the precise expression of local densities.This concept is referred to as the gauge invariance of transport coefficients [47,73].
The explicit expression of the energy flux depends on the Hamiltonian of the system, which, for a classical system of N particles, takes the general form: where V denotes the potential energy.Given the gaugeinvariance principle, any local breakdown of the energy density into atomic contributions is suitable for computing the energy flux [47,73].A valid choice is thus where each atomic energy, ϵ ℓ , is concentrated on the respective atom and can be expressed as: with the first term representing atomic kinetic energy and the second indicating the portion of potential energy assigned to the ℓth atom, such that V = ℓ V ℓ .These atomic energies are phase-space variables dependent on time through atomic positions and momenta.Consequently, the energy flux in Eq. (4) becomes This equation is well-defined in PBC, as it solely depends on interatomic distances calculated using the minimumimage convention.As such, the formula proves suitable for MD simulations involving bulk systems.The concepts discussed thus far are applicable to both empirical force fields (FFs) and MLPs.In the case of the latter, it is possible to generate MD simulations of ab initio quality to sample the relevant fluxes, including the energy flux of Eq. ( 8), thus leading to the computation of κ.This approach becomes especially significant for glasses, where, as previously mentioned, challenges arise in lattice-dynamical calculations involving MLPs.While for semi-empirical FFs assigning the atomic energies in Eq. ( 7) might seem reasonable and simple, especially in the case of pairwise potentials, it is less so in the case of MLPs.MLPs typically rely on a decomposition of the global target quantity into local, atom-centered contributions.This approach offers several advantages, such as computational scalability and the possibility of retrospective interpretation of patterns of local contributions, provided that local predictions are sufficiently "rigid", i.e., that their value is robust under perturbations of the training set and/or model parameters [74].

II. RESULTS
In this Section we first develop, for demonstration purposes, a dataset for the paradigmatic amorphous solid vitreous silica (v-SiO 2 ) based on a referece semi-empirical FF.We train different MLPs on this dataset, and analyze their ability in reproducing thermal conductivity calculations.We then discuss the current limitations of lattice dynamics in thermal conductivity calculations with MLPs, and how these are absent in EMD simulations.We finally apply GKMD to an amorphous alloy useful in the context of solid-state electrolytes.

A. Toy model of vitreous silica
The reference potential is a Tersoff FF developed in [75] as a short-range potential able to describe the amorphous phase of silica with relatively good accuracy.We choose a short-range reference potential in order to keep the model simple; moreover, the MLP to be trained is rigorously short-range, even if strategies to incorporate bona fide long-range interactions exist and are being developed [76][77][78][79][80][81].The dataset is generated starting from 100 independently quenched glassy samples, equilibrated at 10 different temperatures ranging from 300 K to 7000 K. From each of the 100 simulations, 500 uncorrelated configurations are drawn once every 10 ps.The dataset thus comprises 50,000 glassy and molten configurations.Further details can be found in the Supplementary Material.
We trained a committee of four deep neural network DeePMD potentials (DPs) [82][83][84] smooth edition [85] and a neuroevolution potential (NEP) [86] to reproduce energies, forces and virials of the dataset.All the MLPs have a cutoff of 3 Å, and are trained for more than 400 epochs to ensure the proper minimization of the loss function.The optimization of the loss function for the DPs is performed through the Adam stochastic gradient descent [87], while for the NEP via a genetic algorithm [86].Each model within the DP committee differs from the others based on the initial random seed used in the minimization.The NEP features instead a different architecture and slightly different descriptors with respect to DPs [86].
We assess the performance of these models on a validation set, which consists of an additional 2500 configurations produced using the same methodology as the training set.Fig. 1 presents parity plots comparing the reference and predicted potential energies, forces, and virials for the four DP models in the committee.All the DPs exhibit comparable accuracy.The root mean square errors (RMSEs) and relative error with respect of the standard deviation (in parentheses), are 2.86 ± 0.09 meV/atom for the energy (1%), 237.04 ± 0.14 meV/Å for the force (8%), and 6.96 ± 0.04 meV for the virial (5%).In Fig. 1 the performance of only one model of the DP committee is shown (red) for simplicity.The NEP features RMSEs and relative errors of 0.9 meV/atom for the energy (< 1%), of 311 meV/Å for the force (11%), and of 13 meV for the virial (8%).

B. Lattice dynamics vs molecular dynamics
Subsequently, we test the models' performance in replicating the reference thermal conductivity.We compute the thermal conductivity of vitreous silica as determined by Eq. ( 1).To do so, we sample the energy flux and the atomic mass fluxes through EMD simulations on a system with 648 atoms.Given that the Tersoff FF is a many-body potential, particular care needs to be taken when computing the energy flux [88].A correct implementation of the Tersoff energy flux can be found in gpumd [89].As reference, we carried out a 1 ns-long EMD simulation at 300 K using gpumd, collecting the energy and mass fluxes every 5 fs.Further computational details can be found in the Supplementary Material.We employ the same trajectory to sample these quantities using the committee of DPs and the NEP, utilizing the DeePMD and NEP implementations in lammps [89][90][91].Using the same trajectory to test the models allows us to isolate the effects of gauge invariance, eliminating subtle differences that may arise when employing independent EMD runs to sample different trajectories.The thermal conductivity is finally obtained through cepstral analysis of the fluxes' time series [70,92], as implemented in SporTran [93].The total energy decomposition into local contributions is in general different for different models, since the latter is a task arbitrarily performed by the ML algorithm.Thus, also the fluxes differ among the pool of models.Nevertheless, what must remain consistent is the physical observable, namely, the thermal conductivity, as commanded by gauge invariance [73].
In Fig. 2, we show the distributions (normalized histograms) of the deviation of the atomic energies per species with respect to their average value, sampled along the trajectory, as computed by the available MLPs.The DPs feature similar distributions of such quantity, so only one member of the committee is shown, while the NEP results are quite different.Nonetheless, even within the DP committee, the average values of the atomic energies are significantly different.In particular, the set of mean energies for the Tersoff FF, the DPs, and the NEP is reported in Tab.I.
We stress again that such a difference is expected, as the atomic energies are not target quantities in the learning scheme, and only the total energy is physically observable.The difference in atomic energies predicted by the different models is translated into a difference in en- ergy fluxes, as shown in the top panel of Fig. 3. Crucially, the thermal conductivity, as represented by the zero-frequency value of the fluxes' power spectral density [70,92], remains unchanged, as reported in the bottom panel of Fig. 3.We now assess the models' performance in reproducing quantities relevant for lattice-dynamical calculations.We compute the second and third-order interatomic force constants with lammps and obtain normal-mode frequencies and linewidths using κALDo [63].We employ these quantities to compute the QHGK thermal conductivity with the Tersoff FF, a DP, and the NEP.Fig. 4 showcases the results, including the vibrational density of states (VDOS), normal-mode linewidths, γ, and QHGK thermal conductivity.While the VDOS outcomes in the upper panel exhibit agreement, the same cannot be said for anharmonic normal-mode linewidths (central panel), which appear to be accurately captured only by the NEP.In contrast, the DP results considerably overestimate the linewidths, even surpassing the quasi-harmonic regime, denoted by the dashed black line.These disparities have direct implications for the computed QHGK thermal conductivities, as evident in the lower panel of Fig. 4. Specifically, while the NEP and Tersoff results align with each other, the DP results do not share this agreement.
It is essential to recognize that apparently satisfactory performance of the NEP leads to potentially misleading outcomes.Upon closer examination of the central panel of Fig. 4, a notable issue emerges: at low frequencies (in this instance, the first available finite frequency) there are disproportionately large linewidths.This stark contrast with the expected decrease in linewidths as frequencies decrease, as dictated by the hydrodynamics of solids [64,94], precludes the hydrodynamic extrapola- tion the NEP results to obtain a size-converged value for bulk thermal conductivity [60,64].As a further proof that MLPs still present some issues with lattice dynamics, we compute the anharmonic linewidths of a more complex NEP for vitreous silicon (v-Si) developed in [58] and based on the dataset of [95].As shown in Fig. 5, increasing the size of the sample exacerbates the issue of the low-frequency peak in the linewidths, which hinders the possibility of carrying out size-converged thermal conductivity calculations [60,64].This issue persists even when manually imposing acoustic sum rules due to global translational invariance and symmetry under Cartesiandirection permutation.
Due to these limitations, lattice-dynamical calculations prove to be of limited utility when employing MLPs.Consequently, the alternative approach is to resort to GKMD simulations, albeit at the cost of neglecting the quantum BE occupation.It is worth noting that in specific instances, like vitreous silica, this omission remains notable even at room temperature [96,97].Conversely, for other simple glasses like vitreous silicon, this concern does not apply [58,64].

C. Thermal conductivity of Li-intercalated silicon
To exemplify the GKMD methodology within the context of MLPs, we conduct a thermal conductivity assessment of a PBE-accurate model describing amorphous LiSi, as developed by [98].This MLP was employed to investigate the properties of silicon, encompassing both crystalline and amorphous phases, as an anode material for Li-based electrolytes [98].Thermal transport plays a pivotal role in the design of electrolytes, especially for solid-state systems [99].Indeed, an exceedingly low thermal conductivity can lead to excessive heat generation, especially during rapid charging processes, thereby posing the risk of critical incidents such as material melting or explosions.Moreover, the management of thermal dissipation is crucial to optimizing energy conservation and utilization, requiring a delicate equilibrium between minimizing heat dissipation and maximizing electric flux throughout the charging cycle [99].In addition, in situations involving materials where ionic diffusion is not only expected but also a desired attribute, such as solid-state electrolytes, one is forced to use methods that do not rely on the presence of well-defined atomic equilibrium positions [100].
We conduct EMD simulations on Li x Si 1 − x with varying concentration of intercalated lithium, denoted as x, under room temperature conditions (T = 300 K, p = 0 bar).Amorphous samples are prepared through a meltquench-anneal procedure as outlined in [98].The samples so obtained are used as initial configurations for EMD simulations that sample the relevant fluxes in canonical runs [101] carried out for 3 ns at 300 K.The simulations are performed with gpumd using a NEP trained on the dataset of [98].Further details can be found in the Supplementary Material.
Results for the thermal conductivity as a function of Li concentration are presented in Fig. 6.Size effects are important for low Li concentrations, where the system is close to pure vitreous silicon.v-Si is known for being severely affected by size effects [64] due to its high local order [95].Thus, samples with 3,000 atoms are not converged in size, while 70,000 atoms appear to be enough, when compared to calculations on systems with 150,000 atoms.
At higher Li concentrations, even 3,000 atoms are sufficient to achieve convergence.The behavior of κ(x) features a minimum x ≈ 0.5, followed by an increase for larger values of x.Pure v-Si features κ = 1.24 ± 0.03 W m −1 K −1 .This value can be compared with calculations [58] done on structure of similar size that employ a NEP trained on a dataset relying on PW91 DFT calculations [95], rather than PBE.The results are compatible at high quenching rate, while we report a value which is ≈ 30% lower than the one of [58] at the lowest quenching rate.Experimental data on the thermal conductivity of v-Si severely depend on the experimental sample size, due to the high relevance of propagons in this material.At 300 K, they can range from less than 1 W m −1 K −1 [102] to around 4 W m −1 K −1 [103,104].We are not aware of existing experimental data on the thermal conductivity of Li-intercalated amorphous silicon.
The qualitative behavior of the thermal conductivity as a function of Li concentration aligns with calculations on crystalline alloys [105], where the distinctive U-shape of κ(x) is interpreted in terms of phonon scattering due to isotopic mass disorder [106] at the perturbative level.Conversely, in glasses the thermal conductivity reduction with increasing concentration of different-species atoms is mostly due to the increased localization of low-and midfrequency vibrational modes, together with the broadening of the respective linewidths, which hinders their ability to transport heat [107].

III. CONCLUSIONS
In this work, we have reviewed the theory of thermal transport in amorphous solids, focusing on the role of MLPs as a tool to expedite and, in some cases, allow for thermal-transport characterization of glasses from atomistic simulations.We have build an example of MLP for vitreous silica able to reproduce the potential energy surface of a Tersoff FF, and used it to highlight some challenges still present when dealing with lattice dynamics, leaving Green-Kubo molecular dynamics as a valuable alternative to carry out calculations with MLPs.We then applied such methodologies to the Liconcentration dependence of the thermal conductivity of lithium-intercalated amorphous silicon, finding that κ features a minimum at half concentration, coherent with analogous calculations made on other amorphous alloys and interpreted in terms of the localization of propagating modes due to chemical disorder.tum Computing (grant number CN00000013).FG acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Action IF-EF-ST, grant agreement number 101018557 (TRANQUIL).

FIG. 1 .
FIG. 1. Parity plots of the committee of MLPs trained on the vitreous silica Tersoff dataset.

FIG. 3 .
FIG. 3. Thermal conductivity of the MLPs committee.Simulations are carried out at 300 K and ambient pressure.Top panel: the reference (Tersoff) energy flux and the DP committee ones are different, in general.Bottom panel: cepstral analysis of the fluxes' time series.The thermal conductivity, given by the zero-frequency value of the curve, is the same for different models, as commanded by gauge invariance.

FIG. 4 .
FIG. 4. Results based on lattice dynamics of vitreous silica computed with the Tersoff FF, the DP, and the NEP.Upper panel: vibrational density of states computed with the three models.Central panel: normal-mode linewidths.Lower panel: QHGK thermal conductivity.EQ stands for the classical equipartition, while BE for the Bose-Einstein occupation.

FIG. 5 .
FIG. 5. Anharmonic linewidths of a NEP model for v-Si at 300 K, for samples of different sizes.

FIG. 6 .
FIG.6.Thermal conductivity of v-LixSi1 − x as a function Li concentration at room-temperature conditions for different system sizes.