Diffusion in periodic, correlated random forcing landscapes

Transport in random media is extensively studied due to its practical and fundamental importance [1–3]. In many cases, the dynamics is modelled as a Langevin process, with a drift generated by a quenched disordered potential. In theoretical analysis, the potential landscape is taken to be either infinitely extended or periodic in space. Stochastic dynamics in a periodic potential, both random and deterministic, is commonly encountered in many different contexts, including modulated structures [4], superionic conductors [5], colloids in light fields [6, 7], diffusion on regular [8–11] and disordered [12–16] solid surfaces, molecular motors on disordered tracks like DNA/RNA [17–19], and motion in a tilted potential due to a random polymer [20].

Theoretical approaches often assume that the dynamics in a periodic potential reproduces the behavior in an infinitely extended potential. This is implemented by setting the period in the final result, e.g. for the velocity (if any) or the diffusion coefficient, to infinity [21, 22]. It is crucial to investigate how far such an assumption holds. Especially in the context of numerical simulations carried out for periodic systems, one may ask how reliably their results may be extrapolated to infinite systems.

In this work, we address these fundamental questions for a Langevin dynamics x(t) in a periodic, quenched random potential $V(x(t))$ (with $V(x(t)+L)=V(x(t))$):

$$\begin{eqnarray}&&\eta \;\frac{{\rm d}x}{{\rm d}t}=-\frac{{\rm d}V(x)}{{\rm d}x}+\xi (t),\end{eqnarray} \tag{ 1 }$$

with η the friction coefficient, $\xi (t)$ a Gaussian white noise with zero mean and correlations $\overline{\xi (t)\xi (t^{\prime} )}=2\eta T\delta (t-t^{\prime} )$, the overbar being an average over the noise, and the temperature T is in units of the Boltzmann constant. We consider two cases:

• the ratchet case where V(x) is a fractional Brownian motion (fBm) in time $x\in [0,L]$, with $V(0)=0$ and V(L) arbitrary. Thus, V(x) is a Gaussian process with zero mean, $\langle V(x)\rangle =0$, and variance
  $$\begin{eqnarray}&&\langle {{[V(x)-V(y)]}^{2}}\rangle =\frac{V_{0}^{2}}{{{l}^{2\,H}}}|x-y{{|}^{2\,H}};\;\;x,y\in [0,L],\end{eqnarray} \tag{ 2 }$$
  where $H\in (0,1)$ is the Hurst exponent, V₀ and l define respectively the typical amplitude of V(x) and its scale of variation over x. In (2) the angular brackets denote averaging with respect to V(x). Figure 1 (left) shows a realization of V(x), with a jump at x = L.
• the translationally invariant case, where V(x) is a stationary Gaussian process, which at short length scales $|x-y|\ll L$ has the variance (2), and satisfies $V(0)=V(L)$, so that all points are statistically equivalent. The particle in this case diffuses on a ring. A realization of such a V(x) is shown in figure 1 (right).

Zoom In Zoom Out Reset image size

The dynamics (1) involves a combination of two paradigmatic situations: random motion in a periodic potential and the generalized Sinai dynamics in presence of a force $F(x)=-{\rm d}V(x)/{\rm d}x$ that is a time-independent stochastic variable with spatial correlations (except for $H=1/2$ when V (x) is the trajectory of a Brownian motion itself so that (1) is the periodic Sinai model [23]). While the latter produces an archetypal subdiffusion with logarithmically-confined trajectories, the periodicity of the random potential enforces a long–time diffusive behavior with a diffusion coefficient D_L . Here, we show that the trade–off between these two competing trends results in a rich statistical behavior of D_L . In particular, D_L is strongly non-self-averaging, with both negative and positive moments exhibiting an anomalous dependence on the temperature, period L and the order of the moment, and being supported by atypical realizations of V (x). For the ratchet case, we obtain exact analytical results, relying on exact bounds, for both positive and negative moments of D_L . We also discuss the full form of the probability distribution of D_L , and show that it is characterized by a log–normal left tail and a highly singular log–stable right tail, reminiscent of a Lifshitz singularity. We finally highlight the issue of sample–to–sample fluctuations of D_L . From standard scaling arguments and physical intuition, one expects that our exact results for the ratchet case also hold for the translationally invariant situation, which is harder to analyze analytically. This is confirmed below by thorough numerical simulations [24].

