Non-monotonic temperature dependence of chaos-assisted diffusion in driven periodic systems

The general idea of controlling transport is associated with the concept of implementing actions to ensure that a system exhibits desired characteristics. For example, Nature is well known in presenting sophisticated mechanisms that regulate phenomena which take place in all scales of time and space [1]. A careful and systematic exploration of these mechanism reveals that sometimes it is more reliable to control the setup by letting it to fluctuate and change its dynamics as little as possible to drive it to the desired state without applying intense external forces. A typical scenario Nature uses is the resource to tailored non-equilibrium forces, as it is the case for molecular and Brownian motors [2–6]. Brownian machinery presents an archetypal scheme in which, even in the absence of an externally applied bias, directed motion emerges from unbiased environmental non-equilibrium perturbations via the mechanism of breaking the spatial symmetry of the setup [3]. The origin of directed transport stems from the fact that the combined action of non-equilibrium perturbations of stochastic or deterministic nature takes the system out of thermal equilibrium and breaks the detailed balance symmetry. This working principle can be seen as a key for understanding various processes occurring in numerous disciplines of science, including e.g. intracellular transport [7]. In view of the widespread applications of Brownian motors, directed motion controllability has become a focal point of research in non-equilibrium statistical physics which inspired a plethora of new microscale devices displaying unusual transport features [8–16].

Irregular behaviour generally implies also unpredictability, which, when it comes to applications one typically attempts to avoid or at least to minimise. A particular aspect of unpredictability is the appearance of diffusive spreading constituting yet another salient transport quantifier. Beyond doubt, the problem of diffusion has played a central role in the development of both the foundations of thermodynamics and statistical physics [17] and has led to numerous applications within a diversity of scientific disciplines. Inspired by the physics of Brownian motors, we ask to what extent it is possible to influence the diffusive behaviour of the system in presence of thermal noise and deterministic external perturbations and elucidate the mechanisms behind such characteristics. The conventional diffusive spreading characterised by the diffusion coefficient D increases with the temperature of the medium surrounding the system [18–20]. Here, we find that the diffusion coefficient D assumes a non-monotonic behaviour as a function of temperature, possessing a bell-shaped maximum at some intermediate temperature. This form of peculiar diffusive behaviour is markedly distinct from the known case of anomalous diffusion for which the mean square displacement of the particle $\langle {\rm{\Delta }}{x}^{2}(t)\rangle$ grows asymptotically according to a power law $\langle {\rm{\Delta }}{x}^{2}(t)\rangle \sim {t}^{\alpha }$ with $\alpha \ne 1$ [21, 22]. The studied non-monotonic diffusion here must also be distinguished from the phenomenon of giant diffusion. The latter has been observed in the diffusive motion of a Brownian particle in a tilted periodic landscape [23–27], in disordered potentials [28, 29] or for an adiabatically driven inertial Brownian particle moving in a periodic one-dimensional geometry [30].

To investigate our objective of a non-monotonic dependence of diffusion on temperature we consider the non-equilibrium dynamics of a massive Brownian particle moving in a periodic potential under the influence of a time-periodic force. The class of models described by such a dynamics includes periodically driven pendulums [31], super-ionic conductors [32], Josephson junctions [33], dipoles rotating in external fields [34], phase-locked loops [35], dislocations in solid state physics [36], solitons described by the Sine-Gordon equation [37], the Frenkel-Kontorova lattices [38], dynamics of adatoms subjected to a time-periodic force [39], charge density waves [40] and cold atoms in optical lattices [41], to name but a few. The rich physics contained in this model has become evident over the last decades with numerous studies. In particular, with present periodic driving this class of systems comprises operational regimes that are deterministically chaotic. The sensitive dependence on initial conditions and the abundance of unstable periodic attractors are the most salient characteristics of chaotic behaviour [42]. The combination of these features, possibly assisted in addition with noise driving, make chaotic systems one of the most flexible setups enabling the emergence of the discussed peculiar diffusive behaviour.

We formulate the problem in terms of a classical particle of mass M which is (i) subjected to a spatially periodic potential $U(x)=U(x+L)$ of period L possessing the mirror symmetry, i.e. there exists x₀ such that the relation $U({x}_{0}-x)=U({x}_{0}+x)$ holds for all x, (ii) additionally being agitated by an external unbiased time-periodic deterministic force $A\cos ({\rm{\Omega }}t)$ with angular frequency Ω and of amplitude strength A, and (iii) coupled to a thermal bath at temperature θ. The overall dynamics of such a Brownian particle is determined by the following Langevin equation

