Dynamic phase coexistence in glass–forming liquids
One of the most controversial hypotheses for explaining the heterogeneous dynamics of glasses postulates the temporary coexistence of two phases characterized by a high and by a low diffusivity. In this scenario, two phases with different diffusivities coexist for a time of the order of the relaxation time and mix afterwards. Unfortunately, it is difficult to measure the single-particle diffusivities to test this hypothesis. Indeed, although the non-Gaussian shape of the van-Hove distribution suggests the transient existence of a diffusivity distribution, it is not possible to infer from this quantity whereas two or more dynamical phases coexist. Here we provide the first direct observation of the dynamical coexistence of two phases with different diffusivities, by showing that in the deeply supercooled regime the distribution of the single-particle diffusivities acquires a transient bimodal shape. We relate this distribution to the heterogeneity of the dynamics and to the breakdown of the Stokes-Einstein relation, and we show that the coexistence of two dynamical phases occurs up to a timescale growing faster than the relaxation time on cooling, for some of the considered models. Our work offers a basis for rationalizing the dynamics of supercooled liquids and for relating their structural and dynamical properties.
Glass forming systems have a spatially and temporally heterogeneous dynamics1 as revealed, for instance, by the time evolution of the Van Hove (vH) distribution function. This is the probability distribution that a particle has moved of a distance along a fixed direction at time , and is a Gaussian with variance if particles move with a constant diffusion coefficient . Conversely, in glass formers the vH distribution has a temporary non-Gaussian shape2; 3; 4; 5, that indicates the temporary coexistence of particles with different diffusion coefficients. It has been suggested6; 7; 8; 9; 10; 11 that this dynamical heterogeneity reflects the transient coexistence of two phases with different dynamical features, commonly indicated as the ‘fast’ and as the ‘slow’ phase. However, in equilibrium systems it has not yet been identified a dynamical order parameter with a transient bimodal probability distribution, which would support the existence of two coexisting phases; indeed, up to now a dynamical order parameter with a bimodal distribution has only been identified in structural glasses driven out of equilibrium introducing a field pinning some of the particles9, and thus inducing the two phases, or more complex constraints on the relaxation dynamics10. Because of this, in equilibrium supercooled liquids the ‘fast’ and the ‘slow’ phase are usually empirically defined, for instance by considering as ‘fast’ particles, chosen to be the ones with the largest displacement12. These empirical criteria are used because the vH distribution cannot have a bimodal or multimodal shape allowing for the clear identification of different coexisting phases. Indeed, if phases with different diffusivities coexist, then the vH distribution will be the weighted sum of different Gaussian functions, all centered in , and will thus have a single maximum. This clarifies that, in order to investigate whereas two or more dynamical phases coexist, one should investigate the diffusivity distribution, not the vH distribution. Unfortunately, the investigation of the diffusivity distribution is difficult, and only stationary diffusivity distributions can be obtained via a direct inversion of the vH distribution 13; 14.
Here we report the first measure of the time evolution of the single particle diffusion coefficient, for different model systems: the standard Kob–Andersen Lennard–Jones (KALJ) binary mixture15; 16; 17, a binary mixture of soft-spheres in two dimensions18, and the Kob–Andersen lattice gas model19; 20. In the deeply supercooled regime, we find this distribution to temporarily acquire a bimodal shape, thus proving the transient coexistence of two distinct dynamical phases. In the long–time limit the two phases mix and the diffusivity distribution acquires the expected Gaussian shape, with a variance to mean ratio we show to be related to the breakdown of the Stokes–Einstein relation.
We succeeded in measuring the single particle diffusion coefficient by exploiting the intermittent nature of the single particle motion in structural glasses21; 22; 23; 24; 25; 26; 27. Indeed, particles in a glass spend most of their time confined within the cages formed by their neighbors, seldom hopping to different cages. This allows to describe the dynamics through the continuous time random walk (CTRW) formalism, reviewed in the Appendix. In this framework, the diffusivity of each particle at a given time is proportional to its number of jumps, and the distribution of diffusivities is equivalent to the distribution of the number of jumps per particle. Accordingly, in the following we first show that the CTRW approach quantitatively describes the dynamics of the considered systems, when cages and jumps are identified using a recently developed parameter–free algorithm18. We discuss in detail the KALJ system to show that the CTRW approach quantitatively describes the relaxation dynamics of atomistic systems, not only of kinetic lattice models 28; 29. Then, we use this approach to measure the diffusivity distribution and to investigate its time evolution.
CTRW description of the dynamics
In order to prove that the CTRW approach provides a quantitative description of the dynamics of the KA mixture, we have performed a careful analysis of the single particle cage–jump intermittent motion, for temperatures slightly above the mode–coupling one15; 16; 17, . Fig. 1a–c illustrate the distribution of the persistence time and of the jump length , that fix the temporal and spatial features of the system in the CTRW approach, as well as the distribution of the time particles wait in their cages before making a jump, . No correlations between the persistence time and the jump length have been found, in agreement with the CTRW scenario. Panel d illustrates the decay of the persistence. At short times all jumps contribute to the decay of the persistence; we therefore observe , as is the rate at which particles jump, and . At long times () the persistence is found to decay with a stretched exponential, . This implies as verified in Fig.1a.
The temperature dependence of the main quantities characterizing the cage–jump motion is illustrated in Fig.2. We observe the time scales and to have an Arrhenius and a super–Arrhenius behavior, respectively, and the average squared jump length to decrease on cooling. The temperature dependence of these quantities can be used to rationalize those of the diffusion coefficient and of the structural relaxation time at different wavelength (wavevector ), which are commonly accessed experimentally. Indeed, in the CTRW approach it is easy to verify that . Fig.3a illustrates that this relation is verified at the highest temperatures. Deviations emerge on cooling as subsequent jumps of a same particle becomes spatially correlated, as clarified by the subdiffusive dependence of the mean square displacement versus the number of jumps illustrated in panel b. The relaxation time scales as the average time a particle needs to move a distance . Since in the CTRW approach subsequent jumps of a same particle are spatially uncorrelated, this time is that particles need to perform, on average, jumps, and is fixed by the average time particles wait before making the first jump, , and the subsequent ones, , as well as by the average jump duration, :
The last term is actually negligible at low temperatures, where 18.
Fig.3 shows that this prediction agrees very well with the measured data in the investigated range of , with a coefficient of proportionality of the order of . As in the case of the diffusivities, small deviations are observed at the lowest temperatures. These results clearly demonstrate that and respectively correspond to the and to the relaxation time scales of structural glasses28; 30; 31, and confirm that the breakdown of the Stokes–Einstein (SE) relation, which is the increase of the product on cooling, is mainly due to the increase of the ratio, as in lattice model, but it is also affected by the temperature dependence of the jump length. Indeed, we find the length scale28 below which the breakdown of the SE relation occurs, , that is estimated equating the first two terms of the r.h.s. of Eq.1, to also depend on the spatial features of the jumps.
The quantitative description of the relaxation dynamics through the statistical features of the cage–jump motion allows to exploit the features of the CTRW approach to measure the distribution of the single particle diffusivities. Indeed, within the CTRW the diffusivity of particles that have performed jumps at time is . The diffusivity is therefore simply proportional to the number of jumps per unit time. Fig.s4 illustrates the time evolution of the number of jumps per particle rescaled by the average number of jumps , which coincides with the distribution of the single particle diffusion coefficient normalized by the average diffusion coefficient, . The inset and the main panel show results obtained at a high and at a low temperature, respectively. For , is peaked around zero as most particles have not jumped; conversely, in the infinite time limit the distributions have a Gaussian shape with average value . We observe that, at high temperature, the distribution gradually broadens in time, and its maximum move from to . At low temperature, conversely, the distribution acquires a temporary bimodal shape before reaching the asymptotic Gaussian one. The bimodal shape proves the existence of a time window in which two phases of particles with different mobilities coexist. The two phases emerge because of the existence of two well separated timescales and . Indeed the slow timescale, , controls the value of the peak at , that equals the persistence correlation function, . Conversely, the fast timescale, , controls the average value of the distribution, as the position of the second maximum asymptotically occurs at . We stress that the phases with an high and with a low diffusion coefficient cannot be uniquely associated to particles that have moved over a small or over a large distance, respectively, as the average displacement of each particle is zero. This is why a bimodal distribution is not observed in the vH distribution function.
The time evolution of the distribution of the diffusivities gives further insights into the dynamics of the system. Indeed, Fig.s5a,b show that at long times the variance to mean ratio of reaches a plateau value , that grows on cooling. This plateau value can be related to the ratio of the two timescales and . In fact, within the CTRW framework32; 33 and , where . Given the relation between the persistence time and the waiting time distributions34 (see Appendix), it follows and thus . We have verified this prediction considering, beside the KA model, also a binary mixture of harmonic spheres18 and the kinetically constrained Kob–Andersen three dimensional lattice gas model19; 20, as illustrated in Fig.5(inset). The lattice model confirms our predictions. The molecular dynamics simulations reproduce the asymptotically proportionality between and , even though there are small deviations with respect to the CTRW prediction, suggesting the emergence of correlations between successive waiting times at low temperature.
A physical interpretation of the proportionality between and is obtained considering that a distribution with the long–time features of , i.e. a Gaussian distribution with variance , is obtained by randomly assigning the jumps to the particles, in group of elements. Consistently, at high temperature and corresponds to that obtained by randomly assigning the jump to the particles, i.e. a Poisson distribution. The increase of on cooling indicates that at low temperature one might observe, in the same time interval, some particles to perform jumps, and other particles to perform no jumps at all, which clearly suggests .
The CTRW approach does not make any assumption about the spatial correlations between the jumps of different particles. However, particularly in a facilitation scenario in which the jump of a particle facilitates the jumps of nearby particles, one expects these correlation to exist, and hence the two dynamical phases to be spatially segregated. Previous investigations of the spatio–temporal heterogeneities of structural glasses 1 also suggest that this should be the case. Here we investigate these spatial correlation focusing on two correlation functions, both of them related to a scalar field associated to the number of jumps, . Note that is proportional to the average diffusion coefficient of the particles in the volume element .
First, we consider the spatial correlations between the particles that have not jumped at time ,
This equals the correlation function of the particles that have jumped, and thus of the particles that have moved of a distance greater than the jump length. is therefore close to the commonly investigated four-point correlation function, at a wavelength related to the inverse jump length. Then we focus on the spatial correlations between the number of jumps, which is the spatial correlation of the diffusivity, considering the correlation function
We find both correlations functions to decay exponentially, with correlation and , respectively. Their time dependence is illustrated in Fig. 6, for selected temperatures. Both correlation lengths have a maximum as a function of time, We indicate with and , and with and , the time of occurrence and the value of the maxima of the two correlation lengths. As apparent form Fig. 6, both correlations length are small, as usual in structural glasses, and increases on cooling, being much more temperature dependent than . We characterize the temperature dependence of and investigating their scaling with respect to the average persistence time, . Fig. 7 shows that , in agreement with previous results suggesting that the time of the maximum of the dynamical heterogeneities scale as the relaxation time. Conversely, we approximately find . We note that the relation between and is model dependent, as for instance we observe in the Kob–Andersen lattice gas model. Since controls the diffusivity correlations, we expect it to also control the approach of the diffusivity distribution to its asymptotic Gaussian shape, and thus to be the time scale at which the variance to mean ratio reaches its asymptotic value , as in Fig. 5a. Indeed, the data of Fig. 5a are successfully rescaled when normalized and plotted versus , as in Fig. 7(inset).
The study of the time evolution of the diffusivity distribution and of the correlation between the single particle diffusivities allows to identify a new relaxation timescale, . This grows faster than the persistence correlation time on cooling. The emerging physical scenario is as follows: the relaxation time of the system, as measured from the decay of scattering correlation functions, is essentially determined by , as in Eq.1. However, on this time scale the diffusivities of the particles are spatially correlated, and the two dynamical phases are still coexisting, as the diffusivity distribution has not acquired its asymptotic normal shape. It is only on a time of the order of that all correlations are lost. On this time scale the diffusivity distribution has a Gaussian shape, and the diffusivities are spatially uncorrelated.
Our results show that the dynamics of supercooled liquids is characterized by the temporary coexistence of two phases with different diffusivities, one can reveal by describing the intermittent particle motion within the CTRW approach. The presence of these two phases is related to breakdown of the SE relation, that also fixes the variance to mean ratio of diffusivity distribution in the long time limit. The dynamical phase transition is characterized by a time scale , which is that after which the single particle diffusivities are both temporally and spatially uncorrelated. The temperature dependence of this time scale is model dependent, and we have observed it to scale as in the KALJ mixture. This result indicates that the mean squared displacement grows linearly in time for , while the displacement distribution becomes Gaussian on a larger timescale, for . Accordingly, in between these two time scales the dynamics of the system is Fickian but not Gaussian 13; 35; 12; 5.
The clear identification of different dynamical phases might also allow to clarify the debated existence of correlation between the structural and the dynamical properties of supercooled liquids 11; 36; 37; 38; 39; 40. Indeed, these correlations have been looked for arbitrarily dividing the particles in a slow and in a fast phase, introducing a threshold on the particles’ displacements, and then considering how these phases are related to structural properties, such as Vöronoi volume, local order parameters, local elastic constants, or excess entropy. Our results suggest that the slow and the fast phase should correspond to phases with a high and a small diffusion coefficient, we have shown to be unambiguously identified.
We have performed NVT molecular dynamics simulations41 of a standard Kob–Andersen : (:) binary Lennard–Jones (LJ) mixture15. Particles of species and interact via a LJ potential with energy scale and length scale . Values are set as follow: . Particles have the same mass . , and are our units of energy, length and mass. For each temperature, we have first performed simulations to obtain a smooth mean square displacement, from which we have extracted the Debye–Waller factor as in Ref.42. We have then performed other simulations to investigate the statistical features of the cage–jump motion as in Ref.18: we associate to each particle, at each time , the fluctuations of its position computed over the interval , with ballistic time. The trajectory of each particle is segmented in cages and jumps, considering a particle to exit (enter) a cage as becomes smaller (larger) than . This procedure gives access to the duration of each cage, , and to duration and length of each jump. An analogous study has been performed for a two dimensional mixture of particles interacting via a Harmonic potential18. In the case of the KA lattice kinetically constrained lattice glass model19, each particle movement is considered to be a jump.
Appendix – CTRW
The Continuous Time Random Walk (CTRW) approach describes particle motion in supercooled liquids as a stationary isotropic walk process43. The temporal features of this process are fixed by the distribution of the persistence time , which is the time particles wait before making their first step as measured from an arbitary reference time. is related to the distribution of the time particles spend in their cages through the Feller relation34; 44, . The spatial features are fixed by the distribution of the step size . The walk is assumed to be separable as no correlations between and are considered. The relaxation dynamics is monitored by the persistence correlation function20; 28; 29; 31; 45 , that equals the fraction of particles that has not moved up to time . Accordingly the relaxation time , , scales as ; conversely, the diffusivity scales as the number of steps per unit time, . While the CTRW approach assumes the waiting times of different particles to be uncorrelated, this assumption can be easily relaxed to capture the temporal heterogeneities of the dynamics. Indeed, if temporal correlations involve groups of particles, then the fluctuation of the persistence of a particle systems scales as , while its maximum value scales as .
- 1 Berthier, L., Biroli, G., Bouchaud, J-P., Cipeletti, L. & van Saarloos, W. eds. Dynamical Heterogeneities in Glasses, Colloids, and Granular Media, (Oxford University Press, Oxford, 2011).
- 2 Kob, W., Donati, C., Plimpton, S.J., Poole, P.H. & Glotzer, S.C. Dynamical heterogeneities in a supercooled Lennard-Jones liquid. Phys. Rev. Lett. 79, 2827 (1997).
- 3 Donati, C. et al. Stringlike cooperative motion in a supercooled liquid. Phys. Rev. Lett. 80, 2338 (1998).
- 4 Donati, C., Glotzer S.C., & Poole, P.H.. Growing spatial correlations of particle displacements in a simulated liquid on cooling toward the glass transition. Phys. Rev. Lett. 82, 5064 (1999).
- 5 Chaudhuri, P., Berthier, L., & Kob, W. Universal nature of particle displacements close to glass and jamming transitions. Phys. Rev. Lett. 99, 060604 (2007).
- 6 Schmidt–Rohr, K. & Spiess, H.W. Nature of nonexponential loss of correlation above the glass transition investigated by multidimensional NMR. Phys. Rev. Lett. 66, 3020 (1991).
- 7 Glotzer, S.C. Spatially heterogeneous dynamics in liquids: insights from simulation. J. Non–Cryst. Solids 274, 342 (2000).
- 8 Ediger, M.D. Spatially heterogeneous dynamics in supercooled liquids. Annu. Rev. Phys. Chem. 51, 99 (2000).
- 9 Hedges, L.O., Jack, R.L., Garrahan J.P. & Chandler, D. Dynamic Order-Disorder in Atomistic Models of Structural Glass Formers. Science 323, 1309 (2009).
- 10 Cammarota, C. et al. Phase-Separation Perspective on Dynamic Heterogeneities in Glass-Forming Liquids. Phys. Rev. Lett. 105, 055703 (2010).
- 11 Tanaka, H., Kawasaki, T., Shintani, H., & Watanabe, K. Critical-like behaviour of glass-forming liquids. Nature Materials 9, 324 (2010).
- 12 Weeks, E.R., Crocker, J.C., Levitt, A.C., Schofield, A. & Weitz D.A. Three-Dimensional Direct Imaging of Structural Relaxation Near the Colloidal Glass Transition. Science 287, 627 (2000).
- 13 Wang, B., Kuo, J., Bae, S.C., & Granick, S. When Brownian diffusion is not Gaussian. Nature Materials 11, 481 (2012).
- 14 Sengupta, S., & Karmakar, S. Distribution of diffusion constants and Stokes-Einstein violation in supercooled liquids. J. Chem. Phys. 140, 224505 (2014).
- 15 Kob, W. & Andersen, H.C. Scaling Behavior in the –Relaxation Regime of a Supercooled Lennard-Jones Mixture. Phys. Rev. Lett. 73, 1376 (1994).
- 16 Kob, W. & Andersen, H.C. Testing mode–coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function. Phys. Rev. E 51, 4626 (1995).
- 17 Kob, W. & Andersen, H.C. Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture. II. Intermediate scattering function and dynamic susceptibility. Phys. Rev. E 52, 4134 (1995).
- 18 Pastore, R., Coniglio, A. & Ciamarra, M.P. From cage–jump motion to macroscopic diffusion in supercooled liquids. Soft Matter 10, 5724 (2014).
- 19 Kob, W. & Andersen, H.C., Kinetic lattice–gas model of cage effects in high–density liquids and a test of mode-coupling theory of the ideal-glass transition. Phys. Rev. E 48, 4364 (1993).
- 20 Pastore, R., Ciamarra, M.P., de Candia A. & Coniglio, A. Dynamical Correlation Length and Relaxation Processes in a Glass Former. Phys. Rev. Lett. 107, 065703 (2011).
- 21 Appignanesi, G.A., Fris, J.A.R., Montani, R.A. & Kob, W. Democratic Particle Motion for Metabasin Transitions in Simple Glass Formers. Phys. Rev. Lett. 96, 057801 (2006).
- 22 Vogel, M., Doliwa, B., Heuer, A. & Glotzer, S.C. Particle rearrangements during transitions between local minima of the potential energy landscape of a binary Lennard-Jones liquid. J. Chem. Phys. 120, 4404 (2004).
- 23 Vallee, R.A.L., van der Auweraer, M., Paul, W. & Binder, K. Fluorescence Lifetime of a Single Molecule as an Observable of Meta–Basin Dynamics in Fluids Near the Glass Transition. Phys. Rev. Lett. 97, 217801 (2006).
- 24 Fris, J.A.R., Appignanesi, G.A. & Weeks, E.R. Experimental Verification of Rapid, Sporadic Particle Motions by Direct Imaging of Glassy Colloidal Systems. Phys. Rev. Lett. 10, 065704 (2011).
- 25 Vollmayr-Lee, K. Single particle jumps in a binary Lennard-Jones system below the glass transition. J. Chem. Phys. 121, 4781 (2004).
- 26 Carmi, S., Havlin, S., Song, C., Wang K. & Makse, H.A. An energy landscape network approach to the glass transition. J. Phys. A: Math. Theor. 42, 105101 (2009).
- 27 Pastore, R., Pica Ciamarra, M., Pesce, G. & Sasso, A. Connecting short and long time dynamics in hard–sphere–like colloidal glasses. Soft Matter 11, 622 (2015).
- 28 Berthier, L., Chandler, D. & Garrahan, J.P. Length scale for the onset of Fickian diffusion in supercooled liquids. Europhys. Lett. 69 320 (2005).
- 29 Jung, Y., Garrahan, J.P. & Chandler D., Excitation lines and the breakdown of Stokes–Einstein relations in supercooled liquids. Phys. Rev. E 69, 061205 2004.
- 30 Debenedetti, P.G. & Stillinger, F.H. Supercooled liquids and the glass transition. Nature 410, 259 (2001).
- 31 Hedges, L.O., Maibaum, L., Chandler, D. & Garrahan, J.P. Decoupling of exchange and persistence times in atomistic models of glass formers. J. Chem. Phys. 127, 211101 (2007).
- 32 Feller, W., Fluctuation theory of recurrent events. Trans. Amer. Math. Soc. 67, 98 (1949).
- 33 Feller, W. An Introduction to Probability Theory and its Applications (Wiley, New York, 1971), Vol. 2, 2 ed.
- 34 Lax, M. & Scher, H. Renewal Theory and ac Conductivity in Random Structures. Phys. Rev. Lett. 39, 781 (1977).
- 35 Guan, J., Wang, B. & Granick, S. Even Hard–Sphere Colloidal Suspensions Display Fickian Yet Non-Gaussian Diffusion. ACS nano 8, 331 (2014).
- 36 Widmer-Cooper, A., Harrowell, P. and Fynewever, H., How Reproducible Are Dynamic Heterogeneities in a Supercooled Liquid?, Phys. Rev. Lett. 93, 135701 (2004).
- 37 Widmer–Cooper A. & Harrowell, P. Free volume cannot explain the spatial heterogeneity of Debye–Waller factors in a glass-forming binary alloy. J. Non-Cryst. Solids 352, 5098 (2006).
- 38 Hocky, G.M., Coslovich, D., Ikeda, A., & Reichman, D.R. Correlation of Local Order with Particle Mobility in Supercooled Liquids Is Highly System Dependent Phys. Rev. Lett. 113, 157801 (2014).
- 39 Speck, T., Malins, A., & Royall, C.P. First–Order Phase Transition in a Model Glass Former: Coupling of Local Structure and Dynamics. Phys. Rev. Lett. 109, 195703 (2012).
- 40 Malins, A. Williams, S.R., Eggers, J., & Royall, C.P. Identification of structure in condensed matter with the topological cluster classification, J. Chem. Phys. 139, 234506 (2013).
- 41 Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys. 117, 1 (1995).
- 42 Larini, L., Ottochian, A., De Michele C. & Leporini, D. Universal scaling between structural relaxation and vibrational dynamics in glass-forming liquids and polymers. Nature Physics 4, 42 (2007).
- 43 Montroll, E.W. & Weoss, G.H. J. Math. Phys. 6, 167 (1965).
- 44 Tunaley, J.K.E. Theory of ac Conductivity Based on Random Walks. Phys. Rev. Lett. 33, 1037 (1974).
- 45 Chaudhuri, P., Sastry, S. & Kob, W. Tracking Heterogeneous Dynamics During the Relaxation of a Simple Glass Former. Phys. Rev. Lett. 101, 190601 (2008).
We thank H. Tanaka, A. de Candia, A. Fierro and A. Piscitelli for discussions, and acknowledge financial support from MIUR-FIRB RBFR081IUK, from the SPIN SEED 2014 project Charge separation and charge transport in hybrid solar cells, and from the CNR–NTU joint laboratory Amorphous materials for energy harvesting applications.
R.P., A.C. and M.P.C. conceived the project, R.P. and M.P.C carried out simulations and analyzed the data, R.P., A.C. and M.P.C. wrote the paper.
The authors declare no competing financial interests.