The transport of small penetrants through disordered materials with glassy dynamics is encountered in applications ranging from drug delivery to chemical separations. Nonetheless, understanding the influence of the matrix structure and fluctuations on penetrant motions remains a persistent challenge. Here, we use event-driven molecular dynamics to investigate the transport of small, hard-sphere tracers embedded in matrices of square-well particles. Short-range attractions between matrix particles give rise to reentrant dynamics in the supercooled regime, in which the liquid’s relaxation time increases dramatically upon heating or cooling. Heating results in a “repulsive” supercooled liquid where relaxations are frustrated by steric interactions between particles, whereas cooling produces an “attractive” liquid in which relaxations are hindered by long-lived interparticle bonds. Further cooling or heating, or compression, of the supercooled liquids results in the formation of distinct glasses. Our study reveals that tracer transport in these supercooled liquids and glasses is influenced by the matrix structure and dynamics. The relative importance of each factor varies between matrices and is examined in detail by analyzing particle mean-square displacements, caging behavior, and trajectories sampled from the isoconfigurational ensemble. We identify features of tracer dynamics that reveal the spatial and temporal heterogeneity of the matrices and show that matrix arrest is insufficient to localize tracers.

The transport of dilute small molecules or particles within disordered media affects the delivery of drug molecules encapsulated in hydrogels,1,2 the efficacy of polymeric gas separation membranes in capturing carbon dioxide or purifying natural gas,3,4 and the movement of DNA through the crowded cytoplasm during transformation and transcription.5–7 In these processes, penetrant dynamics may couple to the structure and/or to the slow relaxations of the surrounding matrix. As an example, increasing the density of an arrested, disordered matrix leads to anomalous diffusion at a critical density or localization at higher densities.8,9 Likewise, matrix mobility affects gas diffusion in polymeric membranes10–12 and the transport of cytoskeletal and cytoplasmic constituents within the cell.5,13 Understanding the effects of the matrix structure and dynamics on penetrant transport, however, remains a persistent challenge.

Recent progress has been made in understanding tracer transport in complex matrices using well-controlled colloidal models. Anomalous tracer dynamics and localization have been observed in model disordered media consisting of colloidal particles fixed in gel-like14 and liquid configurations.15,16 Similarly, tracer dynamics have also been shown to be coupled to matrix motion, crossing over from localized to diffusive behavior as the matrix relaxes.17,18 This coupling, however, depends on the relative time scales between the tracer and matrix dynamics and also the nature of the matrix relaxations. It is unclear, for example, how the onset of nonergodic, glassy dynamics19,20 may influence this coupling. Moreover, it is also unclear how tracer coupling is affected by the nature of the matrix relaxations, which can be qualitatively altered by modulating the interactions between matrix particles.21,22

In this study, we use event-driven molecular dynamics (MD) to investigate the transport of tracer particles in a model colloidal glass-former consisting of a square-well fluid with short-ranged attractions.21,23 In the supercooled liquid regime, this system exhibits reentrant dynamics characterized by a marked increase in the liquid’s relaxation time upon heating or cooling. Whereas heating produces a “repulsive” glassy liquid in which relaxations are hindered by steric interactions between particles, cooling results in an “attractive” glassy liquid where relaxations are frustrated by long-lived interparticle bonds.21,23 Distinct glasses can be prepared by further heating or cooling of the repulsive and attractive liquids, respectively, and by compression. Here, we examine the dynamics of tracer particles of a critical size, known to exhibit anomalous transport,17,18 embedded in these liquid and glass matrices. We find that tracer transport is affected both by intermediate- and long-time matrix dynamics and by the matrix structure. Intriguingly, sufficiently large local fluctuations in arrested matrices that do not relax on long time scales can allow tracers to escape cages and recover diffusive behavior. In strongly quenched matrices, however, tracer dynamics are primarily determined by the structure of matrix cages. Our results identify the relative contributions of the matrix structure and dynamics on tracer motions in attractive and repulsive glassy matrices and thus provide a framework to understand transport processes in slowly relaxing materials.

Event-driven MD simulations were performed to investigate the transport behavior of tracer particles in glassy matrices. The matrices were modeled using a well-studied equimolar binary (AB) glass-forming mixture that exhibits reentrant dynamics.21,24,25 Following Refs. 21, 25, and 26, the matrix species were assigned unit masses (m = 1) and a hard-core diameter ratio of σAA:σBB = 1.2:1 to frustrate crystallization (Fig. 1). The matrix particles interacted through a short-range square-well potential with depth u0 = 1 and width Δij = 0.03(σij + Δij) for each pair type i, j ∈ A, B, where σij=12(σii+σjj). In the discussion to follow, we adopt conventional simulation units in which Boltzmann’s constant kB = 1 and σBB, u0/kB, and σBB(m/u0)1/2 are the fundamental measures of length, temperature, and time, respectively. To account for the influence of temperature on particle dynamics, we also define a thermal time scale τ = tD0/Dref, where t is the nominal simulation time, D0=σBBkBT/m is the thermal diffusivity at temperature T, and Dref is the reference value of D0 at T = 1.23 

FIG. 1.

(a) Rendering of a representative configuration with the tracers (red) embedded in a matrix of A (green) and B (blue) particles, which have a hard-core diameter ratio of σAA:σBB = 1.2:1. (b) Rendering of tracer trajectories within G0.20 (an attractive glass) over a duration of 360τ. (c) Two-dimensional projection of a tracer trajectory 800τ in duration illustrating tracer cage rearrangement. Time evolution of the trajectory is indicated by color scale.

FIG. 1.

(a) Rendering of a representative configuration with the tracers (red) embedded in a matrix of A (green) and B (blue) particles, which have a hard-core diameter ratio of σAA:σBB = 1.2:1. (b) Rendering of tracer trajectories within G0.20 (an attractive glass) over a duration of 360τ. (c) Two-dimensional projection of a tracer trajectory 800τ in duration illustrating tracer cage rearrangement. Time evolution of the trajectory is indicated by color scale.

Close modal

We examined tracer dynamics in matrices with N = 1372 particles at six different state points specified by {ϕ, T}, where ϕ is the volume fraction of matrix particles (Fig. 2; see the supplementary material). For notational convenience, we refer to these samples as LT or GT (liquid or glass, respectively), where T is the sample temperature. We considered two ergodic liquid states at ϕ = 0.610 with approximately equal long-time diffusion coefficients Di/D0: a high-temperature repulsive glassy liquid (L1.05) and a low-temperature attractive glassy liquid (L0.35), where Di is the nominal diffusion coefficient.18 Two glasses were also prepared at the same temperatures (G1.05, G0.35) and increased matrix volume fraction ϕ = 0.635. Similarly, a glass (G0.20) with ϕ = 0.610 was prepared at T = 0.20, a lower temperature than the attractive glassy liquid. Finally, we also studied a hard-sphere glass at ϕ = 0.610 and T = 1.00 (HSG1.00), which is equivalent to a square-well glass in the high-temperature limit (Fig. 2).

