Modifying continuous-time random walks to model finite-size particle diffusion in granular porous media

The continuous-time random walk (CTRW) model is useful for alleviating the computational burden of simulating diffusion in actual media. In principle, isotropic CTRW only requires knowledge of the step-size, $P_l$, and waiting-time, $P_t$, distributions of the random walk in the medium and it then generates presumably equivalent walks in free space, which are much faster. Here we test the usefulness of CTRW to modelling diffusion of finite-size particles in porous medium generated by loose granular packs. This is done by first simulating the diffusion process in a model porous medium of mean coordination number, which corresponds to marginal rigidity (the loosest possible structure), computing the resulting distributions $P_l$ and $P_t$ as functions of the particle size, and then using these as input for a free space CTRW. The CTRW walks are then compared to the ones simulated in the actual media. In particular, we study the normal-to-anomalous transition of the diffusion as a function of increasing particle size. We find that, given the same $P_l$ and $P_t$ for the simulation and the CTRW, the latter predicts incorrectly the size at which the transition occurs. We show that the discrepancy is related to the dependence of the effective connectivity of the porous media on the diffusing particle size, which is not captured simply by these distributions. We propose a correcting modification to the CTRW model -- adding anisotropy -- and show that it yields good agreement with the simulated diffusion process. We also present a method to obtain $P_l$ and $P_t$ directly from the porous sample, without having to simulate an actual diffusion process. This extends the use of CTRW, with all its advantages, to modelling diffusion processes of finite-size particles in such confined geometries.


INTRODUCTION
Diffusion plays a key role in a wide range of natural and technological processes. A textbook modelling of such processes is the consideration of the diffusion of a single memory-free particle in a given medium. The nature of such a random walk is determined by three probability densities (PDFs): of the step size l, P l (l); of the step directionn, P n (θ, φ); and of the waiting time between steps t, P t (t). These PDFs are, in principle, position dependent, but it is standard practice to derive (or postulate) them assuming position-independence and that P n (θ, φ) is uniform. The diffusion is then modelled as a continuous time random walk (CTRW) in free space, mainly because this alleviates finite size errors due to finite samples. Specifically, the CTRW is constructed by adding vectors of uniformly random orientations, whose lengths are chosen from P l , with waiting times chosen from P t . Averaging over sufficiently many such independent processes, the dependence of the mean square distance (MSD) on time is determined through MSD = DT α . In normal diffusion α = 1 and D is the standard diffusion coefficient. When P t has a slowly decaying algebraic tail and P l does not, the random walk is sub-diffusive (α < 1) [1]. Alternatively, if P l has a slowly decaying algebraic tail and P t does not, the random walk is super-diffusive (α > 1), resembling a Lévy flight [2]. Diffusion processes that have the same value of α are said to be in the same universality class [3].
Using CTRW to model diffusion in confined geometries, such as porous media or different biological systems, is attractive [4][5][6] because it alleviates the need to simulate directly the dynamics of particles within the pore space, reducing significantly the computational burden. This practice is based on the common assumption that the three aforementioned distributions alone control the random walk's universality class. The common procedure is to find first the forms of these distributions in a specific medium, using either small simulations or analytic derivation under some assumptions, and then use these to carry out a many-step CTRW in free space. It is then presumed that the CTRW yields the same universality class as the diffusion in the confined geometry.
The aim of this paper is to demonstrate that this is a misconception when the diffusing particles have finite sizes. We do so by comparing a simulation in an actual porous sample with a CTRW in free space, based on statistical information of the porous media. We show that the medium's connectivity (topology) depends crucially on the particle's size and this affects directly the nature and universality class of the diffusion process. We then propose a method to correct for this effect, which makes it possible to still use CTRW, with its advantages, to model diffusion in confined geometries.
The structure of this paper is the following. In section 2 we describe the simulated porous samples. In section 3 we give details about the diffusion process and present results for a point-like particle. In section 4 we discuss the effects of particle size on the diffusion process and present explicit results. In section 5 we describe the equivalent CTRW simulations and show that these yield a different universality class in spite of having the same step length and waiting time distributions. We propose an explanation of this discrepancy. In section 6 we propose an improved CTRW modelling to alleviate this problem and provide better description of diffusion of particles in confined geometries. We conclude in section 7 with a discussion of the results.