The dynamics (1) in an infinite system for arbitrary H, where $H\gt 1/2$ ($H\lt 1/2$) implies positively (negatively) correlated increments and superdiffusive (subdiffusive) V (x), respectively, was discussed in [25] where it was shown that ${{{\rm lim} }_{t\to \infty }}\langle \overline{{{x}^{2}}(t)}\rangle \sim {{{\rm ln} }^{2/H}}(t)$ (see also [26]). In contrast, in a periodic system, the long–time motion is diffusive for any given realization of the potential V (x), so that we have the diffusion coefficient ${{D}_{L}}\equiv {{{\rm lim} }_{t\to \infty }}\overline{{{x}^{2}}(t)}/(2t)$, with D_L given by [27–31] (see also [7, 8, 12, 32])

$$\begin{eqnarray}&&\frac{{{D}_{L}}}{{{D}^{(0)}}}={{(\int _{0}^{L}\frac{{\rm d}x}{L}\int _{0}^{L}\frac{{\rm d}y}{L}{{{\rm e}}^{\beta [V(x)-V(y)]}})}^{-1}},\end{eqnarray} \tag{ 3 }$$

where β is the inverse temperature, and ${{D}^{(0)}}=T/\eta$. Clearly, D_L is a random variable that fluctuates between realizations of V (x), and has support on $[0,{{D}^{(0)}}]$. The inverse of D_L may be regarded as a product of partition functions in potentials V (x) and $-V(x)$, respectively. The Brownian version of this quantity finds applications in disordered systems and has been extensively studied, while our results for $H\ne 1/2$ apply to more general situations (note that the marginal case H = 0, when V (X) is log-correlated, was studied in [33]). The expression in equation (3) also describes the ground state energy in a toy model of localization, and its average value was studied in [34] for $H=1/2$.

Turning to the discussion of the behavior of D_L , we first reduce the number of parameters. In the following, we measure L in units of l [see equation (2)], absorb V₀ into β, and measure D_L in units of ${{D}^{(0)}}$, so that D_L has support on $[0,1]$ ⁵ . Now, the typical behavior of D_L is easy to estimate as $D_{L}^{{\rm typ}}\propto {{L}^{2}}/{{\tau }_{{\rm typ}}}$, where the dimensionless L sets the scale of an inhomogeneous region, and ${{\tau }_{{\rm typ}}}$ defines the typical (dimensionless) time a particle spends in this region. The random potential V (x) being a fBm, the typical height of the potential barrier over a length L scales like L^H , for both the ratchet and the translationally invariant case. Assuming Arrhenius–type activation, one expects ${{\tau }_{{\rm typ}}}\sim {{{\rm e}}^{\beta {{L}^{H}}}}$, which implies

$$\begin{eqnarray}&&D_{L}^{{\rm typ}}\sim {{L}^{2}}{\rm exp} (-\beta {{L}^{H}}).\end{eqnarray} \tag{ 4 }$$

The average behavior of D_L is a much more delicate question because, as can be seen from equation (3), computing the statistics of D_L is a highly non-trivial task that involves the study of an exponential functional of fBm, for which standard methods like the Feynman–Kac formula are of little use for $H\ne 1/2$.

Let us first summarize our main analytical results obtained for the ratchet case: We find that the average of the logarithm of D_L is given, to leading order in L, by

$$\begin{eqnarray}&&\langle {\rm ln} {{D}_{L}}\rangle \propto -2\;m\beta {{L}^{H}},\end{eqnarray} \tag{ 5 }$$