FIG. 2.

State diagram from Ref. 23 for the model square-well glass former. Symbols denote the locations of the liquid (L) and glass (G) matrices investigated in this study. Arrows denote hypothetical heating/cooling (up/down vertical arrows) and compression (horizontal arrows) protocols for preparing glasses from the attractive and repulsive ergodic liquids (L0.35 and L1.05, respectively). These hypothetical protocols are used to illustrate the relationships between matrices; the actual protocol used to prepare each sample is described in the text.

FIG. 2.

State diagram from Ref. 23 for the model square-well glass former. Symbols denote the locations of the liquid (L) and glass (G) matrices investigated in this study. Arrows denote hypothetical heating/cooling (up/down vertical arrows) and compression (horizontal arrows) protocols for preparing glasses from the attractive and repulsive ergodic liquids (L0.35 and L1.05, respectively). These hypothetical protocols are used to illustrate the relationships between matrices; the actual protocol used to prepare each sample is described in the text.

Close modal

We embedded Nt = 10 tracers into each supercooled liquid or glass matrix. Interactions between the matrix and tracer particles were modeled using purely repulsive, hard-sphere collisions. The diameter of the tracer σtt was chosen such that δt = σtt/σAB = 0.35, which is approximately the size ratio where tracers exhibit anomalous dynamics and couple to matrix fluctuations and relaxations.17,18 Larger tracers exhibit dynamics similar to matrix particles, whereas smaller tracers can traverse through interstitial voids and thus are largely unaffected by matrix fluctuations.17,18

The supercooled liquid matrices were equilibrated at their respective T and ϕ. The glass matrices, by contrast, were prepared by first incrementally compressing the system to a volume fraction of ϕ = 0.610 by increasing the particle radii. The systems were then equilibrated in the NVT ensemble at T = 0.55. Following equilibration, the samples were either instantaneously thermally quenched to their final temperature by rescaling the particle velocities (G0.20 and HSG1.00) or compressed to ϕ = 0.630 in increments of Δϕ = 0.010, followed by additional compression to ϕ = 0.635 in a single step of Δϕ = 0.005, and then thermally quenched (G0.35 and G1.05). After each compression increment, the glasses were simulated for 10τ at constant ϕ and T to relax compression-induced stresses. This protocol was employed in earlier studies,23,26 where it was shown to not qualitatively affect dynamics beyond the microscopic regime.26 For the HSG1.00 sample, the attractive square-well interactions were removed after the thermal quench. All glasses were subsequently aged for a waiting time twτmax (see the supplementary material), producing tw-invariant trajectory data up to the maximum observation time (τmax ≈ 105).23 Statistical properties were calculated by averaging over trajectories computed for ns = 5–20 independent samples prepared for each type of matrix using the protocols described above.

To characterize the restriction of tracer particle motions by the matrices, we performed cage analysis using the method of Doliwa and Heuer.27,28 Their method is based on the assumption that the sequential displacements of caged particles will be directionally anticorrelated. Consider an initial displacement of a caged particle Δr01=r(Δt)r(0) over a time interval Δt on which Δr01 is comparable to the characteristic cage size, where r is the particle’s position vector. Because the neighbor cage restricts further motion along Δr01, the displacement over the next time interval Δr12=r(2Δt)r(Δt) should, on average, be anticorrelated. The displacement Δr12 can be projected onto the unit vector of the preceding displacement Δr01, yielding the caging displacement projection (CDP) x12Δr12Δr01|Δr01|, which is parallel to Δr01. For a caged particle, the ensemble-averaged CDP ⟨x12⟩ is negative and its magnitude grows with |Δr01|. In hard-sphere supercooled liquids, x12=c|Δr01| for small displacements where the cage has not been broken. Larger magnitudes of the proportionality constant c indicate greater displacement memory during caging.27 For larger, cage-breaking displacements, by contrast, ⟨x12⟩ is independent of |Δr01|.29 

In analyzing tracer dynamics, we first compute the non-Gaussian parameter for particle displacements α2t) (see the supplementary material). The maximum in α2t) signifies the time scale τ* on which the per-particle variance in tracer dynamics due to caging and matrix rearrangement is greatest.30,31 We use Δt = τ* as the time interval for computing tracer displacements Δr01 and Δr12.29 The magnitudes of tracer displacements at τ* vary across different matrices. To account for this variation, we report normalized quantities x̃12C1x12 and r̃01C1|Δr01|, where C is the square-root of the tracer mean-square displacement (MSD) Δr2 at lag time τ*.

We performed simulations in the isoconfigurational ensemble (IE) to isolate the effects of the matrix structure on particle dynamics. In this approach, an ensemble of separate simulations is run over a fixed time interval. Each simulation starts from the same initial particle configuration but with a different set of randomly assigned momenta.32 For each matrix system, we analyzed nc,iso = 50–100 configurations (see the supplementary material) extracted from independently prepared samples prepared using the procedures described in Sec. II A. Each of the nc,iso configurations was used to initialize nt,iso = 80 short MD trajectories. Initial particle momenta for the MD trajectories were randomly drawn from the Maxwell-Boltzmann distribution at set temperature T.32,33 Isoconfigurational averages ⟨⋯⟩iso were computed from statistics collected from the MD trajectories. Specifically, to characterize the mobility of individual particles, we calculated the dynamic propensity DPi(t)=(ri(t)ri(0))2Δri2iso, where Δri2 is the ensemble-averaged mean-square displacement (MSD) of the ith particle. This quantity is the second moment of the particle displacement distribution, computed by averaging over the trajectories of particle i. When each particle’s mobility is equal to the average mobility Δri2iso, DPi will be unity.34,35 Thus, examination of this quantity for all particles provides insight into the spatial distribution of dynamic heterogeneity.

To characterize the shapes of tracer rearrangements, we calculated the mass M of a trajectory as a function of its radius of gyration Rg. This analysis is performed by overlaying the trajectory on a cubic lattice composed of unit cells with edge length σtt. The trajectory mass M is evaluated by assigning each cell unit mass and summing over the unique cells visited by the trajectory.36,37 To remove the effects of the initial ballistic motion, we coarse-grain the trajectories over a time scale τcg such that Δr2(τcg)=σtt2 for tracers within a given matrix. The trajectories are then resampled to ensure that successive frames are separated by a time interval τcg. From the coarse-grained trajectories, we compute the radius of gyration Rg=1nngi=1ncg(xixavg)2, where ncg is the number of coarse-grained points in the trajectory, xi is the position of the ith coarse-grained point, and xavg=1ncgi=1ncgxi is the mean position.

