Suppression of crystalline fluctuations by competing structures in a supercooled liquid

TL;DR

This paper introduces a geometric measure called crystal affinity to predict how local structures in supercooled liquids influence crystallization, linking structural covariance to nucleation rates and glass-forming ability.

Contribution

It proposes the crystal affinity as a new quantitative measure to predict crystallization suppression or promotion in supercooled liquids based on local structural covariances.

Findings

Crystal affinity correlates with nucleation rate.

Anticorrelation with crystalline order reduces crystallization tendency.

Predicts glass-forming ability based on local structure.

Abstract

We propose a geometrical characterization of amorphous liquid structures that suppress crystallization by competing locally with crystalline order. We introduce for this purpose the crystal affinity of a liquid, a simple measure of its propensity to accumulate local crystalline structures on cooling. This quantity is explicitly related to the high temperature structural covariance between local fluctuations in crystal order and that of competing liquid structures: favouring a structure that, due to poor overlap properties, anticorrelates with crystalline order reduces the affinity of the liquid. Using a lattice model of a liquid, we show that this quantity successfully predicts the tendency of a liquid to either accumulate or suppress local crystalline fluctuations with increasing supercooling. We demonstrate that the crystal affinity correlates strongly with the crystal nucleation rate and the crystal-liquid interfacial free energy of the low-temperature liquid, making our theory a predictive tool to determine easily which amorphous structures enhance glass-forming ability.

Equations

$$E=\sum _{\mathrm{structures}i}{c}_i{ϵ}_i$$

$$c_{\mathrm{XS}}(T)=c_{\mathrm{XS},\infty }+\frac{Q}{T}+O(T^{-2})$$

$$Q={\frac{\partial c_{\mathrm{XS}}}{\partial \beta }|}_{T=\infty }=-\sum _{\mathrm{structures}i}{C}_{\mathrm{XS},i}{ϵ}_i$$

$$C_{i,j}=N{\mathrm{Cov}}_{T=\infty }(c_i,c_j)$$

$${\tau }_{\mathrm{nucl}}={\tau }_{\alpha }N\exp \left(\frac{{\sigma }^3}{k_{B}T\Delta F^2}\right)$$

$$\sigma \approx {\left[k_{B}T\Delta F^2\log ({\tau }_{\mathrm{nucl}}/{\tau }_{\alpha })\right]}^{1/3}$$

Full text

The ultimate fate of a liquid upon slow cooling – whether it crystallizes or arrests into an amorphous glass – depends on how easily the crystal can nucleate in the supercooled liquid. It has long been speculated that this kinetic stability with respect to crystallization can be related to the geometrical properties of the local structures which, at the molecular scale, are particularly stable. More specifically, such favoured local structures that are geometrically adverse to crystallinity – for instance due to non-crystalline symmetries – would enhance the stability of the liquid. This idea can be traced back to Frank’s 1952 proposal frank_supercooling_1952 that the stability of icosahedral coordination shells in pure metallic liquids might impede crystallization of close packed cubic crystals. From this starting point, a substantial literature has developed royall_role_2015 ; ma , exploring a variety of approaches to the nature and influence of liquid structure. This includes the study of geometrical frustration in liquids frust , the descriptive study of the distribution of local coordination structures in liquids via computer simulations sims and, recently, nano-focused electron scattering nano , and the search for correlation between specific local structures and the local relaxation rates dyn .

In spite of this considerable activity, explicit evidence of the essential thesis, that there is a correlation between liquid structure and crystallization kinetics, remains sparse. Taffs and Royall paddy2016 have reported on the crystallization of a hard sphere liquid subjected to a bias that favours 5-fold common neighbour coordination. They established a clear dependence of the reduced crystallization time on the magnitude of the bias field, confirming that the liquid structure does indeed influence the rate of crystallization. The clear result of ref. paddy2016 is something of an exception: generally, any adjustment of the Hamiltonian to vary the liquid structure also changes the stability of the crystal structure, along with that of any polymorph favoured by the perturbation. For example, Molinero et al molinero found that adjusting the strength $\lambda$ of the 3-body contribution in a model silicon potential resulted in a maximal glass forming ability that coincided with the value of $\lambda$ corresponding to the crossover in stable crystal phases. As this point coincides with the maximal depression of the freezing point, it is difficult to disentangle the role of liquid structure from that of crystal stability. A similar issue has arisen in experimental efforts to confirm the influence of liquid structure. Lee et al lee studied the kinetics of crystal nucleation in a Ti-Zr-Ni alloy in which the degree of local icosahedral order could be varied with composition. Here the increase in liquid icosahedral order coincided with stabilization of a quasicrystal whose nucleation rate was significantly greater than that of the cubic crystal, effectively concealing whatever influence the liquid structure had on the crystallization of the latter crystal.

