A Case for Dynamic Percolation Underlying Mechanistic Crossovers in the Relaxation of Liquids

TL;DR

This paper explains how dynamic percolation in liquids causes changes in relaxation behavior at specific temperatures.

Contribution

The paper introduces dynamic percolation as the underlying mechanism for mechanistic crossovers in liquid relaxation.

Findings

Dynamic percolation explains mechanistic changes at temperatures T_A and T_B in liquid relaxation.

Analysis of the KA system supports the role of mobile and immobile dynamic environments.

Percolation of dynamic environments accounts for key dynamic signatures in glass-forming liquids.

Abstract

One of the outstanding
questions regarding liquid dynamics is the
cause of the apparent mechanistic changes in relaxation at the characteristic
temperatures T
A and T
B, where, as the temperature is lowered, α relaxation
times become super-Arrhenius, and the βJG relaxation
bifurcates from α, respectively. Based on system-averaged picosecond-time
scale dynamic signatures in five molecular liquids and a Kob-Andersen
(KA) model system, we propose that these mechanistic changes arise
from the percolation of distinct dynamic environments in the liquid,
where the dynamic environment of a particle is defined by the number
of structural excitations in its first solvation shell. Analysis of
the KA system trajectories supports this idea and suggests that the
most prominent effects can be understood in terms of environments
that are mobile or immobile on a picosecond time scale. Further, the
existence and percolation of these dynamic environments can account
for many of the characteristic dynamic signatures of glass-forming
liquids.

Figures

Relaxation and transport data for propylene carbonate. (Left Axis) Red circlesα and black trianglesβJG relaxation times from dielectric spectroscopy. (Inset) The stretching exponent βKWW as a function of Temperature. (Right Axis) E A values obtained from relaxation data as [d log­(τα)/d(1/T)]1/2. Dashed lines are guides to the eye. The inset is reproduced from ref with permission from the Royal Society of Chemistry.

(A) Fits of eq to F s(q,t = 1 ps) of propylene glycol at temperatures in the range (60 to 425) K. (B) (1 – Φ0) from the KA system at T = [0.4, 0.5, 0.6, 0.7]. (C) Instantaneous excitation populations, Φ0 as a function of inverse temperature for propylene carbonate (pc), propylene glycol (pg), glycerol (glyc), sorbitol (sorb), and ortho terphenyl (otp). Solid lines are guides to the eye.

Estimated populations of environments with 0, 1, and ≥2 excitations in the first shell as a function of temperature for liquids indicated and the KA model. Green, red, and blue lines indicate P 2, P 1, and P 0 populations, respectively, identified from excitations defined at 1 ps. Green, red, and blue points in the bottom panel are environmental populations calculated from simulation trajectories. T A and T B values, where established, are marked by solid black circles and occur where P 1 and P 2, respectively, have a population of 0.24 ± 0.02. Temperatures marked with T A? or T B? fall where P 0 = 0.24 ± 0.02, and we suspect that it has been erroneously assigned. The low-temperature limits of all graphs except glycerol correspond roughly to the Vogel temperature, T 0.

(A) Populations of dynamic environments in the KA system that persist for at least 0.3 ps after being first defined, for a total persistence time of at least 1.6 ps. The open circles with plus signs indicate T A and T B. The observed percolation threshold p c = 0.14 is indicated as a gray dashed line. (B) Black, blue, and green circles give the probability that (P 0+P 1), P 0 alone, and P 2 domains that persist for ≈1.6 ps in KA will form percolating clusters. The solid lines are sigmoidal fits and are intended as guides to the eye. (C, D) Percolating and nonpercolating environments. The red particles belong to P 2 environments, while the cyan particles are a collection of those in P 0 and P 1. (C) Example percolating environment for P 2, taken from a z-axis slice of a T = 0.5 trajectory. (D) Example nonpercolating environment for P 2, taken from a z-axis slice of a T = 0.4 trajectory.

(A–C) Temporal persistence profiles of dynamic environments and temperatures indicated for the KA model. (D) Environment lifetimes, defined as the first moment of the distributions in panels A–C.

Tables

1. Characteristic Temperatures

molecule	TA [K]	TB [K]	TA? [K]	TB? [K]
propylene	291	209
carbonate	294	200
		187 −
		182
		196
		210
average	292 ± 2	195 ± 10
propylene glycol	321			280
	292			239
				251
				248
average	305 ± 20			255 ± 18
glycerol		310	338
		288	330
		285	288
		295
		283
average		294 ± 13	334 ± 6
ortho terphenyl	455	290 ,,,	341
		285	357
		293	340
		291
		279
average		288 ± 5	346 ± 10
sorbitol		340
		351
average		345 ± 8
Kob-Andersen	0.84 ± 0.09	0.435
average	0.84 ± 0.09	0.435

2. Coordination Numbers

molecule	T [K]	z
propylene carbonate	200	12.5
	220	12.4
	240	12.4
	260	12.4
	280	12.5
	300	12.4
	Value Used	12.4
propylene glycol	298	13
	Value Used	13
glycerol	300	13.1
	350	13.6
	400	13.0
	450	13.2
	Value Used	13.3
ortho terphenyl (LW)	275	10.9
	300	10.8
	400	10.5
	500	10.2
ortho terphenyl (atomistic)	267	13.3
	Value Used	11.5
sorbitol	373	12.9
	Value Used	12.9
Kob-Andersen	1.0	11.3
	Value Used	11.3

Full text

Introduction

