Excitations are localized and relaxation is hierarchical in glass-forming liquids

For several atomistic models of glass formers, at conditions below their glassy dynamics onset temperatures, ${T_\mathrm{o}}$, we use importance sampling of trajectory space to study the structure, statistics and dynamics of excitations responsible for structural relaxation. Excitations are detected in terms of persistent particle displacements of length $a$. At supercooled conditions, for $a$ of the order of or smaller than a particle diameter, we find that excitations are associated with correlated particle motions that are sparse and localized, occupying a volume with an average radius that is temperature independent and no larger than a few particle diameters. We show that the statistics and dynamics of these excitations are facilitated and hierarchical. Excitation energy scales grow logarithmically with $a$. Excitations at one point in space facilitate the birth and death of excitations at neighboring locations, and space-time excitation structures are microcosms of heterogeneous dynamics at larger scales. This nature of dynamics becomes increasingly dominant as temperature $T$ is lowered. We show that slowing of dynamics upon decreasing temperature below $T_\mathrm{o}$ is the result of a decreasing concentration of excitations and concomitant growing hierarchical length scales, and further that the structural relaxation time $\tau$ follows the parabolic law, $\log(\tau / \tau_\mathrm{o}) = J^2(1/T - 1/T_\mathrm{o})^2$, for $T<T_\mathrm{o}$, where $J$, $\tau_\mathrm{o}$ and $T_\mathrm{o}$ can be predicted quantitatively from dynamics at short time scales. Particle motion is facilitated and directional, and we show this becomes more apparent with decreasing $T$. We show that stringlike motion is a natural consequence of facilitated, hierarchical dynamics.