In summary, to establish the influence of liquid structure on crystallization kinetics we would like to be able to vary the liquid structure as freely as possible while still ensuring that the equilibrium crystal state remains unchanged. The aim of this article is to present both a model that satisfies these conditions, and a theoretical framework relating the local geometrical properties of stable liquid structures to the crystal nucleation rate. To this end, we utilize a simple lattice model, the Favoured Local Structures (FLS) model ronceray_variety_2011 ; ronceray_2012 ; ronceray_favoured_2015 . We consider binary spins, representing some local conformational degree of freedom of the liquid (e.g. composition), on a face-centered cubic lattice. We define the local structure of the liquid at a given site as the geometrical arrangement of the 12 spins surrounding it. There are 218 rotationally distinct structures at this nearest-neighbor level, and we associate an energy $\epsilon_{i}$ to each site with local structure $i$. Our model therefore has 218 variable energy levels corresponding to the possible local structures. The set of energies $\{\epsilon_{i}\}_{i=1\dots 218}$ constitutes the local energy landscape and entirely characterizes the system’s Hamiltonian. The energy per site is thus

[TABLE]

where $c_{i}$ is the fraction of sites in local structure $i$. We study this model by canonical ensemble simulations using the Monte-Carlo Metropolis algorithm with single spin flips, in a periodic box with typical dimensions $30^{3}$.

In a previous study ronceray_favoured_2015 , we have studied the simplest case where a single $\epsilon_{i}=-1$ for a selected favoured local structure, while all other structures have zero energy. In this case, the energy per site is simply minus the concentration of the FLS, and the ground state is the densest packing of the FLS. Because structures at nearby sites overlap and thus exert constraints on each other, these ground states can be highly frustrated (i.e. the selected FLS cannot fully “tile” the lattice) and the maximum FLS fraction varies between $1/4$ and $1$. As we have shown ronceray_favoured_2015 , this frustration does not prevent crystallization on moderately slow annealing of the system, leading us to conclude that that some kind of competition between structures – in other words, multiple minima in the local energy landscape – is required to stabilize the liquid. Indeed, we have established more recently ronceray_liquid_2016 that favouring simultaneously two or more of the highly frustrated structures can prevent crystallization even at very slow cooling rates, resulting in the dynamic arrest of an amorphous state at low temperature.

In this article, we use a variant of the FLS model adapted to study the role of competing structures in the crystallization of a supercooled liquid, as illustrated in Figure 1a. A crystalline local structure (denoted by “XS”) competes with a small number $n=2\dots 10$ of frustrated competing structures. All these structures are given an energy $\epsilon=-1$, while other structures have zero energy. For simplicity, the crystalline structure is chosen to be the all-up structure, whose ground state is unambiguously a uniform configuration of up spins, with energy $-1$ per site. The mix of competitors is selected such that it does not lead to the formation of alternate metastable crystals, but continuously arrests (in the absence of the XS) into a high energy amorphous state with high energy $E>-0.6$ (see Methods).

An example demonstrates that the geometry of the frustrated competitors can substantially influence the rate of crystallization of the liquid. In Figure 1b. we present two apparently similar sets of $n=6$ frustrated structures (Sets 1 and 2) that we choose as competitors to the all-up crystalline structure. In Figure 1c. we compare the behaviour of $E(T)$ of these two systems on slow cooling with that of a reference liquid in which only the crystalline structure is favoured. In this latter case (black curve) the system exhibits little ordering in the liquid (less than $0.1\%$ of XS) prior to a sharp first-order transition to the uniform ground state (lower inset). Adding the Set 1 of competing structures, the liquid is much more ordered and accumulates up to $12\%$ of FLS, before freezing into the uniform ground state too. In contrast, Set 2 results in similar liquid energetics, but an absence of crystallization: at low temperature, the system arrests into an amorphous state with no trace of long-range ordering (upper inset of Figure 1c.) in which $43\%$ of the sites are in a favoured structure. We estimate the fastest nucleation rate of each of these systems and find that the reference liquid (i.e. that with only the ’XS’ favoured local structure) freezes in $\tau_{\mathrm{nucl}}\approx 7.10^{4}$ MC steps; with Set 1 favoured, freezing is actually slightly faster with $\tau_{\mathrm{nucl}}\approx 4.10^{4}$, while with Set 2 we find that $\tau_{\mathrm{nucl}}>10^{12}$ steps: we never observe crystallization for this system in our simulations.