Relaxation phenomena in liquids reveal the presence of dynamic heterogeneity, temperature-dependent mechanistic crossovers, and, ultimately, ergodicity breaking at the glass transition. Perhaps the most striking of these is the strongly super-Arrhenius behavior of primary (α) relaxation or viscosity. Goldstein? proposed that neighbor-induced dynamic constraints in liquids could produce this behavior below some onset temperature. Johari and Goldstein later identified a secondary (β_JG_) relaxation process? that distinctly separates from α relaxation at a second, lower temperature and appears to signal yet another mechanistic onset.?

Figure exemplifies dynamic behavior for propylene carbonate (PC) at these onset temperatures. The higher onset temperature is often associated with the point (T A) where the primary relaxation time (τ_α_) crosses over from Arrhenius behavior, τ_α_ ∝ Exp­[E A/k B T], E A = constant for T > T A, to super-Arrhenius behavior. Below T A, liquids evince nonergodicity on the time scale of τ_α_, in other words, dynamic heterogeneity,? which can be seen through the time dependence of the rotational relaxation function c R(t) ∝ Exp­[(−t/τ)^β_KWW_ ^]. As shown in the inset to Figure, β_KWW_ generally begins to decrease from its high-temperature asymptotic value near T A, signifying a broadening in the distribution of dynamic environments that persist on the time scale of τ_α_. At the same time, E A takes on a distinct temperature dependencequadratic for PC and other liquids, ?,? evinced by the linear behavior of E A ^1/2^ vs 1/T in Figure.

Relaxation and transport data for propylene carbonate. (Left Axis) Red circlesα and black trianglesβJG relaxation times from dielectric spectroscopy. (Inset) The stretching exponent βKWW as a function of Temperature. (Right Axis) E A values obtained from relaxation data as [d log­(τα)/d(1/T)]1/2. Dashed lines are guides to the eye. The inset is reproduced from ref with permission from the Royal Society of Chemistry.

Evidence for another mechanistic onset is found at a lower temperature, ?,? often denoted as T B. As shown in Figure, E A in PC takes on a stronger quadratic temperature dependence? at T B, and relationships between transport coefficients such as τ_α_ and the diffusion coefficient (D T) change abruptly. ?,? Additionally, β_KWW_ values approach their asymptotic minimum, and there is evidence for the emergence of a dynamically inactive, solid-like phase. ?,?,? Finally, at T B, the β_JG_ relaxation? appears to bifurcate from the α process.

Stillinger? and others elaborated on Goldstein’s dynamic constraints proposal,? creating a potential energy landscape (PEL) formalism that has been used to frame considerations of liquids and glass. In this framework, dense, amorphous packing leads to locally favored particle configurations referred to as inherent states (IS). These lie in energy basins with distributions of width, depth, and connectivity that encode the thermodynamic properties of the liquid. Dynamic properties of the liquid arise through interbasin (IB) transitions. Chandler proposed that, below T A, IB transitions must be associated with structural excitations,? localized, activated particle configurations that facilitate particle rearrangements and correspond to population at the top of interbasin barriers on a PEL. In simulations, ?,? these excitations are identified by ∼1 ps time scale interbasin barrier crossing events (hops) that occur as excited particle configurations pass over a dividing surface between two basins and rapidly rearrange into a new inherent state. It is important to note that, while excitations are identified by dynamic events, they are themselves not dynamic but thermodynamic, structural entities.

If all dynamic and thermodynamic information is encoded in the PEL, then the underpinnings of the mechanistic crossovers at T A and T B must have their origin in structures and dynamics that can be identified on a picosecond time scale since that is the shortest time on which basin occupation can be determined? and the natural time scale for interbasin barrier crossing events. ?,?,? In support of this notion, it has long been recognized that picosecond motions encode slower dynamics.? Below, we find evidence from picosecond-time scale excitations (populations at the top of interbasin barriers) at T A and T B leading us to postulate that the changes in relaxation mechanism at those temperatures are due to the percolation of distinct dynamic environments that are defined by the number of excitations in a particle’s first solvation shell. We then test this hypothesis by analyzing trajectories of the Kob-Andersen (KA) model system, finding support for and refining this idea.

Methods

QENS Data Acquisition and Reduction

The acquisition of QENS data analyzed here has been reported previously. ?,? In summary, the QENS data was acquired with the Disk Chopper Spectrometer (DCS) installed at the NG4 guide at the NIST Center for Neutron Research.? The neutron wavelength was set at 4 Å, corresponding to a momentum transfer (q) range of (0.2 to 2.8) Å^–1^, and data were binned at steps of 0.1 Å^–1^. The instrumental resolution was 0.19 meV, and the maximum energy transfer was 4.5 meV. Data at each temperature was obtained within (3 to 6) h.

Experiments were performed on cooling from the high temperature to avoid potential artifacts from crystal nucleation. We ensured no sign of crystallization even on reheating when crystal nuclei would have grown if they had been seeded at low temperatures. The absence of crystallization was confirmed by noting an absence of Bragg peaks in plots of scattered intensity vs q.

The data were corrected for (i) background from an empty can, (ii) dark count background with no neutron flux, and (iii) detector efficiencies by measuring a vanadium standard. All the data files were reduced using the DAVE software available at www.ncnr.nist.gov/DAVE. Instrument resolution was estimated as a Gaussian from sample scattering at 30 K.

QENS data were analyzed in the frequency domain and in the time domain after transformation. S(q,ω) was directly transformed to the time domain F(q,t) using DAVE software. Noise is propagated in the Fourier transform operation, and additional errors can arise through truncation and course sampling intervals. These factors are considered and standard uncertainties are calculated by the DAVE software. Errors arising were less than 0.1% of F(q,t) values. The real part of the Fourier transformed data was used to calculate F(q,t), and no filtering options were used as they were not found necessary.

