Rounding of the localization transition in model porous media

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.


Introduction
][9][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. 8In 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. 14 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 tracermatrix 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.6][17][18][19][20] Concomitantly, transport becomes anomalous as manifested in a non-linear, power-law growth of the mean-squared displacement, δr 2 (t) ∼ t 2/z with a universal dynamic exponent z.In the case of soft interactions, the experimental and simulation results 14 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 21 or polymer blends with monomer-size disparity. 22,23Recent work [24][25][26][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, 28 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 tracermatrix interactions, many facets of the localization transition are well understood. 8,19,29Most importantly, the localization transition is due to an underlying continuum percolation transition of the accessible void space, 17 which is accompanied by a series of scaling laws familiar from the theory of critical phenomena of continuous phase transitions. 4Above 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,30Correlations 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. 31Second, 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.

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 2 )πσ 2 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 R R(t) which are piece-wise straight lines.The tracer particles are disks of finite diameter σ T , contrarily to the original overlapping Lorentz model (??a,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 2 )σ 2 , can be defined.The velocity of the tracer is of fixed magnitude v and defines a time scale t o = σ/v.Transport is controlled by variation of n * at fixed η.
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 32 .In the displacement moves, a random particle i with position r r r i is displaced to a new position r r r i + δ δ δ, where the vector δ δ δ is randomly chosen such that |δ δ δ| < σ core .This move is accepted according to a standard Metropolis criterion. 33luster 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 1 ,C 2 ), C 1 = C 2 , defined as two sets of disks satisfying the following condition: When all disks in C 1 are reflected at the pivot, each of them overlaps with at least one disk in C 2 , but none overlaps with disks not in C 2 , 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 1,000 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,35To check whether the system was sufficiently equilibrated we monitored the structure factor S (q) and the pair correlation function g(r) and compared it to the Percus-Yevick approximation. 36,37Additionally, 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 q q|, see ??.The angled brackets represent an ensemble average over the disorder, and {r r 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.
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σ 2 c /L 2 of the reduced number density by varying the distance σ of the tracer-obstacle interaction.First, we have computed the edges of a Voronoi tesselation of the matrix using the free voro++ software library. 38After removal of the edges with a distance smaller than σ to an obstacle center, the obtained network represents the volume accessible to the tracer particle. 39Upon 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) 19 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 ??.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 give an estimate of the width of the distribution of the percolation thresholds.This gives an estimate for the percolation density, 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.
The overlapping Lorentz model corresponds to η = 0, here the percolation threshold is known accurately: 40 n * c (0) = 0.359 081 0± 0.000 000 6 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 (η) = Mean-squared displacement δr 2 Mean-squared displacement δr 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 η hcp = (π/6) √ 3 ≈ 0.9, 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 ) 2 respectively.Note that at η = 0.26 the structure factor closely resembles that of the WCA system discussed later on (??d).

Tracer dynamics
The mean-squared displacements (MSD) 2 were obtained as time-and ensembleaverage over 20 obstacle configurations containing N = 516, 468 obstacles each in production runs up to times 10 9 t o .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,42he MSD are shown in ?? 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 2 (t) 4Dt for n * < n * c with diffusion coefficient D, or saturates, δr 2 (t) 2 for n * > n * c with localization length .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 n * > n * c , the spanning cluster disappears and all tracers are confined.This implies that the localization length 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, ( 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, 46 random walkers on percolation lattices 30 and in the overlapping 2D Lorentz model. 20,47The value was confirmed only recently also for the overlapping 2D Lorentz model with ballistic tracers. 43Our data for the MSD in the cherry-pit model suggest anomalous transport with effective exponents slightly lower than the universal one (??).
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. 30However, the dynamics is extremely sensitive to the density near the percolation threshold.Also, 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 0 in ??.Additionally, the overlapping 2D Lorentz model with ballistic dynamics exhibits strong, non-universal corrections to scaling, which modify the effective exponent over long periods of time. 43In particular, γ(t) slowly approaches 2/z from below.It is thus entirely expected that the cherry-pit exhibits a similarly slow convergence.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, 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 twodimensional standard percolation, ν = 4/3 and β = 5/36 hold exactly, 4 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 ??.Plotting D 1/µ vs. n * (see inset) rectifies the critical law, ??, 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. 43

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, 49 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 B T/ε WCA core = 1.0.To improve numerical stability the potential is multiplied with a smoothing function Ψ(r) := (r − r cut ) 4 /[h 4 + (r − r cut ) 4 ] with the width h = 0.005 σ WCA core .The particle configurations are equilibrated using a simplified Andersen thermostat 50 by randomly drawing their velocities from a Maxwell distribution every 100 steps with thermal velocity v := (k B T/m) 1/2 .We use the Lennard-Jones time as basic unit of time.Newton's equations of motion are integrated numerically with the velocity-Verlet algorithm 51 using a numerical timestep of ∆t = 7.2 • 10 −4 t o .
We generated 100 statistically independent host structures for particle numbers N = 500, 1 000, 2 000, 4 000, and 16 000 at fixed number density n := N/L 2 = 0.278 (σ WCA core ) −2 , 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 tesselation of the particle positions of the host structure, in the same way as for the cherry-pit systems, see ??.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σ 2 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][53][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 ) 2 /4 = 0.225.At this packing fraction, the cherrypit model exhibits a similar percolation threshold, see ??.
The structure factor of the WCA system is included in ??, 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, ??, 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 ) 2 .In the experiment by Skinner et al. 14 , 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. 55Subsequently, 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 micronanonical 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 6 t o .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(β) −1 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 ??, 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.exponent, the local exponent exhibits γ(t) ≈ 0.55 at n * WCA = 0.35 over almost three orders of magnitude in time.

