Persistent Sinai type diffusion in Gaussian random potentials with decaying spatial correlations
Igor Goychuk Vasyl O. Kharchenko Ralf Metzler





igoychuk@uni-potsdam.de, corresponding author
Institute of Physics and Astronomy, University of Potsdam, 
Karl-Liebknecht-Str. 24/25, 14476 Potsdam-Golm, Germany

Institute of Applied Physics, Natl. Acad. Sci. Ukraine, 
58 Petropavlovskaya str., 40030 Sumy, Ukraine

Institute of Physics and Astronomy, University of Potsdam, 
Karl-Liebknecht-Str. 24/25, 14476 Potsdam-Golm, Germany




Logarithmic or Sinai type subdiffusion is usually associated with random force
disorder and non-stationary potential fluctuations whose root mean squared
amplitude grows with distance. We show here that extremely persistent,
macroscopic ultraslow logarithmic diffusion also universally emerges at
sufficiently low temperatures in stationary Gaussian random potentials
with spatially decaying correlations, known to exist in a broad range of
physical systems. Combining results from extensive simulations with a scaling
approach we elucidate the physical mechanism of this unusual subdiffusion. In
particular, we explain why with growing temperature and/or time a first
crossover occurs to standard, power-law subdiffusion, with a time-dependent
power law exponent, and then a second crossover occurs to normal diffusion with
a disorder-renormalized diffusion coefficient. Interestingly, the initial,
nominally ultraslow diffusion turns out to be much faster than the universal
de Gennes-Bässler-Zwanzig limit of the renormalized normal diffusion, which
physically cannot be attained at sufficiently low temperatures and/or for
strong disorder. The ultraslow diffusion is also non-ergodic and displays
a local bias phenomenon. Our simple scaling theory not only explains our
numerical findings, but qualitatively has also a predictive character.


05.40.-a, 82.20.Wt, 87.10.Mn, 87.15.Vv, 87.15.hj

Persistent Sinai type diffusion in Gaussian random potentials with decaying 
spatial correlations
    Ralf Metzler
    December 30, 2023
=================================================================================================






§ INTRODUCTION


Systems with Gaussian energy disorder characterized by spatially decaying
correlations are ubiquitous in physics thanks to the central limit theorem
<cit.>. Such systems include, for instance, disordered organic
photoconductors with long-range electrostatic interactions <cit.>, supercooled liquids <cit.>, as well as naturally
occurring DNA macromolecules encoding biological information in living systems
<cit.>. Colloidal systems in quenched, random
laser-created potentials have also recently become experimentally available
<cit.>.

A common line of thinking  <cit.> treats diffusion and transport
phenomena in such systems as normal on experimentally relevant scales, with a
potential disorder-renormalized diffusion coefficient D_ ren=D_0exp(-[
βσ]^2), where D_0 is the free diffusion coefficient of the diffusing
particle in absence of the potential, σ=⟨δ U^2(x)⟩^1/2
is the root-mean-squared amplitude of the potential fluctuations δ U(x),
and β=1/(k_BT) is the Boltzmann factor proportional to inverse temperature.
This famed renormalization result by De Gennes <cit.>, Bässler
<cit.>, and Zwanzig <cit.> corresponds to a common
measurable temperature dependence of transport coefficients in disordered glassy
systems which is almost indistinguishable from the Vogel-Fulcher-Tammann law
<cit.>, another commonly used temperature dependence used to fit
experimental data.

It came as a surprise when extensive simulations of stochastic Langevin dynamics
in Gaussian random potentials with decaying spatial correlations by Romero and Sancho
<cit.> demonstrated that the emerging diffusion is anomalously slow over the 
whole time range of simulations and
characterized by a power law scaling ⟨δ x^2(t)⟩∝ t^α (0<α
<1) of the mean squared displacement.
This result is very remarkable indeed because within the mean field approximation it
is easy to show that the Gaussian energy disorder yields residence time distributions
ψ(τ) of a log-normal type in finite spatial trapping domains, all moments
of which are finite. This behavior is fundamentally different from the situation
of annealed exponential energetic disorder yielding residence time distributions
of the form ψ(τ)∝1/τ^1+α with α=k_BT/σ
<cit.>: Here subdiffusion emerges for k_BT<σ because of the divergent
mean residence times, which provides the standard model of subdiffusion in the
continuous time random walk framework <cit.>.

Hence, Gaussian energy disorder cannot yield subdiffusion in the absence of spatial
correlations. To neglect such correlations is a usual procedure when they are short
ranged <cit.>. Whether organic photoconductors and other disordered
materials are better described by Gaussian or exponential models of the energy
disorder is a subject of continuing controversy. On one hand, the Gaussian model is
generally accepted for organic photoconductors <cit.>. On the other
hand, some recent experiments <cit.> seem to be more consistent with
the model of exponential disorder, leading to a similar time dependence of transient
currents as in amorphous semiconductor films <cit.>. Models of Gaussian
disorder with significant spatial correlations can be a key to resolve this
controversy. Physically, for instance, such short range correlations may correspond
to the size of a protein molecule diffusing over the potential landscape of a DNA
strand <cit.>, or to the size of a colloidal particle in an
uncorrelated laser field potential <cit.>. 
Morever, long-range electrostatic interactions
give rise to long-range spatial correlations in organic photoconductors <cit.>,
and biological information encoded in the sequence of DNA base pairs manifests itself in 
long-range base pair correlations <cit.>.

The remarkable discovery that subdiffusion persists beyond the mean field
approximation was reaffirmed in a number of recent studies <cit.>. It was explained recently in terms of a
finite non-ergodicity length L_ erg of the spatial random process
exp[βδ U(x)] <cit.> related to the presence of
spatial potential correlations. In this respect, the origin of a weak
ergodicity breaking—the fundamental difference between ensemble and time
averaged physical observables even in the limit of long observations times
<cit.>—is very different from the one in
systems described by continuous time random walks with annealed, uncorrelated
residence time distributions with divergent mean residence times <cit.>. For systems with Gaussian
disorder, the non-ergodic behavior occurs on the transient spatial scale L_
erg which, however, rapidly grows with βσ for spatially correlated
potential fluctuations δ U(x). Indeed, this spatial scale can readily reach
macroscopic sizes already for disorder strengths σ∼4÷5k_BT, which is
a typical value for disordered organic photoconductors already at room temperatures
<cit.>.

