text:
# Absence of self-averaging in the complex admittance for transport through disordered media

## I INTRODUCTION

Disordered systems like amorphous semiconductors have a large number of remarkable properties inherent in their disorder, like the dispersive transport Scher91 . Hence, the disordered systems have wide areas of application like the photoreceptors of solar batteries and photocopiers Scharfe70 . So, they attract continuous research interests. In investigation of such remarkable properties, two kinds of statistics are to be considered. One is of course statistics of thermal fluctuation, which is considered in traditional statistical physics. The other is statistics of random structure of disordered systems. In order to treat this type of statistics, random structure is described with random variables and the probability distribution function of the random variables characterizes the substances.

It is usually assumed in studies of disordered systems that a sample used in experiments is sufficiently large. Hence, a measurement of any observable in such a system corresponds to an average over the ensemble of all realizations of the disorder. Quantities for which this assumption is valid are said to be self-averaging. The assumption means no sample dependence in the quantities. So, if the assumption is valid, reproducibility of experimental results for different samples is guaranteed. In addition, if the assumption is valid, in theoretical analyses only a disorder-averaged quantity is to be calculated since it should coincide with experimental results for any sample.

The expectation that the assumption of self-averaging is valid is based on the law of large numbers. The theorem says that the mean of independent random variables is equal to the expectation value with probability one. However, disordered systems such as amorphous semiconductors are strongly disordered, i.e., the expectation value of the random variables that describe structural disorder diverge. For example, the dispersive transport has been successfully explained by strongly disordered random variables Scher91 . In the theory of the dispersive transport, the transport is modeled by a random walk model on the disordered lattice where the dwell times of carriers at lattice sites are strongly disordered random variables.

The distribution of strongly disordered random variables is broad. Hence, properties of limited areas of samples like the maximal dwell time govern the macroscopic properties like the macroscopic mobility. One simple example is the finite contribution of the maximal term of a set of strong disordered random variables to the sum of the set (see Eq. (84) in appendix.). The properties of the limited areas with principal contribution fluctuate from one sample to another. Hence, the macroscopic properties of such materials may fluctuate. Thus, the assumption of self-averaging should be tested for its validity in strongly disordered systems. In fact, there are several reports on absence of self-averaging in transport phenomena in disordered media: Ref. 3 for mean square displacement, ref. 4 for the mean first passage time in the Sinai model.

In this paper, we study theoretically self-averaging properties of physical quantities obtained in a boundary perturbation experiment. The experimental technique was introduced recently to measure the optoelectrical properties of amorphous semiconductors Jongh96 . In this experimental method, an oscillatory perturbation with laser light is applied at one end of a sample and a response of the photocurrent from the other end is measured. We use as a generic model of the boundary perturbation experiment a random walk model in a one dimensional lattice with a boundary condition oscillating in time. The potential surface that carriers in the amorphous semiconductors feel is rugged due to structural disorder and it is called a rugged energy landscape. Consequently, the jump rates of the carriers are random variables and its probability distribution characterizes the disorder. It is known that the random jump rates obeys the power-law distribution with a negative exponent Pfister78  and it implies that the energy landscapes of the disordered media are extremely rugged.

It is shown that the complex admittance obtained in the experiment is not self-averaging when the jump rates $`w_i`$ are random variables distributed by the power-law distribution with the negative exponent. More precisely, the frequency-dependence of the disorder-averaged admittance $`\chi `$ disagrees with that of the admittance $`\chi `$ of any sample. It implies that the Cole-Cole plot of $`\chi `$ shows a different shape from that of the Cole-Cole plots of $`\chi `$ of each sample. The condition for absence of self-averaging is investigated with a toy version of the complex admittance. The complex admittance for higher dimensional media is investigated and shown to be non-self-averaging.

The rest of this paper is organized in the following way: In Sec. II, the random walk model used in the present study is explained. In terms of random walk, the complex admittance is related to a more physically relevant quantity to describe the transport in a disordered medium in Sec. III. In Sec. IV, absence of self-averaging for the complex admittance of one-dimensional media is shown. In Sec. V, absence of self-averaging is displayed in the Cole-Cole plot impressibly. In order to understand the condition for absence of self-averaging, a toy model is studied with the extended central limit theorem in Sec. VI. Absence of self-averaging for higher dimensional disordered media is shown in the similar way for the one-dimensional case in Sec. VII. Discussions and conclusions including consideration of possibilities that non-self-averaging properties are observed in real experiments are given in Sec. VIII. A part of this work has been reported in the brief note Kawasaki00  without consideration based on the probability distribution of the normalized mean first passage time, study of the toy model and the results for higher dimensional media.

## II A generic model of the boundary perturbation experiment

Response of a system to an external field is a standard tool to be utilized in the study of condensed matter. Recently, the intensity modulated photocurrent spectroscopy Jongh96  has been introduced to investigate transport properties in amorphous semiconductors. In the experiment, the electrode is illuminated with a laser and the incident light intensity is harmonically modulated at frequency $`\omega `$. The light absorbed by an optical transition generates the photocurrent. The light intensity absorbed in the electrode consists of a background intensity and an oscillating component with small amplitude $`\varphi (\omega )`$, which respectively gives rise to a steady state photocurrent and a harmonically varying photocurrent $`\mathrm{}i(\omega )`$. Since the harmonically varying photocurrent $`\mathrm{}i(\omega )`$ may show a phase shift with respect to the absorbed light flux $`\varphi (\omega )`$, the optoelectrical admittance $`\mathrm{}i(\omega )/\varphi (\omega )`$ becomes a complex number.

Since in this experiment an oscillatory perturbation is applied at one end of a system and a response is measured, the experimental technique belongs to a generic method called the boundary perturbation method (BPM). In the presence of a periodically forced boundary condition on one end of the system, the output from the other end of the system is in proportion to the perturbation in the linear response regime and the proportionality constant is called the admittance. The frequency dependence of the admittance contains various information of the dynamics of the system. The BPM is expected to provide useful information on transport properties of disordered media.

We introduce as a generic model of the BPM a random walk model in one-dimensional disordered lattice segment of $`N+1`$ sites with an oscillatory input current (see Fig. 1).

The lattice sites are denoted by integers, $`n=0,1,\mathrm{},N`$. The dynamics of a random walking particle can be described by a master equation for the probability $`P_n(t)`$ that the particle is at the site $`n`$ at time $`t0`$. The master equation for the model is written as

$$\frac{dP_n(t)}{dt}=w_{nn1}P_{n1}(t)(w_{n1n}+w_{n+1n})P_n(t)+w_{nn+1}P_{n+1}(t),$$
(1)
where $`w_{mn}`$ denotes the random jump rate of a particle from the site $`n`$ to the site $`m`$. Here, we assume that the random walking particle jumps only to the nearest neighbor sites. The probability distribution of jump rates characterizes the disordered medium. We introduce a perturbation at the left end, so that the equation for $`P_0(t)`$ is given by

$$\frac{dP_0(t)}{dt}=w_{\mathrm{1\; 0}}P_0(t)+w_{\mathrm{0\; 1}}P_1(t)+J(t).$$
(2)
$`J(t)`$ is the oscillatory current perturbation at the site $`0`$. The right end of the system is assumed to be a sink and the equation for $`P_{N1}`$ is given by

$$\frac{dP_{N1}(t)}{dt}=w_{N1N2}P_{N2}(t)(w_{N2N1}+w_{NN1})P_{N1}(t).$$
(3)
Since the current perturbation $`J(t)`$ oscillates in time around a positive average with the amplitude $`\mathrm{}J`$, the response of the output current $`w_{NN1}P_{N1}(t)`$ oscillates around its stationary state with the amplitude $`w_{NN1}\mathrm{}P_{N1}`$ at the same frequency with a phase-shift. The admittance is defined by the ratio of these two amplitudes

$$\chi _N(\omega )\frac{w_{NN1}\mathrm{}P_{N1}}{\mathrm{}J}.$$
(4)
The system of equations to determine the amplitudes $`\mathrm{}P_n`$ of oscillation of the probabilities is derived from the master equation Eqs. (1), (2) and (3);

$`i\omega \mathrm{}P_n`$ $`=`$ $`w_{nn1}\mathrm{}P_{n1}(w_{n1n}+w_{n+1n})\mathrm{}P_n`$
$`+w_{nn+1}\mathrm{}P\text{for }n=2,3,\mathrm{},N2,`$
$`i\omega \mathrm{}P_0`$ $`=`$ $`w_{\mathrm{1\; 0}}\mathrm{}P_0+w_{\mathrm{0\; 1}}\mathrm{}P_1+\mathrm{}J,`$
$`i\omega \mathrm{}P_{N1}`$ $`=`$ $`w_{N1N2}\mathrm{}P_{N2}(w_{N2N1}+w_{NN1})\mathrm{}P_{N1}.`$ (5)
## III The complex admittance and the first passage time distribution

We show below that the admittance can be related to the first passage time distribution $`F_{N0}(t)`$ (FPTD), which is the probability density that the particle which starts at the site $`0`$ at time $`0`$ arrives for the first time at the site $`N`$ at time $`t`$. Since the site $`N`$ is a sink, the FPTD is given by the output current from the site $`N`$ when there is no input current and the particle starts from the site $`0`$ at time $`0`$, i.e.

$$F_{N0}(t)=w_{NN1}P_{N1}(t).$$
(6)
The Fourier-Laplace transform of the master equation for the case is written as

$`i\omega \stackrel{~}{P}_n(i\omega )`$ $`=`$ $`w_{nn1}\stackrel{~}{P}_{n1}(i\omega )(w_{n1n}+w_{n+1n})\stackrel{~}{P}_n(i\omega )`$
$`+w_{nn+1}\stackrel{~}{P}(i\omega )\text{for }n=2,3,\mathrm{},N2,`$
$`i\omega \stackrel{~}{P}_0(i\omega )`$ $`=`$ $`w_{\mathrm{1\; 0}}\stackrel{~}{P}(i\omega )_0+w_{\mathrm{0\; 1}}\stackrel{~}{P}_1(i\omega )+1,`$
$`i\omega \stackrel{~}{P}_{N1}(i\omega )`$ $`=`$ $`w_{N1N2}\stackrel{~}{P}_{N2}(i\omega )(w_{N2N1}+`$ (7)
$`w_{NN1})\stackrel{~}{P}_{N1}(i\omega ).`$
This set of equations is identical to Eq. (5) divided by $`\mathrm{}J`$. It means that $`\mathrm{}P_n/\mathrm{}J=\stackrel{~}{P}_n(i\omega )`$ and hence the Fourier-Laplace transform of the FPTD $`\stackrel{~}{F}_{N0}(s=i\omega )`$ is identical to $`w_{NN1}\mathrm{}P_{N1}/\mathrm{}J`$. Thus, from Eqs. (4) and (6), we conclude that the admittance is equal to the Fourier-Laplace transform of the FPTD;

$$\chi _N(\omega )=\stackrel{~}{F}_{N0}(i\omega ).$$
(8)
We make the low-frequency expansion of the admittance to see the behavior near the static limit. Since the admittance is given by the Fourier-Laplace transform of the FPTD, the admittance at zero frequency equals to unity due to normalization of the FPTD;

$$\stackrel{~}{F}_{N0}(0)=_0^{\mathrm{}}𝑑tF_{N0}(t)=1.$$
(9)
The mean first-passage time (MFPT) $`\overline{t}_N`$, which is the first moment of the FPTD, is given by $`\overline{t}_N=id\chi (\omega =0)/d\omega `$. Thus, the low-frequency expansion of the admittance is given by

$$\chi _N(\omega )=1i\overline{t}_N\omega +O(\omega ^2).$$
(10)
The MFPT is given by Murthy89

$$\overline{t}_N=\underset{k=0}{\overset{N1}{}}\frac{1}{w_{k1k}}+\underset{k=0}{\overset{N2}{}}\frac{1}{w_{k+1k}}\underset{i=k+1}{\overset{N1}{}}\underset{j=k+1}{\overset{i}{}}\frac{w_{j1j}}{w_{j+1j}}.$$
(11)
## IV Absence of self-averaging for the complex admittance

