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

The spreading of a cloud of independent Brownian particles typically proceeds more effectively at higher temperatures, as it derives from the commonly known Sutherland-Einstein relation for systems in thermal equilibrium. Here, we report on a non-equilibrium situation in which the diffusion of a periodically driven Brownian particle moving in a periodic potential decreases with increasing temperature within a finite temperature window. We identify as the cause for this non-intuitive behaviour a dominant deterministic mechanism consisting of a few unstable periodic orbits embedded into a chaotic attractor together with thermal noise-induced dynamical changes upon varying temperature. The presented analysis is based on extensive numerical simulations of the corresponding Langevin equation describing the studied setup as well as on a simplified stochastic model formulated in terms of a three-state Markovian process. Because chaos exists in many natural as well as in artificial systems representing abundant areas of contemporary knowledge, the described mechanism may potentially be discovered in plentiful different contexts.


Introduction
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,3,4,5,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,9,10,11,12,13,14,15,16].
Irregular behaviour generally implies also unpredictability, which, when it comes to applications one typically attempts to avoid or at least to minimize. 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 characterized by the diffusion coefficient D increases with the temperature of the medium surrounding the system [18,19,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 ∆x 2 (t) grows asymptotically according to a power law ∆x 2 (t) ∼ t α with α = 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,24,25,26,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.

Model
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 0 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 (Ω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 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 ξ(t) of zero mean and unit intensity, i.e., The parameter Γ denotes the friction coefficient and k B is the Boltzmann constant. The noise intensity factor 2Γk B θ 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 Eq. (1), i.e., when A = 0 and in presence of periodic boundary conditions. We focus on the archetypical non-linear situation, namely, a sinusoidal potential with a barrier height 2∆U , Upon introducing the period L and the parameter combination τ 0 = ΓL 2 /∆U as units for length and time, respectively, Eq. (1) is written in a dimensionless form [44], reading, wherex = x/L andt = t/τ 0 . The dimensionless potentialÛ (x) = U (x)/∆U = U (Lx)/∆U =Û (x + 1) assumes the unit period L = 1 and the potential amplitude ∆Û = 1. The remaining re-scaled parameters are the mass m = (1/Γτ 0 )M , the amplitude a = (L/∆U )A and the angular frequency ω = τ 0 Ω. The dimensionless thermal noise readsξ(t) = (L/∆U )ξ(t) = (L/∆U )ξ(τ 0t ) and assumes the same statistical properties, namely it is Gaussian with ξ (t) = 0 and ξ (t)ξ(ŝ) = δ(t −ŝ). The dimensionless noise intensity Q = k B θ/∆U is given by the ratio of thermal energy k B θ and half of the activation energy the particle needs to overcome the original potential barrier 2∆U . From now on, we will use only the dimensionless variables and therefore shall skip the "hat" for all quantities appearing in Eq. (4). The observable of foremost interest in this study is the diffusion coefficient D which characterizes the spread of trajectories and fluctuations around the average position of the particle, namely [45,46], where the averaging · is over all realizations of thermal fluctuations as well as over initial conditions for the position x(0) and the velocityẋ(0). The latter is necessary because in the deterministic limit of vanishing thermal noise intensity Q → 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 (ωt) as well as the friction termẋ in Eq. (4) the particle velocityẋ(t) approaches a unique asymptotic state in which it is characterized 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 ẋ(t) assumes for a long time the form where v is the directed transport velocity. Note that all three terms entering the right hand side of Eq. (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 ξ(t) vanishes according to Eq. (2). Moreover, all three elements are symmetric and in consequence there is no directed transport in the long time stationary regime, i.e. v = 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 ) T 0 dt v(t) = 0. Therefore it is useful to consider also the period averaged velocity v(t) as where T = 2π/ω denotes a period of the external driving a cos (ω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 v(kT ) = (1/T ) where we assumed that the particle starts at 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 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 Eq. (9) the left hand side approaches x 2 (nT ) = 2nT D and hence the diffusion coefficient D can be expressed in terms of the temporal period averaged velocity autocorrelation function, yielding An equivalent expression is given in [50]. The infinite time limit in Eq. (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 ) × T [51]. Because Eq. (4) is a second-order differential equation, we need to specify two initial conditions, namely x(0) andẋ(0). We choose x(0) andẋ(0) to be equally distributed over the intervals [0, 1] and [−2, 2], respectively. Our quantities of interest were averaged over 10 3 − 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 3 times as compared to a common present-day CPU method [51].

Results
The dynamical system described by Eq. (4) exhibits an extremely rich behaviour as a function of the four dimensionless parameters {m, a, ω, Q} with dimensionless friction at Γ = 1. A general overview is provided, e.g. with Ref. [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 Fig. 1. The corresponding parameters of the dynamical system (4) are chosen as {m = 0.9, a = 8.7, ω = 0.275}. At each depicted noise strength Q ∝ θ the diffusion constant D was estimated from 10 5 trajectories by means of Eq. (5). At low noise intensity D increases with Q until it reaches a local maximum at Q ≈ 2 · 10 −5 . From there it decreases to a minimum at Q ≈ 5 · 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 Fig. 1. The decrease of the diffusion constant with increasing temperature Q ∝ θ is counter-intuitive. It stays in clear contrast with the Einstein relation D ∝ θ 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,54,55]. Because relevant structural elements of the stochastic dynamics described by Eq. (4) are determined by its deterministic properties as the first step of our analysis we consider the noiseless case.

Noiseless dynamics: Q = 0
For the set of parameters presented in Fig. 1, the deterministic system (when Q = 0) exhibits chaotic behaviour with a dense set of unstable periodic orbits. For Q > 0, the system 'feels' some unstable orbits which play an important role in controlling diffusion  7). The bifurcation diagrams reveal three qualitatively different dynamical regimes. -For the smallest displayed values of a < a 0 = 8.6867 as well as for a > 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 Eq. (5) has a finite value larger than zero. In the window between a 0 and a c the motion is phase locked with the corresponding winding number w = ±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 v(kT ) ≈ ±0.4 visible in panel (b) in the corresponding parameter window. In panel (a) two period doubling cascades start at the parameter value a 1 = 8.6878 and terminate at a ∞ = 8.6886. In the small parameter interval between a ∞ 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 v(kT ) corresponding to integer winding numbers characterizing unstable periodic orbits. Particularly large is the probability to   Fig. 1. In all three cases the velocity stays predominantly close to v(kT ) = ±0. 4 and v(kT ) = 0. In the majority of cases a trajectory does not stay for a longer period of time close to one of the dominant velocities but rather jumps in the vicinity of another one, without following any obvious rule.
find temporal period averaged velocities v(kT ) ≈ ±0.4, corresponding to the winding numbers w = ±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.

Influence of thermal noise: Q > 0
In Fig. 3 we present the influence of thermal noise on velocity of the Brownian particle. We depict time series of period averaged velocities v(kT ) resulting from a single trajectory in the large time limit for the same set of parameters specified in Fig.  1. Panel (a) exemplifies the deterministic case, with Q = 0. Regions close to the dominant winding numbers w = ±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  Figure 4. The probabilities p R and p L for the particle to be in the running |v(kT )| ≈ 0.4 and in the locked v(kT ) ≈ 0 state, respectively, are plotted against temperature of the system Q ∝ θ. The diffusion coefficient D is maximal when the difference between probabilities p R − p L has a peak, c.f. Fig. 1. When p R and p L intersect the diffusion coefficient D is minimal. Other parameters are the same as in By introducing the thresholds at v(kT ) = ±0.2 and counting the number of period averaged velocities which are between and outside these values one may estimate the probabilities p R = prob(|v(kT )| > 0.2) and p L = prob(|v(kT )| ≤ 0.2) for the occurrence of the running and locked states, respectively. These probabilities are displayed in Fig. 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 v(kT ) ≈ ±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 ) = n−1 k=0 v(kT ) with the diffusion constant D = 0.16T p 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 Fig. 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.

Approximate three-state model description of chaos induced non-monotonic diffusion
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 = ±9 corresponding to the velocities v ± ≈ ±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 ± . The sequence of period averaged velocities v(nT ) is then modelled by a Markov process in terms of transitions between these states. Because of the assumed symmetry the matrix M describing these transitions has the form where q denotes the conditional probability to remain staying in the running state v + , v + → v + , and likewise in v − , v − → v − ; further k, the conditional probability to remain staying in the resting state 0, 0 → 0; and r, the conditional probability of a transition between opposite running states v + → v − and v − → v + . Due to the assumed symmetry of the running states the transitions from the locked into a running state (0 → v ± ) is (1 − k)/2 and from a running into the locked state (v ± → 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 ≤ r ≤ 1 − q. The transition probabilities can be estimated from simulations of the Langevin equation (4) Figure 6. The directly estimated probability p R /2 (red) of a particle to be found in the running state as a function of the noise strength Q is compared with the respective stationary probability (15) (blue) following from the three-state model defined by the matrix (11) of transition probabilities. The agreement is excellent for a not too large noise intensity. For larger noise intensity the agreement is still qualitatively good.
The three independent rates k, q, r are displayed in Fig. 5. Finally, in the present model the states are "decorated" with the respective velocities {v ± , 0} such that the particle changes its position within a period by vT with v ∈ {v ± , 0}.
Collecting the probabilities of finding the particle at the time t = nT in either of the three states as a vector, p(n) = (p + (n), p 0 (n), p − (n)) one can write its dynamics as p(n + 1) = Mp(n) (12) having the formal solution with the initial probability p(0). The stationary distribution p st is invariant under transitions, and hence is the solution of given by In Fig. 6 the stationary probability p 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 Eq. (10) which is based on the stationary auto-correlation function of period averaged velocities. This correlation function follows from the three-state model as where α and β label the three velocity states and the matrix element (M n ) α,β determines the conditional probability to find the state labelled by α after n periods provided the system has started in the state β. Accordingly, p st β denotes the probability (15) to find β in the stationary state. In order to evaluate the n th power of the matrix M it is convenient to use its spectral representation reading where λ κ are the eigenvalues of M and Π (κ) are matrices projecting onto the respective eigenvectors. The eigenvalue λ 1 = 1 results as a consequence of the conservation of total probability. The corresponding matrix Π 1 projects any vector q onto the stationary probability p st , i.e. it acts as Π 1 q = p st α q α . Explicit expressions of the other eigenvalues and eigenprojectors are given in the Appendix Appendix A. Inserting the spectral representation (17) into the expression (16) one finds that the first eigenvalue λ 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 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 Π (κ) as given in the Appendix A. Note that it does not depend on λ n 3 . This is a consequence of symmetry of the matrix M and the resulting symmetry of the eigen-projector Π (3) , see equation (A6). These symmetries reflect the symmetry of the model with respect to two states with non-zero velocities v + and v − .
The sum on n is readily performed to yield for the diffusion matrix where we used the Eqs. (15) and (A.3). In Fig. 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.

Discussion
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 characterizing 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 v = 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 modeling along the reasoning put forward here might as well provide satisfactory insight into the diffusion behavior there.
In contrast, the present setup exhibits chaotic motion in the deterministic limit leading already to a finite diffusion constant D = 0 at θ = 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 Eq. (4) we deal with a minimal model for the peculiar diffusive behaviour.
Since the latter is readily implemented in many diverse physical situations [31,32,33,34,35,36,37,38,39,40,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 −U (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 ∝ −ẋ. The random photon absorption and re-emission events can be modelled by the Gaussian white noise ξ(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 Ref. [63] and Eq. (6) therein.
Finally, in order to generate the external driving a time dependent phase modulation ϕ(t) = a cos (ω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 ϕ(t) = a cos (ω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].