Simulations of the Kob-Anderson

Model

Kob-Andersen (KA) binary LJ simulations were performed with 10,000 particles at the standard choice of 80%/20% A/B particles and reduced density ρ = 1.2. In the KA mixture, the A/B particles interact with LJ pair potentials given by σ_BB_/σ_AA_ = 0.88, σ_AB_/σ_AA_ = 0.8, ϵ_BB_/ϵ_AA_ = 0.5, and ϵ_AB_/ϵ_AA_ = 1.5. KA simulations were run with the OpenMM software, utilizing a Langevin thermostat and pairwise interactions cutoff at r cutoff ≈ 3σ_AA_. The integration time step was set as Δt = 0.0025τ, with τ being the reduced time unit of the KA model.? Simulations were performed at eight different temperatures spanning the range T = 0.4–1.1 in reduced temperature units; this range spans temperatures just below T B to temperatures above T A (Table). Each system was equilibrated for ∼ 10^3^τ, which is sufficient given they were started from the target system density ρ = 1.2. Because the primary analysis of each simulation is the computed self-intermediate scattering function (vide infra), simulation trajectory lengths and output frequency were chosen to provide sufficient resolution/statistics at both short/long time scales. Trajectory snapshots were saved at every 0.05τ time point. The specific trajectory length for each temperature simulation is given in the Supporting Information; trajectories spanned in length from 10^3^–10^4^τ, with longer trajectories generated for the lower temperature simulations.

1: Characteristic Temperatures

The self-intermediate scattering function F s(q,t) was computed from the KA simulations at the different temperatures, using the formula

$$[eqn]$$

The self-intermediate scattering function computed from the KA simulations, in tandem with the corresponding experimental scattering functions obtained from QENS experiments on the molecular liquids, are the primary dynamical data upon which the subsequent analysis of excitations is based.

Determining Characteristic Temperatures

Table displays characteristic temperatures T A and T B obtained from literature sources and calculations herein. The T A values from literature are all determined through a change from Arrhenius to super-Arrhenius temperature dependence of τ_α_. Some literature values listed under T B were obtained from a change in the temperature dependence of τ_α_, usually by plotting its derivative as in Figure. ?,?,?−? ?,? The remaining T B entries from published literature ?−? ? ? ?,?,?,?−? ? were determined as the mode-coupling critical temperature (T c) by fitting material relaxation to a form: τ_α_(T) ∝ (T – T c)^−γ^, or finding the temperature at which nonzero solutions for the mode-coupling nonergodicity factor f(k) first emerge.? Values found as T B or T c are statistically indistinguishable.

Entries marked with letter superscripts (a–c) were determined by us. In cases marked with a superscript a, characteristic temperatures were determined using the method of Stickel et al.? and published relaxation data? as shown in Figure. For entries with a b superscript, we estimated T B as the temperature where τ_α_ and τ_β,JG_ bifurcate. For sorbitol, this was done using τ_β,JG_ from Fujima et al.? to generate a second T B data point. As a check on this approach, we performed a similar analysis for PC and OTP, where T B values are well-established. T B determined in this way agreed reasonably well with T B and T c from literature. The single entry marked with a c superscript is T A for the KA model, which we determined from the onset of super-Arrhenius temperature dependence of α relaxation times. We found τ_α_ as the time of the 1/e point of F _ s _(q,t) at the peak of the structure factor, q _ max , for A particles. The τ_α vs 1/T data from which T A is derived is shown in the Supporting Information.

We now call attention to the last two columns in Table, labeled T A? and T B?. Based on an analysis we will present in Figure below, these values for T A and T B, taken from literature, fall at a what we suspect is a temperature intermediate between actual T A and T B temperatures where we see evidence for percolation of microscopic regions in the liquid that are solid-like on a picosecond time scale.

Coordination

Numbers

Coordination numbers (z) used to construct Figures and ? are given in Table and were calculated from simulated center-of-mass radial distribution functions (g(r)) according to

$$[eqn]$$

where ρ is the average density, and r min is the first minimum in g(r). The KA simulations utilized to compute g(r) have been described in the Methods Section, whereas simulations utilized to compute g(r) for glycerol, OTP, sorbitol, and propylene carbonate are described in the Supporting Information. We note that all other data/analyses of these molecular liquids are based on the experimental QENS data, and the simulations are used for no purpose other than computing coordination numbers.

2: Coordination Numbers

Identifying

Excitations

We are interested in quantifying the population of particles involved in structural excitations at the tops of the interbasin barriers. Excitations are identified in simulations ?,? by picosecond-time scale particle rearrangements over length scales on the order of 20% of a particle radius.

Signatures of inherent state and interbasin motion are also found in experimental quasielastic neutron scattering (QENS) data. A decade ago, we showed? that, at 1 ps, when diffusion is not important, the self-intermediate scattering function (F s(q,t)), the Fourier transform of the scattering function S(q,ω) from QENS data of five molecular liquids could be fit without significant residuals with a two mode function (See FigureA)

$$[eqn]$$

where σ_IS_ and σ_IB_ are distinct length scales of localized motion and Φ­(t) is the time and temperature-dependent fraction of particles executing the larger length scale motion. The strong temperature dependence of σ_IS_ and weak temperature dependence of σ_IB_ found in that work are respectively consistent with increased length scale intrabasin exploration at elevated temperatures and a constant length scale for particle rearrangements associated with interbasin barrier crossing. At all temperatures, σ_IB_ ≥ 3σ_IS_ for each material. The well-separated length scales make it trivial to separate these variables in the fitting process. The σ_IB_ value dominates the slope in FigureA at q ^2^ < 1, below the structure factor peak at q ≈ 1.4 Å^–1^, whereas σ_IS_ dominates the slope at q ^2^ > 2.