Subdiffusion occurs on the corresponding non-ergodic and generally mesoscopic
spatial scales. Somewhat surprisingly, diffusion retains its asymptotically normal
form for any decaying correlations, no matter how slow they decay
<cit.>.  However, for mesoscopic and even macroscopic systems this
limit can be out of reach and become physically completely irrelevant, leading
us to consider this type of subdiffusion as one of the fundamental processes in
a large range of applications. Some of these remarkable features have recently
been confirmed experimentally in colloidal systems <cit.> and, more
tentatively, also for diffusion of regulatory proteins on DNA strands
<cit.>. Namely, although in Ref. <cit.> the results are
interpreted in terms of normal diffusion, the single-trajectory diffusion
coefficient is scattered over three orders of magnitude, which likely implies
a lack of ergodicity. The real physical mechanism of anomalous diffusion in
spatially correlated Gaussian energetic disorder remains elusive, and it is
the main goal of this paper to clarify its physical origin. The whole setup is
very different from the original scenario of Sinai diffusion <cit.>.
It is also different from Sinai diffusion emerging from a
power law energy disorder <cit.>, or logarithmic tails of the waiting time distribution <cit.> in the mean field approximation.  Notwithstanding,
we will show below that a generalized Sinai subdiffusion indeed emerges at
sufficiently low temperatures with ⟨δ x^2(t)⟩∝[ln(t/
t_0)]^a, with a general scaling exponent a, which is a=4 in the special
case of Sinai diffusion <cit.>. 
The latter one is usually
associated with the model of random force disorder <cit.> yielding potential fluctuations whose root mean square amplitude
exhibits an unlimited growth with distance, that is, the potential presents a
free Brownian motion in space <cit.>. A 
generalized Sinai diffusion with a≠ 4 emerges also in the case of
fractional Brownian motion in space <cit.>.

In this work, we reveal that a generalized Sinai type diffusion  universally
emerges for a sufficiently low temperature and/or strong disorder when βσ exceeds a typical value of 5 to 10, even for the case of very short range
exponential or linearly decaying correlations. This fundamental result is fully
confirmed by extensive simulations. Based on the arguments developed herein it
it should be absolutely possible to observe macroscopic Sinai diffusion at
experimentally relevant time scales, in systems governed by spatially correlated
Gaussian energetic disorder.

The remainder of this manuscript is structured as follows. In Section II we
develop the general concepts of our approach and introduce a scaling theory for
the resulting diffusive dynamics. Section III presents results from extensive
simulations for various relevant cases of the spatial correlations decay of the
potential fluctuations. In Section IV we present a detailed discussion of our
results and draw our conclusions.




§ BASIC CONCEPTS AND SCALING THEORY





 §.§ Growing potential fluctuations


We start from an observation, which is central to the theory of extreme
events <cit.> and pivotal for the present study. Namely, even if
a stationary Gaussian process U(x) has a finite root mean squared fluctuation 
σ,
the maximal amplitude of its fluctuations δ U_ max (x) grows slowly
with the distance x, at odds with intuition <cit.>. The
law of this growth can be found from the following consideration. We consider
the potential as a random process with zero average, ⟨ U(x)⟩=0,
and ask the question of how often an energy level U_ max is crossed with
growing distance x→∞, in a stationary limit. For any stationary,
differentiable Gaussian process with normalized autocorrelation function g(x)=
⟨ U(x+x')U(x')⟩/σ^2 and asymptotically decaying correlations,
the answer is known: The averaged number of level-crossing events ⟨ n(x)
⟩ grows linearly with distance: ⟨ n(x)⟩=Γ x, at the
rate Γ=√(-g”(0))exp(-U_ max/2σ^2)/π <cit.>.
Thus, the level U_ max is typically crossed at the distance x found
from ⟨ n(x)⟩=1. Hence,

    U_ max(x)=σ [2ln (x/x_0) ]^1/2,

where x_0=π/√(-g”(0)). This is a very general and universal result
valid for any correlations with finite g”(0), in agreement with the
results of Ref. <cit.>. Note that this result holds only
asymptotically, given a stationary regime is established. Moreover, the longer
the correlations range, the later this asymptotic result is established.

For example, for a Gaussian decay of the correlations we will have g(x)=exp
(-x^2/λ^2) with a correlation length scale parameter 
λ, and in this case x_0=πλ/√(2).
Note that in a scaling sense such a correlation length parameter λ
always exists, even if the correlation length formally defined as λ_
corr=∫_0^∞|g(x)|dx may diverge. The length λ plays a
fundamental role in the presence of correlations. It can be, for instance,
of the size of the protein-DNA contact, or the size of the colloidal particle
in an uncorrelated laser field.

For power law decaying correlations with g(x)=1/[1+x^2/λ^2]^γ/2,
and g(x)∝ 1/x^γ for x≫λ, and thus x_0=πλ/√(γ). This case is very interesting in applications to diffusion of regulatory
proteins on DNA strands, 
where the real base pair sequence arrangement can play a very profound role <cit.>.
In particular, longe-range correlations in the base pair sequence,
which encode biological information, can indeed decay algebraically slow, as shown,
for instance, in Refs. <cit.> with γ∼ 0.6÷0.8.
This should yield corresponding correlations in the protein binding energy with
λ being of the size of the protein-DNA contact, or larger. Namely, we here
have a case for which λ_ corr=∞, which suggests that the
subdiffusive regime should hold much longer than in the case of short ranged
correlations. This is indeed confirmed in our detailed analysis below. Another
important example is provided by some organic photoconductors, for which γ=1
and σ=0.1÷0.125 eV or 4÷5 k_BT at room temperature
<cit.>.