THE POROUS SAMPLE
The three-dimensional porous sample is an opencell structure, generated using Surface Evolver [7] as follows. Initially, N seed points are distributed randomly and uniformly within a cube, and the cube's space is Voronoï-tessellated to determine the cell associated with each point. A cell around a point consists of the volume of all spatial coordinates closest to it. The resulting cellular structure is then evolved with Surface Evolver to minimise the total surface area of the cell surfaces. This procedure is used commonly to model dry foams and cellular materials whose dynamics are dominated by surface tension. The result is an equilibrated foam-like structure, comprising cells, membranes, edges and vertices. A membrane is a surface shared by two neighbouring cells, an edge is the line where three membranes coincide, and a vertex is the point where four edges coincide.
Next, we construct a tetrahedron around every vertex by connecting the mid-points of the four edges emanating from it [8]. Neighbouring tetrahedra are in contact in the sense that they share the midpoint of an edge. This construction results in a pseudo-granular structure, in which every tetrahedron represents a pseudo-grain in contact with exactly four others [9] (see figure 1). Since neighbouring pseudo-grains share the mid-point of the edge between them, the tetrahedra structure is topologically homeomorphic to the original structure. The void space surrounded by the pseudo-grains is still cellular, but a cell surface now consists of triangular facets -the faces of the pseudo-grains surrounding it -and throats -the skewed polygons remaining of the original cell membranes. The membranes over the throats are disregarded, resulting in an opencell porous structure, in which the throats are the openings between neighbouring cells.

THE DIFFUSION PROCESS AS A RANDOM WALK
Modelling of a diffusing particle in the porous medium as a random walk requires knowledge of a number of parameters and distributions: the particle size, the cell volume distribution, the throat area distribution and correlations between these. We start by considering particles that are considerably smaller than the smallest throat in the structure. The particle cannot enter the tetrahedral pseudograins, but only move from cell c to cell c through their shared throat. A step, l c,c , is a vector extending between the centres of these cells. We define a waiting time, t c , which is the number of time steps spent in cell c before a jump occurs. Inside a cell, the particle is assumed to undergo Brownian motion and t c is proportional to: (i) the square of the effective cell radius, R c ≡ 3vc 4π 1/3 , where v c is the cell volume; (ii) the inverse of the fraction of the cell's open surface through which the particle can pass to neighbouring cells, S c ≡ . This is since the particle, on average, has to make 1/S c journeys of length R c until it goes through a throat rather than hits a facet of a pseudo-grain. Here A (t/f) c is the total area of throats/facets that make the surface of the cell. Thus, the cell waiting time is where d is the local diffusion coefficient.
The probability to exit cell c into c is proportional to the area of the throat between them. To reduce finite size effects, we wrap the sample around itself with periodic boundary conditions and let the particle travel larger distances by re-entering the sample. This means that the same cell may occur at different locations along the random walk. To avoid distorting the statistics by using the same cell too many times, we stop the process once a cell has occurred at more than five different locations.
For each walk we calculate the particle's position at time T , relative to the origin, R(T ) ≡ N (T ) l n , where N (T ) is the number of steps before time T . We also calculate MSD(T ) ≡ R(T ) 2 , where the angular brackets denote average over 1000 walks. Figure 2 shows the time dependence of the MSD for one such simulation. As expected, we find that the MSD of the point-like particle is linear with time, with a diffusion coefficient of D = (0.65 ± 0.03)d.