How can slight geometrical changes in the competing structures slow down crystallization by more than seven decades? The origin of this difference does not lie in the thermodynamic stability of the liquid: they have comparable liquid energies and identical ground states. They also have similar liquid relaxation times ($\tau_{\alpha}<10$ MC steps per site). Their crucial difference is revealed by considering the liquid structure. In Figure 1d. we plot the fraction $c_{\mathrm{XS}}$ of sites in the crystalline structure for each of these systems as a function of $1/T$. On cooling, we find that the Set 1 liquid accumulates significantly more crystalline order than the reference liquid: the competitors actually help the liquid accumulating XS, a behaviour that we label as agonist towards crystallization. The Set 2 liquid, in contrast, antagonizes crystallization by suppressing local crystalline order: in spite of the XS being favoured with energy $-1$, their concentration actually decreases when cooling a liquid from the random infinite-temperature limit. This observation intuitively explains the dramatic difference in crystal nucleation rates in these two systems: crystalline local structures are a necessary precursor to freezing, hence their suppression can be expected to make crystallization extremely difficult.

We now propose a general statistical measure of the liquid structure that can quantitatively differentiate liquids like our Set 1 and Set 2 examples. A high-temperature expansion of the concentration of local crystalline structures

[TABLE]

already captures the distinction between crystal-agonist and -antagonist systems through the value of the first-order coefficient $Q$, that we call crystal affinity of this liquid. Denoting by $Q_{0}$ the affinity of the reference liquid with only the XS favoured, agonist systems have $Q>Q_{0}$ and antagonists $Q<Q_{0}$ (dotted lines in Figure 1d). In a previous work ronceray_liquid_2016 , we have developed an approach, that we termed structural covariance, to compute exactly such high-temperature coefficients in terms of geometrical overlap of structures. Indeed, we have shown the following fluctuation-response-like relation

[TABLE]

where $\beta=1/T$, and the coefficients $C_{\mathrm{XS},i}$ quantifies the geometrical interactions between the XS and all structures $i$ (including the XS itself). These coefficients are equal to the covariances of structural concentrations,

[TABLE]

where $N$ is the system size (making $C_{i,j}$ size-independent), and the covariance is computed at infinite temperature, i.e. on completely random spin configurations. In this regime, correlations between structures only occur when structures overlap and share sites, exerting athermal constraints on each other, as sketched in Figure 2a. These covariances can thus be either measured by statistical analysis of high-temperature liquid configurations, or computed exactly by analyzing such overlaps ronceray_2012 ; ronceray_liquid_2016 . Intuitively, similar structures will tend to overlap well and thus have a positive covariance, while very different structures will tend to exclude one from the other’s vicinity. In the FLS model, we can define a geometrical distance between two structures as the minimum number of spin flips required to convert one local structure into the other. As shown in Figure 2b, the value of $C_{\mathrm{XS},i}$ decreases with geometrical distance between competing structure $i$ and the XS, a trend that is not specific to this choice of crystalline structure (Appendix). What is also clear from Figure 2b is that the geometrical distance does not fully determine the structural covariance. Structures in Figure 2c and d both are $5$ spins flips away from the XS, but the former is a strong agonist while the latter is an antagonist (respectively cyan and red squares in Figure 2b). A glimpse at their geometry explains this difference: the first has all its down spins on one side, attracting XS on the other. The antagonist, in contrast, has its down spins scattered over the structure, such that no XS can overlap with it in any of the two ways shown in Figure 2a. The crystal affinity $Q$ defined in Equation 3 therefore encodes, at a pair interaction level, the specific geometric interactions between the crystalline local structure and its frustrated competitors.