with $m=\langle {{{\rm max} }_{s\in [0,1]}}V(s)\rangle$. The result in equation (5) is consistent with the logarithmic growth of the disorder–averaged mean–square displacement in an infinite system: ${{{\rm lim} }_{t\to \infty }}\langle \overline{{{x}^{2}}(t)}\rangle \sim {{{\rm ln} }^{2/H}}(t)$, and the argument leading to equation (4). Further on, we obtain sharp bounds for the positive moments $(k\gt 0)$ of the random variable D_L :

$$\begin{eqnarray}&&{{A}_{k}}(L)\leqslant \langle D_{L}^{k}\rangle \leqslant {{B}_{k}}(L),\end{eqnarray} \tag{ 6 }$$

where, in the limit $L\to \infty$, the bounds satisfy

$$\begin{eqnarray}\begin{array}{rcl} {{A}_{k}}(L) & = & {\rm exp} [-(1+H){{(2k\beta )}^{\frac{1}{1+H}}}{{(\frac{C}{H}L)}^{\frac{H}{1+H}}}], \\ {{B}_{k}}(L) & = & {\rm exp} [-(1+H){{(k\beta )}^{\frac{1}{1+H}}}{{(\frac{c}{H}L)}^{\frac{H}{1+H}}}], \\ \end{array}\end{eqnarray} \tag{ 7 }$$

with constants c and C, $0\lt c\leqslant C\lt \infty$, being independent of L, k and β. Both bounds exhibit the same dependence on $L,k,\beta$, from which we infer that the exact asymptotic result has the same functional form. Finally, for the negative moments of D_L , we find

$$\begin{eqnarray}&&\langle D_{L}^{-k}\rangle \sim {\rm exp} (\frac{{{a}^{2}}{{k}^{2}}{{\beta }^{2}}{{L}^{2H}}}{4}),\end{eqnarray} \tag{ 8 }$$

where a is a constant (independent of $L,k$ and β). Our exact results, equations (6) and (8), thus show that the positive and negative moments are dominated by atypical realizations of V (x), in contrast to $\langle {\rm ln} {{D}_{L}}\rangle$, see equation (5).

We now turn to a derivation of our results. Using equation (3), the logarithm of D_L can be formally written as ${\rm ln} {{D}_{L}}={\rm ln} {{J}_{+}}(L)+{\rm ln} {{J}_{-}}(L)+2{\rm ln} L$, where ${{J}_{\pm }}(L)$ are stationary currents through a finite sample of length L with potentials $\pm V(x)$ [26]: ${{J}_{\pm }}(L)={{[\int _{0}^{L}{\rm d}x\;{\rm exp} [\pm \beta V(x)]]}^{-1}}$. Statistical properties of these currents for the Sinai problem ($H=1/2$) are known [35–40]. Using the results of [41], and noting that with $\langle V(x)\rangle =0$, ${{J}_{+}}(L)$ and ${{J}_{-}}(L)$ have equal moments, we have for arbitrary H and to leading order in L, $\langle {\rm ln} {{J}_{+}}(L)\rangle =\langle {\rm ln} {{J}_{-}}(L)\rangle \propto -m\beta {{L}^{H}}$, which yields equation (5).

The proof of the result in equation (6) is based on a Theorem due to Monrad and Rootzén [42] on the probability that a fBm V (x), with $V(0)=0$, remains within a strip of width for the time $x\in [0,L]$. Defining ${{M}_{L}}\equiv {{{\rm max} }_{0\geqslant x\geqslant L}}|V(x)|$, the Monrad–Rootzén theorem, in our notation, states that $P({{M}_{L}}\leqslant \epsilon )$ satisfies

$$\begin{eqnarray}&&{\rm exp} (-CL\;{{\epsilon }^{-\frac{1}{H}}})\leqslant P({{M}_{L}}\leqslant \epsilon )\leqslant {\rm exp} (-cL\;{{\epsilon }^{-\frac{1}{H}}}),\end{eqnarray} \tag{ 9 }$$