The next case of relevance we consider here is the Ornstein-Uhlenbeck
process with g(x)=exp(-|x|/λ). This case may at first appear problematic
as it is not differentiable, and the corresponding force f(x)=-∂ V(x)/
∂ x has the infinite root mean squared amplitude ⟨ f^2(x)⟩^
1/2=∞. However, in physical applications we consider in fact a regularized
Ornstein-Uhlenbeck process on a lattice with a finite grid size Δ x. This
allows also for a numerical treatment of the corresponding continuous dynamics
(which otherwise would not be possible due to the infinite force root mean squared
amplitude). We note that this point is similar to the subtleties arising when we
consider the model of white Gaussian processes f(x) in the theory of continuous
space Sinai diffusion <cit.>. Indeed, we always generate realizations
of Gaussian processes with some finite resolution Δ x, following a well
established algorithm <cit.>. The corresponding g(x) is getting, in
fact, thereby regularized. An appropriate formal regularization reads g(x)=exp(
-√(x^2+(Δ x/2)^2)/λ+Δ x/(2λ)). It yields the same force
root mean squared ⟨ f^2(0)⟩^1/2=√(2)σ/√(Δ x
λ) as the Ornstein-Uhlenbeck process on a lattice with grid size Δ x.
In this case, x_0=π(λΔ x/2)^1/2. An important further remark is
that, when Δ x→0, then x_0→0, and U_ max(x) grows accordingly
for any fixed x. However, this growth is logarithmically slow, and for realistic
variations of Δ x the effect is small. Nevertheless, the case of
exponentially decaying correlations is especially interesting in the present
context. Notice that such short range correlations appear also in the case
of a biologically meaningless and completely random DNA because of a finite
size λ of the protein-DNA contact, which is typically from 5 to 30
bp <cit.>.

Another pertinent correlation model for the particular application of diffusion
on DNA is the model of linearly decaying correlations g(x)=(1-|x|/λ) for
|x|≤λ, and zero otherwise <cit.>. It emerges when the protein
interacts locally with only one corresponding bp, in the immediate contact. The
corresponding random process U(x) is also non-differentiable, and with respect
to the U_ max(x) behavior its regularized version behaves rather similar
to that of the Ornstein-Uhlenbeck process (considered with the same Δ x).

 

Let us check the above exact analytical results against the results from
extensive simulations: Fig. <ref> shows how the amplitude of the potential
fluctuations grows with distance for several models of Gaussian disorder. The
universal result in Eq. (<ref>) agrees well with the numerics for all four
models considered here. A similar scaling was also found numerically for a
different model of correlations <cit.>. The agreement is better for
differentiable (in the limit Δ x→ 0) processes, and the model of linearly
decaying correlations behaves indeed rather similar to the Ornstein-Uhlenbeck model.
Notice that in the case of power law decaying correlations with γ=0.8, the
convergency is slower than for the Gaussian decay, and for much smaller γ,
for instance, γ=0.2, it is is still far from being achieved for the largest
x in Fig. <ref> (not shown). The results in Fig. <ref> are obtained
from a well-known spectral method <cit.> given the corresponding
correlation functions, on a lattice with grid size Δ x (Δ x=0.02 in
Fig. <ref>) and a maximal spatial interval L_ max (L_ max=NΔ
x, N=2^20 in Fig. <ref>). Notice also that the exact behavior of the
function g(x) is used to generate the realizations of the Ornstein-Uhlenbeck
process and the process with linearly decaying correlations. The shown results
are obtained from averaging over 10^3 uniformly distributed points (particles)
in each potential realization, and 10^2 such potential realizations for
each model of correlations were taken.

It is clear that the strongest potential fluctuations occur in the case of
exponential and linear correlations. They are further increasing for smaller
values of the resolution Δ x due to the singularity of these models.
Another salient feature is that another scaling can be observed transiently.
Namely, we initially find the behavior log x rather than the √(log
x) scaling for the exponential and linear correlations, and intermediately
for the case of power law decaying correlations. Moreover, we notice an initial
(log x)^δ scaling with δ>1 for the power law and Gaussian decay.
To account for these various behaviors we henceforth consider the following
model

    |δ U_ max(x)|=σ_ eff[ln (x/x_ in) ]^δ

for the potential fluctuations, where σ_ eff, δ, and x_
in are fitting parameters, which are generally different from the theoretical
values σ_ eff=2√(2)σ≈ 2.83 σ, δ=1/2,
and x_ in=x_0 in Fig. <ref>. We will now use this generic empirical
formula and pursue the scaling argumentation used originally in the case of
continuous space Sinai diffusion <cit.>.




 §.§ Scaling theory


Let us estimate the time t a particle needs to travel the distance x starting
at x_0, limited in an Arrhenius manner by the largest barrier met on its way,
t=t_0exp[β|δ U_ max(x)|], where t_0 is a proportionality factor of
physical unit of time. From this scaling ansatz, in combination with relation
(<ref>) we immediately obtain our central result for the mean squared
displacement of the diffusing particle,

    ⟨δ x^2(t) ⟩ = x_ in ^2 
    {e^[(k_BT/σ_ eff)ln(t/t_0)]^1/δ-1 }^2.

This in turn yields the asymptotic form

    ⟨δ x^2(t) ⟩≈ 
     x_ in^2 
     [(k_B T/σ_ eff)ln(t/t_0)  ]^2/δ

for small temperatures and/or strong disorder, k_BT/σ_ eff≪1. For
δ=1/2, this is precisely the Sinai type diffusion with a=2/δ=4.
However, for δ =1, ⟨δ x^2(t)⟩≈ x_ in^2[(k_B
T/σ_ eff)ln(t/t_0)]^2.

Next, with increasing temperature and time, when the unity becomes negligible in
Eq. (<ref>), we find the intermediate behavior

    ⟨δ x^2(t)⟩≈ x_ in^2 [t/t_0]^α(t)

with the time dependent scaling exponent

    α(t)=2(k_B T/σ_ eff)^1/δln^1/δ-1(t/t_0).

This result predicts a power-law subdiffusion with a gradually changing anomalous
diffusion exponent. It should be mentioned here that this result follows from the
approximation exp[f(t)]-1≈exp[f(t)], for a logarithmically growing f(t),
which is not fully accurate in our numerics. Nevertheless, it predicts the correct
law for α(t), namely, α(t)∝log t, for δ=1/2, see below.
This result is especially insightful for δ=1 (or log x scaling): here
α(t)= 2k_BT/σ_ eff=α_ min, and we observe subdiffusion
of power law form with a constant anomalous diffusion exponent. Moreover, given 
the results depicted in Fig. <ref>, where initially δ>1, intermediately
δ=1, and asymptotically δ=1/2, one can predict that initially
α(t) will diminish and reach a minimum, and then logarithmically increase
in time. This is what we actually see in the simulations with σ_ eff≈(1.24÷1.52)σ within the temperature range σ/4<k_BT<σ/2,
see below. Such a behavior may indeed mistakingly be attributed to continuous
time random walk subdiffusion.




 §.§ Transient lack of ergodicity


