Rounding of the localization transition in model porous media

Abstract

The generic mechanisms of anomalous transport in porous media are investigated by computer simulations of two-dimensional model systems. In order to bridge the gap between the strongly idealized Lorentz model and realistic models of porous media, two models of increasing complexity are considered: a cherry-pit model with hard-core correlations as well as a soft-potential model. An ideal gas of tracer particles inserted into these structures is found to exhibit anomalous transport which extends up to several decades in time. Also, the self-diffusion of the tracers becomes suppressed upon increasing the density of the systems. These phenomena are attributed to an underlying percolation transition. In the soft potential model the transition is rounded, since each tracer encounters its own critical density according to its energy. Therefore, the rounding of the transition is a generic occurrence in realistic, soft systems.

---

1 Introduction

Molecular transport in strongly heterogeneous media is fundamental for a wide range of disciplines such as molecular sieving,¹ catalysis^1–4 and ion-conductors,^5,6 but also for protein motion in the interior of “crowded” cells.^7–10 Common to all of these systems is the occurrence of a quasi-immobilized host structure which restricts transport to ramified paths through the medium.

The generic features of transport in heterogeneous media^2,11 can be elucidated by studying simplified model systems such as the Lorentz model.⁸ In its simplest variant,^12,13 a point tracer explores the space between randomly distributed hard-disk obstacles, which may overlap and are uncorrelated. Recently, Skinner et al.¹⁴ presented a colloidal realization of a two-dimensional Lorentz model, which differs from the original model with respect to the matrix structure and the interactions. In detail, the experiment uses a slightly size-disparate binary mixture of superparamagnetic colloidal spheres. The larger particle species is immobilized by the cover slides, while the smaller one serves as tracers. If an external magnetic field is applied, magnetic dipoles are induced which lead to a soft repulsion between the particles; the range of the tracer-matrix interaction thus can be tuned by the strength of the magnetic field.

A striking observation in these models is a localization transition with respect to the diffusive motion of the tracer particle. In the hard-core model, long-range transport ceases to exist as a critical obstacle density is approached.^15–20 Concomitantly, transport becomes anomalous as manifested in a non-linear, power-law growth of the mean-squared displacement, δr²(t) ∼ t^2/z with a universal dynamic exponent z. In the case of soft interactions, the experimental and simulation results¹⁴ indicate that the transition is rounded, i.e. the critical singularities are seemingly avoided. Evidence for anomalous transport has been found also in a variety of systems with a strong dynamic asymmetry, e.g. in computer simulations for alkali-doped silica melts²¹ or polymer blends with monomer-size disparity.^22,23 Recent work^24–27 suggests that this holds generically in size-disparate mixtures. Although it seems plausible to expect an analogy between these more realistic models with soft interactions and the theoretical idealization,²⁸ a direct and quantitative link is missing.

The goal of this work is to provide intermediate steps from the hard-disk idealization to more realistic systems, thereby testing the key ingredients leading to anomalous transport. For hard tracer–matrix interactions, many facets of the localization transition are well understood.^8,19,29 Most importantly, the localization transition is due to an underlying continuum percolation transition of the accessible void space,¹⁷ which is accompanied by a series of scaling laws familiar from the theory of critical phenomena of continuous phase transitions.⁴ Above a critical obstacle density, the network of the void space falls apart into a hierarchy of finite-sized pores. At criticality, the void space is a self-similar fractal in the statistical sense, which entails subdiffusion for tracers exploring these structures.^4,30 Correlations in the host matrix modify the geometry of the void space with potential implications on the critical behavior. Moreover, soft interactions smear out the boundaries of the accessible space and change the topology of the percolation network by introducing a potential energy landscape with finite barriers between the pores.

To investigate these issues, we compare simulations of two different models, which represent modifications to the original Lorentz model. We focus on two-dimensional systems which are amenable to colloid experiments. First, we introduce spatial correlations in the host matrix by using an extended tracer in frozen-in configurations of equilibrated hard disks. The resulting host structure is equivalent to the cherry-pit model.³¹ Second, we relax the assumption of hard-core repulsion between both obstacles and tracers by introducing soft interactions. Now, the host structures are generated from snapshots of an equilibrated fluid of soft particles. By this, the transition becomes rounded and we demonstrate that the rounding originates naturally from the energy distribution of the tracers. We find that an effective interaction distance can be assigned to each tracer according to its energy, thereby providing a mapping to the hard-core case.

2 Cherry-pit model

Host structures

In the cherry-pit model, the matrix of obstacles is formed by equilibrium configurations of N non-overlapping disks of diameter σ_core, packing fraction η = (N/L²)πσ²_core/4, and the centers are confined to a square of edge length L. The remaining space is explored by a “ballistic” tracer undergoing specular scattering from the obstacles, yielding trajectories R(t) which are piece-wise straight lines. The tracer particles are disks of finite diameter σ_T, contrarily to the original overlapping Lorentz model (Fig. 1a and b). For comparison to the latter, we introduce the interaction distance σ = (σ_core + σ_T)/2, with which a dimensionless control parameter, the reduced number density n* = (N/L²)σ², can be defined. The velocity of the tracer is of fixed magnitude v and defines a time scale t₀ = σ/v. Transport is controlled by variation of n* at fixed η.

Fig. 1 (a–c) Illustrations of the relevant models. (a) Lorentz model: overlapping obstacles (grey) and point tracer (red). (b) Cherry-pit model: obstacles (dark grey) and extended tracer (red). The area inaccessible to the tracer center is marked in light grey. (c) WCA-system: soft obstacles (grey) and soft tracer (red). (d) Static structure factor of the matrix in the cherry-pit model. For comparison, the matrix structure factor of the WCA-disk system (dashed line) with the effective diameter σ^BH_core is included.