$$\begin{eqnarray}&&M\ddot{x}+{\rm{\Gamma }}\dot{x}=-U^{\prime} (x)+A\cos ({\rm{\Omega }}t)+\sqrt{2{\rm{\Gamma }}{k}_{B}\theta }\,\xi (t).\end{eqnarray} \tag{ 1 }$$

Here, the dot and the prime denote differentiation with respect to the time t and the particle coordinate x, respectively. Thermal fluctuations due to the coupling of the particle with the thermal bath are modelled by δ-correlated Gaussian white noise $\xi (t)$ of zero mean and unit intensity, i.e.,

$$\begin{eqnarray}&&\langle \xi (t)\rangle =0,\quad \langle \xi (t)\xi (s)\rangle =\delta (t-s).\end{eqnarray} \tag{ 2 }$$

The parameter Γ denotes the friction coefficient and k_B is the Boltzmann constant. The noise intensity factor $2{\rm{\Gamma }}{k}_{B}\theta$ with temperature θ of the heat bath follows from the fluctuation-dissipation theorem [43]. The latter ensures the stationarity of the thermal canonical Gibbs state for vanishing driving under the dynamics governed by equation (1), i.e., when A = 0 and in presence of periodic boundary conditions. We focus on the archetypical nonlinear situation, namely, a sinusoidal potential with a barrier height $2{\rm{\Delta }}U$,

$$\begin{eqnarray}&&U(x)={\rm{\Delta }}U\sin \left(\displaystyle \frac{2\pi x}{L}\right).\end{eqnarray} \tag{ 3 }$$

Upon introducing the period L and the parameter combination ${\tau }_{0}={\rm{\Gamma }}{L}^{2}/{\rm{\Delta }}U$ as units for length and time, respectively, equation (1) is written in a dimensionless form [44], reading,

$$\begin{eqnarray}&&m\ddot{\hat{x}}+\dot{\hat{x}}=-\hat{U}^{\prime} (\hat{x})+a\cos (\omega \hat{t})+\sqrt{2Q}\,\hat{\xi }(\hat{t}),\end{eqnarray} \tag{ 4 }$$

where $\hat{x}=x/L$ and $\hat{t}=t/{\tau }_{0}$. The dimensionless potential $\hat{U}(\hat{x})=U(x)/{\rm{\Delta }}U=U(L\hat{x})/{\rm{\Delta }}U=\hat{U}(\hat{x}+1)$ assumes the unit period L = 1 and the potential amplitude ${\rm{\Delta }}\hat{U}=1$. The remaining re-scaled parameters are the mass $m=(1/{\rm{\Gamma }}{\tau }_{0})M$, the amplitude $a=(L/{\rm{\Delta }}U)A$ and the angular frequency $\omega ={\tau }_{0}{\rm{\Omega }}$. The dimensionless thermal noise reads $\hat{\xi }(\hat{t})=(L/{\rm{\Delta }}U)\xi (t)=(L/{\rm{\Delta }}U)\xi ({\tau }_{0}\hat{t})$ and assumes the same statistical properties, namely it is Gaussian with $\langle \hat{\xi }(\hat{t})\rangle =0$ and $\langle \hat{\xi }(\hat{t})\hat{\xi }(\hat{s})\rangle =\delta (\hat{t}-\hat{s})$. The dimensionless noise intensity $Q={k}_{B}\theta /{\rm{\Delta }}U$ is given by the ratio of thermal energy ${k}_{B}\theta$ and half of the activation energy the particle needs to overcome the original potential barrier $2{\rm{\Delta }}U$. From now on, we will use only the dimensionless variables and therefore shall skip the 'hat' for all quantities appearing in equation (4).

The observable of foremost interest in this study is the diffusion coefficient D which characterises the spread of trajectories and fluctuations around the average position of the particle, namely [45, 46],

$$\begin{eqnarray}&&D=\mathop{\mathrm{lim}}\limits_{t\to \infty }\displaystyle \frac{\langle {[x(t)-\langle x(t)\rangle ]}^{2}\rangle }{2t},\end{eqnarray} \tag{ 5 }$$

where the averaging $\langle \cdot \rangle$ is over all realisations of thermal fluctuations as well as over initial conditions for the position x(0) and the velocity $\dot{x}(0)$. The latter is necessary because in the deterministic limit of vanishing thermal noise intensity $Q\to 0$ the dynamics may possess several coexisting attractors thus being non-ergodic and implying that the corresponding results may be affected by a specific choice of those selected initial conditions [46, 47].

