Constrained dynamics of localized excitations causes a non-equilibrium phase transition in an atomistic model of glass formers

Constrained dynamics of localized excitations causes a non-equilibrium phase transition in an atomistic model of glass formers

Thomas Speck Institut für Theoretische Physik II, Heinrich-Heine-Universität, D-40225 Düsseldorf, Germany    David Chandler Department of Chemistry, University of California, Berkeley, California 94720, USA

Dynamical facilitation theory assumes short-ranged dynamical constraints to be the essential feature of supercooled liquids and draws much of its conclusions from the study of kinetically constrained models. While deceptively simple, these models predict the existence of trajectories that maintain a high overlap with their initial state over many structural relaxation times. We use molecular dynamics simulations combined with importance sampling in trajectory space to test this prediction through counting long-lived particle displacements. For observation times longer than the structural relaxation time exponential tails emerge in the probability distribution of this number. Reweighting trajectories towards low mobility corresponds to a phase transition into an inactive phase. While dynamics in these two phases is drastically different structural measures show only slight differences. We discuss the choice of dynamical order parameter and give a possible explanation for the microscopic origin of the effective dynamical constraints.


Since crystallization occurs through nucleation virtually any liquid can be supercooled below its melting temperature. But some liquids never become crystals. Their viscosity increases dramatically and at some point internal relaxation cannot keep up with the cooling anymore and they fall out of equilibrium, reaching a state we call a glass (for reviews see Refs. Angell (1995); Debenedetti and Stillinger (2001)). While in principle protocol-dependent, the temperature range at which the transition occurs is narrow and nearly a material property. The idea that this glass transition is a, or determined by a, thermodynamic transition has influenced the theoretical studies of glasses for decades. However, experimentally determined structure factors of supercooled liquids show little to no change while approaching the glass transition. If not global structure, what is the origin of slow dynamics?

The arguably most striking feature of supercooled liquids is the emergence of dynamical heterogeneity (see Ref. Berthier et al. (2010) for a review) below an onset temperature, i.e., while large regions of the liquid are jammed structural relaxation continues through regions which are less rigid. This phenomenon has been observed directly in, e.g., colloidal glasses Weeks et al. (2000) and granular systems Keys et al. (2007). While simple liquids above the onset temperature are well described by a mean-field theory Hansen and McDonald (2006), dynamical heterogeneity leads to different environments for different particles. This observation forms the foundation for dynamical facilitation theory Garrahan and Chandler (2003), a theory of the glass transition as a dynamical phenomenon, see Ref. Chandler and Garrahan (2010) for a review and references. While the interplay between structure and particle dynamics is complicated, the glass transition is independent of much of these details and dynamics is dominated by effective constraints restricting the accessible regions in space-time. The glass transition is controled by the number of excitations marking locally weak or soft regions able to reorganize, and not by a thermodynamic variable.

Crucial for dynamical facilitation theory is the notion of a mobility field coarse-grained both in time and space. A convenient caricature of mobility fields are spin-like excitations on a lattice Fredrickson and Andersen (1984). The effect of the crowded environment on particle motion is mimicked by kinetic constraints, i.e., for a spin to change its state it must be facilitated by one (or more) neighboring excited spin(s). These dynamical rules suffice to give rise to dynamical heterogeneity and a dramatic slow-down of relaxation in the absence of thermodynamic transitions.

Structural relaxation is the decorrelation of particle positions with their initial positions over time. In a dense liquid particle motion is strongly hindered by surrounding particles leading mainly to vibrational motion, short-lived excursions, and rare, collective, long-lived particle displacements that could be described as “cage breaks” Vollmayr-Lee (2004). One approach to glassy dynamics is to find strategies to predict long-time motion from short-time vibrations probing the local structure Widmer-Cooper and Harrowell (2006); Widmer-Cooper et al. (2009). In this paper we pursue a different route and focus on the long-lived displacements as recorders of excitations in the sense of the simple lattice models. We introduce a binary mobility field but use molecular dynamics to determine its time evolution instead of postulated, and necessarily idealized, dynamics.

