Subdiffusion via dynamical localization induced by thermal equilibrium fluctuations

We reveal the mechanism of subdiffusion which emerges in a straightforward, one dimensional classical nonequilibrium dynamics of a Brownian ratchet driven by both a time-periodic force and Gaussian white noise. In a tailored parameter set for which the deterministic counterpart is in a non-chaotic regime, subdiffusion is a long-living transient whose lifetime can be many, many orders of magnitude larger than characteristic time scales of the setup thus being amenable to experimental observations. As a reason for this subdiffusive behaviour in the coordinate space we identify thermal noise induced dynamical localization in the velocity (momentum) space. This novel idea is distinct from existing knowledge and has never been reported for any classical or quantum system. It suggests reconsideration of generally accepted opinion that subdiffusion is due to broad distributions or strong correlations which reflect disorder, trapping, viscoelasticity of the medium or geometrical constraints.

where the averaging is over all thermal realizations as well as over initial conditions. The diffusion process can be classified through the scaling function 5 σ ∼ .
x 2 The normal diffusion corresponds to the scaling index α = 1. Any deviation from this linear time dependence is classified as anomalous diffusion. For the superdiffusive case, σ t ( ) increases over time faster while for the subdiffusion it grows slower than for normal diffusion. An example of the former is ballistic diffusion with the scaling index α = 2. The hallmark of the latter is famous Sinai subdiffusion which follows the logarithmic law t t ( ) ln σ ∼ 6 . This ultraslow process can be observed for a Brownian particle moving in a static random Gaussian force field imitating quenched disorder in heterogeneous media. "Quenched" means that random traps, barriers or comb-like structures do not evolve with time. This is usually the model considered to describe the dynamical properties of materials containing impurities, defects, or intrinsic randomness like it is the case for amorphous systems 7 . However, recent progress in single particle tracking techniques 8 has allowed to probe transport processes occurring in more complex setups. For instance, the diffusive motion of macromolecules and organelles inside living cells is typically subdiffusive 9,10 . This behaviour is commonly attributed to macromolecular crowding of their interior, summarizing their densely packed, heterogeneous and fluctuating environment 11, 12 . In this new class of systems subdiffusion may be only a transient effect while normal diffusion is observed in the asymptotic long time regime 13 . Nevertheless this anomaly lasts sufficiently long for experimental detection [14][15][16][17] . Despite the outlined fundamental differences researchers try to exploit the well known mechanisms to capture the essence of anomalous diffusion (in the biological context, see e.g. ref. 18 ). They are either broad distributions or strong correlations in diffusive motion. Among others there are two major subdiffusive models. The first one is a system where the particle dynamics is found to be governed by a sequence of trapping-untrapping events, e.g. in energy wells. It is described in the framework of continuous time random walk 5 consisting of jump models where the particle undergoes a series of displacement given in terms of distribution of waiting times with power law tails. The second one is a system where the particle does not simply move in a fixed potential as before, but is a part of an interacting setup exhibiting viscoelastic behaviour meaning that dynamics of different components of the system are correlated. This situation is described in terms of fractional Brownian motion or a fractional Langevin equation [19][20][21] by relaxing the white noise assumption in simple Brownian motion and a priori inclusion of a power law for time correlations of thermal noise. A more interesting situation occurs when this diffusion anomaly is induced by dynamics of the problem itself and is not related to genuine disorder or is not introduced ab initio with broad distributions or strong correlations. We can mention models of anomalous diffusion (both subdiffusion and superdiffusion) in discrete deterministic systems 22,23 . In contrast, here we use a bottom up approach: by analysing the Langevin equation for a Brownian ratchet driven by thermal white noise, we reveal a new mechanism how subdiffusion may emerge from deterministic dynamics superimposed with thermal noise.

