Hopping dynamics of a tracer particle confined in a fluctuating lattice

Abstract

We investigate the hopping dynamics of a tracer particle in a two-dimensional fluctuating lattice using Brownian dynamics simulations. Conventional analyses based on the mean square displacement and the self-intermediate scattering function reveal subdiffusive behavior characteristic of hopping transport, but these averaged quantities cannot resolve individual hopping events. To overcome this limitation, we apply hop functions originally developed for glassy and supercooled liquids, which enable the identification of discrete hopping events and provide microscopic insight into tracer dynamics. In this fluctuating lattice, direct tracking of tracer trajectories between interstitial sites offers a clear reference for validating hop-function definitions. We systematically evaluate different formulations by assessing their accuracy in identifying hopping events, reproducing residence-time distributions, and predicting free-energy barriers from hop-function probability distributions. Our results demonstrate that, among the tested formulations, the newly proposed h₂(t) yields the most accurate identification of hopping events and achieves excellent quantitative agreement with both direct trajectory tracking and umbrella-sampling free-energy barriers across a wide range of fluctuation regimes. These findings establish hop-function analysis as a robust and quantitatively reliable tool for characterizing tracer dynamics in confined, thermally fluctuating environments.

---

1 Introduction

Tracer diffusion in polymer networks has long served as a model system for studying transport in confined environments.^1–8 Previous studies have shown that thermal fluctuations of the network mesh critically determine whether a tracer exhibits normal diffusion or subdiffusion, particularly when the mesh size is comparable to the tracer diameter.^9–14 In this regime, tracers are effectively confined in transient cages formed by crosslinking and entanglement. Large fluctuations generate frequent and sufficiently long-lived pore openings, enabling continuous diffusive motion and resulting in normal diffusion. When fluctuations are small, pore openings become too rare or short-lived to support such motion, so tracers remain confined for extended periods and achieve long-time transport primarily through intermittent hopping between cages, leading to subdiffusive behavior.

The cage-jump picture has provided a widely used framework for understanding intermittent dynamics in glassy and supercooled systems. In this picture, a particle rattles within a transient cage formed by its neighbors and occasionally escapes through a rapid displacement, or jump, to a new cage.¹⁵ Various methods have been developed to identify such jump events and to characterize their statistical signatures and their relation to transport. These include displacement-threshold and trajectory coarse-graining approaches within the continuous-time random walk formalism,^15,16 classification of cage and jump trajectory segments based on the self-part of the van Hove correlation function,¹⁷ cage-jump models defined through hydrogen-bond rearrangement criteria,¹⁸ analysis of single-particle displacement distributions to distinguish hopping from continuous glassy dynamics,¹⁹ and recent statistical approaches for inferring characteristic jump scales while avoiding ad hoc displacement thresholds.²⁰

Building on this broader effort to resolve intermittent particle motion, the hop function has proven particularly effective for identifying discrete hopping events along individual particle trajectories.^21–26 In its conventional form, the hop function detects abrupt positional changes by comparing particle displacements from time-averaged positions in two adjacent time windows before and after a given time point. This analysis has successfully revealed molecular-level dynamical features that are otherwise obscured in conventional ensemble-averaged quantities such as the mean square displacement (MSD), the non-Gaussian parameter (NGP), and the van Hove correlation function.^21,27 Indeed, in molecular dynamics simulations of nanoparticle diffusion in polymer networks, inspection of individual particle trajectories has revealed distinct hopping signatures – abrupt jumps in the squared displacement interspersed with periods of localization – that are entirely masked in the ensemble-averaged MSD.⁷ However, such visual identification remains qualitative and does not provide a systematic criterion for detecting hopping events or extracting the associated free-energy barriers. To overcome this limitation, we apply hop-function analysis to tracer diffusion in network-like environments, where hopping dynamics similarly emerge from structural confinement and thermal fluctuations.

To this end, we perform Brownian dynamics (BD) simulations of a tracer particle in an idealized two-dimensional fluctuating lattice. In this model, lattice particles fluctuate around fixed average positions, and these thermal fluctuations transiently open and close pathways between interstitial regions. Depending on the fluctuation amplitude, the tracer particle may either move diffusively along these pathways or intermittently hop between regions, thereby mimicking transport processes in confined environments. The primary role of this reduced model is not to provide a complete physical description of hopping in polymer networks, but to establish a quantitatively validated hop-function framework for analyzing hopping-mediated transport in complex network environments. This setup offers two key advantages for this purpose. First, the tracer's position can be mapped unambiguously to a specific interstitial region at every time step, allowing transitions between neighboring interstitial regions to be directly tracked and used as a reference for evaluating different hop functions. Second, because lattice sites are fixed on average, the potential of mean force (PMF) for a hop can be calculated straightforwardly via umbrella sampling. This provides a reliable reference for validating free-energy barriers obtained from hop-function statistics via Boltzmann inversion and allows comparison with predictions from Kramers’ rate theory.^28–30

In this study, we apply and validate hop-function analyses within the fluctuating lattice model. We examine the conventional definition h(t) used in studies of glassy and supercooled liquids and introduce two variants, h₁(t) and h₂(t), tailored to interstitial-to-interstitial transitions. We compare the performance of different hop functions by examining their consistency with direct tracer tracking in terms of the number of inconsistent hopping identifications, the probability distributions of residence times, and the free-energy barriers obtained from Boltzmann inversion benchmarked against umbrella sampling. These comparisons provide a set of validation criteria, within which h₂(t) emerges as the most reliable measure, minimizing misidentifications and reproducing both residence-time statistics and free-energy barriers with high accuracy. Our analysis demonstrates that hop-function-based approaches establish a robust and quantitatively reliable tool for characterizing hopping behavior in confined, thermally fluctuating environments. This validated framework also opens the way to applying hop-function analysis to more complex network structures, such as polymer gels and biological systems, where collective network fluctuations make the identification of hopping events less straightforward.