We first characterize the dynamics of the different glassy matrices through the ensemble-averaged mean-square displacement (MSD) Δr2 (Fig. 3), focusing on the intermediate-time and long-time dynamics of the large matrix particles. The intermediate-time dynamics are influenced by cage rattling and interactions between matrix particles, whereas long-time dynamics are controlled by the ability of the matrix to relax when particles escape their local cages. Generally, the MSDs exhibit the expected behavior for glassy matrices:21,23,38 intermediate relaxations are suppressed in attractive matrices relative to those in comparable repulsive matrices and long-time relaxations are suppressed in vitrified samples. Similar behaviors can be observed in the matrices’ self-intermediate scattering functions, Fs(q,τ)=1Nik=1Niexp[jq(rk(τ)rk(0))], where q=|q| is the wavevector magnitude, j=1, Ni is the number of particle species i, and the brackets indicate an ensemble average (see the supplementary material). Detailed comparisons between different matrices reveal further insights into relaxation processes of glassy matrices.

FIG. 3.

Mean-square displacement Δr2 for matrix species A. (a) Δr2 for all matrices with ϕ = 0.610. (b) Δr2 for glasses with ϕ = 0.635 compared to the corresponding liquids with ϕ = 0.610. The two matrices (L0.35 and L1.05) exhibit similar long-time diffusivities D/D0 ≈ 6.5 ± 0.8 × 10−7. The black solid line indicates a power-law slope of one.

FIG. 3.

Mean-square displacement Δr2 for matrix species A. (a) Δr2 for all matrices with ϕ = 0.610. (b) Δr2 for glasses with ϕ = 0.635 compared to the corresponding liquids with ϕ = 0.610. The two matrices (L0.35 and L1.05) exhibit similar long-time diffusivities D/D0 ≈ 6.5 ± 0.8 × 10−7. The black solid line indicates a power-law slope of one.

Close modal

The MSD of the large matrix particles in L1.05 is approximately constant on lag times τ ≈ 100–103, indicating interparticle caging [Fig. 3(a)]. By contrast, Δr2 for L0.35 exhibits a small plateau at τ ≈ 10−1 followed by an increasing, subdiffusive power-law Δr2τβ that extends to τ ≈ 103, where β ≈ 0.32. The small plateau corresponds to the length scale of interparticle bond formation, whereas the power-law region signifies a transition from dynamics dominated by bonding at small τ to dynamics dominated by caging at larger τ.21 The smaller values of Δr2 indicate that particles in L0.35 are more localized than those in L1.05, likely due to the formation of interparticle bonds.39 Thus, on intermediate time scales, the liquids have different relaxation mechanisms. On long time scales, however, the liquids exhibit nearly identical dynamics. The MSD for both L0.35 and L1.05 scales linearly with τ at long times, indicating normal diffusive dynamics. The crossover to diffusive dynamics occurs on similar time and length scales in both liquid matrices and indicates the terminal relaxation of the matrix as particles escape from their local cages.

Next, we examine the dynamics of the glassy matrices (G0.20 and HSG1.00) with the same ϕ as the two supercooled liquids. The G0.20 and HSG1.00 glasses exhibit dynamics on intermediate time scales similar to the corresponding liquids (L0.35 and L1.05) but do not relax on long time scales [Fig. 3(a)]. Δr2 of the hard-sphere glass HSG1.00 exhibits a plateau on intermediate time scales, similar to the one observed for the repulsive liquid L1.05. The smaller plateau height in HSG1.00 relative to that in L1.05 indicates that thermal fluctuations in the repulsive liquid slightly increase the local cage size. Likewise, Δr2 of the attractive glass G0.20 resembles that of the liquid L0.35, exhibiting a small plateau followed by a power-law increase with time. The power-law exponent in this increasing region is β ≈ 0.11, smaller than the exponent β ≈ 0.32 for the corresponding L0.35 liquid matrix. This behavior indicates that matrix rearrangements are restricted on intermediate length and time scales due to the stronger attractions between particles in G0.20.

Two glassy matrices G0.35 and G1.05 can also be produced by compressing from ϕ = 0.610 to 0.635. This increase in ϕ leads to suppressed plateaus in Δr2 on intermediate time scales and prevents the matrix from fully relaxing on long time scales. Δr2 of G1.05 displays an intermediate-time plateau that is suppressed relative to the plateau for the liquid L1.05 [Fig. 3(b)] because of the higher matrix density of G1.05. The MSD of G0.35 exhibits a small shoulder at τ ≈ 10−1, qualitatively similar to but quantitatively lower than the one observed for L0.35, which arises from attractive bond formation. On time scales τ ≳ 102, the MSD of G0.35 exhibits a long-time plateau that contrasts with the extended power law observed on long time scales for Δr2 of L0.35. The plateau indicates that particles in G0.35 remain caged on long time scales, whereas the extended power-law in L0.35 indicates that bond rearrangement and cage escape occur on similar time scales. Thus, comparison of the MSDs for G0.35 and G1.05 with those of L0.35 and L1.05 reveals that increasing matrix ϕ results in caging-driven arrest.

Differences in packing fraction and interparticle interactions influence fluctuations and relaxations in the six matrices, leading to distinct dynamics. The dynamics of tracers of a critical relative size, δt = 0.35,17,18,40 within a slowly relaxing matrix are affected by both caging within matrix voids and matrix relaxations and fluctuations.17,18,40 We thus anticipate that differences in the matrix structure and dynamics will alter the motions of confined tracer particles, providing insight into the nature of the transport processes in these matrices.17,18

After initial ballistic motion, tracers of relative size δt = 0.35 relax differently on intermediate and long time scales in distinct matrices (Fig. 4). The tracers embedded in the L0.35 and L1.05 matrices exhibit subdiffusion on intermediate time scales 10−1τ ≲ 103. Diffusive dynamics are recovered on long time scales τ ≳ 103, with tracers in L1.05 exhibiting a higher diffusivity. The transition from ballistic motion to subdiffusive dynamics reflects the onset of caging by the matrix particles. The value of the MSD at the crossover to subdiffusive behavior is greater for L0.35 than in L1.05, indicating that tracers explore larger voids in L0.35.18 The logarithmic slope β of the tracer MSD in the subdiffusive regime (Δr2tβ), however, is smaller for L0.35 than for L1.05, reflecting smaller matrix fluctuations on these time scales (Fig. 3). The larger voids and slower tracer relaxations in L0.35 are due to subtle changes in the matrix structure arising from the strong attractive bonds between the matrix particles.18,41,42

FIG. 4.

