Glass Transition in Supercooled Liquids with Medium Range Crystalline Order

TL;DR

This study uses molecular dynamics simulations to explore how medium range crystalline order influences the glass transition, revealing distinct properties and validating a robust method for measuring static length scales.

Contribution

It demonstrates that glasses with medium range crystalline order have different static and dynamic properties and confirms the Point-to-set method as a reliable tool for measuring static correlation lengths.

Findings

Glasses with MRCO exhibit distinct static and dynamic behaviors.

The Point-to-set method is validated as an order-agnostic measure of static length scales.

Differences identified suggest different molecular mechanisms in glasses with and without MRCO.

Abstract

The origins of rapid dynamical slow down in glass forming liquids in the growth of static length scales, possibly associated with identifiable structural ordering, is a much debated issue. Growth of medium range crystalline order (MRCO) has been observed in various model systems to be associated with glassy behaviour. Such observations raise the question about the eventual state reached by a glass former, if allowed to relax for sufficiently long times. Is a slowly growing crystalline order responsible for slow dynamics? Are the molecular mechanisms for glass transition in liquids with and without MRCO the same? If yes, glass formers with MRCO provide a paradigm for understanding glassy behaviour generically. If not, systems with MRCO form a new class of glass forming materials whose molecular mechanism for slow dynamics may be easier to understand in terms of growing crystalline order, and should be approached in that manner, even while they will not provide generic insights. In this study we perform extensive molecular dynamics simulations of a number of glass forming liquids in two dimensions and show that the static and dynamic properties of glasses with MRCO are different from other glass forming liquids with no predominant local order. We also resolve an important issue regarding the so-called Point-to-set method for determining static length scales, and demonstrate it to be a robust, order agnostic, method for determining static correlation lengths in glass formers.

Figures

Equations

$${\psi }_6=\frac{1}{N}\sum _{i=1}^N{\psi }_6^i,\text{mbox}where{\psi }_6^j=\frac{1}{n_j}\sum _{k=1}^{n_j}\exp \left(ı6{\theta }_{jk}\right)$$

$$g_6(r)=⟨{\psi }_6^{\ast }(\vec{r}){\psi }_6(0)⟩,$$

$$S_4(q,t)=N\left[⟨\tilde{Q}(q,t)\tilde{Q}(-q,t)⟩-⟨\tilde{Q}(q,t){⟩}^2\right].$$

$$\tilde{Q}(q,t)=\frac{1}{N}\sum _{i=1}^N\exp (i\mathbf{q}.{\mathbf{r}}_i(0))w\left(|{\vec{r}}_i(t)-{\vec{r}}_i(0)|\right),$$

$${\chi }_4(t)=N\left[⟨Q^2(t)⟩-⟨Q(t){⟩}^2\right].$$

$$Q(t)=\frac{1}{N}\sum _{i=1}^{N}w\left(|{\vec{r}}_i(t)-{\vec{r}}_i(0)|\right),$$

$$S_4(q,\tau )={\chi }_0\mathcal{F}(q{\xi }_4),$$

$$B(N,T)=1-\frac{<Q(N,T{)}_{\tau }^4>}{3<Q(N,T{)}_{\tau }^2{>}^2},$$

Full text

glass transition — static length scale — dynamic length scale — hexatic order — Medium Range Crystalline Order

The existence of growing length scales and their connection to the rapid slowing down of dynamics of glass forming liquids approaching glass transition are still a subject of active research. Based on extensive research on the existence of various static and dynamical length scales in glass formers, it is now understood that multiple growing length scales exist and correlate with increasing relaxation time scales as the glass transition is approached, but whether a single or a small number of such length scales satisfactorily describe the essential physics of dynamical slow down and the glass transition is not fully resolved 11BB ; arcmp ; ROPPSMA ; annurev.physchem.51.1.99 ; 05Berthier ; RFOT2 . For quite sometime in the study of glass formers, there has been a focus on identifying structural motifs which are favoured locally but which can not span the whole system due to frustration arising from the incommensurate nature of the structural motifs. There are many experimental and simulation studies which claim the existence of such structural motifs in a variety of glass forming liquids paddyRoyal1 ; paddyRoyal2 . A strong connection between the growth of such locally favoured structures (also widely known as medium range crystalline order, MRCO) and slow relaxation has also been proposed in some glass formers, with a fair degree of supporting evidence PhysRevLett.99.215701 ; Nat2010 ; J.Phys.Condens.Matter23-194121 . In Figure 1, snapshot of such MRCO (red regions in the figure) is shown for one of the model glass formers studied in this work. While the dynamics of different glass forming liquids are very similar to each other in essential aspects, only some glass formers manifest growing MRCO, whereas such structural motifs can not be clearly identified in many others. The failure to identify growing structural order may either be due to their absence or due to the subtlety of the relevant structural order that may not be captured by the metrics employed. It is thus important to understand whether the physics that govern the slowing down of dynamics in different glass forming liquids are generically the same irrespective of the presence of well characterised MRCO or whether essential differences can be found between those systems with MRCO and those without. If it turns out that indeed the physics governing the slowing down in different glass forming liquids are same and all glass forming liquids eventually will start to have locally preferred structurally ordered domains, then the complexity of glass transition will be reduced to the formation of such domains and efforts must be directed at detecting through suitable metrics the growing structural order. On the other hand if it turns out that the behaviour of glass formers with MRCO is not generic, then they form a distinguishable special class of glass forming systems, whose behaviour is perhaps less complex, and must be analysed differently.