The rest of this paper is organized as follows. We describe the simulation models and methods, including the definitions of the hop functions, in the following section. Simulation results are presented in the section “Results and discussion”, highlighting the relationship between lattice fluctuations and subdiffusive behavior and demonstrating the effectiveness of the hop functions in quantifying hopping dynamics. This work is summarized in the section “Conclusions”.

2 Methods

2.1 Simulation models

A two-dimensional fluctuating lattice, illustrated in Fig. 1(a), is used as an idealized model of a network with molecule-sized confinement. This setup allows us to investigate tracer dynamics in a geometrically well-defined environment where hopping events can be directly tracked. Lattice particles with a diameter of σ are placed at the vertices of a square lattice with a spacing of 3.7σ, which corresponds to the average mesh size, ξ. Each lattice particle is constrained independently by a harmonic potentialwhere r_i is the position vector of the i-th lattice particle, r_i,0 is its average position, and k_f is the force constant (in unit of k_BT/σ²) that controls the amplitude of thermal fluctuations. In the present reduced model, lattice particles are not connected to one another by springs or bead-spring connectors. This simplification provides a well-controlled reference system in which the positions of interstitial regions remain clearly defined while the amplitude of environmental fluctuations can be tuned systematically through k_f.

Fig. 1 (a) Schematic of a tracer particle (diameter 3σ) and a two-dimensional ideal fluctuating lattice composed of lattice particles (diameter 1σ) placed at intervals of 3.7σ. The average positions of the lattice particles are fixed in the lattice. (b) and (c) Probability distributions P(ξ) and standard deviations Δξ, of the mesh size around the average value ξ = 3.7σ for various lattice force constant k_f, indicating the scaling relation . The k_f values are given in unit of k_BT/σ².

A tracer particle with a diameter of 3σ is placed in the void space surrounded by four lattice particles, which is referred to as the interstitial region, hereafter. Interactions between the tracer and lattice particles are modeled using a modified repulsive Lennard-Jones (LJ) potential.

where r is the distance between the tracer and lattice particle, r_c = 2^1/6σ is the cut-off distance corresponding to the purely repulsive regime, and r_s = 1σ shifts the potential so that the tracer (3σ) and lattice particle (1σ) repel when their distance is less than r_s + r_c = (1 + 2^1/6)σ. The energy parameter is set to ε = 1k_BT, and periodic boundary conditions are applied in both dimensions.

In the absence of lattice fluctuations, the tracer particle remains entirely confined within a single interstitial region and cannot move between adjacent regions, as its size exceeds the surface-to-surface distance between neighboring lattice particles. By systematically varying the value of k_f in eqn (1), the amplitude of lattice mesh fluctuations, Δξ, is tuned, while the average mesh size remains fixed at ξ = 3.7σ. For large k_f values, the lattice fluctuations are small, whereas for small k_f values, the fluctuations become substantial, suggesting the formation of transient pathways. In this model, transient pathway formation between neighboring interstitial regions arises from local, statistically independent lattice fluctuations rather than from correlated network deformation or collective pore-breathing modes. When the fluctuation amplitude is small compared to the equilibrium separation, Δξ is approximately proportional to , as shown in Fig. 1(b) and (c). Therefore, all simulation results in this work are presented as functions of either k_f or .

2.2 Simulation methods

The dynamic motions of the tracer particle and the lattice particles were described by the position Langevin equation, eqn (3), in order to mimic the solvent-induced Brownian motions of particles without direct incorporation of solvent molecules.³¹ The BD simulations were performed in two spatial dimensions by numerically solving eqn (3) using GROMACS v. 4.5.4.³² To restrict the system to 2D, all particle motion along the z-direction was frozen, ensuring purely in-plane Brownian motion. For each time step Δt, the position r_iα(t) of particle i (either a tracer particle or lattice particles) in Cartesian direction α ∈ {x,y} was updated viawhere F_iα(t) is the total deterministic force acting on the particle and R_iα(t) is a random displacement with a Gaussian distribution exhibiting zero mean and variance–covariance of 〈R_iα(t)R_jβ(t)〉 = 2D⁰_iΔtδ_ijδ_αβ. The diffusion coefficient, D⁰_i, of each lattice particle was set equal to a value of D₀, which sets the time scale by defining the time unit of τ_BD = σ²/D₀. For the tracer particle, a reduced diffusion coefficient D⁰_i = D₀/3 was used to reflect its larger size. A time step of Δt = 10⁻⁴τ_BD was used for all simulations. Each simulation was run for 10⁵τ_BD (corresponding to 10⁹ time steps), with data recorded every 0.1τ_BD. For each value of k_f, forty independent simulations were performed to calculate statistical properties, amounting to a total simulation duration of 4 × 10⁶τ_BD per k_f. Data on shorter time scales were obtained from additional simulations of 1 × 10²τ_BD in duration, with data recorded every 10⁻⁴τ_BD. The total force F_i(t) was calculated from all the interactions described in eqn (1) and (2).

2.3 Average dynamic properties

Tracer dynamics were first characterized using average dynamic properties. The MSD is defined as

〈Δr²(t)〉 ≡ 〈[r(t) − r(0)]²〉,

where r(t) is the position vector of the tracer at time t. The angular brackets denote averaging over multiple time origins and over an ensemble of 40 independent simulation trajectories.

The self-part of the intermediate scattering function (SISF) is defined as

F_s(k,t) ≡ 〈e^{ik·[r(t)−r(0)]}〉.

Here, k = |k| is the magnitude of the wave vector probing displacements on the length scale 2π/k. In this work, we set k = 2π/3, corresponding to the tracer diameter. The SISF provides information on the relaxation of density fluctuations at this length scale.^33–36