The paper is organized as follows: In Sec. I we give a brief reminder on kinetically constrained models. In Sec. II we combine the technique introduced in Ref. Keys et al. (2011) to record excitations in atomistic models of glass formers with the fundamental ideas outlined in Refs. Merolle et al. (2005) and Jack et al. (2006). We thus confirm predictions made previously by the study of kinetically constrained lattice models. In particular, for observation times much longer than the relaxation time we find exponential tails in the distribution of mobile particles and a phase transition upon varying a field that couples to mobility. Such dynamical phase transitions Garrahan et al. (2009) have been studied analytically and numerically for kinetic constrained lattice models Garrahan et al. (2007); Speck and Garrahan (2011), spin-glasses Jack and Garrahan (2010), and have also been found in atomistic glass formers Hedges et al. (2009). In Sec. III we then study some structural and dynamical properties of the inactive phase in more detail before coming to the conclusions in Sec. IV.

Figure 1: (a) Intermediate scattering function Eq. (2) at for two system sizes displaying the two-step decay typical for supercooled liquids. (b) The mean density of excited particles (symbols) as a function of inverse temperature below the onset temperature . The dashed line is the fit to Eq. (1). (c+d) Logarithm of the probability of the intensive order parameter scaled by the particle number for (c) and (d) particles, and different observation times . The dashed lines with slope indicate the exponential tails.

I Kinetically constrained lattice models

We start with a brief reminder on two popular variants of kinetically constrained lattice models Ritort and Sollich (2003): the Frederickson-Andersen (FA) Fredrickson and Andersen (1984) and the East model Jäckle and Eisinger (1991). Both models have a trivial energy function where sets the energy scale and is either 1 or 0 indicating whether site is excited or not. Excitations correspond to regions of the supercooled liquid where particles are unjammed and mobile. The effects of local jamming are incorporated by kinetic constraints: in the FA model a site can change state only if a neighboring site is excited, i.e., ; in the East model the site must be excited. Evolution of the mobility field is assumed to be Markovian. A simple choice for the rates is the following. Annihilation of an excitation occurs with rate 1 defining the time scale. Detailed balance then implies the rate for the creation of an excitation (if the constraint is fulfilled). The equilibrium concentration of excitations is


with respect to the reduced temperature . The temperature marks the onset of dynamical heterogeneity.

As a result of constraints, excitations move in lines. Fluctuations in excitations lead to coalescing and branching of excitation lines. These changes of state, or “kinks”, provide the means to detect excitations through particle motion as discussed in the next section. However, the absence of detectable motion does not signify the absence of excitations which, without other excitations reaching out, are quiescent as a direct consequence of the kinetic constraints. The phase transition that we detail in this paper is a transition between two phases of markedly different concentrations of kinks or fluctuations.

Ii Particle mobility as space-time order parameter

ii.1 Excited particles

We have performed extensive molecular dynamics simulations on the Kob-Anderson binary mixture Kob and Andersen (1994), a popular model for atomistic glass formers Yamamoto and Kob (2000); Berthier and Tarjus (2010) (see appendix A for details). It is composed of 80% large (A) and 20% small (B) particles. Simulations are run at temperature well below the onset temperature of heterogeneous dynamics. For comparison, is an estimate of the glass transition temperature for this system Kob and Andersen (1994); Berthier and Tarjus (2010). We have chosen the higher temperature and therefore only moderately supercooled state point to be able to run millions of trajectories over a couple of structural relaxation times. For that state point, the structural relaxation time is as measured by the decay of the intermediate scattering function (see Fig. 1a)


Particle motion is measured on the length scale corresponding to the peak position of the pair distribution function. Fig. 1a demonstrates that at the structural relaxation for the system sizes considered here shows no finite size effects (see also Ref. Karmakar et al. (2009)).

Kinetically constrained models are minimal models incorporating the crucial ingredient of hindered motion in dense, nearly jammed liquids. What they did not provide so far is a concrete prescription on how to map a set of particle positions onto a binary mobility field . To establish such a mapping we follow Ref. Keys et al. (2011) and use long-lived particle displacements of a given length as recorders of excitations. To this end we define the single particle indicator function


where is the unit step function and if the th particle has moved further than the distance in the time interval and otherwise. In the following we will call particles with mobile, or excited to emphasize their role as recorders of underlying excitations. To distinguish non-trivial particle displacements from mere vibrations we use the inherent structure positions  Stillinger and Weber (1984). The inherent structure of a configuration is obtained by steepest descent to the nearest minimum in the potential energy landscape. Even though there is a hierarchy of motion on different length scales Keys et al. (2011), here we focus on the single length . The commitment time is then chosen to be large enough so that a particle can commit to a new position but small enough so that we do not count multiple jumps.