Since we are interested in the complex admittance for amorphous semiconductors, we employ the power law probability distribution of the random jump rates Eq. (12), which is valid for amorphous semiconductors footnote1 :

$$\rho (w)=\{\begin{array}{cc}\alpha w^{\alpha 1}\hfill & \text{if }0<w<1\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(12)
There is possibility of the absence of self-averaging when disorder of the random dwell time at a site is strong as suggested in Sec. I. Since the dwell time is proportional to inverse of the jump rate, we are interested in the case where the first inverse moment of the random jump rate $`1/w`$ diverges, i.e. $`0<\alpha 1`$. It is important to note that the inequality $`0<\alpha <1`$ holds for amorphous semiconductors Pfister78 .

### IV.1 The probability distribution of the MFPT

In this subsection, we analyze the probability distribution of the MFPT Eq. (11). Since the MFPT gives the coefficient of the low frequency expansion of the admittance, the results of the analysis on the MFPT make it possible to determine whether the admittance is distributed with finite variance, i.e. non-self-averaging, or not. We show below that the complex admittance Eq. (10) is non-self-averaging when $`0<\alpha <1`$.

In order to analyze Eq. (11), we employ the random trap model general interest , where the jump rates $`w_{ij}`$ depends only on $`j`$. For the random barrier model, where the jump rates have the symmetry such that $`w_{ij}=w_{ji}`$, the quantitatively same results for self-averaging properties are obtained. We discuss below the probability distribution function of the normalized MFPT for the three cases of the value of $`\alpha `$. It is important to note that the three cases exhaust all possibilities of the value of $`\alpha `$. The choice of the normalization constants is suggested by the extended central limit theorem Feller66 . The theorem is summarized in appendix .1.

1. The case when $`0<\alpha <1`$ is considered. It is numerically shown that the probability distribution function of the normalized MFPT defined as $`\overline{t}_N/N^{1/\alpha +1}`$ converges to a distribution function with non-zero variance in the limit $`N\mathrm{}`$. Namely the normalized MFPT, $`lim_N\mathrm{}\overline{t}_N/N^{1/\alpha +1}`$, is a random variable even for the infinitely long chain. In addition, the distribution function behaves asymptotically as $`\rho (x)x^\gamma (1<\gamma <2)`$ for large $`x`$. It is shown in Sec. VII.2 that $`\gamma =\alpha +1`$.

Figure 2 shows the probability distribution function of $`\overline{t}_N/N^{1/\alpha +1}`$ when $`\alpha =1/2`$ in Eq. (12). It is clearly seen that the probability distribution functions for sufficiently long chains lie on a same curve.

Figure 3 shows the tail of the probability distribution which appears in Fig. 2. It behaves asymptotically as $`\rho (x)x^\gamma (\gamma =1.5<2)`$.
2. When $`\alpha =1`$, it is numerically shown that the probability distribution function of the normalized MFPT defined as $`\overline{t}_N/N^2b_N/2`$ converges to a distribution function with non-zero variance in the limit $`N\mathrm{}`$. Here, the centering constant $`b_N`$ is defined as

$$b_NN_1^{\mathrm{}}𝑑x\mathrm{sin}(x/N)x^2.$$
(13)
In appendix .2, we show that $`b_N`$ diverges as $`N\mathrm{}`$.

Figure. 4 shows the probability distribution function of the normalized MFPT when $`\alpha =1`$. It is clearly seen that the probability distribution functions lie on a same curve when the length of lattice $`N`$ is sufficiently large.
3. When $`1<\alpha `$, it is also shown numerically that the probability distribution function of the normalized MFPT defined as $`\overline{t}_N/N^2`$ converges to the Dirac delta function with the support at $`1/w/2`$. Figure 5 shows the probability distribution function of $`\overline{t}_N/N^2`$ when $`\alpha =1.5`$. One sees that the probability distribution concentrates on the point $`1/w/2=1.5`$ when the length $`N`$ is large.

### IV.2 Absence of self-averaging in the low frequency behavior of the complex admittance

By using the results for the probability distribution of the normalized MFPT, we discuss below frequency dependence of the complex admittance Eq. (10) for the three cases of the value of $`\alpha `$: $`0<\alpha <1`$, $`\alpha =1`$ and $`1<\alpha `$.

1. When $`0<\alpha <1`$, it is shown that the complex admittance is non-self-averaging. In order to see frequency dependence of the admittance, the long chain limit ($`N\mathrm{}`$) should be taken with keeping $`\overline{t}_N\omega `$ finite. Otherwise, the admittance becomes a trivial constant $`0`$ or $`1`$. Since it was shown that the normalized MFPT $`\overline{t}_N/N^{1/\alpha +1}`$ is finite with probability one, the long chain limit should be taken with fixing $`\omega N^{1/\alpha +1}`$ finite. Thus, from Eq. (10), the admittance for the infinitely long chain is given by

$$\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega N^{1/\alpha +1}}}{lim}\chi _N(\omega )=1i\stackrel{~}{t}\stackrel{~}{\omega }+\mathrm{}$$
(14)
where $`\stackrel{~}{t}`$ is the normalized MFPT defined as $`\stackrel{~}{t}lim_N\mathrm{}\overline{t}_N/N^{1/\alpha +1}`$ and $`lim_{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega N^{1/\alpha +1}}}`$ means the long chain limit with fixing $`\omega N^{1/\alpha +1}`$ at the finite value $`\stackrel{~}{\omega }`$. Equation (14) indicates that the admittance for the infinitely long chain is a function of the normalized MFPT $`\stackrel{~}{t}`$ and the scaled frequency $`\stackrel{~}{\omega }`$. From our result for the probability distribution of the normalized MFPT when $`0<\alpha <1`$, $`\stackrel{~}{t}`$ is a random variable. Thus, the admittance is also a random variable and hence the admittance is non-self-averaging.
2. When $`\alpha =1`$, it is shown below that the admittance is self-averaging in the low frequency range. By introducing the normalized MFPT $`\stackrel{~}{t}_N`$ defined as

$$\stackrel{~}{t}_N\frac{\overline{t}_N}{N^2}\frac{b_N}{2},$$
(15)
$`\overline{t}_N\omega `$ is rewritten as

$$\overline{t}_N\omega =N^2b_N\omega (\frac{2}{b_N}\stackrel{~}{t}_N+1)/2.$$
(16)
It was shown in Fig. 4 that the normalized MFPT $`\stackrel{~}{t}_N`$ is finite with probability one. Thus, in order to see frequency dependence, the long chain limit should be taken with keeping the scaled frequency $`\stackrel{~}{\omega }N^2b_N\omega /2`$ finite. Since $`b_N`$ diverges as $`N\mathrm{}`$, the complex admittance for the infinitely long chain is given by

$`\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=N^2b_N\omega /2}}{lim}\chi _N(\omega )`$ $`=`$ $`\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=N^2b_N\omega /2}}{lim}1iN^2b_N\omega /2({\displaystyle \frac{2}{b_N}}\stackrel{~}{t}_N+1)+\mathrm{}`$ (17)
$`=`$ $`1i\stackrel{~}{\omega }+\mathrm{}`$
Since Eq. (17) has no sample-dependence, the complex admittance is self-averaging in the low frequency region.
3. When $`1<\alpha `$, it is shown below that the complex admittance is self-averaging. In this case, it was shown that the scaled MFPT $`\overline{t}_N/N^2`$ is equal to $`1/w/2`$ in the limit $`N\mathrm{}`$. Thus, the admittance for the infinitely long chain is given by

$$\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega N^2}}{lim}\chi _N(\omega )=1\frac{1}{2}\frac{1}{w}\stackrel{~}{\omega }+\mathrm{}$$
(18)
It means that the complex admittance is self-averaging.

From Eq. (14), the low-frequency expansion of the disorder-averaged admittance is given by

$$\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega N^\beta }}{lim}\chi _N(\omega )=1i\stackrel{~}{t}\stackrel{~}{\omega }+\mathrm{}$$
(19)
When $`0<\alpha <1`$, the distribution function of the normalized MFPT $`\stackrel{~}{t}`$ behaves asymptotically as $`\rho (x)x^\gamma `$. It is shown later in Sec. VII.2 that $`\gamma =\alpha +1`$ and it means that $`1<\gamma <2`$. Hence, the expectation value $`\stackrel{~}{t}`$ of the normalized MFPT diverges. Since it implies that the coefficient of the first order of low $`\stackrel{~}{\omega }`$ expansion diverges, the disorder-averaged admittance is a non-analytic function of $`i\stackrel{~}{\omega }`$ which must behave as

$$\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega N^\beta }}{lim}\chi _N(\omega )1(i\stackrel{~}{\omega })^\mu $$
(20)
where $`0<\mu <1`$.

In order to test the foregoing observation, we numerically solved Eq. (5) for one-dimensional chain of $`20`$ sites. The jump rate $`w_{ij}`$ depends only on $`j`$ (the random trap model) and the jump rates are generated from the power-law distribution Eq. (12). Figure 6 shows the low-frequency behavior of the real part and the imaginary part of the admittances when $`\alpha =1/2`$.

One clearly sees that the low-frequency behavior of the disorder-averaged admittance is completely different from the admittance for any of the samples. On the other hand, figures 7 and 8 show the low-frequency behavior for $`\alpha =1.0`$ and for $`\alpha =1.5`$.

One clearly sees that the admittances for individual samples and their disorder-averaged admittance are both proportional to $`\stackrel{~}{\omega }`$. Although sample-dependence is observed, the sample-dependence is considered as the finite-size effect or effects from the terms of $`O(\omega ^2)`$.

## V Non-self-averaging in the Cole-Cole plot

In the literature of experiments of the BPM Jongh96 , the Cole-Cole plot of the admittance is employed for the analysis of experimental results. It is a parametric plot of the imaginary part of the admittance against the real part.

Figure 9 shows the Cole-Cole plot of $`\chi `$ and $`\chi `$ obtained by solving Eq. (5) when $`\alpha =1/2`$ in Eq. (12). One sees that the shape of the Cole-Cole plot of the disorder-averaged admittance is completely different from the shape of the Cole-Cole plot for each of samples. The boundary of the region where the Cole-Cole plot of the admittance of each sample scatters is the Debye semi-circle, whose center is at $`(1/2,0)`$ and the radius is $`1/2`$. In Fig. 9, solid circles denote the values of the admittance at a given frequency and a diamond is the average over these values. Since the admittance of each of samples at the same frequency (a solid circle) scatters around the Debye semi-circle, the admittance averaged at the same frequency (the diamond) comes inside the Debye semi-circle.

In order to analyze this non-self-averaging property of the admittance in the Cole-Cole plot, we prove rigorously that the Cole-Cole plots of each samples must lie outside the Debye semi-circle. We first note that the inverse of the admittance satisfies the following recursive relation:

$$\chi _n^1(\omega )=\left(\frac{i\omega }{w_{n+1n}}+1+\frac{w_{n1n}}{w_{n+1n}}\right)\chi _{n1}^1(\omega )\frac{w_{n1n}}{w_{n+1n}}\chi _{n2}^1(\omega ),$$
(21)
where $`\chi _n^1(\omega )`$ is the inverse of the admittance for a chain of length $`n`$. The recursive Eq. (21) is derived from the following two equations for the FPTD obtained from its definition: $`\stackrel{~}{F}_{i+1i}(s)=\stackrel{~}{\psi }_{i+1i}(s)+\stackrel{~}{\psi }_{i1i}(s)\stackrel{~}{F}_{i+1i1}(s)`$ and $`\stackrel{~}{F}_{i+1i1}(s)=\stackrel{~}{F}_{i+1i}(s)\stackrel{~}{F}_{ii1}(s)`$, where $`\stackrel{~}{F}_{ij}(s)`$ denotes the Laplace transform of the FPTD from site $`j`$ to site $`i`$ and $`\stackrel{~}{\psi }_{ij}(s)`$ denotes the Laplace transform of the waiting time distribution for the jump from site $`j`$ to site $`i`$. The recursive relation proves inductively that the inverse of the admittance obeys the following four inequalities in the frequency region $`[0,\omega _n]`$, where $`\omega _n`$ is the positive smallest zero of the real part of $`\chi _n^1`$ and is equal to the positive smallest zero of $`\chi _{n}^{}{}_{}{}^{}`$:

$$\{\begin{array}{c}\chi _{n}^{1}{}_{}{}^{\prime \prime }>0\hfill \\ \chi _{n1}^{1}{}_{}{}^{}>0\hfill \\ \chi _{n}^{1}{}_{}{}^{\prime \prime }>\chi _{n1}^{1}{}_{}{}^{\prime \prime }\hfill \\ \chi _{n}^{1}{}_{}{}^{}<\chi _{n1}^{1}{}_{}{}^{}.\hfill \end{array}$$
(22)
Since the second inequality shows $`\omega _n\omega _{n1}`$, the fourth inequality implies $`\chi _{n}^{1}{}_{}{}^{}<\chi _{n1}^{1}{}_{}{}^{}<\mathrm{}<\chi _{1}^{1}{}_{}{}^{}=1`$ in the region $`[0,\omega _n]`$. Since $`\chi _{n}^{1}{}_{}{}^{}=\chi _{n}^{}{}_{}{}^{}/(\chi _{n}^{}{}_{}{}^{}{}_{}{}^{2}+\chi _{n}^{}{}_{}{}^{\prime \prime }{}_{}{}^{2})`$, the inequality $`\chi _{n}^{1}{}_{}{}^{}<1`$ means $`(\chi _{n}^{}{}_{}{}^{}1/2)^2+\chi _{n}^{}{}_{}{}^{\prime \prime }{}_{}{}^{2}>(1/2)^2`$. Thus, the Cole-Cole plot of the admittance of each sample can exist only outside the Debye semi-circle when the frequency is smaller than the positive smallest zero of $`\chi _{n}^{}{}_{}{}^{}`$.

On the other hand, we can show in the following way that the Cole-Cole plot of the averaged admittance is located inside the Debye semi-circle when $`0<\alpha <1`$. We consider the angle $`\theta `$ between the tangent of the Cole-Cole plot at (1, 0) and the horizontal axis. From Eq. (20), $`\theta =\mu \pi /2`$, which is less than $`\pi /2`$ since $`\mu <1`$, is obtained. However the angle $`\theta `$ for the Debye semi-circle is equal to $`\pi /2`$. It means that the Cole-Cole plot of the averaged admittance exists inside the Debye semi-circle when $`0<\alpha <1`$. Since the Cole-Cole plot of the admittance of any sample exists outside the Debye semi-circle, the Cole-Cole plot clearly shows non-self-averaging property of the admittance.

When $`\alpha 1`$, the angle $`\theta `$ of the Cole-Cole plot of the averaged admittance is $`\pi /2`$ and is equal to that of the Cole-Cole plot of each of samples as seen in Figs. 10 and 11.

This property of the angle $`\theta `$ is derived from the regular behavior as $`\chi 1i\stackrel{~}{\omega }`$ (see Eq. (17) and Eq. (18).).

It is important to note that Fig. 10 suggests that the complex admittance is non-self-averaging when $`\alpha =1`$. One sees that the disorder-averaged complex admittance $`\chi `$ exists at different position from that of the complex admittance for each samples. Although it was shown that the low frequency behavior is self-averaging, it seems that the complex admittance at the finite frequency is non-self-averaging. Consequently, we conclude that the complex admittance is non-self-averaging when $`0<\alpha 1`$.

## VI Understanding of non-self-averaging of the complex admittance with a toy model

In the previous section, we discussed absence of self-averaging for the admittance when disorder of the jump rates is strong. In this section, in order to consider the conditions for absence of self-averaging we introduce a simple model, which is a toy version of the complex admittance and can be analyzed rigorously. This simple model will provide understanding of the reason why the admittance is non-self-averaging, since absence of self-averaging for the simple model is shown below with the same discussion as employed in the analysis of the complex admittance.

Analysis of the admittance is difficult because it does not consist of a mean of independent and identically distributed random variables. Hence, we consider a toy model of the admittance which is a simple function of a mean of independent and identically distributed random variables given by

$$ϵ_N(\omega )\frac{1}{1+i\omega _{i=1}^N(1/w_i)/N}.$$
(23)
We call it the dielectric constant. It is obvious that the Cole-Cole plot of each sample is the Debye semi-circle since the form of Eq. (23) is the same as that of the dielectric constant of the Debye relaxation. In addition, since $`|ϵ_N(\omega )|`$ is less than unity, $`ϵ_N(\omega )`$ is finite for any set of the random jump rates $`w_i`$.

In the analyses presented below, we use the following theorems of probability theory Feller66  summarized in appendix .1; Let $`\{x_i\}`$ be a set of independent and identically distributed random variables.

###### Theorem VI.1 (Law of large numbers)

When the expectation value $`x_i`$ exists, law of large numbers holds, i.e.

$$Prob\left(\left|\frac{_{i=1}^Nx_i}{N}x_j\right|>ϵ\right)0\text{as }N\mathrm{}.$$
(24)
###### Theorem VI.2 (Corollary of the extended central limit theorem)

When the expectation value $`x_i`$ diverges, in order that the probability distribution of $`_{i=1}^Nx_i/a_Nb_N`$ converges as $`N\mathrm{}`$ it is necessary that the probability density $`\rho (x_i)`$ is of the form

$$\rho (x)x^{\alpha 1}L(x)$$
(25)
for $`0<\alpha 1`$. Here, $`L(x)`$ is a slowly varying positive function which is precisely defined by Eq. (73). The normalization constant $`a_N`$ is chosen so that

$$N\frac{L(a_N)}{a_{N}^{}{}_{}{}^{\alpha }}C$$
(26)
as $`N\mathrm{}`$.

* When $`0<\alpha <1`$ in Eq. (25), the probability distribution of $`x_i/a_N`$ converges to the stable distribution with the characteristic exponent $`\alpha `$.
* When $`\alpha =1`$ in Eq. (25), the probability distribution of $`x_i/a_Nb_N`$ converges to the stable distribution with the characteristic exponent $`\alpha =1`$. The centering constant $`b_N`$ diverges as $`N\mathrm{}`$ footnote2 .

By using these theorems, we discuss sample-dependence of the dielectric constant $`ϵ(\omega )lim_N\mathrm{}ϵ_N(\omega )`$. The discussion is divided into two parts: the case where the expectation value $`1/w_i`$ is finite and the case where the expectation value diverges.

When the expectation value $`1/w_i`$ is finite, it is shown from law of large numbers that

$$ϵ(\omega )\underset{N\mathrm{}}{lim}\frac{1}{1+i\omega _{i=1}^N(1/w_i)/N}=\frac{1}{1+i\omega 1/w_i}.$$
(27)
Thus, there is no sample-dependence, i.e. the dielectric constant is self-averaging in the limit $`N\mathrm{}`$. It is obvious that the Cole-Cole plot is the Debye semi-circle, which corresponds to the Cole-Cole plot of the dielectric constant for each sample.

When the expectation value $`1/w_i`$ diverges (strong disorder), we discuss below the two cases of the value of $`\alpha `$ in Eq. (25): $`0<\alpha <1`$ and $`\alpha =1`$. In order to see frequency dependence of $`ϵ(\omega )`$, the macroscopic limit ($`N\mathrm{}`$) should be taken so that $`\omega 1/w_i/N`$ is finite. Otherwise, $`ϵ(\omega )`$ becomes a constant $`0`$ or $`1`$.

* When $`0<\alpha <1`$, from the extended central limit theorem, $`\stackrel{~}{T}lim_N\mathrm{}1/w_i/a_N`$ is a random variable with a limiting probability distribution and hence it is finite with probability one. Hence, the limit $`N\mathrm{}`$ should be taken so that $`\omega a_N/N`$ is finite;

$`ϵ(\omega )`$ $``$ $`\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega a_N/N}}{lim}{\displaystyle \frac{1}{1+i\omega (1/w_i)/N}}`$ (28)
$`=`$ $`\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega a_N/N}}{lim}{\displaystyle \frac{1}{1+i\omega a_N/N(1/w_i)/a_N}}`$
$`=`$ $`{\displaystyle \frac{1}{1+i\stackrel{~}{\omega }\stackrel{~}{T}}}.`$
It implies that the frequency-dependence is scaled with $`\omega a_N/N`$ as shown in Fig. 12.

Since $`\stackrel{~}{T}`$ is a random variable, $`ϵ(\omega )`$ is also a random variable and hence $`ϵ(\omega )`$ is non-self-averaging. In addition, since it is known that the probability density $`\rho (x)`$ of a stable distribution with characteristic exponent $`\alpha `$ behaves asymptotically as $`\rho (x)x^{\alpha 1}`$, the probability density $`\rho (\stackrel{~}{T})`$ behaves as $`\rho (\stackrel{~}{T})\stackrel{~}{T}^{\alpha 1}`$. Hence, the expectation value $`\stackrel{~}{T}`$ diverges. Since it implies that the coefficient of the first order of low $`\stackrel{~}{\omega }`$ expansion diverges, the disorder-averaged dielectric constant is a non-analytic function of $`i\stackrel{~}{\omega }`$ which must behave as $`ϵ(\omega )1(i\stackrel{~}{\omega })^\mu (0<\mu <1)`$. Thus, from the same argument employed in the discussion for the Cole-Cole plot of the admittance in Sec. V the Cole-Cole plot of the disorder-averaged dielectric constant appears inside the Debye semi-circle which corresponds to the Cole-Cole plot of the dielectric constant for each sample (see Fig. 13).
* When $`\alpha =1`$, from the extended central limit theorem, $`1/(w_i)/a_Nb_N`$ is a random variable with a limiting probability distribution and hence it is finite with probability one. Furthermore, $`b_N`$ diverges as $`N\mathrm{}`$ as shown in appendix .2. Hence, the limit $`N\mathrm{}`$ should be taken so that $`\omega a_Nb_N/N`$ is finite;

$`ϵ(\omega )`$ $``$ $`\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega a_Nb_N/N}}{lim}{\displaystyle \frac{1}{1+i\omega (1/w_i)/N}}`$ (29)
$`=`$ $`\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{\stackrel{~}{\omega }=\omega a_Nb_N/N}}{lim}{\displaystyle \frac{1}{1+i\omega a_Nb_N/N\left[\left((1/w_i)/a_Nb_N\right)/b_N+1\right]}}`$
$`=`$ $`{\displaystyle \frac{1}{1+i\stackrel{~}{\omega }}}.`$
Hence, $`ϵ(\omega )`$ is self-averaging.

Consequently, we conclude that only when $`0<\alpha <1`$ the dielectric constant is non-self-averaging.

In order to test the foregoing discussion, we present an example of the above analysis where the random jump rates distributed by the power-law distribution Eq. (12) when $`\alpha =1/2`$. From the extended central limit theorem, it is known that the distribution function of $`\stackrel{~}{T}`$ is given by the stable distribution with $`\alpha =1/2`$. From Eq. (83), the probability distribution $`g(\stackrel{~}{T})`$ of $`\stackrel{~}{T}`$ behaves as $`g(\stackrel{~}{T})\stackrel{~}{T}^{3/2}`$. Hence, when $`\stackrel{~}{\omega }0`$ the disorder-averaged dielectric constant behave as

$$ϵ(\omega )1\sqrt{i\stackrel{~}{\omega }}.$$
(30)
On the other hand, the dielectric constant for each sample behaves as

$$ϵ(\omega )1i\stackrel{~}{\omega }.$$
(31)
Thus, the low-frequency behavior is non-self-averaging. Absence of self-averaging is also seen in the Cole-Cole plot (Fig. 13).

## VII Higher dimensional disordered media

It is important to note that the system treated above is purely one-dimensional and it is still an open problem whether the admittance for higher dimensional media is self-averaging.

In this section, we investigate properties of the admittance for higher dimensional disordered media in the same manner for one-dimensional media. At first, the statistics of the MFPT for higher dimensional media is analyzed analytically. With the results, it is shown that the complex admittance is non-self-averaging. After that, the theoretical results are confirmed numerically.

For its feasibility, we investigate a site-disordered model called the random trap model, which is defined as a random walk model where the jump rate $`w_{ij}`$ from site $`j`$ to site $`i`$ depends only on $`j`$.

We assume that the higher dimensional medium is a $`N\times M`$ square (or cubic) lattice. Here, $`N`$ is the linear dimension of the medium in direction of the current and a number of grid points on the cross section perpendicular to the current direction is denoted as $`M`$. The periodic boundary condition on the side ends is assumed, i.e., the medium is a hyper-cylinder, and the oscillatory input current is introduced uniformly on the left end of the medium. An example of two-dimension is illustrated in Fig. 14.

The probability distribution for the jump rates is assumed power-law distribution Eq. (12) with the exponent $`\alpha `$. Since we are interested in amorphous semiconductors, our analyses are limited to the case $`0<\alpha <1`$.

### VII.1 The analytic expression for the MFPT of higher dimensional media

At first, we introduce a general analytical expression for the MFPT in terms of the backward equation. After that, by using the expression, the results for the one-dimensional MFPT and the extended central limit theorem, the probability distribution of the higher-dimensional MFPT is investigated. Since the analyses are general, results are applicable both to two-dimensional and three-dimensional media.

The backward equation, i.e., the time-development equation for a averaged quantity, is derived in the following way.

Let $`P_{\mathrm{𝐣𝐢}}(t)`$ be the probability that the random walker starting from the site $`𝐢`$ at time $`0`$ is found at the site $`𝐣`$ at time $`t`$. We assume that the master equation for the probability is written as

$$\dot{𝐏_𝐢}(t)=H𝐏_𝐢(t).$$
(32)
Here, $`𝐏_𝐢`$ is the vector whose $`𝐣`$-th component is $`P_{\mathrm{𝐣𝐢}}`$ with fixing the starting site $`𝐢`$. Furthermore, $`H`$ is a linear operator given in terms of the jump rates.

We consider the time-development equation of the average of the site-dependent quantity $`f_𝐣`$. The average $`g_𝐢`$ is defined by

$$g_𝐢(t)\underset{𝐣}{}f_𝐣P_{\mathrm{𝐣𝐢}}(t)$$
(33)
and hence the time-derivative is written as

$$\dot{g}_𝐢(t)=\underset{𝐣}{}f_𝐣\dot{P}_{\mathrm{𝐣𝐢}}(t)=\underset{\mathrm{𝐣𝐤}}{}f_𝐣H_{\mathrm{𝐣𝐤}}P_{\mathrm{𝐤𝐢}}(t).$$
(34)
Taken the limit $`t0`$,

$$\dot{g}_𝐢(0)=\underset{\mathrm{𝐣𝐤}}{}f_𝐣H_{\mathrm{𝐣𝐤}}\delta _{\mathrm{𝐤𝐢}}=\underset{𝐣}{}f_𝐣H_{\mathrm{𝐣𝐢}}=\underset{𝐣}{}\left({}_{}{}^{t}H\right)_{\mathrm{𝐢𝐣}}f_𝐣.$$
(35)
Let $`𝐠(t)`$ be the vector whose component is $`g_𝐢(t)`$, then

$$\dot{𝐠(0)}=^tH𝐟.$$
(36)
By using the above result, $`𝐠(t)`$ is Maclaurin-expanded as

$$𝐠(t)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{n!}𝐠^{(n)}(0)t^n=\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{n!}(^tH)^n𝐟t^n.$$
(37)
Taken the time-derivative,

$$\dot{𝐠}(t)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{n}{n!}(^tH)^n𝐟t^{n1}=^tH\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{n!}^n𝐟t^n=^tH𝐠(t).$$
(38)
The last equality gives the backward equation:

$$\dot{𝐠}(t)=^tH𝐠(t).$$
(39)
As an example of the average $`g_𝐢(t)`$, we consider the probability that the random walker staring at the site $`𝐢`$ on the end illuminated with the laser light is found in the medium at the time $`t(>0)`$. The survival probability $`n_𝐢(t)`$ is defined as

$$n_𝐢(t)\underset{𝐣}{}P_{\mathrm{𝐣𝐢}}(t).$$
(40)
Since the probability can be considered as a averaged quantity, $`n_𝐢(t)`$ obeys the backward equation:

$$\dot{𝐧}(t)=^tH𝐧(t),$$
(41)
where $`𝐧(t)`$ is the vector whose component is $`n_i(t)`$.

In terms of the survival probability, the first passage time distribution (FPTD) $`F(t)`$ is written as

$$F(t)=\frac{d}{dt}\underset{𝐢}{}n_𝐢(t)/M,$$
(42)
where $`M_𝐢1`$ and the sum is taken over the grid points on the end of the medium. Hence, the MFPT $`\overline{t}`$ is also written as

$$\overline{t}_0^{\mathrm{}}tF(t)𝑑t=\frac{1}{M}\underset{𝐢}{}_0^{\mathrm{}}t\frac{d}{dt}n_𝐢(t)𝑑t=\frac{1}{M}\underset{𝐢}{}_0^{\mathrm{}}n_𝐢(t)𝑑t=\frac{1}{M}\stackrel{~}{n}_𝐢(0).$$
(43)
Here, $`\stackrel{~}{n}_𝐢(0)`$ is the limit $`s0`$ of the Laplace transform $`\stackrel{~}{n}_𝐢(s)`$. The Laplace-transform of the backward equation is written as

$$s\stackrel{~}{𝐧}(s)𝐧(t=0)=^tH\stackrel{~}{𝐧}(s).$$
(44)
Taken the limit $`s0`$,

$${}_{}{}^{t}H\stackrel{~}{𝐧}(0)=𝐧(t=0),$$
(45)
where $`𝐧(t=0)=^t(1,1,\mathrm{})`$. By solving Eq. (45), the MFPT is obtained as Eq. (43).

In the rest of this section, we concentrate on the random trap model. For the $`d`$-dimensional random trap model, the analytical expression of the MFPT can be easily obtained by Eq. (45) in the following way.

Each components in the same *column* of the matrix $`H`$ are written only in terms of the jump rates of the same site. It means that each components in the same *line* of the transposed matrix $`{}_{}{}^{t}H`$ are written only with the jump rates of the same site. Hence, divided the each lines of Eq. (45) by the jump rates, we obtain

$`{}_{}{}^{t}H_{0}^{}\stackrel{~}{𝐧}(0)=\left(\begin{array}{c}\mathrm{}\\ 1/w_𝐤\\ \mathrm{}\end{array}\right).`$ (49)
Here, $`{}_{}{}^{t}H_{0}^{}`$ is the transpose of the time-development operator when the jump rates of the all sites are equal to unity. From Eq. (49), $`\stackrel{~}{𝐧}(0)`$ is a homogeneous linear function of the inverse of the jump rates.

We denote the grid points in $`d`$-dimensional medium with two indexes, $`n`$ and $`𝐦`$. One index $`n`$ is the label of the cross sections normal to the current direction and the other index $`𝐦`$ is the $`d1`$ dimensional vector to denote a grid point on the cross section. Let the cross section at the left end, which is illuminated with laser light, be numbered $`0`$ and the cross section at the right end be numbered $`N1`$. Each cross sections have $`M`$ grid points.

Since the boundary condition at the side ends is periodic, each grid points in the same cross section are identical geometrically and each jump rates $`w_{n𝐦}`$ in the same cross section $`n`$ have the same contribution to the MFPT. Hence, the results obtained above means that the MFPT for the $`d`$-dimensional medium, $`\overline{t}_d`$ is written as

$`\overline{t}_d`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{𝐦}{}}{\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \frac{a_n}{w_{n𝐦}}}`$ (50)
$`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{𝐦}{}}\left({\displaystyle \frac{a_0}{w_{0𝐦}}}+{\displaystyle \frac{a_1}{w_{1𝐦}}}+\mathrm{}+{\displaystyle \frac{a_{N1}}{w_{N1𝐦}}}\right),`$
where $`\{a_n\}`$ are constants independent of configuration of disorder. If all the jump rates in a cross section are equal, i.e., $`w_{n𝐦}=w_{n𝐦^{}}w_n`$, the $`d`$-dimensional MFPT $`\overline{t}_d`$ is equal to the one-dimensional MFPT $`\overline{t}_1`$, which is given as

$$\overline{t}_1=\frac{N}{w_0}+\frac{N1}{w_1}+\mathrm{}+\frac{1}{w_{N1}}.$$
(51)
Comparison of Eq. (50) to Eq. (51) shows that

$$a_n=Nn.$$
(52)
Consequently, the MFPT for $`d`$-dimensional medium is given as

$`\overline{t}_d`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{𝐦}{}}{\displaystyle \underset{n=0}{\overset{N1}{}}}{\displaystyle \frac{Nn}{w_{n𝐦}}}`$ (53)
$`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{𝐦}{}}\left({\displaystyle \frac{N}{w_{0𝐦}}}+{\displaystyle \frac{N1}{w_{1𝐦}}}+\mathrm{}+{\displaystyle \frac{1}{w_{N1𝐦}}}\right).`$
With respect to the MFPT, the $`d`$-dimensional medium is a bundle of one-dimensional chain and hence the MFPT is rewritten as

