Small particle dynamics in glassy polymers: diffusion, relaxation, and machine-learned softness

Abstract

In this work, we explore the simulated transport dynamics of small particles through a polymer melt at temperatures spanning from liquid-like behavior down to near and below the simulated glass transition temperature. Using softness, a machine-learned scalar quantity, we connect the structural neighborhood surrounding a particle with its dynamic behavior to relate the probability of a glassy rearrangement to its local environment. An energetic and entropic scale for the rearrangement process of these small particles emerges and is compared across systems of different particle sizes and interaction strengths. Diffusion coefficients and relaxation times both show strong dependence on our tuning parameters, and the barriers to rearrangement show increasing nonlinearity as the particles become smaller. The trends we observe provide some insight into how local structure plays a role in small molecule transport when the surrounding medium is undergoing a glass transition, leading to large changes in system mobility.

---

Introduction

Chemical separations processes like distillation and other thermally-driven methods account for nearly fifteen percent of global annual energy consumption.¹ The gradual reduction of the world's energy budget and the carbon intensity of the energy economy have been identified as key elements of a strategy to combat climate change.^2,3 Increasing concentrations of atmospheric carbon dioxide and other pollutants⁴ have been shown unequivocally to be direct causes of climate change⁵ and could be mitigated through either point-source capture upon emission,^6–8 or direct air capture of legacy emissions from the atmosphere.^9,10 The need for not only a reduction in global energy consumption but also with more advanced and increasingly dilute separations processes points toward polymer membrane, pressure-based separations processes as an attractive alternative. Polymer membranes for gas separations have been in use for decades^11,12 and continue to be an area of intense research to push the limits of performance^13,14 and drive the development of new hydrocarbon materials that exhibit the potential for fast and energy-efficient separations.^15–17

Many of these polymer membranes fall under the umbrella of glassy disordered materials,¹⁸ which are amorphous and lack long-range structural order like liquids yet possess many of the physical properties of solids.^19–21 Dynamical heterogeneities in microscopic behavior of these materials below the glass transition temperature T_g are observed, along with a growing structural relaxation time to around 10² seconds alongside the dynamic viscosity increasing to greater than 10¹³ centipoise.²¹ While ubiquitous in glass-forming materials, these behaviors are not entirely explained by any one theoretical approach or understanding.^22–28

Building a clearer understanding of structure–function relationships for transport in polymers has important implications in materials design and application for gas separations.²⁹ Glassy polymer membranes used for gas separations processes are characterized by high solubility coefficients for gas species, moderate diffusion selectivity, and nonequilibrium free volume distributions that allow glassy polymers to dominate the limits of performance compared to rubbery polymers.¹⁴ The dynamics of small tracer molecules within polymers can be studied experimentally using photobleaching experiments^30,31 and tools like quasi-elastic and inelastic neutron scattering.³² Simulations present an important complement to experimental techniques for characterizing materials that have not yet been synthesized or extensively studied at the lab scale.

Computational studies of gas diffusion through polymers have been performed extensively³³ over the years in both coarse-grained^34–37 and atomistic resolutions,^38–41 inspired by classic theoretical frameworks used to describe diffusion in disordered media.^42,43 Early atomistic simulation of oxygen transport through linear n-alkane systems demonstrated that diffusion in glassy systems is governed by a hopping mechanism where particles spend much of their time in a single region before moving to another potential energy basin.⁴⁴ The hopping rearrangement process in glassy systems related to cooperative motion with the nearby polymer chains and does not necessarily correlate with any gradient in local free volume.^45,46 The transitions between these localized regions of permeation have been further studied with other traditional methods such as kinetic Monte Carlo in both glassy and rubbery polymers,^35,47,48 which have also found utility in other nanoporous materials such as zeolites.⁴⁹

In the last decade or so, the advent of machine learning and increased capacity for computer simulation of large size and longer timescales has unlocked a wealth of new tools and measures to characterize the dynamics and plasticity of glassy materials.⁵⁰ One such measure is a scalar quantity called softness, which relates a glassy material's dynamic behavior to its static configuration at a given point in time by using a support vector machine (SVM) to analyze a high-dimensional space of local structural features.^51–53 Softness has proven to be useful in gaining understanding around a wide variety of phenomena in glassy materials including the impacts of structural aging,⁵⁴ spatial variations in thin-film relaxation,⁵⁵ the behavior of polymer nanocomposites under deformation,^53,56–59 and the relationship between static excess entropy and dynamics in the supercooled regime.⁶⁰ When compared with other methods that relate local structure to free volume, such as Voronoi-based indicators like local density, free volume, or Voronoi anisotropy and its divergence, softness has been shown to outperform them in simulations of both glassy quenched and high-temperature Lennard-Jones systems; indeed, using other metrics to measure local free volume and relate it to dynamics has historically been a challenge.^61–63 Softness has been successfully used both in experimental systems⁵¹ and in simulations of simple models like the Kob–Andersen glass former^52,54,60,64 or coarse-grained bead-spring polymers,^{55–59,65,66} making it a versatile approach for characterizing glassy dynamics.

The importance of local structure to diffusion and hopping dynamics makes it tempting to use softness or other high-dimensional machine learning tools to investigate the dynamics of small penetrants through a glassy matrix. Softness has been applied extensively to the binary Kob–Andersen system and to bead-spring polymers, but to our knowledge softness has not yet been used to explore the relationship between local structure and dynamics of small particles moving though a glassy medium. Recent work characterizing the diffusion behavior of small penetrants in both rubbery and glassy polymer materials has included the application of self-consistent nonlinear Langevin equation (SCNLE) theory⁶⁷ and self-consistent cooperative hopping theory (SCCHT),⁶⁸ which may be used to relate the penetrant–segment size ratio, strength of penetrant attraction to the polymer, and temperature to the diffusion.^69–71 These theories are based on the idea that the dynamics in the glassy state are slowed due to a combination of a local barrier for a segment (polymer monomer or a penetrant) to escape its local cage, and the barrier is augmented at low temperatures by the elastic distortion of the matrix required for the cage to dilate and allow the segment to escape. These theories have been corroborated by not only existing experimental measurements of small molecule diffusion, but also molecular dynamics simulations of the purely repulsive penetrant dynamics of small, gas-like molecules in glassy polymers.^34,36,37 Building on this work by examining diffusion behavior of particles that have attractive interactions with the surrounding polymer matrix may give some insights into how systems with real chemical interactions can be tuned to promote desirable transport properties.

In this work, we take the framework of softness and use it to investigate the dynamics of small particle transport through a model glass-forming polymer in the spirit of classic the Kremer–Grest model.^65,66 We use molecular dynamics (MD) to simulate gas-like, monatomic tracer particles permeating a polymer melt as it nears the glass transition, measuring the self-diffusion and relaxation behavior of the tracers, and compare the tends in mobility from these measures qualitatively with the SCCHT model for activated diffusion in a polymer matrix with cross-interactions. We then calculate softness for these systems and extract the characteristic energetic and entropic scales that describe the structural free energy barriers to tracers’ glassy rearrangement. Finally, we compare these barriers across systems with different tracer particle sizes and interaction strengths with the polymer medium to assess how softness can relate local structure to transport dynamics in glassy systems. We examine the differences between energy barriers computed from softness and from diffusion coefficients or relaxation times and theoretical predictions for activation barriers in penetrating improving the overall understanding of connections between these barriers to activated transport is important in its own right and may be practically helpful in designing new materials that exhibit desirable diffusion behavior for use in separations processes and other important engineering applications.

Methods

System initialization and equilibration

The simulations in this work consider bead-spring polymer chains^65,66 consisting of 20 Lennard-Jones (LJ) interaction sites, each representing a monomer of species A. Gas-like tracer particles, which we refer to as species B, are represented as freely moving single LJ interaction sites. We use the standard Lennard-Jones reduced units throughout for time , temperature (T* = T/(ε/k_B)), pressure (P* = P/(σ³/ε)), volume (V* = V/σ³), and energy (U* = U/ε), where σ, ε, and m are the length, energy, and mass units of the reduced Lennard-Jones potential and k_B is Boltzmann's constant.

The reduced Lennard-Jones potential^72,73 in its cut-and-shifted form,

is used for all non-bonded interactions in the system with a cutoff distance of r^nb_c = 2.5σ and shifted by U_cut to be zero at r = r_cut. The polymer–polymer site interaction size and energy are fixed at σ_AA = 1.0σ and ε_AA = 1.0ε respectively. The tracer–tracer interaction sizes are varied at values of σ_BB = 0.3, 0.5, 0.8 and 1.0σ while the tracer–tracer interaction energies are fixed at ε_BB = 1.0ε. The polymer–tracer interaction sizes σ_AB are set by the Lorentz rule,⁷⁴ , while the polymer–tracer interaction energies are chosen to be either ε_AB = 1.0, 1.2, or 1.5ε. A combining rule is not used for interaction energies to suppress spontaneous phase separation at lower temperatures. A unit Lennard-Jones mass of m = 1 is assigned to both the individual monomer beads and to the tracers.