Recently, in a series of papers Tanaka et. al. PhysRevLett.99.215701 ; Nat2010 ; J.Phys.Condens.Matter23-194121 showed that the glass transition in systems with MRCO is determined by the length scale associated with that of locally preferred structures. For example, for two dimensional systems, the correlation length associated with hexatic order is the important length scale and the glass transition in these systems are completely controlled by the hexatic correlation length. In PNASUSA2015 , a polydisperse two dimensional glass forming liquid was studied in which the hexatic correlation grows rapidly with increasing packing fraction. The dynamic heterogeneity PhysRevLett.79.2827 ; PhysRevE.58.3515 ; science.1166665 length scale obtained from the scaling analysis of the wave-vector ($q$) dependence of the four-point structure factor, $S_{4}(q,t)$ JCP119-14 was compared with the hexatic correlation length obtained from the spatial correlation function of the hexatic order parameter, $\psi_{6}$ J.Phys.Condens.Matter23-194121 ; PhysRevE.75.041503 ; JCM14 defined as

[TABLE]

where $\theta_{jk}$ is the angle made by the vector from particle $j$ to its neughbour $k$ with the $x$-axis and $n_{j}$ is the number of neighbours of particle $j$ (see SI for further details). These two length scales were found to be proportional to each other over the studied density range.

The static length scale obtained from few different methods including the Point-to-set (PTS) Nature2008 ; PhysRevLett.108.225506 method was also compared with the hexatic correlation length. The PTS length scale was found to be very different from the hexatic correlation length in the studied density window. In this system with local hexatic order, the behaviour is completely determined by the hexatic correlation length. So, if the PTS method is “order agnostic”, then it should pick up the same hexatic correlation length. From the large difference between the calculated values of the PTS and hexatic length scales, it was concluded that the PTS method fails to pick up the relevant static length scale for systems with MRCO. These results cast a serious doubt about the usefulness of the PTS method in measuring important static length scales relevant for slow dynamics in all glass forming liquids.

In this work we have addressed the important issue of whether a deep connection between structure and dynamics exists for all glass forming liquids. As noted above, a number of theoretical approaches ROPPSMA ; RevModPhys.83.587 ; PhysRevLett.97.195701 ; ARPC58 ; AIP2003 have explored the relevance of a growing length scale to dynamical slow down. A specific goal in some of these approaches were to study the connection between heterogeneous dynamics (also referred to as dynamic heterogeneity) in glass formers annurev.physchem.51.1.99 ; J.Phys.Soc.Jpn ; PhysRevE.52.1694 ; PhysRevLett.80.2338 with rapid change in viscosity. We show that the slowing down in the dynamics while approaching the glass transition in liquids which exhibit a pronounced increase in structural order is strongly coupled to the growth of dynamic heterogeneity length scales, as well as to the growth of the length scale that captures the local structural order. In two-dimensional glass formers with MRCO, the static length scale determined by growing hexatic order is the same as both the dynamic heterogeneity length scale and a static amorphous length scale which we independently evaluate using the PTS method. On the other hand, for liquids in which such structural ordering is not prominent, the temperature dependence of the dynamic heterogeneity length scale is very different from that of the static amorphous length scale. This suggests that the glass transition in liquids which show a strong tendency to the growth of locally favoured structural order is fundamentally different from the transition in liquids in which such local crystalline structures do not grow very prominently.