The matrix configurations are generated by canonical Monte Carlo simulations, where the particles are initially placed onto a hexagonal grid. We consider systems of N = 10 044 or N = 516 468 disks in a square box with varying size assuming periodic boundary conditions in the two spatial directions. For the equilibration of the systems we combine displacement moves with cluster moves proposed by Dress and Krauth.³² In the displacement moves, a random particle i with position r_i is displaced to a new position r_i + δ, where the vector δ is randomly chosen such that |δ| < σ_core. This move is accepted according to a standard Metropolis criterion.³³

Cluster moves are applied periodically after 10 displacement moves. To this end, a pivot is selected as a random point in the system. By starting with one randomly selected disk and recursively searching for disks overlapping with the disks' mirror image with respect to the pivot, we identify a pair of disk clusters (C₁, C₂), C₁ ≠ C₂, defined as two sets of disks satisfying the following condition: When all disks in C₁ are reflected at the pivot, each of them overlaps with at least one disk in C₂, but none overlaps with disks not in C₂, and vice versa. If clusters are larger than 15 disks, the cluster move is rejected. In this manner, the clusters can be exchanged with their reflected counterparts. In the following, the Monte Carlo time is given in terms of cycles, where each cycle consists of N displacement moves and N/10 cluster moves.

Systems with packing fractions ranging from η = 0.02 to 0.90 for the small systems and η = 0.02 to 0.65 for the large systems were generated. At each value of η, the configurations were first equilibrated for at least 1000 cycles for low packing fractions and up to 50 000 cycles for high packing fractions to ensure proper equilibration, particularly for η ≲ 0.7, i.e. for packing fractions lower than the location of the fluid-to-solid transition in hard disks.^34,35 To check whether the system was sufficiently equilibrated we monitored the structure factor and the pair correlation function and compared it to the Percus–Yevick approximation.^36,37 Additionally, for systems η < 0.7, we required that particles are displaced by L/2 on average. For each equilibrated configuration, a production run was performed to yield 20 independent configurations, each of them separated by the respective equilibration time. These configurations served as matrix configurations for the tracer particle dynamics.

The structural correlations contained in the obstacle matrix are a function of the packing fraction η and are measured by the static structure factor of the obstacles,

as a function of the wave number q = |q|, see Fig. 1d. The angled brackets represent an ensemble average over the disorder, and {r_j} denote the positions of obstacle centers, j = 1, …, N. At low packing fractions the system exhibits the structure of a dilute liquid, as indicated by the low amplitude of the first diffraction peak, e.g. S(q_max) ≈ 1.3 at η = 0.26. As the packing fractions increases, the peak grows in amplitude, e.g. S(q_max) ≈ 2.7 at η = 0.6, and S(q) exhibits pronounced short-range order, indicating the structure of a dense liquid. Also included in Fig. 1d is the structure factor for the WCA system (dashed line) which agrees well with that of the cherry-pit model at η = 0.26 (see Section. 3).

Percolation threshold

For the study of the localization transition, it is crucial to precisely know the percolation threshold of the void space accessible to the tracer particle. For a given obstacle configuration, we have determined the threshold value n^*_c = Nσ_c²/L² of the reduced number density by varying the distance σ of the tracer–obstacle interaction. First, we have computed the edges of a Voronoi tessellation of the matrix using the free voro++ software library.³⁸ After removal of the edges with a distance smaller than σ to an obstacle center, the obtained network represents the volume accessible to the tracer particle.³⁹ Upon increasing σ, this connectivity network is diluted until the critical value σ_c is reached, where the residual network barely spans the entire simulation box.

While there is a unique critical density n^*_c for infinitely large systems L → ∞, at finite system sizes L the percolation thresholds of the individual obstacle configurations follow a distribution with a finite width. For decreasing system size, the mean of the distribution is shifted towards a slightly higher critical density n^*_c(L) according to¹⁹n^*_c(L) − n^*_c ∼ L^−1/ν. Additionally, the width of the distribution, which can be measured with the standard deviation δn^*_c(L) for example, scales as δn^*_c(L) ∼ L^−1/ν.

Over the full range of packing fractions 0 < η < η_hcp from the ideal gas to close packing, we have numerically determined the critical reduced density n^*_c(η) for two different system sizes, shown in Fig. 2. The critical density is calculated from the mean of the percolation distance σ_c of the obstacle configurations. The error bars are calculated with the help of the relative standard deviation of the critical distance Δ(η) = δσ_c(η)/σ_c(η) and thus allow to estimate the width of the distribution of the percolation thresholds. This gives an estimate for the percolation density, n^*_c(η) = (Nσ_c²(η)/L²)[1 ± Δ(η)]². We confirmed exemplarily for the case η = 0.26 that the relative standard deviation Δ(η) is indeed a good approximation to the distribution width δn^*_c(η,L) of the critical distance, as we observed that increasing the number of independent configurations up to 300 did not modify Δ(η) within the specified precision.

Fig. 2 Critical reduced density n^*_c for the cherry-pit model as a function of the packing fraction η of the obstacle cores for two different system sizes, 10⁴ and ≈5 × 10⁵ obstacles. The overlapping Lorentz model corresponds to η = 0. For comparison: WCA system at the effective packing fraction η = 0.225 and n^*_c ≈ 0.272.

