text:
# Survival probability and order statistics of diffusion on disordered media

## I INTRODUCTION

Diffusion in disordered media has been an area of intensive research during the last two decades ReviewHB ; BookBH . Transport in non-crystalline, disordered materials cannot be explained by the classical theories of diffusion since anomalous behaviors (relative to what happens in the Euclidean domain) are the rule here. In particular, the mean-square displacement of a random walker $`r^2(t)`$ is no longer proportional to the time $`t`$ as occurs in uniform Euclidean systems of any dimension (Fick’s Law), but the more general anomalous diffusion law $`r^22Dt^{2/d_w}`$ holds for large times, $`r`$ being the Euclidean distance from the position of the random walker at $`t=0`$ , $`d_w`$ the anomalous diffusion exponent, and $`D`$ the diffusion constant. As the geometrical structure of real disordered media is very complex and varied, it is usually modeled by stochastic fractal lattices. Of these, the incipient percolation aggregate embedded in either two or three dimensions is the most widely used ReviewHB ; BookBH ; BookSA ; WeissEdHB ; Essam .

Statistical problems related to a single random walker have traditionally been the subject of more intensive research than those corresponding to $`N>1`$ interacting or independent random walkers Hughes ; Weiss . Of course, problems in which the walkers interact can not be analyzed in terms of the single-walker theory. However, there also exist other multiparticle problems that can not be analyzed in terms of the single-walker theory even though the walkers are independent. These problems have begun to be the target of in-depth studies during the last decades NAT ; PRAHL ; JCPSA ; SNTEuc ; SNTFrac ; WSL ; PRLYuste ; JPA2000 ; YAL ; JPAKR ; PREDK . Particular attention has been paid to (i) the evaluation of the average number of distinct sites visited up to time $`t`$ by $`N`$ independent random walkers all starting from the same origin in both Euclidean and fractal lattices NAT ; PRAHL ; JCPSA ; SNTEuc ; SNTFrac , and (ii) the description of the order statistic of the diffusion processes. This is the subject we address in this paper and can be stated as follows. A set of $`N`$ independent random walkers all initially placed at a given site (the origin) of a medium start to diffuse at time $`t=0`$. Eventually, at time $`t_{1,N}(z)`$ a random walker of this set reaches for the first time a site separated from the origin by the distance $`z`$. Next, a second random walker reaches at time $`t_{2,N}`$ a site at the same distance $`z`$ from the origin, and so on SingMol . Equivalently, we can understand $`t_{j,N}(z)`$ as the time taken (escape time or lifetime) by the $`j`$th particle out of a set of $`N1`$ to escape from a “spherical” region of radius $`z`$ centered at the starting site of diffusion. Our goal in this paper is to calculate the escape-time moments $`t_{j,N}^p(z)`$ when $`N1`$ random walkers diffuse in a disordered medium.

This order-statistic problem was solved for the one-dimensional lattice in Ref. WSL ; PRLYuste  and for some class of deterministic fractal substrates in Ref. PRLYuste . For $`d`$-dimensional Euclidean lattices, the form of the first moment of $`t_{1,N}`$ was guessed in Ref. JPA2000  and checked using simulations results for $`d=2`$ and $`d=3`$. The full solution of this problem for Euclidean media and for arbitrary $`p`$ and $`j`$ has been obtained recently YAL . For all these media, the $`p`$th moment of the time $`t_{j,N}`$ spend by the $`j`$th random walker in reaching the Euclidean distance $`r`$ was given in terms of an asymptotic series for large $`N`$ of the form

$$t_{j,N}^p\left(\frac{r}{\sqrt{2D}}\right)^{pd_w}(\mathrm{ln}N)^{p(1d_w)}\underset{n=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{n}{}}\tau _{nm}\frac{(\mathrm{ln}\mathrm{ln}N)^m}{(\mathrm{ln}N)^n},$$
(1)
where the coefficients $`\tau _{nm}`$ depend on $`p`$, $`j`$ and the substrate. The question we want to answer is whether this is also true for disordered media.

