Theory for Glassy Behavior of Supercooled Liquid Mixtures

TL;DR

This paper develops a dynamical facilitation-based model to predict the relaxation behavior and glass transition temperatures of supercooled liquid mixtures, moving beyond empirical interpolation methods.

Contribution

It introduces a multi-component kinetically constrained model that captures mixture dynamics and predicts relaxation times using a single effective energy scale.

Findings

Experimental relaxation times collapse onto the parabolic law.

Model accurately predicts glass transition temperatures of mixtures.

Provides a physical basis for mixture behavior beyond empirical models.

Abstract

We present a model for glassy dynamics in supercooled liquid mixtures. Given the relaxation behavior of individual supercooled liquids, the model predicts the relaxation times of their mixtures as temperature is decreased. The model is based on dynamical facilitation theory for glassy dynamics, which provides a physical basis for relaxation and vitrification of a supercooled liquid. This is in contrast to empirical linear interpolations such as the Gordon-Taylor equation typically used to predict glass transition temperatures of liquid mixtures. To understand the behavior of supercooled liquid mixtures we consider a multi-component variant of the kinetically constrained East model in which components have a different energy scale and can also diffuse when locally mobile regions, i.e., excitations, are present. Using a variational approach we determine an effective single component model with a single effective energy scale that best approximates a mixture. When scaled by this single effective energy, we show that experimental relaxation times of many liquid mixtures all collapse onto the 'parabolic law' predicted by dynamical facilitation theory. The model can be used to predict transport properties and glass transition temperatures of mixtures of glassy materials, with implications in atmospheric chemistry, biology, and pharmaceuticals.

Figures

Equations

$$\ln \tau /{\tau }^{\mathrm{o}}\sim J^2(\beta -{\beta }^{\mathrm{o}}{)}^2.$$

$$H(\{n_i,p_i\},x)=\sum _i^N{n}_i(p_i{J}_1+(1-p_i)J_0)+C(\{p_i\}),$$

$$C(\{p_i\})=\{\begin{array}{ll}0,& \text{if }\frac{1}{N}{\sum }_i^N{p}_i=x\\ \infty ,& \text{otherwise.}\end{array}$$

$$\begin{array}{rl}{k}_{0\to 1,i}& \propto (p_i{e}^{-\beta J_1}+(1-p_i)e^{-\beta J_0})n_{i-1}\\ k_{1\to 0,i}& \propto n_{i-1},\end{array}$$

$$r_{i,i\pm 1}\sim n_i{n}_{i\pm 1}((1-p_i)p_{i\pm 1}+(1-p_{i\pm 1})p_i),$$

$$H_{\mathrm{m}}(\{n_i,p_i\},x)=J_{\mathrm{m}}\sum _i^N{n}_i+C(\{p_i\}).$$

$$\begin{array}{rl}Z& =\frac{Z_{\mathrm{m}}}{Z_{\mathrm{m}}}\sum _{n_1}\dots \sum _{n_N}\sum _{p_1}\dots \sum _{p_N}{e}^{-\beta H}{e}^{-\beta H_{\mathrm{m}}}{e}^{+\beta H_{\mathrm{m}}}\\ & =Z_{\mathrm{m}}⟨\exp (-\beta \Delta H){⟩}_{H_{\mathrm{m}}}\ge Z_{\mathrm{m}}\exp (-\beta ⟨\Delta H{⟩}_{H_{\mathrm{m}}}),\end{array}$$

$$J_{\mathrm{m}}=x{J}_1+(1-x)J_0.$$

$$H^{\mathrm{o}}(\{n_{i},p_{i}\},x)=\sum_{i}^{N}n_{i}\Big{(}\frac{\widetilde{\beta}_{1}}{\beta}p_{i}J_{1}+\frac{\widetilde{\beta}_{0}}{\beta}(1-p_{i})J_{0}\Big{)}+C(\{p_{i}\}).$$

$$\begin{split}Z^{\mathrm{o}}&=\sum_{n_{1}}\ldots\sum_{n_{N}}\sum_{p_{1}}\ldots\sum_{p_{N}}\exp\big{\{}-\sum_{i}^{N}n_{i}\{\widetilde{\beta}_{1}J_{1}p_{i}\\ &\hskip 71.13188pt+\widetilde{\beta}_{0}J_{0}(1-p_{i})\}-\beta C(\{p_{i}\})\big{\}},\end{split}$$

$$H_{\mathrm{m}}^{\mathrm{o}}(\{n_i,p_i\},x)=\frac{{\tilde{\beta }}_{\mathrm{m}}{J}_{\mathrm{m}}}{\beta }\sum _i^N{n}_i+C\sum _i^N{p}_i,$$

$$\begin{split}Z_{\mathrm{m}}^{\mathrm{o}}&=\sum_{n_{1}}\ldots\sum_{n_{N}}\sum_{p_{1}}\ldots\sum_{p_{N}}\exp\big{(}-\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}}\sum_{i}^{N}n_{i}\\ &\hskip 128.0374pt-\beta C\sum_{i}^{N}p_{i}\big{)},\end{split}$$