$$\overline{t}_d=\frac{1}{M}\underset{𝐦}{}t_1^{(𝐦)},$$
(54)
where $`t_1^{(𝐦)}`$ is the MFPT for a one-dimensional medium defined as

$$t_1^{(𝐦)}\underset{n=0}{\overset{N1}{}}\frac{Nn}{w_{n𝐦}}.$$
(55)
### VII.2 The statistics of the MFPT for higher dimensional media

In order to consider the probability distribution of $`\overline{t}_d`$, the following postulate, which is expected to be true from the results of analyses for the one-dimensional case, is introduced: Let $`\{u_n\}`$ be a set of independent and identically distributed random variables. The tail of the probability density is assumed to be

$$\rho _u(u_n)u_{n}^{}{}_{}{}^{\alpha 1},$$
(56)
where $`0<\alpha <1`$. Then, we postulate that the random variable $`v`$ defined as

$$v\frac{1}{N^{1/\alpha +1}}\underset{n=0}{\overset{N1}{}}(Nn)u_n$$
(57)
has a normalized probability distribution function and the tail obeys the power law.

By making use of the postulate, we analyze the statistics of $`\overline{t}_d`$ in two different ways and determine the asymptotic form of the probability distribution of $`\overline{t}_d`$.

At first, we take the way where the sum in terms of $`n`$ is done first. From our postulate, the random variable $`\stackrel{~}{t}_𝐦`$ defined as

$$\stackrel{~}{t}_𝐦\frac{1}{N^{1/\alpha +1}}t_1^{(𝐦)}$$
(58)
has the normalized probability distribution function and it obeys the power law:

$$\rho _{\stackrel{~}{t}}(\stackrel{~}{t}_𝐦)\stackrel{~}{t}_{𝐦}^{}{}_{}{}^{\alpha ^{}1},$$
(59)
where $`\alpha ^{}`$ is a positive constant. Since $`\overline{t}_d=M^1N^{1/\alpha +1}_𝐦\stackrel{~}{t}_𝐦`$, from the extended central limit theorem, the random variable $`x`$ defined as

$$x\frac{1}{M^{1/\alpha ^{}1}N^{1/\alpha +1}}\overline{t}_d$$
(60)
has a normalized probability distribution function and its tail asymptotically behaves as

$$\rho _x(x)x^{\alpha ^{}1}.$$
(61)
As the other way, we evaluate the sum in terms of $`𝐦`$ first. The extended central limit theorem shows that the random variable $`y_n`$ defined as

$$y_n\frac{1}{M^{1/\alpha }}\underset{𝐦}{}\frac{1}{w_{n𝐦}}$$
(62)
has a normalized probability distribution function whose asymptotic form is

$$\rho _y(y_n)y_{n}^{}{}_{}{}^{\alpha 1}.$$
(63)
Furthermore, since $`\overline{t}_d=M^{1/\alpha 1}_{n=0}^{N1}(Nn)y_n`$, our postulate means that the random variable $`z`$ defined as

$$z\frac{1}{M^{1/\alpha 1}N^{1/\alpha +1}}\overline{t}_d$$
(64)
has a probability distribution function whose tail behaves as

$$\rho _z(z)z^{\alpha 1}.$$
(65)
Here, we compare the two ways of normalization Eqs. (60), (64). Since the two ways of normalization to obtain the normalized distribution function should be identical,