To complement the MSD and SISF analyses, we further examine the non-Gaussian parameter (NGP), α₂(t), which characterizes deviations of the tracer displacement distribution from a Gaussian form. For two-dimensional diffusion, α₂(t) is defined as

α₂(t) vanishes when the displacement distribution is Gaussian and becomes positive when the distribution develops heavy tails, reflecting an enhanced contribution from large-displacement events. Thus, it provides a useful measure of non-Gaussian displacement statistics associated with heterogeneous and intermittent tracer motion, which may not be fully captured by the MSD alone.^{19,33,37–40}

2.4 Hop function

To characterize the discrete hopping events of the tracer particle, we employ a hop function, h(t), following recent studies on glassy and supercooled liquids.^21,27,28 The function is defined as:

h(t) = (〈[r(t) − 〈r(t)〉_B]²〉_A〈[r(t) − 〈r(t)〉_A]²〉_B)^1/2 (5)

Here, r(t) is the position vector of the tracer particle at time t. The angular brackets denote time averages over two adjacent time windows, A and B. These averages are defined as and , representing the backward- and forward-time averages, respectively. The first term under the square root in eqn (5) quantifies the variance of the tracer's displacement from its forward-time average, evaluated over the preceding interval. The second term similarly corresponds to the variance from the backward-time average, evaluated over the subsequent interval. This symmetric structure allows the hop function to sensitively detect abrupt positional changes, enabling the identification of hopping events.

However, the original hop function h(t) is based on a displacement variance and therefore has units of length squared. This quadratic measure makes it less straightforward to directly relate to the physical distance traveled during a hop. To provide a more intuitive and physically transparent description, we introduce two alternative functions, h₁(t) and h₂(t), both formulated to yield values in units of length (distance). h₁(t) is defined similarly to h(t) but each time-averaged position is replaced by the center of the interstitial region containing that position.

h₁(t) = (〈[r(t) − R_B(t)]²〉_A〈[r(t) − R_A(t)]²〉_B)^1/4 (6)

Here, R_A(t) and R_B(t) denote the centers of the interstitial regions containing the backward- and forward-time average positions, respectively. This construction ensures that h₁(t) becomes large only when the tracer actually moves from one interstitial region to another.

By contrast, h₂(t) adopts an even more direct approach, using a single fixed time window centered at t to identify hops:

h₂(t) = ([〈r(t)〉_C − R_C(t)]²)^1/2 (7)

Here, is the time average of X(t) over a symmetric interval of length Δt, and R_C(t) is the center of the interstitial region in which 〈r(t)〉_C is located. Thus, h₂(t) represents the distance between a tracer's time-averaged position and the nearest interstitial center. This value becomes large only when the average position is significantly displaced from that center, for instance, when the particle is midway between regions during a hop. In the present lattice model, the reference positions R_A(t), R_B(t), and R_C(t) are defined as the centers of interstitial regions, which are unambiguously determined from the underlying lattice geometry. In more complex systems such as polymer networks, where such regions are not explicitly defined, these reference positions may be approximated, for example, from the local average positions of surrounding crosslinkers or network nodes.

By using a linear distance scale, these alternative definitions allow a more direct interpretation of hop distances, providing clearer insight into the spatial extent of the tracer's motion. Specifically, h₁(t) compares the tracer position to the centers of adjacent interstitial regions, whereas h₂(t) measures the deviation of the time-averaged position from the center of its current region. Thus, hopping events can be identified based on a physically meaningful displacement scale. Although the three functions differ in scale, each provides complementary information: h(t) is sensitive to sudden fluctuations in motion, h₁(t) highlights transitions between distinct interstitial regions, and h₂(t) captures transient configurations near the boundary between regions.

It is worth noting that, unlike h(t) and h₁(t), h₂(t) is computed using a time average over a symmetric interval centered at t. At first glance, this might suggest that the cross term between the backward- and forward-time averages does not explicitly contribute to the calculation of h₂(t). However, because h₂(t) involves a squared quantity, the product of the backward and forward components implicitly contributes to its value. In fact, the integral over the full interval in 〈r(t)〉_C can be written as the sum of two parts, , so that the squared form of h₂(t) naturally generates a cross term through the multiplication of these backward and forward segments.

All definitions above involve time-averaging the tracer's position over a finite window. In h(t) and h₁(t), the averaging window of length Δt/2 is applied just before or after the time of interest, whereas in h₂(t) the average is taken over a full interval Δt centered on the time of interest. To calculate the hop functions, an appropriate value of Δt must be determined. If Δt is too large, the calculation of h(t) may miss hopping events within each time interval. Conversely, if Δt is too small, the resulting data may be noisy. To ensure adequate temporal resolution for detecting hopping events, Δt is chosen as the time at which the logarithmic derivative of the MSD, d log〈Δr²(t)〉/d log t, reaches its minimum value, following the work by Saito.²¹ Fig. S1 in the SI shows d log〈Δr²(t)〉/d log t for all values of k_f, and the corresponding Δt selected for each k_f is listed in Table 1.

Table 1 Values of Δt used to compute the time averages in the hop function definitions [eqn (5)–(7)]. Each Δt corresponds to the time at which d log〈Δr²(t)〉/d log t is minimized, as shown in Fig. S1 of the SI. These Δt values have been rounded such that the first decimal place is even, ensuring compatibility with the 0.2-time-unit sampling interval used for averaging

kf (kBT/σ^2)	4	10	20	40	60	100
Δt (τBD)	2.0	2.0	2.2	2.4	3.2	3.4

---

3 Results and discussion