The instantaneous density of mobile particles is


with mean . In Fig. 1b we demonstrate that the temperature dependence of this mean density indeed follows the prediction for excitations Eq. (1). The fitted onset temperature agrees excellently with previous estimates Elmatad et al. (2009); Keys et al. (2011). The fitted energy is  111The method employed here to determine mobile particles is slightly different from Ref. Keys et al. (2011). As a consequence the energy scales differently with length scale . We have checked that dynamics is still hierarchical, , with reported in Ref. Keys et al. (2011)..

ii.2 Low mobility tails

Our objects of interest are trajectories of fixed length . To quantify the amount of mobility in a trajectory we count excited particles through the order parameter


with equally spaced and observation time . In Fig. 1 we show distributions for different system sizes and different observation times obtained through umbrella sampling combined with replica exchange (see appendix B). For trajectories in which motion decorrelates on a time scale much shorter than we would sum over many independent events. Following the central limit theorem the probability distribution then approaches a Gaussian. Indeed, for moderate observation times larger than the relaxation time but below a cross-over time we find such Gaussian distributions (see in Fig. 1). Increasing we observe two effects: for small exponential tails emerge and the shape of becomes non-concave. The physical picture is that highly constrained dynamics facilitates the creation of extended immobile regions, i.e., compared to uncorrelated dynamics it is easier to “remove” mobility from a given space-time volume.

The explanation why the tails of are exponential is as follows Merolle et al. (2005). Assuming that excitations are non-interacting the probability to find at an immobile region of particles is proportional to with a geometric factor independent of and . Combining this probability with for the case that this “bubble” persists leads to plus an offset. In Ref. Keys et al. (2011), evidence is presented that the assumption of non-interacting excitations is indeed a good approximation. Of course, not all bubbles span the entire trajectory connecting the initial with the final state. The temporal extent of a typical bubble grows proportionally to the mean persistence time, i.e., the mean time a particle remains at its initial inherent structure position.

The order parameter Eq. (5) is purely dynamical. Since in a crystal particle mobility is low and motion restricted to defects such an order parameter cannot discriminate between a low activity amorphous and a low activity crystalline phase. Therefore, we also monitor the structure through the orientational order parameter defined in Eq. (15). To allow local fluctuations in structure but prevent global long-range order we use the value of


averaged over the large particles of the final configuration of trajectories. We reject all trajectories with . (see also Fig. 5c below)

ii.3 Active–inactive phase transition

So far we have established that trajectories contributing to the exponential tails in Fig. 1 are those which remember their initial conditions and do not relax within the observation time . At least formally and in computer simulations we can apply a bias in the space of trajectories to stabilize these low mobility trajectories. The fact that the distributions are non-concave implies a phase transition between an active phase corresponding to the liquid in which motion is plentiful, and an inactive phase of low particle mobility. To provide a link with traditional thermodynamics, and the Ising model in particular, imagine the to be spins on a lattice with the order parameter taking the role of the magnetization Chandler (1987). Below the critical temperature the system undergoes a first-order transition between a disordered phase of low magnetization and an ordered phase of high magnetization through applying a field (which we will call in the following). While the statistical treatment is analogous, the underlying physics of a supercooled liquid is of course different from the Ising model. In our case the lattice extends over space and time, and the interactions between “spins” is due to short-ranged forces, geometrical confinement, and the thereof resulting temporal correlations of particle motion.

Figure 2: Mean fraction of excited particles (top) and susceptibility (bottom) vs. the biasing field for selected observation times (from right to left: ) and system sizes (solid lines) and (dashed lines). For clarity we only show error bars for the peak values of the susceptibility . Upper inset: Dependence of the coexistence field on trajectory length and system size . The dashed lines are fits to . The fit parameters are and . From these results we estimate the coexistence field to be (arrow) in the limit . Lower inset: System size dependence of .

In Fig. 2 we plot the mean fraction of excited particles


and the susceptibility vs. the biasing field . The plot shows that the density of mobile particles drops from at to about for . For small we can expand the mean with . The linear behavior around in Fig. 2, therefore, reflects the Gaussian nature of the liquid phase. The coexistence field is obtained from the peak position of maximizing the fluctuations of the order parameter. Increasing either the number of particles or the observation time sharpens the transition. However, space and time are not symmetric. At least in part, this asymmetry reflects that we employ periodic boundary conditions in space whereas trajectories can have quite distinct initial and final states. This leads to temporal boundary effects enhancing the mobility at the beginning and the end of the trajectory Garrahan et al. (2009). For fixed one expects to leading order as shown in the upper inset of Fig. 2. From the fits we estimate the limiting coexistence field to be small but nonzero. Finally, in the lower inset of Fig. 2 we demonstrate the finite-size scaling of the peak values plotted vs. the space-time volume . The dashed line corresponds to , which suggests a first-order transition with a diverging susceptibility and a discontinuous jump of at in the limit of large and/or .