The overlapping Lorentz model corresponds to η = 0, here the percolation threshold is known accurately:⁴⁰n^*_c(0) = 0.3590810 ± 0.0000006 for the infinitely large system. For packing fractions η ≲ 0.45, the percolation threshold decreases from this value and can be fitted with a shifted exponential function f(η) = a exp(−bη) + c. For 0.45 ≲ η ≲ 0.9, the percolation threshold is growing, with a “shoulder” around η ≈ 0.7 indicating the 2D melting transition. At , the system displays a hexagonal closed-packed structure and therefore .

For the following study of the tracer dynamics and how it is affected by the structural correlations contained in the matrix, we consider in detail η = 0.26 and η = 0.60 with percolation thresholds n^*_c = 0.262(1 ± 2 × 10⁻³)² and n^*_c = 0.2442(1 ± 7 × 10⁻⁴)² respectively. At η = 0.26, the structure factor closely resembles that of the WCA system discussed later on (Fig. 1d).

Tracer dynamics

The mean-squared displacements (MSDs) δr²(t) = 〈|R(t) − R(0)|²〉 were obtained as time- and ensemble-average over 20 obstacle configurations containing N = 516 468 obstacles each in production runs up to times 10⁹t₀. Each obstacle configuration was probed by at least 8 tracers, while 32 tracers were used close to n^*_c. For each tracer a random point in the void space was chosen as the initial position. Moving time averages were calculated efficiently with an “order-n” algorithm.^41,42

The MSD are shown in Fig. 3 for reduced densities close to the critical one such that the interaction distance σ is changed in geometric progression with the standard deviation Δ as basic scale. For both values of η, the localization transition is evident and qualitatively similar to the overlapping 2D Lorentz model:^18,20,43 for times t longer than a certain crossover time scale t_x, the MSD either grows diffusively, δr²(t) ≃ 4Dt for n* < n^*_c with diffusion coefficient D, or saturates, δr²(t) ≃ [small script l]² for n* > n^*_c with localization length [small script l]. The transport is highly heterogeneous in space: a fraction of tracers is confined to finite pores, which exist at all densities and have a broad distribution of sizes near n^*_c, but only tracers on the spanning cluster contribute to long-range transport (for a snapshot of the Lorentz gas with finite pores, see Fig. 1 of ref. 17). For n* > n^*_c, the spanning cluster disappears and all tracers are confined. This implies that the localization length [small script l] is the root-mean-square size of the finite clusters. As the critical point is approached, n* → n^*_c, a sub-diffusive regime emerges in a growing time window,

δr²(t) ∼ t^2/z, t₀ ≪ t ≪ t_x. (2)

Fig. 3 Mean-squared displacements in the cherry-pit model. Reduced density n* is varied around the critical density n^* = n^*_c(1 + kΔ)², k = 0, ±2⁻¹, …±2⁶ in geometric progression. Obstacle packing fractions (a) η = 0.26, n^*_c = 0.262 with standard deviation of the percolation threshold Δ = 2 × 10⁻³ and (b) η = 0.6, n^*_c = 0.2442 and Δ = 7 × 10⁻⁴. Data below n^*_c fan out towards the upper left, the ones above n^*_c towards the lower right. The solid line indicates a power-law ∝ t^2/z with the dynamic exponent of the Lorentz model z = 3.036. The horizontal dashed line indicates the size of the simulation box.

The exponent z is believed to be universal for particle transport on 2D percolation clusters^44,45 and may be considered the fundamental dynamic exponent of the problem. It was estimated to z = 3.036 ± 0.001 from studies of the conductivity of random resistor networks,⁴⁶ random walkers on percolation lattices³⁰ and in the overlapping 2D Lorentz model.^20,47 The value was confirmed only recently also for the overlapping 2D Lorentz model with ballistic tracers.⁴³ Our data for the MSD in the cherry-pit model suggest anomalous transport with effective exponents slightly lower than the universal one (Fig. 3).

A more thorough test of the value of the anomalous exponent can be achieved with the local exponent γ(t) of the MSD defined as

At short times, the exponent γ(t) decays quickly from its initial value 2 for ballistic motion due to the scattering from the obstacles. At the lowest densities, γ(t) rapidly converges to 1 corresponding to the linear increase of the MSD. At high densities the local exponents converge to 0, reflecting the localization. Values corresponding to anomalous diffusion are found close to the transition, yet the exponent found here seems to underestimate the universal value, obeyed by a random walker on percolation lattices.³⁰ However, the dynamics is extremely sensitive to the density near the percolation threshold. Probably, we slightly overestimate the critical density n^*_c(L), as obtained from the static determination of the percolation threshold (see discussion above). In view of a small overestimation of n^*_c(L), our data is consistent with the expected exponent z, since the local exponent becomes compatible with the universal value of z at the lower end of the error margin for the percolation threshold, see blue lines for k = −2⁰ in Fig. 4.

Fig. 4 Local exponents of the mean-squared displacements of Fig. 3 of the cherry-pit model for the obstacle packing fractions (a) η = 0.26 and (b) η = 0.6. Reduced density increases from top to bottom. The horizontal line indicates the anomalous exponent 2/z with z = 3.036 of the Lorentz model. The shaded areas correspond to one standard deviation Δ in the interaction distance σ.

In the approach to the percolation threshold, long-time diffusion decreases such that it vanishes at the critical point. Scaling arguments^4,48 predict a power-law singularity,

D(ε↑0) ∼ (−ε)^μ, ε = (n* − n^*_c)/n^*_c, (4)