DIFFUSION OF FINITE SIZE PARTICLES
Real diffusing particles have finite sizes which, as we shall see below, can change significantly the nature of the diffusion process. To study this effect, we next model the diffuser as a sphere of radius r. To compare different systems, it is useful to measure r in units of the average effective throat radius, r 0 . The finite size gives rise to two main ef-fects. Firstly, the larger r is the lower the probability of passing through any specific throat, since the effective area that the particle can pass through is smaller. Thus, the overall probability to escape from a cell decreases, increasing the waiting times spent inside cells. Secondly, as the particle size increases, the probability of passing through some throats vanishes identically, changing the system's connectivity for this particle. This size-driven topological modification of the network is especially significant for large particles, as they could be trapped inside disconnected regions. Below we study these two effects.

The waiting time distribution
The observed PDF P t (t), within the porous medium, changes dramatically with increasing particle size. While P t (t) of a point-like particle is nicely bounded, it develops a long tail for sufficiently large particles, P t (t > t 0 ) ∼ t −β , with β a function of r. Figure 3a shows a log-log plot of P t (t) for a particle of radius r = 1.3r 0 with β = 2.29 ± 0.06. In figure  3b we plot the corresponding time dependence of the MSD. We find that the process is sub-diffusive with α = 0.69 ± 0.004.
Sub-diffusion is often seen in CTRW, when the second moment of the step length distribution is finite and P t has an algebraic tail with β < 2, leading to a power-law dependence of the MSD on time with α = β − 1 [10]. It would be then natural to presume that the source of our observations is the emergence of a power-law tail in P t as r increases. Yet, figure 3 shows sub-diffusion even for β = 2.29 ± 0.06, indicating another cause for the sub-diffusion.
To study directly the effect of the confined geometry on the relation between α and β, i.e. between the universality class and the waiting time distribution, we varied r in the simulations and plotted α(r) vs. β(r) (blue dots in figure 5). For small particles, β is much larger than 2 and α = 1, corresponding to standard diffusion with narrow waiting-time PDFs. As r increases, β decreases and the walks eventually become sub-diffusive with α < 1. In our simulations, the transition occurs at β = 2.53 ± 0.03, in contrast to the theoretical prediction of CTRW, β = 2.

MODELLING THE DIFFUSION AS ISOTROPIC CTRW IN FREE SPACE
Our measurements of the universality class transition clearly disagree with the theoretical CTRW prediction. To pinpoint the reason for the disagreement, we examine a traditional isotropic CTRW model of this process, using the particle size-dependent PDFs, P l and P t , that the particle experiences while diffusing in the confined structure. One way to obtain these PDFs is to measure them empirically during a diffusion process. But more efficient is to compute them directly from the structural statistics of the porous medium. We use these as input to an unrestricted CTRW in free space.
A derivation of P l and P t from structural statistics should be made cautiously because a straightforward histogram of the waiting times ignores the inherent correlation between the probability to visit a cell, P (c), and the waiting time [11]. Cells that are difficult to get out of (long waiting times) are also difficult to get into and are therefore visited less frequently. Some cells are completely inaccessible for particles above a certain size.
To take this into account we need to derive P (c). To this end, we first simulate a diffusion process in the finite sample and plot P (c) against the cell total throat area, A (t) c ( figure 4). We find that, to a good accuracy and across 3.5 orders of magnitude, P (c) is proportional to A (t) c , i.e. to its total open surface area, regardless of the cell volume. Using this observation, we can estimate P (c) for a particular particle size by using the effective A (t) c -a direct structural characteristic of the medium. In principle, one expects the visiting probability to be correlated with the visiting probabilities of neighbouring cells, but figure 4 shows that this correlation is negligible. Having P (c), we can now estimate more accurately P t for every particle size. This is done by manipulating the waiting-time histogram in the following way. For every bin consisting of the waiting times of cells {c 1 , ..., c n }, we multiply the bin's height by n 1 P (c i ). Then we normalise the histogram to still represent a PDF. Generally speaking, this suppresses cells with longer waiting times and makes P t less heavy-tailed. For heavy-tailed PDFs, we obtain β using a log-log plot, as we did for a diffusion process in figure 3a.
Collecting the statistics of all possible steps in the sample, we get a PDF described well by P l (l) ∼ exp{−(l − l 0 ) 2 /2σ 2 }, with σ/l 0 = 0.15 ± 0.02. We note that it is also possible to derive P l analytically, given the cell volume distribution and nearest neighbour volume-volume correlations. This, however, is somewhat downstream from the main thrust of this paper and will be detailed in a later report.
We are now able to obtain accurate step-length and waiting-time PDFs for every particle size. Figure 5 shows results of simulations (green dots) of CTRW processes in free space using these PDFs. The curve following them still differs from the theoretical relation α = β − 1, which we believe is purely a finite-time effect. Our measured transition for CTRW is at β = 2.01 ± 0.02, in very good agreement with the theoretical prediction.
Significantly, the normal-to-subdiffusion transition takes place at a lower value of β in the CTRW than in the actual sample. Since the step-length and waiting-time distributions are the same for both processes, we conclude that the reason for the discrepancy is topological, as we now proceed to discuss.
As mentioned above, the size of the diffusing particle determines the effective throat sizes. It follows that the topology is particle size dependent. Moreover, above some size there is no percolating cluster between the the sample's boundaries. Figure 6 shows the dependence of the percolating accessible space on the particle size, where a range of particle sizes around the percolation threshold is marked by red lines. The same range is marked in figure 5. We see that the transition for diffusion in the sample occurs within this range. This supports the idea that the topology plays an important role in determining the universality class. Around the percolation threshold the incipient cluster assumes a fractal-like structure, further biasing the process towards subdiffusion [12,13].