I. INTRODUCTION
The behaviors of supercooled glass-forming liquids manifest complex correlated particle dynamics. Signature behaviors, which appear below an onset crossover temperature, T o , include super-Arrhenius growth of relaxation times with lowering temperature, stretched exponential time-correlation functions, and transport decoupling. Palmer et al. [1] suggest that these behaviors follow from hierarchical dynamics in which excitations on one scale facilitate dynamics of neighboring excitations [2] thereby creating excitations on larger scales. This idea is encoded by a class of dynamical models, so-called kinetically constrained models (KCMs) [3,4].
Results from one such model agree with experimental ob-For each of five atomistic models at several different densities and temperatures, we ask the question: How does an atom move a distance a between relatively long lived neighboring positions? We answer this question by augmenting extensive molecular dynamics simulations with methods of transition path sampling [8]. In particular, molecular reorganization in a supercooled liquid is a arXiv:1107.3628v4 [cond-mat.stat-mech] 19 Oct 2011 rare event, and transition path sampling can harvest unbiased ensembles of trajectories exhibiting such events. By applying this approach, we find several important results.
First, for a smaller than or not much larger than a particle diameter, particle displacements that stick to a new position for a significant period of time are associated with correlated displacements of only a handful of neighboring particles. These persistent changes in particle position, which serve as indicators of what we call "excitations," are closely related to the micro-strings [9] found in earlier computer modeling studies of heterogeneous dynamics in glass forming liquids. At supercooled conditions, i.e., at temperatures T below the onset temperature, T o , their occurrence is sparse. They arise in localized regions of relatively high mobility whose size is largely independent of temperature, and whose spatial distribution is that of a dilute gas.
Second, for a given displacement length a, we find that the equilibrium concentration of excitations, c a , has a Boltzmann temperature dependence where J a grows logarithmically with a, i.e., J a − J a = γ J σ ln (a/a ) .
Here, σ is the space-filling diameter of a particle, the values of γ and J σ 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 J a . In particular, as the average distance between excitations grows with lowering T , the energy scale for relaxation also grows with lowering T . Specifically, from Eq. (2), we argue that the structural relaxation time, τ , obeys with J = γ/d f J σ , where d f is the fractal dimensionality of heterogeneous dynamics. The reference time, τ o , depends little upon temperature, and it is of the order of the time to reorganize a local arrangement of particles in the presence of an excitation, which in turn is of the order of structural relaxation times of the liquid at temperatures above T o . We show that Eq. (3) holds quantitatively with d f /d ≈ 0.9 and 0.8 for physical dimensions d = 2 and 3, respectively.
B. 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 [6,10]. The logarithmic scaling from which it results, Eq. 2, is that of the d = 1 East model [11] and higher dimensional generalizations [5,12]. Here, by providing a microscopic recipe for detecting excitations, computing c a and predicting τ , we validate East-like KCMs as good models for the dynamics of atomistic glass formers.
These technical and quantitative advances seem especially significant in light of the implied physical picture. In particular, we establish herein that the principal growing length that governs the slowing of dynamics in a structural glass former is the mean distance separating excitations, a ∝ c −1/d f a . Relaxation requires the correlated dynamics of neighboring excitations, and the greater their separation, the longer it takes to coordinate their motions. It is a picture that appears significantly different than an often invoked structural perspective originating with Adam and Gibbs [13,14]. That alternative imagines an underlying mosaic of ordered domains, so-called "cooperative rearranging regions." Relaxation follows from reorganizing these regions [15][16][17], and dynamics slows because mosaic domains grow. Our findings do not preclude a picture based upon growing cooperative rearranging regions, but a link between it and the localized excitations that we document would need to relate a to the size of mosaic elements, and whether an operational definition of the latter can accomplish this task is unclear.

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.
A. We study the molecular dynamics of five distinct simple liquid mixtures.
The atomistic models we have simulated are twocomponent mixtures for which particles of type α interact with those of type γ with pair potentials, u αγ (r), where r is the separation of the pair. These potentials are shifted and truncated Lennard-Jones potentials, where u (LJ) (r ; σ, ) = 4 [(σ/r) 12 − (σ/r) 6 ]. For each, we use N = 10, 000 particles in total, with f N and (1 − f )N being the number of A-particles and B-particles, respectively. The different models are distinguished by different choices of length and energy parameters, of mixing fraction f , and of particle masses m A and m B . For notational ease, we use dimensionless quantities throughout, where Boltzmann's constant is unity, and the unit of time is ( /mσ 2 ) −1/2 , where = AA , m = m A , and σ = σ AA are the units of energy, mass, and length, respectively. The first three models detailed in Table I have been developed and studied by others before us. The last two are modifications of models examined by Harrowell and co-workers [21,22]. In particular, we replace the inversetwelve power potential of Refs. [21] and [22] with Weeks-Chandler-Andersen [23] (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 N , V and T fixed, using the HOOMD-Blue simulation code [24] on graphics processors.
B. 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 microstates. We further comment on this point later in this paper.
Non-trivial displacements -the recorders of excitations -are different than fleeting vibrations. Inherent structures [25,26] help distinguish the former from the latter. The inherent structure of a configuration of an N -particle system is obtained by steepest descent to the nearest minimum in the potential energy landscape [27]. 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 ∆t, defined below. Where r i (t) denotes the instantaneous position of the ith particle in the system at time t, we user i (t) to denote its corresponding position in the time-coarse-grained structure, i.e.,r where δt is the coarse-graining time. On the occasions where we refer to inherent structure coordinates, we use r IS i (t) to denote the position of particle i in the inherent structure at time t. 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 ∆t, 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 t particle i is associated with an excitation with displacement length a: Here, θ(x) = 1 or 0 for x ≥ 0 or < 0, respectively. The products are over time steps of a trajectory that extends for a time t a > ∆t. 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, t a is a plateau or commitment time [28,29] -a time for which the equilibrium average of h i (t, t a ; a) grows linearly with t a . It is typically 3 or 4 times the mean instanton time for that displacement length a, ∆t a . As such, the indicator given by Eq. 5 discards trajectories that exhibit displacements with length scale a 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,   [18]. c Wähnstrom model from Ref. [19]. d Weeks-Chandler-Andersen repulsive force mixture from Ref. [20]. e WCA potentials with the same particle-size ratios and mixing rules as the inverse power potentials of Ref. [21]. f WCA potentials with the same particle-size ratios and mixing rules as the inverse power potentials of Ref. [22].
for reasonable choices of a, 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, V is the volume and N is the total number of particles in the system. This second (approximate) equality in Eq. 6, where h i (0, t a ; a) is replaced by θ (|r i (t a + ∆t a ) −r i (0)| − a), holds provided t a 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 h i -functions with the θ-functions when calculating c a at supercooled conditions. This quantity c a 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 h i (t, t a ; a) times the equilibrium distribution of trajectories, each trajectory having a duration of t a . 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 [8,30] 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 excita-tion dynamics [8,30]. The commitment time t a is allowed to vary in transition path sampling [31] 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, a should be large enough to record significant displacements. In addition, values of a should be smaller than those for which typical trajectories would likely visit intermediate states for long periods of time. Larger values of a would obscure separations in time scales, and the equilibrium average h i would fail to exhibit linear growth with respect to t a [32]. For the systems that we consider here, taking a 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 c a is independent of t a . At supercooled conditions, the typical mean values and fluctuations of ∆t and t a are of the order of 10 2 to 10 3 integration steps, and very much smaller than structural relaxation times.
Transition path sampling [8,30] 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 c a from Eq. 6. But in doing so, it is necessary to identify the duration of a time frame, t a or ∆t. Making that identification with transition path sampling is reasonably easy because the equilibrium distribution for instanton times can then be generated automatically [31]. One can then limit the search for excitations with ∆t values not far from most probable values. For a given displacement length a, especially for a larger than a particle diameter, taking ∆t much larger than the most probable value will lead one to harvest sequences of several temporally separated ex-  Table I). The unit of length is σ and the unit of time is ( /mσ 2 ) 1/2 ; i.e, see text. (a) The upper panel shows the displacement for the unprocessed coordinate of a tagged particle at time t together with that displacement in the inherent structure coordinate and in the time coarse-grained structure, i.e., |r1(t) −r1(−∆t/2)|. Here ∆t is the instanton time and ta is the plateau time. The averaging window used is δt = 0.6. The snapshots are drawn to show the space-filling sizes of particles and are colored to indicate the distance of a particle from its position at t = −∆t/2. The arrows indicate the direction of that motion. (b) The lower panels show similarly rendered snapshots for the full 10,000-particle system over one instanton time for a tagged particle at the center of the box. Enlargements show detailed structures of the excitation dynamics, where arrows indicate the direction of the color-labeled displacements. Blue indicates relative immobility while red indicates relative mobility. The depictions are for particle density ρ = 0.75 at four different temperatures (as indicated). The glassy dynamics onset temperature for this system at this density is To ≈ 2.1.
citations at smaller length scales. With straightforward molecular dynamics, one can be assured of not using too large a value of ∆t 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 [33]. In the context of granular media, these features have been studied with order parameters that focus on distinct changes in cages surrounding tagged particles [34]. 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  Table I) at the lowest equilibrated temperature for three different displacement lengths, a = σ/4, σ/2 and σ. (b) The mean value of the rate function R(ta) demonstrating that a plateau value exists for T smaller than the onset temperature, To. The particular system considered here is the W model at the density ρ = 1.296, where To 0.88. Lines in Panel b are drawn through data points as guides to the eye. 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.
A. Dynamic indicators of excitations are localized with spatial and temporal extents that are temperature independent.
The time durations of the dynamical processes that produce persistent displacements are characterized by the probability distribution of the instanton time, P a (∆t). The behavior of P a (∆t) for a = σ 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 T goes below T o , 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 t a . The range of the plateau region (i.e., the separation of time scales) increases with lowering T . For T close to or above T o , 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 R(t a ) implies that statistically dominant motions are instantonic [28,29], 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., T < T o ). 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 P a (∆t). 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 t = −t a /2 and t = t a /2, and with t a in the plateau regime. The function µ(r, t, t ; a) is the mean displacement density at r for the time frame t to t given an excitation is at the origin. In the limit of large r, µ(r, t, t ; a) → ρµ(t − t ), where ρ is the mean particle density and µ(t − t ) = |r i (t) −r i (t )| . Deviations from this asymptotic limit reflect the degree with which particle displacements at r are correlated to the excitation dynamics at the origin at time 0. This function has some similarity with the χ 4 -functions often used to to characterize dynamic heterogeneity [35][36][37]. Those functions, or susceptibilities, measure meansquare fluctuations of densities of particle-position overlap or displacement after a specified time frame. For example, for the special case of taking t = 0 and t − t = t a , the function µ(r, t, t ; a) would be similar to the distinct particle contributions to a χ 4 -function at a time t a . 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 χ 4 -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 χ 4 -functions are important to the precision of our analysis.  Table II. To use this function as a recorder of excitation size, the relevant time frame should surround time 0 and be of width t a , so that a pertinent correlation function to consider is µ(r, −t a /2, t a /2; a) − ρµ(t a ). This function oscillates as a function of r manifesting local packing of molecules. The same oscillations are found in the radial distribution function, ρ(r) 0 , where ρ(r) is the net density in the time coarse-grained structure of particles at r, and the average · · · 0 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 F (r; a) = µ(r, −t a /2, t a /2; a) ρ(r) 0 µ(t a ) − 1 The quantity F (r; a) is plotted relative to its value at r = σ 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.
B. Excitation densities obey Boltzmann statistics with energy scales that grow logarithmically with displacement length.
We have computed the concentration of excitations c a from the second equality of Eq. (6) and we find that c a obeys the Boltzmann distribution of Eq. (1). Figure 3(d) illustrates this finding. For the length scale considered in that figure, a = σ, the plateau time t a is in the range 20 to 30. The proportionality constant of Eq. (1) is of the order of (t a a d ) −1 , changing slightly from one system to another. The onset temperature, T o , signals the high temperature boundary for where c a is characterized by a single energy scale J a . We will see that this same temperature is the high-temperature end to all signatures of supercooled glassy dynamics. The energy scale J a depends upon a in a hierarchical way, Eq. (2). This finding is illustrated in Fig. 4(a). The specific values of J σ and γ that summarize these results are given in Table II.
We provide two different columns of tabulated data for J σ 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.
C. Relaxation times can be predicted from excitation energy scales.
The logarithmic growth of J a with respect to a 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. [1] 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 a for a time slice with thickness of order t a . The spatial volume V will be occupied by V c a /ν of these excitation displacements, where 1/ν is of the order of t a a d . The rate at which a given excitation will dissipate, 1/τ a , is the rate at which that excitation can connect to neighboring excitations of that same displacement a, and the energy to build that connection is its activation energy. Therefore, is the distance to connect a neighboring pair of excitations. Motions in East-like models are nearly linear [4,5], so that the fractal dimensionality, d f , is expected to be close to the physical dimensionality d. It can be less than d to the extent that paths connecting excitations are not linear. Equations (10) and (11) together with (2) yield At a ≈ σ, this relaxation time is the structural relaxation time for the systems we consider. Hence, we find Eq. (3), 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 1/e = F s (q 0 , τ ), where q 0 is the wave-vector at which the structure factor has its main peak q 0 ≈ 2π/σ. For the two-dimensional models, this definition is unsatisfactory because for those cases initial cage relaxation makes F s (q 0 , t) < 1/e long before structural relaxation sets in; so for d = 2, we use 0.1 = F s (q 0 , τ ). 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, τ o , 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 II, are approximations to the relaxation times of their respective liquids at the onset to supercooled behavior.
The fractal dimension d f can be computed directly [38] by clustering neighboring particles that have displaced a distance a = σ over a time window τ . The value of d f is then d f = d ln(n)/d ln(r), where n is the number of particles contained within a spherical volume of radius r positioned at the cluster centroid. We find that this procedure gives d f between 2.5 and 2.7 for physical dimension d = 3, and between 1.8 and 1.9 for d = 2. Any values of d f within those ranges are equally satisfactory for fitting the transport data. The values of d f used in Fig. 4(b) are d f = 2.6 for the three-dimensional systems and d f = 1.8 for the two-dimensional systems, consistent with 0.8d and 0.9d, respectively. The one-dimensional version of the dynamical scaling we find is that of the East model [11], where d f = d = 1. The remarkable data collapse shown in Fig. 4(b) is not the result of simple curve fitting. As J σ , d f , T o and γ are all determined independent of τ , the collapse of data is a demonstration of the dynamical picture we have derived.

