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.

## I. INTRODUCTION

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^{14} 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.

## II. METHODS

### A. Model systems

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*_{0} = 1 and width Δ_{ij} = 0.03(*σ*_{ij} + Δ_{ij}) for each pair type *i*, *j* ∈ A, B, where $\sigma ij=12(\sigma ii+\sigma jj)$. In the discussion to follow, we adopt conventional simulation units in which Boltzmann’s constant *k*_{B} = 1 and *σ*_{BB}, *u*_{0}/*k*_{B}, and $\sigma 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*_{0}/*D*_{ref}, where *t* is the nominal simulation time, $D0=\sigma BBkBT/m$ is the thermal diffusivity at temperature *T*, and *D*_{ref} is the reference value of *D*_{0} at *T* = 1.^{23}

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*_{0}: 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.^{18} 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.^{26} 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^{5}).^{23} 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.

### B. Cage analysis

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 $\Delta r\u219201=r\u2192(\Delta t)\u2212r\u2192(0)$ over a time interval Δ*t* on which $\Delta r\u219201$ is comparable to the characteristic cage size, where $r\u2192$ is the particle’s position vector. Because the neighbor cage restricts further motion along $\Delta r\u219201$, the displacement over the next time interval $\Delta r\u219212=r\u2192(2\Delta t)\u2212r\u2192(\Delta t)$ should, on average, be anticorrelated. The displacement $\Delta r\u219212$ can be projected onto the unit vector of the preceding displacement $\Delta r\u219201$, yielding the caging displacement projection (CDP) $x12\u2261\Delta r\u219212\u22c5\Delta r\u219201|\Delta r\u219201|$, which is parallel to $\Delta r\u219201$. For a caged particle, the ensemble-averaged CDP ⟨*x*_{12}⟩ is negative and its magnitude grows with $|\Delta r\u219201|$. In hard-sphere supercooled liquids, $\u27e8x12\u27e9=\u2212c|\Delta r\u219201|$ 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, ⟨*x*_{12}⟩ is independent of $|\Delta r\u219201|$.^{29}

In analyzing tracer dynamics, we first compute the non-Gaussian parameter for particle displacements *α*_{2}(Δ*t*) (see the supplementary material). The maximum in *α*_{2}(Δ*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 $\Delta r\u219201$ and $\Delta r\u219212$.^{29} The magnitudes of tracer displacements at *τ*^{*} vary across different matrices. To account for this variation, we report normalized quantities $x\u030312\u2261C\u22121\u27e8x12\u27e9$ and $r\u030301\u2261C\u22121|\Delta r\u219201|$, where $C$ is the square-root of the tracer mean-square displacement (MSD) Δ*r*^{2} at lag time *τ*^{*}.

### C. Isoconfigurational ensemble

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 *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)=\u27e8(r\u2192i(t)\u2212r\u2192i(0))2\Delta ri2\u27e9iso$, where $\Delta ri2$ is the ensemble-averaged mean-square displacement (MSD) of the *i*th 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 $\u27e8\Delta ri2\u27e9iso$, DP_{i} will be unity.^{34,35} Thus, examination of this quantity for all particles provides insight into the spatial distribution of dynamic heterogeneity.

### D. Trajectory shape analysis

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 $\u27e8\Delta r2(\tau cg)\u27e9=\sigma 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\u2211i=1ncg(x\u2192i\u2212x\u2192avg)2$, where *n*_{cg} is the number of coarse-grained points in the trajectory, $x\u2192i$ is the position of the *i*th coarse-grained point, and $x\u2192avg=1ncg\u2211i=1ncgx\u2192i$ is the mean position.

## III. RESULTS AND DISCUSSION

### A. Matrix dynamics

We first characterize the dynamics of the different glassy matrices through the ensemble-averaged mean-square displacement (MSD) Δ*r*^{2} (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,\tau )=\u27e81Ni\u2211k=1Ni\u2061exp[\u2212jq\u2192\u22c5(r\u2192k(\tau )\u2212r\u2192k(0))]\u27e9$, where $q=|q\u2192|$ is the wavevector magnitude, $j=\u22121$, *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^{0}–10^{3}, indicating interparticle caging [Fig. 3(a)]. By contrast, Δ*r*^{2} for L_{0.35} exhibits a small plateau at *τ* ≈ 10^{−1} followed by an increasing, subdiffusive power-law Δ*r*^{2} ≈ *τ*^{β} that extends to *τ* ≈ 10^{3}, 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 Δ*r*^{2} indicate that particles in L_{0.35} are more localized than those in L_{1.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 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*^{2} 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*^{2} 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*^{2} on intermediate time scales and prevents the matrix from fully relaxing on long time scales. Δ*r*^{2} 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^{−1}, qualitatively similar to but quantitatively lower than the one observed for L_{0.35}, which arises from attractive bond formation. On time scales *τ* ≳ 10^{2}, 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*^{2} 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.

### B. Tracer dynamics

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^{−1} ≲ *τ* ≲ 10^{3}. Diffusive dynamics are recovered on long time scales *τ* ≳ 10^{3}, 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}.^{18} The logarithmic slope *β* of the tracer MSD in the subdiffusive regime (Δ*r*^{2} ∼ *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^{−1} ≲ *τ* ≲ 10^{3} and diffusive for long time scales *τ* ≳ 10^{3} [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*^{2} in L_{0.35} and G_{0.20} are subdiffusive and nearly identical for time scales 10^{−1} ≤ *τ* ≤ 10^{1} but diverge for *τ* > 10^{1}. In the subdiffusive regime (10^{−1} ≤ *τ* ≤ 10^{1}), 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^{3}. By contrast, Δ*r*^{2} 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^{3}, 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^{5}. The matrix Δ*r*^{2} 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*^{2} 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*^{2} in G_{0.35} and G_{1.05} are nearly identical up to *τ* ≲ 10^{2}. 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^{2}, 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.

### C. Effects of matrix caging

Tracers in all matrices exhibit subdiffusive behavior on intermediate time scales 10^{−1} ≲ *τ* ≲ 10^{3}, 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\u030312$ as a function of the initial displacement magnitude $r\u030301$ for lag times *τ*^{*}. For a particle trapped in a harmonic well, $x\u030312$ varies linearly with $r\u030301$ with a slope of *c* = 0.5, indicating that the second displacement $\Delta r\u219212$ is anticorrelated with the initial displacement $\Delta r\u219201$.^{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\u030301$.^{27} For matrix particles in glassy liquids, by contrast, $x\u030312$ varies linearly with $r\u030301$ for small displacements but deviates from this initial linear behavior for $r\u030301$ beyond a length scale $r\u0303cage$ 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\u0303cage$ (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\u030301<r\u0303cage$. 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\u0303cage$.

Tracers with displacements larger than $r\u0303cage$ escape their cages and rearrange into new positions. Steeper slopes for the tracer CDP $x\u030312$ in the rearrangement regime indicate enhanced backdragging. For tracers within the repulsive liquid L_{1.05} and glass HSG_{1.00}, $x\u030312$ is approximately constant for $r\u030301>r\u0303cage$, indicating that the extent of backdragging does not increase beyond $r\u0303cage$ [Fig. 5(a)]. For the other matrices, by contrast, $x\u030312$ continues to decrease with increasing $r\u030301$ for $r\u030301>r\u0303cage$ 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\u030312$ for L_{0.35}, G_{0.20}, G_{0.35}, and G_{1.05} indicates heterogeneity in tracer dynamics for length scales beyond $r\u0303cage$ 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\u030312$ 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 *χ*_{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 $\chi 4(q,\tau )\u2261Nt(Fs2(q,\tau )\u2212[Fs(q,\tau )]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 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\u0303cage$ (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.

### D. Structural determinism of tracer dynamics

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 structure^{32} 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 $\Delta c,norm=\sigma DP2+\Delta c,normiso$.^{50–52} The first term $\sigma DP2=\u27e8(DPi)2\u27e9\u22121$ 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 $\Delta c,normiso=\u27e8DP2,i\u2212(DPi)2\u27e9$ is the dynamical variance, where $DP2,i=\u27e8(r\u2192i(t)\u2212r\u2192i(0))4\u27e9iso\u27e8\Delta ri2\u27e9iso2$. This variance quantifies the spread in DP_{i} between different isoconfigurational trajectories.

In L_{0.35} and L_{1.05}, the tracer structural variance $\sigma DP2$ increases up to and peaks at *τ* ≈ 10^{1}, corresponding to the time at which tracers experience maximum structural determinism [Fig. 7(a)]. On longer time scales, $\sigma DP2$ decreases and reaches a value near zero at *τ* ≈ 10^{3}, approximately the lag time at which tracer dynamics become diffusive in the MSD Δ*r*^{2} (Fig. 4). The structural variance $\sigma 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 $\sigma DP2$ arise from particle-to-particle variations in dynamic propensity. Hence, the larger values of $\sigma 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 $\sigma 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, $\sigma DP2$ is greater in attractive than repulsive arrested matrices $\sigma DP,HSG1.002<\sigma DP,G0.202$, and the disparity in $\sigma 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}), $\sigma DP2$ increases steeply for tracers for *τ* ≲ 10^{3}, appearing to plateau or even slightly decrease on longer time scales [Fig. 7(b)]. On time scales exceeding 10^{3}, $\sigma DP2$ increases very gradually or decreases slightly for tracers within G_{0.35} and G_{1.05}, respectively. The strong initial increase in $\sigma 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 $\sigma DP2$ to total tracer fluctuations within each matrix via the structural ratio $Rc=\sigma DP2/\Delta 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^{1} 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^{2}, *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 *χ*_{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 *R*_{c} for G_{0.20}, G_{0.35}, and G_{1.05} continue to increase for *τ* ≳ 10^{2}. 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 $\Delta 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.

### E. Consequences for tracer exploration

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.

## IV. CONCLUSIONS

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.^{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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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).

## REFERENCES

*ı*,