Physica A 205 (1994) 487-496 North-Holland
SSDI 0378-4371(94)00014-K
The velocity autocorrelation function for a random walk on percolation clusters Czestaw Oleksy Institute of Theoretical Physics, University of Wroctaw, Pl. Maksa Borna 9, 50-204 Wroctaw, Poland Received 15 September 1993 Revised manuscript received 22 December 1993 We apply the random Lorentz model to diffusion on percolation clusters on a square lattice. In this model a walker performs a random walk on duster sites and can visit sites occupied by impurities which play the role of reflectors. The velocity autocorrelation function and the mean-square displacement of a random walker are calculated with high accuracy using a moment-propagation numerical technique. Exponents describing the decay of correlation function are calculated at the percolation threshold and above it.
1. Introduction
Diffusion in disordered systems is often described by a model of random walk on disordered structure [ 1,2 ]. Diffusion on critical percolation clusters is anomalous, i.e. the mean square displacement of the walker is not a linear function of time but obeys a power-law : R E ( t ) ,,~ t 2/dW with dw > 2. The exponent dw has been found numerically to high accuracy in an exact enumeration approach [ 1,3 ]. Jacobs and Nakanishi [4] studied the behaviour of the velocity autocorrelation function (VAF) for site percolation problem on 2D and 3D lattices. They calculated the VAF for the model o f " a myopic ant in a labyrinth" as well as for "the blind ant". The VAF was obtained from the second derivative of the mean square displacement calculated by the exact enumeration technique. They showed the anomalous behaviour of the VAF at the percolation threshold and for the model of the myopic ant on a bipartite lattice (square, cubic, etc. ) they found even-odd oscillations in the VAF decaying much slower than the diffusive term. In this paper we are going to calculate the VAF directly and with high precision in order to check the long time behaviour of the VAF. We study this problem in site-percolation on the square lattice. The paper is organized as 0378-4371/94/$ 07.00 (~) 1994-Elsevier Science B.V. All rights reserved
C. Oleksy / VAF for random walk on percolation clusters
488
follows. In section 2 we present the random Lorentz model, which we apply to diffusion on percolation clusters. The method of calculation of the VAF using moment propagation technique [5-7] is presented in section 3. Numerical results for the VAF and the mean square displacement are discussed in section 4. The results of calculation of the diffusion coefficient above the percolation theshold are presented in section 5.
2. The model In a previous paper [7] we introduced the random Lorentz model (RLM) and we showed that this model can be applied to diffusion in site-percolation systems. The model consists of the walker moving on the lattice with a fraction of sites, p, occupied by scatterers (impurities) and a fraction, 1 - p, occupied by particles. The walker jumps once each unit time step from a site to one of neighbouring sites. There are two kinds of jumps: (i) random jump from a site occupied by a particle; (ii) correlated jump, dependent on the previous jump, from a site occupied by a scatterer. In other words, in RLM the walker performs a random walk on sites occupied by particles and correlated random walk on sites occupied by scatterers. When scatterers become reflectors then the model describes diffusion on a percolation cluster, because if the walker jumps from a cluster site to a site occupied by an impurity then it will be reflected backward to the cluster site in the next time step. In this paper we will consider percolation clusters on a square lattice. It is worth noting that our model differs slightly from the most popular model of the random walk of the myopic (or blind) ant, because in our case the walker moves on cluster sites and can also visit sites with reflectors which are nearest neighbours of cluster sites. However, the model presented here is more convenient for numerical calculation.
3. Velocity autocorrelation function In the random walk described above the speed of the particle (walker) v is constant and is assumed to be v = 1. Hence the velocity of the particle v is determined by four lattice vectors: ek = (cos(½zck),sin(½nk)),
k = 0,1,2,3.
The velocity autocorrelation function (VAF) is defined as
(3.1)
C. Oleksy / VAF for random walk on percolation clusters
cI) ( N ) = (Vx (0) Vx (N)),
489
(3.2)
where Vx ( N ) is the x-component of the particle velocity after N time steps and (...) means average over different trajectories and different clusters. Applying the moment propagation (MP) numerical technique [5,7] we calculate exactly the VAF for a given cluster. The basic idea of the MP method is following: instead of tracing the individual particle we propagate the site probability P(r, i, N ) , which means the probability that the particle after N steps is at site r with velocity directed along the lattice vector ei. In order to measure the VAF directly, the site probability is weighted with initial velocity Vx (0) and the "weighted" site probability, W (r, i, N), satisfies the following equation:
W(r,i,N + 1) = E Tij(r) W ( r - e j , j , N )
(3.3)
J with initial value W (r, i, O) = P(r, i, O) Vix. The transition matrix T ( r ) depends on the configuration of the refectors
Tij( ) = ¼, = 3 ( v i + v j),
for a cluster site at r, for a reflector at r,
(3.4)
where A (v) is non-zero ( equal to 1 ) only for v --- 0. Choosing the steady state as initial one for a given cluster, the VAF is calculated exactly from the following formula:
• c(N) = ~_,vixW(r,i,N),
(3.5)
r,i
because all possible trajectories are accounted for. In order to get the VAF for a given density of reflectors we should average q0c over different clusters. It is easy to check that in the steady state, all site probabilities P (r, i ) have the same value, P0, which depends on the cluster topology. Let ns denote the number of sites in a given cluster and nrx, nry the number of links from reflectors to cluster sites along the direction x, y, respectively. Hence the number of all posible initial velocities for a given cluster can be expressed as n = 4ns + nrx + n~y, and this yields P0 = 1In. Using the above notion we can write the values of VAF for two initial time steps: q~c(0) = (2ns + nrx)/n, • c(1) = - n r x / n and ~c(2) = 0. The results for N = 1 and N = 2 were obtained by taking into account that there is no correlation on cluster sites and
490
C. Oleksy / VAF for random walk on percolation clusters
that the reflector changes the sign of the velocity. It is worth noting that the contribution to the VAF at all time steps comes only from the sites occupied by reflectors, or in the other words the VAF depends on the distnoution of reflectors. Having a sequence ~c (0), ~c ( 1 ) . . . . . tgc ( N - 1 ) calculated for a given cluster exactly, we can also get the exact value of the mean-square displacement along the x-direction for a given cluster from the following formula:
N-I XZ(N) = N
~c(O) + 2 Z
) ~c(S)
S=I
N-1 - 2Z
S~c(S).
(3.6)
S=I
4. Numerical results Now we present the details of numerical calculation wich were performed on PC computer with 16MB RAM. To calculate the VAF each cluster was generated on a square lattice of 457 x 457 sites by choosing randomly a fraction of sites, p, to be occupied by reflectors and then finding the sites belonging to the percolation cluster. The positions of the cluster sites and reflectors were stored in an array. Two additional arrays were used to update the weighted site probabilities according to Eq. (3.3). The first step of the MP method was to find the steady state. For each time step up to 4000, the value of ~c was calculated and stored. The VAF was averaged over 100 clusters at the percolation threshold (p = 0.407) and over 20 clusters above the threshold (p < 0.407). As in the case of a myopic ant on a percolation cluster [4] we observed an even-odd oscillation of the VAF (Fig. 1 ). The decay of the VAF cannot be described by a simple power law because the exponents of the upper and lower envelopes of the VAF are different. We used the following formula:
O ( N ) ,-~ A ( - 1 ) N N -a - B N -b,
(4.1)
to fit the numerical data. Exponents and amplitudes were calculated using the Levenberg- Marquardt method of nonlinear fitting [8]. To calculate the exponent at time step N we used 101 successive values of the VAF: ~ ( N 50) . . . . . • (N + 50). As is seen in Fig. 2 the exponents a, b change slowly with N. From extrapolation of a (N), b (N) for N ~ oo we obtained the following results at the percolation threshold a = 0.64 + 0.01 and b = 1.292 + 0.01. Our result is different from the result for the myopic ant on a square lattice [4] where the authors found that the upper and lower envelopes of the VAF decay
C. Oleksy / VAF for random walk on percolation clusters
15
......................................
491
15
-10 -151 . 0
.
.
.
.
.
.
.
.
.
.
.
.
1000 ,
1 0
.
,
,
.
.
.
.
.
.
3000
,,r,,i
10 -1
(b)
-1.
o
10 -z-
o
10
D oa o
o
10
--15 4000
.
2000
o a o a
-2
n
10 -3
-a.
10 -'-
10 -* +
10
-6_
I0
-6
10 -6
,
,
,
, ....
i
10
'
'
' ' ' ' " ,
10 2
'
'
' ' ' ' " 1
10
'
'
10
-6
~
N
Fig. 1. (a) The VAF vs the number of time steps at the percolation threshold. The upper curve is for even time steps and the lower curve for odd time steps. (b) Log-log plot of the absolute value of the VAF. The lower curve (crosses) indicates error bars.
The first, dominating, term in formula (4.1) has no physical meaning and is generated by the geometry of the lattice (bipartite) and dynamics. The second term in (4.1) has a diffusive nature and the exponent b is related to the exponent dw describing the anomalous behaviour of the root-mean-square d i s p l a c e m e n t , R d* ( N ) ,,~ N , at the percolation threshold: dw ~ 2 / ( 2 - b ) for large N. In order to obtain the exponent dw we calculated < X 2 (N) > using formula (3.6) and fitted the result with a power-law. The dw as a function of N is presented in Fig. 3 where we also plotted the expression 2 / [ 2 - b ( N ) ] . The small difference between these quantities arose because we used only the two leading terms in the approximation the the VAF by expression (4.1). The third term would be oscillatory and could affect a little the exponent b for large N but would not change the < X 2 (N) >. From extrapolation we found dw = 2.873 + 0.012, which is nearly the same as obtained by Majid et al. [3]: dw = 2.86 + 0.02 in the exact enumeration approach. We calculated the VAF above the percolation threshold for the following values of the density of reflectors: p = 0.05, 0.1,0.15, 0.2, 0.25, 0.3, 0.35. Far a s N -0-69.
492
C. Oleksy / VAF for random walk on percolation clusters
..........
1.5
1.5 1.4
1.4 ~. . . . . . .
'
1.3
'~- 1.3
1.2
1.2
(-31.1
i.i
~1.0
1.0
©o.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
, , , , ,
. . . .
i ~ - r - ~ ,
i000
,
,
,
,
i
. . . .
,
. . . .
2000
i , , , , , , , , ,
3000
0 . 5
4000
N Fig. 2. Plot of the exponents a (the lower curve) and b (the u p p e r curve) vs N. The error bars are drawn for several points.
3.00
T 3-3.00 .00
......................
2.95
~2.95
290
......
2.85
2.85
2.80
/
/
2.80
tt
2.75 2.70
~ . . . .
2.75 ,~ .................. 1000 2000
T........ 2.70 3000 4000
N Fig. 3. Comparing the exponent dw (square) with the expression 2 / ( 2 - b) (circle).
above the threshold p < 0.15 we calculated the VAF for 2000 time steps and for remaining densities N was up to 4000. The even-odd oscillations occur in the VAF's (see Fig. 4) and we see that the correlations become weaker far above the threshold. Approximation of the data by formula (4.1) yields the exponents a and b above the threshold (Fig. 5). Omitting the results for
C Oleksy / VAFfor random walk on percolation clusters
10
-1.
10
-c
l O -1
I
1 0 -3
1 0 -~
10
493
-7
~,.ln . . . . . . . . . . . .1. .0. 2
1
10
-6
10
-7
i5
N Fig. 4. The absolute value of the VAF vs N above the percolation threshold. The triangle, square and cross represent data for the density of reflectors 0.05, 0.2 and 0.35, respectively. Note the strong even-odd oscillations of the VAF.
2.0
2.0
m $
..Q 1.5
(3
1.5
1.0
1.0
0.5
0.5
0.0
........
0.0
, .........
0.1
, ...................
0.2
0.3
,,,,
0.0
0.4
P Fig. 5. The exponents a (circle) and b (square) vs the density of reflectors above the percolation threshold. The full symbols represent data at the threshold. Near the threshold (p > 0.3) the value N = 4000 is not enough to estimate the exponents correctly.
C. Oleksy / VAF for random walk on percolation clusters
494
p > 0.30, where the exponents do not saturate after 4000 time steps, we found that a is close to 1 and b is close to 2. Hence the value of b is the same as in the 2-D Lorentz gas [6]. Our results at the percolation threshold and above it lead to the conlcusion that the exponents a, b seem to be related via the following formula b = 2a,
(4.2)
which is fulfilled in our case within statistical accuracy.
5. The diffusion coefficient The mean-square-displacement, above the percolation threshold, is linear in the number of time steps (for a large number of steps). Hence the time dependent diffusion coefficient can be defined as D(N) -
(X2(N)) 2N
(5.1)
Using the data for the VAF above the threshold we calculated < X 2 (N) > and then the D (N) for several values of the density of reflectors. The diffusion coefficient D = lims--.~ D (N) was calculated using the extrapolation method. We observed that for p < 0.2 the value of D was nearly the same as D (2000) whereas close to the percolation threshold the extrapolated value was significantly reduced (for example, the reduction for p = 0.35 is about 21%). The results for D presented in Fig. 6 were very well fitted by the following formula: D ( p ) = 0.25(1 - 3.107p + 1.584p2).
(5.2)
Taking into account that D(0) = 1/4 we notice that our result is very close to the effective medium result of Watson and Leath [9] for the conductivity in site-percolation on a square lattice,
a ( p ) = cr(O) (I - ~ r p + ½~rp2).
(5.3)
There is only one difference between our model and walk of myopic ant because in the latter case the diffusion coefficient is related to the conductivity in the following way a ( p ) / a ( O ) = (1 - p ) D ( p ) / D ( O ) . This difference can easily be explained. In our model the walker is moving over cluster sites and sites occupied by reflectors. Hence, above the threshold, all sites are practically accessible for the walker.
C. Oleksy / VAF for random walk on percolation clusters
495
0.30
0.25 0.20 D 0.15 0.10 0.05 0.00
0.0
0.i
0.2
0.3
......... 0.4 0.5
P Fig. 6. The diffusion coefficient vs the density of reflectors above the percolation threshold. The symbols denote the extrapolated values and the line is a parabolic fit.
6. Discussion We have presented a very accurate method of direct calculation of the VAF in site-percolation system on a square lattice. For a single cluster the method gives the exact value of the VAF. The statistical errors are due to the averaging of this quantity over different clusters. The method can be applied also for other models - the random walk of a myopic ant, blind ant- as well as to other lattices. The VAF was calculated at the percolation threshold and above it and we found that (i) the decay of the VAF on a square lattice for large number of steps can be described by the sum A ( - 1 ) N N - a -- BN-b; (ii) at the percolation threshold the VAF exhibits anomalous behaviour the correlations become much stronger than they are above the threshold. The even-odd oscillations in the VAF are due to the geometry of the lattice and the dynamics and were also observed in other discrete kinetic models, i.e. in the lattice Lorentz gas on a square lattice [6 ]. However, in site-percolation model on a square lattice these oscillations dominate the long-time behaviour of the VAF. The method allows us obtain to within high accuracy the mean-square displacement using calculated values of the VAF at successive time steps. We found the exponent dw, which is the same as in the exact enumeration technique [3 ]. We hope that the method applied to a model of a random walk
496
C. Oleksy / VAF for random walk on percolation clusters
on percolation clusters, in which the VAF does not possesses even-odd oscillations, will allow us to find logaritmic correction to the long-time tail of the VAF [10]. Finally we want to add that the modification of RLM to biased diffusion on percolation clusters allows us to calculate the drift velocity, trapping probabilities [ 11 ] as well as the VAF and the mean square displacement.
Acknowledgements This work was supported by KBN Grant No 2 0126 91 01.
References [1] S. Havlin and D. Ben-Avraham, Adv. Phys. 36 (1987) 695. [2] D. Stauffer and A. Aharony, Introduction to Percolation Theory ( Taylor and Francis, London, 1992). [3] I. Majid, D. Ben-Avraham, S. Havlin and H.E. Stanley, Phys. Rev. B 30 (1984) 1626. [4] D. Jacobs and H. Nakanishi, Phys. Rev. A 41 (1990) 706. [5] M.A.van der Hoef and D. Frenkel, Phys. Rev. A 41 (1990) 4277. [6] P.M. Binder and D. Frenkel, Phys. Rev. A 42 (1990) 2463. [7] J. Kortus and Cz. Oleksy, J. Phys. A 25 (1992) 1093. [8] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipies in Fortran (Cambridge University Press, Cambridge, 1992). [9] B.P. Watson and P.J. Leath, Phys. Rev. B 9 (1974) 4893. [10] Th.M. Nieuwenhuizen, P.F.J. van Velthoven and M.H. Ernst, Phys. Rev. Lett. 57 (1986) 2477. [1 l] C. Oleksy, The drift velocity in percolating systems, preprint IFT UWr 853 (1993).