Δr2 for tracer particles as a function of normalized time scale τ. (a) Δr2 for tracers in matrices with ϕ = 0.610. The arrows indicate the time regime during which Δr2G0.20Δr2L0.35. (b) Δr2 for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. The arrow indicates the time scale at which the behaviors of Δr2G0.35 and Δr2G1.05 qualitatively diverge. The black solid line indicates a power-law slope of one.

FIG. 4.

Δr2 for tracer particles as a function of normalized time scale τ. (a) Δr2 for tracers in matrices with ϕ = 0.610. The arrows indicate the time regime during which Δr2G0.20Δr2L0.35. (b) Δr2 for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. The arrow indicates the time scale at which the behaviors of Δr2G0.35 and Δr2G1.05 qualitatively diverge. The black solid line indicates a power-law slope of one.

Close modal

To obtain insight into the effects of temperature-induced matrix arrest on tracer dynamics, we compare tracer MSDs in L0.35 and L1.05 to those in G0.20 and HSG1.00. Tracer dynamics in HSG1.00 are subdiffusive for time scales 10−1τ ≲ 103 and diffusive for long time scales τ ≳ 103 [Fig. 4(a)]. In the subdiffusive regime, the tracer MSD of HSG1.00 has a slightly smaller slope and magnitude for a given τ relative to L1.05. Furthermore, diffusive dynamics are recovered later in time in HSG1.00 and the terminal diffusivity is smaller. These behaviors suggest that tracer dynamics in L1.05 are faster and less localized than in HSG1.00 and that tracer cage escape occurs on shorter time scales in L1.05 due to its liquidlike relaxations.

Tracer dynamics within the attractive matrices L0.35 and G0.20 exhibit distinct features not present in L1.05 and HSG1.00, reflecting the effects of attractive bonds between the matrix particles. The tracer MSDs Δr2 in L0.35 and G0.20 are subdiffusive and nearly identical for time scales 10−1τ ≤ 101 but diverge for τ > 101. In the subdiffusive regime (10−1τ ≤ 101), the slopes of the MSDs of L0.35 and G0.20 are smaller than those of L1.05 and HSG1.00. The MSD of tracers in L0.35 transitions from subdiffusive to diffusive scaling, becoming fully diffusive for τ > 103. By contrast, Δr2 of tracers in G0.20 scales subdiffusively, appearing to tend toward recovering diffusive behavior on time scales longer than those readily accessible in simulation. This extended subdiffusion arises from the increased role of matrix interparticle bonding, which reduces matrix fluctuations and increases the residence time of tracers in matrix cages.

A comparison of the matrix and tracer MSDs reveals dynamic coupling on intermediate time scales and suggests that these processes can facilitate tracer transport. Tracer dynamics in HSG1.00 are diffusive on time scales exceeding τ ≳ 103, even though the matrix itself does not relax. This behavior sharply contrasts with the lack of tracer diffusion in G0.20 even on much longer time scales of τ ∼ 105. The matrix Δr2 of HSG1.00 is an order of magnitude larger than that of G0.20 [Fig. 3(a)], indicating that intermediate time scale fluctuations are larger in the repulsive glass. Hence, thermal “cage rattling” in HSG1.00 allows tracers to escape and diffuse but is suppressed in G0.20 due to the presence of attractive bonds. Tracer dynamics begin to recover diffusive behavior in G0.20 on long time scales, when rare fluctuations in the matrix occur that allow tracers to escape. The long-time tracer dynamics in the arrested glasses G0.20 and HSG1.00 qualitatively differ from those in the completely frozen attractive and repulsive matrices examined in Ref. 18. Whereas tracers are localized in cages at long times in frozen matrices, the intermediate-time fluctuations of arrested glasses allow for long-time tracer relaxation and cage escape. This relaxation occurs on longer time scales within G0.20 than in HSG1.00 because G0.20 exhibits smaller matrix fluctuations and hence more closely approximates the environment encountered in completely frozen matrices.

A comparison of transport within attractive G0.20 and repulsive HSG1.00 glasses with the corresponding liquids (L0.35 and L1.05, respectively) reveals how changes in matrix relaxation processes due to temperature-induced vitrification affect tracer dynamics. Insights into how compression-induced vitrification influences tracer dynamics can be obtained by examining the glasses G0.35 and G1.05. Densification reduces the magnitude of the tracer MSDs in glasses at the onset of subdiffusion, relative to all other matrices, reflecting the smaller cages formed at higher ϕ. In addition, the logarithmic slopes of the tracer Δr2 in G0.35 and G1.05 for a given τ are smaller than those of the other matrices and diffusive dynamics are not recovered on long time scales. This reduction in tracer mobility on intermediate and long time scales appears to be primarily a trivial consequence of increased matrix density.

The nature of the interactions between matrix particles in G0.35 and G1.05 also affects tracer dynamics. The tracer Δr2 in G0.35 and G1.05 are nearly identical up to τ ≲ 102. This result is in sharp contrast with the marked differences in short- and intermediate time tracer dynamics in the lower-density glasses (G0.20 and HSG1.00). One possible explanation is that increasing matrix density may lead to more uniform cages in attractive and repulsive glasses that are accessible to tracers of this size on short time scales. For τ ≳ 102, however, the tracer MSDs remain subdiffusive but begin to diverge, with tracers in G0.35 exhibiting slower dynamics than in G1.05. This result suggests that tracer dynamics are still sensitive to differences in cage rattling in G1.05 and G0.35 [Fig. 3(b)] but only on longer time scales.

Tracers in all matrices exhibit subdiffusive behavior on intermediate time scales 10−1τ ≲ 103, suggesting that they are transiently caged by the matrix (Fig. 4). To further investigate the effect of matrix caging on tracer dynamics, we calculate the tracer CDP x̃12 as a function of the initial displacement magnitude r̃01 for lag times τ*. For a particle trapped in a harmonic well, x̃12 varies linearly with r̃01 with a slope of c = 0.5, indicating that the second displacement Δr12 is anticorrelated with the initial displacement Δr01.43 The slope c is the extent to which a particle is “dragged back” (this phenomenon is henceforth referred to as backdragging) as a fraction of its initial displacement r̃01.27 For matrix particles in glassy liquids, by contrast, x̃12 varies linearly with r̃01 for small displacements but deviates from this initial linear behavior for r̃01 beyond a length scale r̃cage identified as the characteristic cage size.27,28,44,45

For all tracers, the CDP linearly decreases for small displacements, indicating that they are caged by the matrices up to a length scale r̃cage (Fig. 5). The cage length scale depends on the matrix and varies between 0.3 and 0.8 (in units of σBB) (see the supplementary material). The slopes c, which reflect the extent of backdragging for tracers in each matrix, are ≳0.4 near the harmonic limit of c = 0.5 and are similar for all tracers for r̃01<r̃cage. The similar values of c indicates that tracers vibrate nearly harmonically within their cages, suggesting that collisions of the tracers with the matrix dominate dynamics up to r̃cage.