Adjacent monomer beads along a single polymer chain are bonded together and do not experience the U^nb potential with one another. The bonded interactions in our system are defined by the finitely-extensible nonlinear elastic (FENE) bonding potential,^65,66 which has the form

with bond stiffness K = 30ε/σ² and maximum extensibility of R₀ = 1.5σ. The first term describes the restorative force of the bond; the second term is the complete non-bonded LJ potential between two monomers; the third term shifts the magnitude of the LJ potential to be zero at its minimum. The second and third terms collectively form the purely repulsive Weeks–Chandler–Anderson (WCA) potential,⁷⁵ which is truncated at its minimum value; for this system, this corresponds to r^nb_c = 2^1/6σ.

All simulations are performed in the LAMMPS^76,77 package with data files initialized by a custom Python script. Each simulation box in our work contains 500 tracer particles and 475 polymer chains, or 10 000 total particles. Since all particles have a unit mass, the mass fraction of tracers in each simulation box is equal to the number fraction of 5%. The total volume fractions of tracers, taken to be the Lennard-Jones volume of the tracers divided by the total Lennard-Jones volume of all beads in the box, correspond to ϕ_B = 0.142%, 0.654%, 2.62%, or 5% for tracers of size σ_BB = 0.3σ, 0.5σ, 0.8σ, and 1.0σ, respectively. For the larger tracer sizes, these are outside of the dilute limit.

To ensure that our starting configuration is not too dense, polymer chains are initialized using a random walk in a large cubic simulation box, after which the tracer particles are inserted randomly. The soft push-off method is then used to relax the overlapping contacts, followed by a final energy minimization.^78,79 After the minimization procedure we perform MD using an isothermal–isobaric (NPT) time integrator^80,81 with a timestep of dt = 0.002τ and damping coefficients of T_damp = 100dt and P_damp = 1000dt for the temperature and pressure, respectively. Velocities are randomly initialized at T* = 1.0 and the system is run at a fixed temperature and pressure of T* = 1.0 and P* = 0.0 for 5 × 10⁶dt or 10⁴τ to equilibrate, at which point the mean-squared displacement of the monomers is several times the mean end-to-end distance of the polymer chains. The equilibrated snapshot is then saved for further use.

Cooling and production runs

An initialized configuration at T* = 1.0 and P* = 0.0 is decreased gradually to zero temperature at a dimensionless cooling rate of Γ = 10⁻⁴T*/τ and constant P* = 0.0. Snapshots are saved along the way at incremental temperatures. The density of the system is recorded as a function of time to identify the temperature at which a glass transition occurs, which we estimate using a nonlinear fitting function used previously,⁸²Here, V* is the system volume (inverse density), w is the width of the transition window in temperature units, M and G are the thermal expansion coefficient of the melt and glass, respectively, and c is the value of log V* at the glass transition temperature .

Simulation snapshots from temperatures of interest are then used for production runs at constant temperature and pressure, using the T* corresponding to their ramping time and a constant P* = 0.0 for a total of 5 × 10⁷dt or 10⁵τ. During these quiescent simulations, several trajectories of 101 frames each are recorded with linear spacing Δt between frames corresponding to each decade from 10⁻²τ to 10³τ to calculate relevant dynamic properties, as well as the radial distribution function g(r). For example, a trajectory is recorded during the entire duration of the simulation with frame spacing of Δt = 10³τ. Following the production run, a trajectory spaced by Δt = 10⁻¹τ with 10⁴ frames is generated spanning 10³τ or 1% of the production time; this is used for evaluating softness in the postprocessing steps, discussed below.

We characterize the structure of our systems using the standard radial distribution function,⁸³ and the dynamic properties are characterized using two metrics. First, we compute the mean-squared displacements⁸⁴ as

where r_i(t) is the position of particle i at some time t. For a trajectory with many samples of displacements over a time lag Δt (i.e. from 0 to Δt, from Δt to 2Δt, etc.) the ensemble averaging uses the moving time origin method, where M is the total number of frames in the trajectory spaced evenly by time Δt, ℓ is the integer number of frames between which we are looking at positional changes, and r_i,m is the position of the particle i in frame m. This expression yields MSD values for all lag times in ℓΔt, with smaller times being averaged over many more independent displacements than larger times. In our analysis using trajectories of 101 frames spaced by a time interval which is a power of 10 ranging from −2 to 3, we keep only the data calculated for values of ℓ up to 10 for all but the largest interval of Δt = 10³τ. The second dynamic measure we use is the intermediate scattering function,where k is the wavevector associated with the calculation, M is the total number of frames in the trajectory spaced evenly by time Δt, ℓ is the integer number of frames between which we are looking at positional changes, and r_i,m is the position of the particle i in frame m. F_s(k,ℓΔt) measures the structural relaxation time on length scales of the order ∼2π/k.

Dynamical indicator function and softness

Following the production run and the generation of a 1000τ trajectory with frames spaced by Δt = 0.1τ for each system, we evaluate the softness of each particle in each frame to quantify the structure of the neighborhood around a particle with a single parameter. To describe the local environment around a particle, we employ a set of two-body Gaussian radial structure functions G inspired by Behler and Parrinello⁸⁵ that have been extensively used for softness,^51–60where we are probing the density of the j particles of species X at some radius r away from our central particle i. The magnitude of the pairwise separation of the two particles is |r_ij| and the Gaussian smearing has some standard deviation which we fix at δ = 0.1σ. We use a set of 60 structure functions, 30 each probing the density of species A and B in the radial direction out to a distance of 3.0σ from the central particle; this encapsulates the force cutoff of 2.5σ and multiple peaks of g(r) in the system.

Using the set of structure functions, a 60-dimensional feature vector F_n is constructed for the particles, the elements of which are used as training data. The training data are divided into an equally sized set of rearranging and non-rearranging particles, which correspond to some binary variable y with values 1 and 0 respectively. We choose an indicator function to assess whether particles in our training set are currently undergoing a rearrangement. The dynamic quantity p_hop has been extensively studied not only as a way of monitoring glassy rearrangements during simulation^86,87 but also as a useful indicator function for evaluating softness.^{52,54,55,57,60} This function has the form,

where i is the particle index, r_i(t) is its position at time t, and the angle bracketed quantities 〈〉_A and 〈〉_B are averages over the time intervals A and B. These intervals are determined by a rearrangement timescale t_R,

A = [t − t_R/2, t] (8)

B = [t, t_R/2] (9)

and correspond to the times before or after rearrangement occurs. The value of p_hop for a particle tends to be very small during times between rearrangements (caging behavior) and peak during a rearrangement from one region to another.^86,87 The value of t_R is chosen to be smaller than the average relaxation time of the tracer particles system. Particles which are rearranging are chosen to have a p_hop value greater than p_cut,R = 0.25 while non-rearranging particles are those whose value of p_hop remains below p_cut,NR = 0.02 for longer than at least the rearrangement timescale t_R. These choices are motivated further in Appendix S1.

Using the rearranging and non-rearranging training data, a linear support vector machine (SVM) is used to construct a hyperplane that best separates the two classes in the high-dimensional space defined by the F_n feature vectors. The hyperplane, defined by,

w_n·F_n − b = 0 (10)

gives us a vector of weights w_n which correspond to each structure function and b is an offset value. We then define the softness, S,

S = w_n·F_n − b (11)

which is a scalar value that contains the data captured by the feature vector F_n for each particle in each frame. To quantify the accuracy of the hyperplane, we usewhere μ_S is the mean of the distribution of all softness values and P(S|R)dS is the normalized probability distribution softness for all rearranging particles. This is the fraction of all rearranging particles’ softness values that lie above the mean of the distribution for the entire trajectory from which the training data is drawn.

Softness is known to provide a relationship between its conveyance of local structure and the probability of rearrangement through a temperature-independent free energy of rearrangement. The rearrangement barrier ΔF(S) is taken from the Arrhenius definition of the rearrangement probability,

ΔF(S) = ΔE(S) − TΣ(S) (14)

where ΔE(S) and Σ(S) are structure-dependent energetic and entropic barriers to the rearrangement process. Prior work on softness^51,52,58 has shown these two quantities to be approximately linear in softness, though we make no prior assumptions about their form for the gas-like tracers in our system.

Results

Glass transition temperature