$$\alpha ^{}=\alpha .$$
(66)
Summarized our results, the MFPT $`\overline{t}_d`$ is normalized as

$$\overline{t}_dM^{1/\alpha 1}N^{1/\alpha +1}z$$
(67)
and the probability distribution function has the power-law tail as

$$\rho _z(z)z^{\alpha 1}.$$
(68)
Furthermore, from Eq. (66), the tail of the probability distribution function of $`v`$ defined by Eq. (57) behaves as

$$\rho _v(v)v^{\alpha 1}.$$
(69)
Since the results are based on the postulate, we confirm numerically them. Figure 15 shows that the probability distribution function of $`z`$ converges to a unique distribution when $`N`$ and $`M`$ are increased. The power-law tail of the distribution is shown in Fig. 16. It is clearly seen that the exponent is equal to $`1.5`$ when $`\alpha =0.5`$.

### VII.3 Non-self-averaging complex admittances for higher dimensional media

By using the results obtained in the previous subsection, Eqs. (67), (68), we show in the same way as for the one-dimensional case that the complex admittance for higher dimensional media is also non-self-averaging when $`0<\alpha <1`$.

The normalization of the MFPT, Eq. (67), means that the complex admittance for an infinitely large higher dimensional medium is given as

$$\underset{N,M\mathrm{}}{lim}\chi _{N,M}(\omega )=1iz\stackrel{~}{\omega }+\mathrm{}$$
(70)
in terms of the scaled frequency $`\stackrel{~}{\omega }`$ defined as

$$\stackrel{~}{\omega }\omega M^{1/\alpha 1}N^{1/\alpha +1}.$$
(71)
Here, the infinite volume limit $`N,M\mathrm{}`$ is taken with fixing $`\stackrel{~}{\omega }`$ at a finite value. Since the normalized MFPT $`z`$ is a random number, the complex admittance is dependent on each samples even for the infinitely large medium. That is, *the complex admittance is non-self-averaging.*

Furthermore, from Eq. (68), it is shown that $`z=\mathrm{}`$. It means non-analyticity of the complex admittance about $`\stackrel{~}{\omega }=0`$ and hence

$$\underset{N,M\mathrm{}}{lim}\chi _{N,M}(\omega )1(i\stackrel{~}{\omega })^\mu ,$$
(72)
where $`\mu `$ is a constant such that $`0<\mu <1`$ (see Fig. 17).

The results obtained above are all concerned with the low frequency behavior. We investigate numerically the non-self-averaging characteristics in the all frequency range with the Cole-Cole plot. Figure 18 shows the Cole-Cole plots of the admittances for one hundred realizations of $`4\times 4,10\times 10,30\times 30`$ square lattices. These numerical results show that the sample-dependence of the admittance seems to exist even when both the length of the lattice and the width are large and hence it implies that the admittance is non-self-averaging even in higher dimensional cases.

## VIII Discussion and conclusions

In conclusion, we presented a generic stochastic dynamical model in one-dimensional medium of the BPM and showed that the admittance has important information about transport in disordered media, e.g. the FPTD. Thus, the BPM will be a powerful technique to investigate the dynamical process in random media. We also showed that the admittance is non-self-averaging when the jump rates are random variables with the power-law distribution as $`\rho (w_i)w_{i}^{}{}_{}{}^{\alpha 1}(0<\alpha 1)`$. It is important to note that this power-law distribution is realized in amorphous semiconductors.

On the other hand, when $`1\alpha `$, the low frequency behavior of $`O(\omega )`$ has no sample dependence as shown in Eqs. (17), (18). However, with respect to the terms of $`O(\omega ^2)`$, the second moment of the inverse of the jump rate $`1/w^2`$, which diverges when $`1\alpha 2`$, is concerned. It implies that some anomalous statistical properties may be observed. Properties for the full frequency range were analyzed with the Cole-Cole plots. In the marginal case when $`\alpha =1`$, the non-self-averaging property appears as shown in Fig. 10. However, Fig. 11 shows that such effects from terms of $`O(\omega ^2)`$ are small when $`1<\alpha `$.

When the complex admittance is non-self-averaging, the behavior of the averaged admittance over infinite number of samples does not coincide with that of the admittance observed in experiments. In order to see the non-self-averaging behavior described above one needs to work in one-dimensional system or the case where the translational invariance is violated in only one direction.

Furthermore, in order to investigate possibility that the non-self-averaging properties of the admittance are observed in real experiments on higher dimensional media, we performed theoretical and numerical analyses of the admittance on higher dimensional disordered media. The analyses on anomalous statistics of the MFPT show that the non-self-averaging characteristics are observed in the low-frequency behavior of the complex admittance of higher-dimensional media. The numerical results confirm the conclusion that the admittance is non-self-averaging. Hence, it suggests possibility that the non-self-averaging properties of the admittance of the BPM may be observed in the real experiments on higher dimensional media.

On the other hand, non-self-averaging behavior, i.e. strong sample-dependence may be observed in other experiments. For example, it is known that anomalous system-size dependence of the mobility, $`\mu L^{11/\alpha }`$, is observed in the time-of-flight experiment of amorphous semiconductors Pfister78 . The anomalous system-size dependence is due to the fact that the smallest value of random jump rates depends on the length of a chain Kehr98 . Since the smallest value dominates the mobility, the mobility is also a random variable. Thus, the anomalous system-size dependence of the mobility implies that the mobility is not self-averaging.

###### Acknowledgements.

I wish to acknowledge valuable discussions with Akira Yoshimori, Hiizu Nakanishi, Kazuo Kitahara, Osamu Narikiyo and Miki Matsuo. One of us (M. K.) is thankful to Institute für Festkörperforschung, Forschungszentrum Jülich for its hospitality, where part of the present work was done. This work was supported in part by the Grant-in-Aid for Scientific Research of Ministry of Education, Culture, Sports, Science and Technology.

## Law of large numbers and the extended central limit theorem

In this appendix, we summarize the notions of probability theory, partly rewritten for our convenience, used to analyze the complex admittance and the “dielectric constant” in Secs. V and VI. After that, we show that $`b_N`$ Eq. (13) diverges as $`N\mathrm{}`$.

### .1 The extended central limit theorem

Here, we summarize theorems on probability distribution of normalized sum of random variables Feller66 .

The following law of large numbers tells that a mean converges to an expectation value even when a variance does not exist. In this section, $`S_n`$ denotes the sum $`_{i=1}^nX_i`$.

###### Theorem .1 (Strong law of large numbers)

Let the $`X_k`$ be independent random variables with a common distribution $`F`$. If they have an expectation $`\mu `$ then $`S_n/n\mu `$ with probability one.

In order to describe the extended central limit theorem, some concepts are defined.

###### Definition .1 (Domain of attraction)

The distribution $`F`$ of the independent random variables $`X_k`$ belongs to the domain of attraction of a distribution $`R`$ if there exist the norming constant $`a_n>0`$ and the centering constant $`b_n`$ such that the distribution of $`a_n^1S_nb_n`$ tends to $`R`$.

###### Definition .2 (Slow variation)

A positive function $`L`$ defined on $`(0,\mathrm{})`$ varies slowly at infinity if

$$\frac{L(st)}{L(s)}1$$
(73)
for each $`t>0`$ as $`s\mathrm{}`$.

It implies that the leading term of $`L(st)`$ is $`L(s)`$ when $`s`$ is large. Constant functions and $`\mathrm{log}(x)`$ are examples.

In addition, the function $`F(x)`$ is assumed to be related to the probability density $`\rho (x)`$ as

$$F(x)_{\mathrm{}}^x𝑑y\rho (y).$$
(74)
Using definitions presented above, the extended central limit theorem is stated as follows:

###### Theorem .2 (The extended central limit theorem)

In order that $`F`$ belongs to some domain of attraction it is necessary that $`F`$ is of the form

$$1F(x)+F(x)x^\alpha L(x)$$
(75)
for some $`0<\alpha 2`$ as $`x\mathrm{}`$.

* If $`\alpha =2`$ then $`F`$ belongs to the domain of attraction of the normal distribution.
* If $`\alpha <2`$ and

$$\frac{1F(x)}{1F(x)+F(x)}p,\frac{F(x)}{1F(x)+F(x)}q$$
(76)
as $`x\mathrm{}`$ then $`F`$ belongs to the domain of attraction of the stable distribution with the characteristic exponent $`\alpha `$. If Eq. (76) fails then $`F`$ belongs to no domain of attraction.

Here, Eq. (76) only means that the limits exist. If positive random variables are considered, $`p=1`$ and $`q=0`$ since $`F(x)=0`$. The limits $`p`$ and $`q`$ give parameters of a stable distribution. In addition, the centering parameter $`b_n`$ is given as followings;

* If $`\alpha >1`$ then $`b_n`$ is given by the expectation value.
* If $`\alpha =1`$ then $`b_n`$ is given by

$$b_n=n_{\mathrm{}}^{\mathrm{}}𝑑x\mathrm{sin}\left(\frac{x}{a_n}\right)\rho (x).$$
(77)
* If $`\alpha <1`$ then $`b_n=0`$.

The norming constant $`a_n`$ is chosen so that

$$\frac{n}{a_{n}^{}{}_{}{}^{\alpha }}L(a_n)C,$$
(78)
where $`C>0`$ is a constant.

This theorem implies that when the disorder is strong and the probability distribution function of the normalized sum $`S_n/a_nb_n`$ converges in the limit $`n\mathrm{}`$ the limiting probability distribution is the stable distribution with the characteristic exponent $`0<\alpha 1`$ given by Eq. (75). The inequality $`0<\alpha 1`$ is due to divergence of the expectation value (the strong disorder).

For our convenience, we rewrite the necessary condition Eq. (75) for the convergence of the probability distribution of $`S_n/a_nb_n`$ in terms of the probability density $`\rho (x)`$. Since $`\rho (x)=F^{}(x)dF(x)/dx`$, the asymptotic behavior of the probability density is given by

$$\rho (x)+\rho (x)x^{\alpha 1}L(x)[1\frac{1}{\alpha }\frac{xL^{}(x)}{L(x)}].$$
(79)
The second term on the right hand side is evaluated from Eq. (73). By differentiation with respect to $`t`$ of Eq. (73) and setting $`t=1`$, it is shown that

$$\underset{s\mathrm{}}{lim}\frac{L^{}(s)s}{L(s)}=0.$$
(80)
It implies that

$$\rho (x)+\rho (x)x^{\alpha 1}L(x)$$
(81)
as $`x\mathrm{}`$.

For the asymptotic behavior of the stable distribution, the following theorem is known.

###### Theorem .3 (Asymptotic behavior of a stable distribution)

If $`G`$ is a stable distribution obtained as the limiting distribution of the extended central limit theorem, then as $`x\mathrm{}`$

$$x^\alpha [1G(x)]Cp\frac{2\alpha }{\alpha },x^\alpha G(x)Cq\frac{2\alpha }{\alpha },$$
(82)
where the parameters $`C,p,q`$ are given by Eqs. (76) and (78).

It implies that the asymptotic behavior of the probability density $`g(x)=G^{}(x)`$ is given by

$$g(x)Cp(2\alpha )x^{\alpha 1},g(x)Cq(2\alpha )x^{\alpha 1}.$$
(83)
Hence, a stable distribution with the characteristic exponent $`\alpha `$ has absolute moments of all orders $`<\alpha `$ and absolute moments of all orders $`\alpha `$ do not exist.

Next, we present the relation holding when $`0<\alpha <1`$. It tells about the contribution of the maximal term to the sum. Let $`X_1,X_2,\mathrm{}`$ be independent random variables with the common distribution $`F`$ satisfying Eq. (75) with $`0<\alpha <1`$, i.e., belonging to the domain of attraction of the stable distribution with the characteristic exponent $`\alpha `$. Put $`M_n=\mathrm{max}(X_1,\mathrm{},X_n)`$. Then,

$$\frac{S_n}{M_n}\frac{1}{1\alpha }$$
(84)
as $`n\mathrm{}`$. It implies that the maximal term is of the same order as the sum with probability one when $`0<\alpha <1`$.