$${\tilde{\beta }}_{\mathrm{m}}{J}_{\mathrm{m}}=x{\tilde{\beta }}_1{J}_1+(1-x){\tilde{\beta }}_0{J}_0.$$

Full text

This paper presents a model for glassy dynamics in mixtures of liquids. Understanding and predicting glassy behavior of liquid mixtures and solutions has been important in investigating the fundamental properties of glassy dynamics and the glass transition Blochowicz and Rössler (2004); Wang et al. (2010); Foffi et al. (2003), as well as in areas of research where solutions of liquids undergo vitrification into amorphous solids, and where crystallization needs to be avoided. These areas include chemistry of solutions Angell et al. (1978), and pharmaceutical and food industries Hancock and Zografi (1997); Yu (2001); Bhandari and Howes (1999). Insight into properties of glassy mixtures is also important in more complex systems such as atmospheric organic aerosols that show glassy behaviour Virtanen et al. (2010); Vaden et al. (2010); Cappa and Wilson (2011); Koop et al. (2011), and in the preservation of biological cells where vitrification may assist the recovery of cells after freezing and/or desiccation Crowe et al. (1998); Green and Angell (1989). The typical approach to predict the glass transition temperature of liquid mixtures is an empirical linear interpolation between the glass transition temperatures of individual components, such as the Gordon–Taylor equation Gordon and Taylor (1952), without a physical basis in a microscopic theory for glassy dynamics. In this work we use the perspective of dynamical facilitation theory to develop a model for glassy dynamics of binary mixtures, which can be used to predict relaxation behavior and glass transition temperatures of supercooled liquid mixtures.

When a glass-former is cooled below its onset temperature ($T^{\mathrm{o}}$) its relaxation time exhibits a super-Arrhenius increase with decreasing temperature Ediger et al. (1996); Angell et al. (2000), and its microstructure exhibits dynamical heterogeneity, i.e., transient but distinct mesoscopic regions of particle mobility and immobility. The increase in relaxation time is explained by dynamical facilitation theory, the idea that regions of particle mobility (‘excitations’) facilitate the motion of neighboring regions in a hierarchical manner, with small motions leading to larger motions Palmer et al. (1984); Sollich and Evans (1999, 2003); Garrahan and Chandler (2002); Chandler and Garrahan (2010); Keys et al. (2011); Speck (2019). The theory predicts a super-Arrhenius form, quadratic in inverse temperature, for the increase in relaxation time $\tau$ with decreasing temperature. This quadratic form collapses experimental data for relaxation times of a wide array of glass-forming liquids onto one universal curve, the parabolic law, given by Speck (2019); Elmatad et al. (2009)

[TABLE]