The glass transition temperature T_g of our neat polymer is computed from the thermal expansion coefficient definition and plotted below in Fig. 1 on the left-hand side. The glass transition temperature is identified to be about in the neat system. The uncertainty and error bars here and elsewhere in this work is defined as the standard deviation of a measured value from three independent configurations. The right-hand side of Fig. 1 shows the calculated values of for different combinations of σ_BB and ε_AB, and compares them with that of the neat system. For all the systems considered, falls within about 4% of the value of the neat system, with the lowest values for a given ε_AB occurring at the σ_BB = 0.5σ size. For tracer particles of size σ_BB = 0.5σ and above, the value of increases with the strength of the cross-interaction ε_AB. The changes in due to the presence of tracers indicate that their concentration in our systems is beyond the dilute limit.

Fig. 1 (a) Plot of log V* vs. T* for neat polymer system with indicated. (b) Gass transition temperatures for various combinations of σ_BB and ε_AB, compared with that of the neat system indicated by the purple horizontal line.

Temperature-dependent diffusion

The ability of tracer particles to diffuse through the polymer system is strongly influenced by temperature, as diffusion of small molecules in porous media is typically an Arrhenius process. Fig. 2(a) shows the MSD of tracers with varying size at two temperatures. At a high temperature of T* = 0.95, over twice , the tracers all transition rapidly from the MSD ∝ τ² ballistic regime to MSD ∝ τ¹ diffusion. For a lower temperature of , there is a reduction in the mobility related to the polymer glass transition. While the tracers of the smallest size are readily diffusive, the larger particles see a precipitous decrease in their mobility, with the σ_BB = 1.0σ tracers moving on average about 10 particle diameters over 10⁵ time units. Tracers of the same particle size as a single monomer bead have only a few times greater mobility than the smallest tracers at temperatures near the glass transition temperature. The rigid, glassy structure of the polymer at low temperatures prevents the larger particles from being able to move through the interstices of the glassy matrix without overcoming exceedingly high energy barriers to rearrangement. Changes in the interaction potential ε_AB between the tracers and the surrounding polymer at a constant size σ_BB also impact the MSD, as shown in Fig. 2(b). The increase in the depth of the potential well impedes the motion of the tracers, though the impact is minimal for the particles above T_g.

Fig. 2 MSD of tracer particles with (a) varying size σ_BB and constant ε_AB = 1.0ε at high and low T*, and (b) varying interaction strength ε_AB and constant size σ_BB = 1.0σ at high, moderate, and low T*. Error bars are smaller than the symbols’ size, and connecting lines are drawn to guide the eye.

The self-diffusion coefficient D is measured using the long-time slope of the MSD. The Einstein relation for diffusion⁸⁸ in three dimensions gives the dependence of D,

where the limit t → ∞ is where MSD ∝ t¹. Values of D are found by taking a linear fit with a slope of 1 to the data, log MSD = b + log τ, and transforming its y-intercept b back to a diffusion coefficient, D = e^b/6, to get the result in units of σ²/τ. A higher relative value of D indicates that under those conditions one type of particle has on average greater mobility. The values of D are shown in Fig. 3.

Fig. 3 Diffusion coefficients D for tracer particles as a function of for varying σ_BB and ε_AB symbols indicate σ_BB, while subplots indicate ε_AB values of (a) 1.0ε, (b) 1.2ε, and (c) 1.5ε. VFT fit is shown with solid connecting lines and extrapolated with dotted lines. Arrhenius fits to high-temperature region are shown with dashed lines. Error bars are smaller than the symbols’ size.

For the smallest of the particles, σ_BB = 0.3σ, we see that the relationship between D and temperature is Arrhenius, until just above , when the value of log D no longer decays linearly with inverse temperature and the activation energy increases. The increase in activation energy is further exacerbated by increases in ε_AB. The larger sizes of particles, σ_BB ≥ 0.5σ, which are closer in size to a segment of the polymer chain exhibit a more pronounced drop in diffusion that becomes stronger with increasing σ_BB and ε_AB. This makes sense intuitively and is analogous to the predicted increase in the relaxation time τ_α for penetrant molecules with cross-interactions in a polymer matrix by SCCHT, which demonstrated that either increasing tracer size (σ_BB) or attraction strength (ε_AB) would lead to an increase in both τ_α and the activation energy. The diffusion coefficient of small penetrants is generally known to be size-dependent to some extent, and we examine this relationship below.

For each system, the results for D may be fit to the empirical expression of the Vogel–Fulcher–Tammann (VFT) type,

where T_VF is a divergence temperature which in all these cases lies below T_g and B is a material property of the system.^20,89 This is consistent with our identification of a glass transition in our polymer system and the tracer particles experiencing a dynamic slowdown as they approach those conditions.

Experiment and theory has demonstrated that, in polymeric supercooled liquids, the diffusion coefficient obeys a dependence such that log(1/D) ∼ (T_g/T)^λ, where the exponent λ is related to system chemistry and is roughly 2 for nonpolar liquids and in the range of 2.2–2.6 for some typical polymers.⁷⁰ SCCHT predictions for this dependence note that in the low-size-ratio limit, the λ dependence should hold, while in the high-size-ratio limit the exponential scaling should be 3λ, which is due to differences in the local cage vs. elastic barrier being the limiting factor. We can examine our own data in this lens by performing a linear fit to extract the scaling exponent λ and the prefactor K in the expression , where we de-dimensionalize the diffusion coefficient by the prefactor from (16). To do this we simply plot , and find K = e^b and λ = a. Linearized plots listing the fit parameters when the data above the glass transition region (0.45 ≤ T* ≤ 0.75) can be found in Fig. S2.1. The fit parameters are shown in Fig. 4.

Fig. 4 Diffusion coefficient fit parameters as they vary with σ_BB and ε_AB. Scaling coefficient λ is shown in (a) with lines to guide the eye, while the prefactor K is shown in (b) with the logarithmic fit to the data indicated.

We identify that, for our smallest particles of size σ_BB = 0.3σ, the scaling exponent is around λ = 1.2, which increases with tracer size to a value around λ = 2 or just greater. The prefactor K appears to be a logarithmic function of the tracer size, as has been predicted for the large-size-ratio limit in the SCCHT. While generally the value of the scaling exponent increases with tracer size, there is no obvious functional form. It is important to note that in the SCCHT model, the value of the scaling exponent is not a fit parameter – it is taken from experiments or theoretical predictions. While we can draw qualitative and some quantitative comparisons to what one might expect from theory, it is important to note that we are limited by the dynamics in the supercooled and glassy regime that can be reasonably observed in simulations (this is further discussed at the end of the next section). It is important to consider how we might quantify the barrier to the diffusion process based on our calculated values of D, so we might compare them with other energy scales. Recall that for an Arrhenius process, the activation energy E_a is determined from the slope of a line of the logarithm of the measured property versus inverse temperature,

which is constant across temperatures. We have identified that diffusion is a non-Arrhenius process as we approach the glass transition, so the activation barrier from the instantaneous slope on an Arrhenius-like plot will be dependent on temperature. We can approximate this temperature-dependent energy barrier by substituting (16) into (17) to get a general expression,which gives an energy barrier associated with a process that has a VFT-like temperature dependence.

Tracer relaxation behavior

Tracer particle relaxation times in the polymer melt are calculated by evaluating the SISF for the wavevector of magnitude k = 7.14σ⁻¹. This value of k approximately corresponds to the first peak of the static structure factor S(k) for the system and serves as a characteristic inverse length scale of the polymer that is constant across all combinations of tracer σ_BB and ε_AB values. As previously with the diffusion coefficients, Fig. 5(a) shows the effect of varying particle size at a high and low temperature while Fig. 5(b) shows the impact of changing the interaction strength for a constant particle size.

Fig. 5 SISF of tracer particles with (a) varying size σ_BB and constant ε_AB = 1.0ε at high and low T*, and (b) varying interaction strength ε_AB and constant size σ_BB = 1.0σ at high, moderate, and low T*. Connecting lines are drawn to guide the eye.

As cooled towards , we see the emergence of a two-step relaxation behavior in the decay of the SISF, one of the hallmarks of glassy dynamics.⁸⁹ The shorter β-relaxation takes place on a time scale ∼1τ, while the longer α-relaxation decays on longer time scales as either σ_BB and ε_AB increase or with decreasing temperature. In the least mobile of particles, a caging plateau like that observed in the MSD of particles near emerges. The relative kinetic arrest of the polymer at means that the motion of polymer chains restricts particles whose sizes are comparable to one of the segments along the chain. To estimate the relaxation time τ_α from the decay of F_s, the data are fit to a sum of two Kohlrausch–Wiliams–Watts (KWW) stretched exponential functions,^90,91

where the first term represents the fast (beta, β) relaxation with its own timescale and stretching exponent, A captures the fraction of the decay due to the slower α process, and b_β and b_α are stretching/compression exponents. When a system exhibits dynamic heterogeneity, typically β_α < 1. We use the fitted function to calculate an effective relaxation time τ^eff_α, by integrating the KWW function over the time domain,⁹²where Γ(x) is the gamma function.