Results
We consider an archetype of the Brownian motor 24,25 consisting of an inertial Brownian particle moving in a one-dimensional periodic potential of a ratchet type (i.e. with broken reflection symmetry) and driven by an external unbiased time-periodic force. Its dynamics is determined by the following dimensionless Langevin equation 26m We refer the reader to the section Methods where we describe in detail the scaling procedure. The dot and the prime denote differentiation with respect to time t and the Brownian particle coordinate x ≡ x(t), respectively. The dimensionless friction coefficient is 1 and the parameter m is the dimensionless mass of the particle moving in a spatially periodic potential U(x) = U(x + L) of period L and driven by the external unbiased time-periodic deterministic force F(t) = a cos (ωt) of amplitude a and angular frequency ω. Thermal fluctuations due to coupling of the particle with the thermal bath of dimensionless temperature θ are modelled by δ-correlated Gaussian white noise ξ(t) of zero mean and unit intensity, i.e., The potential is assumed to be in the ratchet (asymmetric) form 24 of the period L = 2π. Due to the presence of the external driving F(t) = a cos (ωt) and the friction term  x the particle velocity approaches for t → ∞ a unique nonequilibrium stationary state which is characterized by a temporally periodic probability density. Then the mean velocity  〈 〉 x t ( ) takes the form of a Fourier series over all possible harmonics 27 where 〈v〉 is the directed (time independent) velocity while v nω (t) denote harmonic functions of vanishing average over the fundamental period T = 2π/ω. Due to this particular decomposition it is useful to study the period averaged velocity v(t) defined as which may be exploited to evaluate the directed velocity as . A sufficient and necessary condition for the emergence of the directed transport 〈v〉 ≠ 0 is breaking of the mirror symmetry of the potential U(x) which is the case for the form described by Eq. (5) 24 . This operating principle can be seen as a key for understanding the intracellular transport 28 . Despite of simplicity of this system, it exhibits a number of notable and unusual features [29][30][31][32][33][34] including the anomalous diffusion 26,35 . The origin of observed superdiffusion has already been satisfactorily explained in ref. 26 . Unfortunately, the mechanism standing behind subdiffusion still lacks a rewarding illumination. Within this work we aim to fill this major gap.
In panel (a) of Fig. 1 we depict the time-dependent "diffusion coefficient" D(t) defined as shaded with the cyan colour and the asymptotic long time regime of normal diffusion for t > τ 2 . The crossover times τ 1 and τ 2 separating these domains can be controlled by temperature and are monotonically decreasing function of thermal noise intensity 26 . For example, when θ = 0.0007 (the blue curve in Fig. 1), τ 1 ≈ 3.2 · 10 3 and Scientific REPORts | 7: 16451 | DOI:10.1038/s41598-017-16601-0 τ 2 ≈ 10 8 . If temperature is lowered to θ = 0.00016 (the red curve in Fig. 1) the superdiffusion lifetime is extended to τ 1 ≈ 3.2 · 10 6 . It is difficult to numerically determine τ 2 due to limited stability of the utilized algorithm leading to uncontrolled propagation of roundoff and truncation errors. However, if we adopt a more reserved extrapolation from the case θ = 0.0007, the time τ 2 is at least of order 10 11~1 0 13 . So, these two times τ 1 and τ 2 are many, many orders longer than characteristic time scales of the system which in the dimensionless form are of the order of unity! Therefore from the experimental point of view for tailored parameter regimes these diffusion anomalies may be safely treated as nearly persistent effects.
To explain the mechanism standing behind the subdiffusion in such a periodic system let us now study the period averaged velocity v(t) of the Brownian particle, c.f. Eq. (7).
characterises fluctuations of the actual period averaged velocity around its mean value and is depicted in Fig. 1(b) for selected values of temperature θ. On its basis, we can make two observations. The first is that the crossover time τ 1 of superdiffusion is related to the relaxation of the period averaged velocity v(t) to its stationary state which for low to mid temperature may be ultraslow 26 . The second is that when the particle subdiffuses then the variance of the period averaged velocity is much smaller than it is for superdiffusion. This fact suggests that during subdiffusion the coherence of motion is enhanced not only in the coordinate space but also in the velocity space. It is better visualized in Fig. 2(b) where we present the asymptotic value of the period averaged velocity variance σ v computed for the final moment of time t f = 10 6 in our numerical simulations as a function of temperature θ. The interval corresponding to subdiffusion is shaded with the cyan colour, c.f. panel (a) of the same figure where we depict the scaling index α fitted to the asymptotic parts of the mean square deviation of the particle coordinate t ( ) evolved up to the same t f . The variance is more than three times smaller in this temperature window than for values describing superdiffusion and normal diffusion. Moreover, the scaling index α can be controlled by temperature of the system in the non-monotonic way. There is a window of thermal noise intensity θ ∈ [7 · 10 −4 , 8.5 · 10 −4 ] for which it is very small but still non-zero α ≠ 0 indicating nearly coherent transport or the ultraslow subdiffusion. For slightly lower and higher temperature the asymptotic evolution of σ t ( ) is essentialy subdiffusive, e.g. α = 0.5 for θ = 1.175 · 10 −3 , c.f. Fig. 2(a). These findings suggest the prominent role which is played by thermal noise in the observed subdiffusive behaviour.
We further examine it by constructing a rough-and-ready three state stochastic process. In the deterministic limit θ = 0 the system is non-chaotic and possesses three coexisting attractors for the period averaged velocity: v − ≈ −0.4, v 0 ≈ 0 and v + ≈ 0.4 26 . Therefore we now analyse transitions between these states induced by thermal fluctuations. We introduce the following notation: p ++ stands for the (conditional) probability to remain staying in the plus state v + → v + , p +− denotes the probability for a jump between the opposite states v + → v − and p +0 is the probability of a transition between the plus and the zero state v + → v 0 . This convention is analogous for the remaining six probabilities {p 00 , p 0+ , p 0− , p −− , p −0 , p −+ }. By introducing the threshold ν = 0.2 the above probabilities can be estimated from simulations of the Langevin equation (3) as the relative frequencies with which the period averaged velocity v(t) visits the three coarse grained regions V + = {|v(t) − 0.4| < ν}, V 0 = {|v(t)| < ν} and V − = {|v(t) + 0.4| < ν}. A chosen value of the threshold ν = 0.2 follows naturally from inspection of the probability distribution of the period averaged velocity where for low to moderate temperature regimes there are three peaks corresponding to the deterministic coexisting attractors which are approximately ν = 0.2 apart (not depicted). Figure 2(c) presents all transition probabilities as a function of temperature θ. In the low temperature regime the probabilities for the particle to survive in the opposite running states p ++ and p −− corresponding to the velocities v + ≈ 0.4 and v − ≈ −0.4, respectively, are particularly large and almost equal. On the other hand, the probability for staying in the locked state p 00 is twice lower. Consequently, the spread of trajectories is large and ballistic motion is observed. The temperature interval for which subdiffusion is developed is again marked by the cyan colour. There thermal noise induces dynamical localization of the period averaged velocity (momentum), i.e. it resides in the positive running state v + with the probability p ++ very close to unity. It means that once the ensemble of is presented for three values of thermal noise intensity θ proportional to temperature. The region corresponding to the subdiffusive behaviour for θ = 0.0007 is indicated with the cyan colour. Parameters are m = 6, a = 1.899, ω = 0.403. At zero temperature θ = 0, the system is nonchaotic.
particles enter this state in the velocity space it moves almost coherently, i.e. with marginal fluctuation allowed in the chosen threshold |v(t) − 0.4| < ν. These fluctuations are responsible for the ultraslow subdiffusion where the observed scaling index is very small but still nonzero α ≠ 0. We also note that in this interval the probability for surviving in the locked state p 00 corresponding to the particle locked in one of the potential wells is enlarged. Each stay of the particle in the locked state hampers the increase of the spread of trajectories leading to slowing down of diffusion. Similar mechanism contributes to subdifussion observed in systems with disorder, e.g. random barriers, where the particle dynamics is governed by trapping-untrapping events in energy wells. However, certainly the dominant factor which is responsible for the detected subdiffusion is dynamical localization in the positive running state as it is depicted in Fig. 2(c). Further increase of temperature immediately delocalizes the particle. However, since then the probability for being in the locked state p 00 is at the level of p −− and noticeable smaller than p ++ we observe the occurrence of normal diffusion. In the high temperature limit all transition probabilities are almost the same.
To better visualize the explained mechanism we now focus on the single representative Brownian particle trajectory for temperature θ = 0.0004 corresponding to subdiffusive behaviour. The result is shown in Fig. 3. In panel (a) and (b) we present the evolution of the coordinate x(t) and the period averaged velocity v(t), respectively. Each red dot in panel (b) depicts the latter quantity for a given t. The time intervals corresponding to motion in each of the observed states v − , v 0 , v + can be easily identified. Since in this regime p ++ ≈ 1 once the particle enters the plus state it stays in its vicinity v(t) = v + ≈ 0.4 (c.f. panel (b)) for very long time. After the dynamical localization all particles travel with almost identical velocity and the spread of their trajectories changes very little over time implying the subdiffusive behaviour. We illustrate this mechanism in the inset of Fig. 3(a) where we depict an ensemble of 100 trajectories of Brownian particle dynamics for the same temperature. However, we stress that the magnitude of the estimated probability for surviving in the plus state p ++ naturally depends on the chosen threshold ν of the velocity space coarse graining procedure. We restricted ourselves to the simplest stochastic model capturing the essence of the mechanism standing behind the observed subdiffusion. Summing up this part, for the fixed temperature θ and small times t < τ 1 the dynamics of each trajectory is mostly deterministic giving rise to the (ballistic) superdiffusive behaviour provided that it is averaged over unbiased initial conditions for the particle position and velocity. The escape rates from the minus and zero attractors are related to the crossover time τ 1 of the superdiffusive stage of motion. For intermediate time regime τ 1 < t < τ 2 thermal fluctuations induce transitions from the minus and the zero states onto the plus velocity. The observed dynamical localization in the plus state means that the typical escape rate from this solution -which is related to the crossover time τ 2 of subdiffusion -is much larger than the corresponding quantity for two other attractors thus eventually luring there almost all trajectories and hindering the diffusion rate. We numerically verify that for temperature θ ≈ 0.0007 for which the power exponent α ≈ 0 in the interval τ 1 < t < τ 2 the probability for the particle to be in the plus state is very, very close to one p + ≈ 1 meaning that almost all particles are localized in this solution. Finally, for sufficiently long times t > τ 2 random dynamics induced by thermal fluctuations activates escape from all attractors and jumps between various trajectories leads to normal diffusion.

