Thermal Conductivity of Glass-Forming Liquids

TL;DR

This study investigates how thermal conductivity in glass-forming liquids and glasses varies with temperature, aging, and preparation methods, revealing links to structural relaxation and vibrational modes through extensive simulations.

Contribution

It provides new insights into the dependence of thermal conductivity on aging, cooling rate, and inherent structure properties in glass-forming systems.

Findings

Thermal conductivity decreases with aging near the glass transition.

Lower cooling rates produce glasses with reduced thermal conductivity.

Reduced conductivity correlates with less extended high-frequency vibrational modes.

Abstract

Thermal conductivity of a model glass-forming system in the liquid and glass states is studied using extensive numerical simulations. We show that near the glass transition temperture, where the structural relaxation time becomes very long, the measured thermal conductivity decreases with increasing age. Secondly the thermal conductivity of the disordered solid obtained at low temperatures depends on the cooling rate with which it was prepared, with lower cooling rates leading to lower thermal conductivity. Our analysis links this decrease of the thermal conductivity with increased exploration of lower-energy inherent structures of the underlying potential energy landscape. Further we show that the lowering of conductivity for lower-energy inherent structures is related to the high frequency harmonic modes associated with the inherent structure being less extended.

Equations

$${\mathcal{J}}^{\alpha }(\vec{r},t)=\sum _l\delta (\vec{r}-{\vec{r}}_l)[{ϵ}_l{v}_l^{\alpha }+\frac{1}{2}\sum _{n\ne l}(r_l^{\alpha }-r_n^{\alpha })j_{l,n}]$$

$$\kappa =\frac{⟨{\mathcal{J}}^x⟩\times L_x^{\prime }}{\Delta T}.$$

$$\kappa =\frac{1}{3k_B{T}^2}\underset{\tau \to \infty }{\lim }\underset{V\to \infty }{\lim }\frac{\rho }{N}{\int }_0^{\tau }dt\text{mspace}1.0mu⟨{\vec{\mathcal{J}}}_{tot}(0)\cdot {\vec{\mathcal{J}}}_{tot}(t)⟩$$

Full text

Introduction – Thermal conductivity is an important material property with wide ranging applications and understanding the mechanism of heat transport in different materials is crucial from both theoretical and applied perspectives. It is well known that crystalline materials transport heat more efficiently than glassy and disordered materials, for which the thermal conductivity ($\kappa$) can be orders of magnitude smaller Cahill:88 ; Pohl:71 .

A well-studied feature of glassy systems is the temperature $(T)$-dependence of their thermal conductivity, which is very distinct from crystals. One typically finds $\kappa\sim T^{2}$ at very low temperatures, followed by a pronounced plateau and then an eventual gradual increase with temperature. The low temperature features are expected to be of quantum-menchanical origin and a range of mechanisms have been proposed to explain them Anderson:72 ; Buchenau:92 ; Clare:87 . One scenario relates the plateaue to the presence of excess modes, the so-called boson peak, in the low-frequency regime of the vibrational spectrum of glasses Lubchenko:03 ; schirmacher:2006 ; Nagel:10 , which has thereafter been connected with the presence of elastic heterogeneities in such materials Maruzzo:2013 ; Barrat:16 .

More recently, a number of studies have looked at higher temperatures where classical physics dominates. This includes regimes much below the glass transition temperature where the system gets stuck in the basin of a low-energy inherent structure (IS) (local minimum of the potential energy), as well as intermediate temperatures where the system evolves slowly while undergoing transitions between different basins. In the former case the system can be effectively described by a disordered harmonic model obtained by expanding the many-body potential about the minimum. Heat transport is then related to the diffusivity associated with the normal modes of the harmonic solid Allen:89 ; Allen:93 . In the context of the jamming transition, the change in the form of the diffusivity in a soft-sphere packing was studied in Nagel:09 ; Nagel:10 ; wyart:10 , and a crossover was observed near the boson-peak frequency as the transition was approached. In Barrat:13 , it was shown that the introduction of polydispersity leads to an amorphisation transition, whereby the thermal conductivity of the system decreases considerably.

So far the question as to whether thermal transport in glasses depends upon the age of the material, or on the cooling rate used in its preparation has not been explored. Here we address this question. It is well-known that the material properties of glasses depend on their age and history of preparation, e.g. the cooling rate by which they were quenched from a fluid phase Kob:97 ; Binder:96 ; glassbook . This dependence on the history can be seen, for example, in the transient response of glasses to applied shear Shi:06 ; Moorcroft:11 ; Chaudhuri:16 . Surprisingly, in the case of thermal transport, there is not much systematic study. Using extensive numerical simulations of a model structural glass former, we show that the thermal conductivity does depend on the age of the glass and on the rate of cooling during preparation. Further, we demonstrate that the decrease of thermal conductivity in glasses with growing age or slower cooling rate is linked to the exploration of ISs with lower energy. By calculating the participation ratio of the normal modes associated with different ISs, we see that the decrease in $\kappa$ is related to an increase in the degree of localization of high frequency vibrational modes.