The crystal affinity $Q$ is only a first-order term in Equation 2, and so does not quantitatively capture the crystallite concentration at low temperature in Figure 1d. It does very well, however, in capturing both the degree of antagonism (or agonism) of the liquid and the variation in different liquids in $c_{\mathrm{XS}}$ at low temperature, as shown in Figure 3a. This simple geometrical quantity, therefore, provides a very useful insight into the complex structural fluctuations of the supercooled liquid. Perhaps most striking, the crystal affinity also correlates strongly with the nucleation time of the crystal in this liquid. Indeed, in Figure 3b we show that varying $Q$ in similar liquids leads to variations of $\tau_{\mathrm{nucl}}$ over eight orders of magnitude. We remind the reader that $Q$ contains only information about structural pair correlations obtained in the high temperature limit.

The mechanism by which the crystal affinity controls nucleation rates can be interpreted within the framework of Classical Nucleation Theory (CNT) cnt , in which the bulk free energy gain $\Delta\!F$ in crystallization is opposed by a surface free energy cost $\sigma$ for the interface between the growing crystal nucleus and the supercooled liquid. In this theory, the kinetically limiting step to crystallization is assumed to be the thermally activated formation of a critical nucleus, predicting an average nucleation time

[TABLE]

where $N$ is the number of potential nucleation sites and $\tau_{\alpha}$ is the microscopic relaxation time of the supercooled liquid. There are therefore three possible ways to make crystallization slower: by slowing down the kinetics of the liquid (i.e. increasing $\tau_{\alpha}$), by stabilizing thermodynamically the liquid (i.e. lowering $\Delta\!F$), or by increasing the liquid-crystal surface tension $\sigma$.

We have seen in Figure 3b that the crystal affinity $Q$ strongly influences the nucleation rate. In order to understand which of these aspects are affected by this affinity, we study nucleation in a large number of supercooled liquids with $n=2$ to $10$ competing structures, and selected to sample uniformly the available range of $Q$ values. As shown in Figure 4a., the key factor influencing nucleation times in this data set is the crystal affinity, rather than the number of competing structures. In particular, it is worth noting that for negative values of $Q$ (a situation indicating that crystalline structures, in spite of being favoured, are linearly suppressed in $1/T$ as the liquid is cooled down), crystallization is robustly suppressed and never observed in our simulations. For this whole data set, we analyze in Figure 4b-d. the correlations between $Q$ and each of the three aspects of the CNT nucleation rate. Strikingly, we find that $Q$ exhibits no significant correlation with either the relaxation time $\tau_{\alpha}$ (Figure 4b.) or the free energy difference $\Delta\!F$ (Figure 4c.). It does, however, correlate very strongly with the surface tension $\sigma$ (Figure 4d.): with an $r$ coefficient of $-0.92$, the surface tension is almost fully determined by $Q$. This effect is remarkable as it implies that, within our model, the complex problem of estimating the low-temperature surface tension $\sigma$ – a quantity that is somewhat ill-defined in CNT, as it is not a property at thermodynamic equilibrium, and the concept of surface is ambiguous for microscopic clusters – can be bypassed by measuring the crystal affinity $Q$, a simple, high-temperature structural quantity.

In this article, we have introduced the crystal affinity $Q=\left.\partial c_{\mathrm{XS}}/\partial\beta\right|_{T=\infty}$ as a measure of the propensity of the liquid to accumulate crystalline order on cooling, using the local crystal structure concentration $c_{\mathrm{XS}}$ as an inherent probe of crystallinity. This number encapsulates, in an intelligible way, the geometrical interactions of the crystal with non-crystalline stable local structures that compete with crystalline order in the liquid. Indeed, we have shown with Equation 3 that, in a fluctuation-response-like relation, $Q$ can be computed by analyzing the covariances in the number of crystalline and non-crystalline local structures at high temperature. Furthermore, we have demonstrated that $Q$ provides an excellent predictor of the qualitative trends in rate of crystal nucleation: low affinity implies slow nucleation. More precisely, we have established a clear correlation between $Q$ and the effective crystal-liquid interfacial energy (as obtained from the nucleation time data through the assumption of classical nucleation). This result raises the interesting prospect of a treatment of crystal nucleation in terms of the statistics of structural fluctuations, rather than relying on the awkward imposition of a macroscopic interface. Since the only requirement for Equations 2-4 to hold is that the Hamiltonian of the system can be formulated in terms of a short-range local energy landscape, as in Equation 1, the analysis presented here should provide a quite general framework for the systematic study of structural fluctuations in off-lattice supercooled liquids, and the influence of these fluctuations on the rate of crystal nucleation. In particular, structural covariances could be used as a guide to engineer local energy landscapes that kinetically facilitate or hinder the formation of a target crystal structure.

