Exponential distributions of collective flow-event properties in viscous liquid dynamics

We study the statistics of flow events in the inherent dynamics in supercooled two- and three-dimensional binary Lennard-Jones liquids. Distributions of changes of the collective quantities energy, pressure and shear stress become exponential at low temperatures, as does that of the event"size"$S\equiv\sum {d_i}^2$. We show how the $S$-distribution controls the others, while itself following from exponential tails in the distributions of (1) single particle displacements $d$, involving a Lindemann-like length $d_L$ and (2) the number of active particles (with $d>d_L$).

We study the statistics of flow events in the inherent dynamics in supercooled two-and threedimensional binary Lennard-Jones liquids. Distributions of changes of the collective quantities energy, pressure, and shear stress become exponential at low temperatures, as does that of the event ''size'' S P d 2 i . We show how the S distribution controls the others, while itself following from exponential tails in the distributions of (1) single particle displacements d, involving a Lindemann-like length d L and (2) the number of active particles (with d > d L ). DOI Many complex systems, including turbulent flows [1], plastically deforming crystals [2], and financial markets [3], are characterized by intermittent, stochastic dynamics. The statistics of dynamical events are important, with particular consequences for the probability of extreme events: In Gaussian statistics, values more than 5 standard deviations from the mean account for less than 10 À6 of the distribution. For a symmetric exponential distribution, on the other hand, they account for about one thousandth. Examples of exponential distributions (which may be considered intermediate between the Gaussian and power-law cases) include bursts of protein production from gene expression [4] and even returns in financial time series (at least for not-too-large values) [5].
Rearrangements of molecules known as ''flow events'' are the fundamental processes giving rise to structural relaxation and flow in disordered systems ranging from the equilibrium case of highly viscous liquids [6][7][8] to nonequilibrium systems such as glasses [9], granular materials [10], and foams [11]. For viscous (or ''glassforming'') liquids, it has been accepted since Goldstein's 1969 paper [6] that the dynamics may be understood in terms of a division into fast vibrational motion around a particular energy minimum, and relatively rare transitions between neighboring minima. It is the latter that are responsible for the slow dynamics [7,12]; thus, a detailed understanding of their nature is essential, particularly since many theoretical approaches start by making assumptions about the flow events [13][14][15]. Computer simulations provide the means to go beyond assumptions, and in recent years have provided much insight into the kinds of molecular motions that occur in a supercooled liquid. Most workers have concentrated on particle displacements. But of more direct relevance to experiments are collective properties of the dynamics-changes in quantities such as potential energy E, pressure p, and shear stress s . Knowledge of their distributions, and how those depend on temperature, is sure to be relevant for understanding the slowing down of dynamics in viscous liquids. Vogel et al. [16] studied the relation between particle displacements and the size of energy changes during flow events, and while they reported exponential tails in both, they did not examine their relationship in detail, while Schrøder et al. also reported such tails in flow-event displacements [7]. Exponential distributions of forces in a simulated Lennard-Jones fluid were observed already in 1987 [17], although their relevance to flow-event properties is unclear. It is now becoming possible to study individual particle displacements experimentally, as Schall et al. have done for a colloidal system [18]; as yet the analysis is limited to a few events. The exponential tails recently reported for the van Hove correlation function [19] reflect the cumulative effect of many flow events, whereas this work focuses on statistics on single flow events.
In this Letter we present simulation data from two-and three-dimensional binary Lennard-Jones (BLJ) model fluids, brought to a viscous state by cooling. By studying the so-called ''inherent dynamics'' [7]-the trajectory obtained by mapping configurations to local energy minima-we identify flow events and study their statistics. Our main results are (1) at low temperatures T the distributions of changes of E, p, and s are exponential; (2) the sum of squared displacements S, a geometrical measure of the size of an event, is the controlling quantity; (3) the mean event size decreases with decreasing temperature; (4) the exponential distribution (ED) of S can be traced to the existence of an exponential tail in the distribution of particle displacements during events, characterized by a Lindemann-like length scale, which defines ''active'' particles. The number of active particles is broadly distributed at high T, typically comprising a large fraction of the particles in the system, whereas at low T it is also exponentially distributed, with a mean of a few tens of particles. This crossover coincides with the increasing relevance of minima transitions to the dynamics at lower temperature.
The parameters for the BLJ potential are [where and are the energy and length scales for interactions between large (L) and small (S) particles] LL ¼ 1, All particles have the same mass m ¼ 1. These parameters are identical to those of the BLJ introduced by Kob and Andersen [20]. The potential was truncated using an interpolating poly- The American Physical Society nomial between 2:4 and 2:7 (; 2 fL; Sg). All results reported here are from constant volume simulations with periodic boundary conditions, with N p ¼ 700 (1372) particles in 2D (3D), of which 60% (80%) were of type L, at a density of 1:2 À2 LL (1:2 À3 LL ), using a time step of 0.01 LL ffiffiffiffiffiffiffiffiffiffiffiffiffiffi m= LL p (0.005 at T ! 1:0). The system size was chosen to allow a large range of event sizes. From now on, all quantities will be reported in the ''natural'' units defined by LL , LL , and m. Two-step relaxation, a signature of landscape-influenced dynamics, first appears at T $ 0:50 (2D) and T ¼ 0:60 (3D).
To identify flow events we carried out Stillinger's procedure [21] of quenching configurations to the ''nearest'' local minimum. For T ! 0:5, we quenched every time step. For lower T we quenched every tenth time step, and if a change of minimum was observed, the simulation was ''backtracked'' and quenched at each of the intervening time steps. We detected events using the sum of squared where d i is the magnitude of the displacement of the ith particle between successive inherent structures: A value greater than 10 À3 is sufficient to distinguish a genuine change of minimum from numerical noise. This criterion is consistent with one based on changes in inherent energy or stress. Care was exercised in the minimization process [22]. In the following we place somewhat more emphasis on the 2D data because better statistics were obtained; due to a larger number of particles in 3D (though the linear size is more than a factor of 2 smaller) and the larger number of neighbors per particle, the minimization is more time consuming (by about a factor of 6). Thus ð10-20Þ Â 10 3 events per temperature were obtained in 3D compared to 10 5 in 2D. Figure 1 shows PðÁ N XÞ as a function of jÁ N Xj, where X is E, p, or s and the Á N denotes the change normalized by the rms equilibrium fluctuations [ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðX À hXiÞ 2 i p ] calculated from the time series of inherent states (T ¼ 0:33). The tails are very close to exponential for jÁ N Xj * 0:25. Though the distributions differ at small values, the tails for p and E are very similar, with the same decay length $ 0:16 (possibly a reflection of the strong-pressure correlations recently reported for van der Waals systems [25]), while that for shear is larger, $ 0:28. To quantify how close to Gaussian a distribution is we use a non-Gaussian parameter NGP hðxÞ 4 i=3ðhðxÞ 2 iÞ 2 À 1 for any quantity x. This is zero for a Gaussian distribution and unity for a (symmetric) ED. The inset of the figure compares distributions of shear stress changes at the highest and lowest temperatures, now normalized to unit width. That from T ¼ 1:5 is clearly more Gaussian, albeit with an exponential tail a few standard deviations from the mean.
By symmetry, changes of shear stress are uncorrelated with those of energy and pressure, meaning that none of these quantities can be considered a meaningful measure of the ''size'' of an event. A natural measure of the size of an event is the quantity used to detect them, S ¼ P i d 2 i . The S distributions are shown in Fig. 2(a). Going from high to low temperature, there is a dramatic reduction in the average size, while the shape of the distribution becomes increasingly exponential [26]. Examining the distributions of ÁX for a narrow range of S, we find they are Gaussian even at low T when the full distributions are exponential [ Fig. 2(b)]. For not-too large values of S (S & 10), the variance for these fixed-S distributions rises approximately linearly with S [right-hand panels of Fig. 2(b)]. This shows that for a given S, ÁX is essentially a sum of random contributions, whose number is proportional to S, the size of the event. This relationship is independent of T for Á s (lower right-hand panel), but this holds only for small S for Áp and ÁE (upper right-hand panel). It is not clear why this is. In all cases the variance tends to saturate as S exceeds 20 (50 for E at T ¼ 1:0).
The above can be made mathematically explicit by writing PðÁXÞ ¼ R 1 0 PðÁXjSÞPðSÞdS. At low T if PðSÞ ¼ ð1=hSiÞ expðÀS=hSiÞ and the conditional probability PðÁXjSÞ is a Gaussian with variance X S, then integration gives PðÁXÞ ¼ 1=ð2 X hSiÞ expðÀ2jÁXj= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 X hSi p Þ. Linear fits to the S < 10 data for T ¼ 0:33 give for X the values 0.32, 1:5 Â 10 À3 , and 5:8 Â 10 À3 for E, p, and s , respectively. The decay lengths of the corresponding EDs should be ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi X hSi=2 p . For T ¼ 0:33, hSi ¼ 1:83, we get values 0.54, 0.012, and 0.023. After normalizing by the rms fluctuations as was done in Fig. 1, these become 0.16, 0.15, 0.26, in reasonable agreement with the decay lengths determined in Fig. 1.
To understand the distributions of S we consider the distributions of individual (large) particle displacements d in flow events, shown in Fig. 3. As found in Refs. [7,16], clear exponential tails appear for d larger than about 0.2. The decay lengths vary surprisingly little over the simulated temperature range, 0.08-0.15. Most of the variation is in the number of particles in the tail region (the average . We find it useful to take the length scale $0:1 suggested by the decay length seriously, and define the ''active'' particles as those with d > d L ¼ 0:1 (independent of T for simplicity).
Distributions of N A , the number of active particles, are shown in the right-hand panel of Fig. 3. There is a striking shift in both the shape and the mean value from high to low T: A typical event at high T involves most of the system; indeed, larger events are more probable than smaller ones. At low T the distribution becomes exponential with a mean of order a few tens of particles. In this regime, the events are relatively localized though not necessarily compact (see the example in Fig. 3)-stringlike spatial correlation is present in varying degrees in small and intermediatesized events, while the largest events tend to have a more compact structure. The spatial structure of flow events has been investigated by different authors. While Schrøder et al. found stringlike correlations for transitions between minima, Appignanesi et al. used the distance matrix method to identify transitions between so-called metabasins [27]. Their analysis highlights events involving a large fraction of the system, which was relatively small (150 particles). In fact, their data (Fig. 2 of Ref. [27]) are consistent with an exponential tail with characteristic size of order 10-20 particles; the ''democratic events'' would simply be the largest events in the tail. On the other hand , we observe stringlike features in smaller events, but not so much in larger ones. Thus our results encompass both those of Schrøder et al. and Appignanesi et al. Figure 4 shows data for a 3D system. The basic features are similar to the 2D case, and, in particular, the S distribution controls the others. Also shown in the third panel of Fig. 4 is the distribution of N A for a smaller system, N P ¼ 110, which is clearly cut off by the system size-this could well explain a factor of 4 difference in relaxation time between the two sizes that we observe.  The notation d L is meant to suggest a Lindemann-like interpretation. Several authors have emphasized the importance of such a length scale in the context of mechanical instabilities associated with structural relaxation in liquids [15,28,29], also in experiments [18]. Its appearance in the present context can be interpreted as that events are in a loose sense ''discretized'' in units of d L : each event involves some number of particles each of which is displaced some number of d L units. Exponential distributions may be interpreted as the statement that there is no preferred number of basic units (this is analogous to a simple model for polymer size distributions, where the probability to attach a new monomer does not depend on the current length of the chain [30]).
The crossover to exponential distributions is associated with one to well-defined localized events, whose mean size decreases with T. The decrease indicates increasing material stability. A similar result was found by Fabricius and Stariolo when perturbing equilibrium configurations and comparing the corresponding inherent states [31]. If event sizes decrease with decreasing temperature, this presumably has consequences for relaxation times, because a given change takes more events, although decreasing size (in particular ÁE) also suggests that energy barriers for individual events also decrease (we will study energy barriers in upcoming work). On the other hand, when there is a distribution of barriers, there should be increased tendency for repeated reverse transitions. This could be examined by studying correlations between events, which we have ignored here, but will address in upcoming work. One way to account for forward-backward correlations, proposed by Doliwa and Heuer [32], involves grouping minima into larger structures called metabasins and studying transitions between these. Their analysis indicates growing effective barriers as T decreases, consistent with non-Arrhenius slowing down of dynamics. Our analysis focuses on the more fundamental units of dynamics, individual minima transitions, because these are (in principle) unambiguously defined and it is worth understanding their statistics in detail before attempting to coarse grain. The analysis suggests the somewhat paradoxical result that individual event barriers decrease as T does. This emphasizes even more the role of correlations.
Center for viscous liquid dynamics ''Glass And Time'' is sponsored by The Danish National Research Foundation.