ii.4 Choice of order parameter

Figure 3: Mean intensive activity vs. the biasing field for and (solid line). For comparision, the mean fraction of excited particles (right axis) is shown for both the ensembles defined through (, dashed line) and (, dash-dotted line). Inset: Probability distribution of . The dashed line indicates the exponential tail.

To check whether the transition into an inactive phase is robust with respect to the way we measure mobility in trajectories, we have also considered the dynamical order parameter


This order parameter was used in Ref. Hedges et al. (2009) to demonstrate for the first time a transition between a high activity and a low activity phase in an atomistic model. It sums over the short-time mean-square displacements of particles, which obscures the separation of vibrations from reorganization events that lead to structural relaxation. Nevertheless, as demonstrated in Fig. 3, an abrupt transition from high to low activity can still be observed. To be specific, we define a new ensemble


where we denote the biasing field coupling to by and is any observable. As an illustration for and the mean is plotted in Fig. 3. It resembles the curves shown in Fig. 2 and drops abruptly around . Moreover, we find exponential tails in the probability distribution of plotted in the inset of Fig. 3. In addition to we determine the density of mobile particles , which closely follows the former. However, while the activity drops by a factor of less than two, the number of excited particles is reduced by more than one order of magnitude. The comparision with the curve obtained in the -ensemble shows that the mean fraction of mobile particles approaches the same value (within uncertainties) in both ensembles. We, therefore, conclude that both measures and prepare the same inactive phase through applying a biasing field in trajectory space. The transition is more pronounced in since vibrational motion, which does not cease in the inactive phase evolving at the same temperature as the active phase, contributes significantly to Eq. (8).

Iii Properties of the inactive phase

iii.1 Nucleation of activity

An important consequence of first order phase transitions is nucleation: a system crossing the transition line, although the new phase is favored, has to pay a penalty for interfaces and one has to wait for a large enough nucleus to appear spontaneously. Translated to the present case one might ask what happens if at time for we turn off the field . In the picture of facilitated dynamics a new excitation can only appear close to an existing excitation. Since both density of excitations and kinks, as recorded through , are drastically reduced in the inactive phase we have to wait a certain time before excitations “percolate” through the system and it returns to the liquid phase. The nucleation of activity, therefore, is conceptually different from, say, the nucleation of a crystal which is determined by a single large barrier.

Figure 4: Melting of the inactive phase: (a) Fraction of mobile particles vs. time for an initial configuration prepared with the -ensemble at . The systems remains inactive up to time (arrow), when it suddenly unjams. The dashed line indicates the average density . (b) Another melting trajectory showing a more gradual “thawing”. (c) The distribution of for 10,000 trajectories starting out of the same initial configuration. (d) A different initial configuration showing a much broader distribution of waiting times.

To demonstrate this “melting” of jammed configurations we prepare trajectories at . From these trajectories we then take a single configuration, randomize velocities, and run 10,000 unbiased trajectories out of this initial configuration (see also the isoconfigurational ensemble Widmer-Cooper et al. (2004)). In Fig. 4a a single trajectory is shown. Clearly, the system remains inactive for many structural relaxation times (in this example ) as measured by and then, suddenly, becomes active again with large fluctuations of . In Fig. 4b a different trajectory out of the same initial configuration shows a more gradual transition from inactive to active. In Fig. 4c and d we show distributions of waiting times for two different initial states. These distributions are clearly non-exponential, which is consistent with a variable step process as excitations reach out and reconnect. A detailed comparision of these distributions to predictions from kinetically constrained models is left for a future study.

iii.2 Local structure

Figure 5: Different measures of the structure in the liquid () and the inactive phase () for and : (a) Radial distribution function for the large A particles and (b) for the small B particles. The arrow marks the peak corresponding to a B-B bond with 3 common A neighbors, see sketch and main text. (c) Distribution of the structural order parameter measuring long range order. A particle in a perfect crystal would have . (d) Scatter plot of potential energy (PE) per particle versus the concentration of excited particles for both actual positions and inherent states. Red points are from the ensemble of active states and black points are from the ensemble of inactive states. The harmonic contribution to the potential energy is indicated.