Acknowledgements. PR is supported by the Prix des Jeunes Chercheurs de la Fondation des Treilles. PH is supported by the Australian Research Council.

Appendix A Methods

Here we provide the details of the rationale and procedures associated with the specific choices of model parameters used in the paper and the properties used to characterise structural fluctuations and nucleation rates. In a previous publication ronceray_favoured_2015_SI , we have presented a detailed description of the FLS model, including a comprehensive survey of the phase behaviour of the model in the cases where only a single structure is favoured. Readers interested in the model are referred to Ref.  ronceray_favoured_2015_SI .

Calculating the crystal affinity.

The high-temperature covariances between concentrations of local structures, as introduced in Equation 4, are displayed in Figure 5. These coefficients are computed exactly using an enumerative algorithm testing for all possible overlaps between pairs of structures, as described in Ref. ronceray_cluster_2016_SI ; ronceray_liquid_2016_SI .

We show in Figure 6 that the coefficients $C_{i,j}$ depend on geometrical distance between structures, regardless of the nature of $i$ and $j$ (not only for the all-up XS selected in the article).

Generating the liquids.

We discuss here the way we generate and select the 107 liquids that are used for the statistical study in Figure 4. These are selected randomly from a large set of possible combinations of the $74$ structures depicted in Figure 7. These structures are those of the $218$ local structures in the FLS model that have the following properties:

• •

a ground state energy higher or equal than $-1/3$, ensuring that the structure is frustrated and does not permit formation of crystal polymorphs that would compete with the all-up ground state;

• •

a number of “up” spins between $4$ and $8$, so that the interactions between structures rely on geometrical properties rather than on mere ferromagnetic attraction or repulsion.

There are $19$ distinct types of local structures (columns in Figure 7) fulfilling these conditions, coming in up to four variants when including the mirror- and/or spin-inverted variants. For each value of the number of competing local structures $n=2...10$, we generate $10^{4}$ different Hamiltonians by picking randomly one structure in $n$ randomly selected columns of Figure 7 (We don’t use two structures of the same column as we observed empirically that this tends to increase crystallization of composite metastable crystals.). We then pick samples that are regularly spaced in values of $Q$, so as to investigate the correlation between $Q$ and thermodynamic observables. These samples are selected by binning the $10^{4}$ systems by values of $Q$, with interval $1.5\times 10^{-5}$. We then pick one system in one bin out of three, thus obtaining a first batch of 107 Hamiltonians (not all values of $(Q,n)$ are possible).

We are interested in nucleation of the target crystal structure (the all-up system) within a liquid where the competitors are present too. We thus run a pre-selection simulated annealing, with only the frustrated competing structures favoured (no XS), to eliminate those of these systems that tend to crystallize into a polymorph crystal made of a combination of the frustrated competing structures, and thus wouldn’t remain liquid at temperatures near the point at which we study nucleation, $T=0.6$. We thus eliminate and re-pick the system to obtain another Hamiltonian with similar $Q$ value, if, on slow annealing:

• •

the system exhibit a sharp first-order phase transition apparent as a discontinuity in energy per site $E(T)$; or

• •