Another puzzling issue that we address relates to the question of whether the thermal conductivity of the disordered harmonic system that describes the dynamics of small displacements from an IS is infinite. Rayleigh scattering of low frequency ($\omega$) phonons is expected to lead to mean free paths $\ell\sim 1/\omega^{4}$ ziman , which can cause a divergence of $\kappa$ and this is seen in some recent studies Nagel:09 ; kundu:10 ; chaudhuri:10 . In our study we compute the thermal conductivity $\kappa$ of a glass former using three numerical methods (see below), all of which indicate a finite conductivity as $T\to 0$, where the harmonic approximation should be valid. However, we can not rule out a divergence because the system sizes used in our simulations are modest (number of particles $\leq 20,000$). A recent study lerner2016 has shown that the low-frequency part of the vibrational spectrum for such values of $N$ is dominated by the boson peak modes and the effects of long-wavelength phonon modes are observed only in systems with much larger values of $N$. Since the divergence of $\kappa$ is expected to arise from these phonon modes, simulations of much larger systems would be required to provide a clear answer to the question of whether $\kappa$ diverges as $T\to 0$. Our work shows that the contribution of the boson-peak modes to the thermal conductivity remains finite as $T\to 0$.

We consider the well-known model glass former - the Kob-Andersen binary Lennard-Jones mixture Kob:94 (details in section I of Supplementary Information (SI)). We consider $N$ particles in a box of dimensions $L_{x}\times L_{y}\times L_{z}$. Our study is done for a number density of $\rho=N/(L_{x}L_{y}L_{z})=1.2$, for which the Vogel-Fulcher-Tamman ($VFT$) transition temperature of this glass former is $T_{0}\approx 0.3$ in reduced Lennard-Jones units. We compute the thermal conductivity at different temperatures above and below the glass transition. The supercooled liquid and glassy states we have studied were obtained via a well-defined cooling protocol Sastry:98 . The system is well-equilibrated at a high temperature $(T=2.50)$ in its liquid state and then cooled to a low temperature through a number of intermediate temperature steps. The number of intermediate steps and the amount of time spent in each step determine the cooling rate Sastry:98 .

Thermal conductivity – Denoting the position and the velocity of the $l$-th particle by ${\vec{r}}_{l}=\left\{r^{\alpha}_{l}\right\}$ and ${\vec{v}}_{l}=\left\{v^{\alpha}_{l}\right\}$ respectively (where $\alpha=x,y,z$), the energy current density at point ${\vec{r}}$ is given by Hansen:06 :

[TABLE]

where $\epsilon_{l}={\mathbf{v}_{l}^{2}}/2+(1/2)\sum_{n\neq l}V(r_{ln})$ is the energy of the $l$-th particle, and $j_{l,n}=\frac{1}{2}\sum_{\nu}(v_{l}^{\nu}+v_{n}^{\nu})f_{ln}^{\nu}$, with $f_{ln}^{\alpha}=-\partial{V(r_{ln})}/\partial{r_{l}^{\alpha}}$. We evaluated the thermal conductivity by several methods, which we now describe.

Nonequilibrium method (NEMD) – We consider a system with mean temperature $T_{m}$, impose a small temperature difference $\Delta T$ in the $x$-direction across length $L_{x}^{\prime}<L_{x}$, and measure the average steady-state heat current density $\langle\mathcal{J}^{x}\rangle=\frac{1}{V}\int_{V}d{\vec{r}}\langle\mathcal{J}^{x}({\vec{r}},t)\rangle$ where $V$ is the volume of the region in the sample where the measurement is made. The thermal conductivity is then given by

[TABLE]

We consider a rectangular geometry, as illustrated in the schematic diagram of Fig. 1. Periodic boundary conditions were maintained in all directions. We first prepare the initial equilibrium states (at mean temperature $T_{m}$) via the cooling protocol discussed above. A temperature gradient was then applied by keeping a heat source at temperature $1.05T_{m}$ in a region of width $L_{x}/10$ in the middle of the system, and two heat sinks of width $L_{x}/20$ at the two ends at $0.95T_{m}$. The thermostating of the hot and cold reservoirs was done by drawing the velocities of the particles in each reservoir from a Maxwell-Boltzmann distribution at the appropriate temperatures, at intervals of time 0.1. After an initial transient period ($\delta{t}\sim 10^{4}$) following the imposition of the thermal gradient, the system reaches a non-equilibrium steady state, wherein the spatial profiles of local temperature and heat current density typically look like those shown in Fig. 1. The spatial variation of the local temperature $T_{i}$ and the heat current density $\langle\widetilde{\mathcal{J}}^{x}_{i}\rangle$ was obtained by binning $L_{x}$ into $40$ segments, labelled $i=1,\ldots 40$, and doing the averages separately for each segment with $\delta{t}\sim 5\times{10^{4}}$, after the transient time (details in section II(a) of SI). Using the steady state current value from simulations Eq. (2) with $L^{\prime}_{x}=0.4\times L_{x}$ gives $\kappa$ at temperature $T_{m}$. We have done averages over $32$ independent MD trajectories.