In the introduction we have emphasized that global structural differences between liquid and glass are at most minuscule. However, the fact that the system remains jammed for particle configurations taken from trajectories prepared at indicates that there is a structural difference between these configurations and configurations typically visited in the liquid phase. To make this more quantitative we sample trajectories at fixed and compare the structures as measured by three different methods, see Fig. 5. The pair distribution function for the large (A) particles shown in Fig. 5a demonstrates that liquid and inactive phase are globally indistinguishable. Small differences are seen for the small (B) particles in Fig. 5b. Beyond the simple two-point functions we also consider the histogram of the bond-order parameter as defined in appendix D plotted in Fig. 5c. This order parameter is a convenient measure for long-range order. For every particle it quantifies its local order with for a particle in a perfect crystal. All measures clearly show that the inactive phase is amorphous.

In Fig. 5d we show that the potential energy per particle and the density of mobile particles are uncorrelated in both phases. While in the inactive phase the potential energy of particles is typically lower, the mean difference is much less than the vibrational contribution separating real space potential energies from the inherent state energies. Moreover, as demonstrated in Fig. 5d, there is still an overlap of potential energies between both phases. Hence, we conclude that particles are not trapped energetically but rather due to geometrical constraints.

Differences in structure are picked up by the pair distribution function for the small (B) particles plotted in Fig. 5b. It has been shown that the peaks of for binary mixtures can be assigned to certain local structures: two bonded B particles sharing common A neighbors Fernández and Harrowell (2004). Of particular interest is the second peak corresponding to since it indicates icosahedral coordination shells. Fig. 5b shows that in the inactive phase this local structure occurs more often compared to . This is consistent with recent observations of short-ranged structures in supercooled binary mixtures Coslovich and Pastore (2007); Pedersen et al. (2010). Slow relaxation is attributed to reorganization of particles bound in these structures. Moreover, the drop of in the potential energy of inherent states agrees quantitatively with the drop associated to the formation of these structures Pedersen et al. (2010) (albeit for a slightly different model).

iii.3 Dynamical facilitation

We finally study the behavior of facilitation when going from the active to the inactive phase. First, we note that the fraction of excited particles as plotted in Fig. 6 is independent of temperature in the inactive phase. This indicates that the dynamics in the inactive phase is decoupled from the externally fixed temperature. Second, we study the degree to which particle motion is facilitated. Different methods have been reported in the literature including a mobility transfer function Vogel and Glotzer (2004) and the facilitation volume Keys et al. (2011). In the spirit of a mobility transfer we consider the set of newly excited particles for which the binary indicator function


is . We follow a single particle along a trajectory and through


we count the number of excited particles that have been created in a sphere with radius around the tagged particle under the condition that the tagged particle itself had been excited in the preceding time slice. We define the transfer function


The ratio is plotted in the inset of Fig. 6 using roughly corresponding to the first coordination shell. It shows that the probability that a particle becomes excited close to an already mobile particle increases in the inactive phase.

Figure 6: Mean fraction of excited particles at three different temperatures . Inset: The normalized transfer function Eq. (12) showing an increase of facilitation at large . (all data for and )

Putting all observations together the following picture of dynamics in the inactive phase emerges: The persistence time exceeds the observation time and most particles maintain a high overlap with their initial position. However, some activity continues in isolated regions. The fraction of these mobile particles is decoupled from the temperature. Particles do not become mobile (or immobile) at random but are facilitated through existing mobile particles in their vicinity. Turning off the -field these remaining mobile particles are the seeds from which excitations can reconnect before the system returns to its fluid state.

Iv Conclusions

The separation between fast inter-basin vibrations and slow, activated transitions between inherent structures (or meta-basins Heuer (2008)) is the essence of the energy landscape paradigm Goldstein (1969); Debenedetti and Stillinger (2001). It implies a time evolution that is dominated by rare thermal fluctuations that carry the system from one minimum over a barrier into a neighboring minium Fris et al. (2011). However, there is mounting evidence that such a mechanism competes with, or is even shadowed by, relaxation that occurs through channels that present only low energetic barriers but which are rare and found through “surging” particle motion: examples are string-like motion Gebremichael et al. (2004) and participation maps of low-frequency modes Widmer-Cooper et al. (2004); Chen et al. (2011). Dynamical facilitation theory postulates “excitations” to be the fundamental objects describing such structural weaknesses, and facilitation to be the dominant mechanism at low concentrations of excitations.