In order to explain the discussed peculiar behaviour of the diffusion coefficient we will also investigate velocity of the Brownian particle. Due to the presence of the external time-periodic driving $a\cos (\omega t)$ as well as the friction term $\dot{x}$ in equation (4) the particle velocity $\dot{x}(t)$ approaches a unique asymptotic state in which it is characterised by a temporally periodic probability density. This latter density function has the same period T as the driving [48, 49]. In particular, the first statistical moment of the instantaneous particle velocity $\langle \dot{x}(t)\rangle$ assumes for a long time the form

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{t\to \infty }\langle \dot{x}(t)\rangle =\langle {\bf{v}}\rangle +v(t),\end{eqnarray} \tag{ 6 }$$

where $\langle {\bf{v}}\rangle$ is the directed transport velocity. Note that all three terms entering the right hand side of equation (4) are unbiased: the average of the potential force over a spatial period L vanishes as well as that of the time-dependent driving over a temporal period T, and also the average of the random force $\xi (t)$ vanishes according to equation (2). Moreover, all three elements are symmetric and in consequence there is no directed transport in the long time stationary regime, i.e. $\langle {\bf{v}}\rangle =0$. The deviation v(t) is periodic with period T and has a vanishing time average, i.e. $v(t+T)=v(t)$ and $(1/T){\int }_{0}^{T}{\rm{d}}{t}\,v(t)=0$. Therefore it is useful to consider also the period averaged velocity ${\bf{v}}(t)$ as

$$\begin{eqnarray}&&{\bf{v}}(t)=\displaystyle \frac{1}{T}{\int }_{t}^{t+T}{\rm{d}}{s}\,\dot{x}(s),\end{eqnarray} \tag{ 7 }$$

where $T=2\pi /\omega$ denotes a period of the external driving $a\cos (\omega t)$.

The particle positions $x({nT})$ at integer multiples n of the period T can be exactly expressed by the sequence of period averaged velocities ${\bf{v}}({kT})=(1/T){\int }_{{kT}}^{(k+1)T}{\rm{d}}{s}\,\dot{x}(s)$ defined in equation (7) as

$$\begin{eqnarray}&&x({nT})=T\displaystyle \sum _{k=0}^{n-1}{\bf{v}}({kT}),\end{eqnarray} \tag{ 8 }$$

where we assumed that the particle starts at $x(0)=0$ (for $x(0)\ne 0$ one can re-define $x({nT})\to \delta x({nT})=x({nT})-x(0)$). The period averaged velocities at multiple integers of the period T have a vanishing mean value. Therefore the average position remains zero and the variance of the position is given by its second moment reading

$$\begin{eqnarray}&&\langle {x}^{2}({nT})\rangle ={T}^{2}\displaystyle \sum _{k,l}^{n-1}\langle {\bf{v}}({kT}){\bf{v}}({lT})\rangle .\end{eqnarray} \tag{ 9 }$$

In the asymptotic limit of large times and with the temporally periodic probability density the velocity auto-correlation function depends only on the time difference. In view of equation (9) the left hand side approaches $\langle {x}^{2}({nT})\rangle =2{nTD}$ and hence the diffusion coefficient D can be expressed in terms of the temporal period averaged velocity autocorrelation function, yielding

$$\begin{eqnarray}&&D=T\displaystyle \sum _{k=0}^{\infty }\langle {\bf{v}}({kT}){\bf{v}}(0)\rangle .\end{eqnarray} \tag{ 10 }$$

An equivalent expression is given in [50]. The infinite time limit in equation (5) implies that the auto-correlation function has to be determined for the stationary ensemble of the period averaged velocities.

Because neither the Langevin equation (4) nor the corresponding Fokker-Planck variant can be solved within analytical means we performed comprehensive numerical simulations of the model. We did so by employing a weak version of the stochastic second-order predictor-corrector algorithm with a time step typically set to about $({10}^{-3}-{10}^{-2})\times T$ [51]. Because equation (4) is a second-order differential equation, we need to specify two initial conditions, namely x(0) and $\dot{x}(0)$. We choose x(0) and $\dot{x}(0)$ to be equally distributed over the intervals $[0,1]$ and $[-2,2]$, respectively. Our quantities of interest were averaged over ${10}^{3}\mbox{--}{10}^{5}$ sample trajectories. All numerical calculations have been performed by use of a CUDA environment as implemented on a modern desktop GPU. This procedure did allow for a speedup of a factor of the order 10³ times as compared to a common present-day CPU method [51].