Due to the decaying correlations of the Gaussian potential fluctuations in our
model, with increasing temperature for a given root mean squared amplitude
σ (decreasing βσ) and for times longer than a certain
crossover time, a transition to normal diffusion will gradually emerge. 
The corresponding characteristic anomalous diffusion length L_ erg can be
found from the condition that the variance of the random process w(x)=exp(-β
U(x)) equals the mean, which yields the following implicit equation for this, so
far unknown, characteristic length L_ erg <cit.>,

    ∫_0^1 (1-y )e^β^2σ^2g(L y)dy=1.

The transition to normal diffusion starts at ⟨δ x(t)^2⟩⪆
L_ erg^2. However, this transition may last extremely long, and anomalous
diffusion features may persist for appreciably long times for |δ x|>L_
erg.

For the case of linearly decaying correlations, Eq. (<ref>) can be solved
exactly, and we obtain (L_ erg>λ)

    L_ erg   =   λ/(βσ)^2{(e^β^2σ^2-1)-(βσ)^2+(2(βσ)^4+4(βσ)^2
       -   2(βσ)^2e^(βσ)^2+3-4e^(βσ)^2+e^2(βσ)^2
    )^1/2}

which for βσ⪆2.5 is approximated by the simple expression

    L_ erg≈2λ/(βσ)^2e^(βσ)^2

with a very good accuracy. For other models of g(x), Eq. (<ref>) is solved
numerically, the results being listed in Table <ref> for several values of βσ.
One can see that the shortest length L_ erg occurs for the linearly decaying
correlations and the longest one for the power law correlations. Somewhat
surprisingly, for the Gaussian decay the non-ergodicity length is larger than for
exponentially decaying correlations.



As mentioned above, there exist several crossover times in the dynamics. The
first one corresponds to the transition from Sinai like diffusion to the power
law diffusion regime. The corresponding transition time τ^* can be roughly 
estimated from the condition that the argument of the exponential function in
Eq. (<ref>) reaches unity. From this, τ^*∼ t_0exp[σ_
eff/k_BT]. It grows exponentially fast with σ_ eff/k_BT. The second
transition time τ_ erg can be estimated from the condition of how long
the power law diffusion regime with anomalous diffusion exponent α_
min≈ 2k_BT/σ_ eff will last.


Thus we find the
conditions ⟨ x^2(t)⟩∼ x_ in^2(τ_ erg/t_0)^α_ min∼ L_ erg^2, 
from which with α_ min≈2k_BT/σ_ eff and Eq. (<ref>),
we obtain the estimate

    τ_ erg∼ t_0(4λ^2/x_ in^2(k_BT)^4/σ^4)^ σ_ eff/2k_BTexp[σ^2σ_ eff/
    (k_BT)^3].

 Notice the super-exponential growth of τ_ erg with σ/
(k_BT). This is precisely the reason why the power law subdiffusion regime can
last so long already for moderate values σ/(k_BT)∼4÷5, and no
transition to the regime of normal diffusion was revealed in the concrete cases
studied in Refs. <cit.>.




§ RESULTS OF NUMERICAL SIMULATIONS


We now proceed by checking our theoretical predictions against results of
extensive numerical simulations, finding remarkable agreement. To this end,
let us consider a continuous-space Brownian dynamics governed by the
overdamped Langevin equation  

    ηdx(t)/dt=-dU(x)/dx+√(2k_BTη)×ζ(t),
 