RANDOM WALK WITH MEMORY
The above results suggest that, to usefully model diffusion in porous media, the CTRW needs to be modified. The failure of conventional CTRW is in modelling the connectivity of the structure. As dis- Figure 5. Two sets of simulations -diffusion in a porous sample (blue dots) and CTRW in unconfined space (green dots). Each simulation (dot) represents a particular particle size, and is described by the powerlaw of the waiting-time distribution, β, and the anomaly parameter, α. The blue and green lines are fits of the form α = 1 − 1 2 exp{−(β − β0)/τ }, with (β0 = 2.03, τ = 0.44) and (β0 = 1.51, τ = 0.40), respectively. The theoretical CTRW prediction of a universality class transition at β = 2 is denoted by the black dashed line. The red lines mark the range of β's that corresponds to the percolation threshold of the sample. This range matches the universality-class transition for confined diffusion.
cussed above, the process is not entirely Markovian and there are correlations between steps. Moreover, there are correlations between the waiting time and the step direction distributions. Figure 7 shows the PDF for the angle between successive steps of a diffusion process in the porous sample. In particular, it shows a finite probability to make a backward step with P back ≡ P θ (θ = π) = 0.087 ± 0.001. It also shows that P back is definitely larger than 1/13.7 ≈ 0.073, which is the inverse of the average number of throats per cell. This is since the throat, through which a particle enters a cell, is generically larger than average. The enhanced backward step probability can be seen as a correlation between successively used throats. As r increases the total available throat area decreases and this correlation (and P back ) increases.
We examine two methods to modify the CTRW, both introducing a backward bias. In the first, the particle 'remembers' the previous vector step and the next one is chosen relative to it. The relative angle is chosen according to a non-uniform distribution  P θ (θ) which, given the particle radius r, can be calculated directly from the porous structure, similarly to P l and P t . To uniquely define the direction of the new vector step we randomly choose an azimuthal angle, using the previous step as z axis. The addition of anisotropy and a finite P back corrects for the changing connectivity in the porous structure.
Running a set of simulations with the addition of P θ brings the universality-class transition at β = 2.1 ± 0.03, closer to the one measured in the sample, improving somewhat the CTRW model.
In the second method, we take into account both the direction of the last step and the throat area. This is done by simulating an entire cell at every step. We choose the cell volume, the number of throats, and their areas, from the relevant distributions of the porous sample. We then calculate the waiting time for the chosen cell, and choose randomly an exit throat. As in the porous sample, the probability to go through a throat is proportional to its effective area. The step length is the sum of two effective cell radii -the one exited and the one entered -and its direction random. The throat through which the particle exits the cell is always one of the throats for the next cell. If this throat is chosen again, the last step is retraced. In the following we refer to this method as DA, denoting the 2 memory parameters -step Direction and throat Area.
Using DA, once a large throat has been chosen randomly, the particle may pass through it back and forth several times before it moves on and loses memory about this throat. Thus, while the first method includes correlation between successive stepdirections, DA enhances it. In other words, not only does it simulate the correlation between successive steps, it also simulates the correlation between successive backward steps.
With the DA method, we obtain a universalityclass transition at β = 2.50±0.03, in excellent agreement with the original simulation results. We conclude that this model captures much better the particle size-driven topological change.
An important feature of DA is that P back is higher than in the porous sample. To understand this, consider two cells, connected by a relatively large throat, and connected to the rest of the structure only by one other narrow throat. Once the particle enters this structure, it is likely to move back and forth many times before it emerges. But the smaller the throat leading to this sub-structure, the less probable this situation is in the first place. Namely, the more a sub-structure potentially increases P back , the more suppressed it is. In comparison, consider a cell in DA with a relatively large throat. Recall that, other than the entrance throat, a new cell realisation is created at each step. Hence, here the particle is also likely to oscillate between the cells, but the probability of this to occur is no longer small. It is this that allows the DA method to compensate for other missing topological information.

