Flashing subdiffusive ratchets in viscoelastic media

We study subdiffusive ratchet transport in periodically and randomly flashing potentials. Central Brownian particle is elastically coupled to surrounding auxiliary Brownian quasi-particles which account for the influence of viscoelastic environment. Similar to standard dynamical modeling of Brownian motion, the external force influences only the motion of central particle not affecting directly the environmental degrees of freedom (see video). Just a handful of auxiliary Brownian particles suffice to model subdiffusion over many temporal decades. Time-modulation of the potential violates the symmetry of thermal detailed balance and induces anomalous subdiffusive current which exhibits a remarkable quality at low temperatures, as well as a number of other surprising features such as saturation at low temperatures, and multiple inversions of the transport direction upon a change of the driving frequency in nonadiabatic regime. Our study generalizes classical Brownian motors towards operating in sticky viscoelastic environments like cytosol of biological cells or dense polymer solutions.


Introduction
Physics of noise-assisted driven transport presents currently a well-established area of research [1,2]. Most papers are devoted to classical transport in neglecting non-Markovian and even inertial effects. The corresponding stochastic nonlinear dynamics can be described by a Langevin equation in time-dependent potentials and/or by associated Fokker-Planck equation for the noise-averaged dynamics of the probability density of an ensemble of moving particles. Key ingredients are nonlinear dynamics in periodic potentials unbiased on average, friction and thermal noise related by the fluctuation-dissipation theorem (FDT), and an external driving which violates the symmetry of thermal detailed balance ensured by the FDT at thermal equilibrium. The emerged dissipative out-of-equilibrium directed transport is necessarily accompanied by an entropy production and related heat dissipation.
Such profoundly out-of-equilibrium transport should be distinguished from other possibilities such as transport in a running potential, where the particle remains bound to a potential well which moves in space. Basically, in this later case one does not need even a periodic potential. Any trapping confining potential U(x) with a deep minimum at x = x 0 will convey a transport if this minimum is moving with velocity v, i.e. U(x) is replaced by U(x − vt). Such a "peristaltic" transport is clearly not related in essence to any entropy production and for a periodic running potential it can be symbolized by the Archimedean screw pump. This is a purely mechanical system, where both the friction and the noise are not principal for the transport occurrence. Quantummechanical counterparts of the Archimedean pump are also well-known [3]. A peristaltic pump can be also realized with highly overdamped isothermic thermodynamic systems operating close to thermal equilibrium, where dissipation keeps the particle fluctuating near to the potential minimum, and the energy lost due to friction is perpetually restored due to stochastic impact of the environment -the physical content of FDT. Friction plays here in fact a constructive role. Moreover, neglect of the energy gain from hot environment, at odds with FDT, is one of the common mistakes in the literature, leading to erroneous belief that one must always minimize dissipation and cool the environment in order to achieve at the highest efficiencies of isothermal engines possible. This is not necessarily so. For example, the so-called dissipationless Hamiltonian ratchets [4] cannot do any useful work at all against a load, i.e. are pseudo-ratchets with zero efficiency. A dissipative system can remain locally close to thermal equilibrium, but the form of the potential is cyclically and adiabatically slow changed by a driving force so that excess, uncompensated heat exchange between the particle and the thermal reservoir can be minimized. In this quasi-equilibrium scenario, the work required to change the form of potential (to move the minimum) can be transformed with minimal heat losses into the work against the load which opposes this potential modulation. This is why such an isothermal machine considered as a free energy transducer can in principle operate with the efficiency, defined as the ratio of the work against a load to the free energy spent to drive a working cycle, close to one [5][6][7][8]. For example, the efficiencies of highly a) xx xx xx xx xx xx x Figure 1. Snapshots of numerical simulations of the dynamics in Eq. (13) at different instants of time a) t = 0.1, b) t = 0.4, c) t = 1.6, d) t = 2.0, e) t = 2.5. One Brownian particle (big filled circle) is coupled to auxiliary particles (small empty circles) at temperature T = 0.1 and moves in a periodically flashing with frequency ν = 1.0 potential (15) with amplitude U 0 = 0.75. optimized biological ionic pumps as high as 0.75 are common [9] and the efficiency of ATP synthase can approach the theoretical maximum of one [10][11][12].
The primary focus of this paper is different, on the thermal noise assisted transport, where the thermal noise is assumed to play a profound and constructive role. A paradigmatic example is provided here by transport in flashing potentials [13], such as one in figure 1. In a standard Markovian overdamped setup, when the potential is off, an initially localized particle diffuses with the position variance growing linearly, δx 2 (t) ∝ t, which corresponds to normal diffusion. When the potential is on, the particle relaxes to a minimum of the potential and the probability distribution becomes at least bimodal with a larger peak corresponding to sliding down a less steep side of the potential (as a larger basin of attraction corresponds to a larger distance between the minimum and maximum of a spatially asymmetric but periodic potential). When flashing repeats, either periodically or stochastically, but sufficiently slow, the net transport is expected to emerge in the left, natural direction, with the averaged particle's position growing linearly in time. Such normal transport is characterized by mean velocity. This natural direction is opposite to one in the fluctuating tilt potential ratchets of the same potential form [14,15]. For a fixed unbiased on average potential, the total heat exchange between the particle and its environment is zero, FDT holds, and directed transport is forbidden by the symmetry of thermal detailed balance. Out-ofequilibrium potential fluctuations violate this symmetry and induce directed transport. There emerges an overall uncompensated excess heat flow to the environment associated with the corresponding entropy production. A part of the energy put in the repeating potential flashing is dissipated as excess heat and a part can be used to do useful work against a load. If the load is absent all the consumed energy is dissipated as excess heat (futile motor) because the mechanical energy of the motor particle remains on average not changed. This is a well-established by now physical picture [16]. This basic model explains e.g. operating of single-headed kinesin motors [9], where the energy to drive stochastic cycles is provided in effect by the energy of ATP hydrolysis.
Generalization of this approach to account for a non-Markovian friction with memory is not trivial. The memory effects appear e.g. due to viscoelasticity of the environment [17][18][19][20][21][22][23]. Even if a corresponding Generalized Langevin Equation (GLE) [24][25][26][27] is well known for any linear model of friction with the memory, the corresponding non-Markovian Fokker-Planck equation (NMFPE) [28][29][30][31][32] remains simply unknown for general nonlinear force-fields. The cases of constant force, or a linear in coordinate force, where the corresponding NMPPEs are known [28], are not especially useful in the present context. Especially interesting is the case of a power law decaying memory kernel [18] which corresponds to subdiffusion and is associated with the Cole-Cole dielectric response of viscoelastic media [33,34]. Such viscoelastic subdiffusion has recently been found relevant also for transport in polymer networks [35,36], cytosol of biological cells [37] and akin crowded fluids [38]. The corresponding fluctuating tilt or rocking ratchets in viscoelastic media have recently been proposed and studied in Refs. [39,40] with a number of quite unexpected and surprising properties revealed. The interest to ratchet effect in glass-like environments is growing [41]. Flashing ratchet transport in dynamically disordered potentials was proposed and studied earlier in Ref. [42]. Transport in disordered potentials is known to be equivalent within the mean field effective medium approximation to a continuos time random walk (CTRW) approach with independent residence times in traps [43][44][45]. This places our research into a general context of anomalous transport processes in complex media [43][44][45][46]. However, our approach to anomalous transport in viscoelastic media is very different from one based on uncorrelated CTRW. This work is devoted to flashing viscoelastic ratchets which turn out to be no less surprising than rocking viscoelastic ratchets [39,40] and we shall take the advantage of the approach to anomalous subdiffusive transport developed recently in Refs. [22,23,39].