with the exponent fixed by μ = (z − 2)(ν − β/2). The universal exponents ν and β characterize the underlying geometry of the spanning cluster, namely its correlation length ξ ∼ |ε|^−ν (the scale up to which it is self-similar) and its weight P_∞ ∼ |ε|^β. For two-dimensional standard percolation, ν = 4/3 and β = 5/36 hold exactly,⁴ and one computes μ = 1.309 ± 0.002 from the above value of z.

The diffusion coefficients obtained from the long-time behavior of the MSDs are reduced for larger tracers (at fixed η, σ_core), and the suppression of diffusion is compatible with the percolation threshold as determined above, see Fig. 5. Plotting D^1/μvs. n* (see inset) rectifies the critical law, eqn (4), and would yield a straight line if the power law was an accurate description over the full range n* < n^*_c. From the data one infers that the scaling law becomes valid for ε ≲ 0.05, which is similar to the situation in the overlapping case.⁴³

Fig. 5 Diffusion coefficients for the cherry-pit model at packing fractions η = 0.26 and 0.6 as a function of the reduced density n* divided by the respective percolation thresholds n^*_c(η) = 0.262 and 0.2442. Open symbols mark data points which were obtained at densities where the MSD had not quite become diffusive and thus potentially overestimate D. Inset: rectification plot of the same data testing eqn (4) with conductivity exponent μ = 1.309. The straight line serves as guide to the eye.

3 WCA-disk systems

Host structure

Next, we move towards possible experimental realizations and relax the idealization of a hard-core exclusion between the particles, considering soft interactions. To serve as frozen host structures, we take snapshots of an equilibrated liquid of polydisperse particles at moderate densities. The particles interact via a Weeks–Chandler–Andersen (WCA) potential,⁴⁹ which is a truncated and shifted Lennard-Jones potential such that the interaction is purely repulsive.with a cutoff r_cut = 2^1/6σ_αβ. To avoid crystallization, a polydisperse mixture is necessary. To this end, the diameters of the N particles are chosen to be additive, σ_αβ = (σ_α + σ_β)/2, and are taken equidistantly from an interval, σ_α = (0.85 + 0.3α/N)σ^WCA_core with α, β = 1, …, N. The units of length and energy are fixed by σ^WCA_core and ε^WCA_core, respectively. The temperature is set to k_BT/ε^WCA_core = 1.0. To improve numerical stability the potential is multiplied with a smoothing function Ψ(r) = (r − r_cut)⁴/[h⁴ + (r − r_cut)⁴] with the width h = 0.005σ^WCA_core. This smoothing function provides continuity of the force at the cutoff r_cut. As a result, we have not encountered any problems with a drift of the total energy in all the microcanonical runs of the WCA systems.

The particle configurations are equilibrated using a simplified Andersen thermostat⁵⁰ by randomly drawing their velocities from a Maxwell distribution every 100 steps with thermal velocity v_th = (k_BT/m)^1/2. We use the Lennard-Jones time t₀ = σ^WCA_core/v_th = [m(σ^WCA_core)²/ε^WCA_core]^1/2 as basic unit of time. Newton's equations of motion are integrated numerically with the velocity-Verlet algorithm⁵¹ using a numerical timestep of Δt = 7.2 × 10⁻⁴t₀.

We generated 100 statistically independent host structures for particle numbers N = 500, 1000, 2000, 4000, and 16 000 at fixed number density n = N/L² = 0.278(σ^WCA_core)⁻², corresponding to system sizes L/σ^WCA_core = 42.4, 60, 84.8, 120, and 240.

With each generated structure we associate a percolation threshold relying on a Voronoi tessellation of the particle positions of the host structure, in the same way as for the cherry-pit systems, see 2. Averaging over all 100 snapshots at the largest system size yields a critical effective interaction distance σ_c/σ^WCA_core = 0.990 ± 0.009 or equivalently, a critical reduced obstacle density n^*_c = nσ_c² = 0.272 ± 0.005.

It is instructive to structurally compare the WCA system to the cherry-pit model, employing an effective hard-core diameter. Here we use the Barker–Henderson diameter σ^BH_core, originally developed in the context of thermodynamic perturbation theory,^52–54

Numerical evaluation of the integral for σ_αβ = σ^WCA_core yields σ^BH_core ≈ 1.02σ^WCA_core, corresponding to an effective packing fraction of η = nπ(σ^BH_core)²/4 = 0.225. At this packing fraction, the cherry-pit model exhibits a similar percolation threshold, see Fig. 2.

The structure factor of the WCA system is included in Fig. 1, with the wave numbers measured in units of 1/σ^BH_core. It compares well to the one of the slightly denser cherry-pit system at η = 0.26. The positions of the first diffraction peak coincide, while the amplitude in the WCA system is slightly lower by ≈9%. Thus, the percolation threshold, the effective reduced density, and the structure factor of the WCA-system matrix can be mapped consistently onto the cherry-pit model.

The frozen matrices are explored by an ideal gas of tracers. The tracers interact with all matrix particles identically via the smoothly truncated WCA potential, eqn (5), with coefficients ε_WCA = 0.1ε^WCA_core and σ_αβ = σ_WCA. The interaction range σ_WCA is used as the control parameter and defines a reduced number density n^*_WCA = n(σ_WCA)². In the experiment by Skinner et al.,¹⁴ the tuning of the analogous tracer–matrix interaction is achieved by varying the external magnetic field. The tracer particles are inserted and equilibrated in the host structure by grand-canonical Monte-Carlo moves in combination with successive umbrella sampling.⁵⁵ Subsequently, the tracers are equilibrated using the simplified Andersen thermostat. Since the equilibration is performed in the canonical ensemble the average energy of each system is fluctuating. For the microcanonical production runs the systems are brought to the same average energy at the end of the equilibration period by rescaling all tracer velocities in the same system with the same constant, leaving the relative distribution of energies unchanged.