Fig. 2(a) and (b) represent the MSD and SISF, respectively. The MSD results show that the diffusive behavior changes systematically with the fluctuation amplitude of the lattice, controlled by the force constant k_f. For small k_f, where lattice fluctuations are large, the MSD grows linearly with time over the entire range 10⁻¹τ_BD to 10³τ_BD, indicating normal diffusion. For large k_f, the MSD initially grows linearly at short times, then exhibits a slower growth or a plateau at intermediate times, and finally recovers a linear growth at longer times. This intermediate-time slowdown corresponds to subdiffusive behavior, where 〈Δr²(t)〉 ∼ t^α with α < 1, while α = 1 for normal diffusion. The transition from predominantly normal diffusion to pronounced subdiffusion becomes evident as k_f increases from 20 to 40k_BT/σ².

Fig. 2 (a) Mean square displacement (MSD) of the tracer particle for various values of k_f (in unit of k_BT/σ²). (b) Self-part of the intermediate scattering functions (SISF) calculated for k = 2π/3 (corresponding to the diameter of the tracer particle), showing the dynamic relaxation behavior of the tracer within the network. The color codes for the data at different k_f values in panel (b) are consistent with those in panel (a). (c) Non-Gaussian parameter α₂(t) of the tracer for different values of k_f. A peak in α₂(t) emerges at intermediate times and becomes larger and shifts to longer times as k_f increases, indicating increasingly heterogeneous and intermittent tracer motion under stronger confinement.

The SISF in Fig. 2(b), calculated for k = 2π/3 corresponding to the diameter of the tracer particle, show a consistent trend. For small k_f, F_s(k,t) decays in a single step over the full time window, consistent with the normal diffusion observed in the MSD. As k_f increases, the decay develops into a clear two-step process: a rapid initial decay at short times, followed by a pronounced plateau at intermediate times, and a final slow decay at long times. The appearance of the plateau indicates transient confinement of the tracer, with the height of the plateau increasing as k_f increases and lattice fluctuations weaken. The transition from single-step to two-step decay is particularly evident when k_f increases from 4 to 10k_BT/σ². At long times, all SISF curves decay to zero, but the onset of the final decay is progressively delayed with increasing k_f, mirroring the extended subdiffusive regime observed in the MSD.

To further characterize the heterogeneous nature of the tracer motion, we evaluated the NGP, α₂(t), for different k_f values (Fig. 2(c)). Non-Gaussian displacement statistics have been widely used in glassy and soft-matter systems to characterize deviations from Gaussian diffusion and variations in particle mobility, often reflecting intermittent or heterogeneous dynamics.^{19,33,37–40} In the present fluctuating lattice, α₂(t) develops a pronounced peak at intermediate times. As k_f increases, this peak becomes higher and shifts to longer times, mirroring the progressive lengthening of the subdiffusive regime observed in the MSD and the delayed final decay of the SISF. For k_f = 4k_BT/σ², the peak is nearly absent, whereas for k_f = 100k_BT/σ², it reaches α₂ ≈ 2.5 near t ∼ 50τ_BD. This behavior indicates that, as lattice fluctuations weaken, tracer motion becomes increasingly heterogeneous, with localized residence within interstitial regions coexisting with intermittent hopping events. Thus, the NGP provides an ensemble-level signature of hopping-mediated heterogeneity, complementing the MSD and SISF results while not directly resolving individual hopping events.

The MSD, SISF, and NGP together demonstrate the occurrence of subdiffusive and heterogeneous tracer dynamics arising from confinement and intermittent hopping. The MSD and SISF reveal the slowdown of tracer motion and delayed dynamic relaxation, while the NGP captures non-Gaussian displacement statistics associated with heterogeneous and intermittent motion. However, as ensemble-averaged quantities, these measures do not allow direct identification of individual hopping events along tracer trajectories. To resolve such microscopic transitions, we therefore evaluate hop functions within the lattice model. Because the average positions of lattice particles are fixed, tracer displacements between interstitial regions can be tracked unambiguously, allowing rigorous validation of hop functions, as shown in Fig. 3.

Fig. 3 Validation of the hop functions for accurate identification of hopping events. (a)–(c) Comparisons of h(t), h₁(t), and h₂(t) (red curves, plotted against the left y-axis) with the directly tracked tracer trajectory (black curves, plotted against the right y-axis) over the simulation period from 900τ_BD to 2100τ_BD for k_f = 40k_BT/σ². Horizontal dashed lines indicate the optimized threshold values for each hop function (see Table S1 of the SI), used to identify hopping events. (d)–(f) The number of inconsistent identifications of hopping events as a function of threshold values for h(t), h₁(t), and h₂(t). The total number of hopping events directly identified from the entire simulation trajectory (4 × 10⁶τ_BD) for k_f = 40k_BT/σ² is 31 839. The error percentage is defined as the mismatch count divided by this total, multiplied by 100. The optimized threshold values correspond to those that minimize this error percentage.

Panels (a)–(c) of Fig. 3 compare the directly tracked tracer trajectory (black curve) with the hop functions h(t), h₁(t), and h₂(t) (red curve), respectively, over the simulation period, as a representative example, from 900τ_BD to 2100τ_BD for k_f = 40k_BT/σ². (Data for other k_f values are presented in Fig. S2–S6 of the SI.) The horizontal dashed lines indicate the threshold values of the hop functions listed in Table S1 used to identify peaks in the red curves as hopping events. The black curve shows discrete jumps corresponding to individual hopping events, and most red peaks above the thresholds coincide with these jumps. This agreement indicates that the hop functions accurately capture hopping events.

The validity of the hop functions is quantitatively evaluated in terms of the accurate identifications of hopping events. By comparing the hop-function peaks above the threshold with the discrete jumps in the directly tracked tracer trajectory, we count the number of inconsistent identifications for each hop function, as shown for k_f = 40k_BT/σ² in panels (d)–(f) of Fig. 3. In these panels, the number of such mismatches is plotted for different threshold values of each hop function, enabling the determination of optimal thresholds. When each hop function was applied using its own optimal threshold value, on the order of ten thousand mismatches were identified for the original h(t) out of a total of 31 839 hopping events over the entire simulation (4 × 10⁶τ_BD), indicating the need for improvement. The alternative definitions h₁(t) and h₂(t) yield more accurate identification, with h₂(t) achieving the fewest inconsistent identifications among the three and thus serving as the most reliable metric for analyzing hopping dynamics.