Dynamical approach to viscoelastic transport and stochastic modeling
We start with a well-established dynamical approach to the theory of Brownian motion [25,27]. The environment is modeled by a set of N 0 particles with masses m i harmonically coupled with the spring constants κ i to the Brownian particle of mass m: where q i are the coordinates of the environmental particles, x is the coordinate of the Brownian particle, and f (x, t) is an external force which affects only the motion of Brownian particle but not the environment. As usually, using the Green function of harmonic oscillator one can express the bath oscillators coordinates q i (t) in Eq. (2) via the initial values q i (0) and p i (0) for arbitrary x(t): Here, ω i = κ i /m i are the bath oscillators frequencies and the integration by parts has been used to obtain the second equality. By substituting (3) into (1) one immediately obtains where and Eq. (4) is still a purely dynamical equation of motion which is exact, the dynamics of irrelevant degrees of freedom is but excluded. For example, it describes a time-reversible dynamics, if external field f (x, t) does not violate the time-reversal symmetry ‡. Then the trajectories governed by Eq. (4) are also most obviously time-reversal symmetric (see [47] on this point from stochastic perspective). This is contrary to a widespread misperception. The statistical irreversibility comes about on the level of a bunch of trajectories averaged over different initial realizations of q i (0), p i (0). The loss of information due to averaging or coarse-graining is a well known source of irreversibility leading to statistical description of dynamical systems (see e.g. [49]). As a matter of fact, the dynamics described by Eq. (4) is just a projection of a highly dimensional Hamiltonian dynamics onto two dimensional (x, p) phase subspace. Furthermore, one proceeds as in a typical molecular dynamics setup. The initial positions q i (0) and momenta p i (0) of the environmental oscillators in Eq. (6) are sampled (this is the only non-dynamical element in the theory) from a canonical distribution at ‡ Time-reversal symmetry can be dynamically broken by an external time-dependent field. For example, a harmonic mixing driving, f (t) = A 1 cos(Ωt) + A 2 cos(2Ωt + φ) does violate the time-reversal symmetry dynamically for φ = 0 [48].
conditioned on the initial position of the Brownian particle x(0), where Z is the statistical sum of bath oscillators. Here one assumes that initial velocities of the bath oscillators are centered at zero, i.e. the medium is not moving as whole. Likewise, the average q i (0) = x(0), i.e. the medium is initially equilibrated adjusting to the Brownian particle localized at x(0). The corresponding random force ξ(t) becomes a stochastic process which is obviously Gaussian and can be completely characterized by its first two statistical moments. First moment is obviously zero, ξ(t) = 0. Furthermore, with Gaussian averages for any set of bath oscillators. This is the celebrated fluctuation-dissipation relation, or the second fluctuation dissipation theorem by Kubo [24]. The bath oscillators are conveniently characterized by the spectral density [27] It allows to express η(t) as η(t) = (2/π) ∞ 0 dωJ(ω) cos(ωt)/ω and the noise spectral density as S(ω) = 2k B T J(ω)/ω via the Wiener-Khinchin theorem, S(ω) = ∞ −∞ ξ(t)ξ(0) e iωt dt. The choice J(ω) = η α | sin(πα/2)|ω α with 0 < α < 2, yields the fractional Gaussian noise (fGn) ξ(t) introduced by Mandelbrot and van Ness [50]. For 0 < α < 1 (sub-Ohmic thermal bath [27]), where Γ(z) is a standard gamma function. This choice yields subdiffusion asymptotically, δx 2 (t) ∼ t α , for an ensemble of particles [27]. The corresponding GLE (4) is termed also the fractional GLE, or FLE [23,32,51,52] upon the use of the notion of fractional Caputo derivative to shorthand the frictional term, and the coefficient η α is termed then the fractional friction coefficient. Such a dynamical modeling requires a very large N 0 → ∞ number of the bath oscillators with a quasi-dense spectrum.