The dynamical system described by equation (4) exhibits an extremely rich behaviour as a function of the four dimensionless parameters $\{m,a,\omega ,Q\}$ with dimensionless friction at ${\rm{\Gamma }}=1$. A general overview is provided, e.g. with [33]. Our objective is not the systematic exploration of the setup at hand but rather the search for the above mentioned non-monotonic behaviour of the diffusion constant with increasing noise intensity as shown in figure 1. The corresponding parameters of the dynamical system (4) are chosen as $\{m=0.9,a=8.7,\omega =0.275\}$. At each depicted noise strength $Q\propto \theta$ the diffusion constant D was estimated from 10⁵ trajectories by means of equation (5). At low noise intensity D increases with Q until it reaches a local maximum at $Q\approx 2\times {10}^{-5}$. From there it decreases to a minimum at $Q\approx 5\times {10}^{-3}$ turning over to a monotonic function of Q, and, finally, at sufficiently large values of Q becoming strictly proportional to Q, i.e. to the temperature θ of the ambient thermal bath. This high temperature behaviour, however, is not depicted in figure 1. The decrease of the diffusion constant with increasing temperature $Q\propto \theta$ is counter-intuitive. It stays in clear contrast with the Einstein relation $D\propto \theta$ as well as with other known formulas, e.g. Vogel-Fulcher-like laws [52] or Arrhenius-type behaviour for the diffusion of a Brownian particle in periodic potentials [53–55]. Because relevant structural elements of the stochastic dynamics described by equation (4) are determined by its deterministic properties as the first step of our analysis we consider the noiseless case.

Zoom In Zoom Out Reset image size

3.1. Noiseless dynamics: Q = 0

For the set of parameters presented in figure 1, the deterministic system (when Q = 0) exhibits chaotic behaviour with a dense set of unstable periodic orbits. For $Q\gt 0$, the system 'feels' some unstable orbits which play an important role in controlling diffusion properties. Below we explain this. For this aim, in figure 2 we present bifurcation diagrams for the positions $x({kT})$ (panel (a)) and the period averaged velocities ${\bf{v}}({kT})$ (panel (b)) as a function of the driving amplitude a of the periodic driving $a\cos (\omega t)$ for the deterministic system (4) with Q = 0. For any given value of the driving amplitude a the red dots in panel (a) represent the stroboscopic positions of the particle $x({kT})\mathrm{mod}\,L$, i.e. the projected values onto the spatial principal period $[0,1]$, for 1024 trajectories with $k={10}^{4}$. The velocities displayed in panel (b) show the according values ${\bf{v}}({kT})$ obtained from equation (7). The bifurcation diagrams reveal three qualitatively different dynamical regimes. For the smallest displayed values of $a\lt {a}_{0}=8.6867$ as well as for $a\gt {a}_{c}=8.6893$ the motion is chaotic. The stroboscopic positions projected onto the principal spatial period almost densely cover the available interval. In both regions the transport is diffusive [56], i.e. the diffusion constant determined by equation (5) has a finite value larger than zero. In the window between a₀ and a_c the motion is phase locked with the corresponding winding number $w=\pm 9$: the particle proceeds within one temporal period by 9 spatial periods either in positive or in negative direction depending on the initial condition. This is reflected by the two velocities ${\bf{v}}({kT})\approx \pm 0.4$ visible in panel (b) in the corresponding parameter window. In panel (a) two period doubling cascades start at the parameter value a₁ = 8.6878 and terminate at ${a}_{\infty }=8.6886$. In the small parameter interval between ${a}_{\infty }$ and a_c the asymptotic motion is phase locked taking place on two coexisting chaotic attractors. The window of phase locked motion ends at a_c with an attractor merging crisis [57] giving rise to a diffusive attractor covering the full principal period. Also a wide range of velocities suddenly emerges at a_c and continues to be present at larger values of the parameter a. However, one still detects a strong concentration of velocities at several values of ${\bf{v}}({kT})$ corresponding to integer winding numbers characterising unstable periodic orbits. Particularly large is the probability to find temporal period averaged velocities ${\bf{v}}({kT})\approx \pm 0.4$, corresponding to the winding numbers $w=\pm 9$. These two, together with the locked trajectories with w = 0 and a few other phase locked trajectories seem to constitute the backbone of unstable periodic orbits supporting the chaotic motion. As a simplified picture of the chaotic dynamics one may think of a process in which these unstable orbits are visited in a random sequence.

Zoom In Zoom Out Reset image size

3.2. Influence of thermal noise: Q > 0