The calculated values of τ^eff_α are plotted as functions of inverse temperature in an Arrhenius-style manner for each system in Fig. 6, and the stretching exponents b_α in Fig. S2.2. As we can see, relaxation times are super-Arrhenius in nature, increasing by several orders of magnitude in the larger particles as the glass transition is approached. Much like the behavior of the diffusion coefficients, this is to be expected for a system with glassy dynamics and the measured values can be fit to a VFT expression. As with the empirical fits to the diffusion coefficients, a temperature-dependent activation barrier may be found by using (16) and (18) with the replacement B → −B to account for the inverse temperature dependence of τ_α compared to the diffusivity.

Fig. 6 Effective relaxation times τ^eff_α for tracer molecules as a function of . Symbols indicate σ_BB, while subplots indicate ε_AB values of (a) 1.0ε, (b) 1.2ε, and (c) 1.5ε. VFT fit is shown with solid connecting lines and extrapolated with dotted lines. Arrhenius fits to high-temperature region are shown with dashed lines. Error bars are smaller than the symbols’ size.

The stretching exponents shown in Fig. S2.2 are also below a value of 1.0 at temperatures below the onset of glassy dynamics, indicating clearly that there is dynamical heterogeneity in the motion of the particles. Calculated values of β_α larger than unity at high T* are due to the interference of ballistic dynamics with the α relaxation at temperatures where there is not a pronounced caging plateau.

By inspection, the relaxation times increase with particle size and with interaction strength, meaning that smaller particles will be able to overcome the dynamic slowdown at low temperatures, but a decrease in size may not be enough to overcome particularly strong attraction between the chains and the species of interest. For example, in our results, we find some overlap in the relaxation times for the high size ratio particles with σ_BB = 0.8 and 1.0σ when the value of ε_AB is increased. Theoretical predictions for activated penetrant dynamics have previously identified such non-monotonic behavior in the relaxation time when size ratio and interaction strength are increased.⁷¹ We fit the ratio of the relative effective relaxation time τ^eff_α/τ₀, where τ₀ is the prefactor in a fit to eqn (16), to an inverse power-law relationship with a scaling exponent λ. Like the fit for the inverse diffusion coefficients in the prior section, we use data above the glass transition region (0.45 ≤ T* ≤ 0.75). Results for the fit parameters are shown on the linearized plots in Fig. S2.3. While the prefactor K was not found to be a logarithmic function of tracer size for the relaxation time data as it was for the diffusion data, the value of λ increases with both σ_BB and ε_AB as before.

The dynamic slowdown of the penetrants in the system can be quantified relative to the dynamics of the polymer itself. Following the same method outlined above, we calculate the values of the polymer segmental relaxation time τ^eff_α across different systems and temperatures and compare the polymer relaxation times in the presence of penetrants with results from simulation of a neat polymer. Fig. 7 shows how the polymer segmental relaxation time τ^eff_α,poly evolves with , parametrized by either σ_BB or ε_AB of the penetrant species. It is important to note the collapse in τ^eff_α,poly across different systems when temperature is scaled by , showing that the segmental dynamics are perturbed by the presence of the penetrating species insomuch as the penetrants affect the glass transition temperature.

Fig. 7 Polymer segmental relaxation times τ^eff_α plotted against . Panel (a) shows the effect of changing σ_BB for ε_AB = 1.0ε while (b) shows the effect of changing ε_AB for σ_BB = 1.0σ. Stars indicate the segmental relaxation time in a neat polymer system. Note the collapse at constant temperature relative to .

The dynamical decoupling between the mobile penetrants and the polymer matrix can be quantified by plotting the ratio of the relaxation times of the penetrants to the polymer monomers. A ratio of the penetrant-to-segmental relaxation timescales on the order of 1 would imply that the penetrant dynamics track the polymer segmental dynamics, while a smaller value indicates that the tracers are able to resist the slowdown. From the ratios plotted in Fig. 8, we find that tracer dynamics become increasingly decoupled from the matrix as is approached from above. Increasing the interaction strength ε_AB indeed increases the ratio, confirming our prior observation that an increase in the interaction strength leads to reduced mobility of penetrants. For the largest penetrants, increasing the interaction strength above that of the polymer–polymer interaction (ε_AB > 1ε) slows their relaxation times beyond those of the polymer segments. This is consistent with results from a theoretical approach in SCCHT tailored to polymer–penetrant interactions in a homopolymer liquid.⁷¹ That work found a decoupling of the relaxation times with increasing effective packing fraction ϕ_eff of polymer in the system, which approaches a kinetic glass transition in a manner equivalent to our temperature approach to , and that increasing the strength of the interaction at a given value of ϕ_eff (or equivalently T*) reduces the rate of decoupling. In the case of theoretical systems where penetrant–polymer interactions were not considered, packing fractions denser than a glass transition of vitrification point ϕ_glass (or equivalently, temperatures below ) are needed to observe this dynamic decoupling.

Fig. 8 Ratio of tracer-to-polymer relaxation times fitted from eqn (20) plotted against . Symbols indicate σ_BB, while subplots indicate ε_AB values of (a) 1.0ε, (b) 1.2ε, and (c) 1.5ε.

The SCCHT approach also predicts that there is a decrease in the decoupling at equivalent packing fractions with increasing size, meaning that the ratio of relaxation times is larger. Interestingly, our work finds that tracers of size σ_BB = 0.8 and 1.0σ show no dynamic decoupling through their approach to for the time scales accessible in our simulations. Though the particles with the same size as a Kuhn segment modeled using SCCHT had a slight increase in the calculated valued of τ^eff_α,tracer/τ^eff_α,polymer with increasing ϕ_eff, they saw a sharp decrease in the region where ϕ_eff reflects a glass transition or vitrification of the system.⁷¹ One possible explanation for this is that the systems modeled in this work feature a high enough concentration of penetrants that the dilute limit predictions from SCCHT begin to break down; we find that for all particle sizes, the glass transition temperature is perturbed from the neat polymer system, which necessitates being outside the dilute limit. We have, however, identified that the relaxation time of the polymer is roughly constant across systems at the same value of , meaning that relative mobility should be comparable across systems at the same relative temperature. It is also relevant to note that in the SCHHT theory, the interactions between penetrants and polymer beads have a range of only 0.15 times a segment diameter, meaning that the contacts must be very close for them to interact; in this work, the Lennard-Jones potential has an attractive radius 2.5σ, over 16 times that employed in SCCHT. That work found that increasing the standard range of attractions led to an increase in relaxation time, particularly at high penetrant–polymer interaction strengths.⁷¹ Finally, it is known that the accessibility of glassy temperature regimes in equilibrium simulations is subject to greater limitations than in experiments or theory because of the long timescales in real systems and the difficulty in gathering adequate statistics. Simulations on a longer timescale (comparable to that of experiments) that were able to reach an equilibrium closer to or just below the glass transition temperature may demonstrate the beginnings of dynamic decoupling for our largest penetrant particles through the polymer matrix. This would be similar to the predictions in the theoretical approach with the ϕ_eff ≈ 0.64 reported previously.⁷¹ These known limitations in molecular dynamics studies of glassy systems mean that an alternate approach would be warranted in future work.

To directly compare the dynamics in the systems we have studied and to draw analogy to existing experimental and theoretical work on small particle transport, we can also plot the decoupling factor for the tracer molecules, or the product of the diffusion coefficient D and the effective alpha-relaxation time τ^eff_α. For systems where there is dynamic decoupling between penetrant and matrix, as demonstrated in Fig. 8, this can be plotted as the product of the penetrant D and the polymer τ^eff_α. Plots of these quantities similar to Fig. 8 are included in Appendix 2 as Fig. S2.4 and S2.5. The decoupling factor remains less than 1 for all cases studied in the tracer–tracer definition, while for the tracer–polymer definition it demonstrates clear decoupling, Dτ ≫ 1, for large particle sizes as is also seem in Fig. 8.

Softness and barriers to rearrangement

Using the training parameters described in the prior section (p_cut,R = 0.25, p_{cut_NR} = 0.02, t_R = 10τ, and t_NR = 10τ), and a set of M = 60 structure functions surrounding each tracer particle, we trained a hyperplane to describe the softness of our tracers for each of the four values of σ_BB at a fixed interaction strength of ε_AB = 1.0ε. We trained our systems at different temperatures to ensure that enough non-rearranging (in the case of smaller particles) or rearranging (in the case of larger particles) could be identified for the set. The training temperatures and the calculated accuracies using eqn (12) are shown in Table 1. A sensitivity analysis of prediction accuracy and the size of the training set versus the choice of p_cut,R and p_cut,NR for each value of σ_BB with ε_AB = 1.0ε was performed and results are shown in Appendix 3. The non-rearrangement threshold of p_cut,NR = 0.02 generally corresponds to the value below which there will be less than several thousand examples of both rearranging and non-rearranging particles in all cases. The rearrangement threshold of p_cut,R = 0.25 ensures that we are probing displacements that are at least size 0.5σ. Choosing a value below this would increase the size of the training set but not necessarily the accuracy of the classification.