(A) Fits of eq to F s(q,t = 1 ps) of propylene glycol at temperatures in the range (60 to 425) K. (B) (1 – Φ0) from the KA system at T = [0.4, 0.5, 0.6, 0.7]. (C) Instantaneous excitation populations, Φ0 as a function of inverse temperature for propylene carbonate (pc), propylene glycol (pg), glycerol (glyc), sorbitol (sorb), and ortho terphenyl (otp). Solid lines are guides to the eye.

Vispa et al.? subsequently used Bayesian inference to evaluate models of S(q,E) from QENS data of glycerol to establish the unambiguous presence of three distinct dynamical processes with Lorentzian form; translational diffusion and two additional localized modes. The fast and slow localized modes had time scales of 100 fs and 1 ps and distinct length scales. In their model, they proposed that particles could participate in either fast or slow motion, but not both.

We subsequently performed a Bayesian model analysis for S(q,E) and F s(q,t) of propylene carbonate, also finding a three-Lorentzian model to provide the best fit. In a slight departure from the approach of Vispa et al., we proposed a heterogeneous model where all particles undergo fast localized motion, overdamped vibrations around the inherent state configurations, but only some particles participate in slightly slower and larger-scale interbasin transits, consistent with a PEL picture, giving the following postulated form for S(q,E).

$$[eqn]$$

where L _ i _ = Γ_ i π^–1^(E ^2^ + Γ i _ ^2^)^−1^, ⊗ is the convolution operator, the convolutions are over frequency (energy). Γ_D_ = D T q ^2^, where D T is the diffusion coefficient, and a _ IS _ and a IB are respectively q-dependent scattering amplitudes from vibration of inherent states and transits (hopping) over interbasin barriers. Both terms in eq include diffusive (D) and inherent state (IS) vibrational motion. In agreement with Vispa et al., we found that the two localized modes had Gaussian distributions of displacements following the functional form proposed by Rahman:?

$$[eqn]$$

where the σ_ i _ represent the characteristic length scales of motion for mode i.

The Fourier transform of eq gives the time-domain representation as

$$[eqn]$$

where τ_ i _ = ℏ/Γ_ i . On a picosecond time scale when tΓ_D < ℏ < tΓ_IB_ ≤ tΓ_IS_, eq reduces essentially to eq.

In this later work, we noted that the time-dependence of Φ­(t), σ_IS_, and σ_IB_ were also consistent with motion on a PEL lansdscape.? σ_IS_(t) increased monotonically, then plateaued, consistent with bounded exploration of a basin. σ_IB_(t) first increased, then decreased in time, consistent with hopping over a barrier, and subsequently some particles hopping back to the initial basin. We confirmed this origin for the nonmonotonic behavior of σ_IB_ through atomistic MD simulation.? Additional fits of eq to experimental F s(q, t) data, as well as time-dependent σ_IS_ and σ_IB_ from fits to eq are shown in Supporting Information.

The signature of the excitation population is found in Φ­(t). This parameter gives the fraction of particles that have been involved in an excitation-mediated IB barrier crossing up to time t. As shown in FigureB, (1 – Φ­(t)) drops quickly at first but then decays more slowly. This rapid drop is a result? of the same rapid hopping motion that is used in simulation to identify excitations, and the amplitude of the drop gives us the fraction of particles involved in excitations. The 1 ps time scale at which this rapid drop occurs is the time that σ_IB_ reaches a local maximum approximately equal to 20% of the hydrodynamic radius (see Supporting Information). This time and length scale are consistent with criteria applied in particle tracking for identifying excitations through individual particle displacements. ?,?,?

Given the above, we take Φ­(1 ps) = Φ_0_ at the population of particles in excitations, at the top of interbasin barriers. FigureC shows Φ_0_ derived from QENS data for the five liquids. The solid lines are guides to the eye, but we note a positive deviation from Boltzmann behavior at low temperatures for most of the liquids. Assuming Boltzmann statistics apply, Φ_0_ = Exp­[−E Φ/k B T]/Z, and the deviation from behavior at T > T B could be due to temperature-dependent changes in the partition function, Z. Middleton et al.? found distinct groups of interbasin energy barriers for several model systems. It is possible that one such grouping begins to become unavailable below T B, reducing Z and leading to the non-Boltzmann behavior.

In the analysis of simulation data we follow the approach of others, ?,? assigning particles to excitations when they exhibited a sufficiently large displacement δr ^2^ over a time window δt:

$$[eqn]$$

where r⃗ represents the particle’s position and δt = 1.3 ps. We use Φ_0_ determined from eq constrain the number of particles assigned to excitations through a displacement criterion, δr⃗ min so that the fraction of particles with δr⃗ ≥ δr⃗ min is equal to Φ_0_.

Results

Evidence for Distinct Dynamic

Environments

All the systems we have analyzed have Φ_0_ ≈ 0.05 at T B, leading us to search for a potential mechanistic connection. Initially, we explore the hypothesis that the number of nearby excitations presents distinctly different dynamic environments. Retaining the assumption of facilitated dynamics, we argue that only three distinct environments are relevant: those having zero, one, or at least two excitations in their immediate vicinity. We indicate populations in these environments as P 0, P 1, and P 2.

Why three environments? To keep the argumentation simple, we assume the following: (i) significant interbasin barriers exist and excitations mark population at a saddle point of the barrier, (ii) each saddle point between IS basins has only one primary reaction coordinate (a nonbifurcating minimum surface) and that the curvature at the top of the barrier is symmetric about the maximum. In this case, each excitation presents a binary choice between two inherent states for the system, each with the same probability. For convenience, we will define IS _ A,j _ and IS _ B,j _ respectively as the initial inherent state occupied by and the new inherent state available to the particles involved in excitation j.