Here $\beta=1/T$ is inverse temperature, $\beta^{\mathrm{o}}=1/T^{\mathrm{o}}$ is the inverse onset temperature, $J$ is a property of the liquid and is related to the free energy of creating an excitation, and $\tau^{\mathrm{o}}$ is the relaxation time at $T^{\mathrm{o}}$. With knowledge of the onset temperature $T^{\mathrm{o}}$, the energy scale $J$, and the reference time scale $\tau^{\mathrm{o}}$, one can predict the relaxation time at a given temperature, and the glass transition temperature for a given cooling rate Keys et al. (2013); Hudson and Mandadapu (2018). The parabolic law (Eq. 1) differs from the commonly used empirical Vogel–Fulcher–Tammann relation Ediger et al. (1996), $\ln\tau\sim\mathrm{const}/(T-T_{\mathrm{K}})$, and other proposed theories Adam and Gibbs (1965); Gotze and Sjogren (1992); Mézard and Parisi (1999); Xia and Wolynes (2001), in that it does not show a singularity at non-zero temperatures. We note that there is no microscopic physical theory to understand the onset temperature, which is the temperature at which the relaxation behavior crosses over from an Arrhenius to a super-Arrhenius form. The idea of dynamical facilitation is based on prototypical models with kinetic constraints Ritort and Sollich (2003); Garrahan et al. (2011); Garrahan (2018), such as the East model Jäckle and Eisinger (1991), that exhibit hierarchical glassy relaxation. The East model is a one-dimensional lattice of $N$ spins with variables $n_{i}=0,1$, where $i=1,2\ldots N$. Each lattice site represents a region of particles and the spin variable, $n_{i}$, indicates the mobility of the region. $n_{i}=1$ represents an excitation, i.e., a region of mobile particles, and $n_{i}=0$ represents a region of immobile particles. The model has a noninteracting energy function $H_{\mathrm{East}}=J\sum_{i}n_{i}$, where $J$ represents the energy required to generate particle mobility, i.e., an excitation, in a particular region. The dynamics proceeds via single spin flips with the constraint that a site can flip only if its left neighbor is excited, i.e., $n_{i-1}=1$, based on the idea that particle mobility is facilitated in the vicinity of other mobile particles Garrahan and Chandler (2002); Chandler and Garrahan (2010); Keys et al. (2011). The typical concentration of excitations is $c=\langle n\rangle=\exp(-\beta J)/(1+\exp(-\beta J))\sim\exp(-\beta J)$ at low temperatures.

Relaxation of the system is defined as all of the spins changing their state at least once, and progresses in the direction of facilitation, i.e., from left to right. Sollich and Evans Sollich and Evans (1999, 2003) rationalized the timescale of relaxation by considering domains between excitations, e.g. 1000001. The height of the energy barrier to relaxation is defined as the maximum number of excitations required in a single configuration during the process of relaxing, or flipping, the rightmost excitation. The minimum height of the barrier for a domain of length $l=2^{k}$ is reasoned to be $k$ excitations by iteratively bisecting the domain. For example, a domain of length $4L$ requires twice as many excitations as a domain of length $2L$, which in turn requires twice as many excitations as a domain of length $L$. The minimum energy barrier is therefore $kJ=(\ln l/\ln 2)J.$ For the average equilibrium domain length $l^{\mathrm{eq}}=1/c\sim\exp(\beta J)$, the energy barrier is $\beta J^{2}/\ln 2$. The relaxation time $\tau$ is inversely related to the Boltzmann probability of this energy barrier, giving $\tau\sim\exp(\gamma\beta^{2}J^{2})$, which is the parabolic law. The factor $\gamma$ represents the number of paths available to relax a domain, which are not accounted for in the Sollich–Evans argument, and more rigorous analysis shows that it is bounded as $1/2\ln 2<\gamma<1/\ln 2$ Chleboun et al. (2013).

A Multi-Component East Model. In the presence of more than one component, the relaxation dynamics includes not only dynamical heterogeneity but also component diffusion and mixing. The interplay between heterogeneous dynamics and mixing in supercooled liquid mixtures is not understood. With the perspective of dynamical facilitation theory and the East model we construct a multicomponent lattice model for mixtures of glass-formers that are well-mixed and do not phase separate. We begin by considering a mixture of two components with different energy scales, $J_{0}$ and $J_{1}$. Each lattice site represents a region of molecules of one component or the other, represented by the lattice variables $p_{i}=0,1$. The energy cost for an excitation depends on the type of component, i.e., the value of $p_{i}$ for that site. The Hamiltonian for this two component system is

[TABLE]

where $J_{0}$ and $J_{1}$ are the energies required to excite a spin with component $p_{i}=0$ or 1 respectively. Here $C(\{p_{i}\})$ is a constraint function that ensures that the total number of each component remains fixed:

[TABLE]