Table 1 Training temperatures and resulting accuracies for varying σ_BB systems at ε_AB = 1.0

System (σBB)
0.3σ	0.300	80.68%
0.5σ	0.375	77.25%
0.8σ	0.450	74.63%
1.0σ	0.450	66.56%

---

Overall, the accuracy of rearrangement prediction from softness with this set of training parameters and structure functions decreases with increasing particle size. Typical accuracies reported in other computational studies vary from 70% to 90%.^58,60 While tuning these parameters further may be able to improve the prediction accuracy for a specific set of interactions, since we are interested in trying to directly compare the barriers derived from softness across several systems, it is important that our definition of what a rearrangement is consistent across systems. A constant cutoff value p_cut,R ensures that we are probing rearrangements of at least 0.5σ in all cases. The influence of the choice of p_cut on our results is discussed further at the end of this section.

Following the training procedure, we take the hyperplane constructed for some σ_BB and ε_AB = 1.0 and evaluate S for all our particles, using this particle-size hyperplane for each interaction strength to directly compare the values for similar systems. For each system replica and at each T*, we build two histograms of softness values: one for all particles, and one for rearranging particles, both with a bin size of 0.1 softness units. For a given pair of histograms, the value of the rearrangement probability P_R(S) is calculated as the ratio of the two bins with the same softness. As an example, Fig. 9(a) shows the overlaid softness histograms for three independent replicas of the σ_BB = 0.5σ and ε_AB = 1.0ε system at its training temperature of T* = 0.375. The dynamical decomposition of the same system from which the softness dependent energy and entropy barriers are extracted is illustrated in Fig. 9(b). Equivalent plots to Fig. 9(b) for other combinations of σ_BB and ε_AB may be found in Appendix 4.

Fig. 9 Softness and rearrangement data for the σ_BB = 0.5σ and ε_AB = 1.0ε system. (a) Normalized softness distribution histograms for the training data at the training temperature of T* = 0.375. (b) Dynamical decomposition of softness for the system used to derive the energetic and entropic barriers to rearrangement.

From the histograms, we obtain a normalized probability distribution function (PDF) with some mean μ and standard deviation S. The values of μ_S and μ_R, with subscripts S and R referring to the distribution of all softness values and those of rearranging particles respectively, are plotted in Fig. 10. It is interesting to note that, for the plots of μ in the cases of σ_BB = 0.8σ and σ_BB = 1.0σ, the slope of S(T*) of all particles has a change in value near . In contrast, the smaller particles for which mobility is less impaired below , the temperature dependence remains linear with no change in slope. This gives insight into how smaller particles experience the effects of the glass transition occurring in the medium around them, or rather, how they do not.

Fig. 10 Mean values μ of average softness distributions for (a)–(d) all particles and (e)–(h) rearranging particles as functions of temperature. Error bars indicate the standard deviation s of each average distribution.

It is important to understand the difference between the softness distributions of the rearranging and of all particles. When local structure plays a greater role in determining the rearranging behavior, the distributions will become less similar than at high temperatures. To quantify this, we define a signal-to-noise-like quantity SNR,

where the subscript R indicates that the fit parameters are for the rearranging particles. The difference μ_R − μ is effectively always positive by construction (rearranging particles have higher S), and the difference between s_R and s is typically small but nonzero so we take the geometric mean of the two to ensure that effects of both distribution widths are approximately captured. The value of SNR will tend toward zero where the distributions are very similar and increase as they separate from one another (e.g. the softness and surrounding structure of the rearranging particles are quite different), and a positive SNR is approximately the number of standard deviations separating the rearranging and non-rearranging particles. The temperature dependence for SNR is shown in Fig. 11; for each penetrant size, SNR increases as T* is reduced, though the increase begins at lower temperatures for the smaller penetrants, which tend to be less sensitive to the glassy nature of the surrounding matrix. The data for T* < 0.4 in the σ_BB = 0.8 and 1.0σ case are noisy due to the difficulty in sampling rearranging particles below .

Fig. 11 Signal-to-noise ratios of softness distributions for rearrangement of particles as a function of temperature. Lines are drawn to guide the eye.

Next, to determine the softness-dependent barriers, we perform an Arrhenius decomposition as shown in Fig. 9(b) above to extract the S-dependent functions ΔE(S) and Σ(S). We fit the Arrhenius expression to the P_R(S) data in the range to ensure that the system is in equilibrium with the imposed temperature and that we are below the onset temperature of glassy dynamics. We then directly compare the values of absolute softness across systems with the same σ_BB and varying ε_AB because the same hyperplane was used to assess the local structure. This allows us to visualize the relative impact of tuning the interaction strength on ΔE(S) and Σ(S). Results are shown in Fig. 12.

Fig. 12 Energetic and entropic rearrangement barriers as a function of S. Lines indicate a quadratic best fit.

It becomes challenging to compare the energy barriers to rearrangement between systems, because they have been trained on different hyperplanes and therefore have different weights associated with the structure functions of the particle. However, we have shown previously that a rearrangement threshold in p_hop can translate across systems and corresponds to motion of the about the same absolute distance; definitionally, our “rearranging” particles are all traveling at least some escape distance (in this case, 0.5σ) at the time of observation. By taking advantage of this quality of p_hop, and with our knowledge of the distribution of softness values for rearranging particles, we can explore how the barriers to this rearrangement change as a function of temperature for the different systems.

We use the fits to our barrier functions and integrate the product of ΔE(S) or Σ(S) with the normalized probability distribution function of softness values for the rearranging particles at some temperature T, P(S|R,T)dS. The three quantities of interest with units of energy are then,

〈ΔF(T|R)〉 = 〈ΔE(T|R)〉 − T〈Σ(T|R)〉 (24)

which effectively yields the average energetic, entropic, and free energy barriers of the rearranging particles. Note that this is distinct from the average over all particles, since this analysis will necessarily neglect particles that have not yet rearranged and may have larger barriers.

With the energy barriers evaluated using eqn (22), we can also compare these values with those energy barriers estimated using the diffusion coefficient and relaxation times measured, allowing us to draw a connection between those measures of mobility and the hopping mechanism that gives us softness. These results are shown in Fig. 13. The energy barrier to hopping extracted from softness is of a similar magnitude to the barriers assumed from D and τ^eff_α at temperatures well above ; however, as we approach the glass transition of the system, the energy barriers from D and τ^eff_α increase rapidly, as opposed to that from softness, which increases much more slowly. The difference in the increase of the activation energy of on cooling as extracted from D and τ^eff_α is expected from the well-known breakdown of the Stokes–Einstein relationship.⁹³ In all cases, for each mode of mobility considered, the extracted energy barrier increases with both σ_BB and ε_AB at a constant temperature relative to , and the barriers estimated from the alpha-relaxation process are higher than those from diffusion.

Fig. 13 Energetic barriers to particle mobility plotted as a function of extracted from diffusion, relaxation, and softness. Panels (a)–(d) indicate different penetrant sizes σ_BB, noted on the lower right of the panels.

The increase in the barrier height from softness with increasing ε_AB is a valuable finding, as it is an extension of the trends observed in the barriers from diffusion and relaxation. Directly comparing the magnitudes of energy barriers derived from softness which were calculated from different hyperplanes (i.e. for the different particle size systems) is tempting but may be an inaccurate choice. Because the weights associated with the individual structure functions in the construction of softness are different, its physical meaning is not the same. Our proposed transformation into temperature units allows some comparison with other temperature-based trends but does not eliminate this factor. It is worthwhile, however, to examine how the barriers within a family of particle sizes change relative to one another.

Notably, the softness-dependent barriers to rearrangement show a temperature dependence that decreases with the penetrant size. For particles with a high size ratio, the barrier predicted based on softness is nearly independent of temperature. This indicates that the local structure and variations in it play less of a role in the transport of these penetrants and may point to other barriers to the hopping process that are not captured by softness. We can again draw an analogue to SCCHT, where it is demonstrated that local caging effects and non-local collective elasticity both govern transport in the high-size-ratio regime, but that the collective effects are negligible for small particles.^69,70 For our systems, this may be one potential reason why softness does not fully capture the barrier to rearrangement in larger particles, because the softness analysis does not consider any contributions due to elasticity. One might expect that the effects of elasticity would increase the energy barrier associated with particles of a given softness, as local hopping events lead to the nucleation of cage-expansion effects in the SCCHT framework. Upon cooling the matrix becomes more rigid, and a barrier that takes both the local and nonlocal parts into account should then lead to super-Arrhenius growth of the energy barrier. For the time scales where we can reasonably measure the mobility of the penetrants, we do not observe any growth of energy barriers associated with a particular softness value, and the origins of the differences in the barriers estimated by softness and those from direct mobility measures remain an open question.