IV. DYNAMICAL FACILITATION, DIRECTIONALITY, AND STRINGLIKE MOTION
The hierarchical nature of structural relaxation that we have uncovered in the previous two sections implies  Table II. The inset shows the uncollapsed data. The data set 2D-68:32, ρ = 0.75 is omitted for clarity in the inset, as it scale obscures the other data sets. (b) Structural relaxation times predicted from results shown in Fig. 3d and Fig. 4a, which are compared to measured relaxation times. Parameters from statistics are listed in Table II. Data symbols are the same as those in Fig. 3.  a From fits to Eq 1 shown in Fig. 3d. b The ratio of the maximum to minimum values of ca used in fits to Eq 1. c Error in the fit to Eq 1, given by 1 minus the square of the correlation coefficient. d From fits to Eq 2 shown in Fig. 4a. e Error in the fit to Eq 2, given by 1 minus the square of the correlation coefficient.
that dynamical facilitation [3,7] 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 result-ing in dynamical heterogeneity [39]. To illustrate these features more explicitly, we first consider computer rendered movies of trajectories, and then turn to quantitative analysis.
A. Facilitation and hierarchical dynamics are evident in movies of supercooled liquid dynamics. 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 t. 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 longlived 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 [7].

An example of a trajectory is shown in
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. [40] 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.
B. Dynamical facilitation, directionality and string-like motion are pronounced, and increasingly so as temperature is lowered.
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, v F (t) = µ(r, t a /2, t; a) ρ(r) 0 µ(t − t a /2) − 1 dr.
The denominator dividing into the µ-function in Eq. (13) is the value of the µ-function in the absence of dynamical correlations with the initial excitation at the origin.
We find that the integrand in Eq.13 looks much like the function F (r; a) shown graphically in Fig. 3 b and c, but with a range that varies with time. Thus, v F (t) determines the volume of space where dynamics at time t is correlated to that initial excitation. In the absence of dynamical facilitation, particle displacements are not correlated with the initial excitation at the origin, and v F (t) goes to zero. The qualitative behavior of v F (t) is shown in Fig. 6(a) for the W system at several different supercooled temperatures. The other systems we have studied behave similarly. 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 t, as the initial excitations facilitate the formation of new excitations nearby. v F (t) reaches a maximum value, v max F (T ), at t τ , 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 v F (t) to tend towards zero. Figure 6(b) shows that v max F (T ) 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 d = 1 East model [11] and its d ≥ 2 counterparts [5,12]. For the atomistic models studied herein, a quantitative measure of directionality in dynamical facilitation is given by [41] ω(r, cos θ, t, t ; a) = Displacing particles that lie in the direction of the tagged particle's displacement vector have θ i ≈ 0, whereas particles that lie in the opposite direction have θ i ≈ π.
A measure of excess oriented motion is therefore [41] F (r, cos θ, t, t ; a) = ω(r, cos θ, t, t ; σ) where FIG. 5: (Click to view video clips). Time evolution of particle displacements illustrate hierarchical dynamics in the 2D-50:50 system well below its onset temperature. Specifically, this trajectory runs at a temperature T = 1.1, and the onset temperature is To ≈ 2.0. The movie depicts the displacement of each particle in terms of energy-minimized inherent structure coordinates |r IS i (t) − r IS i (0)|, with blue indicating overlap with initial position and red indicating a displacement of at least on particle diameter. See key. The movie, with frames separated by a reduced time of 25, which is approximately equal to ta for this system, spans a time of about 5 τ , where the structural relaxation time is τ ≈ 500ta. At very early times, non-oscillatory particle motions are sparse and spatially decorrelated, indicating the appearance of the initial excitations. As time progresses, regions of high mobility (red, yellow) grow outward from the initial excitations until eventually, the regions connect on a timescale of about τ . Trajectories for the same system at temperatures above the onset temperature do not exhibit these features of correlated dynamics.
FIG. 6: Dynamics is facilitated, directional, and string-like, and more apparently so as temperature is lowered. (a) Dynamical facilitation volume, vF(t), for a representative W system at several different temperatures. (b) The peak of the facilitation volume as a function of temperature for the 12 atomistic models studied. (c) F (r, cos θ, t, t ; a) for a single temperature with parameters a = σ, t = ta/2, and t = 3ta/2 as a function of both the facilitation direction θ and the distance r from the tagged particle at the origin at time ta/2. (d) F (r, cos θ, t, t ; a) as a function of the temperature with the same parameters as in (c) for a single facilitation angle θ ∼ π. (e) F (r, cos θ, t, t ; a) for a single temperature with parameters a = σ, t = ta/2, and t = τ , where τ is the structural relaxation time, as a function of both the facilitation direction θ and the distance r from the tagged particle at the origin at time ta/2. (f) F (r, cos θ, t, t ; a) as a function of temperature with the same parameters as in (e) for a single facilitation angle θ ∼ π. (g) Mean string length ls versus the concentration of excitations of displacement length σ, cσ at several different temperatures T for the 10 three-dimensional systems studied. Data symbols are the same as those in Fig. 3.
Here, as before, · · · 0 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.
The latter observations are consistent with stringlike motion [9,[41][42][43][44][45][46][47]. On time scales of the order t a , stringlike motion takes the form of short "microstrings," whose character does not change with the degree of supercooling [9]. These motions are essentially synonymous with elementary excitations that we characterize here. On timescales of the order τ , microstrings combine to form longer strings, with a length scale that grows with the distance between the initial excitations. In accordance with Refs. [41], [9] and [47], we construct strings by clustering particles that fully replace the initial position of neighboring particles over a time scale that maximizes the mean string length, L s . Strings constructed in this way have a minimum size of 2 particles; thus we subtract this baseline value from the average when quantifying changes in the mean string length, i.e., l s = L s − 2. Figure 6(g) shows that l s grows proportionally as a function of the mean distance between initial excitations, c For our comparison, we consider excitations associated with particle displacements of length a = σ.
Earlier ideas about dynamic heterogeneity, specifically the concepts of strings [9,41,[43][44][45][47][48][49] and microstrings [9], 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 [9]. 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 [9]. 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 [9,41,44,45,47,48] 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
This paper shows that excitations are localized, with static correlations like those of an ideal gas, and that dynamics in these systems are non-trivial because motions occur in a fashion that is facilitated and hierarchical. Perhaps most importantly, the number of particles involved in a single elementary excitation is independent of temperature. As a result, dynamics slow and dynamical heterogeneity length scales grow with decreasing temperature, not because the fundamental mechanism of particle motion becomes increasingly cooperative, but rather because the distance between excitations increases. These findings would thus seem to conflict with pictures, such as the mosaic picture, based on "cooperatively rearranging regions" (CRRs) in which the regions themselves are considered the fundamental object. Instead, here, we view cooperativity as emerging hierarchically over time through facilitation, with excitations being the fundamental object. As noted in the Introduction, we do not discount the possibility of a yet-to-be discovered link that would show how growing CRR size is related to the mean distances connecting neighboring excitations, However, this would presumably require a reformulation of the definition of a CRR. The connection between string length and mean excitation distance demonstrated here might serve as a starting point [50]. While we can await that development, at this stage we conclude this paper by addressing the most common criticisms of the picture supported by our results.
Perhaps the foremost criticism is that no viable procedure has existed from which localized excitations could be seen as emergent properties of atomistic dynamics [16,17]. This paper has established such a procedure.
Another criticism concerns the measured behavior of heat capacity. Unlike reversible heat capacity in an equilibrium system, this property exhibits significant hysteresis near the glass transition temperature, T g . Heat capacity per molecule at temperatures above the range of hysteresis is higher than that at temperatures below the range of hysteresis by an amount typically larger than Boltzmann's constant. In some theories, including Adam and Gibbs', this change in heat capacity, ∆C, is of central importance because it can be related to an imagined vanishing of configurational entropy at a finite temperature, and this vanishing is supposed to signal loss of ergodicity. Inspection of experimental data discredits perceived correlations between relaxation times and these thermodynamic properties (see Fig. 1 of Ref. [10]). Nevertheless, interest in ∆C persists, and one may wonder whether its typical values are consistent with dynamics in structural glass formers pictured in terms of localized excitations [51].
Because localized excitations emerge from molecular dynamics, it is unlikely that this picture is inconsistent with thermodynamics, at least for the systems studied herein. Still, typical values of ∆C reflect the number of translational degrees of freedom removed upon cooling below the glass transition, and this number seems physically interesting. The removal of translational freedom coincides with the removal of fluctuations in the number of excitations. Therefore, a simple connection between ∆C and excitation concentration is ∆C ∝ c σ (T g ). This proportional relationship was proposed in Ref. [5], where it was shown to be consistent with experimental data. The constant of proportionality must grow with the number of particles correlated with an excitation of unit displacement length. We see from this paper that this number is of order 10 or more. Whether further steps can be taken to predict precise values of the proportionality constant remains to be seen. These steps may involve attempts to understand connections between local vibrational modes and excitations. The former are not causative of dynamics [52], but they are correlated to the latter [53], and the former extends over regions of space that are large compared to the latter [52].
Another criticism is the recent suggestion that facilitation is not the cause of intermittency in dynamical activity (births of so-called "avalanches") [22,34], and further that facilitation diminishes as temperature is significantly lowered below (or density is increased above) the onset to supercooled behavior [54]. 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 [39].
Although excitations will connect throughout spacetime, 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 [22] and [54] 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 [55]. 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 [20], and it becomes more striking as relaxation times grow because the lengths of these quiescent periods grow with increased supercooling [20,55]. Whether this understanding can be enhanced from analysis presented in Refs. [22] and [54] remains to be seen.