FIG. 5.

(a) Normalized CDP x̃12 as a function of previous displacement magnitude r̃01 for tracers in matrices with ϕ = 0.610. (b) Normalized CDP for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Shaded regions indicate statistical uncertainty. The black solid line indicates a slope of −0.5.

FIG. 5.

(a) Normalized CDP x̃12 as a function of previous displacement magnitude r̃01 for tracers in matrices with ϕ = 0.610. (b) Normalized CDP for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Shaded regions indicate statistical uncertainty. The black solid line indicates a slope of −0.5.

Close modal

Tracers with displacements larger than r̃cage escape their cages and rearrange into new positions. Steeper slopes for the tracer CDP x̃12 in the rearrangement regime indicate enhanced backdragging. For tracers within the repulsive liquid L1.05 and glass HSG1.00, x̃12 is approximately constant for r̃01>r̃cage, indicating that the extent of backdragging does not increase beyond r̃cage [Fig. 5(a)]. For the other matrices, by contrast, x̃12 continues to decrease with increasing r̃01 for r̃01>r̃cage but not as sharply as in the small-displacement regime. This behavior indicates that tracers are dragged back but not as far as predicted from linear extrapolation from the harmonic region.29 Similar behavior has also been observed in the rearrangement regime for probe particles in Laponite clay gels.46,47 Analyses of simple model systems show that such deviations from linearity are observed if the matrix is heterogeneous on the scale of the tracers or the relaxation rate of the probe particle is spatially dependent.47 Accordingly, the behavior of x̃12 for L0.35, G0.20, G0.35, and G1.05 indicates heterogeneity in tracer dynamics for length scales beyond r̃cage and within these matrices.18,41,42 A comparison of the slopes for different matrices suggests that backdragging outside of the local cage is enhanced by stronger matrix bonding (e.g., G0.20 vs L0.35). Furthermore, the fact that the x̃12 slopes are similar for G0.35 and G1.05 but different for the other matrices (L0.35, L1.05, and HSG1.00) indicates that matrix density, not matrix interparticle interactions, is the dominant factor controlling the extent of backdragging in the rearrangement regime.

To directly quantify tracer dynamical heterogeneity arising from differing matrix interparticle interactions, we calculated the dynamic susceptibility χ4(q, τ) at the wavevector magnitude q for which the peak in χ4 is maximized [Figs. 6(a) and 6(b)]. The dynamic susceptibility is the variance of tracer self-dynamics, defined as χ4(q,τ)Nt(Fs2(q,τ)[Fs(q,τ)]2). For all matrices, χ4 exhibits a peak whose location and width corresponds to the maximum and persistence of tracer dynamic heterogeneity, respectively. A comparison of the χ4 widths reveals that the persistence of the tracer dynamic heterogeneity varies across the different matrices, increasing such that L1.05 < HSG1.00 < G0.20. Based on the tracer CDPs for these matrices, we posit that the increase in tracer dynamic heterogeneity persistence may be associated with larger cages r̃cage (e.g., larger in HSG1.00 than in L1.05; see the supplementary material) and enhanced backdragging (e.g., in G0.20 vs HSG1.00) in the rearrangement regime.

FIG. 6.

Tracer dynamic susceptibility χ4 for BB for which the peak in χ4 is maximized. (a) χ4 for tracers in matrices with ϕ = 0.610. (b) χ4 for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. The wavevector magnitudes are BB = 2.00, 2.50, 1.60, 2.75, 3.00, and 3.00 for L0.35, L1.05, G0.20, HSG1.00, G0.35, and G1.05, respectively.

FIG. 6.

Tracer dynamic susceptibility χ4 for BB for which the peak in χ4 is maximized. (a) χ4 for tracers in matrices with ϕ = 0.610. (b) χ4 for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. The wavevector magnitudes are BB = 2.00, 2.50, 1.60, 2.75, 3.00, and 3.00 for L0.35, L1.05, G0.20, HSG1.00, G0.35, and G1.05, respectively.

Close modal

The susceptibility χ4 reveals that the dynamics of the embedded tracer particles are temporally heterogeneous. To quantify the structural and dynamical contributions to this heterogeneous behavior, we performed simulations in the isoconfigurational ensemble. Because this analysis allows a particle’s capacity for motion to be characterized for a given initial configuration, it enables the spatial distribution of dynamics in a given system to be linked to the system’s structure32 without requiring correlations to specific structural metrics (e.g., free volume, potential energy, etc.) to be established, which has proven to be extremely challenging.48,49 In this ensemble, total tracer dynamical fluctuations can be expressed as Δc,norm=σDP2+Δc,normiso.50–52 The first term σDP2=(DPi)21 is the structural variance, where ⟨⋯⟩ is the ensemble average over all nc,iso configurations, nt,iso trajectories, and Nt tracer particles. This variance is a measure of fluctuations in the tracer dynamic propensity DPi. The second term Δc,normiso=DP2,i(DPi)2 is the dynamical variance, where DP2,i=(ri(t)ri(0))4isoΔri2iso2. This variance quantifies the spread in DPi between different isoconfigurational trajectories.

In L0.35 and L1.05, the tracer structural variance σDP2 increases up to and peaks at τ ≈ 101, corresponding to the time at which tracers experience maximum structural determinism [Fig. 7(a)]. On longer time scales, σDP2 decreases and reaches a value near zero at τ ≈ 103, approximately the lag time at which tracer dynamics become diffusive in the MSD Δr2 (Fig. 4). The structural variance σDP2 is larger and attains a maximum at a larger lag time τ for tracers within L0.35 compared to L1.05. Strictly positive values of σDP2 arise from particle-to-particle variations in dynamic propensity. Hence, the larger values of σDP2 indicate that tracer relaxation times are more broadly distributed in L0.35, which is a consequence of the structural heterogeneity arising from strong interparticle bonds in this matrix. The existence of greater structural heterogeneity in L0.35 than in L1.05 is evidenced by the shorter first peak in the matrix static structure factor S(q), which indicates reduced translational ordering. Similar evidence can be found by inspecting the matrix radial distribution function g(r) (see the supplementary material).

FIG. 7.

Tracer structural variance σDP2 of the dynamic propensity DPi distribution. (a) σDP2 for tracers in matrices with ϕ = 0.610. (b) σDP2 for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Shaded regions indicate statistical uncertainty.

FIG. 7.

Tracer structural variance σDP2 of the dynamic propensity DPi distribution. (a) σDP2 for tracers in matrices with ϕ = 0.610. (b) σDP2 for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Shaded regions indicate statistical uncertainty.

