Reaction Front in an A+B -> C Reaction-Subdiffusion Process

0.35 MB·English

by  S. B. Yuste

Reaction Front in an A+B → C Reaction-Subdiffusion Process
S. B. Yuste†, L. Acedo†, and Katja Lindenberg‡
†Departamento de F´ısica, Universidad de Extremadura, E-06071 Badajoz, Spain
‡Department of Chemistry and Biochemistry, and Institute for Nonlinear Science,
University of California San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0340, USA
4
0 WestudythereactionfrontfortheprocessA+B→C inwhichthereagentsmovesubdiffusively.
0 Ourtheoreticaldescriptionisbasedonafractionalreaction-subdiffusionequationinwhichboththe
2 motion and the reaction terms are affected by thesubdiffusive character of the process. We design
numerical simulations to check our theoretical results, describing the simulations in some detail
n
because the rules necessarily differ in important respects from those used in diffusive processes.
a
J Comparisons between theory and simulations are on the whole favorable, with the most difficult
quantities to capture being those that involve very small numbers of particles. In particular, we
7
analyzethetotalnumberofproductparticles,thewidthofthedepletionzone,theproductionprofile

] of product and its width, as well as the reactant concentrations at the center of the reaction zone,
h all as a function of time. Wealso analyze theshape of the product profile as a function of time, in
c particular its unusualbehavior at thecenter of thereaction zone.
e
m
PACSnumbers: 82.33.-z,82.40.-g,02.50.Ey,89.75.Da
-
t
a
t I. INTRODUCTION basis of scaling arguments, later made more rigorous, a
s
number of exponents were deduced to characterize this
.
t evolution. The firstexperimentsconfirmingthese results
a It is very well known that diffusion-limited binary re-
m were carried out with species A and B diffusing in a gel
actions in low dimensions may lead to the spontaneous
contained in a thin capillary, with one species initially
- appearance of spatial order and spatial structures, and
d occupying one side, the other occupying the other side,
to associated “anomalous” rate laws for the global den-
n andasharpfrontbetweenthem[24]. Morerecentexper-
sities of the reacting species. For example, the reac-
o iments havebeen carriedout in gel-free systems [25, 26].
c tions A + A → C (some selected references of many
Theevolutionofthereactantfrontsandoftheproductof
[ in the literature are [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
the reactioninthese experimentsboth reflectthe kinetic
and A + B → C (again, some of many references
1 anomalies.
are [1, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23])
v
In this paper we extend the front analysis to subdif-
1 under “normal” circumstances are described by second-
fusive reactions. Subdiffusive motion is characterizedby
7 order rate laws, whereas the asymptotic rate law for the
a mean square displacement that varies sublinearly with
0 former reactionis ofapparentorder(1+2/d)for dimen-
1 sions d<2, and for the mixed reaction it is of apparent time,
0
order (1 +4/d) for d < 4. The slow-down implied by 2K
4 the higher order is a consequence of the rapid deviation hr2(t)i∼ γ tγ, (1)
0 Γ(1+γ)
of the spatial distribution of reactants from a random
/
at distribution. This is in turn a consequence of the fact with 0 < γ < 1. For ordinary diffusion γ = 1, and
m that diffusion is not an effective mixing mechanism in K ≡ D is the ordinary diffusion coefficient. We argue
1
low dimensions. for the importance of this generalization on a number
-
d To design an experiment in a constrained geometry of grounds. First, there exists a huge literature on sys-
n in order to measure these anomalies is not at all sim- temsthatdeviatefromdiffusivebehaviorandareinstead
o
ple, especially for the mixed reaction [24, 25, 26]. It characterizedby motion all the way from subdiffusive to
c
is simpler for the A + A problem because a number superdiffusive [39, 40]. Subdiffusive motion is particu-
:
v of appropriate non-chemical species can be identified larly important in the context of complex systems such
i
X that essentially undergo the simplest annihilation reac- as glassy and disordered materials, in which pathways
tion or variants thereof. Examples include exciton an- are constrainedfor geometric or for energetic reasons. It
r
a nihilation experiments in one-dimensional pores and in is also particularly germane to the way in which exper-
effectively one-dimensional polymer wires [15], excited iments in low dimensions have to be carried out. Such
molecule naphthalene fusion and quenching experiments experiments must avoid any active or convective or ad-
in one-dimensional pores [27], and kink-antikink simula- vective mixing so as to ensure that any mixing is only a
tions in one dimension [28]. Experimental observations consequence of diffusion. To accomplish this usually re-
of the A+B anomalies instead generally involve reac- quirestheuseofgelsubstratesand/orhighlyconstrained
tion fronts. Early on G´alfi and Ra´cz [29] and later oth- geometries(thefirstgel-freeexperimentswerecarriedout
ers [30, 31, 32, 33, 34, 35, 36, 37, 38] recognized that recently [25, 26]). Under these circumstances it is not
the kineticanomaliesinthe homogeneoussystemswould clear whether the motion of the species is actually dif-
be reflected in the evolution of reaction fronts. On the fusive, or if it is in fact subdiffusive. Indeed, a recent
2
detailed discussion on ways to extract accurate parame- Sec.III we discuss our Monte Carlosimulationmethods.
tersandexponentsfromsuchexperimentsconcludesthat Section IV is a compendium and comparison of numer-
at least the experiments presented in that work, carried ical and theoretical results. Some closing comments are
outinagel,reflectsubdiffusiveratherthandiffusivemo- presented in Sec. V.
tion [41].
We recently solved the A + A reaction-subdiffusion
probleminonedimension[42]. Tosolvethisproblem,we II. THE MODEL
generalizedmethodsfirstappliedtothereaction-diffusion
A+Aproblem. FortheA+A→Aproblemthemethodof We start with a system of A particles on one side and
intervalsallowsanexactformulationintermsofintervals B particles on the other of a sharp linear front, defined
onthelinethatareemptyofAparticles[4,6,7,9,10,11]. to lie perpendicular to the x axis. The particles diffuse
The distribution of intervals evolves linearly, and there- and react with a given probability upon encounter. A
fore one can find an exact solution. In the reaction- standard mean-field model for the evolution of the con-
diffusion problem the description involves a diffusion centrations a(x,t) and b(x,t) of A and B particles along
equation, while the reaction-subdiffusion problem in- x is given by the reaction-diffusion equations
volves a subdiffusion equation [42]; both can be solved
∂ ∂2
exactly. For the A+A → C problem one uses instead a(x,t) = D a(x,t)−ka(x,t)b(x,t)
the odd/even parity method [7, 8], whereby one keeps ∂t ∂x2 (2)
trackofthe parityofthe number ofparticlesinaninter- ∂ ∂2
b(x,t) = D b(x,t)−ka(x,t)b(x,t),
val. The associated distribution again satisfies a linear ∂t ∂x2
diffusion or subdiffusion [42] equation.
where D is the diffusion coefficient assumed to be equal
Before these exact methods were developed, it was
for the two species. The initial conditions are that
customary to model these and other binary reactions
a(x,t) = const = a for x < 0 and a(x,t) = 0 for
0
by writing down a reaction-diffusion equation for each
x ≥ 0. Similarly, b(x,t) = const = b for x > 0
0
species in the reaction. Such equations typically contain
and b(x,t) = 0 for x ≤ 0. With these conditions,
a diffusion term and a binary reaction term, the latter
no matter the dimensionality of the system, the sys-
being some truncated form of a two-particledistribution
tem of equations is effectively one-dimensional. The
function. For instance, in the A+B reactionthe typical
front problem was first analyzed via a scaling descrip-
reaction term simply involves the product of the local
tion [29] and later refined by a large number of authors
concentrations of reactants, −ka(r,t)b(r,t) [1]. In the
usingmorerigoroustheoreticalandcarefulnumericalap-
A+Areactiononehastobeslightlymorecarefulbecause
proaches [30, 31, 32, 33, 34, 35, 36, 37, 38]. One upshot
in writing, for example, −ka2(r,t) one must be careful
oftheextensiveworkisthatd=2isacriticaldimension
not to include spurious self-reaction contributions [5, 8].
for the mean field description to be appropriate. Below
Once the A+A exact models were developed that did
d=2 one must take into account fluctuations, neglected
not require one to explicitly write a reaction term, it
in this description, that completely change the outcome
was possible to analyze the accuracy of the approximate
of the analysis. A particularly transparent argument for
truncations [43]. Also, it was not necessary to consider
this critical dimension was provided by Krapivsky [37].
the generalization of the reaction term to the subdiffu-
He argued that the reaction constant in the mean-field
sive case since the exact methods could be generalized
reaction rate r = kab should in general depend on the
directly.
diffusion constant D and the radius R of the reacting
The situationismorecomplicatedforthe A+B prob- particles. Dimensional analysis gives k ∼ DRd−2, but
lem,becausenosuchexactformulationsorsolutionshave on physical grounds one expects the reaction rate con-
beendevelopedinthiscase. Thereisalargeliteratureon stant to be an increasing function of the radius R. The
the reaction-diffusion problem with different truncation conclusionisthatthemeanfieldmodelcanthereforenot
schemestorepresentthereactionterm,buttheliterature be valid for d < 2. While it has been assumed that the
on the reaction-subdiffusion problem is far more recent meanfield model holds for the criticaldimensions d=2,
and relatively unsettled. In particular, at the current Krapivsky finds logarithmic corrections that have also
stage of development of this problem it is necessary to been observed in simulations [31]. In our analysis and
think about how to (approximately) model the reaction simulationswewilltaked=2(whichturnsouttobethe
term. critical dimension for the subdiffusive problem as well)
In Sec. II we present a discussion of the model to be and will therefore not deal with the lower-dimensional
usedforthedescriptionoftheA+Breaction-subdiffusion fluctuation effects. In this first study we will not deal
problem. Having arrived at a particular set of fractional with logarithmic corrections.
equations. We apply a scaling theory to these equations Inordertogeneralizethereaction-diffusionproblemto
akinto thatofG´alfiandRa´cz[29],but nowforasubdif- reaction-subdiffusion,wemustdealwiththesubdiffusive
fusive front. To support the theoretical conclusions,it is motionoftheparticles(generalizationofthefirsttermin
necessarytoperformnumericalsimulations,whichisnot Eq. (2)) and with their reaction rate law (second term).
a trivial matter for a problem involving subdiffusion. In We discuss each separately.
3
Subdiffusion is not modeled in a universal way in the for the concentrations and
literature. Among the more successful approaches to
R (x,t)=t−µRˆ (xt−α) (7)
thesubdiffusionproblemhavebeencontinuoustime ran- γ γ
dom walks with non-Poissonian waiting time distribu-
for the reaction term. The exponents θ, α, and µ are to
tions [44, 45, 46], and fractional dynamics approachesin
bedeterminedfromthreerelations. Thescalingformsare
which the diffusion operator is replaced by a generalized
only valid for x≪W , that is, well within the depletion
d
fractional diffusionoperator[39,45,47,48,49]. Therela-
zone.
tionbetweenthetwohasalsobeendiscussed[39,45,48].
Twoofthe threerelationsneededtofixthe scalingex-
In particular, the fractional dynamics formulation that
ponents do not require further specification of the reac-
leadsto the meansquaredisplacement(1)canbe associ-
tion term. Since the reaction zone increases more slowly
atedwith acontinuoustime randomwalkwithawaiting
thanthewidthofthedepletionzone(anassumptionthat
time distribution between steps which at long times be-
ex post turns out to be correct), we can focus on the
haves as
concentration difference u(x,t) = a(x,t)−b(x,t) to de-
ducethewidthofthelatter. Thereactiontermdropsout
ψ(t)∼t−γ−1. (3)
whenonesubtractstheequationsinEq.(4),anditsform
thereforedoesnotmatter atthis point. Generalizingthe
We adopt the fractional dynamics approach, and com-
procedurein[29],onecanscaletheresultingequationby
ment later on some issues associated with it that must
measuring concentrations in units of a , time in units of
carefullybe consideredinthe contextofnumericalsimu- 0
τ = 1/(ka ), and length in units of l = (K τγ)1/2, so
lations. We thus replaceEq.(2)with the setofreaction- 0 γ
that the equation is simply
subdiffusion equations
∂ ∂2
∂ ∂2 u(x,t)= D1−γ u(x,t) (8)
a(x,t) = K D1−γ a(x,t)−R (x,t) ∂t 0 t ∂x2
∂t γ 0 t ∂x2 γ
(4)
∂ b(x,t) = K D1−γ ∂2 b(x,t)−R (x,t), andthe onlycontrolparameterisq =b0/a0 inthe initial
∂t γ 0 t ∂x2 γ condition:
u(x,0) = 1 for x<0
where Kγ is the generalizeddiffusion coefficient that ap- (9)
pears in Eq. (1), and D1−γ is the Riemann-Liouville u(x,0) = −q for x>0.
0 t
operator,
The solution is
0Dt1−γf(x,t)= Γ(1γ)∂∂tZ0tdτ(tf−(xτ,)τ1)−γ. (5) u(x,t)=−q+ 1+2 qH11,,10"tγx/2 (cid:12) ((01,,1γ2)) #. (10)
(cid:12)
(cid:12)
The reaction term Rγ(x,t) ≡ R will be discussed subse- where H1,0 is the Fox H-function [4(cid:12)9, 50, 51]. When
quently,becausecertainaspectsofthe problemareinde- 1,1 (cid:12)
γ =1 this reduces to the diffusion result [29]
pendent of the specific form of this term.
1+q x
u(x,t)=−q+ erfc . (11)
2 2t1/2
A. Scaling independent of reaction term (cid:16) (cid:17)
FromEq.(10)weseethatthewidthofthedepletionzone
scales as
As the reaction proceeds, a depletion zone develops
around the front. This is the region where the concen- W ∼tγ/2, (12)
d
trations of reactants are significantly smaller than their
initial values. How the width W evolves with time is i.e., ∂a(x,t)/∂x ∼ ∂b(x,t)/∂x ∼ t−γ/2. Then, from
d
one of the measures typically used to characterize the Eq.(6), the followingrelationbetweenscalingexponents
process. Within this depletion zone lies the so-called re- follows immediately:
actionzone,the regionwhere the concentrationc(x,t)of γ
θ+α= . (13)
the product C is appreciable. This concentration profile 2
hasawidthw whosevariationwithtimeisanotherchar-
Thesecondrelationfollowsfromthefactthatthecon-
acteristic of the evolving reaction. The evolution of the
centrationgradientofAandB leadstoafluxofparticles
productionrateofC (whichdeterminestheheightofthe
towardthe reactionregion. The assumption that the re-
profile of c(x,t) in the reaction zone) is a third measure
action is fed by these particle currents then leads to the
of the process. To find these time dependences we adapt
quasistationary form in the reaction zone
the original scaling approach [29, 31] to the subdiffusive
case, and assume the scaling forms ∂2
0 = K D1−γ a(x,t)−R (x,t)
γ 0 t ∂x2 γ
a(x,t) = t−θaˆ(xt−α) (14)
b(x,t) = t−θˆb(xt−α) (6) 0 = Kγ 0Dt1−γ∂∂x22b(x,t)−Rγ(x,t),
4
which requires that memory with the reaction term. Some assume that, as
in the case of ordinarydiffusion, reactionscan simply be
µ=θ+2α+1−γ . (15) modeled by a space-dependent form of the law of mass
action,e.g., by setting R=ka(x,t)b(x,t). Some ofthese
Forthewidthofthereactionzonetogrowmoreslowly
assumptions may be appropriate if the reaction is very
than the depletion zone caused by subdiffusion requires
rapid, but not if many encounters between reactants are
that
required for the reaction to occur.
We adopt the viewpoint put forth in a recent theory
α<γ/2. (16)
developed for geminate recombination [54, 55] but, as
On the other hand, the quasistationarity condition re- the authors themselves point out, much more broadly
quires that applicable. Thistheorygoesbacktothecontinuoustime
random walk picture from which the fractional diffusion
∂2 equation can be obtained, and considers both the mo-
K D1−γ a(x,t) ∼ t−(θ+2α+1−γ)
γ 0 t ∂x2 tion and the reaction in this framework. In the context
∂ of geminate recombination the authors define a reaction
≫ a(x,t)∼t−(θ+1), (17)
zone and argue that a geminate pair within the reaction
∂t
zone will not necessarily react for any finite intrinsic re-
which again leads to Eq. (16). action rate (which they call γ ) because one of the par-
rc
Equations (13) and (15) combined lead to the relation ticles may leave the zone before a reaction takes place.
α − µ = γ/2 − 1 that is easily checked by numerical The dynamics of leaving the reaction zone is ruled by
simulations,sinceitisdeterminedbytheproductionrate the waiting time distribution ψ (t)=ψ(t)e−γrct where
out
ofC. The rateofchangeofthe totalamountofproduct, ψ(t)isthe waitingtime thatregulatesthe restofthe dy-
dNC/dt,is givenbythe integralofthe reactionrateover namics [cf. Eq. (3)], and therefore the reaction rate will
the reaction zone, acquire a memory that arises from the same source as
the memory associated with the subdiffusive motion. In
dN
C = dxR (x,t) ∼ t−µ dxRˆ (x/tα) the continuum limit this model then leads to a reaction-
γ γ
dt reaction reaction subdiffusion equation in which both contributions have
Z zone Z zone
∼ t−(µ−α) ∼tγ/2−1, (18) a memory. Seki et al. obtain a subdiffusion-reaction
equation which at long times corresponds to choosing a
that is, reaction term of the form
NC(t)∼tγ/2. (19) Rγ(x,t)=k 0Dt1−γa(x,t)b(x,t). (22)
We stress that this total amount of product as a func-
Here“longtimes”setinveryquicklyifthereactionzone
tion of time, which is numerically more robust than its
is narrowand the intrinsic reactionrate small. As noted
derivative, is predicted to grow as tγ/2 regardless of the
earlier, although the derivation is specifically for gemi-
specific form of the reaction term.
nate recombination, the arguments can be generalized.
Another accessible quantity that is independent of R
γ
is the location x of the point at which the production Our full reaction-subdiffusion starting equations on
f
rate of C is largest. This should occur where a(x,t) ∼ which the remainder of this paper is based then are
b(x,t),thatis,u(x ,t)∼0. Thetime dependenceofthis
f
equimolar point is found from Eq. (9) to be ∂ ∂2
a(x,t) = D1−γ K a(x,t)−ka(x,t)b(x,t)
∂t 0 t γ∂x2
x (t)=K tγ/2 (20) (cid:26) (cid:27)
f f ∂ ∂2
b(x,t) = D1−γ K b(x,t)−−ka(x,t)b(x,t) .
where K is determined from the equation ∂t 0 t γ∂x2
f
(cid:26) (cid:27)
(23)
2q (1,γ) From the specific reaction term given in Eq. (22) we can
=H10 K 2 . (21)
1+q 11" f(cid:12)(0,1) # now obtain the third relation between the scaling expo-
(cid:12) nents by balancing the terms within the brackets:
(cid:12)
(cid:12)
(cid:12)
B. Choice of reaction term and resultant scaling µ=2θ+1−γ. (24)
Further relations involving the scaling exponents Simultaneous solutionof Eqs.(13), (15), and(24)finally
aimedattheirexpressionintermsofmodelquantitiesre- yields
quirespecificationofthereactionterm. Thereisavaried
literatureonthis subject,basedona numberofdifferent γ γ γ
α= , θ = , µ=1− . (25)
assumptions [52, 53, 54, 55]. Most do not associate a 6 3 3
5
C. Simulated quantities III. SIMULATION DETAILS
It is useful to list here the quantities that willbe com- MonteCarlosimulationmethods forreaction-diffusion
pared with numerical simulations. Each is character- processes are ubiquitous. For a two-dimensional simula-
ized by an exponent explicitly given in terms of γ. The tionone startswith asquarelattice, anddeploysa given
firstandsecondareindependentofthechoiceofreaction number of particles at each site according to the initial
term, but the others are sensitive to this choice. distribution. The particles then perform a random walk
simulatedbythe parallelupdate ofthe coordinatesofall
1. ThetotalamountofproductC producedasafunc-
particles at each time step t = m∆t, m = 1,2,···. The
tion of time, given in Eq. (19), is
entirelattice is exploredatperiodic intervals∆t (which
r
could and often does coincide with ∆t), and reactions
N (t)∼tγ/2. (26)
C take place at each site on which there are A and B par-
ticles, with probability kab. Here k ≪ 1 is the reaction
Thisscalingis independentoftheformofthe reac-
rateconstantandaandbareproportionaltothenumber
tion term.
of particles of type A and B on that site. Clearly, kab
2. We measure the width W ofthe depletion zone as must (in the appropriate sense, since the quantity is not
d
the width of the profile dimensionless)besmall. Therearevariantsofthisproce-
durethatareinconsequentialforouranalysis(e.g., some
U (x,t)≡1−|a(x,t)−b(x,t)|. (27)
P excluded volume effects). A necessary condition to be
in the diffusion-controlled regime described by the usual
The prediction, which is also independent of the
reaction-diffusion equations is that the random walkers
form of the reaction term, is given in Eq. (12),
onaverageperforma largenumber ofsteps before react-
W ∼tγ/2. (28) ing.
d
Adjustments that must be made to this procedure
in order to describe subdiffusion are neither trivial nor
3. We carry out our simulations with an equal initial
straightforward. Firstandmostimportantly,onecannot
unit concentrationofA and B. In this case x =0
f assumethattheparticlesalljumpatthesametime. The
for all time. We monitor the number of C parti-
distributionofjumpingtimesisnowverybroad: onecan
cles produced at this point of maximum produc-
imagineeachparticleoutfittedwithanalarmclock,with
tion of C, N (x=0,t). Since R (0,t)=dN (x=
C γ C a jump to a randomly selected nearest neighbor taking
0,t)/dt∼t−µ =tγ/3−1, we have
place when the alarm goes off, at which time the alarm
isresetaccordingtoa distributionwhoseasymptoticbe-
N (0,t)∼tγ/3. (29)
C haviorgoesasinEq.(3). Jumpingis thereforearenewal
process [56]. An example of a normalized distribution
This is thus a check on the exponent µ.
with this behavior is the Pareto law:
4. The concentration a(0,t)=b(0,t) of each reactant
γ/τ
atthecenterofthereactionzoneisdifficulttomon- ψ(t)= . (32)
(1+t/τ)1+γ
itorbecauseitisverysmallandthereforesubjectto
largefluctuations. Instead,wemonitortheintegral
Theparticlesarelabeled,andjumpingtimesareassigned
of this concentration over time,
tothemaccordingtothisdistribution. Thesetimes,from
t t smallest to largest, must be sorted, and the list must be
a(0,τ)dτ ∼ τ−γ/3dτ ∼t1−γ/3, (30) sorted after each jump or reaction.
Z0 Z0 Since the particles no longer jump at the same time, a
decision must be made about when they are allowed to
with a similar result for the other reactant. This
react. There are at least two alternatives: (1) A reac-
then is a check on the exponent θ.
tion attempt occurs only when a particle first arrives at
5. The width w(t) of the product profile grows, ac- a site; (2) Reaction attempts occur at each site at peri-
cording to the scaling equation (7), as w(t) ∼ tα. odic intervals ∆t , and occur with probability kab. The
r
According to Eq. (25) we then have,as a test of α, first alternative does not seem physically reasonable for
the subdiffusive problem since it implies that a pair of
A and B particles that remain at a given site and that
w(t)∼tγ/6. (31)
did not reactupon firstencounter will not reactno mat-
ter how long they remain at the site, which on average
6. Finally,wemonitortheentireprofile(27)asafunc- is infinite. They can only react if they move apart and
tion of position and time. This is a difficult quan- then encounter one another again. The second alterna-
titytofollowbecauseitinvolvesregionsofverylow tive,whichwechooseformostofoursimulations,canbe
concentration. Inawayitconstitutesacheckofthe associated with a number of physical explanations. On
simulation methodology, as we will see below. the one hand one can think of reactions induced or ac-
6
(cid:3)
tivated periodically by some external agent (a laser, for
example). More in line with our thinking of subdiffusion
(cid:20)(cid:19)(cid:22)
as a way to describe movement in a disordered or glassy
or porous medium is to think of this as a mesoscopic de- 1(cid:11)W(cid:12)F (cid:3)
scription. At a microscopic level small jumps may occur
diffusively,butthemotionfromonemesoscopicregionto (cid:20)(cid:19)(cid:21)
another on a longer time scale is much slower because of
geometricbottlenecksthataffectthislongerrangemove-
(cid:20)(cid:19)(cid:23) (cid:20)(cid:19)(cid:24) (cid:20)(cid:19)(cid:25)
ment. Our “sites” would then correspond to mesoscopic
W
regions in which a walker can spend a long time moving
diffusively from one part of the regionto another. Reac-
FIG.1: Log-logplotsofthetotalnumberofproductparticles
tions can then take place within one of these regions at
vs time for γ = 0.75 (squares), γ = 0.5 (circles), and γ =
regular time intervals.
0.25 (triangles). The linear fit slopes are 0.37, 0.25, and 0.15
Since the subdiffusive process has a long memory, we
respectively. The mean field prediction for the slope is given
must be careful about the initiation of the process. In
in Eq. (28) as γ/2.
particular, it is not appropriate to choose the initial
jumping times as indicated above because that would (cid:3)
(cid:20)(cid:19)(cid:21)
biasthe initialconditionto one inwhichallthe particles
jumped simultaneously at time t=0. Instead, after this
first selection of times we choose another set of jumping
times from the distribution, and repeat this procedure :G (cid:20)(cid:19)(cid:20) (cid:3)
a large number of times. The number of repetitions is
usually chosen as the total number of particles initially
in the system. Only then do we choose t = 0 by taking
thesmallestjumpingtimeasourneworiginoftimefrom (cid:20)(cid:19)(cid:19)
(cid:20)(cid:19)(cid:22) (cid:20)(cid:19)(cid:23) (cid:20)(cid:19)(cid:24) (cid:20)(cid:19)(cid:25) (cid:20)(cid:19)(cid:26)
which the process is launched. W
Finally, it is noted that even in a diffusive process,
reaction events are not really restricted to occur only at
periodic time intervals ∆t . A large literature points to FIG. 2: Width Wd of the depletion zone vs time for γ =0.5
r
(circles) and γ = 0.25 (triangles). The linear fit slopes are
thefactthatthecontinuoustimeprocessunderlyingsuch
0.257 and 0.154 respectively. The mean field prediction for
a step process is one in which times are selected from an
theslope is given in Eq. (28) as γ/2.
exponential distribution [57, 58]
P(t)=κe−κt, (33)
t =1,024,000. Results are averagedover 100 runs.
max
whereκisthereactionrateconstant. Wehavealsotested
this procedure in our subdiffusive system, allowing each
pair(A,B)ofparticlesonone sitetoreactatatime dic- IV. COMPARISONS WITH SIMULATIONS
tated by such an exponential distribution. If a particle
leaves a site before a reaction takes place, the reaction Here we compare simulations of the six quantities
“clock” of each particle is reset. This is also the view- enumerated in Sec. II with the theoretical predictions.
point followed by Seki et al. [54, 55]. Note that whereas Thesimulationsrapidlybecomeincreasinglydifficultand
inthek,∆t formulationofthe reactioneventsonespec- time-intensive with decreasing γ, and it is therefore ex-
r
ifies two parameters,in the exponentially distributed re- pected that agreement with the theory improves with
actioneventsthereactionsarespecifiedbythesinglerate increasing γ. As we will show, the agreement is on the
parameterκ,whichidentifies both the reactionratecon- wholegood,especiallyforthelargervaluesofγ. Wealso
stant κ and the mean time between reaction events 1/κ. stressthatfourofthesixcomparisonsinvolveresultsthat
The parametersusedin oursimulationsareasfollows: decidedly depend on the choice of reaction term. Agree-
we place two particles initially at each site (this corre- ment would not be obtained with the usual memoryless
sponds to unit concentration for each species). The lat- local law of mass action.
tice dimensions are usually (L ,L ) = (20,10), except Figure1showsoursimulationresultsforthetotalnum-
x y
for γ = 3/4 where we use L = 90, and in some cases ber of product particles as a function of time in units of
x
specified later where we use (L ,L ) = (160,20). The τ =1.0(usedthroughout)forseveralvaluesofγ. Thisis
x y
maximum number of particles allowed at a given site is perhapsthemostrobustglobalquantitytobesimulated.
40. The rate coefficient is k = 0.05, and the time be- The linear fit slope is given for each γ, and agrees very
tween reaction events is ∆t =10. For some of the sim- well with the theoretical prediction given in Eq. (26) for
r
ulationswe useexponentiallydistributed reactionevents the two larger values of γ.
with κ = 10−4 or κ = 10−5, which corresponds to a Figure 2 shows our simulation results for the width of
muchlowerreactionrate. The maximumtime per run is the depletion zone as a function for time for two values
7
(cid:3) (cid:3)
(cid:20)(cid:19)(cid:20) (cid:20)(cid:19)(cid:20)
1(cid:11)[ (cid:19)(cid:15)W(cid:12)F (cid:3) Z(cid:11)W(cid:12) (cid:3)
(cid:20)(cid:19)(cid:19)
(cid:20)(cid:19)(cid:19)
(cid:20)(cid:19)(cid:22) (cid:20)(cid:19)(cid:23) (cid:20)(cid:19)(cid:24) (cid:20)(cid:19)(cid:25) (cid:20)(cid:19)(cid:22) (cid:20)(cid:19)(cid:23) (cid:20)(cid:19)(cid:24) (cid:20)(cid:19)(cid:25)
W W
FIG. 3: Log-log plots of the production profile of product C FIG. 5: Log-log plot of the width of the product profile as a
at x=0as afunction oftime for γ =0.75 (squares), γ =0.5 function of timeforγ =0.75 (squares), γ =0.5 (circles), and
(circles), and γ = 0.25 (triangles). The linear fit slopes are γ =0.25(triangles). Thelinearfitslopesare0.129,0.084,and
0.24,0.162,and0.093respectively. Themeanfieldprediction 0.042 respectively. The mean field prediction for the slope is
for the slope is given in Eq.(29) as γ/3. given in Eq. (31) as γ/6.
(cid:3)
(cid:20)(cid:17)(cid:19)
(cid:19)(cid:17)(cid:27)
(cid:3)
(cid:19)(cid:17)(cid:25)
RI(cid:3)D(cid:11)(cid:19)(cid:15)W(cid:12)(cid:15)(cid:3)E(cid:11)(cid:19)(cid:15)W(cid:12)(cid:20)(cid:20)(cid:19)(cid:19)(cid:23)(cid:24) 8(cid:11)[(cid:15)W(cid:12)3(cid:19)(cid:19)(cid:17)(cid:17)(cid:21)(cid:23) (cid:3)
JUDO(cid:3) (cid:3) (cid:19)(cid:17)(cid:16)(cid:19)(cid:20)(cid:19)(cid:19) (cid:16)(cid:24)(cid:19) (cid:19) (cid:24)(cid:19) (cid:20)(cid:19)(cid:19)
PH(cid:3)LQWH(cid:20)(cid:19)(cid:22) [
7L
(cid:20)(cid:19)(cid:21)
(cid:20)(cid:19)(cid:24) (cid:20)(cid:19)(cid:25)W (cid:20)(cid:19)(cid:26) FIG.6: SimulationresultsfortheprofileUP(x,t)forγ =0.75
and t = 946, 3777, 15073, 60149, 240025, and 957828. The
width of the profiles increase with time, as seen in Fig. 5.
The simulation was carried out on a lattice of size (160,20)
FIG.4: Time integrals ofreactant concentrations at thecen-
and averaged over 102 runs. Notice the apparent evolution
ter of the reaction zone. The mean field prediction for the
slopesisgiveninEq.(30)as1−γ/3. Thesteepercurveisfor from a characteristic sharp-pointed profile for short times to
a vaulted profile at longer times. For a discussion of this
γ = 0.25 and the reaction events governed by the exponen-
tial distribution Eq (33) with κ= 10−5. The linear fit slope anomaly, and for thevalues of other parameters, see text.
is 0.814 while mean field theory yields 0.917. The shallower
curve is for γ = 0.5 and κ = 10−4. The linear fit slope is
0.683, themean field slope is 0.833. sion of the profile U (x,t).
P
of γ. The linear fit slope is in good agreement with the
theoryasgiveninEq.(28)forthelargervalueofγ. Later Figure 5 contains our simulation results for the width
we discuss some difficulties, particularly for small values of the product profile, which should be compared with
of γ, in the accurate simulation of the profile UP(x,t) the prediction of Eq. (31). The agreement is very good
whose width is used to obtain these results. for all the values of γ.
Figure 3 shows our simulation results for the produc- Finally, in Figs. 6 and 7 we present perhaps the most
tion profile of the product of the reaction at x = 0 as a difficult quantity to capture accurately, namely, the pro-
function of time, for severalvalues of γ. For the larger γ fileU (x,t) definedinEq.(27). Forthesesimulationswe
p
the linear fit slope agrees well with the theoretical pre- use a lattice of size (L ,L )=(160,20). As pointed out
x y
diction given in Eq. (29). earlier, the difficulty arises from the fact that it involves
Figure 4 presents our simulation results for the time- regions of very low concentration. It is instructive to il-
integral of the concentration of a reactant at the center lustrate the difficulty, and that is why we have included
of the reaction zone, to be compared with the theoret- Fig. 6. The simulation profiles shown at the different
ical prediction of Eq. (30). While the agreement is not times are actually time averagedovera small time inter-
spectacular, the trend is correct. Also of interest here is val around the times shown. The noteworthy feature is
theimprovedagreementwhenthereactionrateisgreatly the evolutionof the sharp-pointedprofile near the origin
reduced, as expected. Even more dramatic effects of the at short times to a more rounded shape at longer times.
reactionrate areseenbelow inthe contextof ourdiscus- The mean field theory presented in this paper does not
8
(cid:3)
(cid:20)(cid:17)(cid:19) theevolutionofareactionfront. Becauseweareworking
withasetofmeanfieldcontinuumequations,ourresults
(cid:19)(cid:17)(cid:27)
are only valid above the critical dimension d=2.
The subdiffusive motionis modeledvia the usualfrac-
(cid:19)(cid:17)(cid:25)
(cid:11)[(cid:15)W(cid:12)3 (cid:19)(cid:17)(cid:23) (cid:3) teiroantaolr,eqEuqa.t(io5n).thTahtiscocnhtoaiicneshtahseaRileomnganhnis-tLoiroyu,vailnledoipts-
8
virtues and shortcomings are clearly understood. Less
(cid:19)(cid:17)(cid:21) clear has been the selection of the local reaction term,
and the question of the way in which the memory in
(cid:19)(cid:17)(cid:19) the Riemann-Liouville operator does (or does not) af-
(cid:16)(cid:27)(cid:19) (cid:16)(cid:25)(cid:19) (cid:16)(cid:23)(cid:19) (cid:16)(cid:21)(cid:19) (cid:19) (cid:21)(cid:19) (cid:23)(cid:19) (cid:25)(cid:19) (cid:27)(cid:19)
fect the way in which the reaction is modeled. While
[
the literature on this subject has presented a number of
viewpoints, we argued, in agreement with [54, 55], that
the reaction term should also be modified from its usual
FIG. 7: Simulation (jagged lines) and theory (dashed lines)
fortheprofileUP(x,t)forγ =0.75andt=3738,8577,19682, simple instantaneous product form, at least for small re-
45165,103638,237815,and545708. Noticetheabsenceofthe action rate constants. Our reaction-subdiffusion model
rounding anomalies in the profile for small x. The reaction is thus given by Eq. (23).
rateismuchsmallerinthiscase,seetext. Thesimulationwas
Following the approach of G´alfi and Ra´cz [29] for the
carriedoutonalatticeofsize(160,20)andaveragedover102
evolutionofa frontin the reaction-diffusionproblem, we
runs. The times are approximately equivalent to those used
assumed scaling solutions for the various quantities that
in Fig. 6 if the total number of C particles in the system is
canbe calculatedfromthe model. Some of these quanti-
used as a measure of time.
tiesdependexplicitlyontheformchosenforthereaction
termwhileothersdonot. Wecomparedtheresultingex-
producethis rounding,sothatitseemedatfirstthatthe ponentswiththoseobtainedfromnumericalsimulations.
theory and simulations differed in some profound way. We found very good agreement between the theory and
However,the simulations in Fig. 6 were carriedout with simulations for the exponents µ and α that characterize
the rate coefficientk =0.05 with reactions occurringpe- the reaction term, Eq. (7). In particular, in terms of the
riodically at time interval ∆t =10, a reaction rate that power γ < 1 that characterizes the subdiffusive process
r
turns outto be too high for comparisonwith our theory. we found that the theoretical values µ = 1− γ/3 and
In Fig. 7 we show both the simulation results (jagged α = γ/6 are recovered in the simulations with greater
curves) and those of our theory (dashed curves), now fidelity for larger γ. The exponent θ, governingthe time
with exponentially distributed reaction events according decay of the reactant concentrations as in Eq. (6), is,
to Eq. (33), with κ = 10−5. The profile near x = 0 theoretically, given by γ/3. Simulation results give val-
remains pointed for all the times shown, as predicted ueswhichcorrectlyfollowthistrend,buttheagreementis
by the theory. The quantitative disagreements for long notquantitative. However,wehavetorememberthatthe
times (t = 237185 and especially t = 545708) are due to quantity a(0,t) is local and, consequently, it is more dif-
finite size effects. Boundary effects are negligible only as ficult to achieve good statistical averages with the small
long as U (x = −L /2,t) = U (x = L /2,t) = 0, and systems and number of particles considered in the sim-
p x P x
for the long times we find that this is not the case. In ulations. Perhaps the most challenging quantity to cap-
our other simulation results we have not included such ture is the profile UP(x,t). The theory predicts a cusp
results in our averages, but have left them in this fig- at x = 0 which we were able to capture by our simula-
ure simply to stress the finite size effects one must be tionswhenthereactionrateconstantissufficientlysmall.
aware of in calculating the behavior of quantities in the The quantitative agreement between the theory and the
depletion zone when the zone extends all the way to the simulations for this profile was ultimately limited by our
boundaries of the system. To provide accurate results finitesystemsize. Wenotethatourresultsareagoodex-
for such long times it is necessary to run simulations on ample of what is sometimes referred to as subordination
larger systems. in that the subdiffusive scaling behavior can be deduced
from the corresponding diffusive behavior with the sub-
stitution t→tγ [59].
V. CLOSING COMMENTS This work can clearly be pursued along a number of
directions. Among them is the description of this same
In this paper we have proposed a set of continuum reaction-subdiffusionproblemwiththeusualuniformini-
fractionaldiffusion equations to describe the behavior of tial condition for the species, to investigate what sorts
a reaction front in the A+B → C reaction-subdiffusion of segregation patterns might evolve on the way to ex-
problem. Subdiffusion may be appropriate to describe tinction, or on the way to equilibration if the reaction
the way reactants move in complex (glassy, disordered, is reversible. Another is the study of the fluctuations
highlyconstrained)geometries,andwewereinterestedin that must be added to the model in order to describe
exploringhowthis constraintonthe motionwouldaffect thisprocessinaone-dimensionalsystemwherethemean
9
field description is no longer appropriate, and the pos- BFM2001-0718 and by the Engineering Research Pro-
sible logarithmic correction in two dimensions that may gram of the Office of Basic Energy Sciences at the U.
explainsome of our small-γ discrepancies. A third is the S. Department of Energy under Grant No. DE-FG03-
effect of different subdiffusion coefficients for the species 86ER13606.
A and B, and even of different exponents γ and γ . In
A B
this latter case subordination would necessarily be more
complicated if valid at all. Work along these directions
is in progress [60].
Acknowledgments
This work was partially supported by the Ministe-
rio de Ciencia y Tecnolog´ıa (Spain) through Grant No.
[1] For a review see E. Kotomin and V. Kuzovkov, Modern (1991).
Aspects of Diffusion-Controlled Reactions: Cooperative [21] F. Leyvraz and S. Redner, Phys. Rev. Lett. 66, 2168
Phenomena in Bimolecular Processes (Elsevier, Amster- (1991).
dam, 1996). [22] G. Zumofen, J. Klafter, and A. Blumen, Phys. Rev. A
[2] D.C.TorneyandH.McConnell,J.Phys.Chem.87,1941 45, 8977 (1992).
(1983). [23] P. Argyrakis, R. Kopelman, and K. Lindenberg, Chem.
[3] J. L. Spouge, Phys. Rev.Lett. 60, 871 (1988). Phys.177,693(1993); K.Lindenberg,P.Argyrakis,and
[4] C. R. Doering and D. ben-Avraham, Phys. Rev. A 38, R. Kopelman, J. Phys.Chem. 98, 3389 (1994).
3035 (1988); M. A. Burschka, C. R. Doering and D. [24] H. Taitelbaum, Y-E. L. Koo, S. Havlin, R. Kopelman,
ben-Avraham, Phys. Rev. Lett. 63, 700 (1989); D. ben- andG.H.Weiss,Phys.Rev.A46,2151(1992);H.Taitel-
Avraham, M. A. Burschka and C. R. Doering, J. Stat. baum,A.Yen,R.Kopelman,S.Havlin,andG.H.Weiss,
Phys.60,695(1990);D.ben-Avraham,Mod.Phys.Lett. Phys. Rev. E 54, 5942 (1996); H. Taitelbaum, B. Vilen-
9, 895 (1995). sky, A. Lin, A. Yen, Y-E. Lee Koo, and R. Kopelman,
[5] K.Lindenberg, P. Argyrakisand R. Kopelman, J. Phys. Phys. Rev.Lett. 77, 1640 (1996)
Chem. 99, 7542 (1995). [25] S. H. Park, S. Parus, R. Kopelman, and H. Taitelbaum,
[6] K. Krebs, M. P. Pfannmuller, B. Wehefritz and H. Hin- Phys. Rev.E 64, 055102(R) (2001).
richsen, J. Stat. Phys78, 1429 (1995). [26] C. N.Baroud, F. Okkels, L.M´en´etrier, and P. Tabeling,
[7] T. O. Masser and D. ben-Avraham, Phys. Lett. A 275, Phys. Rev.E 67, 060104(R) (2003).
382 (2000); Phys. Rev. E 63, 066108 (2001); Phys. Rev. [27] R. Kopelman, Z-Y. Shi, and C-S. Li, J. Lumin. 48/49,
E 64, 062101 (2001). 143 (1991) and references therein.
[8] S.Habib,K.Lindenberg,G.Lytheand C.Molina-Par´ıs, [28] S.HabibandG.Lythe,Phys.Rev.Lett.84,1070(2000).
J. Chem. Phys. 115, 73 (2001). [29] L. G´alfi and Z. R´acz, Phys. Rev A (Rapid Comm.) 38,
[9] C.MandacheandD.ben-Avraham,J.Chem.Phys.112, 315 (1988).
7735 (2000). [30] Z. Jiang and C. Ebner, Phys.Rev.A 42, 7483 (1990).
[10] M.HenkelandH.Hinrichsen,J.Phys.A34,1561(2001). [31] S. J. Cornell, M. Droz, and B. Chopard, Phys. Rev. A
[11] E. Abad, H. L. Frisch and G. Nicolis, J. Stat. Phys. 99, 44, 4826 (1991).
1397 (2000). [32] S. J. Cornell, Phys. Rev. Lett. 75, 2250 (1995); S. J.
[12] A.A.OvchinnikovandY.G.Zeldovich,Chem.Phys28, Cornell, Phys.Rev.E 51, 4055 (1995).
215 (1978). [33] E.ben-NaimandS.Redner,J.Phys.A25,L575(1992).
[13] D. Toussaint and F. Wilczek, J. Chem. Phys. 78, 2642 [34] H. Larralde, M. Araujo, S. Havlin, and H. E. Stanley,
(1983). Phys. Rev A 46, 855 (1992); H. Larralde, M Araujo,
[14] K.KangandS.Redner,Phys.Rev.Lett.52,955(1984). S. Havlin, and H. E. Stanley, Phys. Rev. A 46, R6121
[15] R. Kopelman, Science 241, 1620 (1988) and references (1992); M. Araujo, S. Havlin, H. Larralde, and H. E.
therein. Stanley,Phys.Rev.Lett.68,1791(1992);M.Araujo,H.
[16] S.Kanno,Prog. Theor.Phys.79,721(1988); ibid,1330. Larralde, S. Havlin, and H. E. Stanley, Phys. Rev. Lett.
[17] E.Cl´ement,L.M.Sander,andR.Kopelman,Phys.Rev. 71, 3592 (1993); M. Araujo, H. Larralde, S. Havlin, and
A 39, 6455 (1989); ibid,6466; ibid,6472. H. E. Stanley,Phys. Rev.Lett. 75, 2251 (1995).
[18] P. Argyrakis and R. Kopelman, Phys. Rev. A 41, 2114 [35] B.Chopard,M.Droz,T.Karapiperis,andZ.Racz,Phys.
(1990); ibid,2121. Rev.E47,pR40(1993);B.Chopard,D.Droz,J.Magnin,
[19] K.Lindenberg,B.J.West,andR.Kopelman,Phys.Rev. and Z. Racz, Phys.Rev.E 56 5343 (1997).
Lett.60,1777(1988); ibid,Phys.Rev.A42,890(1990); [36] B. P.Lee and J. Cardy, Phys.Rev.E 50, R3287 (1994);
W-S.Sheu,K.Lindenberg,andR.Kopelman,Phys.Rev. G. T. Barkema, M. J. Howard, and J. L. Cardy, Phys.
A 42, 2279 (1990). Rev. E53, R2017 (1996).
[20] M. Bramson and J. L. Lebowitz, J. Stat. Phys. 62, 297 [37] P. L. Krapivsky,Phys.Rev.E 51, 4774 (1995).
10
[38] Z. Koza, J. Stat. Phys. 85, 179 (1996); Z. Koza and H. [48] R.Metzler,E.BarkaiandJ.Klafter,Europhys.Lett.46,
Taitelbaum,Phys.Rev.E54,R1040(1996);Z.Kozaand 431 (1999); E. Barkai, R. Metzler, and J. Klafter, Phys.
H. Taitelbaum, Phys. Rev. E 56, 6387 (1997); Physica Rev. E61, 132 (2000).
A 240, 622 (1997); Philos. Mag. B 77, 1437 (1998); H. [49] W. R. Schneider and W. Wyss, J. Math. Phys. 30, 134
Taitelbaum and Z. Koza, Physica A 285, 166 (2000); (1989).
Phys.Rev.E 66, 011103 (2002); Europhys.J. B32,507 [50] A. M. Mathai and R. K. Saxena, The H-function with
(2003). Applications in Statistics and Other Disciplines, (John
[39] See e.g. R. Metzler and J. Klafter, Phys. Rep. 339, 1 Wiley and Sons, NewYork,1978).
(2000);B.J.West,M.Bologna,andP.Grigolini,Physics [51] D. ben-Avraham,Mod. Phys. Lett. 9, 895 (1995).
of Fractal Operators (Springer, NewYork,2003). [52] B.I.HenryandS.L.Wearne,PhysicaA276,448(2000);
[40] D.ben-AvrahamandS.Havlin,Diffusion and Reactions B. I. Henry and S. L. Wearne, SIAM (Soc. Ind. Appl.
in Fractals and Disordered Systems (Cambridge Univer- Math.) J. Appl. Math. 62, 870 (2002); M. O. Vlad and
sity Press, Cambridge, 2000); J.Ross,Phys.Rev.E66,061908 (2002);S.Fedotovand
[41] T. Kostol owicz, K. Dworecki, and St. Mr´owczyn´ski, V. M´endez, Phys. Rev.E 66, 030102 (20020).
cond-mat/0309072. [53] J. Sung, E. Barkai, R. J. Silbey, and S. Lee, J. Chem.
[42] S. B. Yuste and K. Lindenberg, Phys. Rev. Lett. 87, Phys. 116, 2338 (2002).
118301 (2001); ibid,Chem. Phys.284, 169 (2002). [54] K.Seki,M.Wojcik,andM.Tachiya,J.Chem.Phys.119,
[43] J-C. Lin, C. R. Doering and D. ben-Avraham, Chem. 2165 (2003).
Phys. 146, 355 (1990); J-C. Lin, Phys. Rev. A 44, 6706 [55] K.Seki,M.Wojcik,andM.Tachiya,J.Chem.Phys.119,
(1991). 7525 (2003).
[44] ForreviewsseeJ.W.HausandK.Kehr,Phys.Rev.150, [56] D. R. Cox and H. D. Miller, The Theory of Stochastic
264 (1987); J.-P. Bouchaud and A. Georges, Phys. Rep. Processes (Wiley, NewYork,1965).
195, 127 (1990). [57] D. Bedeaux, K. Lakatos-Lindenberg, and K. E. Shuler,
[45] E. Barkai, Chem. Phys. 284, 13 (2002). J. Math. Phys. 12, 2116 (1971).
[46] I.M.Sokolov,A.BlumenandJ.Klafter,Europhys.Lett. [58] V. M. Kenkre, E. W. Montroll, and M. F. Shlesinger, J.
56, 175 (2001). Stat. Phys. 9, 45 (1973).
[47] R. Hilfer, Editor, Applications of Fractional Calculus in [59] A.Blumen,J.Klafter,andG.Zumofen,inOptical Spec-
Physics, (World Scientific, Singapore, 2000); V. Balakr- troscopy of Glasses, ed. I. Zschokke (Reidel, Dordrecht,
ishnan, Physica A 132, 569 (1985); C. R. Doering and 1986); A. Blumen, J. Klafter, B. S. White, and G. Zu-
D. ben-Avraham, Phys. Rev. Lett. 62, 2563 (1989); R. mofen, Phys. Rev. Lett. 53, 1301 (1984); J. Klafter, A.
Metzler, E. Barkai and J. Klafter, Phys. Rev. Lett. 82, Blumen, andG. Zumofen,J. Stat.Phys. 36, 561(1984).
3563 (1999); [60] S. B. Yusteet al., in preparation.