We address these issues by performing extensive computer simulations of four different glass forming liquids in two dimensions and measuring all the relevant length scales in each system using a wide range of procedures. The model systems studied are: (i) the two-dimensional Kob-Andersen binary Lennard-Jones mixture (referred here as 2dKA) PhysRevE.51.4626 at fixed density, with the temperature as the variable that determines changes in dynamics; (ii) a two- dimensional system characterized by a repulsive inverse power-law potential, referred to as 2dIPL PhysRevLett.105.157801 , studied with density as the relevant variable. For the above two systems, structural order or MRCO grows significantly while approaching the putative glass transition. We study another system (iii) with a repulsive interaction potential (referred to as 2dR10) PhysRevE.82.031301 . No pronounced crystalline order is found in the studied temperature range for this model. This is consistent with previous work PhysicaA that showed that the 2dR10 model has almost no tendency to exhibit locally favored crystalline order. We also study (iv) the polydisperse system studied by Tanaka et al PNASUSA2015 , for reasons which we discuss later. The details of the models and simulations are given in the SI.

Rest of the paper is organized as follows: first, we briefly discuss the systems with local hexatic order (2dKA, 2dIPL) and compare various length scales in those model systems. We then discuss the model (2dR10) for which there is no prominent hexatic order and compare the findings with those for the previous set of model systems. Then we will discuss the controversy regarding the PTS method and describe results that resolve that controversy. Finally, we conclude with a summary of the perspective regarding the relationship between structural order and dynamics in different classes of glass formers.

For the 2dKA model, hexatic order grows with decreasing temperature as illustrated in Figure 1 where particles are coloured according to the value of the local hexatic order parameter, $\psi_{6}^{i}$. Similarly, hexatic order grows for the 2dIPL model with increasing density. We have extracted the correlation length associated with the hexatic order by fitting the decay profile of the peaks of the hexatic correlation function, $g_{6}(r)$ JOPCMEPP defined as

[TABLE]

where $\psi_{6}(\vec{r})=\sum_{i=1}^{N}\delta(\vec{r}-\vec{r}_{i})\psi_{6}^{i}$. We have normalized it by pair correlation function $g(r)$ (see SI for further details). We have used exponential function and the two dimensional version of the Ornstein-Zernike function, to fit the decay of the peaks of the correlation function to get the correlation length. The correlation length $\xi_{6}$ is plotted in Figure 2 (a), 2 (b) and 2 (c) together with other dynamical and static correlation lengths which are described below.

Previous studies have shown Nat2010 ; J.Phys.Condens.Matter23-194121 a one-to-one correspondence between dynamic heterogeneities and regions of high/low hexatic order. We also find the same correspondence for the model systems studied here where medium range hexatic correlations are prominent. We have calculated the dynamical length scale $\xi_{4}$ from the four point structure factor $S_{4}(q,t)$ JCP119-14 ; PhysRevLett.97.195701 , which is defined as

[TABLE]

with $\tilde{Q}(q,t)$ is defined as

[TABLE]

where the window function, $w(x)=1$ if $x<0.3$ and [math] otherwise.

A dynamical susceptibility $\chi_{4}(t)$ is defined in terms of the fluctuations of two point overlap correlation function $Q(t)$ as

[TABLE]

The function $Q(t)$ (see SI for further details) defined as

[TABLE]

This function measures the overlap of a configuration of particles at a given initial time $(t=0)$ with the configuration at a later time $t$. Moreover $\chi_{4}(t)\equiv\lim_{q\to 0}S_{4}(q,t)$. The $\alpha$ relaxation time $\tau_{\alpha}$ is defined as $Q(t=\tau_{\alpha})=1/e$. We define another time scale $\tau$ as the time at which $\chi_{4}(t)$ peaks - it is proportional to $\tau_{\alpha}$ with a proportionality constant close to

1. We consider $S_{4}(q,t)$ evaluated as $t=\tau$ in our analysis using the Ornstein-Zernike form,

[TABLE]

to obtain the dynamical correlation length, $\xi_{4}$. Here, $\chi_{0}=S_{4}(q=0,\tau)$ is the $q=0$ value of the structure factor and is same as $\chi_{4}(\tau)$ only in the Grand canonical Ensemble ROPPSMA ; PhysRevLett.105.015701 ; PhysRevLett.105.019801 ; PhysRevLett.105.217801 , and $\mathcal{F}(x)$ is a scaling function (see SI for further details of the scaling analysis).

We have also estimated the dynamical correlation length $\xi_{4}$ from detailed finite size scaling analysis of the Binder Cumulant ROPPSMA ; J.Phys.Condens.Matter23-194121 ; PNASUSA2009 ; Zeitschrift , defined as