Discussion
In conclusion, we explained the mechanism of long lasting subdiffusion in a nonequilibrium noisy system, which in the deterministic limit is non-chaotic. As the reason for this anomalous behaviour we identified classical dynamical localization 36 of the particle velocity (momentum) induced by thermal equilibrium fluctuations, i.e. there is a temperature window for which the particle localizes in the positive running state v + with the probability very close to unity. This temperature interval matches exactly the region where the subdiffusive regime is observed in the finite data acquisition time.
Dynamical localization is a well known phenomenon in quantum systems and was introduced in ref. 37 . It manifests itself in quantum suppression of the chaotic classical diffusion in momentum space due to interference effects 38 . A different but closely related phenomenon to dynamical localization in deterministic systems is Anderson spatial localization in disordered systems 39,40 . With this work we established a relation between subdiffusion in the coordinate space and the dynamical (quasi)localization in the momentum space. To the best of our knowledge, a similar idea has never been reported for any other classical or quantum systems. The findings are distinct from existing knowledge and suggest reconsideration of generally accepted views that subdiffusion is due to broad distributions or strong correlations which reflect disorder, trapping, viscoelasticity of the medium or geometrical constraints.
Since we demonstrated the new mechanism of subdiffusion using the model of a Brownian ratchet, the results provide essentially new insight into processes governing transport in complex systems especially of biological origin. Notably, they may explain certain aspects of strange kinetics occurring in living cells, with a particular emphasis on subdiffusive motion of molecular motors which are responsible for the intracellular transport. While we illuminate this phenomenon in the case of the driven Brownian particle moving in the asymmetric (ratchet) potential, in principle there are no restrictions for this universal mechanism to be likely observed also for other nonequilibrium setups: driven or non-driven but tilted symmetric periodic systems 31,32 and including those which operate in strong disspation regime or in overdamped limit. The latter is especially true due to the fact that in the deterministic limit our setup is in a non-chaotic regime. Therefore complex chaotic dynamics which is characteristic in one dimension for driven inertial systems is not needed for this mechanism to occur. Moreover, due to simplicity and universality of the system with its physical clarity as well as appealing strength of Brownian motion with its intrinsic Gaussian noise propagator our findings can open a wide area of studies and may be corroborated experimentally with a wealth of physical systems. One of the most promising setups for this purpose are optical lattices [41][42][43][44] and asymmetric SQUID devices 45,46 .