The prior knowledge of the short-time asymptotic expression of the survival probability $`\mathrm{\Gamma }(z,t)`$, which is defined as the probability that a random walker who starts at an origin site has not arrived at a spherical boundary of radius $`z`$ in the time interval $`(0,t)`$, was a key that allowed to find rigorously the full asymptotic approximation (1) for $`t_{j,N}^p`$ for Euclidean and some deterministic fractal media surviEucDetFrac  in Refs. WSL ; PRLYuste ; YAL . (Here, and henceforth, we use the term “short time” or “short-time regime” to mean that $`z/z^2(t)^{1/2}1`$, $`z^2(t)`$ being the mean square distance $`z`$ travelled by a single random walker by time $`t`$.) However, the use of the asymptotic procedures of Refs. WSL ; PRLYuste ; YAL  for estimating $`t_{j,N}^p`$ for disordered media is impeded by the fact that the value of the short-time survival probability is unknown for these substrates. In spite of this, in Ref. PRLYuste  it was conjectured that Eq. (1) is valid for disordered media too conjecture . However, Dräger and Klafter PREDK , using a scaling approach, have shown recently that

$$t_{1,N}=O\left[\mathrm{}^{d_w^{\mathrm{}}}(\mathrm{ln}N)^{1d_w^{\mathrm{}}}\right]$$
(2)
for ordered and disordered structures, which differs from (1) because, in general, $`d_w^{\mathrm{}}d_w`$. Here $`d_w^{\mathrm{}}=d_w/d_{\text{min}}`$ is the chemical diffusion exponent defined by the relation $`\mathrm{}\text{const}\times t^{1/d_w^{\mathrm{}}}`$, $`\mathrm{}(r)`$ being the chemical or topological distance defined as the length of the shortest path connecting two sites on a substrate that are separated by the Euclidean distance $`r`$, and $`d_{\text{min}}`$ is the fractal dimension of the shortest path on the fractal: $`\mathrm{}(r)\text{const}\times r_{\text{min}}^d`$ ReviewHB ; BookBH . Notice that $`d_{\text{min}}=1`$ for Euclidean lattices and for the deterministic fractals discussed in Ref. PRLYuste , so that, for these cases, Eqs. (1) and (2) agree. Unfortunately, the scaling approach of Dräger and Klafter is not precise enough to lead to a fully correct asymptotic expression of the form of Eq. (1). This can be seen, for example, because the asymptotic series for $`t_{1,N}`$ given by these authors in Ref. PREDK  does not agree with the rigorous asymptotic series reported in Refs. WSL ; PRLYuste  for the one-dimensional lattice.

In this paper we deal with the order statistic of the diffusion process in disordered media following the procedure previously carried out in Euclidean media WSL ; PRLYuste ; YAL  and in (some class of) deterministic fractal lattices PRLYuste ; YAL . Therefore, our first step must be to propose a short-time asymptotic expression for the survival probability. The functional form we take coincides with that obtained for Euclidean and renormalizable fractal lattices. Next, we apply the asymptotic methods already used for the Euclidean and deterministic fractal cases in order to obtain asymptotic expansions for the order-statistics moments $`t_{j,N}^m`$.

The paper is organized as follows. In Sec. II, a short-time asymptotic expression for the survival probability for disordered media is proposed. We check this conjecture in the two-dimensional incipient percolation aggregate and, on the way, estimate for this case the unknown parameters appearing in the proposed relationship. In Sec. III we give the asymptotic expressions of the moments of $`t_{j,N}`$ for stochastic fractals assuming that the survival probability for the short-time regime is given by the expression proposed in Sec. II. These theoretical results are compared with simulation data for the two-dimensional percolation aggregate. The paper ends with some conclusions and remarks.

## II SHORT-TIME SURVIVAL PROBABILITY

From now on, we will denote by $`z`$ either the chemical distance $`\mathrm{}`$ or the Euclidean distance $`r`$. As discussed in Sec. I we have to know the survival probability $`\mathrm{\Gamma }(z,t)`$, or equivalently the (boundary) mortality function $`h(z,t)=1\mathrm{\Gamma }(z,t)`$, for short times \[i.e., for $`z/z^2(t)^{1/2}1`$\] in order to calculate the escape-time moments $`t_{j,N}^m`$. This mortality function is the probability that a given random walker starting at the site $`z=0`$ will have been trapped by an absorbing boundary of radius $`z`$ by time $`t`$. In this paper we will consider media for which this mortality function (averaged over all realizations of the lattice if the lattice is stochastic) grows for $`\xi z/t^{1/d_w^z}1`$ as a stretched exponential with power-law corrective terms:

$$h(z,t)=h(\xi )A\xi ^{\mu v}\mathrm{exp}\left[c\xi ^v\right]\left\{1+h_1\xi ^v+\mathrm{}\right\}$$
(3)
where $`A`$, $`\mu `$, $`v`$, $`h_1`$ are characteristic parameters of the lattice. The anomalous diffusion coefficient $`d_w^z`$ is replaced by $`d_w`$ if the Euclidean distance is used ($`z=r`$) and by $`d_w^{\mathrm{}}`$ if $`z=\mathrm{}`$. There are good reasons to propose the functional form of Eq. (3). To start with, this relation holds for Euclidean lattices. Then, the mortality function was obtained in Ref. PRLYuste ; JPAYuste  for some class of deterministic fractals by using a renormalization procedure which involved only boundaries containing the nearest neighbors of the origin after successive decimations, and the result is in agreement with Eq. (3) too. In these media $`zr`$, $`\mu =1/2`$ and $`v=d_w/(d_w1)`$. (Although the results of renormalization cannot be directly applied to any origin and any arbitrary boundary, there are again good reasons to suspect that the functional form is the same Propag .) Also, a closely related quantity to $`h(z,t)`$, the (site) mortality function $`h_t(𝐫)`$ defined as the probability that a specific site $`𝐫`$ has been visited by a single random walker by time t, follows Eq. (3) for the two-dimensional percolation aggregate SNTFrac . Finally, one expects that the mortality function and the propagator $`P(z,t)`$ decay for $`\xi 1`$ in the same way JPAYuste ; Propag , and it is known that $`P(z,t)`$ decays as $`\mathrm{exp}(c\xi ^v)`$ for stochastic fractals, where $`v=d_w^z/(d_w^z1)`$ ReviewHB ; BookBH ; BookSA . We will devote the remaining part of this section to verify that $`h(\mathrm{},t)`$ and $`h(r,t)`$ behave as conjectured in Eq. (3) for the two-dimensional percolation aggregate.

### II.1 Mortality function for the two-dimensional percolation aggregate: Chemical distance

We start by estimating the mortality function $`h(\mathrm{},t)`$ when the “circular” absorbing boundary is placed at the chemical distance $`\mathrm{}`$ from the starting site. The numerical evaluation of this quantity is performed by the Chapman-Kolmogorov method (also called the exact enumeration method ReviewHB ; BookBH ). The circular boundary of absorbing traps is simulated by a set of special sites belonging to the cluster that absorbs all the probability density that enters them without giving back any probability to their neighbors. In the simulation we locate these boundaries at distances $`\mathrm{}=40,100,160`$, and perform three experiments of absorption for $`t=0,1,\mathrm{},t_{\text{max}}`$ with $`t_{\text{max}}=1000`$ on every aggregate built, one for each boundary. The resulting mortality function is averaged over $`2000`$ realizations of the percolation aggregate, which are generated by the Leath method BookBH ; Leath

In Fig. 1 we plot $`\mathrm{ln}[\mathrm{ln}h(\mathrm{},t)]`$ versus $`\mathrm{ln}\widehat{\xi }`$ with $`\widehat{\xi }\mathrm{}/t^{1/d_w^{\mathrm{}}}`$ and $`d_w^{\mathrm{}}=2.4`$. \[Hereafter we will put the symbol $``$ ( $``$ ) over quantities corresponding to the chemical (Euclidean) distance.\] The value of $`d_w^{\mathrm{}}`$ was taken from one of our previous works SNTFrac  and is in agreement with the values reported by other authors ReviewHB ; BHR ; Majid . If the conjecture in Eq. (3) is right we should observe the linear behavior $`\mathrm{ln}[\mathrm{ln}h(\mathrm{},t)]\mathrm{ln}\widehat{c}+\widehat{v}\mathrm{ln}\widehat{\xi }`$. Certainly the plots seem linear in Fig. 1 except for a portion in the range of large $`\widehat{\xi }`$ where the curves deviate upwards. This is a finite size effect (already analyzed in the case of the two-dimensional Sierpinski gasket in Ref. Propag ) associated with the existence of a minimum arrival time corresponding to a random walker who travels “ballistically” along a chemical path from the origin to the absorbing boundary, which in turn implies a maximum available value of $`\widehat{\xi }`$ in the simulations: $`\widehat{\xi }_{\text{max}}=\mathrm{}^{11/d_w^{\mathrm{}}}`$. This value is around $`\widehat{\xi }_{\text{max}}2.15`$ for $`\mathrm{}=40`$ and $`\widehat{\xi }_{\text{max}}2.96`$ for $`\mathrm{}=160`$. This apparently means that the reliable interval for numerical fitting is larger for $`\mathrm{}=160`$ but we must also take into account that the minimum value of $`\xi `$ attained in the simulations is $`\mathrm{}/t_{\text{max}}^{1/d_w^{\mathrm{}}}`$ which is proportional to $`\mathrm{}`$ (in our simulations $`t_{\text{max}}=1000`$). For this reason, in order to carry out the numerical fit, we have concatenated the simulation data for which the plots are almost linear. The result $`h(\mathrm{},t)=\widehat{A}\widehat{\xi }^{\widehat{\mu }\widehat{v}}\mathrm{exp}(\widehat{c}\widehat{\xi }^{\widehat{v}})`$ with $`\widehat{A}=1`$, $`\widehat{c}=0.9`$, $`\widehat{v}=d_w^{\mathrm{}}/(d_w^{\mathrm{}}1)1.714`$ and $`\widehat{\mu }=0.4`$ is also plotted in Fig. 1. The agreement is excellent.