The energy landscape picture assumes that the way the landscape is sampled is governed by temperature. From the evidence presented here, we arrive at a seemingly different perspective: the concentration of excitations determines relaxation. At constant temperature the system can be forced into an inactive glassy phase through either removing excitations or suppressing fluctuations leading to “frozen” excitations that are quiescent. In this inactive phase particles vibrate around local energy minima, the statistics of which is consistent with the externally fixed temperature. In contrast, transitions between local minima, or inherent states, are rare and decoupled from this temperature. It appears that these jammed states, in which excitations are arrested, can be created rather easily through local particle rearrangements that do not affect the global structure. In Ref. Jack et al. (2011) it is shown that these states are mechanically more stable than fluid states at a lower temperature. Here we have demonstrated that the melting of jammed states is not consistent with a single crossing of a large free energy barrier but that it is rather a multi-step process. We attribute this multi-step process to the “unfreezing” of excitations, an interpretation that is supported by an increased degree of facilitation in these jammed states. These observations, together with the emergence of exponential tails equivalent to those observed in kinetically constrained models, leads us to the conclusion that the transition is indeed caused by local dynamic constraints. The precise pathways and microscopic mechanisms of the particle rearrangements underlying the active-inactive transition are left for future studies.

As a final note we emphasize that the active to inactive transition is reminiscent of the glass transition. The fundamental difference is that the transition demonstrated in this paper is a transition controlled by a field coupling to a dynamical observable, while the experimental glass transition is controlled by rate of temperature decrease. The connection between the two remains to be quantified.

We have profited from discussions with A.S. Keys and U.R. Pedersen. We thank L.O. Hedges for assistence with early stages of the code and J.D. Chodera for help in understanding and implementing MBAR. TS was supported in part by the Alexander-von-Humboldt foundation, and TS and DC were supported in part by the Director, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division and Chemical Sciences, Geosciences, and Biosciences Division of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix A Simulation details

We have performed extensive molecular dynamics simulations on the binary mixture Kob and Andersen (1994). It is composed of 80% large (A) and 20% small (B) particles possessing the same mass . Particles interact through the continuous truncated and shifted Lennard-Jones potentials

where . The parameters read , , , , , and . Simulations are run at constant volume , constant temperature , and constant number density with the number of particles at positions . Throughout the paper we employ reduced Lennard-Jones units with respect to the large particles, i.e., we measure length in units of , energy in units of , time in units of , and we set Boltzmann’s constant to unity.

Newton’s equations of motion are integrated through the velocity Verlet algorithm with time step 0.005 using LAMMPS Plimpton (1995). For energy minimization we employ the FIRE algorithm Bitzek et al. (2006).

Appendix B Importance sampling

Just as standard Monte Carlo simulations Frenkel and Smit (2002) do importance sampling of configurations, the methods we employ do importance sampling of trajectories. Specifically, we harvest new trajectories using moves from transition path sampling Bolhuis et al. (2002); Dellago et al. (2002). These moves preserve the equilibrium weight of trajectories. Hence, accepting or rejecting a move according to the usual Metropolis criterion

generates an ensemble of trajectories with weight

Here, is a dynamical order parameter that calculates a real number from trajectory , see Eqs. (5) and (8); and is the weight function.

For a single trajectory we store configurations at times with . We employ the “massive stochastic collision” thermostat, i.e., all velocities are randomized at these times and the center-of-mass velocity is subtracted. The use of a stochastic thermostat allows us to perform transition path sampling with so-called ’half moves’ as described in detail in Ref. Dellago et al. (2002). In order to efficiently sample probability distributions of the order parameter we use the quadratic form


To speed up the sampling of trajectory space and decrease the correlations of subsequently generated trajectories we use replica exchange between replicas with different values of . Hence, a single cycle consists of two consecutive steps: (i) every replica generates a new trajectory which is either accepted (and replaces the previous trajectory) or rejected, and (ii) trajectories are swapped between replicas. To obtain good mixing we attempt swaps between all of the replicas, not only neighbors. Data has been acquired from two independent runs, i.e., two independent seed trajectories (except , which is from one run). For a single run we let the trajectories relax for cycles and then recorded trajectories. Table 1 shows an overview for the data gathered to produce Figs. 1 and 2.