Green-Kubo method(EMD) – Alternatively, $\kappa$ can be calculated using the Green-Kubo relation Hansen:06 :

[TABLE]

where ${\vec{\mathcal{J}}}_{tot}(t)$ is the integral of the heat current density ${\vec{\mathcal{J}}}({\vec{r}},t)$ over the whole system and the time correlation function is evaluated at equilibrium.

For the Green-Kubo calculation we considered a cubic simulation box, with the same density of particles as in the nonequilibrium case. An initial equilibrium state is prepared following the cooling protocol described earlier, and subsequently the NVE dynamics is switched on. The fluctuations of the heat current are then measured over a time-window of ${\delta}t~{}\sim 10^{4}$ as the system evolves and the thermal conductivity ($\kappa$) is calculated using the Green-Kubo formula (Eq. 3). The estimates of $\kappa$ were obtained via averages over $96$ independent trajectories.

In Fig. 2 we compare the values of $\kappa$ measured via the two schemes outlined above, for a given cooling rate. As expected for a disordered system, the magnitude of the thermal conductivity continuously decreases over the entire temperature range, before nearly saturating at very low temperatures. We observe very good agreement between the two different measurements of $\kappa$ over a very broad range of temperatures. Surprisingly, this is true even in the glassy regime, where the two methods provide the same result for states which have the same preparation history, i.e. produced via the same cooling rate. In the rest of the work we have used these two methods independently, in accordance with the question to be addressed. By performing simulations at various system sizes (keeping $\rho$ constant), we have verified that the results for $\kappa$ become essentially independent of system size for $N\sim 10^{3}$ (EMD) and $N\sim 10^{4}$ (NEMD). In the inset of Fig. 2 we show the system size dependence of $\kappa$ (details in section II of SI).

Effect of cooling rate – It is known that the cooling rate influences the glassy state obtained at low temperatures glassbook . For example, the mean energy of the underlying ISs, $\overline{E}_{IS}$, corresponding to any temperature, depends on the choice of the cooling rate Sastry:98 , with slower cooling leading to exploration of lower energy ISs. Using data from our simulations, this is illustrated in the inset of Fig. 3, where we see that below a certain temperature, the variation of $\overline{E}_{IS}$ with $T$ branches out for two contrasting cooling rates. In the same spirit, we check how $\kappa$ (from non-equilibrium measurements) varies with temperatures along the same two cooling branches. This is shown in Fig. 3, and we observe that the slower cooling corresponds to lower values of $\kappa$. However, compared to data shown for $\overline{E}_{IS}$, in this case, the branching occurs at a lower temperature, typically where the dynamics is known to fall out of equilibrium glassbook .

Effect of aging – The dependence of thermal conductivity on the history of preparation of the glassy state, as demonstrated above, implies that $\kappa$ would also vary with the age of the glass when the system evolves after a thermal quench from high temperature to low temperatures. To investigate this, we quench the system from $T=2.50$ to $T=0.30$ (which is close to the $VFT$ transition temperature) and let the system evolve in a constant temperature environment for a time $t_{w}$. Once the system has reached the age $t_{w}$, the $NVE$ dynamics is switched on. Under such conditions, we measure the total heat current and calculate the corresponding auto-correlation function $\langle{\vec{\mathcal{J}}}_{tot}(t_{w})\cdotp\vec{{\mathcal{J}}}_{tot}(t_{w}+t)\rangle$. In Fig. 4, we show the current auto-correlation function for increasing age $t_{w}$. Note that the correlators relax quickly in time, even though the structural relaxations are extremely slow at this temperature Scheidler:01 . The variation of conductivity $\kappa$ with changing age $t_{w}$, calculated from these correlation functions using the Green-Kubo formula, is shown in the inset of Fig. 4: with increasing age, the thermal conductivity decreases. Now, it is known that glassy systems explore deeper minima in the potential energy landscape with increasing age Sciortino:99 . Therefore, the observed dependence of $\kappa$ on $t_{w}$ and our earlier observation regarding the dependence of $\kappa$ on cooling rate suggest that the thermal conductivity of the glass is closely related to the energy of the relevant ISs.