### II.2 Mortality function for the two-dimensional percolation aggregate: Euclidean distance

We also simulated the mortality function $`h(r,t)`$ in the percolation aggregate when the circular absorbing boundary has a Euclidean radius $`r`$. The analysis of these results parallels those of the previous subsection. In Fig. 2 we have plotted the simulation results for the double logarithm of the mortality function $`\mathrm{ln}[\mathrm{ln}h(r,t)]`$ versus $`\mathrm{ln}\stackrel{~}{\xi }`$ (where $`\stackrel{~}{\xi }=r/t^{1/d^w}`$) for circular boundaries of radii $`r=35,55,75`$, with $`t_{\text{max}}=1000`$, and $`2000`$ aggregates to perform the average. The anomalous diffusion coefficient for that time range was taken as $`d_w=2.8`$ in agreement with that obtained by us in a previous work SNTFrac . In the long time limit, a slightly greater value $`d_w=2.87`$ has been found ReviewHB , but as our simulations are restricted to $`t1000`$ it is more reasonable to consider the (effective) value $`d_w=2.8`$ which better represents the diffusive behavior in this time range.

In contrast with the results for the chemical boundary case, now we get a poorer collapse to a single line for the curves corresponding to different values of $`r`$. This we attribute to the fact that a given Euclidean distance $`r`$ may correspond to many chemical distances $`\mathrm{}r`$ depending on the site and the particular aggregate we are using in the simulation. The result is that, in contrast with their chemical counterparts, the propagator, $`P(r,t)`$ BHR , the mortality function for a single trap site SNTFrac , and the mortality function for an absorbing spherical boundary, $`h(r,t)`$ (see Fig. 3), exhibit a broad distribution over the different percolation aggregate realizations. As a consequence, statistical precision requires more substrate averaging in the case of quantities defined in terms of Euclidean distance than quantities referred to the more natural chemical distance. Anyway, proceeding as in Sec. II.2, we linked the simulation results for $`r=35,55,75`$ where the plot is almost linear in order to perform an overall fit. Assuming the functional form of Eq. (3) we find that the data can be roughly described by $`h(r,t)=\stackrel{~}{A}\stackrel{~}{\xi }^{\stackrel{~}{\mu }\stackrel{~}{v}}\mathrm{exp}(\stackrel{~}{c}\stackrel{~}{\xi }^{\stackrel{~}{v}})`$ with $`\stackrel{~}{A}=1`$, $`\stackrel{~}{\mu }=1.8`$, $`\stackrel{~}{c}=1.65`$ and $`\stackrel{~}{v}d_w/(d_w1)1.56`$ (see Fig. 2).

## III ORDER STATISTICS ON A STOCHASTIC FRACTAL LATTICE

Now, with the functional form for $`h(z,t)`$ found in the preceding section, we can proceed with the evaluation of $`t_{j,N}^m`$ by means of the asymptotic technique already developed for Euclidean and deterministic fractal substrates. Only the key steps in the calculation will be outlined. The interested reader is referred to previous references for the details WSL ; PRLYuste ; YAL .

We will find it convenient in this section to write the mortality function as follows:

$$h(t)h(z,t)a(z)t^{\mu \beta }\mathrm{exp}\left[(t_0/t)^\beta \right]\left\{1+\overline{h}_1t^\beta \mathrm{}\right\}$$
(4)
where $`a(z)=Az^{\mu v}`$, $`\beta =v/d_w^z`$, $`t_0(z)=c^{1/\beta }z^{d_w^z}`$ y $`\overline{h}_1=h_1z^v`$. The generating function of the $`m`$th moment of the $`j`$th passage time, $`𝒰_{N,m}(\zeta )=_{j=1}^Nt_{j,N}^m\zeta ^{j1}`$, can be written as WSL

$$𝒰_{N,m}=\frac{m}{1\zeta }_0^{\mathrm{}}𝑑tt^{m1}\left\{[1h(t)+h(t)\zeta ]^N\zeta ^N\right\}.$$
(5)
Dropping the $`\zeta ^N`$ term we get

$$U_{N,m}\frac{m}{1\zeta }_0^{\mathrm{}}𝑑tt^{m1}\mathrm{exp}\left\{N\mathrm{ln}\left[1h(t)(1\zeta )\right]\right\}.$$
(6)
The point is that $`𝒰_{N,m}(\zeta )`$ and $`U_{N,m}(\zeta )`$ have the same Taylor series expansion up to the term of order $`\zeta ^{N1}`$, so that $`t_{j,N}^m`$ can also be estimated through the evaluation of this pseudo-generating function $`U_{N,m}(\zeta )`$. Proceeding as in Refs. WSL ; PRLYuste ; YAL , one gets:

$`U_{N,m}(\zeta )`$ $`=`$ $`{\displaystyle \frac{t_0^m}{1\zeta }}{\displaystyle \frac{1}{\mathrm{ln}^\alpha \lambda }}\{1+{\displaystyle \frac{\alpha \left(\mu \mathrm{ln}\mathrm{ln}\lambda \gamma \right)}{\mathrm{ln}\lambda }}+{\displaystyle \frac{\alpha }{2\mathrm{ln}^2\lambda }}[(1+\alpha )({\displaystyle \frac{\pi ^2}{6}}+\gamma ^2)+2\mu \gamma 2\overline{h}_1t_0^\beta `$ (7)
$`2\mu (\mu +(1+\alpha )\gamma )\mathrm{ln}\mathrm{ln}\lambda +(1+\alpha )\mu ^2\mathrm{ln}^2\mathrm{ln}\lambda ]+O\left({\displaystyle \frac{\mathrm{ln}^3\mathrm{ln}\lambda }{\mathrm{ln}^3\lambda }}\right)\}`$
where $`\alpha m/\beta `$ and $`\gamma 0.577215`$ is the Euler constant.

Once the generating function is known, the escape times and their moments are calculated straightforwardly because $`t_{j,N}^m`$ is simply the coefficient of $`z^{j1}`$ in the Taylor series expansion of $`U_{N,m}(\zeta )`$. Therefore, the $`m`$th moment of the first passage time of the first out of $`N1`$ diffusing particles is equal to $`U_{N,m}(0)`$, i.e.,

$`t_{1,N}^m`$ $``$ $`{\displaystyle \frac{t_0^m}{\mathrm{ln}^\alpha (\lambda _0N)}}\{1+{\displaystyle \frac{\alpha \left(\mu \mathrm{ln}\mathrm{ln}\lambda _0N\gamma \right)}{\mathrm{ln}\lambda _0N}}+{\displaystyle \frac{\alpha }{2\mathrm{ln}^2(\lambda _0N)}}[(1+\alpha )({\displaystyle \frac{\pi ^2}{6}}+\gamma ^2)+2\mu \gamma 2\overline{h}_1t_0^\beta `$ (8)
$`2\mu (\mu +(1+\alpha )\gamma )\mathrm{ln}\mathrm{ln}\lambda _0N+(1+\alpha )\mu ^2\mathrm{ln}^2\mathrm{ln}\lambda _0N]+O\left({\displaystyle \frac{\mathrm{ln}^3\mathrm{ln}\lambda _0N}{\mathrm{ln}^3\lambda _0N}}\right)\}`$
with $`\lambda _0=at_0^{\mu \beta }=Ac^\mu `$ and $`\alpha =m/\beta =mv/d_w^z`$ as before. Notice that the main term of this expression for $`m=1`$ agrees with the result given in PREDK  if we use $`v=d_w^z/(d_w^z1)`$.