Previous work on softness found that the choice of a different rearrangement threshold changed the quantitative measure of the barriers to rearrangement but did not affect the qualitative result.⁵² It is therefore relevant to explore how a change in the rearrangement definition might impact the results of our analysis. We find that changing the value of p_cut affects the results of the analysis in this section quantitatively, but not qualitatively. A lower threshold above which particles are rearranging means that the average barrier to reach that rearrangement threshold is lower, thus the energetic barrier and entropic barrier to rearrangement as a function of softness are shifted down. This is analogous to how relaxation times would be shorter if quantified using a self-intermediate scattering function evaluated over a longer wavevector k (i.e., a shorter length scale of relaxation). To give an example, when the value of the threshold is cut in half to p_cut = 0.125, we are now probing rearrangements of at least an escape distance of . Retraining a hyperplane to seek rearrangements at this value of p_cut yields a slightly lower accuracy as there are far greater distinct examples of rearranging particles. When the analysis of this section is repeated, we find the same qualitative result: the average barrier still decreases with increasing softness, and when the transformation to temperature is performed, we see the same trend observed in Fig. 13 for the softness barriers. Results of this comparison for varying σ_BB and a constant ε_AB = 1.0ε using p_cut = 0.125 are shown in Appendix 5.

Discussion

In this work, we took the well-known model of a bead-spring polymer and used molecular dynamics to simulate the transport of small particles in a manner analogous to gas molecules moving through a glassy polymer. The particles’ relative size to and interaction strength with the surrounding polymer was varied to understand how the dynamics change in response to more favorable interactions and less steric hindrance with the surrounding chains. Super-Arrhenius behavior in the self-diffusion coefficients and effective relaxation time for the particles was observed, and we were able to fit these measured quantities to an inverse-power-law relationship previously described in theory on activated penetrant dynamics. The pre-exponential factors for the inverse-power-law relationship of the diffusion coefficients was found to be a logarithmic function of the penetrant size, in accordance with SCCHT predictions for the large-size-limit ratio of activated penetrant transport. Scaling of the power-law exponent was consistent with the low-size-limit regime prediction, where the inverse diffusion coefficient scales with (T_g/T)^λ and the exponent is on the order of λ = 2. While the discrepancy in penetrant size regimes suggests that the systems we have modeled here lie in the intermediate size region, the trends observed in these fits are qualitatively consistent with the predictions from theory and with the kind of scaling observed in experiments.⁷⁰ We observed a relationship between the penetrant-to-matrix relaxation time ratio and the inverse temperature that is reminiscent of predictions from SSCHT in a homopolymer liquid indicating a dynamical decoupling in regions near the glass transition.⁷¹

To further compare the limited data we have in the supercooled region with theoretical predictions, it is relevant to consider fitting the logarithm of both the inverse penetrant diffusion coefficients and the penetrant-to-polymer relaxation time ratios to the aforementioned power-law scaling. Shifting the window of temperature data that we fit to the inverse-power-law relationship from the 0.45 < T* < 0.75 previously used in the text to 0.40 < T* < 0.55, closer to and just at the point of glass transition, we can find the relevant values of λ and K and plot them in Fig. S2.6 for the diffusion coefficients and Fig. S2.7 for the relaxation time ratios. The values of the scaling exponent λ increase by about 50% in the case of diffusion of the largest-size penetrants with σ_BB = 1.0σ, increasing from a range of 2.0 < λ < 2.2 to a range of 2.8 < λ < 3.1. While this is not exactly the 3λ scaling of the theoretical prediction in the high-size-ratio limit, the increase demonstrates that close to the supercooled glassy regime the dynamics we observe in simulations follow trends not unlike the theoretical predictions. In the case of the relaxation time ratios, we see a similar phenomenon where the scaling exponents in the high-size-ratio limit increase from a range of 2.2 < λ < 2.3 to a range of 3.3 < λ < 3.5. Further simulations sampling a broader range of size ratios and at lower particle concentrations to more accurately study a dilute limit, as well as experimental studies in the deeply supercooled and glassy regimes, are warranted to make additional quantitative comparison with established theory.

Using the framework of softness, previously applied to other glassy systems to study dynamics, we investigated the relationship between the local structure around the tracer particles in a glassy matrix and their mobility, extracting softness-dependent energetic and entropic barriers to the rearrangement process. These barriers, known to be approximately linear in systems where all particles are similar in size, become nonlinear as the particle size of the tracer shrinks, and the barriers increase in magnitude with the increasing depth of the Lennard-Jones potential well between tracer and polymer. For systems where we delineate between rearranging and non-rearranging penetrants with the same set of local structure functions and training parameters, the accuracy of the hyperplane training decreases with increasing penetrant size. Typical accuracies reported in other simulations are on the order of 70%–90% for either homopolymers^56–59 or Kob–Anderson glass formers^52,54,60 where the binary species have similar dynamics. The decrease in the ability of softness, a proxy for the local structure, to accurately predict rearrangements at the larger penetrant-to-matrix size ratios gives rise to an interesting parallel with theory. For high-size limit penetrants it is known that local cage barriers (presumably captured by softness) and collective matrix elasticity (not captured by softness) both play a role in determining mobility, whereas in low-size-limit systems the collective elasticity becomes significantly less important.^69,70 The increased role of collective elasticity in determining mobility in the high-size limit may explain why the changes in mobility are not adequately captured by softness, which is a local measure.

Comparison between systems where softness is trained with different hyperplanes has remained a challenge in drawing similarities between the machine-learned quantities and more well-understood, experimentally-accessible properties.⁹⁴ We propose a simple change-of-variables method to transform the softness-dependent barriers to a temperature dependence by integrating over the distribution of softness values at each temperature, allowing for an examination of the barriers relative to the glass transition. We find that, consistent with other measures examined in this work, the energetic barrier to rearrangement for particles of a given size σ_BB increases with polymer–penetrant interaction strength ε_AB. For the different tracer sizes studied, each shows that the barrier height increases as the glass transition temperature is approached, with the rate of increase slowing as opposed to the super-Arrhenius behavior we have seen in the barriers predicted from diffusion coefficients and relaxation times. As we have suggested above, this somewhat unspectacular finding may be indicative of the phenomenon predicted by SCCHT which notes that local caging effects play the dominant role in determining mobility for small penetrants while large penetrants are also affected by collective effects in the matrix. To further explore the phenomena discussed here, future work should explore a greater array of particle sizes and interaction strengths to identify transitions between the low- and high-size-ratio regimes of penetrant behavior to enable a deeper comparison with theory.

Finally, computational methods we have introduced here may be used to extend softness methods to more complex systems which it has not previously been applied to. For purposes of exploring model gas separation systems, limiting use to the low-size-ratio regimes would likely ensure that softness can accurately capture rearrangements in systems with penetrants and a matrix that have decoupled dynamics. One possible study could involve multicomponent diffusion in a polymer matrix with two penetrant species, examining the effects of changing both the penetrant–polymer interactions and the penetrant–penetrant interaction. This could provide some additional insights into the fundamental physics of competitive penetrant transport in glassy systems and create a direct comparison to the real-world application of the membrane-based separation of gas mixtures.

Author contributions

S. J. L.: data curation, formal analysis, investigation, resources, visualization, writing – original draft, writing – review and editing. R. A. R.: conceptualization, funding acquisition, methodology, project administration, resources, supervision, writing – review and editing.

Conflicts of interest

The authors have no conflicts of interest to disclose.

Data availability

Appendices referenced in this article are available as part of its supplementary information (SI). These include S1 (Quantifying glassy rearrangements), S2 (Tracer diffusion and relaxation behavior), S3 (Sensitivity analysis of hyperplane accuracy to p_cut,R and p_cut,NR), S4 (Tracer softness training histograms and Arrhenius plots for different systems), and S5 (Impact of changing rearrangement cutoff parameter in softness analysis). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00837a.

Data for this article, including volume vs. temperature data from cooling simulations, mean-squared displacement, self-intermediate scattering function, and softness histogram data, are available at a GitHub repository hosted at https://github.com/sjlayding/small_particle_transport_SI.git.

Acknowledgements

