Tracer transport in attractive and repulsive supercooled liquids and glasses

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 membranes^10–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-like¹⁴ 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 dynamics^19,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 u₀ = 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 k_B = 1 and σ_BB, u₀/k_B, 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 τ = tD₀/D_ref, where t is the nominal simulation time, $D0=σBBkBT/m$ is the thermal diffusivity at temperature T, and D_ref is the reference value of D₀ at T = 1.²³ 

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 L_T or G_T (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 D_i/D₀: a high-temperature repulsive glassy liquid (L_1.05) and a low-temperature attractive glassy liquid (L_0.35), where D_i is the nominal diffusion coefficient.¹⁸ Two glasses were also prepared at the same temperatures (G_1.05, G_0.35) and increased matrix volume fraction ϕ = 0.635. Similarly, a glass (G_0.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 (HSG_1.00), which is equivalent to a square-well glass in the high-temperature limit (Fig. 2).

We embedded N_t = 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 (G_0.20 and HSG_1.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 (G_0.35 and G_1.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.²⁶ For the HSG_1.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 ≈ 10⁵).²³ Statistical properties were calculated by averaging over trajectories computed for n_s = 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 $Δr→01=r→(Δt)−r→(0)$ over a time interval Δt on which $Δr→01$ is comparable to the characteristic cage size, where $r→$ is the particle’s position vector. Because the neighbor cage restricts further motion along $Δr→01$⁠, the displacement over the next time interval $Δr→12=r→(2Δt)−r→(Δt)$ should, on average, be anticorrelated. The displacement $Δr→12$ can be projected onto the unit vector of the preceding displacement $Δr→01$⁠, yielding the caging displacement projection (CDP) $x12≡Δr→12⋅Δr→01|Δr→01|$⁠, which is parallel to $Δr→01$⁠. For a caged particle, the ensemble-averaged CDP ⟨x₁₂⟩ is negative and its magnitude grows with $|Δr→01|$⁠. In hard-sphere supercooled liquids, $⟨x12⟩=−c|Δr→01|$ for small displacements where the cage has not been broken. Larger magnitudes of the proportionality constant c indicate greater displacement memory during caging.²⁷ For larger, cage-breaking displacements, by contrast, ⟨x₁₂⟩ is independent of $|Δr→01|$⁠.²⁹ 

In analyzing tracer dynamics, we first compute the non-Gaussian parameter for particle displacements α₂(Δt) (see the supplementary material). The maximum in α₂(Δt) 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 $Δr→01$ and $Δr→12$⁠.²⁹ The magnitudes of tracer displacements at τ^* vary across different matrices. To account for this variation, we report normalized quantities $x̃12≡C−1⟨x12⟩$ and $r̃01≡C−1|Δr→01|$⁠, where $C$ is the square-root of the tracer mean-square displacement (MSD) Δr² 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.³² For each matrix system, we analyzed n_c,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 n_c,iso configurations was used to initialize n_t,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)=⟨(r→i(t)−r→i(0))2Δri2⟩iso$⁠, 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 $⟨Δri2⟩iso$⁠, DP_i 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 R_g. 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=1nng∑i=1ncg(x→i−x→avg)2$⁠, where n_cg is the number of coarse-grained points in the trajectory, $x→i$ is the position of the ith coarse-grained point, and $x→avg=1ncg∑i=1ncgx→i$ is the mean position.

We first characterize the dynamics of the different glassy matrices through the ensemble-averaged mean-square displacement (MSD) Δr² (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,τ)=⟨1Ni∑k=1Ni⁡exp[−jq→⋅(r→k(τ)−r→k(0))]⟩$⁠, where $q=|q→|$ is the wavevector magnitude, $j=−1$⁠, N_i 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.

The MSD of the large matrix particles in L_1.05 is approximately constant on lag times τ ≈ 10⁰–10³, indicating interparticle caging [Fig. 3(a)]. By contrast, Δr² for L_0.35 exhibits a small plateau at τ ≈ 10⁻¹ followed by an increasing, subdiffusive power-law Δr² ≈ τ^β that extends to τ ≈ 10³, 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 τ.²¹ The smaller values of Δr² indicate that particles in L_0.35 are more localized than those in L_1.05, likely due to the formation of interparticle bonds.³⁹ 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 L_0.35 and L_1.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 (G_0.20 and HSG_1.00) with the same ϕ as the two supercooled liquids. The G_0.20 and HSG_1.00 glasses exhibit dynamics on intermediate time scales similar to the corresponding liquids (L_0.35 and L_1.05) but do not relax on long time scales [Fig. 3(a)]. Δr² of the hard-sphere glass HSG_1.00 exhibits a plateau on intermediate time scales, similar to the one observed for the repulsive liquid L_1.05. The smaller plateau height in HSG_1.00 relative to that in L_1.05 indicates that thermal fluctuations in the repulsive liquid slightly increase the local cage size. Likewise, Δr² of the attractive glass G_0.20 resembles that of the liquid L_0.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 L_0.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 G_0.20.

Two glassy matrices G_0.35 and G_1.05 can also be produced by compressing from ϕ = 0.610 to 0.635. This increase in ϕ leads to suppressed plateaus in Δr² on intermediate time scales and prevents the matrix from fully relaxing on long time scales. Δr² of G_1.05 displays an intermediate-time plateau that is suppressed relative to the plateau for the liquid L_1.05 [Fig. 3(b)] because of the higher matrix density of G_1.05. The MSD of G_0.35 exhibits a small shoulder at τ ≈ 10⁻¹, qualitatively similar to but quantitatively lower than the one observed for L_0.35, which arises from attractive bond formation. On time scales τ ≳ 10², the MSD of G_0.35 exhibits a long-time plateau that contrasts with the extended power law observed on long time scales for Δr² of L_0.35. The plateau indicates that particles in G_0.35 remain caged on long time scales, whereas the extended power-law in L_0.35 indicates that bond rearrangement and cage escape occur on similar time scales. Thus, comparison of the MSDs for G_0.35 and G_1.05 with those of L_0.35 and L_1.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 L_0.35 and L_1.05 matrices exhibit subdiffusion on intermediate time scales 10⁻¹ ≲ τ ≲ 10³. Diffusive dynamics are recovered on long time scales τ ≳ 10³, with tracers in L_1.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 L_0.35 than in L_1.05, indicating that tracers explore larger voids in L_0.35.¹⁸ The logarithmic slope β of the tracer MSD in the subdiffusive regime (Δr² ∼ t^β), however, is smaller for L_0.35 than for L_1.05, reflecting smaller matrix fluctuations on these time scales (Fig. 3). The larger voids and slower tracer relaxations in L_0.35 are due to subtle changes in the matrix structure arising from the strong attractive bonds between the matrix particles.^18,41,42

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

Tracer dynamics within the attractive matrices L_0.35 and G_0.20 exhibit distinct features not present in L_1.05 and HSG_1.00, reflecting the effects of attractive bonds between the matrix particles. The tracer MSDs Δr² in L_0.35 and G_0.20 are subdiffusive and nearly identical for time scales 10⁻¹ ≤ τ ≤ 10¹ but diverge for τ > 10¹. In the subdiffusive regime (10⁻¹ ≤ τ ≤ 10¹), the slopes of the MSDs of L_0.35 and G_0.20 are smaller than those of L_1.05 and HSG_1.00. The MSD of tracers in L_0.35 transitions from subdiffusive to diffusive scaling, becoming fully diffusive for τ > 10³. By contrast, Δr² of tracers in G_0.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 HSG_1.00 are diffusive on time scales exceeding τ ≳ 10³, even though the matrix itself does not relax. This behavior sharply contrasts with the lack of tracer diffusion in G_0.20 even on much longer time scales of τ ∼ 10⁵. The matrix Δr² of HSG_1.00 is an order of magnitude larger than that of G_0.20 [Fig. 3(a)], indicating that intermediate time scale fluctuations are larger in the repulsive glass. Hence, thermal “cage rattling” in HSG_1.00 allows tracers to escape and diffuse but is suppressed in G_0.20 due to the presence of attractive bonds. Tracer dynamics begin to recover diffusive behavior in G_0.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 G_0.20 and HSG_1.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 G_0.20 than in HSG_1.00 because G_0.20 exhibits smaller matrix fluctuations and hence more closely approximates the environment encountered in completely frozen matrices.

A comparison of transport within attractive G_0.20 and repulsive HSG_1.00 glasses with the corresponding liquids (L_0.35 and L_1.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 G_0.35 and G_1.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 Δr² in G_0.35 and G_1.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 G_0.35 and G_1.05 also affects tracer dynamics. The tracer Δr² in G_0.35 and G_1.05 are nearly identical up to τ ≲ 10². This result is in sharp contrast with the marked differences in short- and intermediate time tracer dynamics in the lower-density glasses (G_0.20 and HSG_1.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 τ ≳ 10², however, the tracer MSDs remain subdiffusive but begin to diverge, with tracers in G_0.35 exhibiting slower dynamics than in G_1.05. This result suggests that tracer dynamics are still sensitive to differences in cage rattling in G_1.05 and G_0.35 [Fig. 3(b)] but only on longer time scales.

Tracers in all matrices exhibit subdiffusive behavior on intermediate time scales 10⁻¹ ≲ τ ≲ 10³, 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 $Δr→12$ is anticorrelated with the initial displacement $Δr→01$⁠.⁴³ 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$⁠.²⁷ 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$⁠.

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 L_1.05 and glass HSG_1.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.²⁹ 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.⁴⁷ Accordingly, the behavior of $x̃12$ for L_0.35, G_0.20, G_0.35, and G_1.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., G_0.20 vs L_0.35). Furthermore, the fact that the $x̃12$ slopes are similar for G_0.35 and G_1.05 but different for the other matrices (L_0.35, L_1.05, and HSG_1.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 χ₄(q, τ) at the wavevector magnitude q for which the peak in χ₄ 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, χ₄ exhibits a peak whose location and width corresponds to the maximum and persistence of tracer dynamic heterogeneity, respectively. A comparison of the χ₄ widths reveals that the persistence of the tracer dynamic heterogeneity varies across the different matrices, increasing such that L_1.05 < HSG_1.00 < G_0.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 HSG_1.00 than in L_1.05; see the supplementary material) and enhanced backdragging (e.g., in G_0.20 vs HSG_1.00) in the rearrangement regime.

The susceptibility χ₄ 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 structure³² 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)2⟩−1$ is the structural variance, where ⟨⋯⟩ is the ensemble average over all n_c,iso configurations, n_t,iso trajectories, and N_t tracer particles. This variance is a measure of fluctuations in the tracer dynamic propensity DP_i. The second term $Δc,normiso=⟨DP2,i−(DPi)2⟩$ is the dynamical variance, where $DP2,i=⟨(r→i(t)−r→i(0))4⟩iso⟨Δri2⟩iso2$⁠. This variance quantifies the spread in DP_i between different isoconfigurational trajectories.

In L_0.35 and L_1.05, the tracer structural variance $σDP2$ increases up to and peaks at τ ≈ 10¹, 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 τ ≈ 10³, approximately the lag time at which tracer dynamics become diffusive in the MSD Δr² (Fig. 4). The structural variance $σDP2$ is larger and attains a maximum at a larger lag time τ for tracers within L_0.35 compared to L_1.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 L_0.35, which is a consequence of the structural heterogeneity arising from strong interparticle bonds in this matrix. The existence of greater structural heterogeneity in L_0.35 than in L_1.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).

To understand the effects of matrix arrest on structural determinism in tracer dynamics, we compare the structural variances in L_0.35 and L_1.05 to those in G_0.20 and HSG_1.00 [Fig. 7(a)]. The tracer $σDP2$ in HSG_1.00 exhibits a larger peak at a larger τ than in L_1.05. The structure of the nearly arrested HSG_1.00 matrix obstructs tracer rearrangement, leading to a broader DP_i distribution. Surprisingly, $σDP2$ is greater in attractive than repulsive arrested matrices $σDP,HSG1.002<σDP,G0.202$⁠, and the disparity in $σDP2$ between L_0.35 and G_0.20 is much larger than the disparity between L_1.05 and HSG_1.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 (G_0.35 and G_1.05), $σDP2$ increases steeply for tracers for τ ≲ 10³, appearing to plateau or even slightly decrease on longer time scales [Fig. 7(b)]. On time scales exceeding 10³, $σDP2$ increases very gradually or decreases slightly for tracers within G_0.35 and G_1.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 G_0.20. The common feature of G_0.20, G_0.35, and G_1.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 R_c is the fractional contribution of isoconfigurational fluctuations due to the structure.^50,51 For tracers in all matrices, R_c is approximately constant on short time scales τ ≲ 10¹ 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 R_c 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 τ ≳ 10², R_c for tracers in L_1.05, L_0.35, and HSG_1.00 decay toward zero as tracer dynamics become diffusive (Fig. 4) and χ₄ 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 R_c for G_0.20, G_0.35, and G_1.05 continue to increase for τ ≳ 10². 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 G_0.20, G_0.35, and G_1.05 on intermediate to long time scales (see the supplementary material). Sufficiently strong vitrification, whether through attractive bonding (G_0.20) or increased matrix density (G_0.35 and G_1.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.

To characterize tracer exploration within each matrix, we analyze the scaling of the tracer trajectory radius of gyration R_g 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 L_1.05 and HSG_1.00, M scales approximately as a power law for large R_g (R_g > 0.5, long times). The logarithmic slope d_f 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 R_g) in all matrices in which tracer dynamics are ergodic. For tracers within L_0.35 and G_0.20, by contrast, the instantaneous value of d_f 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., d_f = 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.

Finally, we examine the role of matrix density on the tracer trajectory shape by comparing the fractal scaling of tracer trajectories in G_0.35 and G_1.05 to L_0.35 and L_1.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 G_0.35 and G_1.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 (G_0.35, G_1.05) or enhancing attractions (G_0.20, L_0.35) increased the tracer dynamic heterogeneity within these matrices relative to the repulsive matrices with ϕ = 0.610 (L_1.05, HSG_1.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 media^58,59 and glasses prepared through vapor deposition onto surfaces.⁶⁰ 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 R_c, 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).