Another validation of the hop functions is provided by the residence time, t_R, defined as the interval a tracer spends within an interstitial region before hopping. Fig. 4(a) presents the probability distributions of the residence time, P(t_R), obtained from h(t), h₁(t), and h₂(t), along with that from direct tracer tracking for k_f = 40k_BT/σ² (see Fig. S7 of the SI for other k_f values). In all cases, P(t_R) decays exponentially, consistent across all definitions of the hop function, which suggests that the hopping process follows a Poisson process irrespective of the k_f value.

Fig. 4 Residence time, t_R, defined as the time interval between consecutive hopping events. (a) Probability distributions of the residence time, t_R, obtained from h(t), h₁(t), and h₂(t), compared with that from direct tracer tracking. Data are shown for k_f = 40k_BT/σ²; results for other k_f values are provided in Fig. S7 of the SI. (b) Comparison of the average residence time, t̄_R, obtained from h(t), h₁(t), and h₂(t), with that from direct tracer tracking. The theoretical prediction of t̄_R from Kramers' rate equation in eqn (9) is also included for comparison; the equation and detailed discussion are provided later (after Fig. 7). Error bars, representing statistical uncertainties, are smaller than the symbols.

Despite the large number of inconsistent hopping identifications in h(t), its P(t_R) remains close to that from direct tracking. The agreement improves substantially for h₁(t) and h₂(t), whose distributions nearly coincide with the direct-tracking result. Among the three definitions, h₂(t) provides the most accurate measure of hopping dynamics, yielding far fewer mismatches in identifying hopping events and reproducing the residence-time distribution with high fidelity.

The quality of agreement also depends on the fluctuation amplitude. At large k_f ≥ 40k_BT, the agreement of P(t_R) from all hop functions with direct tracking is strong, but deviations increase as k_f decreases, as shown in Fig. S7 of the SI. At small k_f, the tracer exhibits nearly continuous diffusive motion, making segmentation of trajectories into discrete hopping events inherently ambiguous and thus accounting for the discrepancies. Nevertheless, the overall trend is clear: h₁(t) and h₂(t) consistently outperform h(t) across the full range of k_f. This improvement is also reflected in the average residence time, t̄_R, shown in Fig. 4(b), where estimates from h₁(t) and h₂(t) agree more closely with direct tracking than that from h(t), over the entire range of k_f or .

In our system, hopping events identified by direct tracking are defined by transitions between neighboring interstitial regions. Because h₁(t) and h₂(t) are constructed relative to interstitial-region centers, they can reliably detect such transitions even under small k_f conditions, where hopping becomes less distinct. In contrast, the original hop function h(t), which is based on displacement variance relative to backward/forward time-averaged positions, tends to miss clear crossings under the same conditions, leading to undercounted hops and consequently overestimated residence times.

The hop functions exhibit numerous peaks with varying amplitudes, as illustrated in Fig. 3. The probability distributions of the hop-function values, P(h), P(h₁), and P(h₂), are shown in Fig. 5(a)–(c). Regardless of the definition, the distributions are narrow and centered at small hop-function values for large k_f ≥ 40k_BT/σ², whereas they become broader and shift toward larger values for smaller k_f. This indicates that the tracer's motions are suppressed for large k_f and enhanced for small k_f. Although each hop-function value reflects the translational motion of the tracer particle, a hopping event is defined to occur when the motion is sufficiently large to escape from an interstitial region, that is, when the hop-function values exceed the threshold values h*, , and . These threshold values are indicated by vertical dashed lines in Fig. 5(a)–(c). The numerical values of these thresholds are summarized in Table S1 of the SI.

Fig. 5 Probability distributions and potential of mean force (PMF) profiles corresponding to the hop functions h(t) ((a) and (d)), h₁(t) ((b) and (e)), and h₂(t) ((c) and (f)) are presented for each applied force constant k_f = 4–100k_BT/σ². The results highlight the dependence of hopping behavior on the magnitude of the external constraint. Vertical dashed lines indicate the optimized threshold values (h*, , and ) used to identify hopping events; the numerical values are summarized in Table S1 of the SI.

Although these threshold values were obtained by optimizing the accuracy of hop-event identification against direct trajectory tracking, they can also be estimated without prior knowledge of the exact hopping events from the probability distributions of the hop functions (Fig. S8 of the SI). For h(t), P(h) consists of a Gaussian-like peak from caged thermal fluctuations and an exponential tail enriched in hopping events. The threshold, h*, is therefore identified as the onset of this exponential tail, following the approach employed in previous hop-function analyses.^21,26,28 For h₁(t), P(h₁) exhibits a bimodal structure, with the main peak corresponding to caged motion and the secondary peak corresponding to hopping events. The local minimum between these two populations provides a natural threshold, . For h₂(t), P(h₂) shows an exponential decay from caged fluctuations followed by an approximately plateau-like regime associated with hopping configurations. This plateau-like contribution over a range of h₂ values arises because a hop can occur at any instant within the time-averaging window. As a result, the time-averaged position during a hop can be distributed broadly between the initial and final interstitial centers. The threshold is set at the onset of this plateau-like regime, where the magnitude of the semi-logarithmic slope becomes small. For k_f = 40k_BT/σ², these distribution-based criteria give h* ≈ 0.96, , and , respectively, in good agreement with the direct-tracking optimized values of h* = 0.90, , and . These results indicate that the threshold values can be selected from the hop-function distributions themselves, enabling application of the method to systems where exact hopping events are not known a priori.