Methods
In this work we analyse the generic model of a ratchet system which consists of (i) a classical inertial particle of mass M, (ii) moving in a deterministic asymmetric ratchet potential U(x), (iii) driven by an unbiased time-periodic force Acos(Ωt) of amplitude A and angular frequency Ω, and (iv) subjected to thermal noise of temperature T. The corresponding Langevin equation reads 26 The parameter Γ stands for the friction coefficient and k B is the Boltzmann constant. Thermal equilibrium fluctuations are modeled by δ-correlated, Gaussian white noise ξ(t) of zero mean and unit intensity, i.e., The spatially periodic potential U(x) is assumed to be in a double-sine form of period 2πL and a barrier height ΔU, namely As only relations between scales of length, time and energy are relevant but not their absolute values we next formulate the above equations of motion in its dimensionless form. To do so, we first introduce the characteristic dimensionless scales for the system under consideratioñx so that the dimensionless form of the Langevin dynamics (3) reads Here, the dimensionless potential ) . The dimensionless noise intensity Q = k B T/ΔU is the ratio of thermal and the activation energy the particle needs to overcome the nonrescaled potential barrier.
The system described by Eq. (9) possesses four characteristic time scales~~τ The quantity τ 0 is used as a unit of time in our scaling procedure, see Eq. (12). It is a characteristic time for an overdamped particle to move from the maximum of the potential U(x) to its minimum. The scale τ 1 is a relaxation time of the velocity of the free Brownian particle (when U(x) = A = 0). Note that the dimensionless mass m = τ 1 /τ 0 is a ratio of the two characteristic times. The time τ 2 is a characteristic scale for the conservative system (when Γ = A = 0). The last scale τ 3 is a period of the external time-periodic force. Thermal fluctuations are modelled here approximately as white noise so its correlation time is zero and there is no characteristic time scale associated with it. However, in real systems it is non-zero but much, much smaller than the other characteristic time scales.
The Langevin equation (9) has been originally derived for the asymmetric superconducting quantum interference device (SQUID) which is composed of a loop with three capacitively and resistively shunted Josephson junctions, see Eq. (14) in ref. 45 . The particle coordinate x and velocity v corresponds to the Josephson phase and the voltage drop across the device, respectively. The particle mass stands for the capacitance of the SQUID, the friction coefficient translates to the reciprocal of the SQUID resistance. The time-periodic force corresponds to the external current. The asymmetry parameter ϕ of the potential (11) can be controlled by the external magnetic flux which pierces the device.
The Langevin equation (13) is solved by employing a weak version of the stochastic second-order predictor corrector algorithm and using a CUDA environment implemented on a modern desktop GPU. This procedure allows for a speedup of a factor of the order 10 3 times as compared to a common present-day CPU method. Details on this implementation can be found in 47 . Since Eq. (13) is a second-order differential equation we need to specify two initial conditions: one for the position x(0) and the other one for the velocity  x(0) of the Brownian particle. For some regimes of system parameters the studied dynamics may depend on a specific choice of these initial conditions and therefore to avoid this unwanted behaviour we chose x(0) and x(0)  to be equally distributed over the intervals [0, 2π] and [−2, 2]. However, we note that the effects discussed in the paper are robust with respect to variation of these initial conditions. In particular, in ref. 26 we have considered three various forms of the probability distribution for the initial velocity of the particle: the Dirac delta P(v) = δ(v), uniformly distributed P(v) = U(−2, 2) on the velocity interval [−2, 2] and normally distributed (Gaussian) P(v) = N(0, 1) of zero mean and standard deviation equals to one. All considered forms lead to the consistent results proving that the observed behaviour is robust with respect to the probability distribution for the initial velocity of the particle.