In our idealized scenario, particles in P 0 environments with no excitations in their vicinity have no opportunities to reorganize and are temporarily trapped in their present IS. Single excitations give the system two equally probable options, occupation of IS _ A,1_ or IS _ B,1_. With no other reorganization events in the same vicinity, IS _ A,1_ and IS _ B,1_ both remain open to the system (the rearrangement is reversible), and no real relaxation has occurred. If two excitations coincide in the same vicinity we have four equally probable outcomes: {IS _ A,1_; IS _ A,2_}, {IS _ B,1_; IS _ A,2_}, {IS _ A,1_; IS _ B,2_}, {IS _ B,1_; IS _ B,2_}. Without considering coupling, each possibility is equally weighted, and there is a 75% chance that at least one new IS will be occupied, but the system can still return to its original state. If, in our simple model, we allow a particle rearrangement to induce a strain field? that biases inherent state potentials for nearby particles, then we must consider coupling between excitations. In this case, a IS _ A,1_ → IS _ A,2_ rearrangement in excitation 1 biases the potential of the nearby particles involved in excitation 2, favoring either the IS _ A,1_ or IS _ B,2_ choice of that excitation. Once this biased choice is made for excitation 2, the IS _ B,1_ state is stabilized, and a relaxation has occurred. By this argumentation, P 0 and P 1 are fundamentally distinct from P 2, in that the former two cannot relax, but the latter can. Environments with two nearby excitations behave qualitatively the same as environments with three or more nearby excitations.

To quantify environment populations, we must define a “vicinity” as used in the paragraph above. Assuming that dynamics are influenced primarily by the structure of the first solvation shell, ?,? we classify local dynamic environments by the number of molecules within the first shell involved in excitations.

Using measured Φ_0_ values, coordination numbers (z), and mild assumptions, we can estimate populations P 0, P 1, and P 2. In Figure, we present these population estimates for all systems studied. We arrive at these population estimates by assuming that excitations are randomly distributed in space, that P 0 = 1 at T = 0, and that z is temperature-independent. Under these assumptions we can write a set of ordinary differential equations that are exact in the limit Φ_0_ → 0: dP 0/dΦ 0 = −(z + 1) P 0, where z is the coordination number, dP 1/dΦ 0 = (z + 1) (P 0 – P 1), and dP 2/dΦ 0 = (z + 1) P 1. These equations do not account for correlation effects that come into play at higher Φ_0_ values, but for systems such as those of interest here, with z ≈ 12, the deviations are inconsequential. (See Supporting Information for comparison to exact lattice calculations.) Additionally, the bottom panel in Figure compares the estimated environment populations and those directly tabulated from the simulations, showing only minor deviations. Coordination numbers for each system are calculated from simulation as described in Methods and shown in Table.

Estimated populations of environments with 0, 1, and ≥2 excitations in the first shell as a function of temperature for liquids indicated and the KA model. Green, red, and blue lines indicate P 2, P 1, and P 0 populations, respectively, identified from excitations defined at 1 ps. Green, red, and blue points in the bottom panel are environmental populations calculated from simulation trajectories. T A and T B values, where established, are marked by solid black circles and occur where P 1 and P 2, respectively, have a population of 0.24 ± 0.02. Temperatures marked with T A? or T B? fall where P 0 = 0.24 ± 0.02, and we suspect that it has been erroneously assigned. The low-temperature limits of all graphs except glycerol correspond roughly to the Vogel temperature, T 0.

The T A and T B values in Figure and their uncertainties are derived from literature values listed in Table. We find it striking that one of the populations P _ i _ has a value 0.24 ± 0.02 at each of these temperatures. In most cases, P 1 and P 2 = 0.24 within uncertainty at T A or T B respectively. In this figure, vertical error bars at T A and T B are obtained by projecting temperature uncertainty onto the associated population line.

The colored data points in the bottom panel of Figure show environment populations calculated directly from simulation trajectories. Although excitations are inhomogeneously distributed in the KA model,? our P _ i _ population estimations are quite close to measured values in this system, suggesting that the assumption of randomly distributed excitations is not too severe. Aside from the assumption of random distribution, the populations are estimated with no free parameters.

The probability that a single characteristic temperature for any individual liquid would randomly fall at P _ i _ = 0.24 is approximately equal to the uncertainty in that P _ i _ value (2σ_ i ) at that temperature, divided by the total range that P _ i _ takes. Accordingly, we estimate the log probability of the null hypothesis for N such temperatures as ln(H 0) ≈ N ln[2σ i _/range­(P _ i _)], where we conservatively assign range­(P _ i _) = 0.35 (although P 0 and P 2 range from 0 to 1). H 0 values for individual liquids are modest, but H 0 ≈ 10^–8^ for all systems combined. This strongly supports our initial hypothesis that there is a mechanistic tie between populations of dynamic environments and the mechanistic crossovers at T A and T B. A closer investigation of the data allows us to refine this initial hypothesis.

The data of Figure indicates that mechanistic changes to relaxation occur for each liquid at a specific fraction (0.24) of a given dynamic environment. To us, this suggests that the percolation of dynamically distinct environments underlies these mechanistic changes. Consistent with this idea, we note that for OTP, glycerol, and PG, characteristic temperature values thought to have been T A or T B have been reported at precisely the temperatures where P 0 = 0.24. If the percolation hypothesis is correct, this would correspond to a solidification temperature where, on a picosecond time scale, a solid-like domain would first percolate. These temperatures are reported in Table under column headings “T A?” and “T B?”.