for $0\lt \epsilon \leqslant {{L}^{H}}$, c and C being L-independent constants [see equation (6)].

Consider the lower bound in equation (6). Suppose we average the positive definite quantity D^k _L by considering instead of the entire set Ω of all possible paths V (x) only a subset $\Omega ^{\prime} \subset \Omega$ of paths such that ${{M}_{L}}\leqslant \epsilon$. This gives the lower bound ${{\langle D_{L}^{k}\rangle }_{\Omega }}\geqslant {{\langle D_{L}^{k}\rangle }_{\Omega ^{\prime} }}P({{M}_{L}}\leqslant \epsilon )$. However, for paths in $\Omega ^{\prime}$, we have ${\rm exp} (\beta [V(x)-V(y)])\leqslant {\rm exp} (2\beta \epsilon )$, and hence, ${{(\int _{0}^{L}\frac{{\rm d}x}{L}\int _{0}^{L}\frac{{\rm d}y}{L}{{{\rm e}}^{\beta [V(x)-V(y)]}})}^{-1}}\geqslant {{{\rm e}}^{-2\beta \epsilon }}$. Therefore, we obtain ${{\langle D_{L}^{k}\rangle }_{\Omega }}\geqslant {{{\rm e}}^{-2k\beta \epsilon }}\;P({{M}_{L}}\leqslant \epsilon )$. Making the inequality more stringent by choosing the lower bound in equation (9), we get

$$\begin{eqnarray}&&{{\langle D_{L}^{k}\rangle }_{\Omega }}\geqslant {{{\rm e}}^{-2k\beta \epsilon }}\;{\rm exp} (-CL{{\epsilon }^{-\frac{1}{H}}}),\end{eqnarray} \tag{ 10 }$$

which holds for any with $0\lt \epsilon \leqslant {{L}^{H}}$.

The function on the right-hand side (rhs) of equation (10) is a non-monotonic function of , attaining its maximum at $\epsilon ={{\epsilon }_{{\rm opt}}}={{(CL/k2\beta H)}^{H/(1+H)}}$. Clearly, the best lower bound corresponds to the choice $\epsilon ={{\epsilon }_{{\rm opt}}}$, leading to the lower bound in equation (6). Note that to satisfy the conditions of validity of the Monrad–Rootzén theorem, we require that ${{\epsilon }_{{\rm opt}}}\leqslant {{L}^{H}}$, that is, $2k\beta L\geqslant C/H$, which is easily realized for sufficiently large L. The derivation of the lower bound is an example of the Lifshitz optimal fluctuation method [43], which has been used to bound the survival probability of particles diffusing in the presence of randomly scattered immobile traps (see, e.g., [44]).

We now discuss the derivation of the upper bound. To this end, we discretize x, and write the rhs of equation (3) as

$$\begin{eqnarray}&&\int _{0}^{L}\frac{{\rm d}x}{L}\int _{0}^{L}\frac{{\rm d}y}{L}{{{\rm e}}^{\beta [V(x)-V(y)]}}\sim \mathop{\sum }\limits_{i,j=1}^{N}{{{\rm e}}^{\beta [V(j)-V(i)]}}.\end{eqnarray} \tag{ 11 }$$

Given that the fBm starts at $V(0)=0$, at least one term in the double sum on the rhs takes the value ${\rm exp} (\beta {{M}_{L}})$, corresponding to the point x = 0 and the point where $|V(x)|$ attains its maximal value M_L . Now, as all the other terms are positive, we have the bound $D_{L}^{-1}\geqslant {\rm exp} (\beta {{M}_{L}})$, and, thus,

$$\begin{eqnarray}&&\langle D_{L}^{k}\rangle \leqslant \langle {{{\rm e}}^{-k\beta {{M}_{L}}}}\rangle =k\beta \int _{0}^{\infty }{\rm d}\epsilon P({{M}_{L}}\leqslant \epsilon ){{{\rm e}}^{-k\beta \epsilon }}.\end{eqnarray} \tag{ 12 }$$