the system exhibits a peak of heat capacity at $T>0.45$ (sign of either a weakly first order phase transition, or a second order transition that would result in (near-)critical slowdown of the kinetics, and hence affect the nucleation rate in an undesired way; or

• •

the system exhibits obvious long-range crystalline ordering in the final, low-temperature state.

Once these pre-selection tests and repicking are performed, we have obtained a large set of 107 Hamiltonians composed of a mixture of the structures depicted in Figure 7, regularly spaced in $Q$ and $n$ values, that do not, when only these structures are favoured, crystallize on cooling. This phenomenology is further described in Ref.ronceray_liquid_2016_SI in the case of mixtures of $2$ frustrated FLS’s. Note that the $13$ systems presented in Figure 3 are simply the $n=7$ line in Figure 4.

Measuring nucleation times.

We can next measure nucleation times for each of these Hamiltonians. Our protocol is the following: we prepare the liquid in a system of $20^{3}$ sites with periodic boundary conditions, by annealing a similar system, but without favouring the XS so as to prevent premature nucleation. Having such an equilibrated liquid, we check that it hasn’t crystallized into a frustrated polymorph (i.e. that the system’s energy is consistent with the corresponding liquid), then add the XS to the set of favoured structures. We then perform MC moves at a fixed temperature $T=0.6$ (chosen to be the temperature of fastest nucleation, in a quite robust way) until crystallization is detected, i.e. when the energy falls below a threshold energy of $-0.9$, from which the system always falls into the ground state. This protocol is however spoiled by the presence of ballistic relaxation, both before nucleation (relaxation of the liquid after adding the XS) and after nucleation (crystal growth). For short nucleation times, this makes the measured time before hitting the threshold strongly reproducible. We thus obtain an upper bound to the nucleation time by repeating the simulation several times (up to 20 times or $10^{12}$ MC steps) and measuring the standard deviation of the crystallization time, which gives us an estimate of $\tau_{\mathrm{nucl}}$, assuming that nucleation times are exponentially distributed (i.e. no aging or memory in the liquid prior to nucleation).

Liquid relaxation times.

The relaxation times $\tau_{\alpha}$ of the liquid (as presented in Figure 4b.) are the autocorrelation times of the energy of the liquid, expressed in Monte-Carlo steps per site. They are obtained using the batch autocorrelation method.

Free energies.

Liquid-crystal free energy differences presented in Figure 4c. are computed by approximating the crystal’s free energy by its energy value $E\approx-1$, neglecting its entropy. An approximation of the liquid free energy is obtained by thermodynamic integration from infinite temperature of the system without the XS favoured, so that it does not crystallize in the annealing run performed to integrate the free energy. Neglecting the XS leads to a minor error in the free energy of the liquid (less than $1\%$ in systems for which crystallization is slow enough to allow measurement of the free energy with the XS favoured).

Surface tensions.

We do not measure directly the surface tensions presented in Figure 4d.; rather, we infer it through Classical Nucleation Theory assumption, using

[TABLE]

where $k_{B}=1$ in our unit system, $T=0.6$ for all nucleation runs presented in this article, and $\tau_{\mathrm{nucl}}$, $\tau_{\alpha}$ and $\Delta\!F$ are obtained as discussed above. Note that we do not include the volume factor $N=20^{3}$ in Equation 6; doing so would result in undefined values for $\sigma$ for many of the high-affinity liquids, as nucleation is often faster than $N\tau_{\alpha}$ in these cases. In any case, $\sigma$ in this article should not be interpreted too strictly as a surface tension – as discussed at the end of the main text, this notion is somewhat ill-defined, on top of being particularly difficult to measure directly. Rather, it quantifies the aspect of the nucleation rate that is not a consequence of the bulk free energy gain $\Delta\!F$ in crystallizing, as illustrated in Figure 8 by showing the near absence of correlation between $\sigma$ and $\Delta\!F$.

Influence of the number of competitors.

In Figure 4, we present a large data set of liquids with varying number of competing structures $n$ and crystal affinity $Q$. Panels b-d. of this figure analyze the influence of $Q$ on the different parameters $\tau_{\alpha}$, $\Delta\!F$ and $\sigma$ that affect the nucleation time, showing that $Q$ strongly correlates with $\sigma$ while being essentially independent of the other parameters. As a complement, we present in Figure 9 the same correlations, but with the other parameter of our data set, the number of frustrated competing structures $n$. It presents a strong correlation with $\Delta\!F$ ($r=0.94$), while showing no significant correlation to $\sigma$. These results can be intuitively understood by noting that adding more competitors to the liquid increases the number of ways for the liquid to lower its energy: it will thus increase its entropy, at a given energy, and thus lower the free energy of the liquid, while leaving the crystal unaffected. On the other hand, these additional competitors are equally likely to act as agonist or antagonist with the crystalline structure, hence they do not significantly affect the crystal affinity $Q$ – and, therefore, do not affect the surface tension either.