In figure 3 we present the influence of thermal noise on velocity of the Brownian particle. We depict time series of period averaged velocities ${\bf{v}}({kT})$ resulting from a single trajectory in the large time limit for the same set of parameters specified in figure 1. Panel (a) exemplifies the deterministic case, with Q = 0. Regions close to the dominant winding numbers $w=\pm 9$ and w = 0 are significantly more frequently visited than others. Yet, the trajectory almost never dwells near any of these states for a longer period of time but typically leaves the state within the period after it had arrived there. Panels (b) and (c) illustrate the influence of noise on the period averaged velocity time series for the two values $Q=2.16\times {10}^{-5}$ and $Q=5.25\times {10}^{-3}$ which correspond to the positions of the maximum and the minimum of the diffusion constant D displayed in figure 1, respectively. At the smaller noise strength displayed in panel (b) still some of the fine details of the deterministic time series are visible while others are washed out. At the larger noise intensity (panel (c)) all fine details have disappeared apart from the fact that positive and negative velocities near ${\bf{v}}({kT})\approx \pm 0.4$ still occur with relatively high probability as well as small velocities ${\bf{v}}({kT})\approx 0$.

Zoom In Zoom Out Reset image size

By introducing the thresholds at ${\bf{v}}({kT})=\pm 0.2$ and counting the number of period averaged velocities which are between and outside these values one may estimate the probabilities ${p}_{R}=\mathrm{prob}(| {\bf{v}}({kT})| \gt 0.2)$ and ${p}_{L}=\mathrm{prob}(| {\bf{v}}({kT})| \leqslant 0.2)$ for the occurrence of the running and locked states, respectively. These probabilities are displayed in figure 4 as a function of the noise intensity Q. We note that the running states occur significantly more frequent at weaker noise than at larger noise. Considering as a very rough model an independent sequence $\{{v}_{k}\}$ of velocities ${\bf{v}}({kT})\approx \pm 0.4,0$ occurring with the respective probabilities ${p}_{R},{p}_{L}$ one finds a diffusive behaviour for the spreading of the positions $x({nT})={\sum }_{k=0}^{n-1}{\bf{v}}({kT})$ with the diffusion constant $D=0.16{{Tp}}_{R}$. With the observed noise-dependence of the probability p_R one already finds a qualitative agreement with the non-monotonic behaviour of the diffusion constant displayed in figure 1. However, in the next section we will present a slightly more realistic model based on the same three velocity states but with a more adequate dynamics.

Zoom In Zoom Out Reset image size

In order to describe the non-monotonic behaviour of the diffusion constant D in dependence of the noise strength Q we consider a simplified model only, containing transitions between the most relevant unstable orbits. As such we consider the phase locked trajectories with winding numbers $w=\pm 9$ corresponding to the velocities ${v}_{\pm }\approx \pm 0.4$ and those with winding number w = 0 corresponding to ${v}_{0}=0$. The model is assumed to be symmetric with respect to the two running states ${v}_{\pm }$. The sequence of period averaged velocities ${\bf{v}}({nT})$ is then modelled by a Markov process in terms of transitions between these states. Because of the assumed symmetry the matrix ${\bf{M}}$ describing these transitions has the form

$$\begin{eqnarray}{\bf{M}}:= \left(\begin{array}{ccc}q & \tfrac{1}{2}(1-k) & r\\ 1-q-r & k & 1-q-r\\ r & \tfrac{1}{2}(1-k) & q\end{array}\right),\end{eqnarray} \tag{ 11 }$$

where q denotes the conditional probability to remain staying in the running state ${v}_{+}$, ${v}_{+}\to {v}_{+}$, and likewise in ${v}_{-}$, ${v}_{-}\to {v}_{-};$ further k, the conditional probability to remain staying in the resting state 0, $0\to 0;$ and r, the conditional probability of a transition between opposite running states ${v}_{+}\to {v}_{-}$ and ${v}_{-}\to {v}_{+}$. Due to the assumed symmetry of the running states the transitions from the locked into a running state ($0\to {v}_{\pm }$) is $(1-k)/2$ and from a running into the locked state (${v}_{\pm }\to 0$) is $1-q-r$. Probabilities q and k may take non-negative values less than or equal to 1, while r is restricted to $0\leqslant r\leqslant 1-q$. The transition probabilities can be estimated from simulations of the Langevin equation (4) as the relative frequencies with which the period averaged velocity $v={\bf{v}}({nT})$ visits the three coarse grained regions ${V}_{+}=\{v| v\geqslant 0.2\}$, ${V}_{0}=\{v| \,-\,0.2\lt | v| \,\lt \,0.2\}$ and ${V}_{-}=\{v| v\,\leqslant \,-0.2\}$. The three independent rates k, q, r are displayed in figure 5. Finally, in the present model the states are 'decorated' with the respective velocities $\{{v}_{\pm },0\}$ such that the particle changes its position within a period by ${vT}$ with $v\in \{{v}_{\pm },0\}$.