CONCLUSIONS
We compared between two numerical models of diffusion: a direct simulation of the diffusion process in a computer-generated porous sample, and the other an unconfined isotropic CTRW. We presented a method to construct genuine step-length and waiting-time PDFs, P l and P t , based on statistical information of the porous media. Thus, we confirmed that, for every particle size, the two models use the same P l and P t . For each model we presented the dependence of the universality class on the waiting-time distribution. We demonstrated that the two models show different behaviours: while the CTRW simulations follow the theoretical prediction (up to finite-time effects) with transition to subdiffusion at β = 2, the diffusion simulations in the confined geometry exhibit a transition at β = 2.5. We showed that the difference stems from a particle size-driven change of accessible topology in the actual sample, which is not taken into consideration in the CTRW model. This has been supported by showing that the universality class transition and the percolation transition occur at the same range of particle sizes.
We emphasise that, in contrast to naive expectation, the discrepancy in the location of the transition is not due to different waiting-time distribution -the same distribution can lead to different universality classes, as we demonstrated.
It is interesting to note that a similar process arises when probe particles diffuse in biological networks of entangled filaments, e.g. F-actin [6]. A sub-diffusive motion was found when the particles size was comparable to the typical network mesh size due to long time trapping inside 'pores'. Indeed, it was found that the in-pore waiting times were distributed with a power-law tail, which was offered as the explanation for this effect. The experimental study in [6] concluded that the universality class of the diffusion process depends only on the ratio of the probe radius to mesh size, in agreement with our results. However, a scrutiny of their data shows that there is no size-driven topological change in their experiments and their values of α and β were in good agreement with the CTRW model. This could be due to the partial mechanical flexibility of F-actin filaments, which allows trapped particles to eventu-ally escape by deformation of the pores. In contrast, our rigid porous structures do not allow such events, i.e. throats cannot expand momentarily to let a large particle out.
As mentioned, the main advantage of the CTRW is that it overcomes potential finite-size problems. However, this is achieved at the expense of ignoring topological information about local connectivity. To improve the modelling of diffusion in porous media by CTRW, we presented two anisotropic models. One model includes memory of the last step direction and a non-uniform distribution of the step direction. The second model, DA, remembers, in addition, the area of the last throat, adding correlation between backward steps. DA shows a universalityclass transition at β = 2.50 ± 0.03, in good agreement with the one measured in confined geometry. We conclude that DA is a better alternative to the traditional CTRW for modelling diffusion in porous media. It combines the CTRW advantages -overcoming the finiteness of the sample and convenience of application -with a better capturing of the universality class.