Stochastic modeling
Alternatively, one can use just a handful of N auxiliary Brownian particles clouding around the central particle [23,53] while modeling the rest as friction and noise acting on these representative ones. Such Brownian quasi-particles can serve to model sticky viscoelastic media. A particle clouded by other particles reminds conceptually polaron picture in condensed matter physics. Then, one replaces Eqs. (1,2) with where x i are the coordinates of auxiliary particles, k i are the corresponding coupling constants, √ 2η i k B T ζ i (t) are thermal Gaussian forces, with zero range of correlation, , and η i are the corresponding friction coefficients. Moreover, the overdamped limit for these auxiliary stochastic medium oscillators, m i → 0, yieldṡ where . The last equation for u i is similar to the Maxwell's relaxation equation for viscoelastic force in macroscopic theory of viscoelasticity [17] which is augmented by the corresponding Langevin force in accordance with FDR. Such a description was introduced in Refs. [22,23,39,53] to model anomalous Brownian motion in complex viscoelastic media within a generalized Maxwell model. For a particular case of the only one auxiliary particle (Maxwell model of viscoelasticity), the earlier description in Refs. [54,55] leading to GLE with thermal Ornstein-Uhlenbeck noise is readily reproduced. Excluding the dynamics of auxiliary variables u i (t) and assuming that initially forces u i (0) are thermally Gaussian-distributed with zero mean and dispersion u 2 i (0) = k i k B T yields again the GLE (4) with the corresponding memory kernel presented by a sum of exponentials, The corresponding noise ξ(t) presents the sum of Ornstein-Uhlenbeck noises. The power-law memory kernel (10) can be reliably approximated by such a sum [22] Physical meaning of ν 0 is a high-frequency (short-time) cutoff of the stochastic process ξ(t), t l = ν −1 0 . The low-frequency (long-time) cutoff corresponds to t h = b N −1 t l . Similar scaling and approximation to a power law are well known in the theory of anomalous relaxation [43,56]. By adjusting b and N one can approximate power law over about r = N log 10 b−2 time decades. Weak dependence of r on b and N ensures a very powerful and accurate numerical approach to FLE dynamics. Just a handful of auxiliary particles can suffice for all practical purposes. Of course, such a Markovian embedding of FLE is not unique [57]. However, it appeals by a clear physical interpretation. Recently, the set (13) was formally generalized to numerically integrate also superdiffusive FLE [58].