Close modal

To understand the effects of matrix arrest on structural determinism in tracer dynamics, we compare the structural variances in L0.35 and L1.05 to those in G0.20 and HSG1.00 [Fig. 7(a)]. The tracer σDP2 in HSG1.00 exhibits a larger peak at a larger τ than in L1.05. The structure of the nearly arrested HSG1.00 matrix obstructs tracer rearrangement, leading to a broader DPi distribution. Surprisingly, σDP2 is greater in attractive than repulsive arrested matrices σDP,HSG1.002<σDP,G0.202, and the disparity in σDP2 between L0.35 and G0.20 is much larger than the disparity between L1.05 and HSG1.00. These observations suggest that arrest by bonding in attractive matrices results in local environments that are highly structurally heterogeneous compared to those in the repulsive matrices.

In the higher-density matrices (G0.35 and G1.05), σDP2 increases steeply for tracers for τ ≲ 103, appearing to plateau or even slightly decrease on longer time scales [Fig. 7(b)]. On time scales exceeding 103, σDP2 increases very gradually or decreases slightly for tracers within G0.35 and G1.05, respectively. The strong initial increase in σDP2 indicates heterogeneous local environments and a broad distribution of tracer relaxation times. This increase is similar to that observed in the lower-density, strongly arrested glass G0.20. The common feature of G0.20, G0.35, and G1.05 is that matrix rattling is suppressed on intermediate to long time scales, indicated by the low values of the matrix MSD (Fig. 3). Thus, reduced matrix fluctuations appear to lead to an increase in structural determinism in the tracer dynamics.

Finally, to characterize the relative importance of the structure for tracer dynamics, we examine the relative contribution of σDP2 to total tracer fluctuations within each matrix via the structural ratio Rc=σDP2/Δc,norm, where Rc is the fractional contribution of isoconfigurational fluctuations due to the structure.50,51 For tracers in all matrices, Rc is approximately constant on short time scales τ ≲ 101 but varies depending on the matrix [Figs. 8(a) and 8(b)]. This behavior suggests that the relative contributions from the matrix structure are fixed on these time scales but depend on matrix density and interactions between matrix particles. The values of Rc for matrix particles are larger than those of the corresponding tracers (see the supplementary material) for all matrices, indicating that structural contributions are more important for the matrix dynamics than for tracer dynamics. For τ ≳ 102, Rc for tracers in L1.05, L0.35, and HSG1.00 decay toward zero as tracer dynamics become diffusive (Fig. 4) and χ4 peaks [Fig. 6(a)]. Collectively, these observations suggest that relatively large dynamical matrix fluctuations in these matrices allow tracers to escape their local cage so that the matrix structure no longer affects tracer dynamics. By contrast, the values of Rc for G0.20, G0.35, and G1.05 continue to increase for τ ≳ 102. The absence of tracer diffusion within these matrices on these time scales suggests that tracers are partially or fully localized, and thus their dynamics are strongly influenced by the matrix structure. This interpretation is consistent with the nearly constant dynamical variances Δc,normiso for G0.20, G0.35, and G1.05 on intermediate to long time scales (see the supplementary material). Sufficiently strong vitrification, whether through attractive bonding (G0.20) or increased matrix density (G0.35 and G1.05), reduces dynamical matrix fluctuations and thereby hinders the ability of tracers to escape their cages. As a result, tracer dynamics in highly vitrified matrices are more strongly influenced by the matrix structure.

FIG. 8.

Fraction of tracer fluctuations due to initial structure Rc. (a) Rc for tracers in matrices with ϕ = 0.610. (b) Rc for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Shaded regions indicate statistical uncertainty.

FIG. 8.

Fraction of tracer fluctuations due to initial structure Rc. (a) Rc for tracers in matrices with ϕ = 0.610. (b) Rc for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Shaded regions indicate statistical uncertainty.

Close modal

To characterize tracer exploration within each matrix, we analyze the scaling of the tracer trajectory radius of gyration Rg as a function of its mass M. A comparison of the trajectory shapes for tracers within attractive and repulsive matrices reveal that the nature of interactions between matrix particles influences tracer exploration (Fig. 9). For tracers within L1.05 and HSG1.00, M scales approximately as a power law for large Rg (Rg > 0.5, long times). The logarithmic slope df in this region defines the fractal dimension, which is approximately 2.0 in both matrices. This value corresponds to the limit of free diffusion that is expected for tracers at long times (large Rg) in all matrices in which tracer dynamics are ergodic. For tracers within L0.35 and G0.20, by contrast, the instantaneous value of df is larger than 2.0. This result indicates that the trajectories in attractive matrices are more compact than those within repulsive matrices, with a fractal dimension that approaches that of a geometric solid (i.e., df = 3). The nearly indistinguishable mass scaling of the arrested and liquid matrices suggests that it is matrix interparticle interactions and not dynamical arrest that leads to the difference between attractive and repulsive matrices.

FIG. 9.

Mass M of tracer trajectories as a function of trajectory radius of gyration Rg as calculated by box-counting (a) M for tracers in matrices with ϕ = 0.610. (b) M for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Black dotted and dashed lines indicate fractal scalings of 3 and 2, respectively.

FIG. 9.

Mass M of tracer trajectories as a function of trajectory radius of gyration Rg as calculated by box-counting (a) M for tracers in matrices with ϕ = 0.610. (b) M for tracers in glasses with ϕ = 0.635 and in the corresponding liquids. Black dotted and dashed lines indicate fractal scalings of 3 and 2, respectively.

Close modal

Finally, we examine the role of matrix density on the tracer trajectory shape by comparing the fractal scaling of tracer trajectories in G0.35 and G1.05 to L0.35 and L1.05 [Fig. 9(b)]. The slopes in the glasses are larger than those in the corresponding liquid, illustrating that tracer trajectories within the glasses have larger fractal dimensions than those within the liquids. Thus, tracer trajectories in G0.35 and G1.05 are less extended in space, indicating that tracers explore cages for longer periods of time and rearrange over smaller distances as matrix density is increased.