The PMF values were estimated from Boltzmann inversion of the hop-function statistics, i.e., −log(P), as a function of the hop-function values, h, h₁, and h₂, as shown in Fig. 5(d)–(f). Each curve was shifted so that its minimum value was set to zero, and the threshold values h*, , and are again marked with vertical dashed lines. The value at each threshold corresponds to the free-energy barrier for hopping from one interstitial region to another.

It is worth noting that the slight peak in −log(P) around h₁ ≈ 1.2–1.3, visible particularly at large k_f in Fig. 5(e), corresponds to the local minimum in P(h₁) discussed above: the discrete transition from R_A = R_B to R_A ≠ R_B during a hop makes h₁ values in this range inherently less probable, producing a barrier-like feature in the PMF representation.

The PMF was also calculated using umbrella sampling, a technique widely used in molecular simulations to obtain quantitative insight into thermodynamic and kinetic processes.^41–43 In particular, umbrella sampling has been applied to compute tracer PMFs in polymer networks, revealing how tracer–polymer interactions modulate the hopping energy barrier.⁸ For the present model system, the reaction coordinate was defined as the tracer displacement along the x-direction between neighboring interstitial regions, as illustrated in Fig. 6(a). The center-to-center distance between the two adjacent interstitial regions was divided into 75 windows with a spacing of 0.05σ, and in each window the tracer's x-position was restrained using a harmonic potential of strength k_x = 500k_BT/σ². The resulting biased distributions were then combined to reconstruct the PMF using the weighted-histogram analysis method (WHAM).⁴⁴ To restrict sampling to the intended reaction pathway, the tracer's motion along the y-direction was weakly restrained toward the mean y-position of the reaction coordinate using a harmonic potential of strength k_y = 3k_BT/σ². This restraint prevented the tracer from drifting into neighboring rows of interstitial regions while still allowing local fluctuations along the y-direction (see Fig. S9 of the SI for the y-displacements at various k_y).

Fig. 6 (a) Schematic illustration of the reaction coordinate for the calculation of the potential of mean force (PMF). (b) PMF obtained via umbrella sampling along the reaction coordinate for various k_f values. The tracer position in the y-direction is weakly restrained near the center of the interstitial regions by a harmonic potential with a force constant of 3k_BT/σ².

Fig. 6(b) shows the resulting PMF profiles. Starting from the center of one interstitial region (x = 0), the PMF increases, reaches a maximum at the boundary between regions (x = 1.85σ), and then decreases toward the center of the neighboring region (x = 3.7σ). As k_f increases (equivalently, as fluctuations weaken), the barrier height rises, reflecting stronger confinement. These umbrella-sampling PMF profiles serve as benchmarks for validating free-energy barriers obtained from hop-function statistics via Boltzmann inversion (Fig. 7). We note that the PMF profiles shift slightly depending on the choice of the y-directional force constant k_y (see Fig. S10 of the SI). In this study, we chose k_y = 3k_BT/σ², which allows mild y-directional fluctuations while preventing tracer migration into neighboring rows along the y-direction.

Fig. 7 Energy barriers E_b (in unit of k_BT) estimated using umbrella sampling (US) and hop functions, h(t), h₁(t), and h₂(t), for each force constant k_f (in unit of k_BT/σ²).

Assuming that the PMF calculated via umbrella sampling provides the correct free energy barrier, it can be used to validate the hop functions. Fig. 7 compares the free energy barriers obtained from umbrella sampling with those estimated using the Boltzmann inversion method from the probability distributions of the hop functions shown in Fig. 5(d)–(f). Among the three, the agreement between h₂(t) and the umbrella-sampling results is excellent. Except for a small value of k_f = 4k_BT/σ² (or ), the barriers are nearly identical, indicating that h₂(t) provides a highly accurate description of the hopping dynamics. The agreement for h₁(t) is moderate: deviating more at large k_f ≥ 40k_BT/σ² (or ) and becoming more accurate at smaller k_f.

In contrast, h(t) performs the least accurately, matching umbrella-sampling results only when the confinement is strong (k_f ≥ 40k_BT/σ²) and deviating substantially at small k_f. This outcome is consistent with the significant errors in residence-time estimates shown in Fig. 4, and with the fact that at small k_f the tracer undergoes nearly continuous diffusive motion, which makes segmentation into discrete hops inherently ambiguous.

Overall, our results demonstrate that h₂(t) is the most robust and reliable measure of hopping dynamics in this confined, fluctuating-lattice system, with h₁(t) also performing well under most conditions. At the same time, it is worth noting that the relative performance of the hop functions may depend on the physical system: hop functions like h₁(t) and h₂(t) are expected to work especially well when confinement creates well-defined spatial compartments, while h(t) may still provide advantages in systems where boundaries are more diffuse and hops cannot be cleanly separated.

It is also useful to clarify how the present reduced model differs from more realistic flexible polymer-network models.^{9–14,45–47} In flexible gels and polymer networks, crosslinks, bead-spring connectivity, and collective network fluctuations can strongly influence tracer transport. Recent three-dimensional simulations of flexible model gels and polymer networks have demonstrated that such effects can modify both the accessibility of transient pathways and the resulting long-time diffusivity. In particular, concerted breathing modes of a flexible network can transiently widen pores and facilitate the passage of large tracers, as highlighted in Godec et al.¹⁰ These collective modes are absent in the present uncoupled lattice by construction. Thus, transient pathway formation in our model is governed by local, statistically independent lattice fluctuations rather than by correlated pore-breathing motions.

This difference should be regarded as a deliberate reduction rather than as a complete representation of flexible gel transport. The advantage of the present model is that it provides a geometrically well-defined setting in which transitions between neighboring interstitial regions can be directly tracked. This makes it possible to benchmark hop-function definitions quantitatively. The resulting framework provides a basis for future applications to more realistic three-dimensional polymer-network models, where the influence of bead-spring connectivity and collective network fluctuations on hopping dynamics can be examined explicitly.