Despite the considerable circumstantial evidence for percolation, 0.24 ± 0.02 is much higher than expected for a site percolation critical fraction (p c) for amorphous systems with coordination numbers in the range (11–13) for as for systems considered here (see Table). For fcc and hcp lattices, which have z = 12, p c ranges from (0.195 to 0.199). Random-packed systems are expected to have p c values equal to or slightly smaller than regular lattices with the same z.? We turn to an analysis of the KA model trajectories to test whether the patterns in Figure actually signify percolation events, and if so, what precisely is percolating. Of course, the KA model does not represent the wide range of intermolecular interactions represented by the other liquids, but all these liquids exhibit what appear to be very similar or the same mechanistic crossovers at T A and T B so it is reasonable to assume that the detailed topology of the PEL landscape may not matter much, and that what matters is there is a landscape with significant barriers. Additionally, we have previously demonstrated that on time scales and length scales associated with interbasin hops, motion in molecular systems can be mapped onto the behavior of simple systems such as the KA model. ?,? Below, we test the percolation hypothesis in the KA system by directly calculating percolation thresholds for dynamic environments. As we will see below, analysis of the KA model provides evidence that the high p c values inferred from the data in Figure are an artifact of counting both transient and longer-lived dynamic environments when percolation of only the longer-lived environments might be relevant.

Percolation of Mobile and Immobile Domains

We identify dynamic environments for each A particle in the KA system as follows: First, all A particles are given a label value, l = 1 if the particle is involved in an excitation or l = 0 if not. First coordination shells are then defined as all A particles within an interaction radius, R AA = 1.42 σ, and the environment of the center particle is assigned as P 0, P 1, and P 2 if the sum of l values over the central particle and its first shell is 0, 1, or ⩾2, respectively.

Once environments are assigned, we compute the percolation probability for each dynamic environment P 0, P 1, and P 2. Percolation for a given type of environment (e.g., P 0) is defined as when the particles form a system-spanning network composed of continuous connections between particles, fully spanning the periodic simulation box. Connections are contacts between nearest neighbor particles, using the criteria above R AA = 1.42 σ to define nearest neighbors. For each simulation snapshot, a value of “1” is assigned if the environment type (e.g., P 0) percolates as defined above, or “0” otherwise; the percolation probability is then the average over the simulation trajectory. Percolation of each frame is determined using a recursive search algorithm of Edvinsson et al.? It is known that percolation behavior follows a “step function”, exhibiting a sharp transition at the percolation threshold “p c”, with systems percolating/nonpercolating for concentrations above/below p c. This sharp transition in percolation behavior is only rigorous in the infinite system size limit,? and for finite system sizes (i.e., with periodic boundary conditions) p c must be estimated, as the observed percolation transition occurs over a concentration range. For the KA system sizes investigated here, we estimate a few percent uncertainty in determining p c due to finite size effects,? which is acceptable for our purposes given the corresponding uncertainties in determining populations P 0, P 1, and P 2.

As a control, we also compute percolation for a random distribution of particles. We randomly assigned excitation labels to varying fractions of A particles and found that each environment percolated at p c = 0.2 as expected. When excitations were identified based on particle mobility as described above, we found that domains P 0, P 1, or P 2 percolate with p c = 0.14, indicating that excitations are not randomly distributed but form local clusters, consistent with the results of others.? This confirmed that the apparent p c = 0.24 value from Figure was aberrant.

Excitations and dynamic environment populations were assigned using a rolling 1.3 ps window, applied at 0.1 ps intervals. We filter the most unstable environment assignments by considering only those dynamic environments that retain the same designation over three consecutive assignment intervals. Applying this criterion, we obtain environment percolation at the expected temperatures. P _ i _ values for persistent environments are shown in FigureA. p c is exceeded at T A for immobile particles, those that persist in P 0 or P 1 domains. Likewise, the population of mobile particles, those that persist in P 2 environments, drops below p c at T B.

*(A) Populations of dynamic environments in the KA system that persist for at least 0.3 ps after being first defined, for a total persistence time of at least 1.6 ps. The open circles with plus signs indicate T A and T B. The observed percolation threshold p c = 0.14 is indicated as a gray dashed line. (B) Black, blue, and green circles give the probability that (P 0+P 1), P 0 alone, and P 2 domains that persist for ≈1.6 ps in KA will form percolating clusters. The solid lines are sigmoidal fits and are intended as guides to the eye. (C, D) Percolating and nonpercolating environments. The red particles belong to P 2 environments, while the cyan particles are a collection of those in P 0 and P

1. (C) Example percolating environment for P 2, taken from a z-axis slice of a T = 0.5 trajectory. (D) Example nonpercolating environment for P 2, taken from a z-axis slice of a T = 0.4 trajectory.*

From top to bottom, the three panels in FigureB show percolation probabilities for all slow environments (particles in either P 0 or P 1 environments), solid-like environments P 0, and fluid environments (P 2). We see that (P 0 + P 1) percolates at T A, rather than simply P 1 as suggested in Figure. The percolation of solid-like environments occurs at a temperature intermediate between T A and T B, and the bottom panel shows that mobile P 2 domains cease percolating at T B.

FigureC and D provide examples of percolating and nonpercolating environments for P 2, above and below T B ≈ 4.3 at T = 0.5 and T = 0.4, respectively. The images are 2D slices along the z-axis of ∼1.3σ in depth from our 3D trajectories. All particles from the slice were projected onto a single plane for visualization.