system size trajectory length total per replica
2000 50000
2000 40000
2000 40000
2000 50000
2000 50000
2000 50000
2000 50000
3000 60000
6000 60000
222Low activity umbrellas might not have fully equilibrated. 10000 30000
Table 1: System sizes studied and number of harvested trajectories.

Appendix C Mbar

To calculate distributions and expectation values from raw data we use Shirts’ and Chodera’s multistate Bennett acceptance ratio (MBAR) method Shirts and Chodera (2008) and its extension to path ensembles Minh and Chodera (2009). For completeness we briefly summarize the method. We solve

self-consistently for the set of “free energies” . Here, is the weight function corresponding to replica , and is the value of the order parameter along the th trajectory of replica . While in principle MBAR provides an error estimate it requires independent samples, whereas consecutive trajectories obtained using the method described in the previous section are highly correlated. The statistical inefficiency of replica is , where is the correlation time (in unit samples) of some representative observable (here we use ). One possibility to obtain independent samples is to subsample the data series with stride . Here we use all data but weigh replicas according to their relative statistical inefficiencies. Errors are estimated by splitting the data into chunks of trajectories and calculating the standard error of expectation values.

Expectation values of an observable are calculated through

Here, and can correspond to one of the replicas, i.e., and . The advantage of MBAR is that we can employ an in principle arbitrary weight function (given sufficient statistical weight of the sampled data) with


Probabilities are calculated as the expectation value of an indicator function which is 1 if the value falls into bin and 0 otherwise. In particular, the distributions shown in Fig. 1 correspond to the unbiased ensemble with .

The curves shown in Fig. 2 for the mean value are obtained through using the weight function in Eq. (14). In order to sample trajectories at fixed we use this weight function instead of Eq. (13) for a chain of replicas with different values ranging from to . To obtain a set of independent trajectories we keep only every 1000th trajectory for analysis.

Appendix D Orientational order

To quantify orientational order we follow Ref. Steinhardt et al. (1983). For each particle a complex vector

is defined, where are spherical harmonics and the angles and describe the orientation of the displacement vector between particles and with respect to a fixed reference frame. The sum is over all neighbors in the first coordination shell with radius 1.42. The normalized scalar product of the -vectors is

The average over neighbors