Dependence on $E_{IS}$ – We now explore this possible link between properties of ISs and the values of $\kappa$. We consider a number of ISs with mean energy $\overline{E}_{IS}$ and standard deviation $\sim 10^{-6}$ and prepare low temperature ($T=0.002$) initial states consistent with small harmonic fluctuations around the potential minimum (details in section IV of SI). In the harmonic theory, the total energy of an equilibrium state with temperature $T$ is $\overline{E}_{IS}+[(6N-3)k_{B}T/2]$. An initial state with this energy was prepared by choosing positions corresponding to the IS and giving each particle a random velocity taken from a Maxwell distribution at temperature $2T$. MD simulations within the NVE ensemble are carried out for $N=1000$. For the states that remain confined to the basin of the initial IS during the course of the simulation, the thermal conductivity is calculated using the Green-Kubo method. In Fig. 5 we show the dependence of $\kappa$ on the energy of the IS – we find $\kappa$ increases linearly with $\overline{E}_{IS}$. This explains why the thermal conductivity decreases with slower cooling, or with increasing age of the system after the thermal quench – the system starts exploring ISs with lower energy.

Properties of harmonic excitations – At very low temperatures, the system mostly remains close to an IS of the potential energy landscape and one can approximate its potential energy to that of a harmonic solid. It is of interest to relate the thermal conductivity to the properties of the harmonic excitations, especially the effect of Anderson localization nagel84 ; laird91 ; kundu:10 . An expansion of the potential energy to quadratic order in the displacements from the IS gives the Hessian matrix which completely describes the properties of the harmonic solid. We ask how the properties of the eigenvalues and eigenfunctions of the Hessian matrix affect thermal transport and the relation between $\kappa$ and $\overline{E}_{IS}$. Following the formulation in Allen:93 , the conductivity is given by $\kappa=(k_{B}/V)\int_{0}^{\infty}d\omega g(\omega)d(\omega)$, where the density of states $g(\omega)$ and the thermal diffusivity $d(\omega)$ are known completely in terms of the eigenvalues and eigenfunctions of the Hessian matrix (details in section III of SI). We expect the diffusivity of modes to depend on their degree of delocalization which can be quantified by the participation ratio ($PR(\omega)$).

In Fig. 6(a), we consider ISs with three different $\overline{E}_{IS}$ and compare the results for $d(\omega),g(\omega)$ and $PR(\omega)$ for the three cases. We see that, $d(\omega)$ of the ISs with lower $\overline{E}_{IS}$ are smaller over most of the allowed frequency range except for a small region near the boson peak, as compared to the ISs with higher $\overline{E}_{IS}$. We recall that for amorphous solids, excess density of states compared to the Debye law, $g(\omega)\sim\omega^{2}$, appears in the low frequency region, and this is known as the boson peak glassbook . Here we find that the boson peak is more pronounced for ISs with higher $\overline{E}_{IS}$ than the ISs with lower $\overline{E}_{IS}$, as shown in Fig. 6(b). Finally from the PR (shown in Fig. 6(c)) we see that both the low and high frequency vibrational modes are more extended for the ISs with higher $\overline{E}_{IS}$. This indicates that the higher thermal conductivity is related to the increased delocalization of these vibrational modes, with the high frequency modes playing a more important role (details in section III(e) of SI).

Conclusions — Using numerical simulations, we have studied thermal transport in a model glass-forming liquid. The values of the thermal conductivity, obtained from non-equilibrium and Green-Kubo calculations are found to agree well over a wide range of temperatures. One of our main findings is that the phenomenon of aging has a significant effect on thermal transport. Near the glass transition temperature, the conductivity drops by about $50\%$ with increasing age. Similarly the conductivity of the low temperature glass can vary by about $10\%$ depending on the coooling protocol used to form it. The lowering of conductivity with growing age or slower cooling can be rationalized in terms of the system exploring lower energy minima and our numerical results demonstrating that the thermal conductivity increases linearly with the energy of the IS. Further analysis of the harmonic solid associated with the ISs helped us to understand the underlying reasons in terms of the diffusivity and the extent of localization of the normal modes. For the disordered solid, we do not find a diverging diffusivity at low frequencies, but we can not rule out the possibility of a divergence in much larger systems for which long-wavelength phonon modes are expected to dominate the low-frequency part of the vibrational spectrum.

We thank J.-L. Barrat and J. Horbach for discussions; D. Banerjee and H. Khare for their help with figures. PJB & RM acknowledges financial support from CSIR, India. We thank SERC IISc, IMSc and ICTS for computing support. AD thanks support from the Indo-Israel joint research project No. 6-8/2014(IC) and from the French Ministry of Education through the grant ANR (EDNHS). CD acknowledges support from DST, India.