Tracer dynamics
The situation changes qualitatively if all tracers are set to exactly the same energy, which restores the critical behavior. 14Then, the system undergoes a localization transition at n * WCA ≈ 0.320, where the MSD exhibits subdiffusion with the expected exponent, see dashed lines in ?? and ??.
The difference between these two systems is also strikingly apparent in the long-time diffusion coefficient, shown in ??.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 ξ. 19 For constant and incrementally increas-ing 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 ?? 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 ??.
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 ??, 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 56 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), ?? (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 ??b.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 tesselation.

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 ?? and that dE/dn * eff is negative.In ??b, the distributions p(n * eff ) are shown for a range of n * WCA , with p(E) directly taken from the simulation, see ??.The diffusion coefficients calculated from the energy-resolved MSDs for the same systems are shown in ??a.To account for the trivial scaling with the microscopic time scale of the particles t o = σ eff /v with the velocity v of the particles, the diffusion coefficients have to be plotted rescaled as D/(σ eff v).The velocity v for each energy was extracted from the short-time behavior of the MSD, δr 2 (t; E) = v 2 t 2 .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.
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, with t x ∼ z ∼ | | −z(ν−β/2) diverging at the transition. 8The scaling function δR 2 − holds on the delocalized side of the transition, and δR 2 + on the localized side.At long times, t x t, regular diffusion is recovered by δR 2 − (t/t x ) ∼ (t/t x ) 1−2/z , while on the localized side δR 2 − (t/t x ) ∼ (t/t x ) 2/z holds.For t o t t x , both scaling functions tend to the same constant, representing the regime of anomalous diffusion.
Both scaling functions are displayed in ??, where the energyresolved 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 17 and can in principle be further improved by considering corrections to scaling. 30rcolating 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 ), p(E) dE, (12)   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 ??, 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).
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.

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 transi-tion 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 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-twodimensional experiment recently. 14The 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.

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

Fig. 1 :
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 cherrypit model.For comparison, the matrix structure factor of the WCA-disk system with the effective diameter σ BH core is included.

Fig. 2 :
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 4 and ≈ 5 • 10 5 obstacles.The overlapping Lorentz model corresponds to η = 0.For comparison: WCA system at the effective packing fraction η = 0.225 and n * c ≈ 0.272.

Fig. 4 :Fig. 5 :
Fig. 4: Local exponents of the mean-squared displacements of ?? 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 σ.

Fig. 6 :
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.

Fig. 8 :
Fig.8: Local exponent of the mean-squared displacements of the WCA system for the same data as in ??, 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.

Fig. 9 :Fig. 10 :
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 meansquared displacements, isolated errorbars at higher densities from finitesize 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. 11 :
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 ??.The red line indicates the critical asymptote from ??. b) Distribution of effective reduced obstacle density in the same systems.
Furthermore, E(n * eff → 0) = +∞ holds.Thus, one findsp perc ≈ A exp −βE(n * c ) = A exp −8βε WCA n * 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 ) 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 ??.Instead of the Lorentz model 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.(Datapublished previously in Ref.14) 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.