Flashing ratchet
To study flashing ratchet transport we take f ( is a deterministic force acting on the particle (Brownian motor) when a potential U(x) is switched on. Furthermore, r(t) switches on/off taking just two values, zero and one, r(t) = {1, 0}. For U(x) we take a typical ratchet potential with amplitude U 0 and period L. For r(t) two models are considered: periodic vs. random switching. In the periodic case, r(t) = 1/2 [1 + sign (sin(2πνt))], where sign(·) is standard signum function and ν is linear frequency of oscillations. The first switching occurs at t I = 1/2ν and the total number of switching is N sw ≡ T tot /t I = 2νT tot , where T tot is a total computation time which equals N s dt; N s is a total number of integration steps of Eq. (4), dt is a time-step. Finally, one can write N sw = 2νN s dt. For random switching r(t), we consider the model of dichotomous Markovian process (DMP) with the probability p = 2νdt to switch within dt. Such DMP is symmetric with equal residence time distributions in two states, ψ(τ ) = 2ν exp(−2ντ ). Here, 2ν is transition rate from one state to another one (and also inverse of the mean residence time in a state). The mean flashing frequency ν equals ν.
In all simulations we scale coordinate in units of L and time in units of τ v = (m/η α ) 1/(2−α) . It is assumed to be temperature-independent in accordance with the underlying Hamiltonian model [25,27]. This is a standard assumption done also in other ratchet models [1,2]. Furthermore, energy is scaled in E = L 2 η 2/(2−α) α /m α/(2−α) and temperature T = E/k B . In this work we choose α = 0.5 and b = 10, C 1/2 (10) = 1.3, and use ν 0 = 100, N = 12. Such a set of parameters gives an excellent approximation to the FLE dynamics over at least ten time decades, until t max = 10 8 . Similar to [22,39] this was checked by comparison with the exact analytical solution for the position variance obtained within GLE and FLE [24,51] in the force-free case. The numerical errors due to the memory kernel approximation are negligible as compare with the overall statistical error in our stochastic simulations (several percents), see for details in [22]. Like in [40], numerical solutions of stochastic differential equations were performed on the graphical processor units (GPUs) with double precision using stochastic Heun method. The use of GPU computing allowed to parallelize and accelerate simulations by a factor of about 100, as compare with conventional CPU computing.
Schematically, viscoelastic dynamics in the flashing ratchet potential is shown in figure 1 (see also video abstract). The auxiliary particles are placed on the level U = 0 and they are not influenced by potential. Initially, when the potential is switched on, the Brownian particle moves in the potential well (see snapshots a and b). When the potential is off (snapshots c and d) the Brownian particle "freely" subdiffuses, interacting only with auxiliary particles. After the potential is switched on again, Brownian particle moves into the next potential well (see snapshot e) with the net motion occurring in the negative direction. During such motion not all auxiliary particles follow immediately to the Brownian particle. Some of them are at large distances from the Brownian particle and are essentially less mobile than the Brownian particle and its nearest environment because they have essentially larger friction coefficient than other auxiliary particles. The weaker is the coupling constant k i of an auxiliary particle to the Brownian one the larger is the corresponding frictional coefficient η i . These very sluggish particles create a slowly fluctuating quasi-elastic force (a temporal biasing force [22,23]) acting on the central particle.

Results and discussion
In all simulations we monitored two main statistical quantities, the mean position x(t) and the position variance δx 2 (t) = x 2 (t) − x(t) 2 . For the ensemble average, N ens = 10 4 trajectories were used. The focus is on such physical quantities as subdiffusive current (subvelocity) v α , subdiffusion coefficient D α and generalized Peclet number [39] Pe α := v α /D α . The latter one surves as a measure for the coherence and quality of anomalous stochastic transport, by analogy with normal one [60]. The subvelocity and the subdiffusion coefficient are defined as usually [23,39,45,53] Here, the limit is understood in the following sense: t is large but still much smaller than the memory cutoff, which can always be made unreachable in numerical simulations and thus is irrelevant for the results presented. Practically, the values of v α and D α were calculated by fitting the numerical dependencies x(t) and δx 2 (t) with a power-law function a {v,D} t α , extracting the corresponding a v and a D within the last time window of simulations. Total time of simulations was varied in the interval [2 × 10 5 , 10 6 ] depending on system parameters to guarantee a good fit with α = 0.5 (convergency is slow); the time step was fixed at dt = 2 × 10 −3 . Typical dependencies for the mean particle position and variance on time for some different values of the potential amplitude U 0 and flashing frequency ν are shown in figure 2(a) and (b), respectively. Here black solid, red dash and blue dash-dot lines correspond to the periodic flashing, whereas green dash-dot-dot line relates to the random one. The direction of transport in figure 2(a) is negative and the transport has clearly subdiffusive character. We indicate in figure 2(a) also the corresponding power-law asymptotics. One can notice that stochastic flashing tends to delay the corresponding transport process in comparison with the periodic one of the same frequency (cf. red dash vs. green dashdot-dot lines), and the transport becomes faster with the increased flashing frequency (this is but not a universal feature, see below). The diffusion is initially always ballistic (universal regime), δx 2 (t) = v 2 T t 2 , where v T = k B T /m is thermal velocity. This is because the Brownian particle's velocity is initially thermally distributed and the action of the medium requires some time to settle in. After some transient, subdiffusion follows a) asymptotically to one and the same universal dependence, δx 2 (t) = 2D α t α /Γ(1 + α) with D α = k B T /η α , as in the absence of potential, independently of the potential height and the presence of driving, cf. in figure 2(b). This universality (or a weak sensitivity in the case of a strong fast driving [39]) is a benchmark of viscoelastic subdiffusion [22,23,39,53]. We elaborate on this fact in more detail below.

Ergodicity
An important issue in anomalous transport is ergodicity [61,62], i.e. whether timeaverage over a single particle trajectory delivers non-random result, same for all identical particles, or a principal randomness and scatter in the single-particle averages emerge. If this is the case, then the behavior and fate of each individual though identical particle is different, even in the limit of infinitely long trajectories, and only ensemble-averaging smears out this principal randomness and delivers non-random result. In the case of ergodic transport, a single-trajectory average should coincide with the result of the ensemble averaging. For example, CTRW subdiffusion with divergent mean residence times in traps is patently nonergodic [63] and single-trajectory averages obey some universal fluctuation laws [63][64][65][66] leading to a universal scaling relating such transport and diffusion in periodic potentials [68]. On the contrary, viscoelastic subdiffusion, both free [67], and in time-independent potentials [22] was shown to be mostly ergodic, though some transient nonergodic features can be present. Whether it remains ergodic also in time-fluctuating fields presents a nontrivial problem. First, we checked the asymptotic ergodicity of ratchet viscoelastic transport which is expected. To show that this indeed is the case we plotted the fluctuating subvelocity v α (t), obtained for a single trajectory as v α (t) = Γ(1 + α)x(t)/t α in figure 3(a)   averages coincide and the transport is ergodic. For any finite t, statistical fluctuations are present, of course, also in ergodic case. Next, a single trajectory average of the squared displacement can be obtained as [22,64,65,67], In the ergodic case, this average should coincide in the limit T → ∞ but for any t with the result of the ensemble averaging. If the agreement holds only for large t, then diffusion is asymptotically ergodic. For viscoelastic subdiffusion in fixed periodic potentials the agreement holds both for small t (ballistic regime) and for large t [22]. However, some deviations from ergodicity occur for intermediate t, on the time scale of diffusion over several potential periods. This is because even if the mean escape time exists the escape kinetics remains anomalous for activation barriers of an intermediate height [22]. Diffusional spreading derived from single trajectories is compared with the corresponding ensemble averages in figure 3(b). In the ballistic regime, using Clearly, if in the limit T → ∞ the time average of v 2 (t) coincides with the ensemble average yielding v 2 T , then ergodicity holds also in the ballistic regime. In the absence of driving this is indeed the case [22]. However, a periodic driving exciting coherent oscillations emerging due to a combined action of the external trapping potential and the viscoelastic cage effect can cause large additional temporal fluctuations in v(t) which are not self-averaged in v 2 (t) T . Accordingly, such a periodically driven subdiffusion ceases to be ergodic in the ballistic regime. Temporally but on a large intermediate time scale ergodicity is generally broken. Nevertheless, asymptotically it remains ergodic as figure 3(b) implies. Convergency to such asymptotically ergodic regime is very slow. However, it improves with the increased frequency of flashing, see the case ν = 5 in figure 3(b). For random driving, ergodicity in the discussed sense holds also in ballistic regime.

Frequency dependence
We proceed further with the dependence of anomalous current v α on frequency ν at a fixed value of temperature for different potential amplitudes U 0 . In the case of fluctuating tilt subdiffusive ratchets [39] this dependence was especially intriguing. That anomalous transport vanished in the limit of vanishing modulation frequency (adiabatic driving limit) and exhibited a maximum at intermediate frequencies.
The origin of this maximum was related to stochastic resonance effect in Ref. [40]. This is in a sharp contrast with normal diffusion fluctuating tilt ratchets where transport optimizes namely in the adiabatic limit. The results for periodic flashing are shown in figure  4 and for stochastic driving in figure 5. For small driving frequencies subdiffusive transport occurs always in the natural negative direction, which is opposite to the natural positive direction of fluctuating tilt ratchet. Our anomalous Brownian motors share these features with their normal diffusion counterpart [1,2]. There are but novel features which are due to a profound role of the inertial effects in our model. We recall that in the used scaling the unit of time corresponds to a scaling time constant τ v of the decay of velocity autocorrelation function. Correspondingly, the ballistic regime occurs until t ∼ 1, when the potential is off. For periodic driving and a small potential height U 0 = 0.25 = T (see red line with squares in figure 4) subvelocity takes negative values for all flashing frequencies considered. The dependence of |v α | on ν has a unique maximum at ν c0 ≃ 0.5, where subdiffusive transport is optimized, with negative v α . For the fast-flashing regime, when ν > 1 subvelocity asymptotically decays to zero with ν. Increase in the ratchet potential height 2U 0 accompanied by an increased role of the inertial effects leads to a profound change in the character of the frequency-dependence v α (ν) (see black curve with circles for U 0 = 0.5 in figure 4). Let us study this dependence in more detail. The first difference from the former case is a multiple (double) current inversion, realized around ν c1 ≃ 0.335 and ν c2 ≃ 0.402. Enlarged picture of this inversion effect is shown in the insert a) placed on the right panel in figure 4. Here, a well defined peak occurs at ν max 1 = 0.375 with positive values for ν c1 < ν < ν c2 . Furthermore, on the left side of ν < ν c1 there is a maximum of |v α | for negative v α . Increase of ν beyond ν c2 leads to a global maximum of |v α | realized at ν max 2 ≃ 0.75. With a further increase in flashing frequency a complicated oscillatory pattern with several more relative minima and maxima in |v α | emerges. With a further increase of the potential amplitude these oscillatory resonance-like features become more pronounced -compare two cases with U 0 = 0.50 and U 0 = 0.75 in figure 4. This is clearly related to an increased role of the inertial effects. First, the amplitude of oscillations in v α (ν) increases. Second, ever more extrema appear, see insets b) and c) in figure 4. The situation is quite different in the case of random flashing. Dependencies of anomalous current versus mean linear frequency ν at fixed temperature T = 0.25 for different values of U 0 are shown in figure 5. Here, contrary to the case of periodic flashing one has more simple dependencies v α ( ν ) with one broad minimum for each considered potential amplitude U 0 . Increase in U 0 leads to increase in absolute values for subvelocity, broadening the minimum and shifting minimum position towards larger values of mean linear frequency. Thus, in the random flashing case the system is not manifestly sensitive to inertial effects due to stochastic nature of flashing. The maximal value of subvelocity becomes also much smaller, compare figure 4 and figure 5, due to the absence of resonances.
A benchmark feature of viscoelastic subdiffusion is that it asymptotically does not depend on the presence of static periodic potential [22,23,53]  weakly dependent on driving in the case of rocking subdiffusive ratchets [39]. The diffusional spread in flashing ratchets is also practically not dependent on the potential amplitude U 0 and flashing frequency as figure 6 illustrates. One can see, that for both periodic and random flashing the anomalous diffusion coefficient does not display any profound dependence on the flashing frequency and the potential amplitude. It takes values around 0.25, which is D α = T in the scaled units, with a deviation which is less than 3% and lies within the statistical error margin of our simulations. Therefore, the corresponding generalized Peclet number Pe α := v α /D α , which measures the coherence and quality of transport, resembles the behavior of the absolute value of subcurrent in figures 4 and 5, which should be multiplied by factor 4 in order to obtain Pe α .