The authors acknowledge the generous support of both the University of Pennsylvania Ashton Fellowship and NSF Research Traineeship Award #2152205. Simulations were performed using the Stampede2, Stampede3, and Ranch resources at the Texas Advanced Computing Center (TACC) through NSF XSEDE⁹⁵ allocation TG-DMR150034 and NSF ACCESS⁹⁶ allocation TG-CHM230003. Additionally, S. J. L. would like to thank Dr. Ian Graham (University of Pennsylvania Department of Physics and Astronomy) for assistance with developing the softness calculation tools and Prof. Zahra Fakhraai (University of Pennsylvania Department of Chemistry) for many useful discussions on glass-forming materials. No generative AI tools were used in the preparation of this manuscript.

References

1. D. S. Sholl and R. P. Lively, Seven chemical separations to change the world, Nature, 2016, 532, 435–437.
2. IEA. Net Zero by 2050. https://www.iea.org/reports/net-zero-by-2050 ( 2021).
3. IEA. Net Zero Roadmap: A Global Pathway to Keep the 1.5 °C Goal in Reach. https://www.iea.org/reports/net-zero-roadmap-a-global-pathway-to-keep-the-15-0c-goal-in-reach ( 2023).
4. P. Friedlingstein, et al., Global Carbon Budget 2023, Earth Syst. Sci. Data, 2023, 15, 5301–5369.
5. K. Calvin et al. IPCC, 2023: Climate Change 2023: Synthesis Report. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Core Writing Team, H. Lee and J. Romero (Eds.)]. IPCC, Geneva, Switzerland. https://www.ipcc.ch/report/ar6/syr/ ( 2023) DOI:10.59327/IPCC/AR6-9789291691647.
6. S. J. Layding and H. S. Caram, Comparing the energy requirements for amine absorption processes for CO2 capture from flue gases calculated using a minimal information simple thermodynamic model with experimental and detailed simulation results, Int. J. Greenhouse Gas Control, 2022, 115, 103595.
7. S. J. Layding and H. S. Caram, A thermodynamic approach to energy requirements for CO₂ capture and a comparison between the minimum energy for liquid and solid sorbent processes, Int. J. Greenhouse Gas Control, 2024, 137, 104227.
8. C. Chao, Y. Deng, R. Dewil, J. Baeyens and X. Fan, Post-combustion carbon capture, Renewable Sustainable Energy Rev., 2021, 138, 110490.
9. N. McQueen, et al., A review of direct air capture (DAC): scaling up commercial technologies and innovating for the future, Prog. Energy, 2021, 3, 032001.
10. J. Young, et al., The cost of direct air capture and storage can be reduced via strategic deployment but is unlikely to fall below stated cost targets, One Earth, 2023, 6, 899–917.
11. L. M. Robeson, Correlation of separation factor versus permeability for polymeric membranes, J. Membr. Sci., 1991, 62, 165–185.
12. L. M. Robeson, Polymer membranes for gas separation, Curr. Opin. Solid State Mater. Sci., 1999, 4, 549–552.
13. L. M. Robeson, The upper bound revisited, J. Membr. Sci., 2008, 320, 390–400.
14. L. M. Robeson, Q. Liu, B. D. Freeman and D. R. Paul, Comparison of transport properties of rubbery and glassy polymers and the relevance to the upper bound relationship, J. Membr. Sci., 2015, 476, 421–431.
15. H. W. H. Lai, et al., Hydrocarbon ladder polymers with ultrahigh permselectivity for membrane gas separations, Science, 2022, 375, 1390–1392.
16. T. H. Lee, M. Balcik, W.-N. Wu, I. Pinnau and Z. P. Smith, Dual-phase microporous polymer nanofilms by interfacial polymerization for ultrafast molecular separation, Sci. Adv., 2024, 10, eadp6666.
17. Z. Yang, et al., Surpassing Robeson Upper Limit for CO2/N2 Separation with Fluorinated Carbon Molecular Sieve Membranes, Chem, 2020, 6, 631–645.
18. G. M. Iyer, L. Liu and C. Zhang, Hydrocarbon separations by glassy polymer membranes, J. Polym. Sci., 2020, 58, 2482–2517.
19. F. H. Stillinger, A Topographic View of Supercooled Liquids and Glass Formation, Science, 1995, 267, 1935–1939.
20. P. G. Debenedetti and F. H. Stillinger, Supercooled liquids and the glass transition, Nature, 2001, 410, 259–267.
21. C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan and S. W. Martin, Relaxation in glassforming liquids and amorphous solids, J. Appl. Phys., 2000, 88, 3113–3157.
22. L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, Rev. Mod. Phys., 2011, 83, 587–645.
23. G. Biroli and J. P. Garrahan, Perspective: The glass transition, J. Chem. Phys., 2013, 138, 12A301.
24. L. Berthier and M. D. Ediger, Facets of glass physics, Phys. Today, 2016, 69, 40–46.
25. V. Bapst, et al., Unveiling the predictive power of static structure in glassy systems, Nat. Phys., 2020, 16, 448–454.
26. J. Langer, The mysterious glass transition, Phys. Today, 2007, 60, 8–9.
27. W. F. Drayer and D. S. Simmons, Is the Molecular Weight Dependence of the Glass Transition Temperature Driven by a Chain End Effect?, Macromolecules, 2024, 57, 5589–5597.
28. D. L. Baker, M. Reynolds, R. Masurel, P. D. Olmsted and J. Mattsson, Cooperative Intramolecular Dynamics Control the Chain-Length-Dependent Glass Transition in Polymers, Phys. Rev. X, 2022, 12, 021047.
29. J. Yang, L. Tao, J. He, J. R. McCutcheon and Y. Li, Machine learning enables interpretable discovery of innovative polymers for gas separation membranes, Sci. Adv., 2022, 8, eabn9545.
30. M. D. Ediger, Spatially Heterogeneous Dynamics in Supercooled Liquids, Annu. Rev. Phys. Chem., 2000, 51, 99–128.
31. M. T. Cicerone and M. D. Ediger, Photobleaching technique for measuring ultraslow reorientation near and below the glass transition: tetracene in o-terphenyl, J. Phys. Chem., 1993, 97, 10489–10497.
32. L. Liu, et al., Quasielastic and inelastic neutron scattering investigation of fragile-to-strong crossover in deeply supercooled water confined in nanoporous silica matrices, J. Phys.: Condens. Matter, 2006, 18, S2261–S2284.
33. N. Vergadou and D. N. Theodorou, Molecular Modeling Investigations of Sorption and Diffusion of Small Molecules in Glassy Polymers, Membranes, 2019, 9, 98.
34. K. Zhang, D. Meng, F. Müller-Plathe and S. K. Kumar, Coarse-grained molecular dynamics simulation of activated penetrant transport in glassy polymers, Soft Matter, 2018, 14, 440–447.
35. M. L. Greenfield and D. N. Theodorou, Coarse-Grained Molecular Simulation of Penetrant Diffusion in a Glassy Polymer Using Reverse and Kinetic Monte Carlo, Macromolecules, 2001, 34, 8541–8553.
36. K. Zhang and S. K. Kumar, Molecular Simulations of Solute Transport in Polymer Melts, ACS Macro Lett., 2017, 6, 864–868.
37. D. Meng, K. Zhang and S. K. Kumar, Size-dependent penetrant diffusion in polymer glasses, Soft Matter, 2018, 14, 4226–4230.
38. M. Brownell, A. L. Frischknecht and M. A. Wilson, Subdiffusive High-Pressure Hydrogen Gas Dynamics in Elastomers, Macromolecules, 2022, 55, 3788–3800.
39. M. Balçık, S. B. Tantekin-Ersolmaz, I. Pinnau and M. G. Ahunbay, CO₂/CH₄ mixed-gas separation in PIM-1 at high pressures: Bridging atomistic simulations with process modeling, J. Membr. Sci., 2021, 640, 119838.
40. G. S. Larsen, P. Lin, K. E. Hart and C. M. Colina, Molecular Simulations of PIM-1-like Polymers of Intrinsic Microporosity, Macromolecules, 2011, 44, 6944–6951.
41. T. R. Cuthbert, N. J. Wagner and M. E. Paulaitis, Molecular Simulation of Glassy Polystyrene: Size Effects on Gas Solubilities, Macromolecules, 1997, 30, 3058–3065.
42. J. W. Haus and K. W. Kehr, Diffusion in regular and disordered lattices, Phys. Rep., 1987, 150, 263–406.
43. S. Havlin and D. Ben-Avraham, Diffusion in disordered media, Adv. Phys., 1987, 36, 695–798.
44. H. Takeuchi, A jump motion of small molecules in glassy polymers: A molecular dynamics simulation, J. Chem. Phys., 1990, 93, 2062–2067.
45. F. Müller-Plathe, L. Laaksonen and W. F. Van Gunsteren, Cooperative effects in the transport of small molecules through an amorphous polymer matrix, J. Mol. Graphics, 1993, 11, 118–120.
46. M. L. Greenfield and D. N. Theodorou, Coupling of Penetrant and Polymer Motions During Small-Molecule Diffusion In a Glassy Polymer, Mol. Simul., 1997, 19, 329–361.
47. I. Cozmuta, M. Blanco and W. A. Goddard, Gas Sorption and Barrier Properties of Polymeric Membranes from Molecular Dynamics and Monte Carlo Simulations, J. Phys. Chem. B, 2007, 111, 3151–3166.
48. G. Dorenbos and K. Morohoshi, Modeling gas permeation through membranes by kinetic Monte Carlo: Applications to H2, O2, and N2 in hydrated Nafion®, J. Chem. Phys., 2011, 134, 044133.
49. D. Paschek and R. Krishna, Diffusion of Binary Mixtures in Zeolites: Kinetic Monte Carlo versus Molecular Dynamics Simulations, Langmuir, 2001, 17, 247–254.
50. D. Richard, et al., Predicting plasticity in disordered solids from structural indicators, Phys. Rev. Mater., 2020, 4, 113609.
51. E. D. Cubuk, et al., Identifying structural flow defects in disordered solids using machine-learning methods, Phys. Rev. Lett., 2015, 114, 108001.
52. S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras and A. J. Liu, A structural approach to relaxation in glassy liquids, Nat. Phys., 2016, 12, 469–472.
53. E. D. Cubuk, et al., Structure-property relationships from universal signatures of plasticity in disordered solids, Science, 2017, 358, 1033–1037.
54. S. S. Schoenholz, E. D. Cubuk, E. Kaxiras and A. J. Liu, Relationship between local structure and relaxation in out-of-equilibrium glassy systems, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 263–267.
55. D. M. Sussman, S. S. Schoenholz, E. D. Cubuk and A. J. Liu, Disconnecting structure and dynamics in glassy thin films, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 10601–10605.
56. R. J. S. Ivancic and R. A. Riggleman, Identifying structural signatures of shear banding in model polymer nanopillars, Soft Matter, 2019, 15, 4548–4561.
57. E. Yang, R. J. S. Ivancic, E. Y. Lin and R. A. Riggleman, Effect of polymer–nanoparticle interaction on strain localization in polymer nanopillars, Soft Matter, 2020, 16, 8639–8646.
58. E. Yang and R. A. Riggleman, Role of Local Structure in the Enhanced Dynamics of Deformed Glasses, Phys. Rev. Lett., 2022, 128, 097801.
59. E. Yang, et al., Understanding creep suppression mechanisms in polymer nanocomposites through machine learning, Soft Matter, 2023, 19, 7580–7590.
60. I. R. Graham, P. E. Arratia and R. A. Riggleman, Exploring the relationship between softness and excess entropy in glass-forming systems, J. Chem. Phys., 2023, 158, 204504.
61. F. W. Starr, S. Sastry, J. F. Douglas and S. C. Glotzer, What Do We Learn from the Local Geometry of Glass-Forming Liquids?, Phys. Rev. Lett., 2002, 89, 125501.
62. A. Widmer-Cooper and P. Harrowell, On the relationship between structure and dynamics in a supercooled liquid, J. Phys.: Condens. Matter, 2005, 17, S4025–S4034.
63. R. A. Riggleman, J. F. Douglas and J. J. De Pablo, Tuning polymer melt fragility with antiplasticizer additives, J. Chem. Phys., 2007, 126, 234903.
64. W. Kob and H. C. Andersen, Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 4626–4641.
65. G. S. Grest and K. Kremer, Molecular dynamics simulation for polymers in the presence of a heat bath, Phys. Rev. A:At., Mol., Opt. Phys., 1986, 33, 3628–3631.
66. K. Kremer and G. S. Grest, Dynamics of entangled linear polymer melts: A molecular-dynamics simulation, J. Chem. Phys., 1990, 92, 5057–5086.
67. R. Zhang and K. S. Schweizer, Statistical Mechanical Theory of Penetrant Diffusion in Polymer Melts and Glasses, Macromolecules, 2016, 49, 5727–5739.
68. R. Zhang and K. S. Schweizer, Correlated matrix-fluctuation-mediated activated transport of dilute penetrants in glass-forming liquids and suspensions, J. Chem. Phys., 2017, 146, 194906.
69. B. Mei and K. S. Schweizer, Activated penetrant dynamics in glass forming liquids: size effects, decoupling, slaving, collective elasticity and correlation with matrix compressibility, Soft Matter, 2021, 17, 2624–2639.
70. B. Mei, G. S. Sheridan, C. M. Evans and K. S. Schweizer, Elucidation of the physical factors that control activated transport of penetrants in chemically complex glass-forming liquids, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2210094119.
71. B. Mei and K. S. Schweizer, Theory of the Effects of Specific Attractions and Chain Connectivity on the Activated Dynamics and Selective Transport of Penetrants in Polymer Melts, Macromolecules, 2022, 55, 9134–9151.
72. J. E. Jones, On the determination of molecular fields.—I. From the variation of the viscosity of a gas with temperature, Proc. R. Soc. London, Ser. A, 1924, 106, 441–462.
73. J. E. Jones, On the determination of molecular fields. —II. From the equation of state of a gas, Proc. R. Soc. London, Ser. A, 1924, 106, 463–477.
74. H. A. Lorentz, Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase, Ann. Phys., 1881, 248, 127–136.
75. J. D. Weeks, D. Chandler and H. C. Andersen, Role of repulsive forces in determining the equilibrium structure of simple liquids, J. Chem. Phys., 1971, 54, 5237–5247.
76. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
77. A. P. Thompson, et al., LAMMPS – a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
78. R. Auhl, R. Everaers, G. S. Grest, K. Kremer and S. J. Plimpton, Equilibration of long chain polymer melts in computer simulations, J. Chem. Phys., 2003, 119, 12718–12728.
79. Y. R. Sliozberg and J. W. Andzelm, Fast protocol for equilibration of entangled and branched polymer chains, Chem. Phys. Lett., 2012, 523, 139–143.
80. D. J. Evans and B. L. Holian, The Nose-Hoover thermostat, J. Chem. Phys., 1985, 83, 4069–4074.
81. G. J. Martyna, D. J. Tobias and M. L. Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys., 1994, 101, 4177–4189.
82. K. Dalnoki-Veress, J. A. Forrest, C. Murray, C. Gigault and J. R. Dutcher, Molecular weight dependence of reductions in the glass transition temperature of thin, freely standing polymer films, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2001, 63, 031801.
83. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University PressOxford, 2017 DOI:10.1093/oso/9780198803195.001.0001.
84. E. J. Maginn, R. A. Messerly, D. J. Carlson, D. R. Roe and J. R. Elliot, Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0], Living J. Comput. Mol. Sci., 2018, 1(1), 6324.
85. J. Behler and M. Parrinello, Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces, Phys. Rev. Lett., 2007, 98, 146401.
86. A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer and D. Chandler, Excitations Are Localized and Relaxation Is Hierarchical in Glass-Forming Liquids, Phys. Rev. X, 2011, 1, 021013.
87. A. Smessaert and J. Rottler, Distribution of local relaxation events in an aging three-dimensional glass: Spatiotemporal correlation and dynamical heterogeneity, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 88, 022314.
88. D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, 2023 DOI:10.1016/C2009-0-63921-0.
89. P. G. Debenedetti, Metastable Liquids: Concepts and Principles, Princeton University Press, 2020 DOI:10.2307/j.ctv10crfs5.
90. R. Kohlrausch, Theorie des elektrischen Rückstandes in der Leidener Flasche, Ann. Phys., 1854, 167, 56–82.
91. G. Williams and D. C. Watts, Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function, Trans. Faraday Soc., 1970, 66, 80.
92. J. S. Shaffer, Effects of chain topology on polymer dynamics: Configurational relaxation in polymer melts, J. Chem. Phys., 1995, 103, 761–772.
93. M. T. Cicerone and M. D. Ediger, Enhanced translation of probe molecules in supercooled o-terphenyl: Signature of spatially heterogeneous dynamics?, J. Chem. Phys., 1996, 104, 7210–7218.
94. A. R. Moore, Effects of Inter-Molecular and Intra-Molecular Factors on the Properties of Simulated Glasses, ProQuest Dissertations and Theses, University of Pennsylvania, United States, Pennsylvania, 2022.
95. J. Towns, et al., XSEDE: Accelerating Scientific Discovery, Comput. Sci. Eng., 2014, 16, 62–74.
96. T. J. Boerner, S. Deems, T. R. Furlani, S. L. Knuth and J. Towns, ACCESS: Advancing Innovation, Practice and Experience in Advanced Research Computing, ACM, New York, NY, USA, 2023, pp. 173–176 DOI:10.1145/3569951.3597559.

---