### .2 The centering constant $`b_n`$ when $`\alpha =1`$

Here, we show that the centering constant $`b_n`$ diverges for the positive random variables $`x`$ satisfying Eq. (75) with $`\alpha =1`$;

$$1F(x)x^1L(x).$$
(85)
Let $`x_0>0`$ denote the lower limit of the distribution and the probability density is given by $`\rho (x)=F^{}(x)`$. From Eq. (78), the norming constant $`a_n`$ is chosen so that

$$\frac{n}{a_n}L(a_n)C,$$
(86)
where $`C>0`$ is a constant. The asymptotic behavior of the probability density $`\rho (x)`$ is derived from Eq. (81);

$$\rho (x)x^2L(x).$$
(87)
We evaluate $`b_n`$ by dividing the region of integration of Eq. (77);

$$b_n=n_{x_0}^{x_n}𝑑x\mathrm{sin}\left(\frac{x}{a_n}\right)\rho (x)+n_{x_n}^{\mathrm{}}𝑑x\mathrm{sin}\left(\frac{x}{a_n}\right)\rho (x)$$
(88)
where $`x_n`$ is a constant satisfying

$$x_na_n,\underset{n\mathrm{}}{lim}x_n=\mathrm{}.$$
(89)
Since $`x_na_n`$, the first term on the right hand side of Eq. (88) is given by

$$\text{the first term}\frac{n}{a_n}_{x_0}^{x_n}𝑑xx\rho (x)>0.$$
(90)
From Eq. (87), the leading term of the second term is evaluated as

the second term $``$ $`n{\displaystyle _{x_n}^{\mathrm{}}}𝑑x\mathrm{sin}\left({\displaystyle \frac{x}{a_n}}\right)x^2L(x)`$ (91)
$`=`$ $`{\displaystyle \frac{n}{a_n}}{\displaystyle _{x_n/a_n}^{\mathrm{}}}𝑑x\mathrm{sin}\left(x\right)x^2L(a_nx),`$
where a numerical constant is ignored. Since $`L`$ is a function slowly varying, the leading term is given by

the second term $``$ $`{\displaystyle \frac{n}{a_n}}L(a_n){\displaystyle _{x_n/a_n}^{\mathrm{}}}𝑑x\mathrm{sin}\left(x\right)x^2`$ (92)
$``$ $`{\displaystyle \frac{n}{a_n}}L(a_n)\mathrm{log}(a_n).`$
From Eq. (86), it is shown that

$$\text{the second term}\mathrm{log}(a_n).$$
(93)
Consequently, since the first term is positive and the second term becomes positive infinity in the limit $`n\mathrm{}`$, the centering constant $`b_n`$ diverges in the limit $`n\mathrm{}`$.


format: text

text:
# Avalanche Dynamics in Wet Granular Materials

## Abstract

We have studied the dynamics of avalanching wet granular media in a rotating drum apparatus. Quantitative measurements of the flow velocity and the granular flux during avalanches allow us to characterize novel avalanche types unique to wet media. We also explore the details of viscoplastic flow (observed at the highest liquid contents) in which there are lasting contacts during flow, leading to coherence across the entire sample. This coherence leads to a velocity independent flow depth at high rotation rates and novel robust pattern formation in the granular surface.

Avalanches and landslides Brabb89 ; Dikau96 ; mudflow  are among the most dramatic of natural catastrophes, and they also provide an evocative metaphor for a wide range of propagating breakdown phenomena from impact ionization in semiconductors Clauss90  to magnetic vortex motion in superconductors Field95 . On the other hand, the existence of avalanches, i.e. the sudden collapse of the system previously frozen into a high energy state, is a fundamental manifestation of the metastable nature of granular materials. Studies of avalanches Rajchenbach90 ; Liu91 ; dryava  and surface flows dryflow  in granular media have largely focused on dry grains. By wetting such media, however, one introduces controllable adhesive forces between the grains Schubert84 ; Craig97  which lead to qualitatively new behavior Hornbaker97 ; Tegzes99 ; Halsey98 ; Bocquet98 ; Fraysse99 ; Nase01 ; Samadani00 . In our previous work Tegzes99  we identified three fundamental regimes for the repose angle of wet granular materials as a function of the liquid content. The granular regime at very low liquid contents is dominated by the motion of individual grains; in the correlated regime corresponding to intermediate liquid contents, a rough surface is formed by the flow of separated clumps; and the repose angle of very wet samples results from cohesive flow with viscoplastic properties.

Here we report investigations of the avalanche dynamics and flow properties of wet granular materials, employing a rotating drum apparatus (a cylindrical chamber partly filled with a granular medium and rotated around a horizontal axis) Rajchenbach90 ; Liu91 . At low rotation rates, the medium remains at rest relative to the drum while its surface angle is slowly increased by rotation, up to a critical angle ($`\theta _{max}`$) where an avalanche occurs, thus decreasing the surface angle to the repose angle ($`\theta _r`$). The flow becomes continuous at high rotation rates, but the transition between avalanching and continuous flow is hysteretic in rotation rate in dry media Rajchenbach90 ; Caponeri95 .

Previous studies of cohesive granular media in a rotating drum Bocquet98 ; Fraysse99 ; Nase01  have focused on the surface angles of the medium before and after avalanches. For example, Quintanilla et al. recently performed a statistical analysis of avalanche size based on these angles, and showed that the average clump size increased with cohesion Viturro01 . In our measurements, we focus instead on characterizing the dynamics of cohesive flow. We quantitatively investigate the flow dynamics during avalanches at different liquid contents by analyzing the time evolution of the averaged surface profile obtained from hundreds of avalanche events, and we also measure surface velocities during continuous flow. In particular, we explore the nature of the viscoplastic flow, which displays unique characteristics associated with coherent motion over the entire granular surface.

We studied glass spheres thoroughly mixed with small quantities of hydrocarbon oil oilprop  varying between $`\tau =`$ $`0.001\%`$ and $`5\%`$ of the void volume. In this regime the flow of oil due to gravity can be neglected. The inner diameter of our drum was $`16.8`$ cm, its width was $`3.2`$ cm, and the granular filling was $`30\%`$. Measurements performed in a thinner ($`2`$ cm) drum verified that wall effects do not modify the qualitative behavior. We studied two different sizes of beads ($`d_1=0.5`$ mm$`\pm 20\%`$ and $`d_2=0.9`$ mm$`\pm 11\%`$), and the results were qualitatively equivalent in the two cases except where noted otherwise. The experiment was recorded by a CCD video camera interfaced to a computer which could analyze the spatiotemporal evolution of the surface profile (height variations in the axial direction were negligible). To determine the position of the surface, we used background illumination and an algorithm based on local contrast which gave us a resolution of $`\delta t=0.03`$ seconds and $`\delta x=0.5`$ mm.

Surface angles. — We first investigated the avalanching – continuous flow transition for various liquid contents. We slowly increased and then decreased the rotation rate and measured the average surface angles, $`\theta _{max}`$ and $`\theta _r`$ in the avalanching regime. We consider the flow continuous when the medium never comes to rest with respect to the drum, which coincides with $`\mathrm{\Delta }\theta \theta _{max}\theta _r<1^{}`$. The results of typical runs are presented in Fig. 1 (a), (b) and (c), where $`\theta _{max}`$ and $`\theta _r`$ are plotted as a function of the rotation rate.

The three basic types of behavior observed previously are also reflected in the present experiment. In the granular regime at low liquid contents, (e.g. $`\tau =0.04\%`$), the behavior is qualitatively similar to the dry case, and we observe a clear hysteretic transition between continuous flow and avalanching (Fig. 1 (a)). At a somewhat higher liquid content ($`\tau =0.125\%`$, see Fig. 1 (b)) the correlated behavior is marked by a lack of hysteresis, and avalanching behavior is observed at a relatively high rotation rate compared to the other regimes. The continuous flow in this regime consists of a stream of separated clumps rather than the constant flux seen in the other two regimes. At $`\tau =0.375\%`$ (Fig. 1 (c)) hysteresis is again observed, reflecting the onset of the viscoplastic continuous flow, which is smooth and coherent over the entire sample.

An interesting feature of the curves in Fig. 1 (a-c), is that $`\theta _{max}`$ decreases with increasing rotation rate, although at slow rotation (where the rotation of the drum during the avalanche is negligible) the rotation rate influences the avalanche process only through changes in the waiting time between successive avalanches. We find that $`\theta _{max}`$ depends logarithmically on this waiting time unpublished , reminiscent of other granular aging phenomena Bocquet98 ; Losert00 . The effect is more pronounced at higher liquid contents, and we attribute it to the motion of liquid flowing towards the contact points rather than condensation effects Bocquet98  since the vapor pressure of the oil is low.

If the rotation rate is sufficiently low, then discrete avalanches occur at all liquid contents. In the following paragraphs we analyze this avalanching regime. Fig. 1 (d) shows $`\theta _{max}`$ (averaged over several hundred avalanches) as a function of liquid content. The distribution is close to Gaussian, and the error bars in the figure indicate the standard deviation. The average avalanche size $`\mathrm{\Delta }\theta =\theta _{max}\theta _r`$ and the width of the avalanche size distribution $`2\sigma _{\mathrm{\Delta }\theta }=2(\mathrm{\Delta }\theta ^2\mathrm{\Delta }\theta ^2)^{1/2}`$ are also plotted in Fig. 1 (d). We interpret these results below in the context of a more detailed investigation of the avalanche dynamics.

Flow Dynamics. — By using the rotating drum apparatus, we can obtain information not just about the medium before and after the avalanche events, but we can also study the details of the grain motion during avalanche events. In order to analyze the dynamics of avalanches we have obtained two dimensional space-time matrices, $`h(x,t)`$, characterizing the sample surfaces throughout the avalanche process which can then be analyzed to produce a variety of information about the individual avalanches. By taking derivatives of the $`h(x,t)`$ data, we obtain the local angle ($`\alpha (x,t)=\mathrm{arctan}[_xh(x,t)`$\]) and the local vertical velocity ($`u(x,t)=_th(x,t)`$) of the surface profile. Furthermore, by then integrating the vertical velocity (using the continuity equation and assuming constant density Douady98 ), we also obtain the local flux in the avalanche, i.e. $`\varphi (x,t)=\rho w_{D/2}^x_th(x^{},t)dx^{}`$ (where $`\rho `$ is the grain density, $`w`$ is the width of the drum, and $`D`$ is the drum diameter), which represents the material flowing through a given vertical plane at position $`x`$.

In Fig. 2 (a) we present snapshots of the progression of single avalanches for several typical liquid contents. Fig. 2 (b) displays $`\alpha (x,t)`$ as a function of space (horizontally) and time (downwards) for the same individual avalanches. In Fig. 2 (c-e), we present the average behavior of 300-500 avalanches at the same liquid contents, and show similar graphs of the time evolution of $`\alpha (x,t)`$, $`u(x,t)`$, and $`\varphi (x,t)`$, where $``$ denotes averaging over avalanches. By obtaining these quantitative measures of the averaged properties, we can separate the robust characteristics of the avalanche dynamics from the large fluctuations which are inherent in avalanche processes.

In the avalanches among dry grains Rajchenbach90 ; Liu91 ; dryava , the avalanches consist of flow in a thin layer near the smooth surface and have a much longer duration and much smaller flux than in the wet media. The size of the avalanches is $`\mathrm{\Delta }\theta 2^{}\pm 0.5^{}`$. At low but nonzero liquid contents (e.g. $`\tau =0.04\%`$) $`\theta _{max}`$ and $`\mathrm{\Delta }\theta `$ become larger due to the onset of intergrain cohesion. In this granular regime the avalanche is initiated by a front traveling downhill (bright line in Fig. 2 (c)), but we can also observe a second front traveling uphill (dark line). The second front is also apparent in the motion of the point of maximum flux, Fig. 2 (e), and it corresponds to the grains reaching a solid barrier at the bottom Makse97 .