Considering only longer-lived dynamic environments leads to reasonable values for p c and percolation at the expected temperatures, but there is a notable difference in the ratio of P 0 to P 1 populations defined over shorter times and slightly longer periods, with P 0/P 1 increasing significantly for longer-lived environments. To gain insight into this, we consider the persistence characteristics of each environment. We define a time-dependent environment population P _ i (Δt) as the fraction of particles that continuously meet the criteria for inclusion in envronment P _ i _ for a time interval Δt, and η i (Δt) = P _ i (Δt)/P _ i (0) where P _ i (0) = P _ i _ as defined above. FigureA–C shows these lifetime plots for A particles of the KA model, and panel D shows average environment lifetimes, calculated as first moments of η i (Δt). P 1 environment assignments are short-lived, with characteristic lifetimes (τ P1) that are only weakly temperature dependent and less than 300 fs at all temperatures. This profile suggests that P 1 environments are found primarily at interfaces between P 0 and P 2 domains. The temperature dependences of τ P0 and τ_ P2_ suggest that these environments form bulk domains. These would be solid-like and fluid-like, respectively. With lifetimes no larger than the inverse relaxation attempt frequency, it is understandable that P 1 environments on their own do not seem to have an impact on dynamic mechanism. By contrast, P 0 and P 2 lifetimes are such that significant fractions of these environments persist long enough to influence local relaxation attempts.

(A–C) Temporal persistence profiles of dynamic environments and temperatures indicated for the KA model. (D) Environment lifetimes, defined as the first moment of the distributions in panels A–C.

Discussion

Above, we have presented an analysis of experimental and simulation data suggesting that the crossovers in relaxation mechanism at T A and T B arise from the percolation of distinct dynamic environments in the liquid on a picosecond time scale. Others have also found that percolation appears to play a role in thermodynamic? and dynamic ?−? ? ? ? ? properties of liquids and glasses. Of these previous reports, the works of Novikov et al.? and Gao et al.? are most directly relevant to the present study.

Gao et al.? recently showed that β_JG_ and α relaxation, measured through the shear modulus, occur respectively on the shortest time scale at which mobile particles begin to percolate and immobile particles cease to percolate for a given temperature. In their simulations, they defined all particles as either mobile or immobile depending on their displacements over the time window in which the material relaxation was probed.

Novikov et al.? showed that experimental values for T B (the mode-coupling T c) and the Vogel temperature, T 0 (far below T g) are consistent with percolation of a mobile domain at T 0 and percolation of a solid-like domain at T B, assuming percolation thresholds of 15%, and classifying all particles as in solid-like or mobile environments depending on whether they exhibit picosecond-time scale mean-squared displacements in the upper half or lower half of an assumed log-normal distribution.

On their face, the Gao? and Novikov? interpretations of mobile particle percolation appear inconsistent. Together, they seem to suggest that β_JG_ relaxation should occur on a picosecond time scale at T 0, whereas τ_β,JG_ is much slower, usually on the order of μs to ms at the much higher temperature, T g. Our analysis resolves this issue and is consistent with results from both papers.

Given that Novikov et al.? considered large excursions on a ps time scale, they would have considered particles in P 1 and P 2 environments as mobile. Since their analysis dealt with low temperatures, P 1 environments would be three to four times more prevalent than P 2. On the other hand, Gao et al.? would have considered only particles in persistent P 2 environments as mobile since they performed their analysis on time scales roughly (3 to 5) orders of magnitude larger, and the lifetime of P 1 environments is so short. Our analysis shows that Gao’s mobile particles (those in P 2), percolate on a ps time scale at T B. Below this temperature, P 2 percolation would occur on increasingly long time scales, consistent with Gao’s results. Likewise, Novikov’s mobile particles in P 1 environments appear to cease percolating near T 0 for most of the systems we have analyzed in Figure. P 1 populations cross the apparent percolation line near the lowest temperatures plotted in these systems, which, in all cases except glycerol, correspond to T 0.

We likewise attribute the incompatibility of the Gao and Novakov analysis of immobile particle percolation to the assignment of immobile particles. Novikov et al. inferred the presence of immobile particles and estimated that they should percolate near T B. As we show in Figures and ?, solid-like domains do percolate near, but at a slightly higher temperature than T B. In other published work, this higher temperature percolation event seems to have been mistaken for T B or T A in OTP, glycerol, and propylene glycol (see Table). Gao et al. assign mobile and immobile particles to be consistent with a two-Gaussian van Hove function, similar to what we have done, but on a longer time scale, so their assignments should correlate closely to the immobile (P 0+P 1) assignments in our work, which, on a picosecond time scale percolate at T A. As Dyre points out,? this percolation, which Gao showed occurs at longer times for lower temperatures, is essentially a facilitation, allowing α relaxation to occur while only sampling the smallest barriers to reorganization.

In addition to the evidence presented above for the proposed percolation scenario, we discuss below how this proposed scenario, along with the assumption that excitations are required for relaxation, leads us to plausible pictures for mechanistic crossovers in liquid dynamics. To reiterate from previous sections, P 0 environments are solid-like and unable to relax, and P 1 environments support only single IB, reversible transitions. In P 2 environments, reversible IB transitions can be entropically trapped by a quasi-synchronous and nearby IB transition, leading to nonreversing relaxation events.