[TABLE]

where $Q(N,T)_{\tau}\equiv Q(t=\tau)$ for system size $N$ at temperature $T$. The Binder cumulant is defined in terms of $4^{th}$ and $2^{nd}$ moments of the distribution of $Q$ so that it becomes zero for a Gaussian distribution and acquires positive values for bimodal distributions. It is a scaling function of only $L/\xi_{4}$ where $L$ is the linear dimension of the system (further details can be found in SI). From Figures 2.(a),  2.(b) and  2.(c) it is clear that this length scale obtained from finite size scaling of $B(N,T)$ is in good agreement with the dynamical length scale calculated from the four-point structure factor $S_{4}(q,\tau)$. For one model (2dKA) we also compute the dynamical length scale from the four-point structure factor of bond breakage correlation function (see SI for details), and as seen in Figure. 2 (a), this length scale also compares very well with the lengths obtained by the other methods.

The growth of the calculated dynamical correlation length, $\xi_{4}$ is found to be in good agreement with the hexatic correlation length, $\xi_{6}$ for models 2dKA and 2dIPL. This is in complete agreement with the observation made in Ref.PNASUSA2015 where it was shown that for systems with growing hexatic order, the hexatic correlation length is the same as the dynamical heterogeneity length scale.

To understand possible connections between the growing hexatic order and static length scales associated with amorphous order, we have calculated the static length scale using three different methods: (a) PTS correlation method, (b) Finite size scaling of the relaxation time $\tau_{\alpha}$ and (c) Finite size scaling of the minimum eigenvalue of the Hessian matrix PhysicaA calculated at local minima of the potential energy (known as Inherent structures) visited by the system at a given temperature. Computing PTS correlations Nature2008 ; PhysRevLett.108.225506 is an elegant, general method to capture the multi-point static structural correlations in viscous liquids. The main idea is to freeze the majority of particles from an equilibrated configuration of the liquid, outside a specified volume, and to measure how the structure of the unfrozen particles in the inner volume are affected. It has been argued that in the geometry in which particles outside a spherical cavity are frozen, the PTS correlation should detect the typical domain size predicted by the random first order transition theory (RFOT) JChemPhys.2013Mar28 ; JCP121-7347 . Here we calculate the point to set length scale (referred to as $\xi_{PTS}$) in the cavity geometry. Bulk equilibrium configurations are generated at desired temperature and density using NVT Molecular dynamics simulations and then cavities are constructed by freezing the particles outside a circular region of radius $R$. Static overlap correlations are calculated as a function of cavity radius $R$ and from the decay profile we obtain the PTS length scale (see SI for further details).

The static length scale obtained from the finite size scaling of the relaxation time is referred to here as $\xi_{\tau}$ arcmp ; PNASUSA2009 . The details of the analysis and the scaling plots are given in the SI. The length scale $\xi_{\lambda}$ obtained from the scaling of the minimum eigenvalue of the Hessian matrix is calculated only at low temperatures due to the harmonic nature of the analysis which will not be a good approximation at higher temperatures (see the SI for further details of the analysis and scaling plots for different model systems).

The comparison of different length scales are summarized in Figures  2.(a),  2.(b) and  2.(c). The data for 2dKA and 2dIPL in which MRCO grows significantly, both dynamic and static length scales obtained using different methods grow very similarly as the hexatic correlation length. On the other hand for the 2dR10 model where no significant hexatic order is observed in the studied temperature range, the dynamic length scales grow much faster than both the static length scale and the hexatic length scale. Note that for the 2dR10 model, the growth of the hexatic correlation length is even smaller than the modest growth of the static length scale.

These results clearly show that the relationship between dynamical, static and hexatic order length scales is significantly different for glass formers with MRCO (for which all length scales are comparable) and those without (for which different length scales are significantly different).