At intermediate liquid contents (the correlated regime, $`0.1\%\tau 0.5\%`$), both $`\theta _{max}`$ and $`\mathrm{\Delta }\theta `$ increase substantially. In this regime, the principal failure mechanism is a fracture along a curved slip plane (approximated by the green lines in Fig. 2 (a)), analogous to the dynamics of a class of geological events known as “slides” Dikau96 . In the figure there is a single slip plane, which leads to large avalanches with a relatively narrow distribution. However, at larger liquid contents the avalanches occur through a succession of local slip events followed by a large avalanche extending over the whole surface unpublished . This leads to large fluctuations in the size of the avalanches (see Fig. 1 (d)) as observed in cohesive powders Viturro01 . The medium becomes more cohesive with increasing liquid content, and the upward traveling front disappears for $`\tau 0.2\%`$ since the material moving as a connected block stops coherently when it hits the bottom.

At the highest liquid contents, $`\tau 0.51\%`$, the onset of coherent viscoplastic flow is apparent through numerous qualitative changes in the dynamics of flow. The coherence of avalanche flow in this regime is demonstrated by the parallel lines during the avalanche in Fig 2 (b) and (c). This behavior is qualitatively similar to a different class of geological event, called a ”debris flow” or ”mudflow” Dikau96 ; mudflow . Since the whole surface moves coherently (rather than breaking apart and evolving separately in different parts of the surface), fluctuations are strongly suppressed Tegzes99  and the avalanche size distribution is rather narrow (Fig. 1 (d)). The small decrease in $`\theta _{max}`$ is probably due to lubrication Tegzes99 .

A novel property of the viscoplastic avalanches is the topology of the top surface which spontaneously forms a nearly periodic pattern (seen in Fig. 2 (a) and (b)). This surface structure is maintained essentially intact during the avalanche, indicating that there are lasting contacts in the flowing layer. Moreover, the same pattern in the surface features is reproduced at the end of each avalanche. This is demonstrated most clearly in Fig 2 (c), where the average of 347 avalanches of the $`\tau =5\%`$ sample has the same features as the typical individual avalanche shown in Fig. 2 (a) and (b). The robust nature of the surface structure of the wettest grains is in sharp contrast to the other regimes where averaging completely smoothes out the smaller surface features. We can understand this behavior as resulting from coherence of the entire flow, which strongly reduces fluctuations in this regime. With minimal fluctuations, the final surface structure after each avalanche is essentially the same, thus setting the same initial condition for the next avalanche. With the same initial conditions for each avalanche, naturally the surface features are reproduced each time. Our experiments with the larger ($`d=0.9`$ mm) beads also revealed pattern formation, but with a smaller characteristic size corresponding to 8-10 grain diameters. The difference is probably due to the smaller ratio of the cohesive forces to the gravitational forces on the grains.

In addition to differences in the avalanching behavior, the viscoplastic regime also displays continuous flow which is rather different from that at lower liquid contents. We compare the viscoplastic and granular continuous flows by measuring their surface velocities $`v`$ (by means of tracer particles) as a function of rotation rate in Fig. 3 (continuous flow is only observed in the correlated regime above $`6`$ rpm, and we therefore do not include data from that regime here). For these measurements we used the larger, ($`d=0.9`$ mm) beads. The viscoplastic flow is slower than the granular flow by a factor of two, and the two curves correspond to rather different functional forms. The curvature in the granular case suggests that the flow depth increases with rotation rate as previously suggested finitedepth . The linear $`v(\omega )`$ function of the viscoplastic flow suggests a constant flow depth which is independent of the flow rate. This observation is consistent with the coherent nature of the flow: the flow depth is not determined by local mechanisms, but is fixed by the geometry of the whole system dryflow .

Discussion. — While the flows we observe appear to have analogies with geological events, it is important to note that real geological materials usually consist of polydisperse irregular particles often with very high ($`\tau 100\%`$) water contents and that the scaling of our system to geological lengthscales is non-trivial. Furthermore, avalanche studies in real soil have demonstrated additional phenomena associated with soil saturation Iverson00 . Interestingly we still recover some of the basic dynamical processes in our model system, which constitutes an important step towards the description of qualitatively different flow behaviors in the framework of a single model Nase01 ; Gray01 . The changes in the dynamical behavior with wetting are associated with the increasingly coherent nature of the flow, i.e. the formation of coherently moving clusters – clumps – due to the increased cohesion and viscous effects. Within a cluster, local velocity fluctuations should be suppressed, and thus the local granular temperature ($`T=v^2v^2`$) should approach zero, but the clusters themselves both form and break apart during an avalanche process in a finite container. An important theoretical question raised by our data is how a length scale describing the size of the clumps may emerge from a granular flow model, and how such a length scale would vary with the type of media, the nature of intergranular adhesion, and the type of granular flow.

###### Acknowledgements.

We gratefully acknowledge helpful discussions with J. Banavar, A.-L. Barabási, and Y. K. Tsui. We are also grateful for support from the Petroleum Research Fund and the NASA Microgravity Fluid Physics Program. P. T. and T. V. are grateful for the partial support from OTKA Grant No. T033104.


text:
# Invariant Imbedding Equations for Electromagnetic Waves in Stratified Magnetic Media: Applications to One-Dimensional Photonic Crystals

## Abstract

We derive the invariant imbedding equations for plane electromagnetic waves propagating in stratified magnetic media, where both dielectric and magnetic permeabilities vary in one spatial direction in an arbitrary manner. These equations allow us to obtain the reflection and transmission coefficients of the waves and the field amplitudes inside the media exactly for any polarization and incident angle of the incoming wave by solving an initial value problem of a small number of ordinary differential equations. We apply our results to one-dimensional photonic crystals, where the periodic variations of both dielectric and magnetic permeabilities create photonic band gaps in the frequency spectrum.

The propagation of electromagnetic waves in inhomogeneous media is an important topic in various branches of physics, including optics, plasma physics, astrophysics, and condensed matter physics yeh ; ginz ; chandra ; sheng . Great attention has been paid to the cases where the inhomogeneity is random or periodic yeh ; sheng ; ishi ; krav ; kly1 ; kush ; jkps1 ; jkps2 . The problem of electromagnetic wave propagation in random dielectric media is analogous to the Anderson localization problem of noninteracting electrons in a random potential sheng ; john1 ; john2 ; crete . In recent years, much effort has been made to demonstrate the localization of light experimentally lagen ; gena  and to use the effect for developing the so-called random laser lawa ; cao . Another important area of research has been created from an analogy with the electronic band structure problem in condensed matter physics. In media where the dielectric permeability varies periodically in space, a series of photonic band gaps can appear in their frequency spectra john2 ; yab1 ; yab2 . Because of the promising possibility that a number of devices based on this phenomenon can be used in various optoelectronic applications, research activity in photonic band gap structures has increased explosively in the last decade crete ; mit .

A number of analytical and numerical techniques have been developed for analyzing the wave propagation in inhomogeneous media yeh ; ginz ; chandra ; sheng ; ishi ; krav ; kly1 ; pendry ; jkps3 . In the most general cases where the inhomogeneity is three-dimensional, the exact solution of the problem is extremely hard to obtain, and one is forced to use some approximation scheme. In the cases where the inhomogeneity is one-dimensional, however, there exist several theoretical methods that allow (at least numerically) exact solutions of the problem in some simple situations. The invariant imbedding method is a powerful and most versatile tool for handling one-dimensional situations kly1 ; kly2 ; kly3 ; ram1 ; ram2 ; kim . It can be used to obtain exact solutions for the reflection and transmission coefficients of incoming waves and the electric and magnetic field amplitudes inside arbitrarily inhomogeneous media. When the inhomogeneity is random, it can be used to obtain the exact disorder-averaged reflection and transmission coefficients and field amplitudes kim . The invariant imbedding method can also be used to study the wave propagation in nonlinear dielectric media kly2 ; ram2  and the propagation of several coupled waves kly3  in an exact manner.

In the present letter, we apply the invariant imbedding method to cases where both the dielectric and magnetic permeabilities vary in one direction in space in an arbitrary manner and derive a set of first-order ordinary differential equations called the invariant imbedding equations. To the best of our knowledge, these equations have never been derived before. We solve the initial value problem of the invariant imbedding equations numerically to obtain the exact reflection and transmission coefficients and field amplitudes. We apply our results to the cases where one-dimensional periodic variations of both the dielectric and magnetic permeabilities create photonic band gaps in the frequency spectrum. This situation can be realized in one-dimensional magnetic photonic crystals made of alternating layers of dielectric and ferrite films sou ; kee2 ; kee3 ; kee4 .

We are interested in the propagation of a plane, monochromatic, and linearly-polarized electromagnetic wave of angular frequency $`\omega `$ and vacuum wave number $`k_0=\omega /c`$, where $`c`$ is the speed of light in vacuum. The wave is assumed to be incident from a homogeneous region on a layered or stratified medium, where the dielectric permeability $`ϵ`$ and the magnetic permeability $`\mu `$ (and therefore the refractive index $`n=\sqrt{ϵ\mu }`$ and the wave impedance $`Z=\sqrt{\mu /ϵ}`$) vary only in one direction in space. We take this direction as the $`z`$-axis and assume the inhomogeneous, but isotropic, medium lies in $`0zL`$. Without loss of generality, we assume that the wave propagates in the $`xz`$-plane. Since the medium is uniform in the $`x`$-direction, the dependence on $`x`$ can be taken as being through a factor $`e^{iqx}`$, with $`q`$ being a constant. When $`q=0`$, the wave passes through the medium normally. If $`q0`$, the wave is said to pass obliquely.

For $`q0`$, we have to distinguish two independent cases of polarization. In the first case, the electric field vector is perpendicular to the $`xz`$-plane. Then one can show easily from Maxwell’s equations that the complex amplitude of the electric field, $`E=E(z)`$, satisfies

$$\frac{d^2E}{dz^2}\frac{1}{\mu (z)}\frac{d\mu }{dz}\frac{dE}{dz}+\left[k_0^2ϵ(z)\mu (z)q^2\right]E=0.$$
(1)
This type of wave is known as an $`s`$ (or TE) wave. In the other case, the magnetic field vector is perpendicular to the $`xz`$-plane. Then the complex amplitude of the magnetic field, $`H=H(z)`$, satisfies

$$\frac{d^2H}{dz^2}\frac{1}{ϵ(z)}\frac{dϵ}{dz}\frac{dH}{dz}+\left[k_0^2ϵ(z)\mu (z)q^2\right]H=0.$$
(2)
This type of wave is called $`p`$ (or TM) wave. We note that there is a simple symmetry between Eqs. (1) and (2): Eq. (2) is obtained from Eq. (1) by exchanging $`ϵ`$ and $`\mu `$ and replacing $`E`$ by $`H`$. From now on, we will consider only the $`s`$ wave. The results for the $`p`$ wave follow from those for the $`s`$ wave in a trivial manner.

We assume that the wave is incident from the region where $`z>L`$ and is transmitted to the region where $`z<0`$. The dielectric and magnetic permeabilities are assumed to be given by

$`ϵ(z)`$ $`=`$ $`\{\begin{array}{cc}ϵ_1\hfill & \text{if }z>L\hfill \\ ϵ_R(z)+iϵ_I(z)\hfill & \text{if }0zL\hfill \\ ϵ_2\hfill & \text{if }z<0\hfill \end{array},`$ (6)
$`\mu (z)`$ $`=`$ $`\{\begin{array}{cc}\mu _1\hfill & \text{if }z>L\hfill \\ \mu _R(z)+i\mu _I(z)\hfill & \text{if }0zL\hfill \\ \mu _2\hfill & \text{if }z<0\hfill \end{array},`$ (10)
where $`ϵ_1`$, $`ϵ_2`$, $`\mu _1`$, and $`\mu _2`$ are real constants and $`ϵ_R(z)`$, $`ϵ_I(z)`$, $`\mu _R(z)`$, and $`\mu _I(z)`$ are arbitrary real functions of $`z`$. When $`\theta `$ is defined as the angle of incidence, the constant $`q`$ in Eqs. (1) and (2) is equal to $`\sqrt{ϵ_1\mu _1}k_0\mathrm{sin}\theta `$.

