Excitations are localized and relaxation is hierarchical in glass-forming liquids
Aaron S. Keys, Lester O. Hedges, Juan P. Garrahan, Sharon C. Glotzer, David Chandler
I Introduction
For each of five atomistic models at several different densities and temperatures, we ask the question: How does an atom move a distance between relatively long lived neighboring positions? We answer this question by augmenting extensive molecular dynamics simulations with methods of transition path sampling Bolhuis et al. 2002. In particular, molecular reorganization in a supercooled liquid is a rare event, and transition path sampling can harvest unbiased ensembles of trajectories exhibiting such events. By applying this approach, we find several important results.
Second, for a given displacement length , we find that the equilibrium concentration of excitations, , has a Boltzmann temperature dependence
where grows logarithmically with , i.e.,
Here, is the space-filling diameter of a particle, the values of and are material-dependent, and is of order unity.
Third, we find that observed super-Arrhenius temperature variation of transport properties follows directly from this logarithmic scaling of . In particular, as the average distance between excitations grows with lowering , the energy scale for relaxation also grows with lowering . Specifically, from Eq. (2), we argue that the structural relaxation time, , obeys
I.2 Results demonstrate the perspective of kinetically constrained models.
Equation (3) is the temperature dependent form for that successfully collapses disparate transport properties of all fragile glass formers Elmatad et al. 2009; Elmatad et al. 2010. The logarithmic scaling from which it results, Eq. 2, is that of the East model Jäckle and Eisinger 1991 and higher dimensional generalizations Garrahan and Chandler 2003; Berthier and Garrahan 2005. Here, by providing a microscopic recipe for detecting excitations, computing and predicting , we validate East-like KCMs as good models for the dynamics of atomistic glass formers.
II Methods
The abundance of results collected in this paper are unprecedented by the current standards of this field. These results establish the behaviors of many distinct large systems, and demonstrate the commonality of these behaviors over a broad range of time scales. Even so, further studies could provide additional documentation and perhaps refinements of the conclusions we are able to draw from the results presented herein.
The atomistic models we have simulated are two-component mixtures for which particles of type interact with those of type with pair potentials, , where is the separation of the pair. These potentials are shifted and truncated Lennard-Jones potentials,
The first three models detailed in Table 1 have been developed and studied by others before us. The last two are modifications of models examined by Harrowell and co-workers Hurley and Harrowell 1995; Candelier et al. 2010a. In particular, we replace the inverse-twelve power potential of Refs. Hurley and Harrowell 1995 and Candelier et al. 2010a with Weeks-Chandler-Andersen Weeks et al. 1971 (WCA) repulsive potentials.
By exploring the dynamics of these models for larger system sizes and longer times than have been considered before, we find that both variants of the 2D-50:50 system appear to coarsen and freeze. As such, we use those systems for illustrative qualitative studies over timescales smaller than coarsening times. For quantitative behavior of reversible transport and related properties, we use the other four models, none of which exhibit coarsening at the conditions we consider. Our molecular dynamics for those systems is carried out with , and fixed, using the HOOMD-Blue simulation code Anderson et al. 2008 on graphics processors.
II.2 Excitations are manifested by particle displacements that persist for significant periods of time.
While glassy systems tend to be locally rigid or jammed, some locations are atypically less rigid and provide opportunities for structural reorganization. These locations coincide with defects or excitations in an otherwise jammed material. Their presence can be inferred or detected by observing non-trivial particle displacements associated with transitions between relatively long-lived configurations. Because these displacements are reliable recorders of excitations, we often treat them as synonymous. Nevertheless, the two are distinct. Displacements refer to dynamics in small parts of trajectories, while excitations refer to underlying configurations or micro-states. We further comment on this point later in this paper.
Non-trivial displacements – the recorders of excitations – are different than fleeting vibrations. Inherent structures Stillinger and Weber 1984; Heuer 2008 help distinguish the former from the latter. The inherent structure of a configuration of an -particle system is obtained by steepest descent to the nearest minimum in the potential energy landscape Bitzek et al. 2006. The inherent structure evolves as dynamics progresses, but most intra-basin vibrations are not apparent in that evolution. Excitations are present when the dynamics produces significant changes in inherent structure. The same effect of removing relatively high-frequency vibrations can be obtained by coarse graining coordinates over a small increment of time. We have generally chosen this latter procedure. But we have checked in a few cases that our principal results for large enough length scales and time scales are unaffected by this choice, provided the coarse-graining time is of the order of the typical period for atomic length scale intrabasin vibrations and this time is of the order of, but generally smaller than, the typical instanton time , defined below. Where denotes the instantaneous position of the th particle in the system at time , we use to denote its corresponding position in the time-coarse-grained structure, i.e.,
Figure 1 renders part of a trajectory from one of the models we have considered. It illustrates typical behavior of these processed coordinates. While inherent structures change discontinuously as the system moves between potential energy basins, the raw and coarse-grained coordinates are continuous. The particular trajectory illustrated in Fig. 1(a) shows the appearance of the type of motion with which we detect an excitation. It is characterized by three periods: During the first, the inherent structure of a tagged particle remains near its initial point; following this sojourn, the system undergoes a rapid transition that lasts up to a time , termed the “instanton” time, after which the tagged particle is found in a distinctly different position for another relatively quiescent period.
This behavior suggests using the following functional of path as the indicator that at time particle is associated with an excitation with displacement length :
Here, or 0 for or , respectively. The products are over time steps of a trajectory that extends for a time . This length of time is at least as long as the time for trajectories to commit to one basin or the other while traversing a transition state region. In the parlance of rare-event dynamics, is a plateau or commitment time Chandler 1978; Hänggi et al. 1990 – a time for which the equilibrium average of grows linearly with . It is typically 3 or 4 times the mean instanton time for that displacement length , . As such, the indicator given by Eq. 5 discards trajectories that exhibit displacements with length scale if they fail to persist for at least as long as the typical instanton time. The instanton time is the shortest time separating the initial and final sojourns. It is a functional of path, and it varies from one trajectory to another.
Fig. 1(b) illustrates the spatial arrangements of these excitations. In particular, at supercooled temperatures, for reasonable choices of , the pictures give the impression that the mean excitation concentration is small. As temperature is lowered, typical sizes of excitations are unchanged, but their concentration clearly decreases. This concentration is recorded by
where the angle brackets denote equilibrium ensemble average, is the volume and is the total number of particles in the system. This second (approximate) equality in Eq. 6, where is replaced by , holds provided is indeed a plateau time. At conditions where the liquid is not supercooled a separation of time scales need not exist, and thus there is no plateau behavior. We replace -functions with the -functions when calculating at supercooled conditions.
This quantity is the average number of displacements per unit space-time. We refer to it as the average excitation concentration. With transition path sampling, without the simplifying approximation of Eq. 6, we collect satisfactory equilibrium statistics for this average by accepting or rejecting trajectories with the weight functional formed from the indicator functional times the equilibrium distribution of trajectories, each trajectory having a duration of . This is done by first equilibrating the system with a trajectory that runs for many structural relaxation times. A few short sections of this trajectory are then chosen to provide first examples of a transition or excitation associated with the tagged particle. From these first trajectories, a series of shooting and shifting moves are then performed Bolhuis et al. 2002; Dellago et al. 2002 generating an ensemble of thousands of independent examples of excitations. The ensemble produced by this Monte Carlo walk through trajectory space is the subspace, with proper statistical weight, of those trajectories in the equilibrium distribution that exhibit excitation dynamics Bolhuis et al. 2002; Dellago et al. 2002. The commitment time is allowed to vary in transition path sampling Bolhuis 2008 so that dynamics determines a distribution of these times without pre-conceived notions of its typical values. In some cases, we have used transition path sampling to sample excitations in very cold systems out of equilibrium. These systems are initialized from configurations of a warmer equilibrated system, but with temperatures chosen from a Maxwell-Boltzmann distribution at a lower temperature.
Considering the range of possible displacement lengths, should be large enough to record significant displacements. In addition, values of should be smaller than those for which typical trajectories would likely visit intermediate states for long periods of time. Larger values of would obscure separations in time scales, and the equilibrium average would fail to exhibit linear growth with respect to Dellago et al. 1999. For the systems that we consider here, taking to be 0.2 to 2 particle diameters proves to be satisfactory, and for this range of displacement lengths there is a range of suitable commitment times where is independent of . At supercooled conditions, the typical mean values and fluctuations of and are of the order of to integration steps, and very much smaller than structural relaxation times.
Transition path sampling Bolhuis et al. 2002; Dellago et al. 2002 is a necessary tool in this study to extract information at very low temperatures. At moderate supercooling conditions, however, where sufficient numbers of excitations are present at any one time frame, much of what we have done can be done with straightforward molecular dynamics. Indeed, this is the procedure we use to evaluate from Eq. 6. But in doing so, it is necessary to identify the duration of a time frame, or . Making that identification with transition path sampling is reasonably easy because the equilibrium distribution for instanton times can then be generated automatically Bolhuis 2008. One can then limit the search for excitations with values not far from most probable values. For a given displacement length , especially for larger than a particle diameter, taking much larger than the most probable value will lead one to harvest sequences of several temporally separated excitations at smaller length scales. With straightforward molecular dynamics, one can be assured of not using too large a value of by checking that this time is not beyond plateau value times, as discussed below.
Many features we highlight in this study of atomic glass forming liquids are also features of granular media Dauchot et al. 2005. In the context of granular media, these features have been studied with order parameters that focus on distinct changes in cages surrounding tagged particles Candelier et al. 2009. Those order parameters might be useful alternatives to those considered in this work.
III Localized excitations and hierarchical dynamics
With the operational definition of excitation dynamics given above, the structure and energetics of this dynamics can be examined. Thus, in this section we establish that the equilibrium statistics of excitation density is that of a dilute gas, and that the energetics of excitations grows logarithmically with the length scale of excitation displacement. This logarithmic growth is the signature of a particular class of hierarchical dynamics, and we show that it predicts the temperature variation of structural relaxation times of a supercooled liquid.
The time durations of the dynamical processes that produce persistent displacements are characterized by the probability distribution of the instanton time, . The behavior of for is plotted in Fig. 2(a) for temperatures ranging from the onset temperature to the lowest accessible temperature in our equilibrium molecular dynamics simulations. The physical meaning of the onset temperature is seen explicitly in Fig. 2(b). As goes below , molecular motions are activated, and the concomitant separation of time scales is apparent from plateau values show in that figure. Beyond the plateau region, the function
will decrease with increasing . The range of the plateau region (i.e., the separation of time scales) increases with lowering . For close to or above , however, motion is not activated and a plateau value does not exist. (The determination of onset temperatures is presented in the next section). The existence of plateau behavior in implies that statistically dominant motions are instantonic Chandler 1978; Hänggi et al. 1990, like that illustrated in Fig. 1.
Transport properties vary by three to four orders of magnitude over the range of temperatures considered in Fig. 3(a), but the distribution of instanton times varies little, and its width is far smaller than structural relaxation times at deeply supercooled temperatures (i.e., ). In this sense, the extent of excitation dynamics in supercooled materials is relatively localized in time. Changes in displacement length, shown in Fig. 2(a), produce significant changes in . These changes relate to the hierarchical nature of the dynamics that we demonstrate in subsequent sections, where we also report on structural relaxation times.
Along with being temporally localized, excitation dynamics are spatially localized. To show this, we consider
for and , and with in the plateau regime. The function is the mean displacement density at for the time frame to given an excitation is at the origin. In the limit of large , , where is the mean particle density and . Deviations from this asymptotic limit reflect the degree with which particle displacements at are correlated to the excitation dynamics at the origin at time 0.
This function has some similarity with the -functions often used to to characterize dynamic heterogeneity Lačević et al. 2003; Berthier et al. 2005; Chandler et al. 2006. Those functions, or susceptibilities, measure mean-square fluctuations of densities of particle-position overlap or displacement after a specified time frame. For example, for the special case of taking and , the function would be similar to the distinct particle contributions to a -function at a time . But the -function is different because it is not limited to that choice of time variables. It is also different because it refers to particle coordinates that are coarse grained over time, while -functions refer to raw coordinates. Further, the -function refers to a displacement density conditioned on an excitation at the origin, and an excitation is distinct from a displacement because an excitation produces a long-lived displacement and not simply a fleeting motion. These features that distinguish the -function from -functions are important to the precision of our analysis.
To use this function as a recorder of excitation size, the relevant time frame should surround time 0 and be of width , so that a pertinent correlation function to consider is . This function oscillates as a function of manifesting local packing of molecules. The same oscillations are found in the radial distribution function, , where is the net density in the time coarse-grained structure of particles at , and the average is taken with an additional particle fixed at the origin. Throughout the supercooled regime, this radial distribution function is essentially temperature-independent, and its correlation length (the coarse-graining length over which its oscillations disappear) is of the order of one particle diameter. Thus, for the purpose of viewing the spatial extent of a single excitation, it is useful to divide out the radial distribution function and consider
The quantity is plotted relative to its value at in Fig. 3(b). The graphs show that, throughout the supercooled regime of the KA model, motions in the same time frame correlate over only a small range of distances – no more than a few atomic diameters – and this correlation range is independent of temperature. This observation holds for all equilibrium state points studied, Fig. 3(b), as well as non-equilibrium state points sampled using transition path sampling, Fig. 3(c). The same holds for all other models we have studied.
The fact that motions in the same time frame correlate over only a small range of distances implies that the underlying configurations, which facilitate motion, must also be small with pair correlation lengths no larger than the motional correlation lengths seen in Panels (a) and (b) of Fig. 3. A large underlying excitation or large correlation length between such excitations would imply large motional correlation lengths, and these are not observed. In contrast, as we will see, much larger length scales that grow with decreasing temperature are associated with structural relaxation and dynamical heterogeneity. In this sense, therefore, elementary excitations are spatially localized. Growing length scales of structural relaxation must arise from correlations between those excitations at different time frames. These dynamical correlations are hierarchical, as we discuss next.
We have computed the concentration of excitations from the second equality of Eq. (6) and we find that obeys the Boltzmann distribution of Eq. (1). Figure 3(d) illustrates this finding. For the length scale considered in that figure, , the plateau time is in the range to . The proportionality constant of Eq. (1) is of the order of , changing slightly from one system to another. The onset temperature, , signals the high temperature boundary for where is characterized by a single energy scale . We will see that this same temperature is the high-temperature end to all signatures of supercooled glassy dynamics. The energy scale depends upon in a hierarchical way, Eq. (2). This finding is illustrated in Fig. 4(a). The specific values of and that summarize these results are given in Table 2.
We provide two different columns of tabulated data for because this quantity can be obtained by either of two ways – by fitting data to Eq. (1) or by fitting data to Eq. (2). It is significant that the two methods give similar values. Those obtained in the former way have the better statistical certainty, so we use those values for predicting relaxation times, which we turn to now.
III.3 Relaxation times can be predicted from excitation energy scales.
The logarithmic growth of with respect to has an important implication with respect to the temperature dependence of transport properties. In particular, and in contrast with excitations simply diffusing as a random walker, the logarithmic growth implies a hierarchical dynamics as imagined by Palmer et al. Palmer et al. 1984 in which excitations on one scale combine to build excitations on larger scales. To see how this implies a specific temperature dependence of transport, consider excitations of displacement length for a time slice with thickness of order . The spatial volume will be occupied by of these excitation displacements, where is of the order of . The rate at which a given excitation will dissipate, , is the rate at which that excitation can connect to neighboring excitations of that same displacement , and the energy to build that connection is its activation energy. Therefore,
To check this prediction we have determined structural times through calculations of the self correlation function
For the three-dimensional models, we define the structural relaxation time, , according to , where is the wave-vector at which the structure factor has its main peak . For the two-dimensional models, this definition is unsatisfactory because for those cases initial cage relaxation makes long before structural relaxation sets in; so for , we use . Figure 4(b), shows the excellent agreement between measured structural relaxation times and theoretical prediction for all supercooled systems studied. The reference structural relaxation time, , is determined by computing a single relaxation time in the moderately supercooled regime and using this value to align the curves. These values, listed in Table 2, are approximations to the relaxation times of their respective liquids at the onset to supercooled behavior.
IV Dynamical facilitation, directionality, and stringlike motion
The hierarchical nature of structural relaxation that we have uncovered in the previous two sections implies that dynamical facilitation Fredrickson and Andersen 1984; Chandler and Garrahan 2010 plays a central role in relaxation phenomena. Dynamical facilitation refers to structural rearrangements or excitations allowing for the birth and death of excitations nearby in space. The superposition of such dynamics leads to the formation of excitations on longer length and time scales, thus resulting in dynamical heterogeneity Garrahan and Chandler 2002. To illustrate these features more explicitly, we first consider computer rendered movies of trajectories, and then turn to quantitative analysis.
An example of a trajectory is shown in Fig. 5 for the 2D-50:50 system. The movie depicts inherent structures, with particles colored according to the length of their displacement vector over a time window . Movies with time coarse-grained structures are similar, but we choose to show those of inherent structures to highlight the underlying physics of a most striking feature in these movies. These features are the ubiquitous low-frequency low-amplitude motions that permeate the system (manifested by aqua-colored particles). These motions are not low frequency harmonic modes, as harmonic motions are absent from the inherent structure. Rather, these motions are soft, anharmonic motions of the disordered system. On time scales that are long compared to the period of these small amplitude soft motions, significant long-lived particle displacements occur (yellow, red particles). These motions are associated with what appear to be defects or excitations in the system. String-like patterns of slight mobility (aqua and green particles) surge outward from, and retract back towards, the initial excitations Chandler and Garrahan 2010.
Eventually, rare, strong surges cause local regions to irreversibly deform on larger length scales and time scales. These processes, whereby ubiquitous surging excitations play a significant role in producing rarer structural relaxations, are an intrinsic feature of East-like models. Such surging was noted in Ref. Widmer-Cooper and Harrowell 2009 but the connection to the East model was overlooked. Activation energies grow as displacement length grows, so that frequencies of displacements decrease as length scales grow. Further larger length scale displacements are built from smaller length scale displacements through facilitation. This is why in the movie of Fig. 5, it is evident that short time local relaxation processes represent a microcosm of relaxation on longer length scales and time scales. As time progresses, larger regions of mobility (red, yellow) grow outward from the initial relaxed regions, giving rise to more collective surges that are characterized by broader strings.
While movies of particle motions are instructive, we obtain quantitative measures of hierarchical facilitation from the calculation of relaxation time presented in the previous section. We can also estimate facilitation volumes,
At low temperatures, the sparsity of excitations makes it difficult to obtain accurate statistics for this facilitation volume. Irregularities in the variations seen in Fig. 6(b) give a sense of the statistical uncertainties in our estimates for this quantity. The data collected is sufficient to reveal the following trends: The facilitation volume initially increases as a function of time , as the initial excitations facilitate the formation of new excitations nearby. reaches a maximum value, , at , and this maximum value grows with decreasing temperature as the distance between the initial excitations increases. Over time scales exceeding , the relaxed regions overlap and the system becomes increasingly dynamically uniform with time, causing to tend towards zero. Figure 6(b) shows that consistently grows with decreasing temperature, thus demonstrating the presence of dynamical facilitation at all supercooled state points.
Dynamical facilitation occurs with directionality. This is a feature of the East model Jäckle and Eisinger 1991 and its counterparts Garrahan and Chandler 2003; Berthier and Garrahan 2005. For the atomistic models studied herein, a quantitative measure of directionality in dynamical facilitation is given by Donati et al. 1998
Displacing particles that lie in the direction of the tagged particle’s displacement vector have , whereas particles that lie in the opposite direction have .
A measure of excess oriented motion is therefore Donati et al. 1998
Here, as before, denotes the average taken with particle 1 exhibiting the excitation at the origin. Figure 6(c) shows that on relatively short time scales, excitations are most likely to occur behind earlier excitations a short distance away. Figure 6(d) shows that this directional preference increases with decreasing temperature. Figures 6(e) and 6(f) show how directional effects dissipate with increasing time.
Earlier ideas about dynamic heterogeneity, specifically the concepts of strings Donati et al. 1998; Donati et al. 1999; Schrøder et al. 2000; Glotzer 2000; Aichele et al. 2003; Gebremichael et al. 2004; Vogel and Glotzer 2004; Keys et al. 2007 and micro-strings Gebremichael et al. 2004, seem to be connected with our results. In particular, these and other earlier works highlighted the presence of correlated stringlike motions in the heterogeneous dynamics. A small subset of these motions involve a few particles that move nearly simultaneously, a process which is called a microstring Gebremichael et al. 2004. Longer strings and then clusters of strings are built from these microstrings, processes that take place over times that are significantly longer than those of a microstring Gebremichael et al. 2004. In the context of this paper, microstrings coincide with dynamics of elementary excitations, and the building of longer strings and clusters of strings coincide with the facilitated hierarchical dynamics of excitations. The exponential distribution of string lengths reported by earlier studies of dynamic heterogeneity in many types of liquids and in granular materials Donati et al. 1998; Donati et al. 1999; Aichele et al. 2003; Gebremichael et al. 2004; Vogel and Glotzer 2004; Keys et al. 2007 may now be understood as a consequence of the ideal gas statistics obeyed by the elementary excitations. Consequently, the cooperative length scale of stringlike motion that grows with decreasing temperature coincides with a growing distance between excitations. It is a plausible connection that is worthy of further study.
V Discussion
Perhaps the foremost criticism is that no viable procedure has existed from which localized excitations could be seen as emergent properties of atomistic dynamics Biroli and Bouchaud 2009; Berthier and Biroli 2011. This paper has established such a procedure.
Another criticism is the recent suggestion that facilitation is not the cause of intermittency in dynamical activity (births of so-called “avalanches”) Candelier et al. 2009; Candelier et al. 2010a, and further that facilitation diminishes as temperature is significantly lowered below (or density is increased above) the onset to supercooled behavior Candelier et al. 2010b. These seemingly contradictory conclusions to our findings are the result of differences in the definition of facilitation. In our usage, facilitated dynamics is where changes in configurations or microstates occur only in the vicinity of excitations. Excitations refer to configurations or microstates. Thus, in facilitated dynamics, the birth or death of excitations occurs only in the presence of neighboring excitations. From this definition, it follows that connected lines of excitations permeate trajectory space in systems dominated by facilitated dynamics Garrahan and Chandler 2002.
Although excitations will connect throughout space-time, changes in states will not necessarily connect in facilitated dynamics. These changes are the so-called “kinks” of kinetically constrained models. Kinks necessarily record the presence of one or more excitations, but excitations can persist for long periods without the appearance of kinks. References Candelier et al. 2010a and Candelier et al. 2010b ascribe a loss of facilitation to growing spatial and temporal separations between changes in state, but these are growing separations between kinks, not excitations. Indeed, the behaviors noted in those papers are predicted from the facilitated dynamics of the East model and its generalizations. Specifically, the motion of particles in the presence of excitation lines is intermittent Jung et al. 2004. Clustered bursts of activity occur at space-time regions where a particle or group of particles intersects excitation lines. Bursts end and disconnect from other bursts when this intersection ends, after which the particle persists in an inactive state for a relatively long period. This behavior is the essential element of decoupling phenomena Hedges et al. 2007, and it becomes more striking as relaxation times grow because the lengths of these quiescent periods grow with increased supercoolingJung et al. 2004; Hedges et al. 2007. Whether this understanding can be enhanced from analysis presented in Refs. Candelier et al. 2010a and Candelier et al. 2010b remains to be seen.