The integral in the rhs is dominated, for large L, by the small region, where we can thus use the upper bound in (9). Performing the remaining integral over by the saddle-point method and omitting the pre-exponential terms lead to the upper bound in equation (6). The result in (6) has several striking features. Namely, the function

$$\begin{eqnarray}&&\mu (k,\beta ,L)=-{\rm ln} \langle D_{L}^{k}\rangle ,\end{eqnarray} \tag{ 13 }$$

(as compared to its typical counterpart given by equation (4) as ${{\mu }_{{\rm typ}}}(k,\beta ,L)\sim k\beta {{L}^{H}}$), (a) grows sub-linearly with k (multifractality), (b) is a non-analytic, sublinear function of β, which implies a rather unusual sub-Arrhenius dependence of the positive moments on the temperature, and (c) exhibits a slower anomalous growth with L as $\sim {{L}^{\frac{H}{1+H}}}$. This means that the disorder-averaged D_L is generically larger than the one expected on the basis of typical realizations of the disorder. In turn, this implies that the behavior of the positive moments of D_L is supported by atypical realizations of disorder, reminiscent of the so-called Lifshitz singularities, as discussed above. In conclusion, one cannot infer the dynamical behavior in an infinite system from the positive moments of D_L . This is surprising at first glance, as $\langle \overline{{{x}^{2}}(t)}\rangle$ is linearly proportional to $\langle {{D}_{L}}\rangle$, and shows that the limits $t\to \infty$ and $L\to \infty$ do not commute in this system.

The behavior of the negative moments $\langle D_{L}^{-k}\rangle$ with $k=1,2,...$ is determined by essentially the same approach as above. Note that both the lower and the upper bound on $D_{L}^{-k}$ are made tighter for a given realization of V (x) by using $D_{L}^{-k}\sim {\rm exp} (k\beta S)$, where S is the span of V (x) (the difference between the maximum and minimum) on the interval $[0,L]$. Therefore, in contrast to the positive moments, the negative moments are supported by realizations of V (x) with a large span. Using the result that for large S, $P({{M}_{L}}=S)\sim {\rm exp} (-{{S}^{2}}/a{{L}^{2\,H}})$, with a a constant, integration of equation (12) gives the result announced in equation (8), which displays a super-Arrhenius dependence on the temperature, a superlinear dependence on k, and a strong dependence on L. A similar result was obtained earlier in [45]. We note that, as for the positive moments, one cannot deduce the behavior in an infinite system from that of negative moments of D_L in a periodic system, as the latter is supported by atypical realizations of V (x) that have anomalously large span scaling as $S\sim {{L}^{2\,H}}$, while the typical behavior is ${{S}_{{\rm typ}}}\sim {{L}^{H}}$.

Based on our results for the moments, we now obtain the probability distribution $P({{D}_{L}})$. As already explained, the behavior of the negative moments is supported by anomalously stretched trajectories of V (x) for which the value of D_L is small. One may thus expect in view of the form of the moments in equation (8) that for small D_L , $P({{D}_{L}})$ is log-normal:

$$\begin{eqnarray}&&P({{D}_{L}})\sim \frac{1}{a\beta {{L}^{H}}{{D}_{L}}}{\rm exp} (-\frac{{{{\rm ln} }^{2}}({{D}_{L}})}{{{a}^{2}}{{\beta }^{2}}{{L}^{2\,H}}}).\end{eqnarray} \tag{ 14 }$$

To analyze the behavior of $P({{D}_{L}})$ for D_L close to 1, we recall the formal definition of the one-sided Lévy distribution ${{}_{\nu }}(z)$, $0\leqslant z\lt \infty$, of order ν (see, e.g., [46]):

$$\begin{eqnarray}&&\int _{0}^{\infty }{\rm d}z{{{\rm e}}^{-pz}}{{}_{\nu }}(z)={{{\rm e}}^{-{{p}^{\nu }}}}.\end{eqnarray} \tag{ 15 }$$