where $x$ is the fraction of sites with $p_{i}=1$.

Extending the idea of dynamical facilitation theory that regions of mobile particles facilitate the movement of neighboring immobile regions, irrespective of the type of particles in each region, we define the dynamics of the mixture model as similar to the single component East model, i.e., a spin can flip only if its left neighbor is excited. The rates of creating and destroying an excitation at site $i$ are

[TABLE]

where $n_{i}=0,1$ again represents regions of immobile and mobile particles respectively. With the reasoning that only regions of mobile particles are able to diffuse, two adjacent regions will be able to exchange components only if they are mobile. We model this as two adjacent spins exchanging $p_{i}$ variables only if they are both excited, as illustrated in Fig. 1. This occurs with rate

[TABLE]

and is based on the assumption that mixing of particles in two adjacent regions occurs over the same timescale as motion of particles.

To illustrate the relaxation behavior of the mixture model we choose a system with two energy scales, $J_{0}=1$ and $J_{1}=0.5$. We define the relaxation time, $\tau$, as the time taken for ninety percent of the system to change its spin variable, $n_{i}$, at least once, which represents ninety percent of particles making non-trivial displacements at least once Keys et al. (2011). In Fig. 2(a), we show relaxation times for mixtures with different fractions of the two components. Similar to a single component system, the curves all appear super-Arrhenius, with the mixture relaxation curves lying in between the relaxation curves for the individual components.

An Effective Single Component Model. Based on the super-Arrhenius relaxation behavior of the mixtures in Fig. 2(a) we hypothesize that a mixture effectively behaves as a single, homogeneous material with an effective, predictable energy scale for excitations. This means that mixing of components gives rise to a new, effective dynamical heterogeneity at different length and time scales compared to the original components. To capture the effective behavior of the two component mixture we postulate a single component East Hamiltonian, $H_{\mathrm{m}}$, with an effective energy scale $J_{\mathrm{m}}$. We attempt to predict $J_{\mathrm{m}}$ as a function of the energy scales of the individual component materials. The postulated single component Hamiltonian $H_{\mathrm{m}}$ is given by

[TABLE]

We use a variational method to calculate $J_{\mathrm{m}}$ that best approximates the partition function for the multi-component Hamiltonian $H$ (Eq. 2). Using Jensen’s inequality (sometimes known as the Gibbs–Bogoliubov–Feynman approach Feynman (2018); Chandler (1987)) we have

[TABLE]

where $Z$ and $Z_{\mathrm{m}}$ are the canonical partition functions for Hamiltonians $H$ and $H_{\mathrm{m}}$ respectively, $\langle\ldots\rangle_{H_{\mathrm{m}}}$ indicates the canonical ensemble average with energy $H_{\mathrm{m}}$, and $\Delta H=H-H_{\mathrm{m}}$.

Calculating the value of $J_{\mathrm{m}}$ that maximizes the right hand side of the inequality we obtain

[TABLE]

The effective energy scale $J_{\mathrm{m}}$ can be interpreted as giving rise to effective excitations that occur with probability $\exp(-{\beta}J_{\mathrm{m}})$. Revisiting the Sollich and Evans argument Sollich and Evans (1999, 2003), the average equilibrium domain length between effective excitations is $l_{\mathrm{m}}^{\mathrm{eq}}=\exp({\beta}J_{\mathrm{m}})$, creating a minimum energy barrier of $\gamma J_{\mathrm{m}}\ln l_{\mathrm{m}}^{\mathrm{eq}}=\gamma{\beta}J_{\mathrm{m}}^{2}$. This gives a relaxation time $\tau\sim\exp(\gamma{\beta}^{2}J_{\mathrm{m}}^{2})$, which we refer to as the parabolic law for mixtures. Using the effective energy $J_{\mathrm{m}}$ from Eq. 7 to rescale the data in Fig. 2(a) we find that all the relaxation curves collapse onto a universal parabolic form (Fig. 2(b)). This indicates that the two component system behaves as a single component system with an effective energy scale of $J_{\mathrm{m}}$, an average of the energy scales of the individual components.