The calculation of $`t_{j,N}^m`$ for $`j>1`$ is more involved. We begin with the identity $`\mathrm{ln}^n(1z)=n!_{i=n}^{\mathrm{}}(1)^iS_i(n)z^i/i!`$, where $`S_i(n)`$ are the Stirling numbers of the first kind Abra . Using this relation, the expansions of $`1/\mathrm{ln}^\alpha \lambda `$ and similar terms in Eq. (7) in powers of $`1/\mathrm{ln}\lambda _0N`$ can be found PRLYuste . After some algebra we have

$$t_{j,N}^mt_{1,N}^m+\frac{t_0^m\alpha }{\mathrm{ln}^{\alpha +1}(\lambda _0N)}\underset{n=1}{\overset{j1}{}}\frac{\mathrm{\Delta }_n(\alpha )}{n}$$
(9)
where $`j=2,3,\mathrm{}`$ and

$`\mathrm{\Delta }_n(\alpha )`$ $`=`$ $`1+{\displaystyle \frac{\alpha +1}{\mathrm{ln}\lambda _0N}}[(1)^n{\displaystyle \frac{S_n(2)}{(n1)!}}+\mu \mathrm{ln}\mathrm{ln}\lambda _0N`$ (10)
$`{\displaystyle \frac{\mu }{\alpha +1}}\gamma ]+𝒪({\displaystyle \frac{\mathrm{ln}^2\mathrm{ln}\lambda _0N}{\mathrm{ln}^2\lambda _0N}}).`$
Finally, we will quote the main term of the asymptotic expression for the variance, $`\sigma _{j,N}^2=t_{j,N}^2t_{j,N}^2`$, which is easily derived from Eqs. (8) and (9) yielding

$$\sigma _{j,N}^2=\frac{t_0^2}{\beta ^2}\frac{d_j}{(\mathrm{ln}\lambda _0N)^{2+2/\beta }}\left[1+O\left(\frac{\mathrm{ln}^3\mathrm{ln}\lambda _0N}{\mathrm{ln}\lambda _0N}\right)\right]$$
(11)
with

$$d_j=\left[\frac{\pi ^2}{6}\left(\underset{n=1}{\overset{j1}{}}\frac{1}{n}\right)^2+2\underset{n=1}{\overset{j1}{}}(1)^n\frac{S_n(2)}{n!}\right]$$
(12)
and $`j=1,2,\mathrm{}`$. We use the convention that the sums are equal to zero when the upper limit is zero. It is clear that the main and first corrective terms of $`t_{j,N}^2`$ are equal to those of $`t_{j,N}^2`$, so that only the difference between their second corrective terms contributes to the main term of the variance. For the sake of comparison with the simulation results it is more convenient to consider the quotient $`t_{j,N}/\sigma _{j,N}`$ whose expression is given by

$$\frac{t_{j,N}}{\sigma _{j,N}}=\beta d_j\mathrm{ln}N\left[1+O\left(\frac{\mathrm{ln}^3\mathrm{ln}N}{\mathrm{ln}N}\right)\right],$$
(13)
where $`d_1=(\pi ^2/6)^{1/2}`$, $`d_2=(\pi ^2/61)^{1/2}`$, $`d_3=(\pi ^2/65/4)^{1/2}`$,….

Let us now compare the theoretical predictions in Eqs. (8), (9) and (13) with simulation results in the two-dimensional percolation aggregate when an absorbing “circular” boundary is placed either at a given chemical distance $`z=\mathrm{}`$ or at a given Euclidean distance $`z=r`$.

### III.1 Absorbing boundary at a given chemical distance

The first passage time of the first few random walkers out of a set of $`N=2^i`$, $`i=2,3,\mathrm{},16`$ to a “circular” boundary of chemical radius $`\mathrm{}=50`$ was simulated on $`2000`$ aggregates. In Fig. 4 we plot the scaled simulation results of $`[t_{1,N}/\mathrm{}^{d_w^{\mathrm{}}}]^\delta `$, with $`\delta =1/(d_w^{\mathrm{}}1)`$, versus $`\mathrm{ln}N`$ for $`\mathrm{}=50`$ and compare them with the theoretical predictions of Eq. (8), namely,