Newton's equations of motion are integrated numerically with the velocity-Verlet algorithm with the same timestep as for the host particles. Between 50 and 10 000 tracers for each host structure configuration are used to obtain ensemble averages for runs of up to nearly 10⁶t₀. For the calculation of time averages, typically 10 moving time origins per run were used and were spaced equidistantly over the whole simulation time.

The probability distribution of the energy per tracer p(E), as defined by p(E) = Z(β)⁻¹ exp[−βE]D(E) with the density of states D(E) and the partition function Z(β), can be directly calculated from the simulation data as the histogram of the tracer energy. For the binning of the energies, a bin width of ΔE/ε^WCA_core = 0.1 was chosen. The distribution p(E) has a peak at small energies, see Fig. 6, and decays exponentially at large energies, see inset. The energy distribution is nearly unchanged for all densities n^*_WCA. Merely slight variations in the peak height are observed, which are probably due to fluctuations in the potential energy frozen into the matrix.

Fig. 6 Energy distribution p(E) in the WCA-disk system for a canonical ensemble of tracers for a range of reduced densities n*. Inset: the same data in semilogarithmic presentation.

Tracer dynamics

The system undergoes a localization transition similarly to the overlapping Lorentz model and the cherry-pit model: At long times, the MSD becomes diffusive for n^*_WCA ≤ 0.4 and saturates for n^*_WCA > 0.4, see solid lines in Fig. 7. This implies a transition point (n^*_WCA)_c ≈ 0.4. At intermediate densities, the MSD is subdiffusive at intermediate times but it never matches the critical subdiffusion of the Lorentz model, with exponent 2/z with z = 3.036. This was already discussed shortly by some of us in ref. 14. This is even more apparent by direct inspection of the local exponent γ(t), see solid lines in Fig. 8. Instead of the Lorentz model exponent, the local exponent exhibits γ(t) ≈ 0.55 at n^*_WCA = 0.35 over almost three orders of magnitude in time.

Fig. 7 Mean-squared displacements of the WCA system for a canonical ensemble of tracers (solid lines) and tracers with exactly one energy (dashed lines) for a range of n^*_WCA. The straight line ∼ t^2/z with the dynamic exponent z = 3.036 of the Lorentz model serves as guide to the eye (data published previously in ref. 14).

Fig. 8 Local exponent of the mean-squared displacements of the WCA system for the same data as in Fig. 7, i.e. for a canonical ensemble of tracers (solid lines) and tracers with exactly one energy (dashed lines). The horizontal line indicates the anomalous exponent 2/z with z = 3.036 of the Lorentz model.

The situation changes qualitatively if all tracers are set to exactly the same energy, which restores the critical behavior.¹⁴ Then, the system undergoes a localization transition at n^*_WCA ≈ 0.320, where the MSD exhibits subdiffusion with the expected exponent, see dashed lines in Fig. 7 and 8.

The difference between these two systems is also strikingly apparent in the long-time diffusion coefficient, shown in Fig. 9. At large densities, the diffusion coefficient could not be directly measured from the MSD but was determined via finite-size scaling. For small separations ε from the critical point, the diffusion coefficient D is expected to vanish as D ∼ (−ε)^μ for ε → 0. If the size of the simulation box L is smaller than the correlation length ξ of the system, then this scaling is replaced by the finite-size scaling D ∼ L^−μ/ν for L ≪ ξ.¹⁹ For constant ε and incrementally increasing L, the diffusion coefficient will first follow the finite-size scaling D ∼ L^−μ/ν before converging to the true value at large-enough L ≫ ξ. This behavior is approximately fulfilled by the fitting function D = aL^−μ/ν + D_lower. Therefore, even if the simulated systems are not large enough to allow determining the true value of D, a fit to the small-L data will return a true lower bound D_lower, see Fig. 10 for an illustration. A true upper bound for D is given by the value obtained in the largest simulated system. With this procedure we calculated the bounds shown as vertical bars in Fig. 9.

Fig. 9 Diffusion coefficient D of the WCA system for a canonical ensemble of tracers and the single-energy case as function of the reduced obstacle density n^*_WCA. Connected symbols are obtained directly from the mean-squared displacements, isolated error bars at higher densities from finite-size scaling, see text. The solid black line ∝ (−ε)^μ with the conductivity exponent μ = 1.309 of the Lorentz model serves as guide to the eye. Inset: rectification plot of the same data.

Fig. 10 (a) Schematic representation of the finite size scaling of the diffusion coefficient D at some finite but small distance ε to the localization transition. Dots represent data obtained from simulations (not actually simulated here) plotted as a function of L^−μ/ν. The true scaling (red line) is some unknown function fulfilling L^−μ/ν at small L. Fitting to this small-L part provides a lower bound for D as L → ∞. (b) Illustration of a channel between two obstacles at distance 2σ_eff with the potential energy in greyscale. Obstacle centers are marked by dots, the equipotential line of the WCA potential at the energy E = 2V_αβ(σ_eff) where the channel closes is given in black, and the corresponding effective hard disks are given by red circles.

While the data of the canonical ensemble of tracers is not compatible with the critical behavior of the Lorentz model, D ∼ (−ε)^μ, the single-energy case is. The difference between the two cases is even starker in the rectification plot given in the inset of Fig. 9, where data following the critical power-law will fall on a straight line. While this holds for the single-energy case near the transition, where the critical asymptote becomes valid for roughly ε ≲ 0.1 as in the cherry-pit model, the canonical ensemble case shows a strong rounding.