We consider a plane wave of unit magnitude $`\stackrel{~}{E}(x,z)=E(z)e^{iqx}=e^{ip(Lz)+iqx}`$, where $`p=\sqrt{ϵ_1\mu _1}k_0\mathrm{cos}\theta `$, incident on the medium from the right. The quantities of main interest are the complex reflection and transmission coefficients, $`r=r(L)`$ and $`t=t(L)`$, defined by the wave functions outside the medium:

$$\stackrel{~}{E}(x,z)=\{\begin{array}{cc}e^{ip(Lz)+iqx}+r(L)e^{ip(zL)+iqx},\hfill & z>L\hfill \\ t(L)e^{ip^{}z+iqx},\hfill & z<0\hfill \end{array},$$
(11)
where $`p^{}=\sqrt{ϵ_2\mu _2}k_0\mathrm{cos}\theta ^{}`$ and $`\theta ^{}`$ is the angle that outgoing waves make with the negative $`z`$-axis.

We generalize the invariant imbedding method developed in Ref. 23 to the cases where both $`ϵ`$ and $`\mu `$ are inhomogeneous. Let us consider the $`E`$ field in Eq. (1) as a function of both $`z`$ and $`L`$: $`E=E(z;L)`$. Starting from the crucial observation that the boundary value problem of the wave equation, Eq. (1), can be transformed into an integral equation

$`E(z;L)=g(z;L)+{\displaystyle \frac{ip}{2}}{\displaystyle _0^L}𝑑z^{}g(z;z^{})`$
$`\times \left\{{\displaystyle \frac{ϵ(z^{})}{ϵ_1}}{\displaystyle \frac{\mu (z^{})}{\mu _1}}+{\displaystyle \frac{q^2}{p^2}}\left[{\displaystyle \frac{ϵ(z^{})}{ϵ_1}}{\displaystyle \frac{\mu _1}{\mu (z^{})}}\right]\right\}E(z^{};L),`$
$`g(z;z^{})=\mathrm{exp}\left[ip\mathrm{s}gn(zz^{}){\displaystyle _z^{}^z}𝑑z^{\prime \prime }{\displaystyle \frac{\mu (z^{\prime \prime })}{\mu _1}}\right],`$ (12)
we derive

$$\frac{E(z;L)}{L}=a(L)E(z;L),$$
(13)
where

$`a(L)`$ $`=`$ $`ip{\displaystyle \frac{\mu (L)}{\mu _1}}+{\displaystyle \frac{ip}{2}}\{{\displaystyle \frac{ϵ(L)}{ϵ_1}}{\displaystyle \frac{\mu (L)}{\mu _1}}`$ (14)
$`+{\displaystyle \frac{q^2}{p^2}}[{\displaystyle \frac{ϵ(L)}{ϵ_1}}{\displaystyle \frac{\mu _1}{\mu (L)}}]\}E(L;L).`$
Using this equation, we obtain exact differential equations satisfied by $`r`$ and $`t`$:

$`{\displaystyle \frac{dr(L)}{dL}}=2i\sqrt{ϵ_1\mu _1}k_0\mathrm{cos}\theta {\displaystyle \frac{\mu _R(L)+i\mu _I(L)}{\mu _1}}r(L)`$
$`+{\displaystyle \frac{i}{2}}\sqrt{ϵ_1\mu _1}k_0\mathrm{cos}\theta [{\displaystyle \frac{ϵ_R(L)+iϵ_I(L)}{ϵ_1}}{\displaystyle \frac{\mu _R(L)+i\mu _I(L)}{\mu _1}}`$
$`+{\displaystyle \frac{ϵ_R(L)+iϵ_I(L)}{ϵ_1}}\mathrm{tan}^2\theta {\displaystyle \frac{\mu _1}{\mu _R(L)+i\mu _I(L)}}\mathrm{tan}^2\theta ]`$
$`\times \left[1+r(L)\right]^2,`$
$`{\displaystyle \frac{dt(L)}{dL}}=i\sqrt{ϵ_1\mu _1}k_0\mathrm{cos}\theta {\displaystyle \frac{\mu _R(L)+i\mu _I(L)}{\mu _1}}t(L)`$
$`+{\displaystyle \frac{i}{2}}\sqrt{ϵ_1\mu _1}k_0\mathrm{cos}\theta [{\displaystyle \frac{ϵ_R(L)+iϵ_I(L)}{ϵ_1}}{\displaystyle \frac{\mu _R(L)+i\mu _I(L)}{\mu _1}}`$
$`+{\displaystyle \frac{ϵ_R(L)+iϵ_I(L)}{ϵ_1}}\mathrm{tan}^2\theta {\displaystyle \frac{\mu _1}{\mu _R(L)+i\mu _I(L)}}\mathrm{tan}^2\theta ]`$
$`\times \left[1+r(L)\right]t(L).`$ (15)
These invariant imbedding equations are supplemented with the initial conditions for $`r`$ and $`t`$, which are obtained using the well-known Fresnel formulas:

$`r(0)`$ $`=`$ $`{\displaystyle \frac{\mu _2\sqrt{ϵ_1\mu _1}\mathrm{cos}\theta \mu _1\sqrt{ϵ_2\mu _2ϵ_1\mu _1\mathrm{sin}^2\theta }}{\mu _2\sqrt{ϵ_1\mu _1}\mathrm{cos}\theta +\mu _1\sqrt{ϵ_2\mu _2ϵ_1\mu _1\mathrm{sin}^2\theta }}},`$
$`t(0)`$ $`=`$ $`{\displaystyle \frac{2\mu _2\sqrt{ϵ_1\mu _1}\mathrm{cos}\theta }{\mu _2\sqrt{ϵ_1\mu _1}\mathrm{cos}\theta +\mu _1\sqrt{ϵ_2\mu _2ϵ_1\mu _1\mathrm{sin}^2\theta }}}.`$ (16)
For given values of $`ϵ_1`$, $`ϵ_2`$, $`\mu _1`$, $`\mu _2`$, $`k_0`$ (or $`\omega `$), and $`\theta `$ and for arbitrary functions $`ϵ_R(L)`$, $`ϵ_I(L)`$, $`\mu _R(L)`$, and $`\mu _I(L)`$, we solve the nonlinear ordinary differential equations in Eq. (15) numerically, using the initial conditions in Eq. (16), and obtain the reflection and transmission coefficients $`r`$ and $`t`$ as functions of $`L`$. The reflectivity $``$ and the transmissivity $`𝒯`$ are given by

$$=|r|^2,𝒯=\frac{\mu _1\sqrt{ϵ_2\mu _2ϵ_1\mu _1\mathrm{sin}^2\theta }}{\mu _2\sqrt{ϵ_1\mu _1}\mathrm{cos}\theta }|t|^2.$$
(17)
If the permeabilities have no imaginary parts, the quantities $``$ and $`𝒯`$ satisfy the law of conservation of energy, $`+𝒯=1`$.

The invariant imbedding method can also be used in calculating the field amplitude $`E(z)`$ inside the inhomogeneous medium. We rewrite Eq. (13) as

$`{\displaystyle \frac{E(z;l)}{l}}=i\sqrt{ϵ_1\mu _1}k_0\mathrm{cos}\theta \{{\displaystyle \frac{\mu _R(l)+i\mu _I(l)}{\mu _1}}`$
$`+{\displaystyle \frac{1}{2}}[{\displaystyle \frac{ϵ_R(l)+iϵ_I(l)}{ϵ_1}}{\displaystyle \frac{\mu _R(l)+i\mu _I(l)}{\mu _1}}`$
$`+{\displaystyle \frac{ϵ_R(l)+iϵ_I(l)}{ϵ_1}}\mathrm{tan}^2\theta {\displaystyle \frac{\mu _1}{\mu _R(l)+i\mu _I(l)}}\mathrm{tan}^2\theta ]`$
$`\times [1+r(l)]\}E(z;l).`$ (18)
For a given $`z`$ ($`0<z<L`$), the field amplitude $`E(z;L)`$ is obtained by integrating this equation from $`l=z`$ to $`l=L`$ using the initial condition $`E(z;z)=1+r(z)`$.

We demonstrate the validity and utility of our invariant imbedding equations by applying them to simple one-dimensional photonic crystals, where both $`ϵ`$ and $`\mu `$ vary periodically in one direction in space. More specifically, we consider periodic arrays of alternating layers, characterized by the dielectric permeabilities $`ϵ_a`$ and $`ϵ_b`$, the magnetic permeabilities $`\mu _a`$ and $`\mu _b`$, and the thicknesses $`d_a`$ and $`d_b`$. We assume that the imaginary parts of the permeabilities are zero and the ambient medium is vacuum. In the limit where the length of the system, $`L`$, diverges to infinity, one can find the exact positions of the photonic band gaps by computing the roots of the analytical expression kee4

$$\mathrm{cos}\beta _a\mathrm{cos}\beta _b\frac{1}{2}\left(\frac{h_a}{h_b}+\frac{h_b}{h_a}\right)\mathrm{sin}\beta _a\mathrm{sin}\beta _b+1=0,$$
(19)
where $`n_i=\sqrt{ϵ_i\mu _i}`$, $`Z_i=\sqrt{\mu _i/ϵ_i}`$ ($`i=a,b`$) and

$`h_i=\{\begin{array}{cc}\sqrt{n_{i}^{}{}_{}{}^{2}\mathrm{sin}^2\theta }/(n_iZ_i)\hfill & \text{for the }s\text{ wave}\hfill \\ Z_i\sqrt{n_{i}^{}{}_{}{}^{2}\mathrm{sin}^2\theta }/n_i\hfill & \text{for the }p\text{ wave}\hfill \end{array},`$ (22)
$`\beta _i=k_0d_i\sqrt{n_{i}^{}{}_{}{}^{2}\mathrm{sin}^2\theta }.`$ (23)
In Fig. 1, we plot the transmittance $`T`$ ($`|t|^2`$) of $`s`$ and $`p`$ waves impinging on a 100-period array with $`ϵ_a=4`$, $`\mu _a=2`$, $`ϵ_b=2`$, $`\mu _b=1`$, and $`d_a=d_b`$ as a function of $`k_0\mathrm{\Lambda }`$ where $`\mathrm{\Lambda }=d_a+d_b`$. In this case, the refractive index alternates periodically between $`n_a=2\sqrt{2}`$ and $`n_b=\sqrt{2}`$, whereas the wave impedance is uniform. We find that there is no photonic band gap when the wave is normally incident. For obliquely incident waves, small gaps are created. Furthermore, we observe that the band gaps for the $`s`$ wave agree precisely with those for the $`p`$ wave in their sizes and positions. When $`\theta =45^{}`$ and $`L\mathrm{}`$, we find from Eq. (19) that the lowest band gap exists for $`1.54<k_0\mathrm{\Lambda }<1.63`$, which is fully consistent with our numerical result obtained using the invariant imbedding method.

In Fig. 2, we plot $`T`$ for the case of a 100-period array with $`ϵ_a=4`$, $`\mu _a=1`$, $`ϵ_b=2`$, $`\mu _b=2`$, and $`d_a=d_b`$. In this case, the wave impedance alternates periodically between $`Z_a=1/2`$ and $`Z_b=1`$, while the refractive index is uniform. We find that large band gaps are formed for all incident angles. Similarly to the previous case, the band gaps for the $`s`$ wave agree with those for the $`p`$ wave in their sizes and positions. When $`\theta =0^{}`$ and $`45^{}`$, the first gaps are predicted to exist for $`1.23<k_0\mathrm{\Lambda }<1.91`$ and for $`1.32<k_0\mathrm{\Lambda }<2.04`$, respectively, in agreement with our result.

From the two examples discussed above, we conclude that the wave impedance, not the refractive index, is a crucial parameter mainly responsible for the creation of photonic band gaps, as has been extensively discussed in Refs. 29 to 31. The equality of the $`s`$ and $`p`$ wave gaps is a special feature that works only when either the refractive index or the wave impedance is uniform. It is easy to verify this by examining Eq. (19).

It is straightforward to apply our equations to more realistic cases where the imaginary parts of the permeabilities are nonzero. They can also be applied to cases where the parameters contain randomly-varying components. This versatility of the invariant imbedding method makes it a very useful and convenient tool in designing practical optical devices for various application purposes.

###### Acknowledgements.

This work has been supported by the Korea Research Foundation through Grant No. KRF-2001-015-DP0167.