Temperature dependence
The temperature dependence of subvelocity v α , and generalized Peclet number Pe α is mostly intriguing, whereas the dependence of subdiffusion coefficient D α on temperature is expected to be trivial, D α (T ) = T . We plotted numerical results for all these three quantities in figure 7 8(a). First, the subvelocity is finite in the limit T → 0 (this is always so for random flashing, but only for certain frequency windows in the case of periodic flashing). Second, it diminishes to zero in the limit of high temperatures. The latter is expected, but the former is not, being highly surprising. Indeed, in a combination with the subdiffusion coefficient following to the expected linear dependence and vanishing in the limit T → 0 [see figures 7(b) and 8(b)], the non-zero value of v α in this limit means that the basic mechanism of operating our subdiffusive viscoelastic Brownian motors is very different from one in the case of a) normal Brownian motion. This is at odds with our expectations (see Introduction). Indeed, in the absence of viscoelastic effects the flashing Brownian motors of the studied kind are known to be essentially based on thermal diffusion, i.e. such transport is assisted by thermal fluctuations [9]. It is therefore expected to vanish in the limit T → 0, where the thermal fluctuations vanish (in neglecting quantum effects). Our results say, however, that viscoelastic flashing ratchets are not based primarily on the thermal diffusion and thermal fluctuations. These are long-range elastic correlations which are at work here in combination with the potential fluctuations. Namely, the external field acting on the Brownian particle changes abruptly (on-off process) the energy stored in the elastic springs. Not only the Brownian particle moves under the influence of f 0 (x), but also the auxiliary particles under the influence of Brownian particle. Some of them rapidly adjust their positions, but there are also slow particles which cannot follow immediately. When the potential is off these extra sluggish particles pull the Brownian particle. This is why its motion is generally not frozen even at T = 0 within our purely classical setup. Because subdiffusion diminishes with lowering temperature, There are also profound differences in the temperature-dependence of anomalous transport in the cases of periodic and random flashing. So, in the case of periodic flashing the transport can occur in the counterintuitive positive direction, see the curve for U 0 = 0.5, ν = 0.375 in figure 7(a). In this case, v α remains finite also at T = 0, and an inversion of the current direction with increasing temperature is possible. This is not so for other three curves in figure 7(a), where transport vanishes jump-like at T = 0. Therefore, it seems obvious that for a periodic flashing there are frequency windows for permitted transport at T = 0 (a detailed investigation of this issue is left for a separate study). On the contrary, for random driving in figure 8, a, the transport remained always finite at T = 0. Moreover, here transport always occurs in the intuitive negative direction. Both in the random case and in the periodic case the transport can be optimized at some temperature T max = 0. It can be but also that T max = 0.