Clearly, the canonical ensemble does not exhibit the critical dynamics of the Lorentz model, while the single-energy case does. This can be explained by an averaging of the dynamics in the case of the canonical ensemble. In contrast to the cherry-pit model, the WCA system contains finite energy barriers. As a consequence, the available void space and its topology are a function of tracer energy, i.e. barriers which can be surmounted by fast tracers may not be passable for slower tracers. It will be shown in the following how this notion can be quantified with the help of a mapping of the system onto hard disks.

Hard-disk mapping

The hard-disk mapping will yield an effective hard-disk interaction diameter σ_eff for each tracer as a function of n^*_WCA and its energy E and will thus show that the dynamics in the WCA-disk system can be understood as an average over a distribution of effective Lorentz models.

What is needed is a mapping of the WCA-disk system onto an equivalent system with (overlapping) hard-disk obstacles and a point-like tracer. In order for it to be useful, the mapping must conserve the topology of the void space: open channels have to stay open and closed channels have to remain closed under the mapping. Otherwise, the percolation transition of the void space would not be correctly mapped. While mappings such as the Barker–Henderson mapping, which was used to estimate the packing fraction of the matrix, or a mapping using the Andersen–Weeks–Chandler approximation⁵⁶ can be very successful for the mapping of glassy systems, for example, they do not guarantee conservation of topology. Greater care is necessary, here.

In two dimensions, a channel in the void space is defined by two obstacles. The potential landscape in such a channel has the shape of a saddle. A tracer is able to pass the channel if its energy matches or surpasses the potential energy on the saddle point of the channel, i.e. at the point exactly between the obstacles. In the presented mapping, the effective hard-core interaction distance σ_eff between a given tracer and the obstacles is then calculated as the closest distance between two obstacles forming a channel through which the tracer is barely able to pass.

For obstacles at a distance 2r, the potential energy in the center of the channel is given by 2V_αβ(r), eqn (5) (the smoothing function Ψ(r) can be neglected), and a tracer of energy E cannot pass the channel if E < 2V_αβ(r). Thus the topology of the accessible space is preserved if the soft obstacles are mapped to hard disks of radius σ_eff (assuming a point tracer) with the condition E = 2V_αβ(σ_eff), see Fig. 10b. Explicitly,

which has two positive solutions for σ_eff. Only one of them respects the cutoff of the potential σ_eff ≤ r_cut,

The reduced effective density of the matrix then reads

For the mapping to be successful, it has to correctly map the critical point as determined by the single-energy dynamics onto the percolation point of the matrix. From the dynamics, the critical point can be read off as n^*_WCA ≈ 0.320 where the simulation was performed at the tracer energy E/ε^WCA_core = 1.143. Via the hard-disk mapping this corresponds to an effective hard-disk critical radius of (σ_eff)_c/σ^WCA_core = 0.982 and a critical hard-disc reduced density (n^*_eff)_c = 0.268. This fully agrees with the percolation threshold determined above via Voronoi tessellation.

Energy-resolved dynamics

The mapping clearly exposes that tracers with different energies experience matrices with different densities n*. Thus, it is useful to consider the tracer dynamics as a function of tracer energy. To this end, the total energy of each tracer was calculated at the beginning of the simulation and tracers with similar energies were grouped into bins of width ΔE/ε^WCA_core = 0.1. The MSD was then calculated for each tracer and averaged over each energy bin. Since the particles of each bin have approximately the same energy E and same interaction range σ_WCA, their state can be uniquely expressed by the interaction diameter σ_eff = σ_eff (σ_WCA, E).

The energy distribution of the tracers p(E) corresponds to a distribution of effective densities p(n^*_eff), which can be directly calculated via

and p(n^*_eff) = 0, else. Note, that E(n^*_eff) is given by the inversion of eqn (9) and that dE/dn^*_eff is negative.

In Fig. 11b, the distributions p(n^*_eff) are shown for a range of n^*_WCA, with p(E) directly taken from the simulation, see Fig. 6. The diffusion coefficients calculated from the energy-resolved MSDs for the same systems are shown in Fig. 11a. To account for the trivial scaling with the microscopic time scale of the particles t₀ = σ_eff/ν with the velocity v of the particles, the diffusion coefficients have to be plotted rescaled as D/(σ_effν). The velocity v for each energy was extracted from the short-time behavior of the MSD, δr²(t; E) = ν²t². Without any further rescaling, the diffusion coefficient as a function of n^*_eff falls onto a single master curve, which is in agreement with the single-energy data. As the percolation transition at n^*_c = 0.268 is approached, the master curve approaches the critical behavior expected for the Lorentz model, D ∼ (−ε)^μ. This demonstrates clearly that the hard-disk mapping is successful and that the long-time dynamics of single WCA-particles are compatible with the Lorentz model dynamics.

Fig. 11 (a) Master plot of diffusion coefficients D in the WCA-disk system for a canonical ensemble of tracers resolved by their energy E. The white line shows D for the single energy case of Fig. 9. The red line indicates the critical asymptote from Fig. 9.(b) Distribution of effective reduced obstacle density in the same systems.

Furthermore, it is possible to demonstrate that not only the diffusion coefficient but also the full dynamics satisfies the critical scaling of the Lorentz model when mapped onto hard-disks. For the Lorentz model, the MSD is expected to follow an asymptotic scaling which incorporates the regimes of regular and anomalous diffusion, as well as the localized regime into a single functional form,

δr²(t) = t^2/zδR²_±(t/t_x), t₀ ≪ t, (11)