where η is the frictional coefficient and ζ(t) is unbiased, white
Gaussian noise with δ-correlation ⟨ζ(t)ζ(t')⟩=
δ(t-t'). Distance is scaled in units of λ, time in units of
τ_0=λ^2η/σ, and temperature in σ/k_B. Initially,
10^4  particles are uniformly distributed in random potentials (10
realizations), and the particle motion is integrated using periodic boundary
conditions with a very large period L_ max. The random potential is
generated on a lattice with spacing Δ x=0.02, and discrete points are
connected by parabolic splines, that is, the potential is locally parabolic,
and the corresponding force entering the Langevin equation (<ref>) is
piece-wise linear. In this respect our setup is similar to the one considered
in Refs. <cit.>, but different from that in
Refs. <cit.>, where a discrete hopping dynamics in both
space and time was studied. In most simulations, we employ L_ max=2^19Δ
x≈1.0485 ×10^4λ. We use the stochastic Heun method with a
time integration step Δ t=2×10^-4 in most simulations.




 §.§ Exponential correlations










We start our numerical analysis with the case of the exponentially decaying
correlations, g(x)=exp(-|x|). From the simulated trajectories, we first
perform an ensemble averaging to obtain the mean squared displacement. The
results are depicted in Fig. <ref> in panels (a) and (b). In Fig. <ref>a,
the fit of the numerical results is based on Eq. (<ref>) using δ=1/2,
σ_ eff∼4÷4.7, and x_ in∼0.59÷5.01. Note that the
fitting values σ_ eff differ only by a factor of less than two from
the theoretical value 2.83 in Fig. <ref>. The corresponding theoretical
value x_0 is x_0≈ 0.31. Here, the agreement is worse, however, the
numerical results in Fig. <ref> in this case are also somewhat different
from the theoretical asymptotics, which is still not reached. The same numerical
data are fitted to expression (<ref>) with δ=1 in Fig. <ref>b,
with σ_ eff∼1.23÷1.69. For T=0.1 and T=0.2, both fits have
almost the same quality, although the fit with δ=1 appears slightly better.
For T=0.1, Sinai like diffusion (<ref>) with a≈3.7, near to
the Sinai value a=4, covers about six decades in time, Fig. <ref>b.
It is worthwhile noting that for T≤ 0.25 the fit with δ=1 works better,
see Fig. <ref>b, due to the initial log x scaling in Fig. <ref> (see
symbols therein). For this reason, α(t) shows a nearly constant behavior,
at intermediate temperatures, for an extended time period, Fig. <ref>a. In
fact, in this regime the power law exponent is nearly proportional to temperature, 
α≈2k_BT/σ_ eff, where the value of σ_ eff is
obtained from fitting the minimal value α_ min in Fig. <ref>a by
this dependence. It turns out that σ_ eff≈1.52, Fig. <ref>.
However, with increasing temperature, when diffusion covers larger distances, the
fit with δ=1/2 becomes much better, compare the case T=0.5 in
Figs. <ref>a and <ref>b. An excellent fit holds over about six to seven
decades in time with δ=1/2, where the explicit time dependence of the
anomalous diffusion exponent α(t) in Eq. (<ref>) becomes apparent
for sufficiently long times. At larger temperatures, for instance, T=1, the
crossover to the asymptotic regime of normal diffusion is already accomplished
up to t_ max=10^6 in our simulations. Such high temperatures, or weak
disorder strengths, are not of interest for anomalous diffusion of the kind
discussed. Generally α_ eff reaches a minimum α_ min,
may stay nearly constant for a certain time period, and then logarithmically
increases, a dependence which can be fitted by Eq. (<ref>) with
δ=1/2 and some values of σ_ eff and t_0, which are different
from those used in Fig. <ref>a, see Fig. <ref>a for T=0.5 and T=1/3.
The reason for this discrepancy in the corresponding fitting values is that a
transition from Eq. (<ref>) to Eqs. (<ref>) and (<ref>)
is still not quite justified numerically. Nevertheless, the prediction of
a logarithmically increasing α(t), as confirmed by the simulations, is a
remarkable success of our simple scaling theory.

Interestingly, our simulations demonstrate that generally the observed subdiffusion
is much faster than the corresponding limit of disorder-renormalized normal
diffusion, which is also shown in Fig. <ref>a for several values of the
temperature. The presence of correlations thus leads to a dramatic increase in
the particle mobility on intermediate but relevant time scales. Note that whereas
for T=0.5 this limit is already gradually approached in Fig. <ref>, the
corresponding asymptotics cannot even be depicted for T=0.2 and T=0.1 in this
figure, as they lie outside of the plotting range, and, therefore, are completely
irrelevant on the corresponding time scale.

As discussed above, the physical origin of the observed subdiffusion is due to
a weak breaking of ergodicity. Hence, single trajectory time averages

    δ x^2(t)^ T_w=1/ T_w-t∫_0^ T_w-t[
    δ x(t|t')]^2dt'

of the mean-squared displacement δ x(t|t')=x(t+t')-x(t') over the time
window T_w  are expected to be very different from the above ensemble
result, and their amplitudes should be broadly scattered <cit.>.
This is indeed so, as evidence in Fig. <ref>a, where at all times t≪
T_w, and hence this scatter is not a trivial statistical effect occurring when
t∼ T_w. Remarkably, some of the particle trajectories are quickly
localized, while others are diffusing very fast. The slow ones start near to and
become trapped in a low potential valley, while the fast particles start from a
relatively high value of U(x) and move downhill to a much lower value of U(x),
thus experiencing a local energy bias. We note that this phenomenon is analogous to
the Golosov localization in standard Sinai diffusion <cit.>,
when particles starting at the same (thermal) initial position are not significantly
separated in the course of time. One may observe two such very close trajectories
in Fig. <ref>a. Thus, particles diffuse similarly and are correlated. This
feature can be very important for the diffusion of proteins on DNA, which may be
locally biased, even if the bias is absent on average. This behavior is very 
different from the scatter of single trajectory averages in the case of annealed
continuous time random walk subdiffusion with divergent mean resident time. In the
latter case, even identical particles starting at the same place will follow very
different, diverging trajectories. This feature can therefore be used for a crucial
experimental test to distinguish between different types of subdiffusion. 
Overall, both ensemble and
time averaged mean squared displacements are by many orders of magnitude faster
than in the limit of renormalized normal diffusion whose diffusion coefficient
is suppressed by the factor exp(-25)≈1.39× 10^-11, as compared
with free normal diffusion in Fig. <ref>a for the same temperature T=0.2.

As can be seen from single trajectory recordings (not shown), particles typically
continue their diffusive motion after being localized for a certain time. This
feature can also been seen in some trajectory averages depicted in Fig. <ref>a,
where the diffusional spread displays a step-like feature. Namely, the diffusional
spread continues after temporally reaching a plateau. In data analyses this
might mistakingly be attributed to a continuous time random walk subdiffusion
behavior. The ensemble average of the single trajectory time averages is also of
interest. This is how experimentalists often proceed to smooth out single trajectory
averages <cit.>. Such an ensemble averaged time average (EATA) is
also shown in Fig. <ref>. It is nicely fitted by a power-law time dependence.
Hence, on the level of this ensemble averaged time averaged mean squared
displacement the power law subdiffusive regime is established earlier.




All moments of the corresponding residence time distribution in any finite spatial
domain are finite. Moreover, the average residence time is much (for k_BT≪σ
orders of magnitude) smaller than the one expected from the normal diffusion
characterized the renormalized result D_ ren <cit.>. In the present case, weak
ergodicity breaking does not rely on infinite mean residence times, what makes it
especially interesting in biological applications. Residence time distributions of
escape times for particles initially located in the center of a symmetric interval
±λ in space for k_BT≤0.5σ are nicely described by a
generalized log-normal distribution reading

    ψ(t)   =   C/t[ e^-|ln(t/t_m)/κ_1|^b_1θ(t_m-t)
          +e^-|ln(t/t_m)/κ_2|^b_2θ(t-t_m)],

where C=b_1b_2/[b_2κ_1Γ(1/b_1)+b_1κ_2Γ(1/b_2)] is a
renormalization constant, b_1,2>1, and κ_1,2>0. Γ(x)
denotes the complete Gamma function. For any b_1,2>1, this distribution
(<ref>) has finite moments ⟨ t^n⟩ of all orders n. 
For b_1=b_2=2 and κ_1=κ_2, this is the well-known log-normal
distribution. However, in our case b_1 and b_2 are different, and κ_1
≠κ_2. Namely, for the exponential correlations b_1≈2.6÷2.7,
and it is weakly temperature dependent. However, b_2 does depend visibly on
temperature. For T=0.5, b_2=2, and it becomes smaller with decreasing
temperature while t_m becomes longer. For instance, for T=0.2, b_2≈
1.54, see Fig. <ref>a, where the distribution of times is plotted for
the log-transformed variable ln t, attaining a maximum at ln t_m.




 §.§ Power law correlations


We proceed with the case of power law correlations of the spatial potential
fluctuations, with γ=0.8. The ensemble averages for the mean squared
displacement are depicted in Fig. <ref>c. Here, the fit with δ=1
works generally better, which can be rationalized from Fig. <ref>, and
for T=0.1 a generalized Sinai diffusion with a≈3 nicely fits the
numerical data over six decades in time. Note in Fig. <ref>b that
α(t) stays nearly constant over five decades in time, up to the end
of the simulations even for T=1/3. For the larger T=0.5, it also grows
logarithmically in time, as in the case of exponential correlations.  Why
α(t) changes slower in time in this case and a fit with δ=1
and σ_ eff≈1.41 works nicely over many time decades for
intermediate temperatures can be rationalized from the fact that the
non-ergodicity length L_ erg in this case is much larger, see Table
<ref>. The case of power law correlations is especially important for
the diffusion of regulatory proteins on DNA strands, where such long-range
correlations with γ∼0.6÷0.8 emerge due to the way biological
information is encoded in the base pair sequence. Our results imply that
such diffusion should be typically anomalously slow with α∼ 1.41 k_BT/
σ in the corresponding range of temperatures, for γ=0.8. It must
be stressed in this respect that diffusion does not become immediately normal
for ⟨δ x^2(t)⟩>L_ erg^2, but rather a very slow
crossover with growing α(t) emerges, see Fig. <ref>c for T=0.5.
Here, the transition to the asymptotic regime of normal diffusion is still far
from being established.

Somewhat surprisingly, for power law correlations
subdiffusion is essentially faster in absolute terms than in the case of
exponentially short correlations, see Fig. <ref>. This may first appear
counter-intuitive. However, one reason for this feature becomes clear from the
potential fluctuations depicted in Fig. <ref>, which for the exponential
correlations are essentially larger (thinking in units of k_BT). Clearly,
exponential correlations and other singular (in the limit Δ x→ 0)
models of static disorder present a preferred case to observe Sinai type
diffusion in such systems. The second explanation is that in the case of
power law correlations a local bias is present on a much longer scale. This
leads to an essential acceleration of diffusion in the ensemble sense, see
Fig. <ref>b, where the ensemble mean squared displacement is in fact
superdiffusive (see also Fig. <ref>b for the initial regime) within several
correlation lengths at T=0.2. Remarkably, it is even faster than the potential
free diffusion. This is namely due to the presence a local bias, which is also
the reason for the Golosov phenomenon in Sinai diffusion. The single trajectory
time averaged mean squared displacement, however, does not show such a striking
feature, see Fig. <ref>b. The reasoning for this feature is that the local
bias is averaged out in the single trajectory time averaging for a sufficiently
large measurement time T_ w. Single trajectory averages spread out slower
than the ensemble averaged result and, nevertheless, their broadening is many
orders of magnitude faster than the result of the renormalized normal diffusion.
This is why even for σ∼ 4÷5, relevant for some important biophysical
situations <cit.>, regulatory proteins can diffuse on
DNA tracks despite of the fact that the classical de Gennes-Zwanzig-Bässler
result would predict that they should be practically localized (on biophysically
relevant time scales) for such a strong disorder. Naturally, single trajectory
averages exhibit a large scatter due to a transient lack of ergodicity. It is
important to mention that a strong single trajectory scatter was indeed observed
experimentally <cit.>. Note also that such a local bias can be functionally
very important, directing the protein toward a specific binding site on DNA.
The ensemble averaged time averaged mean squared displacement in this case also
shows a power law scaling in time, see Fig. <ref>b. Judging from the power
law exponent α_ EATA=0.31, it appears to be slower than in the case
of exponential correlations, where α_ EATA=0.54, compare with
Fig. <ref>a. However, the prefactor in this case is, in fact, much larger,
which makes it in fact faster rather than slower.

The distribution of escape times is also nicely captured by the generalized
log-normal distribution (<ref>), as seen in Fig. <ref>b. Note
that the escape times in this case are essentially shorter than in the case
of exponential correlations, Fig. <ref>a. Particularly, the scaling
exponent b_2 is larger and the half width κ_2 is smaller. Conversely,
the power exponent b_1 is pretty robust, b_1∼2.6 (also for other models
of disorder, see below).




 §.§ Linear correlations


The case of linear correlations presents another important model of singular
disorder with potential fluctuations growing with diminishing Δ x.
Therefore, it is expected to be similar to the case of exponential correlations,
despite the fundamental difference between these two models. The non-ergodicity
length L_ erg here is the shortest within the four considered models, see
Table <ref>. The ensemble averaged mean squared displacement is depicted
in Fig. <ref>d. Qualitatively, it appears very similar to the previous cases,
however, it turns out to be the slowest one, as Fig. <ref> reveals. In this
case, σ_ eff∼4.22÷4.91, and x_ in∼ 0.27÷0.67,
except from T=0.1, where x_ in≈3.65. Now, the agreement with x_0
≈0.314 is much better. It might seem paradoxical that the shortest
correlations, which exactly vanish at distances exceeding λ, yield the
slowest diffusion. However, it must be kept in mind that in the absence of
spatial correlation we have just normal diffusion, which is orders of magnitude
slower than the considered correlation-induced subdiffusion. The time behavior
of the anomalous diffusion exponent α(t) is indeed more similar to the
one in the case of exponential correlations than to the case of power law
correlations, see Fig. <ref>c, and compare with panels (a) and (b) therein.
The scatter of the single trajectory time averaged mean squared displacements
is also quite pronounced, see Fig. <ref>c. The escape times obey the same
distribution (<ref>), as depicted in Fig. <ref>c, but with different
parameters. Similar to previous cases, b_1∼2.6÷2.7, nearly independent
of temperature, while b_2 strongly depends on temperature, decreasing with
decreasing temperature. All temporal moments of the resident times remain,
however, finite with decreasing temperature.




 §.§ Gaussian correlations


Finally, we consider the case of Gaussian spatial correlations. Judging from
Fig. <ref> it should be the second fastest of the four models considered
here. This is indeed the case, as seen in Fig. <ref>. The ensemble averaged
mean squared displacements display very similar generic features as in the other
cases, see Fig. <ref>e. Here, σ_ eff∼4.50÷5.10 and x_
in∼0.41÷4.70, whereas the theoretical is x_0≈ 2.22. Once again,
the discrepancy by a factor of less than two with the theoretical σ_
eff=2.83, and at least the correct order of magnitude for x_ in are
quite impressive, given the simplicity of our scaling argumentation. Sinai type
diffusion with a=4 is featured in the low temperature behavior at T=0.1.
In this case, however, α(t) already starts to slightly increase at the
end of the simulations already at T=0.1, Fig. <ref>d. This type of
correlations indeed presents the worst case to observe a Sinai like diffusion.
The reasons are quite obvious: Namely, (1) this model of disorder is not
singular and (2) the correlations are short-ranged. Single trajectory averages
of the mean squared displacement are strongly scattered, and their ensemble average
yields a power law dependence for sufficiently large times, Fig. <ref>d.
Finally, the distribution of the escape times obeys the same common pattern,
as seen in Fig. <ref>d.






§ DISCUSSION AND CONCLUSIONS


Diffusion in systems characterized by Gaussian fluctuations of the potential
with spatially decaying correlations is typically thought of as normally
diffusive on experimentally relevant time scales, following the renormalization
idea for these systems developed by de Gennes, Zwanzig, and Bässler.
Following recent discoveries of extended anomalous diffusion in such systems,
we here provided a clear physical picture for the origin of subdiffusion in
stationary Gaussian random potentials with decaying spatial correlations. While
the anomalous diffusion is, of course, transient, the scaling theory developed
and confirmed herein demonstrates that this transient may easily reach macroscopic
time scales for a physically relevant strength of the potential
fluctuations.

The primary reason for the origin of this subdiffusion is that the maximal amplitude of
the potential fluctuations grows logarithmically with the distance, in accordance
with the asymptotic theoretical result (<ref>) and the simulations data
depicted in Fig. <ref>. This behavior is universal, and it has a predictive
power. In particular, it allows to predict for which model of the potential correlations
decay the anomalous diffusion will be faster in
absolute terms—and this prediction completely agrees with numerics. The
strongest trapping effects leading to anomalous diffusion occur in the case of
singular models of disorder, such as exponential or linear decays of the
correlations, for which the stationary energy autocorrelation function is not
differentiable at the origin. This leads to unbounded static force root mean
squared fluctuations ⟨ f^2⟩^1/2→∞ in the strict 
Δ x→ 0 limit of the coarse graining. Strictly speaking, such models are not
amenable to any Langevin simulations unless the random potential is considered
on a lattice with a spatial grid Δ x, as incorporated in our simulations,
with discrete potential points connected by parabolic splines (locally piece-wise
linear forces). A maximal time step of the Langevin simulations must be chosen
accordingly, depending on Δ x and the averaged local force root mean
square ⟨ f^2⟩^1/2. Singular disorder models are regularized
accordingly. However, they display the largest potential fluctuations in
Fig. <ref>.

Spatially increasing potential fluctuations may appear strange and contradictory
to the very fact that the random potential is stationary and possesses the well
defined root mean squared magnitude σ. This property of continuous albeit
logarithmic increase with distance is, however, the very cornerstone of the theory
of extreme events <cit.>: An extreme event will always occur, if only
one waits sufficiently long, and even if the probability density of such events
is strongly localized such as in the case of the considered Gaussian distribution.
In our case, the extreme events occur in space rather than time, but have the
same origin. The law (<ref>) of potential fluctuations with a (log x)^δ scaling (δ=1/2) is valid asymptotically. The longer the spatial
correlations of the potential fluctuations reach, the later this asymptotic is
achieved. Numerically, we also observed different scaling exponents such as
δ=1 transiently and δ>1 initially.

Based on these observations we put forward a very simple scaling theory of
anomalous diffusion in the potential landscape U(x), which in essence is
very similar to the scaling theory of Sinai diffusion developed in the
classical Ref. <cit.>. It leads to our major result given in
Eq. (<ref>) and explains why for sufficiently low temperatures we
indeed observe a generalized Sinai diffusion with a power law exponent, which
is generally different from a=4 and can be even smaller as shown in
Figs. <ref>b, c, and d: The diffusion is formally even slower than for the
case of Sinai diffusion in random force fields. Our analytical results were
shown to be valid over five to seven decades in time, at least, which is a
remarkable success of the relatively simple scaling approach. The nominally
ultraslow diffusion is in fact many orders of magnitude faster than the
asymptotical limit of the disorder-renormalized normal diffusion, which for
the physically relevant parameters considered here simply cannot be attained
neither physically nor numerically for such a strong disorder. Hence, it becomes
completely irrelevant in such settings. This provides a striking example of how
the formally valid mathematical result of asymptotic, disorder-renormalized
normal diffusion in the de Gennes-Zwanzig-and Bässler theory may produce
physically inappropriate descriptions on meso- and macroscales.

This result also leads to a number of further predictions, which are confirmed
in numerical simulations and agree with previous simulations of the considered
stochastic dynamics in random potentials. Namely, this kind of anomalous
diffusion is generally characterized by a time dependent power law exponent
α(t), for which we obtained the theoretical result (<ref>).
This result predicts that, after an initial decay (δ>1), the anomalous
diffusion exponent α(t) will reach a minimum for which α_
min=2k_BT/σ_ eff (for δ=1) and then gradually grow in
logarithmic fashion, α(t)=2(k_BT/σ_ eff)^2ln(t/t_0)
(δ=1/2). These two nontrivial predictions are remarkably well confirmed
by extensive numerical simulations. Indeed, α_ min=2k_BT/σ_
eff with σ_ eff∼1.24÷1.52, within a broad range of
temperatures, as demonstrated in Fig. <ref>. Moreover, at sufficiently large
temperatures, α(t) indeed grows logarithmically, see Fig. <ref>. It
is indeed remarkable how a simple scaling theory can be predictive to such a
degree. Of course, the scaling results cannot be fully accurate on a quantitative
level and claim full consistency over the entire range of parameters and times
studied herein. In particular, the values of σ_ eff and t_0 used
to fit the numerical data, for instance, in Figs. <ref> and <ref>
are in fact different. This is due to the fact that the approximation of
Eq. (<ref>) by Eqs. (<ref>) and (<ref>) is, indeed, a bit too coarse.
Nevertheless, the coherence and validity of our scaling theoretical predictions
are compelling.

We note that in the biologically relevant and important case of power law decaying
correlations the scaling δ=1 applies to a much larger range of times and
temperatures, as shown in Figs. <ref>c and <ref>b. In this case, power
law subdiffusion of the form ⟨δ x^2(t)⟩∝ t^α with
α≈1.41  k_BT/σ in fact presents a good approximation over a broad
range of times and temperatures. This behavior may easily be confused with the law
α≈ k_BT/σ produced by the mean field continuous time random
walk with exponential energy disorder <cit.>. In
striking contrast to the latter model, however, the residence times in a finite
spatial domain are not power-law distributed but follow the generalized log-normal
distribution (<ref>) shown in Fig. <ref>, for all the studied models
of the correlation decay. All moments of this distribution are finite, which is a
very attractive physical feature. This subdiffusion is clearly non-ergodic, as it
is demonstrated in Fig. <ref>, where the maximal time is merely 1% of the
time window T_ w used to accumulate the single-trajectory averages. Hence,
a large scatter is a truly non-ergodic effect, even if it is a transient one and
must vanish in the mathematical limit T_ w→∞ <cit.>. It
is, however, not necessarily attainable experimentally.

The origin of this lack of ergodicity was already explained in
Ref. <cit.>: It is due to the absence of self-averaging of the
statistical weight function w(x)=exp[±β U(x)] on the spatial scale
L_ erg, implicitly defined in relation (<ref>), depending on the
disorder autocorrelation function g(x). Solving this equation analytically
in the case of linearly decaying correlations yields the exact result
(<ref>) and a very handy and highly accurate approximate result in
Eq. (<ref>). In this case of linear decay, L_ erg turns out to
be the shortest one within all models of g(x) considered in this paper, see
Table <ref>. This result implies that L_ erg in units of the
correlation scale λ is only by a factor (1/2)(σ/k_BT)^2 smaller
than the factor exp[σ^2/(k_BT)^2] suppressing the asymptotically normal
diffusion coefficient in the renormalization sense. This, in turn, means that
such mesoscopic subdiffusion can readily reach macroscopic scales even for a
moderate disorder strength σ∼ 4÷ 5  k_BT featuring many physical
and biological systems already at room temperatures.

Different from the annealed continuous time random walk subdiffusion,
the single trajectory time averages for the mean squared displacement in
the present case of quenched potential fluctuations are characterized by
different power law exponents α≠1 <cit.>, and not
just a linear scaling with the lag time (erroneously suggesting normal
diffusion) but with a significantly scattered diffusion coefficient
<cit.>. Importantly, especially with respect to
biological applications is that the subdiffusion considered here is in fact
orders of magnitude faster than the normal diffusion predicted by the classical
renormalization result with effective diffusivity D_ ren. For example,
for the diffusion of regulatory proteins on DNA tracks,  σ≈
4.3 k_BT_ room, which can be deduced from the results presented in
Ref. <cit.>. Hence, D_ ren=9.33× 10^-9D_0, and for
a typical experimental value D_0=3 μ m^2/s <cit.>
it would become D_ ren≈ 2.8× 10^-2 nm^2/s. This
would mean that within some hundreds of seconds the diffusional spread would
be merely several nanometers. However, this renormalization result 
underestimates massively, on the relevant time scales of this phenomenon, 
the actual protein mobility. Our results demonstrate that the actually
occurring subdiffusion is orders of magnitude faster than one suggested by
this normal diffusion result of the renormalization approach. Hence,
correlations-induced persistent subdiffusion makes diffusional search
feasible in such a situation despite a strong binding-energy disorder.

An important role in applications may also be played by the presence of a local
bias, which is especially clearly expressed in the case of power law correlations.
This is the analogous reason for the Golosov phenomenon in the case of random
potential exhibiting Brownian motion in space <cit.>. Then,
the particles which were initially localized nearby diffuse similarly, in a
correlated fashion. The distance between them does not grow dramatically in
time, being bounded. One can find similar single trajectory averages also in
our numerical results. Conversely, different particles starting in locally
different environments can move into opposite directions, and this can give
rise to an enhanced diffusivity in the ensemble sense. In fact, due to this
reason an initial regime of superdiffusion can be realized in Fig. <ref>b
and d on the ensemble level in the case of power-law correlations, for which
it can extend over several scaling lengths λ, see also
Fig. <ref>b and d for small times. In the case of single trajectory time
averages, the local bias averages out. Therefore, such a superdiffusive regime
is absent. Diffusion on the level of single trajectories is typically slower
than the one on the ensemble level, see, for instance, Fig. <ref>.
The difference between the ensemble average and the ensemble-averaged time average
of the mean squared displacement becomes smaller with increasing time. The latter
ensemble-time-averages display a power law dependence on time for sufficiently
large times even in the Sinai like regime on the ensemble level, which can be
an important observation with respect to possible experimental manifestations.

To conclude, we elucidated the physical mechanism leading to subdiffusion in
stationary correlated potentials with spatially decorrelating Gaussian disorder,
and we showed that a generalized Sinai diffusion typically emerges at sufficiently
low temperatures and/or strong disorder for various models of decaying
correlations. Our scaling theory also explains how a standard power law
subdiffusion emerges with increasing temperature and in the course of
time. Such subdiffusion is weakly non-ergodic, displays a local bias and
proceeds much faster than de Gennes-Bässler-Zwanzig limit of normal diffusion,
which for sufficiently low temperatures and/or finite size of the system
simply cannot be attained physically on typical mesoscales. We believe that
our results provide a new vista on the old problem of potential disorder, and
that they will be very useful in the context of non-ergodic diffusion processes,
especially in relation to various biologically relevant problems on the cellular
level. Likely they are also important for diffusion of colloidal particles in
laser created random potentials, a conjecture calling for further experimental
studies of such systems.



§ ACKNOWLEDGMENT
 
Funding of this research by the Deutsche Forschungsgemeinschaft (German
Research Foundation), Grant GO 2052/3-1 is gratefully acknowledged.

PRX