At T > T A in the proposed scenario, both P 0 and P 1 exist but neither forms extended domains. FigureD shows that both these environments are transient on the time scale of α relaxation at T > T A, which is also evinced by the Arrhenius temperature dependence and lack of obvious dynamic heterogeneity at the τ_α_ time scale. The existence of slow environments that persist on a ps time scale at T > T A is demonstrated by the bimodal structure of the van Hove function at 1 ps.? We understand these results from P 0 and P 1 existing primarily as individual loci or small droplets surrounded by P 2 at T > T A. Intimate proximity with P 2 environments means that P 1 or P 0 particles will quickly become involved in an excitation, leading to the interconversions P 0 ↔ P 1 ↔ P 2 on the time scale τ_Φ_, where τ_Φ_ < 0.25 ⟨τ_α_⟩ for T ≥ T A.? Interconversion of these environments on such a short time scale would obfuscate dynamic heterogeneity effects, and α relaxation would appear to be homogeneous and collisional.

At T < T A, signatures of dynamic heterogeneity and excitation-facilitated hops become detectable, ?,?,? which we understand as a consequence of increased persistence of slow domains (P 0 + P 1) concomitant with their percolation. The appearance of extended P 0 + P 1 environments is consistent with expectations of a spinodal at this temperature.?

At all temperatures T > T B, where extended P 2 environments exist, multiple pathways will be available to each particle for interbasin transits in those environments. Thus, particles in extended P 2 domains undergoing an interbasin transit will likely have a particle in their first shell involved in a nearly independent but simultaneous interbasin transit. In these cases, the local structure would rearrange somewhat within the period of the barrier transit. If close enough, those nearby rearrangements would tend to alleviate local stress and reduce the extent of displacement required to establish a new inherent state. This accounts for the reduced hop length (σ_IB_) observed at elevated temperatures in actual liquids.? This scenario allows us to define small-step diffusion as instantaneously facilitated relaxation dynamics and understand why small-step diffusion is observed only at T > T B, ?,?

?,? where extended P 2 domains exist. The existence of extended P 2 domains only at T > T B further explains the observation that saddle points are encountered when quenching simulations to identify inherent states from this temperature regime,? and with a localization transition (loss of access to saddle points) at T < T B.

Below the temperature where P 0 percolates, solid-like domains should begin to persist long enough to impact dynamics. Accordingly, there is a shift in the temperature dependence of the Brillouin line at precisely this temperature in PC.? Additionally, this temperature seems to have been mistakenly identified as T A or T B in OTP, glycerol, and PG. Finally, two-level systems are known to disappear in glasses prepared with a sufficiently low equilibrium structure.? This is predicted in the scenario proposed here when P 1 drops below the detection limit in the equilibrated glass at some temperature well below T B. Finally, we note in passing that Figure suggests P 1 environments may drop below a percolation threshold at the Kauzman temperature. Investigating the implications (if any) of this on dynamics will be addressed at a later time.

Summary

We provide evidence that the percolation of mobile and immobile domains on a picosecond time scale drives mechanistic crossovers in relaxation at T A and T B. We also show that such percolation events provide mechanistic explanations for a wide range of the dynamic phenomena documented at these temperatures. Discriminating between three distinct picosecond-time scale dynamic environments also allows us to resolve apparent contradictions between results of previous works that also propose dynamic percolation. ?,? Assuming that this percolation hypothesis is substantially correct, we are left with two compelling questions: (1) Why do we find percolation for mobile and immobile regions only when we define dynamic environments as those that persist for at least 300 fs after they have been identified, and (2) Why do all the experimental liquids seem to behave just as the KA system in Figure, with apparent p c values of 0.24 when dynamic environment populations are estimated from Φ_0_?

The KA model cannot reproduce all the dynamic features that arise from the detailed chemical structures of molecular liquids. On the other hand, we have previously given examples where KA maps onto molecular systems when time scales and length scales are properly averaged, ?,?,? the dynamic features we are concerned with hereexcitations and interbasin hopsappear to provide that averaging to a large extent. Regarding the environment persistence question, ≈1 ps, where Φ­(t) transitions to a slower decay (see Figure) is the shortest time at which it is physically meaningful to define an excitation or a dynamic environment. Likely, some of the excitations identified at this time are merely rare deformations of unusually elastic inherent states. This would lead to a transient misidentification of some dynamic environments, which would be rectified with a slightly extended persistence criterion.

Fully addressing these questions will require analyzing atomistic simulation trajectories for at least each of the liquids in Figure, as we have done for the KA system. This is presently out of reach for two reasons. The first is that T A and T B are not known for the force field Hamiltonians of the liquids other than the Kob-Andersen model. There is no question that T A and T B from simulations of propylene carbonate, glycerol, sorbitol, and OTP will vary from the values obtained from the liquids themselves. For example, in PC we have found that even a very sophisticated, polarizable ab initio force field cannot reproduce the experimental scattering function F(q,t) with enough fidelity for correspondence of excitation populations between experiment and simulation.? This precludes any correspondence between simulation thresholds (e.g., percolation) and experimental thresholds of P 0, P 1, and P 2 as determined in Figure. Determination of T A and T B for each liquid force field/Hamiltonian would require extensive simulation beyond that presented in this work. The second, more compelling reason is that precisely how one defines excitations in a molecular liquid simulation has yet to be worked out.

Finally, the transition between the picture suggested by Figures and ? requires specific lifetimes for the three dynamic environments. If excitations are prone to cluster, then P 0 and P 2 will form bulk domains and P 1 would be found mostly at interfaces as seems to be the case here. If this is the case, then the persistence time scales come down to the size and fractal dimension of the domains, and the fundamental time scale for excitation mobility. The latter should be several multiples of the zero-crossing time of the particle velocity autocorrelation time, and should not depend strongly on temperature or molecular system. If the characteristic temperatures T A and T B are defined by some function of the size and fractal dimension of the P 0 and P 2, then these are constrained to maintain similarity in some gauge for all systems. We look forward to exploring these questions.

Supplementary Material