is the bond order parameter. It is for a particle in a perfect crystal and acquires a broad distribution with mean for particles in a disordered environment.


  • Angell (1995) C. A. Angell, Science 267, 1924 (1995).
  • Debenedetti and Stillinger (2001) P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2001).
  • Berthier et al. (2010) L. Berthier, G. Biroli, J.-P. Bouchaud, and R. L. Jack, arXiv:1009.4765 (2010).
  • Weeks et al. (2000) E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, and D. A. Weitz, Science 287, 627 (2000).
  • Keys et al. (2007) A. S. Keys, A. R. Abate, S. C. Glotzer, and D. J. Durian, Nature Physics 3, 260 (2007).
  • Hansen and McDonald (2006) J. Hansen and I. McDonald, Theory of Simple Liquids (Academic Press, Amsterdam, 2006), 3rd ed.
  • Garrahan and Chandler (2003) J. P. Garrahan and D. Chandler, Proc. Natl. Acad. Sci. U.S.A. 100, 9710 (2003).
  • Chandler and Garrahan (2010) D. Chandler and J. P. Garrahan, Annu. Rev. Phys. Chem. 61, 191 (2010).
  • Fredrickson and Andersen (1984) G. H. Fredrickson and H. C. Andersen, Phys. Rev. Lett. 53, 1244 (1984).
  • Vollmayr-Lee (2004) K. Vollmayr-Lee, J. Chem. Phys. 121, 4781 (2004).
  • Widmer-Cooper and Harrowell (2006) A. Widmer-Cooper and P. Harrowell, Phys. Rev. Lett. 96, 185701 (2006).
  • Widmer-Cooper et al. (2009) A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, J. Chem. Phys. 131, 194508 (2009).
  • Keys et al. (2011) A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Phys. Rev. X 1, 021013 (2011).
  • Merolle et al. (2005) M. Merolle, J. P. Garrahan, and D. Chandler, Proc. Natl. Acad. Sci. U.S.A. 102, 10837 (2005).
  • Jack et al. (2006) R. L. Jack, J. P. Garrahan, and D. Chandler, J. Chem. Phys. 125, 184509 (2006).
  • Garrahan et al. (2009) J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, and F. van Wijland, J. Phys. A 42, 075007 (2009).
  • Garrahan et al. (2007) J. P. Garrahan, L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, and F. van Wijland, Phys. Rev. Lett. 98, 195702 (2007).
  • Speck and Garrahan (2011) T. Speck and J. Garrahan, Eur. Phys. J. B 79, 1 (2011).
  • Jack and Garrahan (2010) R. L. Jack and J. P. Garrahan, Phys. Rev. E 81, 011111 (2010).
  • Hedges et al. (2009) L. O. Hedges, R. L. Jack, J. P. Garrahan, and D. Chandler, Science 323, 1309 (2009).
  • Ritort and Sollich (2003) F. Ritort and P. Sollich, Adv. Phys. 52, 219 (2003).
  • Jäckle and Eisinger (1991) J. Jäckle and S. Eisinger, Z. Phys. B 84, 115 (1991).
  • Kob and Andersen (1994) W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376 (1994).
  • Yamamoto and Kob (2000) R. Yamamoto and W. Kob, Phys. Rev. E 61, 5473 (2000).
  • Berthier and Tarjus (2010) L. Berthier and G. Tarjus, Phys. Rev. E 82, 031502 (2010).
  • Karmakar et al. (2009) S. Karmakar, C. Dasgupta, and S. Sastry, Proc. Natl. Acad. Sci. U.S.A. 106, 3675 (2009).
  • Stillinger and Weber (1984) F. Stillinger and T. Weber, Science 225, 983 (1984).
  • Elmatad et al. (2009) Y. S. Elmatad, D. Chandler, and J. P. Garrahan, J. Phys. Chem. B 113, 5563–5567 (2009).
  • Chandler (1987) D. Chandler, Introduction to Modern Statistical Mechanics (Oxford University Press, Oxford, 1987).
  • Widmer-Cooper et al. (2004) A. Widmer-Cooper, P. Harrowell, and H. Fynewever, Phys. Rev. Lett. 93, 135701 (2004).
  • Fernández and Harrowell (2004) J. R. Fernández and P. Harrowell, J. Phys. Chem. B 108, 6850 (2004).
  • Coslovich and Pastore (2007) D. Coslovich and G. Pastore, J. Chem. Phys. 127, 124504 (2007).
  • Pedersen et al. (2010) U. R. Pedersen, T. B. Schroder, J. C. Dyre, and P. Harrowell, Phys. Rev. Lett. 104, 105701 (2010).
  • Vogel and Glotzer (2004) M. Vogel and S. C. Glotzer, Phys. Rev. Lett. 92, 255901 (2004).
  • Heuer (2008) A. Heuer, J. Phys.: Condens. Matter 20, 373101 (2008).
  • Goldstein (1969) M. Goldstein, J. Chem. Phys. 51, 3728 (1969).
  • Fris et al. (2011) J. A. R. Fris, G. A. Appignanesi, and E. R. Weeks, Phys. Rev. Lett. 107, 065704 (2011).
  • Gebremichael et al. (2004) Y. Gebremichael, M. Vogel, and S. C. Glotzer, J. Chem. Phys. 120, 4415 (2004).
  • Chen et al. (2011) K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, Z. Zhang, A. J. Liu, and A. G. Yodh, Phys. Rev. Lett. 107, 108301 (2011).
  • Jack et al. (2011) R. L. Jack, L. O. Hedges, J. P. Garrahan, and D. Chandler, Phys. Rev. Lett. 107, 275702 (2011).
  • Plimpton (1995) S. Plimpton, J. Comp. Phys. 117, 1 (1995), available at
  • Bitzek et al. (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Phys. Rev. Lett. 97, 170201 (2006).
  • Frenkel and Smit (2002) D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications (Academic Press, San Diego, 2002), 2nd ed.
  • Bolhuis et al. (2002) P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Annu. Rev. Phys. Chem. 53, 291 (2002).
  • Dellago et al. (2002) C. Dellago, P. G. Bolhuis, and P. L. Geissler, Adv. Chem. Phys. 123, 1 (2002).
  • Shirts and Chodera (2008) M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008).
  • Minh and Chodera (2009) D. D. L. Minh and J. D. Chodera, J. Chem. Phys. 131, 134110 (2009).
  • Steinhardt et al. (1983) P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description