Conclusions
With this work we put forward anomalously slow Brownian motors based on flushing subdiffusion in viscoelastic media. Both anomalous transport and subdiffusion were shown to be asymptotically ergodic. However, such driven subdiffusion can be transiently nonergodic on appreciably long time scales. The transport exhibits a remarkably good quality at low temperatures. It remains finite even at zero temperature in many cases (but not always in the case of periodic flashing) with zero dispersion. In this limit, the transport is clearly induced by long-range elastic correlations in a combination with potential flashing and not by the thermal noise of environment. This is surprising and at odds with a popular explanation of the origin of noise-induced flashing transport in the case memory-free stochastic dynamics [1,2]. Moreover, in the case of periodic driving anomalous current exhibits multiple resonance-like features, can flow in the counter-intuitive direction, and invert its direction both with the change of flashing frequency and with change of temperature. The inertial cage effects are essential for many observed features. In particular, they are crucial for the resonance-like features manifested for high potential barriers in the case of periodic flashing. We hope that our research will stimulate further cross-fertilization of the field of fluctuation-induced transport and the field of anomalously slow transport and subdiffusion which flourish at present with little interaction. There emerges an increasing experimental support for the occurrence of subdiffusion in cytoplasm of biological cells and we believe that there might be also place for operating sluggish molecular motors of the kind described.

Acknowledgments
Support of this research by the Deutsche Forschungsgemeinschaft, Grant GO 2052/1-1 is gratefully acknowledged.