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

0.63 MB·English

by  Shahar Amitai

Noname manuscript No.
(will be inserted by the editor)
Modifying continuous-time random walks to model
finite-size particle diffusion in granular porous media
Shahar Amitai · Raphael Blumenfeld
6
1
0
2

c
e Received:date/Accepted:date
D

5
Abstract Thecontinuous-timerandomwalk(CTRW) good agreement with the simulated diffusion process.

model is useful for alleviating the computational bur- We also present a method to obtain P and P directly
] l t
h den of simulating diffusion in actual media. In prin- fromtheporoussample,withouthavingtosimulatean
c
ciple, isotropic CTRW only requires knowledge of the actualdiffusionprocess.ThisextendstheuseofCTRW,
e
m step-size, P , and waiting-time, P , distributions of the withallitsadvantages,tomodellingdiffusionprocesses
l t
random walk in the medium and it then generates pre- of finite-size particles in such confined geometries.
-
at sumablyequivalentwalksinfreespace,whicharemuch
Keywords Anomalous diffusion · Granular pore
t faster. Here we test the usefulness of CTRW to mod-
s space · Continuous-time random walk · Anisotropy
. ellingdiffusionoffinite-sizeparticlesinporousmedium
t
a generatedbyloosegranularpacks.Thisisdonebyfirst
m
simulatingthediffusionprocessinamodelporousmedium
1 Introduction
- of mean coordination number, which corresponds to
d
n marginal rigidity (the loosest possible structure), com- Diffusionplaysakeyroleinawiderangeofnaturaland
o putingtheresultingdistributionsPlandPtasfunctions technological processes. A textbook modelling of such
c of the particle size, and then using these as input for a processesistheconsiderationofthediffusionofasingle
[
freespaceCTRW.TheCTRWwalksarethencompared memory-free particlein a givenmedium. The natureof

5 to the ones simulated in the actual media. such a random walk is governed by three probability
v
Inparticular,westudythenormal-to-anomaloustran- densityfunctions(PDFs):ofthestepsize,P (l );ofthe
8 l i
9 sition of the diffusion as a function of increasing par- stepdirection,P (nˆ );andofthewaitingtimebetween
n i
9 ticle size. We find that, given the same Pl and Pt for steps,Pt(ti).ThesePDFsare,inprinciple,positionde-
3 the simulation and the CTRW, the latter predicts in- pendent, but it is standard practice to derive (or pos-
0
correctly the size at which the transition occurs. We tulate) them assuming position-independence and that
.
1 show that the discrepancy is related to the dependence P (nˆ ) is uniform. The diffusion is then modelled as a
0 n i
oftheeffectiveconnectivityoftheporousmediaonthe continuous-time random walk (CTRW) in free space.
5
1 diffusing particle size, which is not captured simply by Specifically, the CTRW is constructed by adding vec-
: these distributions. tors of uniformly random orientations, whose lengths
v
i WeproposeacorrectingmodificationtotheCTRW are chosen from Pl, at time intervals chosen from Pt.
X
model – adding anisotropy – and show that it yields Averagingoversufficientlymanysuchindependentpro-
r
cesses, the dependence of the mean square distance
a
ShaharAmitai (MSD)ontimesatisfies(cid:104)x2(cid:105)=Dtα.Innormaldiffusion
ESE,ImperialCollegeLondon,LondonSW72AZ, UK
α = 1 and D is the standard diffusion coefficient. But
E-mail:[email protected]
when P and/or P are very wide, the diffusion might
l t
RaphaelBlumenfeld
become anomalous (α(cid:54)=1). In particular, when P has
ESE,ImperialCollegeLondon,LondonSW72AZ, UK t
a slowly decaying algebraic tail and P does not, the
CollegeofScience, NUDT,Changsha,Hunan,China l
CavendishLaboratory,JJThomsonAvenue,CambridgeCB3 random walk is sub-diffusive (α < 1) [1, 2]. Alterna-
0HE,UK tively, if P has a slowly decaying algebraic tail and P
l t
2 ShaharAmitai,RaphaelBlumenfeld
does not, the random walk is super-diffusive (α > 1), that the CTRW yields the same universality class as
resembling a L´evy flight [3]. Diffusion processes that the diffusion in the confined geometry.
have the same value of α are said to be in the same The first aim of this paper is to demonstrate that
universality class [4]. thisdoesnotapplywhenthesizeofthediffusingparti-
cleiscomparabletothroatsizes.Wedosobyanalysing
Anomalousdiffusioncanarisefromdifferentsources,
trajectories of individual particles diffusing in a porous
whichcanonlybeidentifiedbygoingbeyondtheMSD.
sampleandshowstatisticaldeviationsfromCTRWpre-
Whensingleparticletrackingispossible,themovement
dictions. We also compare these simulations with an
canbeevaluatedbythetime-averagedMSD(TAMSD),
equivalent CTRW model. We show that the effective
δ2(t,T) [5].While the MSD is the ensemble average of
change in the medium’s connectivity with varying par-
the squared distance, made during a time interval t,
ticle size affects directly the nature and universality
over different realisations, the TAMSD, δ2(t,T), is the
class of the diffusion process. We conclude that the
average of the same quantity along a single trajectory
sub-diffusion is the result of CTRW on a percolation
of length T. Within the model of sub-diffusive CTRW,
cluster. Indeed, a combination of underlying mecha-
the TAMSD satisfies (cid:104)δ2(cid:105) ∼ t·Tα−1, where the angu-
nisms, leading to sub-diffusion, has also been observed
lar brackets denote a further ensemble average. In con-
in [12, 13, 14, 15].
trast, the MSD is sub-linear in t, which makes CTRW
Thesecondaimofthepaperistoproposeamethod
non-ergodic – the time-average and ensemble-average
tocorrectforthetopologicaleffect,whichmakesitpos-
differ. In particular, the dependence of the TAMSD on
sible to still use CTRW, with its advantages, to model
T points to the ageing nature of CTRW [5].
diffusion of any finite size particle in confined geome-
Akeyfeatureofsub-diffusiveCTRWistherandom-
tries. To maximise the range of validity of our results
ness of its TAMSD. Since P is scale free, the longest
t (see discussion below), we consider very high porosity
waitingtimeseachindividualtrajectoryencountersvary
porousmedia.Thesecorrespondtomarginallyrigidas-
significantly, as do the amplitudes of the individual
semblies of frictional particles, whose mean coordina-
TAMSDs. To quantify this, we define the amplitude
tion number is four [16]. The least confined of these
scatter, ξ =δ2/(cid:104)δ2(cid:105). For ergodic processes (e.g. α=1)
are model systems whose each particle has exactly four
its PDF is P(ξ) = δ(ξ − 1) for sufficiently long tra-
contacts.
jectory times. But within CTRW this PDF broadens
The structure of this paper is the following. In sec-
as α decreases. Defining the ergodicity breaking (EB)
tion 2 we describe the simulated porous samples. In
parameter, EB=(cid:104)ξ2(cid:105)−(cid:104)ξ(cid:105)2, it can be derived analyt-
section 3 we describe the diffusion process, and dis-
ically for CTRW processes as a monotonically decreas-
cuss the effects of particle size. We perform statistical
ing function of α.
analysis of the particle trajectories and show disagree-
Anothercauseforsub-diffusioniswalkinginafractal- ments with some predictions of the CTRW model. In
like environment [6, 7]. Such environment is charac- section4wedescribetheequivalentCTRWsimulations
terised by a network of narrow passages and dead ends andshowthattheyyielddifferentbehaviourinspiteof
atdifferentlengthscales,whichhinderthewalk.Unlike havingthesamestep-lengthandwaiting-timedistribu-
CTRW,thisprocessisstationaryandthereforeergodic. tions. We propose an explanation for this discrepancy.
TheTAMSD,liketheMSD,issub-linearint,indepen- In section 5 we propose a modification to the conven-
dent of T and its EB parameter vanishes. tional CTRW model to alleviate this problem, making
itmoresuitableformodellingdiffusionoffinitesizepar-
UsingCTRWtomodeldiffusioninconfinedgeome-
ticles in confined geometries. We conclude in section 6
tries, such as porous media formed by either sintered
with a discussion of the results.
orunconsolidatedgranular materials,isvery attractive
[8, 9, 10, 11] because it alleviates the need to simu-
late directly the dynamics of particles within the pore 2 The porous sample
space,reducingsignificantlythecomputationalburden.
In addition, it alleviates finite-size errors due to finite To simulate a three-dimensional porous granular as-
samples.Thispracticeisbasedonthecommonassump- sembly of coordination number four, we first generate
tionthatP ,P andP alonecontroltherandomwalk’s an open-cell structure, using Surface Evolver as follows
l t n
universalityclass.Thecommonprocedureistofindfirst [17, 18]. Initially, N seed points are distributed ran-
the forms of these distributions in a specific medium, domly and uniformly within a cube, and the cube’s
usingeithersmallsimulationsoranalyticderivationun- space is Vorono¨ı-tessellated to determine the cell as-
der some assumptions, and then use these to carry out sociatedwitheachpoint.Acellaroundapointconsists
a many-step CTRW in free space. It is then presumed of the volume of all spatial coordinates closest to it.
Modifyingcontinuous-timerandomwalks... 3
The resulting cellular structure is then evolved with 3 Diffusion in the porous sample
Surface Evolver to minimise the total surface area of
the cell surfaces. This procedure is used commonly to We model the diffuser as a sphere of radius r, mea-
model dry foams and cellular materials whose dynam- sured in units of the average effective throat radius, r .
0
ics are dominated by surface tension. The result is an We start by considering particles that are considerably
equilibratedfoam-likestructure,comprisingcells,mem- smaller than the smallest throat in the structure. The
branes, edges and vertices. A membrane is a surface particlecannotenterthetetrahedralpseudo-grains,but
shared by two neighbouring cells, an edge is the line only move from cell c to cell c(cid:48) through their shared
where three membranes coincide, and a vertex is the throat. The simulation progresses by moving the par-
point where four edges coincide. ticle from one cell, c, to a neighbouring cell, c(cid:48). Each
suchaneventisastep,l ,namelyavectorextending
c,c(cid:48)
between the centres of these cells. We define a wait-
ing time, t , which is the number of time steps spent
c
in cell c before a jump occurs. Inside a cell, the parti-
cle is assumed to undergo Brownian motion and t is
c
proportional to (i) the square of the effective cell ra-
dius, Rc ≡ (cid:0)34vπc(cid:1)1/3, where vc is the cell volume, and
(ii) the inverse of the fraction of the cell’s open surface
through which the particle can pass to neighbouring
cells, Sc ≡ A(ctA)+(ctA)(cf). This is since the particle, on av-
erage, has to make 1/S journeys of length R until
c c
it goes through a throat rather than hits a facet of a
pseudo-grain. The effective area of a throat is the area
through which a particle of radius r can pass and A(t)
c
(A(f))isthetotalareaofthroats(facets)thatmakethe
c
surface of the cell. A(t) is then the sum of the effective
c
Fig.1:Pseudo-grainsaroundthe foam vertices. throat areas accessible for the particle to go through.
Thus, the waiting time within a cell is
Next, we construct a tetrahedron around every ver- R2 1 (cid:18)3v (cid:19)2/3(cid:32) A(f)(cid:33)
tex by connecting the mid-points of the four edges em- t ≡ c = c 1+ c , (1)
anating from it [19]. Neighbouring tetrahedra are in c 2dSc 2d 4π A(ct)
contact in the sense that they share the mid-point of
anedge.Thisconstructionresultsinapseudo-granular where d is the local diffusion coefficient, which, using
structure of volume fraction φ = 34%, in which ev- the Stokes-Einstein relation, is inversely proportional
ery tetrahedron represents a pseudo-grain in contact to the particle radius, d=(r0/r)d0.
with exactly four others [20] (see fig. 1). Since neigh- Theprobabilitytoexitcellcintoc(cid:48),P(c(cid:48) |c),ispro-
bouring pseudo-grains share the mid-point of the edge portional to the area of the throat between them. To
betweenthem,thetetrahedrastructureistopologically reduce finite-size effects, we wrap the sample around
homeomorphictotheoriginalstructure.Thevoidspace with periodic boundary conditions and let the particle
surrounded by the pseudo-grains is still cellular, but travel larger distances by re-entering the sample. This
a cell surface now consists of triangular facets – the means that the same cell may occur at different loca-
faces of the pseudo-grains surrounding it – and throats tions along the random walk. To avoid distorting the
– the skewed polygons remaining of the original cell statistics by using the same cell too many times, we
membranes. The membranes over the throats are dis- stop the process once a cell has occurred at more than
regarded, resulting in an open-cell porous structure, in five different locations.
which the throats are the openings between neighbour- Weemphasisethatwedonotsimulatethediffusion
ing cells. The pseudo-grains volumes are smaller than withinthecell–eachstepinoursimulationcorresponds
those of real convex grains, which curve out into the to a transition of the particle from one cell to another,
cells of this structure. This increases the pore volume and the waiting time associated with the step is calcu-
and forms a limiting case, which establishes the valid- lated from the cell properties. This process constitutes
ity of our results for any porous medium, as will be a random walk on a graph, whose nodes are the cell
discussed in the concluding section. centres. After waiting for a period of time t in cell c,
c
4 ShaharAmitai,RaphaelBlumenfeld
determined by eq. (1), the particle makes a step to a Assumingthatallcellshavethesameeffectivepotential
neighbour cell, according to P(c(cid:48) |c). U, we get:
Before continuing, it is instructive to put the prob-
∆E v2/3
lem in thermodynamic context. The cell can be re- c
=Const.−lnT +ln . (4)
gardedasapotentialwellofheight∆E,andtheproba- kT Sc(r)
bility to escape from it is P(t)=P e−t/τ. We can then
0 Thisequationestablishestheheightoftheeffectivebar-
use Kramer’s escape rate formula,
rier in terms of the cell volume and the fraction of its
2πkT surface through which the particle can escape.
τ = e∆E/kT , (2)
(cid:112) Notethattheparticle’smeanfreepathwithinacell
d U”(a)U”(b)
is assumed to be well smaller than the cell size, regard-
where k is the Boltzmann constant, T is the tempera-
lessoftheparticlesize,andthereforethatKnudsendif-
ture and U”(a) and U”(b) are the second derivative of
fusion[21,22]neednotbeconsidered.However,evenif
thepotentialatthebottomandtopofthewell,respec-
thisassumptionisnotborneout,thiswouldonlymod-
tively. Interpreting t as the half-life of the particle in
c ifythecoefficientdineq.(1),whichisarbitraryanyway
the cell, t = τln2, we can combine eq. (1) and (2) to
c in our simulations. Also note that the above assump-
get: tion, P(t)∼e−t/τ, means that typical escape times do
2ln(2)πkT 1 (cid:18)3v (cid:19)2/3 not deviate much from the mean or half-life time. This
e∆E/kT = c . (3) justifies our choice of taking t as a representative.
(cid:112)U”(a)U”(b) 2Sc(r) 4π c
103 3 100
β 2.18
2x1100›fi12 P(t)0 0t 0 P(l)li12 CP(t)ti·111000−−−11550 ∼ 2x 111›fi000123 α∼1.00
α∼1.00 10−20 100 α 0.63
β 4.05 ∼
100 10−2 10t−1 100 0.0 0.4 0l.8 1.2 1.6 10−25101∼102 10t3 104 105 10−1 102 103t 104 105
i i
(a) (b) (c) (d)
2δ(t,T)10111−›fi0001012 ∼t1.00 ∼t0.65 2δ(t,T)1011−›fi00101 ∼T0 P(ξ)1231233 rrr===011...9235r50rr0 Bparameter0000....1234
10−2 10−2 ∼T0.1 12 0 E0.0
101 102 t 103 104 103 T104 105 0 1 ξ 2 0.6 0.8α 1.0
(e) (f) (g) (h)
Fig. 2: Statistical analysis of diffusion in the porous sample. All averages are over 1000 walks. (a) MSD vs. time for a small
particle (r = r0/100), for which α ∼ 1, indicating normal diffusion. Inset: overall waiting-time distribution of the walks. (b)
Step-length distribution, Pl(li), calculated directly from the sample (see text), with a normal fit. Pl is nearly constant for
all particle sizes, as the cells positions are fixed. P is narrow, and thus does not affect the universality class of the walk. (c,
l
d) Waiting-time distributions (left) and MSD vs. time (right) for large particles (r/r0 = 1,...,1.35). Both α and β decrease
with increasing r. All power-law relations hold for at least 1.5 orders of magnitude. (e, f) Ensemble-averaged TAMSD, (cid:104)δ2(cid:105),
vs. t (left) and vs. T (right) for the same large particles. For t > 102.5, (cid:104)δ2(cid:105) becomes sub-linear in t, like the MSD, and is
nearlyindependentofT.(g)PDFoftheamplitudescatter,ξ,from2000individualwalksofthreedifferentparticlesizes.The
distribution broadens dramatically with particle size, which indicates the walk is non-ergodic. (h) The EB parameter for the
same three particle sizes, as well as for r=r0/100, vs. the anomaly parameter, α. The solid line is the expected relation for
CTRW, which goes through all points. Error bars denote 2σ ranges, and were established by the variance of 10 independent
measurements.
Modifyingcontinuous-timerandomwalks... 5
For each walk we calculate the particle’s position
at time t, relative to the origin, x(t)=(cid:80)nN=(t1)ln, where 1.1
r=1.1r
N(t)isthenumberofstepsmadebeforetimet.Wethen 0
calculate the MSD, (cid:104)x2(cid:105), as a function of t, where the 1.0
angularbracketsdenoteaverageover1000walks.Figure
2a shows (cid:104)x2(cid:105) vs. t for a particle of size r = r0/100. 0.9 r=0.95r0
The linear relation indicates normal diffusion, with a
r=1.25r
diffusioncoefficientofD =(0.65±0.03)d (d=100d0). α0.8 0
The inset figure shows a narrowly bounded P .
t
0.7
r=1.35r
Largeparticles Wenextconsiderparticlesofsizescom- 0
parable to r . Such particles diffuse differently due to 0.6
0
two effects. One is delay and trapping inside cells. The Percolationthreshold
larger r is the lower the probability of passing through 0.51.5 2.0 2.5 3.0 3.5 4.0 4.5
anyparticularthroat,sincetheeffectivearea,whichthe β
particle can pass through, is smaller. This reduces the
Fig.3:Twosetsofsimulations–diffusioninaporoussample
overall probability to exit a cell, increasing the waiting (bluedots)andCTRWinunconfinedspace(greentriangles).
times spent inside cells. As a result, while the waiting Each simulation (dot) represents a particular particle size,
times are narrowly bounded for a small particle, P (t ) and is described by the power-law of the waiting-time dis-
t i
tribution, β, and the anomaly parameter, α. Small (large)
developsapower-lawtailforsufficientlylargeparticles,
particles appear at the top right (bottom left) corner of the
Pt(ti >t(0))∼t−i β, with β a function of r. The second graph – they experience a narrow (wide) waiting-time dis-
effect is that the topology changes with particle size; tributionandundergonormal(sub-)diffusion.Theblueand
as it increases, the probability of passing through some greenlinesarefitsoftheformα=1−21exp{−(β−βt−21)/τ},
throatsvanishesidentically,changingthesystem’scon- with (βt(sample) = 2.53, τ = 0.42) and (βt(CTRW) = 2.01,
τ =0.40), respectively. The theoretical CTRW prediction of
nectivity for this particle.
a universality class transition at β = 2 is denoted by the
Fig. 2c shows Pt for a few large particles. We see black dashed line. The red lines mark the range of β’s that
that β decreases with r, corresponding to longer wait- corresponds to the percolation threshold of the sample. This
ingtimes.Fig.2dshowstheMSDvs.timeforthesame range matches the universality-class transition for confined
diffusion. Error bars denote 2σ ranges, and were established
particles. We see that beyond a certain particle size
bythevarianceof20independentmeasurements.
(r ∼1.2r )αstartstodecrease–thediffusionbecomes
0
anomalous.Toquantifythisrelation,wechoosealarger
set of radii and plot α(r) vs. β(r) (blue dots in fig. 3). broadening of the PDF as r increases indicates that
We see that for smaller particles β (cid:29) 2 and α = 1, the process is non-ergodic and that there is ageing [5].
correspondingtonormaldiffusionwithnarrowwaiting- Fig. 2h shows that the corresponding EB parameters,
time PDFs. As r increases, β decreases and the walks i.e. the variance of these PDFs, follow the theoretical
eventually become sub-diffusive with α < 1. We mea- expectation from CTRW.
sure a transition at β(sample) =2.53±0.03.
t
A short comment on scaling windows is due – α is
evaluated along t ∈ (103.5,105) (see fig. 2d). At much 4 Modelling the diffusion as isotropic CTRW
longertimes,thediffusionisnormalforallparticlesizes.
This is merely a consequence of the finite system size – Next we show that an attempt to simulate this pro-
therearenocellswithlongerwaitingtimesthan∼105. cess with straightforward isotropic CTRW fails. This
To investigate the diffusion process further, we cal- may not come as a surprise, as the statistical analysis
culatetheTAMSD,δ2,anditsaverageover1000walks, showed some deviations from the traditional CTRW,
(cid:104)δ2(cid:105). Fig. 2e and 2f show (cid:104)δ2(t,T)(cid:105) vs. the time lag, t, but it is still constructive to describe the CTRW simu-
andvs.theoveralltrajectorytime,T,respectively.(cid:104)δ2(cid:105) lationtobetterunderstanditsadjustmentsinsection5.
issub-linearint,deviatingfromthelineart-dependence Forsuchasimulationweusether-dependentPDFs,P
l
intheCTRWmodel,anditisnearlyindependentofT. and P , that the particle experiences while diffusing in
t
This behaviour is the same as for a random walk on the confined structure. Recall that these PDFs refer to
a fractal. In contrast, the amplitude scatter, ξ, follows the transition between cells, rather than the movement
the CTRW prediction: Fig. 2g shows the PDF of ξ for within a cell. One way to obtain these PDFs is to mea-
three different particle sizes of order r . The dramatic sure them empirically during a diffusion process. How-
0
6 ShaharAmitai,RaphaelBlumenfeld
ever, more efficient is to compute them directly from Collecting the statistics of all possible steps in the
the structural statistics of the porous medium, as we sample, we get a PDF described well by the Gaussian
outline next. P (l )∼exp{−(l −¯l)2/2σ2},withσ/¯l=0.15±0.02(see
l i i
A derivation of Pl and Pt from structural statistics fig.2b).AsPl isnarrow,itdoesnotaffecttheuniversal-
should be made cautiously because a straight-forward ityclassofthewalk.Pl staysnarrow,andindeednearly
histogram of the waiting times of all cells in the sam- constant, for all particle sizes, as the cells positions are
ple ignores the inherent correlation between the prob- fixed. Note that it is also possible to derive Pl analyt-
ability to visit a cell, P(c), and the waiting time [23]. ically, given the cell volume distribution and nearest
Cellsthataredifficulttogetoutof(longwaitingtimes) neighbour volume-volume correlations. This, however,
tend to have a lower probability of getting into and are is somewhat downstream from the main thrust of this
thereforevisitedlessfrequently.Inparticular,somecells paper.
arecompletelyinaccessibleforparticlesaboveacertain Wearenowabletoobtainaccuratestep-lengthand
size. waiting-time PDFs for every particle size, which could
To this end we use the observation that, to a very be used as input into unconstrained CTRW. Fig. 3
good accuracy in our diffusion process, P(c) is linear shows results of CTRW simulations (green triangles)
in the cell’s total effective throat area, A(t). This ob- usingthesePDFs.Thefittothiscurve(solidgreenline)
c
differs from the theoretical prediction [24], α = β−1,
servation, which is independent of the cell volume, can
duetofinitetimeandsizeeffects.Ourmeasuredtransi-
be seen over 3.5 orders of magnitude in fig. 4. Further-
tion for CTRW is β(CTRW) =2.01±0.02, in agreement
more, this holds for all particle sizes, both well smaller t
with the theoretical prediction.
and comparable to r . This allows us to estimate P(c)
0
for a particular particle size by using the effective A(t) A key observation is that the CTRW exhibits a
c
normal-to-anomalous transition at a lower values of β
– a direct structural characteristic of the medium. In
than in the actual sample. Since both processes have
principle,oneexpectsthevisitingprobabilitytobecor-
the same step-length and waiting-time distributions,
related with the visiting probabilities of neighbouring
but one is performed on a graph and the other in free
cells, but fig. 4 shows that this effect is negligible.
space,thediscrepancymuststemfromtheconnectivity
of the sample, which the CTRW model cannot account
for. This is supported by the statistical analysis of the
10−2 diffusion in the porous sample (section 3), that show
behaviours typical to random walks on a fractal.
Asmentionedabove,thesizeofthediffusingparticle
10−3 P(c)∼Ac(t) determinestheeffectivethroatsizes,andhencethecon-
nectivity of the porous sample. Moreover, above some
size there is no path percolating between the sample’s
)
(c10−4 boundaries. Fig. 5 shows the dependence of the per-
P
colating accessible volume on r, where a range of radii
aroundthepercolationthresholdismarkedbyredlines.
10−5 Thesamerangeismarkedinfig.3.Weseethattheuni-
versality class transition in the sample occurs within
this range. This is more evidence that the connectivity
10−6 10−3 10−2 10−1 100 plays an important role in determining the universal-
A(t) ity class. Specifically, around the percolation threshold
c
the incipient cluster assumes a fractal-like structure,
Fig. 4: Visiting probability, P(c), vs. the cell’s total effective further inducing sub-diffusion [6, 7]. We conclude that
throatarea,A(ct),inadiffusionprocesswith r=r0. our simulations of diffusion in porous media are best
described as CTRW on a percolation cluster.
We can now estimate more accurately P , and, par-
t
ticularly,thepower-lawβ,foreveryparticlesize,byma- 5 Anisotropic CTRW
nipulating the waiting-time histogram as follows. For
everybinconsistingofthewaitingtimesofcells{c ,...,c }, To retain the usefulness of the CTRW model, it would
1 n
wemultiplythebin’sheightby(cid:80)nP(c )andnormalise bedesirabletomodifyittocapturetheeffectofconnec-
1 i
thehistogramintoaPDF.Thissuppresseslongwaiting tivity.Wenextproposesuchamodification,inspiredby
times in P . the PDF of the angle between successive steps in the
t
Modifyingcontinuous-timerandomwalks... 7
1.0 1.0
0.8 0.8
)
)
100 s(θ0.6 ack0.6
%) 90 co0.4 Pb0.4
( 80 P(
e 0.2 0.2
siz 70 0.0 0.0
1.0 0.0 1.0 0.4 0.8 1.2
er 60 − cos(θ) r/r
t 0
s 50
u
cl 40 Fig. 6: Left: the PDF of cos(θ), with θ the angle between
t
n 30 successive steps of a diffusion process in the porous sample
e
pi 20 for r=r0/100. The singular point at cos(θ)=−1 marks the
ci finiteprobabilityofabackwardstepP =0.087±0.001.In
n 10 back
I addition,P(−1<cos(θ)<−0.88)=0.Right:thedependence
0 of P on r. For large particles, backward steps dominate
back
thewalk.
0.5 1.0 1.5 2.0 2.5
r/r
0
Fig. 5: The percentage of cells belonging to the incipient- thesample.Testingthismethodbyasetofsimulations,
clustervs.theparticleradius(inunitsofr0).Thepercolation we find that the universality class transition occurs at
rangeismarkedinred.Thecorrespondingβ rangeismarked
β(anisotropic) = 2.1±0.03 – closer to the one measured
infig.3. t
in the porous sample.
The second method is more drastic: at every step
oftheCTRWweconstructanewcell,accordingto the
porous sample (fig. 6). There is a finite probability to
structural statistics of the porous sample. We choose
make a backward step, i.e. go back through the throat
thecell’svolume,thenumberofthroats,andthethroats’
the particle entered a cell, P ≡P (θ =π). For very
back θ
areasfromthecorrespondingdistributions,derivedfrom
small particles P = 0.087±0.001. This is in con-
back
thesample.Usingtheseandtheparticlesize,wecalcu-
trasttotheconventionalCTRWmodel,wherethenext
late the waiting time for the new cell. Then we choose
step direction is uniformly distributed. Moreover, we
randomly an accessible exit throat. The probability to
see that P >1/13.7≈0.073, which is the inverse of
back
exit through the throat is proportional to its effective
the average number of throats per cell, and is the ex-
area. Theparticle thenmakesa stepin the directionof
pected value for P when all steps are equiprobable.
back
theexitthroat.Anothercellisthenconstructedaround
This is because the mean size of an entrance throat is
it. The step length is the sum of the radii of these two
larger than the mean size of all the throats. The en-
cells. This process is then iterated as many times as
hanced backward step probability can be regarded as
required.
a correlation between successively visited throats. As r
increases, the total available throat area decreases and A key feature of this method is that the particle
P increases, as can be seen in the right panel of fig. ‘remembers’ the exit direction, the last used cell and
back
6. For very large particles, P dominates the walk. the exit throat. The latter is used as one of the throats
back
This can be seen as the ‘lowest order’ effect of connec- in constructing the next cell. If this throat is chosen
tivity, which we next try to capture within CTRW. again then the last step is retraced. We refer to this
method as DA for the two variables that the particle
We examine two methods to modify the CTRW
‘remembers’ – step direction and throat area.
model, both introducing a backward bias. In the first,
Within the DA process, the particle may pass back
weaddtoP andP athirddistribution,P ,ofstepdi-
l t θ
and forth several times through a large throat before
rectionrelativetothepreviousstep.Inthismethod,the
it moves on and loses memory of this throat. This pro-
particle ‘remembers’ the direction of the previous step,
cesscorrelatesnotonlysuccessivesteps,astheprevious
andarelativeangleischosenfromthenon-uniformdis-
method does, but also successive backward steps.
tribution P (θ), e.g. the one in fig. 6. P can be calcu-
θ θ
lateddirectlyfromtheporousstructureforanyparticle Using the DA model, we obtain a universality-class
radius r, similarly to P and P (see section 4). A ran- transitionatβ(DA) =2.50±0.03,inexcellentagreement
l t t
domstepisthenmadeattheangleθrelativetothepre- with the original simulation results. Thus, this model
vious step direction. The waiting time and step length captures much better the particle size-driven topologi-
at each step are chosen as before. This method intro- cal change.
duces a correlation between consecutive steps, which is AnimportantfeatureoftheDAmodelisthatP
back
presentunavoidablyinthediffusionoflargeparticlesin ishigherthanintheporoussample.Tounderstandthis,
8 ShaharAmitai,RaphaelBlumenfeld
considertwocellsintheporoussample,connectedbya
large throat, and connected to the rest of the structure 1.0
bysmallerthroats.Oncetheparticleentersoneofthese Sample
cells,itislikelytomovebackandforthseveraltimesbe-
0.8 CTRW
fore it emerges. However, the smaller the throats lead-
Method1
ing to this pair of cells, the less likely is the particle
n
to enter in the first place. This means that the occur- tio0.6 Method2(DA)
a
rence frequency of such sub-structures, which increase el
P , is suppressed. In contrast, once a particle passes rr0.4
back o
C
througharelativelylargethroatintheDAmodel,then,
other than this throat, an entire new cell is generated 0.2
for each step. This results in a higher probability that
theparticleoscillatesacrosssuchathroat.Thisfeature
0.0
appears to compensate for other, more complex, miss-
0 2 4 6 8 10 12 14
ingtopologicalfeatures,makingtheDAabettermodel Stepdistance
for the diffusion process.
Fig. 7: Step-length correlation function for the four types of
Asafurtherinvestigationoftheproposedmethods,
simulations discussed in the paper, all with r =1.2r0. Solid
we present the step-length correlation function for the lines are decaying exponentials. The correlation function of
different models, all using r =1.2r (fig. 7). The corre- thefirstadjustmenttoCTRWdecaysexponentially,andthe
0
lationinthesampleismainlyduetothefactthateach secondadjustment(DA) –moreslowly.
two consecutive steps enter and exit a certain cell, c.
If c is small, then the two steps will tend to be short,
from regular to anomalous diffusion. We showed that
and vice versa, leading to positive correlation. In addi-
the the two models give different results – while the
tion, a high P means that many consecutive steps
back CTRWsimulationsfollowthetheoreticalprediction,up
are of exactly the same length. This further increases
to finite-time effects, with a transition to sub-diffusion
the correlation for larger particles. As expected, the
at β ≈ 2, the diffusion simulations in the confined ge-
traditional CTRW exhibits no correlations. Our first
ometry exhibit a transition at β ≈ 2.5. We established
adjustment to CTRW, introducing the possibility of a
that the difference stems from changes to the effective
backward step, adds correlation. However, since this is
connectivity available to the particle with increasing
onlyaone-stepcorrelationitdecaysexponentially.The
size. This particle size-driven change in connectivity is
increaseincorrelationwithintheDAprocessisbecause
not taken into consideration in the CTRW model. We
oftheincreasedprobabilityforalongsequenceofback-
supported this conclusion by investigation of the time
wardsteps,asdiscussedabove.Asaresult,theDAcor-
average MSD and by showing that the transition in
relationfunctiondoesnotdecayexponentially,agreeing
universality class occurs at the same range of parti-
better with the diffusion simulation in the sample. The
cle sizes that corresponds to the percolation transition.
correlationoftheDAatastepdistanceofoneishigher
Our findings show unequivocally that the discrepancy
thanthatinthesimulateddiffusionprocessbecauseits
betweenthetwomodelsisnot duetodifferentwaiting-
P is higher, yet the DA correlation decays faster
back timedistributions,sinceusingidenticaldistributionsdo
with the number of steps since it lacks the more in-
lead to different universality classes.
volved topological correlations.
It is important to comment on the range of valid-
ityofourresults.Increasingtheparticlesizecanbere-
6 Conclusions garded,alternatively,asshrinkingtheporousstructure,
while keeping the particle size unchanged. Evidently,
Toconclude,wecomparedbetweentwonumericalmod- the smaller the pore space the more restricted the dif-
elsofdiffusionofafinitesizeparticleinporousmedia:a fusingparticleisandthelargerthediscrepancybetween
directsimulationofthediffusionprocessinacomputer- thesimulateddiffusionandtheCTRW.Thus,bystart-
generatedsample,andwhatiscommonlybelievedtobe ing from a medium with a very large pore space, we
an equivalent CTRW. We first presented a method to established that our results hold true for a wide range
construct the representative PDFs of the step-length of porous media with lower porosity.
andwaiting-time,P andP ,giventheparticlesizeand Wong et al. [10] studied experimentally a related
l t
statisticalinformationaboutthestructureoftheporous process of trace particles diffusing in biological net-
media alone. Using the same particle size dependent works of entangled F-actin filaments. There, the dif-
P and P in both models, we analysed the transition fusion of particles, of size comparable to the typical
l t
Modifyingcontinuous-timerandomwalks... 9
network mesh size, is sub-linear and P decays alge- 5. R. Metzler, J.H. Jeon, A.G. Cherstvy, E. Barkai,
t
braically. The universality class they observe is a func- Phys. Chem. Chem. Phys., 16(44), 24128 (2014)
tion of the particle-to-mesh size ratio, in agreement 6. Y. Gefen, A. Aharony, S. Alexander, Phys. Rev.
with our results. However, they do not observe size- Lett. 50, 77 (1983)
driven topological effects and their values of α and β 7. R.Pandey,D.Stauffer,A.Margolina,J.Zabolitzky,
are in good agreement with the CTRW model. This is J. Stat. Phys. 34(3-4), 427 (1984)
because of the flexibility of the gel-like network, which 8. B. Berkowitz, A. Cortis, M. Dentz, H. Scher, Re-
allowstrappedparticlestoeventuallyescapebydeform- views of Geophysics 44(2) (2006)
ing the filament network, a phenomenon also modelled 9. B. Bijeljic, M.J. Blunt, Water Resources Research
recentlyin[25].Therigidityofthestructureconsidered 42(1) (2006)
hereprecludesthisparticleescapemechanismandisthe 10. I.Y. Wong, M.L. Gardel, D.R. Reichman, E.R.
reason for this difference. A potentially related process Weeks, M.T. Valentine, A.R. Bausch, D.A. Weitz,
of large particles, diffusing in rigid porous media, was Phys. Rev. Lett. 92(17), 178101 (2004)
studied experimentally in [26]. That study focused on 11. P. De Anna, T. Le Borgne, M. Dentz, A.M. Tar-
hydrodynamic in-pore effects, which might be interest- takovsky, D. Bolster, P. Davy, Phys. Rev. Lett.
ing to eventually include in our model. 110(18), 184502 (2013)
The main advantages of the CTRW model are that 12. S.M.A. Tabei et al., PNAS 110 4911 (2013)
it overcomes potential finite-size problems and is less 13. A.V. Weigel, B. Simon, M.M. Tamkun, D. Krapf,
demandingcomputationally.However,aswehavedemon- PNAS 108 6438 (2011)
strated here, this is achieved at the expense of ignor- 14. J.H. Jeon et al., Phys. Rev. Lett. 106, 048103
ingtopologicalinformationaboutlocalconnectivity.To (2011)
preserve the advantages of CTRW, these need to be 15. E. Yamamoto, T. Akimoto, M. Yasui, K. Yasuoka,
takenintoconsideration.Tothisend,weintroducedtwo Scientific Reports 4 4720 (2014)
anisotropicCTRWmodels.Oneincludesmemoryofthe 16. R. Blumenfeld, S.F. Edwards, S.M. Walley, Gran-
last step direction and a non-uniform distribution of ular systems, in The Oxford Handbook of Soft
stepdirection.Theother,theDAmodel,addsmemory Condensed Matter, Eds. E.M. Terentjev and D.A.
of the area of the last throat visited, effectively corre- Weitz, (Oxford University Press, Oxford, UK,
lating successive backward steps. The DA model shows 2015)
auniversality-classtransitionatβ(DA) =2.50±0.03,in 17. K.A. Brakke, Experimental Mathematics 1(2), 141
t
good agreement with the one measured in our simula- (1992)
tions of the diffusion process in the confined geometry 18. Y. Wang, S. Neethling, Minerals Engineering
ofaporousmedium.WeconcludethattheDAmodelis 19(10), 1069 (2006). Selected papers from Compu-
a better alternative to the traditional CTRW for mod- tational Modelling 05, Cape Town, South Africa
elling diffusion of finite size particles in such media. It 19. G. Frenkel, R. Blumenfeld, P.R. King, M.J. Blunt,
combinestheCTRWadvantages,overcomingthefinite- Adv. Eng. Materials 11(3), 169 (2009)
nessofthesampleandconvenienceofapplication,with 20. R. Blumenfeld, S.F. Edwards, The Euro. Phys. J.
a better capturing of topological effects. E: Soft Matter and Biological Physics 19(1), 23
(2006)
Acknowledgement: SA is grateful for support from 21. M. Knudsen, Ann. Phys. 4F 28, 75 (1909)
the Alan Howard Scholarship. 22. P. Clausing, Ann. Phys. 5F 4, 567 (1930)
23. Y. Edery, A. Guadagnini, H. Scher, B. Berkowitz,
Conflict of interest: The authors declare that they Water Resources Research 50(2), 1490 (2014)
have no conflict of interest. 24. E.W. Montroll, G.H. Weiss, J. Math. Phys. 6(2),
167 (1965)
25. A. Godec, M. Bauer, R. Metzler, New J. Phys. 16,
References
092002 (2014)
26. M.JSkaug,L.Wang,Y.Ding,D.K.Schwartz,ACS
1. H.Scher,E.W.Montroll,Phys.Rev.B 12(6),2455
Nano 9 2148 (2015)
(1975)
2. H. Scher, M.F. Shlesinger, J.T. Bendler, Physics
Today 44(1), 26 (1991)
3. B.B. Mandelbrot, The fractal geometry of nature,
vol. 173 (Macmillan, 1983)
4. L.P. Kadanoff, Physics 2(6), 263 (1966)