Comparison with Experimental Data. With this understanding from the East model we attempt to validate the idea of an effective mixture energy scale with available experimental data for supercooled mixtures. We analyze dielectric relaxation measurements for binary mixtures of various glass-forming liquids Duvvuri and Richert (2004); Wang and Richert (2005); Wang et al. (2010); Puzenko et al. (2005). Relaxation times for thirty three mixtures of five pairs of liquids are shown in Fig. 3(a). We first fit parabolas to relaxation time data for single liquid components to determine $J$ values for the individual components. (We use the low temperature data for these fits because the higher temperature data may contain significant contributions from the non-supercooled liquid.) We then calculate $J_{\mathrm{m}}$ for each mixture using Eq. 7. $\beta^{\mathrm{o}}$ and $\tau^{\mathrm{o}}$ for each mixture are treated as fitting parameters, and $\gamma$ is treated as absorbed into the $J$ value for all experimental data. On scaling the inverse temperature axis by $J_{\mathrm{m}}$, we find that all the single component as well as mixture data collapse onto the parabolic law (Fig. 3(b)). Consistent with the model, it appears that mixtures of liquids effectively behave as a single liquid with an energy scale that is an average of the energy scales of the individual components as per Eq. 7.

Recalling that the onset temperature is the temperature below which spatially and temporally heterogeneous dynamics and super-Arrhenius relaxation behavior occur, it is reasonable to expect that the onset temperature of a mixture lies in between, and proportional to the fraction of, the two individual components. Examining the estimated inverse mixture onset temperatures, $\beta^{\mathrm{o}}$, we find that they are bounded by the onset temperatures of the two individual components and change systematically depending on the fraction of each component, with the exception of sorbitol–glycerol mixtures. In the latter case we find that all the mixture data are best fit using onset temperatures that are lower than both the individual components, even as the curvature of the parabola, $J_{\mathrm{m}}$, is well predicted by Eq. 7. We revisit this pair of liquids in the discussion below.

Modeling Onset Temperatures. In the previous section we determined relaxation behavior by predicting the effective energy scale, $J_{\mathrm{m}}$, of the mixture. We now attempt to determine mixture relaxation times by predicting both the effective energy scale as well the onset temperature. The East model shows super-Arrhenius relaxation behavior for all $T>0$, i.e., it has no meaningful onset temperature. To model materials more realistically the temperature field is modified to $\widetilde{\beta}=\beta-\beta^{\mathrm{o}}$, where $\beta^{\mathrm{o}}$ is the inverse onset temperature of the material Keys et al. (2013). To more effectively predict relaxation behavior for mixtures, where each component can have a different onset temperature, we further modify our multicomponent East model so that depending on the component present each lattice site feels an effective temperature field $\widetilde{\beta}_{0}=1/T-1/T_{0}^{\mathrm{o}}$ or $\widetilde{\beta}_{1}=1/T-1/T_{1}^{\mathrm{o}}$, where $T_{0}^{\mathrm{o}}$ and $T_{1}^{\mathrm{o}}$ are the onset temperatures of the two components. To achieve this we modify the mixture Hamiltonian (Eq. 2) to

[TABLE]

The corresponding partition function, $Z^{\mathrm{o}}$, is given by

[TABLE]

where it can be seen that $\widetilde{\beta}_{0}$ and $\widetilde{\beta}_{1}$ act as effective temperature fields depending on the type of species at lattice site $i$. Postulating again that mixtures of components behave as a single component with effective inverse onset temperature $\widetilde{\beta}^{\mathrm{o}}_{\mathrm{m}}$ and energy scale $J_{\mathrm{m}}$, we propose an effective single component Hamiltonian, $H_{\mathrm{m}}^{\mathrm{o}}$ of the form

[TABLE]

where $\widetilde{\beta}_{\mathrm{m}}=\beta-\widetilde{\beta}^{\mathrm{o}}_{\mathrm{m}}$. The partition function, $Z_{\mathrm{m}}^{\mathrm{o}}$, is given by

[TABLE]

where $\widetilde{\beta}_{\mathrm{m}}$ acts as the effective temperature field for every lattice site $i$. Again using Jensen’s inequality to now find $\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}}$ that best approximates $Z^{\mathrm{o}}$, we obtain

[TABLE]

This equation gives an effective temperature–energy scale for a mixture given the energy scales and onset temperatures of the individual components, and reduces to our earlier result (Eq. 7) at low temperatures.