We next consider the relationship between the dynamical length scales and relaxation times. In Figure.3 we plot the relaxation time (scaled by its value in the infinite size limit) against the Binder cumulant for the three different models, and for different temperatures (or densities for the 2dIPL model) and system sizes in each case. If the system studied exhibits dynamical finite-size scaling, the data for different system sizes and temperatures should fall on a master curve. Further, since we expect the relaxation time to be governed by a static length scale (see above), such a data collapse would also indicate the equality of static and dynamical length scales, which is not generally expected to hold. Figure. 3 shows $\tau_{\alpha}(N,T)/\tau_{\alpha}(N\rightarrow\infty,T)$ vs. $B(N,T)$ for 2dKA and 2dR10 model and $\tau_{\alpha}(N,\rho)/\tau_{\alpha}(N\rightarrow\infty,\rho)$ vs. $B(N,\rho)$ for the 2dIPL model systems we study. The top two panels, for 2dKA and 2dIPL, show convincing data collapse, indicating that a unique length scale, which we may attribute to growing hexatic order, controls the dynamics, as well as the static order. In contrast, the bottom panel, for 2dR10, shows no data collapse, and confirms our previous observation in Figure. 3 that the growths of dynamic and static length scales are different and decoupled from each other in this system.

However, our conclusion about the equality of static, dynamic and structural length scales in systems exhibiting prominent hexatic order has recently been questioned. In Ref. PNASUSA2015 , it was argued that the PTS method PhysRevLett.116.098302 ; PhysRevLett.111.165701 fails to capture the growth of the relevant static length scale for systems which exhibit growing MRCO. The system studied in Ref.PNASUSA2015 is a polydisperse one with $11\%$ polydispersity. In the PTS method one calculates the static overlap as a function of different cavity radii (R) for a given temperature or density and then from the decay of the static overlap with increasing cavity radius one extracts the static length-scale.

For randomly chosen small cavities, the packing fraction and polydispersity will be distributed around the mean values, with the variance being larger for smaller cavities. It may be expected, however, that for a sufficiently large number of realisations, such that the average packing fraction and polydispersity lie within a reasonably small tolerance, the estimated properties (in this instance, the static correlations to obtain the PTS length scale) will converge to the correct values. We thus estimate the number of realisations required for the average polydispersity and packing fraction to lie within $2\%$ of the bulk values, different values of the PTS length are obtained. In the inset of the left panel of Figure. 4 we show such an estimate for a system of $N=10000$ polydisperse particles in two dimension with $11\%$ polydispersity (henceforth we will refer to this model as 2dPoly).

One can clearly see that for small cavities, one needs averaging over as large as $250$ different cavities in order to reach the bulk polydispersity. The distributions of the polydispersity and packing fraction obtained for different cavity sizes are shown in the left and middle panels of Figure. 4.

We then obtained the PTS length scale by averaging over the estimated number of cavities for each radius. The results obtained this way are shown in the right panel of Figure.4 (blue triangles). Similarly to the results in Ref.PNASUSA2015 , we find that the PTS length scale so estimated falls well below the hexatic length scale and other static and dynamical length scales that track the hexatic length scale. This result is not expected to change if we average over a larger number of cavities for each $R$.

If, on the other hand, we adopt, following Nature2008 ; PhysRevLett.108.225506 , the procedure of evaluating the PTS length scale by selecting only those cavities for any radius for which the packing fraction and polydispersity within the cavity fall within a cutoff of $2\%$ from the bulk values. The PTS length scales obtained from this procedure (magenta triangle) match remarkably well the hexatic order length, the dynamical heterogeneity length scale and the static length scale estimated from finite size scaling of the relaxation time. We believe this procedure to be more appropriate because thermodynamic parameters such as temperature, density and polydispersity of the particles inside the cavity are assumed to be fixed in the theoretical arguments Nature2008 ; PhysRevLett.108.225506 that relate the PTS length scale to the mosaic scale of RFOT. We thus conclude that the difference between the PTS length scale and the hexatic order length scale observed in Ref. PNASUSA2015 arises from the method of calculation, rather than from any fundamental limitations of the PTS method in determining static length scales. The PTS method appears to be a robust, order agnostic approach, and while the issues surrounding what constitutes a proper way of calculating $\xi_{PTS}$ need to be better understood, there may be no need to modify it, as recently proposed PhysRevE.94.032605 .

In conclusion, we have shown that the growth of different measures of static and dynamic length scales in supercooled liquids with medium range crystalline order (MRCO) is fundamentally different from the generic behaviour of glass forming liquids where local crystalline order is not prominent. For liquids with MRCO, the growth of dynamic and static length scales tracks very closely the MRCO length scale. The growth of relaxation times in these systems are apparently controlled by the growth of the MRCO length scale. On the other hand, for generic glass formers where the growth of MRCO as temperature or density is varied is not prominent, the dynamic length scales grow more rapidly than static length scales and the increase in relaxation times is controlled by the static length scale. Finally we have conclusively shown that the point-to-set correlation method is indeed order agnostic and captures the relevant length scale for all model systems including those with prominent MRCO.