The asymptotic behavior of ${{}_{\nu }}(z)$ is well-known [46], and, in particular, one has ${{}_{\nu }}(z)\sim {{z}^{-\sigma }}{\rm exp} (-b/{{z}^{\tau }})$ for $z\to 0$, where b is a computable constant, $\sigma =(2-\nu )/(2(1-\nu ))$, and $\tau =\nu /(1-\nu )$. It is important to note that this precise asymptotic form is responsible for the stretched-exponential behavior in equation (15), which is immediately verified by substituting the form in equation (15), and performing the integration by the saddle-point method. Moreover, one realizes by making in equation (15) a change of the integration variable $z={\rm ln} (1/{{D}_{L}})/\beta {{L}^{H}}$, choosing $\nu =1/(1+H)$, and setting $p=k/\beta {{L}^{H}}$ that equation (15) becomes identical to the result in equation (6), up to numerical factors. It follows that for D_L close to 1 (i.e., z close to 0), the distribution function behaves as

$$\begin{eqnarray}&&P({{D}_{L}})\sim \frac{1}{\beta {{L}^{H}}{{D}_{L}}}{{}_{\frac{1}{1+H}}}[{\rm ln} (D_{L}^{-1})/\beta {{L}^{H}}].\end{eqnarray} \tag{ 16 }$$

Using the asymptotic ${{}_{\nu }}$ given above, we get that $P({{D}_{L}})$ is highly singular near the right edge, $P({{D}_{L}})\sim {\rm exp} [-b{{\beta }^{1/H}}L/{{(1-{{D}_{L}})}^{1/H}}]$, similar to the Lifshitz singularity.

We now consider the translationally invariant case. Here, we expect our above analysis, in particular, result (6) to hold, up to possible numerical factors. To demonstrate this, we now present results of extensive numerical simulations: figure 2 for $-{\rm ln} (\langle {{D}_{L}}\rangle /{{L}^{2}})$, $-\langle {\rm ln} ({{D}_{L}}/{{L}^{2}})\rangle$, and $\mu (k,\beta ,L)$, and figure 3(a) for $P({\rm ln} {{D}_{L}})$ indeed show a very good agreement that supports our expectations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To close, we ask: if we have two different realizations of V (x), and correspondingly, two different values, D_L and $D{{^{\prime} }_{L}}$, of the diffusion coefficient, how likely are these values equal? We introduce a random variable $\mathcal{D}\equiv \frac{{{D}_{L}}}{{{D}_{L}}+D{{^{\prime} }_{L}}}\;,\;\;\mathcal{D}\in [0,1]$, and analyze its distribution $P(\mathcal{D})$ via numerical simulations. Clearly, $\mathcal{D}=1/2$ maximizing $P(\mathcal{D})$ implies that the two values of D_L are most likely very close to one another. Variables such as $\mathcal{D}$ play a key role in various scale-independent hypothesis testing procedures, in classical problems in statistics, in signal processing (see, e.g., [47]), and in the analysis of chaotic scattering in few-channel systems [48]. Such variables are used to characterize the effective width of narrow distributions possessing moments of arbitrary order [49, 50].

In figure 3(b), we present numerical results for $P(\mathcal{D})$ for different values of L and $H=1/3$, for the translationally invariant case. We observe an interesting phenomenon of a change in the form of the distribution as L is increased. For relatively small L, the distribution is bell-shaped and centered at $\mathcal{D}=1/2$. However, on increasing L, $P(\mathcal{D})$ broadens, becomes almost flat at a certain critical L (whose value depends on H), and then changes its shape so that $\mathcal{D}=1/2$ minimizes the distribution. This implies that for sufficiently large L, two values of the diffusion coefficients obtained for two different realizations of V (x) are most likely very different, and the event ${{D}_{L}}=D{{^{\prime} }_{L}}$ is the least probable. Note that a similar dependence in the distribution of the prefactor in the Sinai law on the strength of disorder was recently observed in [51].

Acknowledgments