We used event-driven MD to investigate the dynamics of dilute, hard-sphere tracers in attractive and repulsive liquid matrices with similar long-time relaxations and in analogous attractive and repulsive glass matrices prepared via thermal quenching or compression. A comparison of the matrix and tracer MSDs revealed the effects of matrix dynamics on tracer dynamics. Although tracers within glasses were less mobile than tracers within the corresponding liquids, matrix arrest was insufficient to guarantee tracer localization. Through cage analysis, we determined whether tracers are caged for small displacements and characterized the heterogeneity of tracer cage rearrangements for large displacements. This analysis revealed that increasing matrix density ϕ from 0.610 to 0.635 (G0.35, G1.05) or enhancing attractions (G0.20, L0.35) increased the tracer dynamic heterogeneity within these matrices relative to the repulsive matrices with ϕ = 0.610 (L1.05, HSG1.00). The tracer dynamic susceptibility revealed that tracer dynamics were spatiotemporally heterogeneous, as also shown through cage analysis. By performing simulations in the isoconfigurational ensemble and calculating the dynamic propensity, we quantified the structurally-induced heterogeneity of tracer dynamics and the extent to which tracer dynamics were determined by the matrix structure. This analysis revealed that strong arrest of the matrix, driven by attractive bonding or high density, enhanced structural determinism. Finally, the mass scaling of tracer trajectories revealed that increasing matrix attractions or matrix density leads to more compact tracer trajectories. These results collectively reveal how the spatial and temporal heterogeneity in matrices is reflected in the dynamics of embedded tracers.

Our simulations also demonstrate that tracers are able to diffuse on long time scales through glass matrices if the matrix fluctuations are sufficiently large. This result has interesting implications for understanding the ability of tracers to penetrate dense, slowly relaxing matrices, suggesting that fluctuations above a critical size can facilitate transport even in matrices that do not fully relax on long time scales. The findings from our study also motivate future work in a number of different directions. Whereas our investigation focused on understanding equilibrium tracer dynamics, penetrant transport in most practical settings is driven by a chemical potential gradient and hence occurs under nonequilibrium conditions. Although much work has been done toward understanding nonequilibrium transport through rigid matrices,53–55 it remains unclear how these processes are influenced by structural fluctuations and slow matrix relaxations. Additionally, fluctuations within the matrices studied here are isotropic due to the bulk nature of the samples imposed through the use of periodic boundary conditions. Experimental studies of supercooled liquids and glasses show, however, that fluctuations in these systems can become anisotropic by imposing different boundary conditions.56,57 This scenario has been encountered, for example, when examining the dynamics of confined supercooled liquids in porous media58,59 and glasses prepared through vapor deposition onto surfaces.60 An intriguing future line of inquiry would be to investigate how anisotropic fluctuations in these systems affect the dynamic coupling of the tracer and matrix particles. Finally, we investigated small particles in the isolated tracer limit. For asymmetric hard-sphere binary mixtures, both reentrant melting of the large-particle glass and vitrification of the small particles were observed upon increasing small-particle density.61,62 How matrix interactions affect and, in turn, are affected by tracers at high concentrations has not yet been studied. We anticipate that these outstanding questions can be addressed using computational approaches similar to those employed in this study.

See the supplementary material for a table of matrices and simulation parameters, tracer non-Gaussian parameters, matrix self-intermediate scattering functions, estimates of the tracer caging length scales, matrix structure factors and radial distribution functions, matrix Rc, and tracer dynamical variances.

The authors acknowledge support from the Welch Foundation (Grant Nos. E-1869 to J.C.C. and E-1882 to J.C.P.) and the National Science Foundation (Grant No. CBET-1705968).