The effective scale $\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}}$ can again be interpreted as giving rise to effective excitations with probability $\exp(-\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}})$. The average distance between these excitations, i.e., the average equilibrium domain length, is $l_{\mathrm{m}}^{\mathrm{eq}}=\exp(\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}})$, giving rise to an energy barrier of $J_{\mathrm{m}}\ln l^{\mathrm{eq}}_{\mathrm{m}}=\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}}^{2}$. The expected time taken to overcome this barrier, and therefore the relaxation time of the average equilibrium domain length, is $\tau\sim\exp(\beta\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}}^{2}).$ With reference to the relaxation time at the onset temperature, $\tau^{\mathrm{o}}$, we have $\tau/\tau^{\mathrm{o}}\sim\exp(\widetilde{\beta}_{\mathrm{m}}^{2}J_{\mathrm{m}}^{2})$, which is the parabolic form for relaxation of mixtures that takes into account both the energy scales, $J$, and the onset temperatures, $\beta^{\mathrm{o}}$, of the individual components. This form reduces to our earlier result, $\tau\sim\exp({\beta}^{2}J_{\mathrm{m}}^{2})$, at low temperatures.

To compare the new parabolic equation for mixtures with experimental data we use the individual values of $J$ and $\beta^{\mathrm{o}}$ in Eq. 12 to calculate the combined energy–temperature scales. Here we only need one fitting parameter, the reference time scale $\tau^{\mathrm{o}}$. In Fig. 3(c) we find that the relaxation behavior for all liquid mixtures agrees well the parabolic form when the data are scaled by $\widetilde{\beta}_{\mathrm{m}}J_{\mathrm{m}}$ as per Eq. 12, with the exception of mixtures of sorbitol and glycerol. When using Eq. 7 in the previous section, we noted that even though the $J_{\mathrm{m}}$ for sorbitol–glycerol mixtures agreed well with the data, the estimated values for $\beta^{\mathrm{o}}$ were unexpectedly lower than $\beta^{\mathrm{o}}$ of the individual liquids. This, along with the lack of agreement in Fig. 3(c), leads us to suggest further measurements for this pair of liquids.

Predicting Glass Transition Temperatures of Mixtures. Empirical linear interpolations, such as the Gordon–Taylor equation Gordon and Taylor (1952), are typically used to predict the glass transition temperature for mixtures Gordon and Taylor (1952). This approach is not based on a microscopic physical description of glassy dynamics or glass formation, and also does not take into account the protocol dependence of the glass transition temperature. The work presented here provides a prediction for the behavior of supercooled mixtures across various temperatures, from which the glass transition temperature of a mixture can be calculated for a given cooling rate Keys et al. (2013); Hudson and Mandadapu (2018). For a given cooling rate, $\nu$, the glass transition temperature, $T_{\mathrm{g}}$, is the temperature at which the relaxation time becomes slower than the time available to equilibrate at that temperature, i.e., the temperature at which $\nu^{-1}\sim|{d\tau}/{dT}|_{T=T_{\mathrm{g}}}$. Using the parabolic expression $\tau\sim\tau^{\mathrm{o}}\exp(\widetilde{\beta}_{\mathrm{m}}^{2}J_{\mathrm{m}}^{2})$ in the above equation gives a prediction for the glass transition temperature of a mixture.

In summary, we find that a mixture of well-mixed glass-forming liquids behaves as a single material whose dynamics is governed by an effective energy scale, an average of the energy scales of the individual components. We anticipate the use of our effective model in predicting transport properties and glass transition temperatures of mixtures especially for compounds, and under conditions, where making measurements is challenging. Extending this model from binary to polydisperse mixtures may provide insight into the effects of polydispersity on the dynamics of glass-formers Berthier et al. (2016). We also anticipate the use of this model in studying the interplay between dynamical heterogeneity of the system and the diffusion and phase separation of individual components, possibly involving interfacial fronts between two dynamical phases bearing a dynamical interfacial tension Katira et al. (2018). This is especially relevant in multi-component systems such as atmospheric organic aerosols, where melting (and formation) of amorphous solids is hypothesized to occur via an increase (or decrease) in water content Koop et al. (2011).