Zoom In Zoom Out Reset image size

Collecting the probabilities of finding the particle at the time t = nT in either of the three states as a vector, ${\bf{p}}(n)=({p}_{+}(n),{p}_{0}(n),{p}_{-}(n))$ one can write its dynamics as

$$\begin{eqnarray}&&{\bf{p}}(n+1)={\bf{M}}p(n)\end{eqnarray} \tag{ 12 }$$

having the formal solution

$$\begin{eqnarray}&&{\bf{p}}(n)={{\bf{M}}}^{n}{\bf{p}}(0)\end{eqnarray} \tag{ 13 }$$

with the initial probability ${\bf{p}}(0)$. The stationary distribution ${{\bf{p}}}^{\mathrm{st}}$ is invariant under transitions, and hence is the solution of

$$\begin{eqnarray}&&{{\bf{p}}}^{\mathrm{st}}={{\bf{Mp}}}^{\mathrm{st}}\end{eqnarray} \tag{ 14 }$$

given by

$$\begin{eqnarray}&&{{\bf{p}}}^{\mathrm{st}}=\left(\begin{array}{c}\tfrac{1-k}{2(2-k-q-r)}\\ \tfrac{1-q-r}{2-k-q-r}\\ \tfrac{1-k}{2(2-k-q-r)}\end{array}\right).\end{eqnarray} \tag{ 15 }$$

In figure 6 the stationary probability ${p}_{+}^{\mathrm{st}}$ to find the particle in the running state is compared with ${p}_{R}/2$ estimated from simulations of the Langevin equation (4). The diffusion constant can be determined in the framework of the Markovian three-state model by means of equation (10) which is based on the stationary auto-correlation function of period averaged velocities. This correlation function follows from the three-state model as

$$\begin{eqnarray}&&\langle {\bf{v}}({nT}){\bf{v}}(0)\rangle =\displaystyle \sum _{\alpha ,\beta }{v}_{\alpha }{v}_{\beta }{({{\bf{M}}}^{n})}_{\alpha ,\beta }{p}_{\beta }^{\mathrm{st}},\end{eqnarray} \tag{ 16 }$$

where α and β label the three velocity states and the matrix element ${({{\bf{M}}}^{n})}_{\alpha ,\beta }$ determines the conditional probability to find the state labelled by α after n periods provided the system has started in the state β. Accordingly, ${p}_{\beta }^{\mathrm{st}}$ denotes the probability (15) to find β in the stationary state. In order to evaluate the nth power of the matrix ${\bf{M}}$ it is convenient to use its spectral representation reading

$$\begin{eqnarray}&&{\bf{M}}=\displaystyle \sum _{\kappa =1}^{3}{\lambda }_{\kappa }{{\rm{\Pi }}}^{(\kappa )},\end{eqnarray} \tag{ 17 }$$

where ${\lambda }_{\kappa }$ are the eigenvalues of ${\bf{M}}$ and ${{\rm{\Pi }}}^{(\kappa )}$ are matrices projecting onto the respective eigenvectors. The eigenvalue ${\lambda }_{1}=1$ results as a consequence of the conservation of total probability. The corresponding matrix ${{\rm{\Pi }}}^{1}$ projects any vector ${\bf{q}}$ onto the stationary probability ${{\bf{p}}}^{\mathrm{st}}$, i.e. it acts as ${{\rm{\Pi }}}^{1}{\bf{q}}={{\bf{p}}}^{\mathrm{st}}{\sum }_{\alpha }{q}_{\alpha }$. Explicit expressions of the other eigenvalues and eigenprojectors are given in the appendix. Inserting the spectral representation (17) into the expression (16) one finds that the first eigenvalue ${\lambda }_{1}$ does not contribute to the auto-correlation function because it yields a factor that is proportional to the equilibrium expectation value of the velocity and hence vanishes. Therefore the period averaged velocity auto-correlation function simplifies to read

$$\begin{eqnarray}\langle {\bf{v}}({nT}){\bf{v}}(0)\rangle & = & \displaystyle \sum _{\kappa =2}^{3}{\lambda }_{\kappa }^{n}\displaystyle \sum _{\alpha ,\beta }{v}_{\alpha }{v}_{\beta }{{\rm{\Pi }}}_{\alpha ,\beta }^{(\kappa )}{p}_{\beta }^{\mathrm{st}}\\ & = & 2{v}^{2}{\lambda }_{2}^{n}{p}_{1}^{\mathrm{st}},\end{eqnarray} \tag{ 18 }$$