As an additional time-resolved measure of intermittent tracer dynamics, we also evaluated the time-local diffusion coefficient, D_t(t), also referred to as the temporal diffusion coefficient,⁴⁸ according to

where Δ is the lag time and t_w is the measurement time window. Because D_t(t) measures the local diffusivity over a finite trajectory segment, it can capture transient bursts of mobility associated with hopping events. Representative traces obtained with various values of Δ and t_w are shown in Fig. S11 of the SI. For Δ = 5τ_BD and t_w = 6τ_BD, sharp peaks in D_t(t) coincide well with abrupt jumps in the directly tracked tracer trajectory. A threshold analysis based on this condition is provided in Fig. S12, where the optimal value was determined by minimizing the number of inconsistent identifications relative to direct trajectory tracking. Although this analysis confirms that D_t(t) is useful for detecting hopping-like intermittent motion, the number of inconsistent identifications remains substantially larger than that obtained using h₂(t). We therefore use D_t(t) as a complementary indicator of dynamic heterogeneity, while h₂(t) remains the primary and more accurate criterion for identifying discrete hopping events.

In a historical paper, Kramers related the reaction rate constant to the free energy barrier.^29,30 Using Kramers' rate theory and the relation between the rate constant and the residence time, Ma et al. recently determined the free energy barrier from residence times obtained in experiments (eqn (9)).²⁸

where D⁰_t is the diffusion coefficient of a tracer particle diffusing freely in the absence of lattice particles, and are the second derivatives of U(x), the free-energy profile, at the positions of the minimum and maximum, respectively, and E_b is the free-energy barrier for hopping.

Using the PMF determined from umbrella sampling, we estimated the average residence time t̄_R using eqn (9) and compared it with values obtained from direct tracer tracking and from h(t), h₁(t), and h₂(t), as shown in Fig. 4(b). and in eqn (9) were taken from the minimum and maximum of the umbrella-sampling PMF, as shown in Fig. S13 of the SI, and D⁰_t was obtained from the MSD of a freely diffusing tracer particle.

The average residence time predicted by Kramers' theory is approximately 2.5 times larger than that obtained from simulation across the entire range of k_f, as shown in Fig. 4(b) and Fig. S13 of the SI. This quantitative deviation appears to arise mainly from the simplified nature of the theoretical model, which assumes a quasi one-dimensional energy landscape, whereas the tracer in simulation diffuses in two dimensions. In such higher-dimensional systems, the tracer can escape along multiple pathways, leading to a shorter mean residence time than that predicted by the one-dimensional model. Nevertheless, the theoretical and simulation results maintain an almost constant ratio across all k_f values, indicating that Kramers' theory correctly reproduces the relative variation of residence time with confinement strength. Thus, even though the absolute values are overestimated by roughly a factor of 2–3, the theory accurately captures the underlying scaling relationship between the residence time and the energy barrier, providing a physically consistent and predictive framework for hopping-dominated diffusion.

4 Conclusions

We investigated the hopping dynamics of a tracer particle confined within a two-dimensional fluctuating lattice using BD simulations. By systematically varying the lattice fluctuation amplitude, we revealed a clear transition from normal diffusion at small k_f (large fluctuations) to subdiffusive behavior at large k_f (small fluctuations). This transition, observed consistently in both the MSD and the SISF, highlights the central role of thermal fluctuations of the network mesh.

To resolve microscopic hopping events that are obscured in averaged measures, we applied and validated several hop functions within this controlled lattice model. Direct tracking of tracer motion between interstitial regions provided an unambiguous benchmark, showing that the newly defined h₂(t) yields the most accurate identification of hopping events, reproduces residence-time statistics with high fidelity, and delivers free-energy barriers in close agreement with umbrella sampling results.

Overall, our findings establish hop-function analysis – particularly using h₂(t) – as a robust and quantitatively reliable methodology for characterizing tracer dynamics in confined, thermally fluctuating environments. More broadly, this validated framework offers a general approach for analyzing hopping-mediated transport in complex heterogeneous systems. In many material and biomolecular systems, association and dissociation processes occur at spatially localized traps or binding sites, and the h₂(t) introduced here provides a useful descriptor for analyzing the dynamics between such localized sites across a wide range of systems.

In more realistic polymer-network models, such as flexible gels with crosslinks and bead-spring connectivity, hopping dynamics can be affected by network architecture and collective fluctuation modes.^{9–14,45–47} In such systems, the local reference positions required by h₁(t) and h₂(t) are not defined a priori as they are in the present ideal lattice, but may be approximated from structural features such as crosslinks or network nodes. Extending the present hop-function analysis to three-dimensional polymer networks, together with robust criteria for defining local reference positions and threshold values, remains an important direction for future work.

Finally, although our results demonstrate that h₂(t) performs best – and h₁(t) also works reliably – across regimes where hopping is well defined, we note that h(t) may still provide value in situations where compartment boundaries are less sharply defined, and discrete hops are inherently more difficult to identify.

Author contributions

J. S. K., S. S., and B. J. S. conceived, designed, and supervised the work. S. K. performed the computer simulations and visualization. S. K. and Y. S. performed the data analysis. S. K. and J. S. K. wrote the manuscript. S. S. and J. S. K. reviewed and edited the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Table S1 and Fig. S1–S13. See DOI: https://doi.org/10.1039/d6sm00250a.

Acknowledgements

This research was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT and MOE) under Grant No. NRF-2020R1A5A2019210, RS-2024-00338551, RS-2025-00519673, and RS-2025-16063688.

Notes and references