with t_x ∼ [small script l]^z ∼ |ε|^{−z(ν−β/2)} diverging at the transition.⁸ The scaling function δR²₋ holds on the delocalized side of the transition, and δR²₊ on the localized side. At long times, t_x ≪ t, regular diffusion is recovered by δR²₋(t/t_x) ∼ (t/t_x)^1−2/z, while on the localized side δR²₋(t/t_x) ∼ (t/t_x)^2/z holds. For t₀ ≪ t ≪ t_x, both scaling functions tend to the same constant, representing the regime of anomalous diffusion.

Both scaling functions are displayed in Fig. 12, where the energy-resolved MSD at n^*_WCA = 0.35 is divided by the critical asymptote t^2/z and time is rescaled appropriately. For this, the separation parameter ε was calculated from the effective interaction distance, ε = (n^*_eff − n^*_c)/n^*_c with n^*_c = 0.268. The collapse of the data onto the scaling functions is roughly as successful as in the overlapping Lorentz model¹⁷ and can in principle be further improved by considering corrections to scaling.³⁰

Fig. 12 Master plot of energy-resolved mean-squared displacement in the WCA-disk system for a canonical ensemble of tracers at n^*_WCA = 0.35. The MSD is divided by the critical asymptote ∼ t^2/z and shown as a function of rescaled time using the separation parameter ε =(n^*_eff − n^*_c)/n^*_c with n^*_c = 0.268. The energy E of the MSDs increases from bottom to top and the MSDs with smallest and largest energy are annotated.

Percolating fraction

To quantify the rounding of the transition, it is instructive to calculate the fraction of tracers with an energy sufficiently high to allow for long-range transport. At a given n^*_WCA, this fraction corresponds to the percolation probability, p_perc, of the effective hard-disk system. Provided that n^*_WCA is large enough that some tracers are on the localized side of the transition, p_perc is obtained as the integral over all subcritical states of p(n^*_eff),and p_perc = 1, otherwise. At large densities n^*_WCA, only tracers with the largest energies are delocalized. Then, it is reasonable to assume that p(E) ≈ Aβ exp(−βE), e.g. from inspection of the inset of Fig. 6, with A ≳ 1 (If the approximation were meant to hold for all E, then due to normalization A = 1 would hold exactly, but this would underestimate the probability distribution at large E). Furthermore, E(n^*_eff → 0) = +∞ holds. Thus, one finds

The approximation of the energy distribution by an exponential overestimates p_perc at small densities, but the approximation should become exact for large n^*_WCA. Therefore, p_perc > 0 holds for all finite n^*_WCA, but becomes exponentially suppressed at large densities.

4 Summary and conclusion

We have performed simulations in two dimensions of particle transport in two models of porous media which represent systematic steps away from the standard overlapping Lorentz model towards realistic systems. In the Lorentz model, a percolation transition in the void space entails a localization transition in the dynamics with anomalous transport being a key signature. Our systems allow testing which properties of porous media are necessary for anomalous transport and a localization transition.

In the cherry-pit model, the host matrix contains structural correlations which modify the structure of the void space. In the WCA model, interactions between tracer particles and the host particles are modeled with a purely repulsive and soft potential. This changes the topology of the void space by introducing a potential energy landscape with finite barriers. The dynamics have been analyzed in terms of the mean-squared displacement and quantities derived from it.

In the cherry-pit model, we have determined the percolation threshold as a function of the obstacle packing fraction, which is a measure of structural correlations contained in the host matrix. For both a weakly and a strongly correlated system, the localization transition is observed coinciding with the percolation threshold and the dynamics is found to be compatible with the critical predictions for the Lorentz model. However, the convergence to the universal predictions is poor, with a possible origin being corrections to scaling.

In the WCA model with a canonical ensemble of non-interacting tracer particles, a localization transition is observed, but the critical predictions do not apply, i.e. the transition is rounded. The behavior is similar to what has been observed in a quasi-two-dimensional experiment recently.¹⁴ The situation is clarified by a mapping of the WCA matrix onto hard-disks which reveals that the dynamics of each tracer can be fully mapped onto the Lorentz model as a function of its diameter and energy. This is confirmed by a scaling analysis of the dynamics as a function of tracer energy. The dynamics of the full system thus represent an energy-average over a distribution of effective Lorentz models. As a consequence, in systems with soft potentials like WCA-disks, one can only observe the idealized Lorentz model scenario in a simulation where the energy of tracers can be precisely controlled and held constant over the whole simulation, i.e. only for Newtonian dynamics. In Brownian dynamics or for interacting tracers, we expect the rounding of the transition to become more pronounced, as each tracer samples the full energy distribution over time.

Acknowledgements

This work has been supported by the Deutsche Forschungsgemeinschaft DFG via the Research Unit FOR1394 “Nonlinear Response to Probe Vitrification”.

References