where we inserted the velocities ${v}_{1}=v$, ${v}_{2}=0$ and ${v}_{3}=-v$ and made use of the explicit form of the projection matrices ${{\rm{\Pi }}}^{(\kappa )}$ as given in the appendix. Note that it does not depend on ${\lambda }_{3}^{n}$. This is a consequence of symmetry of the matrix M and the resulting symmetry of the eigen-projector ${{\rm{\Pi }}}^{(3)}$, see equation (A.6). These symmetries reflect the symmetry of the model with respect to two states with non-zero velocities ${v}_{+}$ and ${v}_{-}$.

Zoom In Zoom Out Reset image size

The sum on n is readily performed to yield for the diffusion matrix

$$\begin{eqnarray}D & = & 2{v}^{2}T\displaystyle \frac{{p}_{1}^{\mathrm{st}}}{1-{\lambda }_{2}}\\ & = & {v}^{2}T\displaystyle \frac{1-k}{(2-k-q-r)(1-q+r)},\end{eqnarray} \tag{ 19 }$$

where we used the equations (15) and (A.3). In figure 7 the resulting expression is compared with the diffusion constant obtained from simulations of the Langevin equation (4). It turns out that (19) reproduces the non-monotonic behaviour of the diffusion constant with its maximum and minimum in good agreement with the simulation result. However, the absolute magnitude of D resulting from (19) is too large by approximately a factor of 1.5.

Zoom In Zoom Out Reset image size

In this work we studied the spreading of a cloud of inertial, massive Brownian particles independently moving in a periodic potential. As Brownian particles they are subject to a friction as well as to a fluctuating force, both being imposed by the interaction with a thermal heat bath held at a temperature θ. Moreover, they are driven by a force periodically varying in time. In the present investigation we considered a parameter regime in which the system behaves chaotic in the absence of thermal noise, i.e., if formally the temperature is set to zero. If the temperature is gradually increased the diffusion constant characterising the spreading rate of the particle cloud first increases until it reaches a local maximum and then decreases again until a minimum is reached. A further increase of the temperature causes a continuous growth of the diffusion constant which eventually becomes proportional to temperature.

As the reason for the counter-intuitive non-monotonic behaviour of the diffusion coefficient versus temperature we identified the dynamics of the populations of certain regions in phase space. The latter contain certain deterministically unstable periodic orbits. Unstable periodic orbits are known to constitute the backbone of deterministic chaos [58]. They allow one to reconstruct the chaotic dynamics of a system in a hierarchical way. Using a simplified model we considered three unstable periodic orbits. Two of them move with opposite temporal period averaged velocities and a third one that is locked. The frequencies of transitions change with temperature, thereby leading to a non-trivial dependence of the diffusion coefficient on temperature. Considering Markovian transitions between these states with transition rates estimated from numerical simulations of the Langevin equation (4) then leads to a qualitatively agreeing description, however, with a too large diffusion coefficient in the considered range of temperatures. This discrepancy may have different reasons. First of all, non-Markovian memory effects might lead to a smaller diffusion constant. Additionally, a more systematic partitioning of the phase space [59, 60] and the inclusion of further unstable periodic orbits likely will improve the agreement. Here, however, we have limited ourselves to the generic simplest stochastic model with three states in discrete time. The presented analysis has been restricted to one set of the system parameters. It is not the only one and exceptional region where a non-monotonic dependence of the diffusion coefficient D on temperature θ is detected. We have found other regimes but are not presented here because the temperature dependence of D and its mechanisms are similar.

Non-monotonic diffusion was also recently found in a rocked ratchet [45]; i.e. a driven Brownian particle moving in a periodic potential lacking any mirror symmetry. Due to the absence of the latter property a non-trivial directed transport effect $\langle {\bf{v}}\rangle \ne 0$ may emerge in this system. Contrary to the case considered in this paper, the deterministic limit of the mentioned regime is non-chaotic and possesses three stable coexisting states. The normal diffusion of the particle is caused by thermal noise produced by a heat bath in the small-to-moderate temperature limit. In the limit of weak noise a three-state modelling along the reasoning put forward here might as well provide satisfactory insight into the diffusion behaviour there.