$`\left[{\displaystyle \frac{t_{1,N}}{\mathrm{}^{d_w^{\mathrm{}}}}}\right]^\delta `$ $``$ $`\left[{\displaystyle \frac{\mathrm{ln}N}{\widehat{c}}}\right]^{\frac{\delta d_w^{\mathrm{}}}{\widehat{v}}}\{1`$ (14)
$`{\displaystyle \frac{\delta d_w^{\mathrm{}}}{\widehat{v}}}{\displaystyle \frac{\widehat{\mu }\mathrm{ln}\mathrm{ln}N\gamma \mathrm{ln}\widehat{A}\widehat{c}^{\widehat{\mu }}}{\mathrm{ln}N}}\}`$
where $`\delta =1/(d_w^{\mathrm{}}1)`$. Notice that $`\left(t_{1,N}/\mathrm{}^{d_w^{\mathrm{}}}\right)^\delta `$ depends linearly on $`\mathrm{ln}N`$ if $`\widehat{v}=d_w^{\mathrm{}}/(d_w^{\mathrm{}}1)`$. Figure 4 shows that this is indeed the case.

From Eq. (9) one gets

$$R\frac{t_{1,N}}{t_{j,N}t_{1,N}}\left(\underset{n=1}{\overset{j1}{}}\frac{1}{n}\right)^1\frac{\widehat{v}}{d_w^{\mathrm{}}}\mathrm{ln}N.$$
(15)
so that $`R`$ should be linear in $`\mathrm{ln}N`$ and independent of the radius of the chemical boundary $`\mathrm{}`$. This is confirmed in Fig. 5 where this quotient is plotted versus $`\mathrm{ln}N`$ for $`j=2,3`$ and $`\mathrm{}=50,80`$. A good superposition of the simulation data for $`\mathrm{}=50`$ and $`\mathrm{}=80`$ is observed. The simulation results are in good agreement with the prediction of Eq. (15) with $`\widehat{v}=d_w^{\mathrm{}}/(d_w^{\mathrm{}}1)`$ and $`d_w^{\mathrm{}}=2.4`$.

Finally, in Fig. 6 we show the quotient between the average lifetime of the $`j`$th particle and its variance, $`t_{j,N}/\sigma _{j,N}`$, versus $`\mathrm{ln}N`$ for the first, second and third random walkers to arrive at the chemical boundary with $`\mathrm{}=50`$. This is a remarkable figure that shows the crucial importance of the corrective terms on the order statistics formulae. Except for the first walker ($`j=1`$) we find striking discrepancies between the zeroth-order (main term) of the asymptotic expression in Eq. (13) and the simulation results. We attribute these discrepancies to the $`O\left[(\mathrm{ln}^3\mathrm{ln}N)/\mathrm{ln}N\right]`$ corrective terms not considered in the zeroth-order asymptotic expression. Let us explain why. In Fig. 7 we plot the function $`\mathrm{ln}^n\mathrm{ln}N/\mathrm{ln}^mN`$ for several values of $`n`$ and $`m`$. One sees that $`\mathrm{ln}^3\mathrm{ln}N/\mathrm{ln}N`$ cannot be neglected versus 1 even for values of $`N`$ larger than the Avogadro’s number. Therefore, it is not strange to find a poor agreement between the zeroth-order asymptotic expression and the simulation results. (In fact, what is truly unexpected is the good agreement of the zeroth-order theoretical prediction with simulation results for $`j=1`$.) Notice that if we had calculated another corrective term in the asymptotic expansion (7), then the term neglected in Eq. (13) would be of order $`\mathrm{ln}^4\mathrm{ln}N/\mathrm{ln}^2N`$, which, as shown in Fig. 7, would likely be much smaller than the $`\mathrm{ln}^3\mathrm{ln}N/\mathrm{ln}N`$ corrections. Unfortunately, corrective terms of order $`\mathrm{ln}^3\mathrm{ln}N/\mathrm{ln}N`$ involve parameters such as $`h_1`$ \[see Eq. (3)\] which are very difficult to estimate relying only on simulation results. One also observes in Fig. 7 that the relative errors of the zeroth- and first-order approximations to $`t_{j,N}`$ committed by ignoring terms of the form $`\mathrm{ln}\mathrm{ln}N/\mathrm{ln}N`$ and $`\mathrm{ln}^2\mathrm{ln}N/\mathrm{ln}^2N`$, respectively, will likely be small. This explains the good agreement between theory and simulation for lifetimes shown in Figs. 4 and 5.

### III.2 Absorbing boundary at a given Euclidean distance