1. B. Amsden, Macromolecules, 1998, 31, 8382–8395.
2. L. Masaro and X. X. Zhu, Prog. Polym. Sci., 1999, 24, 731–775.
3. F. Höfling and T. Franosch, Rep. Prog. Phys., 2013, 76, 046602.
4. Z. E. Dell and K. S. Schweizer, Macromolecules, 2014, 47, 405–414.
5. L.-H. Cai, S. Panyukov and M. Rubinstein, Macromolecules, 2015, 48, 847–862.
6. K. Saalwächter and S. Seiffert, Soft Matter, 2018, 14, 1976–1991.
7. V. Sorichetti, V. Hugouvieux and W. Kob, Macromolecules, 2021, 54, 8575–8589.
8. B.-R. Zhao and B. Li, J. Chem. Phys., 2022, 157, 104901.
9. H. W. Cho, H. Kim, B. J. Sung and J. S. Kim, Polymers, 2020, 12, 2067.
10. A. Godec, M. Bauer and R. Metzler, New J. Phys., 2014, 16, 092002.
11. H. Zhou and S. B. Chen, Phys. Rev. E, 2009, 79, 021801.
12. N. Kamerlin and C. Elvingson, J. Phys.: Condens. Matter, 2016, 28, 475101.
13. P. Kumar, L. Theeyancheri, S. Chaki and R. Chakrabarti, Soft Matter, 2019, 15, 8992–9002.
14. Y. Chen, R. Ma, X. Qian, R. Zhang, X. Huang, H. Xu, M. Zhou and J. Liu, Macromolecules, 2020, 53, 4172–4184.
15. M. P. Ciamarra, R. Pastore and A. Coniglio, Soft Matter, 2016, 12, 358–366.
16. R. Pastore, T. Kikutsuji, F. Rusciano, N. Matubayasi, K. Kim and F. Greco, J. Chem. Phys., 2021, 155, 114503.
17. S. Dueby, V. Dubey and S. Daschakraborty, J. Phys. Chem. B, 2019, 123, 7178–7189.
18. T. Kikutsuji, K. Kim and N. Matubayasi, J. Chem. Phys., 2019, 150, 204502.
19. V. Sposini, C. N. Likos and M. Camargo, Soft Matter, 2023, 19, 9531–9540.
20. S. Pigolotti and S. Roldán-Vargas, Phys. Rev. E, 2025, 112, 035415.
21. S. Saito, J. Chem. Phys., 2024, 160, 194506.
22. R. Candelier, O. Dauchot and G. Biroli, Europhys. Lett., 2010, 92, 24003.
23. S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras and A. J. Liu, Nat. Phys., 2016, 12, 469–471.
24. R. Candelier, A. Widmer-Cooper, J. K. Kummerfeld, O. Dauchot, G. Biroli, P. Harrowell and D. R. Reichman, Phys. Rev. Lett., 2010, 105, 135702.
25. R. Candelier, O. Dauchot and G. Biroli, Phys. Rev. Lett., 2009, 102, 088001.
26. A. Smessaert and J. Rottler, Phys. Rev. E, 2013, 88, 022314.
27. B. Kang, J. Yu, S. Saito, J. Jang and B. J. Sung, Adv. Sci., 2025, 12, e09205.
28. X. Ma, Z. S. Davidson, T. Still, R. J. S. Ivancic, S. S. Schoenholz, A. J. Liu and A. G. Yodh, Phys. Rev. Lett., 2019, 122, 028001.
29. P. Hänggi, P. Talkner and M. Borkovec, Rev. Mod. Phys., 1990, 62, 251–341.
30. H. A. Kramers, Physica, 1940, 7, 284–304.
31. D. L. Ermak and J. A. McCammon, J. Chem. Phys., 1978, 69, 1352–1360.
32. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
33. J. Kim, C. Kim and B. J. Sung, Phys. Rev. Lett., 2013, 110, 047801.
34. W. Kob and H. C. Andersen, Phys. Rev. E, 1995, 52, 4134–4153.
35. T. Sentjabrskaja, E. Zaccarelli, C. De Michele, F. Sciortino, P. Tartaglia, T. Voigtmann, S. U. Egelhaaf and M. Laurati, Nat. Commun., 2016, 7, 11133.
36. N. Kuon, A. A. Milischuk, B. M. Ladanyi and E. Flenner, J. Chem. Phys., 2017, 146, 214501.
37. I. Y. Wong, M. L. Gardel, D. R. Reichman, E. R. Weeks, M. T. Valentine, A. R. Bausch and D. A. Weitz, Phys. Rev. Lett., 2004, 92, 178101.
38. M. Levin, G. Bel and Y. Roichman, J. Chem. Phys., 2021, 154, 144901.
39. J. M. Miotto, S. Pigolotti, A. V. Chechkin and S. Roldán-Vargas, Phys. Rev. X, 2021, 11, 031002.
40. S. Song, S. J. Park, M. Kim, J. S. Kim, B. J. Sung, S. Lee, J.-H. Kim and J. Sung, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 12733–12742.
41. J. Kästner, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 932–942.
42. N. Kim, J. H. Lee, Y. Song, J. H. Lee, G. C. Schatz and H. Hwang, J. Phys. Chem. B, 2023, 127, 6061–6072.
43. L. Zhao, J. Yuan, Y. Xing, J. Qi, P. Peng, Z. Liu and A. Zheng, Chem. Soc. Rev., 2026, 55, 210–253.
44. A. Grossfield, WHAM: the weighted histogram analysis method, https://membrane.urmc.rochester.edu/wordpress/?page_id=126, Accessed 2025, Version 2.0.10.
45. Y. Kim, S. Joo, W. K. Kim and J.-H. Jeon, Macromolecules, 2022, 55, 7136–7147.
46. P. Kumar and R. Chakrabarti, Phys. Chem. Chem. Phys., 2023, 25, 1937–1946.
47. K. Goswami, A. G. Cherstvy, A. Godec and R. Metzler, Phys. Rev. E, 2024, 110, 044609.
48. E. Yamamoto, T. Akimoto, A. Mitsutake and R. Metzler, Phys. Rev. Lett., 2021, 126, 128101.

---