1. H. Gleiter, Acta Mater., 2000, 48, 1–29.
2. H. Brenner and D. Edwards, Macrotransport Processes, Butterworth-Heinemann, 1993.
3. O. Bénichou, C. Chevalier, J. Klafter, B. Meyer and R. Voituriez, Nat. Chem., 2010, 2, 472–477.
4. D. Ben-Avraham and S. Havlin, Diffusion and Reactions in Fractals and Disordered Systems, Cambridge University Press, Cambridge, 1st edn, 2000.
5. A. Bunde, Solid State Ionics, 1998, 105, 1–13.
6. T. Voigtmann and J. Horbach, Europhys. Lett., 2006, 74, 459–465.
7. M. Weiss, in New Models of the Cell Nucleus: Crowding, Entropic Forces, Phase Separation, and Fractals, ed. R. Hancock and K. W. Jeon, Academic Press, 2014, vol. 307, ch. 11, pp. 383–417.
8. F. Höfling and T. Franosch, Rep. Prog. Phys., 2013, 76, 046602.
9. I. M. Sokolov, Soft Matter, 2012, 8, 9043.
10. M. J. Saxton, Biophys. J., 2012, 103, 2411–2422.
11. P. Adler, Porous Media: Geometry and Transports, Butterworth-Heinemann Limited, 1992.
12. H. Lorentz, Proc. R. Acad. Sci. Amsterdam, 1905, 7, 438–453.
13. H. V. Beijeren, Rev. Mod. Phys., 1982, 54, 195–234.
14. T. O. E. Skinner, S. K. Schnyder, D. G. A. L. Aarts, J. Horbach and R. P. A. Dullens, Phys. Rev. Lett., 2013, 111, 128301.
15. W. Götze, E. Leutheusser and S. Yip, Phys. Rev. A, 1981, 23, 2634.
16. W. Götze, E. Leutheusser and S. Yip, Phys. Rev. A, 1982, 25, 533–539.
17. F. Höfling, T. Franosch and E. Frey, Phys. Rev. Lett., 2006, 96, 165901.
18. F. Höfling and T. Franosch, Phys. Rev. Lett., 2007, 98, 140601.
19. F. Höfling, T. Munk, E. Frey and T. Franosch, J. Chem. Phys., 2008, 128, 164517.
20. T. Bauer, F. Höfling, T. Munk, E. Frey and T. Franosch, Eur. Phys. J.: Spec. Top., 2010, 189, 103–118.
21. J. Horbach, W. Kob and K. Binder, Phys. Rev. Lett., 2002, 88, 125502.
22. A. J. Moreno and J. Colmenero, J. Phys.: Condens. Matter, 2007, 19, 466112.
23. A. Moreno and J. Colmenero, Phys. Rev. Lett., 2008, 100, 126001.
24. A. J. Moreno and J. Colmenero, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 021409.
25. A. J. Moreno and J. Colmenero, J. Chem. Phys., 2006, 125, 164507.
26. T. Voigtmann and J. Horbach, Phys. Rev. Lett., 2009, 103, 205901.
27. T. Voigtmann, Europhys. Lett., 2011, 96, 36006.
28. J. Horbach, T. Voigtmann, T. Franosch and F. Höfling, Eur. Phys. J.: Spec. Top., 2010, 189, 141–145.
29. M. Spanner, S. K. Schnyder, F. Höfling, T. Voigtmann and T. Franosch, Soft Matter, 2013, 9, 1604.
30. A. Kammerer, F. Höfling and T. Franosch, Europhys. Lett., 2008, 84, 66002.
31. S. Torquato, Random Heterogeneous Materials, Springer-Verlag, 2002.
32. C. Dress and W. Krauth, J. Phys. A: Math. Theor., 1995, 597, L597–L601.
33. N. Metropolis, A. Rosenbluth, M. Rosenbluth and A. Teller, J. Chem. Phys., 1953, 21, 1087–1092.
34. E. P. Bernard and W. Krauth, Phys. Rev. Lett., 2011, 107, 155704.
35. S. C. Kapfer and W. Krauth, Soft-disk melting: from liquid-hexatic coexistence to continuous transitions, arXiv: 1406.7224v1 [cond-mat.stat-mech].
36. J. K. Percus and G. J. Yevick, Phys. Rev., 1958, 110, 1–13.
37. M. Adda-Bedia, E. Katzav and D. Vella, J. Chem. Phys., 2008, 128, 184508.
38. C. H. Rycroft, Chaos, 2009, 19, 041111.
39. A. R. Kerstein, J. Phys. A: Math. Gen., 1983, 16, 3071–3075.
40. J. A. Quintanilla and R. M. Ziff, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 051115.
41. D. Frenkel and B. Smit, Understanding Molecular Simulation. From Algorithms to Applications, Academic Press, London, 2nd edn, 2002.
42. P. H. Colberg and F. Höfling, Comput. Phys. Commun., 2011, 182, 1120–1129.
43. F. Höfling, in preparation.
44. B. I. Halperin, S. Feng and P. N. Sen, Phys. Rev. Lett., 1985, 54, 2391–2394.
45. J. Machta, R. A. Guyer and S. M. Moore, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 4818–4825.
46. P. Grassberger, Phys. A, 1999, 262, 251.
47. F. Höfling, K.-U. Bamberg and T. Franosch, Soft Matter, 2011, 7, 1358–1363.
48. Y. Gefen, A. Aharony and S. Alexander, Phys. Rev. Lett., 1983, 50, 77–80.
49. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237.
50. H. C. Andersen, J. Chem. Phys., 1980, 72, 2384.
51. K. Binder, J. Horbach, W. Kob, W. Paul and F. Varnik, J. Phys.: Condens. Matter, 2004, 16, 429.
52. J. A. Barker and D. Henderson, J. Chem. Phys., 1967, 47, 4714–4721.
53. D. Henderson, Mol. Phys., 1977, 34, 301–315.
54. J. Hansen and I. McDonald, Theory of simple liquids, Academic Press, London, 3rd edn, 2006.
55. P. Virnau and M. Müller, J. Chem. Phys., 2004, 120, 10925–10930.
56. M. Schmiedeberg, T. K. Haxton, S. R. Nagel and A. J. Liu, Europhys. Lett., 2011, 96, 36010.

---