In contrast, the present setup exhibits chaotic motion in the deterministic limit leading already to a finite diffusion constant $D\ne 0$ at $\theta =0$. At low temperatures the combination of chaotic motion with very weak noise leads to an increase of the diffusion. With increasing noise the conditional probability to stay in the locked state eventually increases whereas that of the running state decreases going hand in hand with a decrease of the diffusion constant. For this particular mechanism to work all terms of the Langevin equation (4) are relevant. Without the inertial term, i.e. in the overdamped limit, the dynamics take place in a two-dimensional state space and therefore cannot display chaos. Likewise, the absence of the periodic driving also would restrict the system to a two dimensional phase space without the option of exhibiting chaos. Moreover, the periodic forcing drives the system out of thermal equilibrium in which the diffusion constant is a strictly monotonic function of temperature. For the case with vanishing friction the model becomes simple; then a dose of finite noise corresponds to infinite temperature rendering any bounded potential ineffective. In this sense with our setup in equation (4) we deal with a minimal model for the peculiar diffusive behaviour.

Since the latter is readily implemented in many diverse physical situations [31–41] we expect that a similar non-monotonic behaviour of the diffusion may be observed in various contexts recurring in the most diverse areas of science such as physics, chemistry, biology, engineering, computer science or even sociology.

We expect that our findings can be experimentally corroborated with any of the physical systems mentioned in the introductory part of the article. As a possibility we propose a promising setup for this purpose, namely the dynamics of cold atoms dwelling optical lattices [61]. These systems are known for their high tunability, providing a precise control of the amplitude and period of a defect-free symmetric, spatially periodic optical potential. We consider the simplest and the most common model of a dissipative optical lattice consisting of atoms of mass m with a two level structure illuminated by the optical molasses, i.e. by counter-propagating light fields with the same frequency tuned slightly below the electronic transition of the atoms. This trap generates the potential force $\mbox{--}U^{\prime} (x)$ with precisely adjusted period and amplitude. The red-detuning of the laser fields leads to the Doppler cooling mechanism [62] which can be described as a classical damping term $\propto -\dot{x}$. The random photon absorption and re-emission events can be modelled by the Gaussian white noise $\xi (t)$ with noise intensity Q. Thus, the light beams play a role of the bath at an effective temperature to which atoms are coupled to [63]. Temperature of the atoms is mainly controlled by the molasses frequency or by heating of the optical lattice due to the noisy electronic phase shifting of the beam, see [63] and equation (6) therein.

Finally, in order to generate the external driving a time dependent phase modulation $\varphi (t)=a\cos (\omega t)$ can be applied to one of the lattice building beams. In the laboratory reference frame such a laser configuration generates a moving optical potential. In the co-moving frame the potential is static and the atoms experience an inertial force which is proportional to the time dependent phase modulation $\varphi (t)=a\cos (\omega t)$ [64, 65]. In this way, the predicted non-monotonic behaviour of the diffusion constant as a function of temperature should manifest itself in an optical lattice experiment. Last but not least, we stress that the detected peculiar diffusive behaviour manifesting in non-monotonic temperature dependence of the diffusion coefficient can be observed in deep optical lattices. This non-intuitive diffusion behaviour is therefore quite distinct from the well explored regime with shallow potentials where typically anomalous diffusion processes occur [66, 67].

This work was supported in part by the MNiSW programme via the Diamond Grant (JS) and the NCN grant 2015/19/B/ST2/02856 (J Ł). PT thanks the Foundation for Polish Science (FNP) for granting him an Alexander von Humboldt Honorary Research Fellowship.

The spectral representation (17) of the matrix ${\bf{M}}$ defined in equation (11) is given by the eigenvalues ${\lambda }_{\kappa }$ reading

$$\begin{eqnarray}&&{\lambda }_{1}=1\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{\lambda }_{2}=q-r\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{\lambda }_{3}=q+k+r-1\end{eqnarray} \tag{ A.3 }$$

and the projection matrices

$$\begin{eqnarray}{{\rm{\Pi }}}^{(1)}=\displaystyle \frac{1}{2-q-k-r}\left(\begin{array}{ccc}\tfrac{1}{2}(1-k) & \tfrac{1}{2}(1-k) & \tfrac{1}{2}(1-k)\\ 1-q-r & 1-q-r & 1-q-r\\ \tfrac{1}{2}(1-k) & \tfrac{1}{2}(1-k) & \tfrac{1}{2}(1-k)\end{array}\right),\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}{{\rm{\Pi }}}^{(2)}=\displaystyle \frac{1}{2}\left(\begin{array}{ccc}1 & 0 & -1\\ 0 & 0 & 0\\ -1 & 0 & 1\end{array}\right),\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}{{\rm{\Pi }}}^{(3)}=\displaystyle \frac{1}{2(2-q-k-r)}\left(\begin{array}{ccc}1-q-r & -1+k & 1-q-r\\ -2(1-q-r) & 2(1-k) & -2(1-q-r)\\ 1-q-r & -1+k & 1-q-r\end{array}\right).\end{eqnarray} \tag{ A.6 }$$