Flooding dynamics of diffusive dispersion in a random potential

We discuss the combined effects of overdamped motion in a quenched random potential and diffusion, in one dimension, in the limit where the diffusion coefficient is small. Our analysis considers the statistics of the mean first-passage time T(x) to reach position x, arising from different realisations of the random potential. Specifically, we contrast the median T¯(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{T}}(x)$$\end{document}, which is an informative description of the typical course of the motion, with the expectation value ⟨T(x)⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle T(x)\rangle $$\end{document}, which is dominated by rare events where there is an exceptionally high barrier to diffusion. We show that at relatively short times the median T¯(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{T}}(x)$$\end{document} is explained by a ‘flooding’ model, where T(x) is predominantly determined by the highest barriers which are encountered before reaching position x. These highest barriers are quantified using methods of extreme value statistics.


Introduction
There are many situations where particles move under the combined influence of thermal diffusion and a static (or quenched) random potential [1].The particles might be electrons, holes or excitons diffusing in a disordered metallic or semiconductor sample [2], or molecules diffusing in a complex environment such as the cytoplasm of a eukaryotic cell [3].The state of knowledge of this problem is surprisingly under-developed, and in this work we present new results on the simplest version of this problem, in one dimension, where the equation of motion is Here V (x) is a random potential, D is the diffusion coefficient, and η(t) is a white noise signal with statistics defined by ( • denotes expectation value throughout).We assume that V (x) is a smooth random function, defined by its statistical properties, which are stationary in x, and independent of the temporal white noise η(t).The one and two-point statistics of this potential are The correlation function C(∆x) is assumed to decay rapidly as |∆x| → ∞.We also assume that the tails of the distribution of V are characterised by a large-deviation 'rate' (or 'entropy') function J(V ), so that when |V | is large, the probability density function of V is approximated by where throughout we shall use P X to denote the probability density function (PDF) of a random variable X.If P V is a Gaussian distribution, then the entropy function is quadratic, J(V ) ∼ V 2 /2C(0).
It has been proposed that the behaviour of this system is diffusive, with an effective diffusion coefficient which vanishes very rapidly as D → 0: Zwanzig [4] gave an elegant argument which implies that, when V (x) has a Gaussian distribution, the effective diffusion coefficient is ( An earlier work by De Gennes [5] proposes a similar expression.We discuss the origin of this result, and present a generalisation of it to non-Gaussian distributions, in section 2. When D is small, this estimate for the diffusion coefficient depends upon rare events where the potential is unusually large, and it is very difficult to verify equation ( 5) numerically.In addition, numerical experiments show that the model exhibits sub-diffusive behaviour and it has been suggested that there is anomalous diffusion, in the sense that x 2 ∼ t α , with 0 < α < 1 [6,7,8].It is desirable to achieve an analytical understanding of the sub-diffusive behaviour which is observed in numerical simulations of (1).
We should mention that there are also exact results [9,10,11] on a closely related model (motion in a quenched velocity field, which is not the derivative of a potential with a well-defined probability distribution) showing that x 2 ∼ (ln t) 1/4 .This 'Sinai diffusion' process is fundamentally different, because the particle becomes trapped in successively deeper minima of the potential, from which it takes ever increasing time intervals to escape.
We will argue that, while equation ( 5) and its generalisation to non-Gaussian distributions describes the long-time asymptote of the dispersion of particles, the diffusive behaviour only emerges at very long times.At intermediate times, the dynamics of typical realisations is not diffusive.We show that it is determined by the time taken to diffuse across the largest potential barrier which must be traversed to reach position x.The diffusion process is able to traverse a barrier of height ∆V after a characteristic time T ∼ exp(∆V /D) [12], and as time increases the height of the barriers which can be breached, leading to 'flooding' of the region beyond, increases.According to this picture, the dispersion distance x is determined by a problem in extreme value statistics: how large must x be before we reach a barrier of height ∆V ≈ D ln t?By considering the solution of this problem in extreme value statistics, we argue that the median T (with respect to different realisations of the potential V (x)) of the mean-first-passage-time (averaged over η(t)) satisfies ln where x is a lengthscale which characterises the typical distance between extrema of the potential.In the case where the potential has a Gaussian distribution, this implies that the dispersion is sub-diffusive, satisfying ln which is quite distinct from the usual anomalous diffusion behaviour, characterised by power-laws such as x 2 ∼ t α .After a sufficiently long time, the dynamics becomes diffusive, with a diffusion coefficient given by (5).Our arguments will depend upon making estimates of sums of the form where f j are independent identically distributed (i.i.d.) random variables, and ǫ is a small parameter, which we identify with the diffusion coefficient D. We term the S N 'extreme-weighted sums', because the largest values of f j make a dominant contribution to S N as ǫ → 0. In section 2 we show how the mean-first-passage time is related to sums like (8), and in section 3 we analyse some of their statistics, which are used in section 4 to justify our principle result, equation (6).Section 5 describes our numerical investigations, and section 6 is a summary.

The mean first passage time
Our discussion of the dynamics of (1) will focus on the mean first passage problem: what is the mean time T (x) at which a particle released from the origin reaches position x.First passage problems are discussed comprehensively in the book by Redner [13].The result that we require can be found in multiple sources: [14] is the earliest reference that we are aware of and the key formula, equation (9) below, was already applied to equation (1) by Zwanzig [4].
In this section we first quote the general formula for the mean first passage time T (x), as a functional of the potential V (x).If particles are released at x 0 = 0, the mean first passage time to reach position x is given by where the averaging is with respect to realisations of the noise η(t) in the equation of motion (1), with V (x) frozen, so that T (x) is a functional of V (x).
We then (subsection 2.1) discuss the result obtained by Zwanzig [4] for the expectation value T (x) (averaged with respect to realisations of V (x)).Zwanzig gave the result for a potential with Gaussian fluctuations, which we extend to the case of a general form for the large-deviation entropy function (as defined by equation ( 4)).The result obtained by Zwanzig suggests that the dispersion is diffusive, with a diffusion coefficient D eff which vanishes in a highly singular fashion as D → 0. We shall argue that this result is a consequence of the expectation value of T (x) being dominated by very rare large excursions of the potential V (x), and that for typical realisations of the potential the dispersion is much more rapid than the value of T (x) suggests.This requires a more delicate analysis of the structure of the integrals in the expression for T (x), equation (9).In subsection 2.2, we discuss how these integrals may be approximated by sums, analogous to (8), in the limit as D → 0.

Expression for expectation value of mean first passage time
We can make an additional average of (9), with respect to different realisations of the potential, which leads to where in the second line we consider the leading order behaviour as x → ∞.If the motion were simple diffusion, with V = 0, equation ( 9) would evaluate immediately to T = x 2 /2D, so that it is reasonable to identify x 2 /2 T as the effective diffusion coefficient.Hence, assuming that the PDF of V (x) is symmetric between V and −V , we have where V * is the stationary point of the exponent, satisfying From this we obtain In the Gaussian case, where equation ( 14) agrees with (5).

Summation approximations
In order to understand the implications of equation ( 9), we should consider the behaviour of the integral in the limit as D → 0. When D is small this quantity may be estimated from the minima of the potential: where V − j are the values of the N minima between 0 and x, occurring at positions x − j , and where we have defined Note that and consider how to estimate T (x) in the limit as D → 0. Note that S(y) is determined by the values of the minima of V (y) in the interval [0, y], jumping by an amount exp occurring at positions x + j , then T (x) jumps at local maxima.The evolution of S(x) and T (x) are therefore determined by a pair of coupled maps: where we have defined again . These equations are difficult to analyse in the general case, but in the next section we discuss an approach which can be used to treat the limit where D is small.

Statistics of extreme-weighted sums
We have seen that when D is small the integrals defining the mean first passage time are approximated by sums over extrema of the potential, as described by equation ( 17).Accordingly, we study properties of random sums of the form (8) where ǫ is a small parameter and where the f j are drawn from a distribution for which the probability for f j being greater than f is Q(f ).In the case where f has a Gaussian distribution, ( 8) is a sum of log-normal distributed random variables.There is some earlier literature on these sums which shows very little overlap with our results, see [15] and references therein, also [16], which discusses a phase transition which arises in a limiting case.We also consider sums of the form where g j are drawn from the same i.i.d.distribution as the f j .This is a model for the summation which approximates the integral T (x) defined by equation ( 19).When ǫ is sufficiently small, these sums are determined by the largest values of f j and g j , and for this reason we shall refer to S N and T N as extreme-weighted sums.We write the distribution function for f in the form where J (f) is a large deviation rate function.We are interested in the asymptotic behaviour of statistics of the sums S N and T N for small ǫ and large N .The sums vary wildly in magnitude and the mean is dominated by the tail of the distribution of f .Unless N is sufficiently large, values of f j which determine the mean are unlikely to be sampled.This suggests that it will be useful to characterise the distribution of the S N by the median, rather than the mean.We denote the median of X by X and its expectation by X .

Estimate of median of S N
The sum S N may be well approximated by its largest term, which is where f is the largest of the N realisations, f j , with index j = ĵ.We write where If F is close to unity, we can estimate SN by sN .Let us first estimate sN and return to consider F later.Note that where f is the median of the largest value of N samples from the distribution of f .This is determined by setting the probability for N samples to be less than f to be equal to one half: When N ≫ 1, this is determined by the tails of the distribution, where Q(f ) is approximated using (22): An important special case is where the f have a Gaussian distribution, so that in the case where f = 0 and f 2 = 1, implying that In the limit where N is extremely large, we can approximate f by and consequently the median of sN is approximated by Next consider how to estimate the quantity F in equation ( 24), when ǫ ≪ 1.
When N ≫ 1, either F is close to unity or else it is the sum of a large number of terms which make a comparable contribution.The value of F depends upon f.
The f j which contribute to F are i.i.d.random variables, each with a PDF which is the same as that of the general f j , except that there is an upper cutoff at f: the adjustment of the normalisation due to the loss of the tail, f > f , can be neglected.
If the PDF of f is then the expectation value of F is obtained as follows where Noting that we have Hence, we obtain a rather simple approximation for S N , depending upon the extreme value f of the sample of N realisations of the f j : The median of S N is therefore approximated by where f is the solution of equation (32).

Interpretation and generalisation to T N
We shall see that equation (41) gives a quite precise approximation for the median SN , but it is not immediately clear when either of the two terms is dominant.In order to clarify the structure of equation ( 41), we consider an approximate form of the equation determining f, and transform to logarithmic variables.As well as leading to a transparent understanding of equation ( 41), this facilitates making an estimate for TN in the limit where ǫ ≪ 1 and N ≫ 1.We define Note that η and τ are logarithmic measures of, respectively, distance and time, so that a plot of η versus τ gives information about the dispersion due to the dynamics.
Let us consider the limit where the first term in (41), exp[ f/ǫ], is dominant.Note that the condition (32), determining the extreme value of a set of N samples, can be approximated by the requirement that the PDF of f is approximately equal to 1/N.That is 1 ∼ N exp[−J( f)].For the purposes of considering the N → ∞ and ǫ → 0 limit, we can therefore approximate f by a solution of the equation If the second term in equation ( 41) is negligible, as might be expected when f − f * ≪ 1, equations ( 42) and ( 43) then yield a simple implicit equation for σ: If f − f * ≫ 1, and if the second term in (41) is dominant, then SN ∼ S N = N exp(f/ǫ) , and using the Laplace principle we find Note the (45) indicates that dσ dη = 1.Let us compare this with the value of dσ dη obtained from (44), which predicts dη dσ = ǫJ ′ (f).The approximation (44) therefore becomes sub-dominant when 1 = ǫJ ′ (f), which is precisely the equation for f * , equation (37).If we define η * and σ * by writing then assembling these results and definitions, the relationship between η and σ can be summarised in the following equation (47) Note that η(σ) and its first derivative are continuous functions.In the foregoing we defined x as the median of x, but it should be noted that our arguments will lead to equations ( 45) and (46) as N → ∞ if SN denotes any fixed percentile of S N .
Thus far we have considered the behaviour of η as a function of σ rather than of τ , but it is the function η(τ ) which describes the dynamics of the dispersion.Consider the form of the sum T N defined in equation ( 22).When σ < σ * , the value of Sn is almost always determined by f the largest value of f j , and similarly, one the factors exp(g j /ǫ) corresponding to ĝ, the largest of the g j , will predominate over the others.In one half of realisations, those where k(N) > ĵ(N), the largest value of f j contributes to the sum which is multiplied by exp(ĝ/ǫ), and we have T N ∼ exp(ĝ/ǫ) exp( f/ǫ).In cases where ĵ > k, T N is expected to be small in comparison to this estimate.Noting that exp( f/ǫ) and exp(ĝ/ǫ) are independent and both have probability one half to exceed exp( f/ǫ) and exp( ḡ/ǫ) respectively, there are one quarter of realisations where exp[( f+ ĝ)/ǫ] exceeds exp[( f + ḡ)/ǫ] and in half of these realisations T N ≪ exp[( f + ĝ)/ǫ].If we now use the overbar to represent the upper octile of the distribution of T N , rather than the median, we have Using the assumption that the f j and g j have the same PDF, we can conclude that TN ∼ S2 N and hence that τ = 2σ.The equation describing the dispersion as a function of time is therefore where τ * is determined by the condition that dη/dτ = 1 2 when τ = τ * .When τ > τ * , we have TN (50) Equations ( 49) and ( 50) are a description of the logarithm of the typical dispersion η as a function of the logarithm of the time, τ .Usually the function J(V ) has a quadratic behaviour for small values of V , so that the initial dispersion, described by (49), is sub-diffusive.The factor of one half in (50) indicates that the longtime limit is diffusive.Writing η 2 ∼ τ + ln D eff , we see that the effective diffusion coefficient is which is consistent with (14).

Flooding dynamics model for dispersion
In section 2, we showed that the integrals which are used to compute the meanfirst-passage time may be approximated by sums when D is small.In section 3, we considered the statistics of these sums, S N and T N , defined by equations ( 8) and ( 21) respectively.In terms of the calculation discussed in section 2, our estimate of TN corresponds, for N < N * , to the value of T (x) being determined by the difference between the lowest minimum of the potential and its highest maximum, provided the minimum occurs before the maximum.We can therefore think of T (x) being determined by a 'flooding' model, according to which the probability density for locating the particle occupies a region which is constrained by a potential barrier which can trap a particle for time T .As T increases, higher barriers are required.
In terms of the original problem, discussed in section 2, N is the number of extrema of the potential before we reach position x.The arguments of section 3 imply that the upper octile of the mean-first-passage time, T (x), satisfies an equation similar to (48).We define logarithmic variables where x is the mean separation of minima of V (x).In terms of these logarithmic variables, the dispersion is described by which is valid up to τ * , which is defined by the condition Equation ( 53) is our principal result.It applies to any percentile of the distribution which remains fixed when we take the limits N → ∞ and ǫ → 0. When τ is large compared to τ * , equation ( 53) is replaced by a linear relation, with an effective diffusion coefficient An important example is the case where V has a Gaussian distribution, so that J ∼ V 2 /2C(0).In terms of the diffusion coefficient D, equations ( 53)-(55) give and using (55) we find D eff ∼ D exp[−C(0)/D 2 ], in agreement with (5).A sketch of the dependence of η upon τ for the Gaussian case is shown in Fig. 4.

Numerical studies
We performed a variety of numerical investigations, using Gaussian distributed random variables f j to test the theory of extreme-weighted sums, and a Gaussian random function V (x) to test the analysis of continuous potentials.In both cases the Gaussian variables had zero mean and unit variance.In the case of the random potential, we also used a Gaussian for the correlation function, with a correlation length of order unity: Fig. 1 The dynamics of a typical realisation, characterised by the median T of the mean-firstpassage time, shows a crossover from sub-diffusive 'flooding' dynamics to slow diffusion.

Discrete sums
We characterised the statistics of the discrete sum ( 8) by making a careful estimate of its median, equation (41).In order to evaluate equation ( 41), we need a solution of the implicit equation (32), which determines f.By substituting (33) into (32), we find f The expression for the median approaches that for the mean at large values of N when f exceeds f * = 1/ǫ.For very large N and very small ǫ, the medians of S N and T N are estimated by simplified expressions, relating σ = ln SN and τ = ln TN to η = ln N .In the Gaussian case, these equations (47), (49) and (50) give where These equations imply that, in the limit as ǫ → 0, if we plot y = σ/η * as function of x = η/η * , the numerical data for SN should collapse onto the function Similarly, y ′ = τ /η * plotted as a function of x = η/η * should collapse to y ′ = 2f(x).
Figure 2 plots ln SN , and ln S N as a function of η = ln N , for different sample sizes, for ǫ = 1/4 (a) and ǫ = 1/6 (b).We compare with the theoretical prediction, obtained from (41) and (58) (for the median) and (59) (for the mean).The agreement with the theory for the median is excellent.Note that the convergence of the mean value for different sample sizes is very poor when η < η * = 1/2ǫ 2 (this is especially apparent for smaller values of ǫ).In figure 3, we plot y = σ/σ * (a) and y ′ = τ /τ * (b) as a function of x = η/η * , for all of the values of ǫ in our data set, using the largest sample size (M = 1000) in each case, comparing with the theoretical scaling function (63).We see convergence towards the function (63) as ǫ → 0. In panels (c) and (d), we make the same comparison using the upper octile rather than the median.

Continuous potentials
In the case of a continuous potential, we require the mean separation of maxima or minima, x, in order to make a comparison with theory.The density D of zeros of V ′ (x) may be determined by the approach developed by Kac [17] and Rice [18].If P(V, V ′ , V ′′ ) is the joint PDF of V (x) and its first two derivatives, evaluated at the same point, we find that By noting that the vector (V, V ′ , V ′′ ) has a multivariate Gaussian distribution, and expressing P(V, V ′ , V ′′ ) in terms of the correlation function of the elements of this vector, we obtain D and hence the separation of minima x for the potential satisfying (57): First we investigated whether the mean-first passage time can be accurately represented by sums over maxima and minima of the potential.In figure 4, we compare the numerical evaluation of the integrals S(x) (a) and T (x) (b), given by Eq. ( 16) and Eq. ( 19), respectively, with the approximations which estimate the integrals using maxima and minima, Eq. (20).We evaluated the median T (x) and mean T (x) of the mean first passage time T (x) for 1000 realisations of the potential V (x), up to xmax = 10 4 , for D ∈ {1/3, 1/4, 1/5, 1/6, 1/7}.According to the discussion in section 4, we expect that τ = ln T (x) and η = ln(x/x) are related by η = J(Dτ /2) = D 2 τ 2 /8, up to a maximum value of η, given by η * = 1/2D 2 .In figure 5 we plot Y med = 2D 2 ln T (x) and Yav = 2D 2 ln T (x) as a function of X = 2D 2 ln(x/x) for different values of D, and compare with the theoretical scaling function, given by equation (63).

Conclusions
In his analysis of equation ( 1), Zwanzig considered the mean-first-passsage time T (x) to reach displacement x.Computing the expectation value T (x) over different realisations of the random potential, he showed [4] that T (x) ∼ x 2 , which is consistent with a diffusive dispersion, with an effective diffusion coefficient D eff .The effective diffusion coefficient vanishes in a highly singular manner as D → 0, and numerical studies have suggested that equation (1) exhibits anomalous diffusion [6,7,8].It seems evident that the discrepancy between these two pictures of the dynamics results from the expectation value T (x) being dominated by rare events, where an unusually large fluctuation of the potential V (x) acts as a barrier to dispersion.The central limit theorem is applicable to this problem, and at sufficiently large values of x the ratio T (x)/ T (x) is expected to approach unity, for almost all realisations of V (x).However, at values of x which are of practical relevance, most realisations T (x) will be much smaller than T (x) .
In order to give a description of the dynamics of (1) which is both empirically useable and analytically tractable, we considered the median (with respect to different realisations of the potential) of the mean-first-passage time.In the limit where the diffusion coefficient D is small, the integrals which appear in the expression for the first passage time, equation (9), are dominated by maxima and minima of the potential, described by equations ( 17) and (20).This observation led us to consider the statistics of sums of exponentials of random variables, equations ( 8) and (21).We gave a quite precise estimate, equation (41), for the median of ( 8) and also derived simple relations describing the asymptotic behaviour of these sums, equations (47) and (49).
It is these expressions which enable us to formulate a concise asymptotic description of the dynamics of (1) in the limit as D → 0, in terms of the large deviation rate function of the potential, J(V ).We argued that at very long length scales T (x) approaches the expectation value T (x) , and that the dispersion is diffusive, in accord with the theory of Zwanzig [4].On shorter timescales T (x) is determined by a 'flooding' model, according to which the probability density for locating the particle occupies a region which is constrained by a potential barrier which can trap a particle for time T .As T increases, higher and higher barriers are required.For a Gaussian distribution of barrier heights, equation (56) implies that the dispersion is described as sub-diffusive, of the form x ∼ x exp D 2 (ln T ) 2 8C(0) (66) which is distinctively different from the power-law anomalous diffusion which has been reported by some authors [6,7,8].Our numerical investigations of the dynamics of equation ( 1) for different values of D, illustrated in figure 5, show a data collapse which is in excellent agreement with equation (63), verifying (66).

Fig. 4
Fig. 4 Plots of S(x) (a) and T (x) (b), for one realisation of the smooth random potential, with D = 1/8.The numerically evaluated integrals (solid curves) are in good agreement with approximations based on sums over maxima and minima (dashed lines).

Fig. 5
Fig.5Numerical results on the median T (x) of the mean first-passage time to reach x, for different values of the diffusion coefficient D. We plot logarithmic variables Y med = D 2 ln T (x)/2 (a) and Yav = D 2 ln T (x) (b) as a function of X = 2D 2 ln(x/x) (where x is the mean distance between maxima of V (x)).The numerical results converge towards the theoretically predicted scaling function, equation (63), as D → 0.