We have also simulated the order-statistics of the arrival of $`N`$ random walkers at a circular boundary of Euclidean radius $`r`$. The analysis of these results parallels those of the previous subsection. The numerical results for $`[t_{1,N}/r^{d_w}]^{1/(d_w1)}`$ and $`t_{1,N}/t_{j,N}t_{1,N}`$ versus $`\mathrm{ln}N`$ are plotted in Figs. 8 and 9, respectively. The radii of the circular boundary are $`r=50`$ and $`r=75`$, and we used $`2000`$ aggregates to perform the substrate average. The theoretical predictions of Eqs. (14) and (15) \[replacing $`\mathrm{}`$ by $`r`$, $`d_w^{\mathrm{}}`$ by $`d_w`$, $`\delta `$ by $`1/(d_w1)`$ and $``$ by $``$\] lead to an agreement with numerical results that is not as good as that obtained when the chemical distance was used (see Figs. 4 and 5), although the linear behavior of $`T`$ and $`R`$ versus $`\mathrm{ln}N`$ and the superposition of the simulation data for different values of $`r`$ is still there. This slightly worse prediction is in accordance with the fact, discussed in Sec. II.2, that the collapse of the numerical mortality function for different distances $`r`$ to a single curve is not so nearly perfect as when chemical distances are used.

## IV CONCLUSIONS AND REMARKS

In this paper we studied the order statistics of the survival times (lifetimes or exit times) of $`N1`$ independent random walkers that are put initially on a stochastic lattice at the center of a “spherical” absorbing boundary of radius $`z`$. We found the moments $`t_{j,N}^m`$ of the lifetime of the $`j`$th random walker of the set of $`N`$ in terms of an asymptotic series that decays mildly in powers of $`1/\mathrm{ln}N`$. These theoretical results were compared with numerical simulations for the two-dimensional incipient percolation aggregate with the absorbing boundary placed either at a given chemical distance $`z=\mathrm{}`$ or at a given Euclidean distance $`z=r`$. The agreement found was reasonable and in accordance with the large size of the asymptotic corrective terms.

The theoretical approach relies on the knowledge of the mortality function $`h(z,t)`$ for short times. We assumed that the functional form of the short-time mortality function for a stochastic lattice averaged over all realizations coincides with that of the mortality function for Euclidean lattices and for the class of deterministic fractal lattices whose sites are isolated by their nearest neighbors. This is the main assumption of this paper and we showed that it is compatible with our simulation results for the two-dimensional percolation aggregate with Euclidean and chemical absorbing circular boundaries. This allows us to describe the order statistics of the diffusion process on many kinds of substrates (namely, Euclidean media, and deterministic and stochactic fractals) with the same approach and formulae.

A useful feature of the order statistics description is that it allows one to infer properties of the diffusive system (diffusion constant, number of diffusing particles, concentration of traps, effective dimension of the diffusive substrate,…) from the analysis of the behavior of those particles that are trapped first. This could be an advantage when it is impractical or impossible to wait until the entire reaction is over.

Finally, it should be noted that our simulation results for the stochastic fractal we have studied (the incipient percolation aggregate) are insufficient to perform a completely reliable numerical fit to the general form of the mortality function $`h(z,t)`$ with its dominant exponential and subdominant power-law terms simultaneously. Indeed, even to find the true exponent value $`v`$ of the dominant term for the propagator is not at all easy as past controversies about that true value reveal. Moreover, the value of a hypothetical power-law subdominant term \[equivalent to the term $`\xi ^{\mu v}`$ of Eq. (3)\] is still a matter of discussion JPAYuste ; Propag ; propuestaspropagador . Hence, one must expect that to find a completely reliable description of the short-time mortality function will be similarly difficult, and this task will therefore require more detailed simulations on larger lattices, with longer simulation times and, of course, averaging over many more realizations of the lattice than those we have used here. This is simply too much for a common workstation, but we think the effort on more powerful systems would be worthwhile: first, because of the interest in the percolation aggregate as a model of disordered media ReviewHB ; BookBH , and, second, because simulation of the mortality function (with circular boundary) is easier than that of the propagator and, since the asymptotic coefficients and exponents of these two functions are related, estimation of these parameters for the mortality function would help to solve the long-standing controversy on the exact asymptotic form of the propagator in fractal media.

###### Acknowledgements.

This work has been supported by the Miniterio de Ciencia y Tecnología (Spain) through Grant No. BFM2001-0718. SBY is also grateful to the DGES (Spain) for a sabbatical grant (No. PR2000-0116) and to Prof. K. Lindenberg and the Department of Chemistry of the University of California San Diego for their hospitality.