1.
N. A.
Hadjiev
and
B. G.
Amsden
,
J. Controlled Release
199
,
10
(
2015
).
2.
J.
Li
and
D. J.
Mooney
,
Nat. Rev. Mater.
1
,
16071
(
2016
).
3.
R. W.
Baker
,
Ind. Eng. Chem. Res.
41
,
1393
(
2002
).
4.
D. R.
Paul
and
L. M.
Robeson
,
Polymer
49
,
3187
(
2008
).
5.
R. J.
Ellis
,
Curr. Opin. Struct. Biol.
11
,
114
(
2001
).
6.
S.
Nakano
,
D.
Miyoshi
, and
N.
Sugimoto
,
Chem. Rev.
114
,
2733
(
2013
).
7.
K.
Regan
,
D.
Wulstein
,
H.
Rasmussen
,
R.
McGorty
, and
R. M.
Robertson-Anderson
,
Soft Matter
15
,
1200
(
2019
).
8.
F.
Höfling
,
T.
Franosch
, and
E.
Frey
,
Phys. Rev. Lett.
96
,
165901
(
2006
).
9.
Y.
Jin
and
P.
Charbonneau
,
Phys. Rev. E
91
,
042313
(
2015
).
10.
M. H.
Cohen
and
D.
Turnbull
,
J. Chem. Phys.
31
,
1164
(
1959
).
11.
F.
Müller-Plathe
,
J. Chem. Phys.
94
,
3192
(
1991
).
12.
A. A.
Gusev
and
U. W.
Suter
,
J. Chem. Phys.
99
,
2228
(
1993
).
13.
J.
Hwang
,
J.
Kim
, and
B. J.
Sung
,
Phys. Rev. E
94
,
022614
(
2016
).
14.
N.
Alcázar-Cano
and
R.
Delgado-Buscalioni
,
Soft Matter
14
,
9937
(
2018
).
15.
S. K.
Schnyder
,
M.
Spanner
,
F.
Höfling
,
T.
Franosch
, and
J.
Horbach
,
Soft Matter
11
,
701
(
2015
).
16.
T. O. E.
Skinner
,
S. K.
Schnyder
,
D. G. A. L.
Aarts
,
J.
Horbach
, and
R. P. A.
Dullens
,
Phys. Rev. Lett.
111
,
128301
(
2013
).
17.
T.
Sentjabrskaja
,
E.
Zaccarelli
,
C.
De Michele
,
F.
Sciortino
,
P.
Tartaglia
,
T.
Voigtmann
,
S. U.
Egelhaaf
, and
M.
Laurati
,
Nat. Commun.
7
,
11133
(
2016
).
18.
R. C.
Roberts
,
R.
Poling-Skutvik
,
J. C.
Palmer
, and
J. C.
Conrad
,
J. Phys. Chem. Lett.
9
,
3008
(
2018
).
19.
W. K.
Kegel
and
A.
van Blaaderen
,
Science
287
,
290
(
2000
).
20.
E. R.
Weeks
,
J. C.
Crocker
,
A. C.
Levitt
,
A.
Schofield
, and
D. A.
Weitz
,
Science
287
,
627
(
2000
).
21.
E.
Zaccarelli
,
G.
Foffi
,
K. A.
Dawson
,
S. V.
Buldyrev
,
F.
Sciortino
, and
P.
Tartaglia
,
Phys. Rev. E
66
,
041402
(
2002
).
22.
K. N.
Pham
,
A. M.
Puertas
,
J.
Bergenholtz
,
S. U.
Egelhaaf
,
A.
Moussa
d
ı,
P. N.
Pusey
,
A. B.
Schofield
,
M. E.
Cates
,
M.
Fuchs
, and
W. C. K.
Poon
,
Science
296
,
104
(
2002
).
23.
E.
Zaccarelli
and
W. C. K.
Poon
,
Proc. Natl. Acad. Sci. U. S. A.
106
,
15203
(
2009
).
24.
K. A.
Dawson
,
Curr. Opin. Colloid Interface Sci.
7
,
218
(
2002
).
25.
E.
Zaccarelli
,
H.
Löwen
,
P. P. F.
Wessels
,
F.
Sciortino
,
P.
Tartaglia
, and
C. N.
Likos
,
Phys. Rev. Lett.
92
,
225703
(
2004
).
26.
E.
Zaccarelli
,
F.
Sciortino
, and
P.
Tartaglia
,
J. Phys.: Condens. Matter
16
,
S4849
(
2004
).
27.
B.
Doliwa
and
A.
Heuer
,
Phys. Rev. Lett.
80
,
4915
(
1998
).
28.
B.
Doliwa
and
A.
Heuer
,
J. Phys.: Condens. Matter
11
,
A277
(
1999
).
29.
E. R.
Weeks
and
D. A.
Weitz
,
Chem. Phys.
284
,
361
(
2002
).
30.
W.
Kob
,
C.
Donati
,
S. J.
Plimpton
,
P. H.
Poole
, and
S. C.
Glotzer
,
Phys. Rev. Lett.
79
,
2827
(
1997
).
31.
M. S.
Shell
,
P. G.
Debenedetti
, and
F. H.
Stillinger
,
J. Phys.: Condens. Matter
17
,
S4035
(
2005
).
32.
A.
Widmer-Cooper
and
P.
Harrowell
,
J. Chem. Phys.
126
,
154503
(
2007
).
33.
A.
Widmer-Cooper
,
P.
Harrowell
, and
H.
Fynewever
,
Phys. Rev. Lett.
93
,
135701
(
2004
).
34.
M.
Fitzner
,
G. C.
Sosso
,
S. J.
Cox
, and
A.
Michaelides
,
Proc. Natl. Acad. Sci. U. S. A.
116
,
2009
(
2019
).
35.
G. C.
Sosso
,
J.
Colombo
,
J.
Behler
,
E.
Del Gado
, and
M.
Bernasconi
,
J. Phys. Chem. B
118
,
13621
(
2014
).
36.
A. A.
Saberi
,
Phys. Rev. E
84
,
021113
(
2011
).
37.
D. A.
Russell
,
J. D.
Hanson
, and
E.
Ott
,
Phys. Rev. Lett.
45
,
1175
(
1980
).
38.
D. N.
Perera
and
P.
Harrowell
,
Phys. Rev. Lett.
81
,
120
(
1998
).
39.
K.
Dawson
,
G.
Foffi
,
M.
Fuchs
,
W.
Götze
,
F.
Sciortino
,
M.
Sperl
,
P.
Tartaglia
,
T.
Voigtmann
, and
E.
Zaccarelli
,
Phys. Rev. E
63
,
011401
(
2000
).
40.
R.
Poling-Skutvik
,
R. C.
Roberts
,
A. H.
Slim
,
S.
Narayanan
,
R.
Krishnamoorti
,
J. C.
Palmer
, and
J. C.
Conrad
,
J. Phys. Chem. Lett.
10
,
1784
(
2019
).
41.
W. P.
Krekelberg
,
V.
Ganesan
, and
T. M.
Truskett
,
J. Phys. Chem. B
110
,
5166
(
2006
).
42.
W. P.
Krekelberg
,
J.
Mittal
,
V.
Ganesan
, and
T. M.
Truskett
,
J. Chem. Phys.
127
,
044502
(
2007
).
43.
A.
Heuer
,
M.
Kunow
,
M.
Vogel
, and
R. D.
Banhatti
,
Phys. Rev. B
66
,
224201
(
2002
).
44.
E. R.
Weeks
and
D. A.
Weitz
,
Phys. Rev. Lett.
89
,
095704
(
2002
).
45.
Y.
Gebremichael
,
M.
Vogel
, and
S. C.
Glotzer
,
J. Chem. Phys.
120
,
4415
(
2004
).
46.
F. K.
Oppong
,
L.
Rubatat
,
B. J.
Frisken
,
A. E.
Bailey
, and
J. R.
de Bruyn
,
Phys. Rev. E
73
,
041405
(
2006
).
47.
J. P.
Rich
,
G. H.
McKinley
, and
P. S.
Doyle
,
J. Rheol.
55
,
273
(
2011
).
48.
A.
Widmer-Cooper
and
P.
Harrowell
,
J. Phys.: Condens. Matter
17
,
S4025
(
2005
).
49.
A.
Widmer-Cooper
and
P.
Harrowell
,
J. Non-Cryst. Solids
352
,
5098
(
2006
).
50.
L.
Berthier
and
R. L.
Jack
,
Phys. Rev. E
76
,
041509
(
2007
).
51.
H.
Tong
and
H.
Tanaka
,
Phys. Rev. X
8
,
011041
(
2018
).
52.
M.
Leocmach
and
H.
Tanaka
,
Nat. Commun.
3
,
974
(
2012
).
53.
S.
He
,
J. C.
Palmer
, and
G.
Qin
,
Microporous Mesoporous Mater.
249
,
88
(
2017
).
54.
K. E.
Gubbins
,
Y.-C.
Liu
,
J. D.
Moore
, and
J. C.
Palmer
,
Phys. Chem. Chem. Phys.
13
,
58
(
2011
).
55.
E. J.
Maginn
,
A. T.
Bell
, and
D. N.
Theodorou
,
J. Phys. Chem.
97
,
4173
(
1993
).
56.
C. R.
Nugent
,
K. V.
Edmond
,
H. N.
Patel
, and
E. R.
Weeks
,
Phys. Rev. Lett.
99
,
025702
(
2007
).
57.
G. L.
Hunter
,
K. V.
Edmond
, and
E. R.
Weeks
,
Phys. Rev. Lett.
112
,
218302
(
2014
).
58.
A.
Faraone
,
L.
Liu
,
C.-Y.
Mou
,
C.-W.
Yen
, and
S.-H.
Chen
,
J. Chem. Phys.
121
,
10843
(
2004
).
60.
M. D.
Ediger
,
J. Chem. Phys.
147
,
210901
(
2017
).
61.
T.
Sentjabrskaja
,
A. R.
Jacob
,
S. U.
Egelhaaf
,
G.
Petekidis
,
T.
Voigtmann
, and
M.
Laurati
,
Soft Matter
15
,
2232
(
2019
).
62.
M.
Laurati
,
T.
Sentjabrskaja
,
J.
Ruiz-Franco
,
S. U.
Egelhaaf
, and
E.
Zaccarelli
,
Phys. Chem. Chem. Phys.
20
,
18630
(
2018
).

Supplementary Material