Anomalous transport in driven periodic systems: distribution of the absolute negative mobility effect in the parameter space

Absolute negative mobility is one of the most paradoxical forms of anomalous transport behaviour. At the first glance it contradicts the superposition principle and the second law of thermodynamics, however, its fascinating nature bridges nonlinearity and nonequlibrium in which these fundamental rules are no longer valid. We consider a paradigmatic model of the nonlinear Brownian motion in a driven periodic system which exhibits the absolute negative mobility. So far research on this anomalous transport feature has been limited mostly to the single case studies due to the fact that this model possesses the complex multidimensional parameter space. In contrast, here we harvest GPU supercomputers to analyze the distribution of negative mobility in the parameter space. We consider nearly $10^9$ parameter regimes to discuss how the emergence of negative mobility depends on the system parameters as well as provide the optimal ones for which it occurs most frequently.


Introduction
Nonlinear systems exhibit rich spectrum of unusual behaviour which is absent in their linear counterparts [1]. It is rooted in the fact that they are exempted from the superposition principle which loosely speaking tells that the response of linear system caused by two or more forces is the sum of the reactions that would have been induced by each of them individually. This property opens a new avenue for the emergence of remarkable phenomena like chaos, in which deterministic evolution of the system is completely not predictable [2] or multistability, when several stable states coexist in the setup dynamics [3,4].
Similarly, when the system is taken out of thermal equilibrium, monumental Thermodynamic Laws and various symmetries such as the detailed balance lose their validity. Solely this remark opens a new landscape of phenomena that despite many years of active research in nonequilibrium statistical physics still remains a terra incognita. Yet, some progress in exploring this fascinating ground has been achieved in the form of understanding effects like, for instance, stochastic resonance [5], anomalous diffusion [6,7,8,9], noise assisted transport [10,11] or deriving various celebrated fluctuation theorems [12,13,14] which bridge physics in and out of equilibrium.
In this work we unite these two worlds of nonlinearity and nonequilibrium that inspired the entire research fields which, as outlined above, have been intensively explored over recent decades. In doing so we investigate the form of anomalous transport behaviour, namely the absolute negative mobility, in which the particle moves in the direction opposite to the net acting force around zero bias [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. This phenomenon rests on the combination of both nonlinearity as well as nonequilibrium and cannot emerge without these two properties in a one-dimensional system [18,19].
Therefore we consider a paradigmatic model of nonequilibrium statistical physics, namely, the nonlinear Brownian motion in a driven periodic system. Since it possesses a multidimensional parameter space which has been too complex to explore systematically, so far research on the absolute negative mobility has been limited mostly to the single case studies. In contrast, in this work we exploit the state of the art computer simulations to analyze the distribution of absolute negative mobility effect in the parameter space by harvesting the power of GPU supercomputers. This innovative method [37] allowed us to consider nearly 10 9 parameter regimes to draw a number of important qualitative and quantitative conclusions about the emergence of absolute negative mobility. In particular we provide parameters of the model which are optimal for the occurrence of this anomalous transport behaviour.
The paper is organized as follows. In Sec. 2 we recall the formulation of the model, introduce the dimensionless quantities and discuss differences between the two most common scaling regimes. In Sec. 3 we briefly review the state of the art of the absolute negative mobility effect. The next Sec. 4 contains the description of employed research methodology. In Sec. 5 we present the results of our simulations. First we discuss the qualitative dependence of the absolute negative mobility on the model parameters. Then we elaborate on the distribution of negative mobility effect in the parameter space. Section 6 provides a summary and conclusions.

Model
In this work we consider the Langevin equation describing the dynamics of a Brownian particle dwelling in a one-dimensional spatially periodic potential [31]. The kinetics of the particle depends on its mass M and friction coefficient Γ. We assume that the potential is symmetric and has a spatial period L, namely The particle is driven by a harmonic force A cos (Ωt) as well as a constant bias F . The system is coupled to a thermostat of temperature T . Thermal fluctuations are modelled by δ-correlated Gaussian white noise ξ(t) of vanishing mean, i.e.
Such a model can be expressed by the following Langevin equation [31] Mẍ where dot means differentiation with respect to time t and x is the position of the particle. The factor √ 2Γk B T follows from the fluctuation-dissipation theorem [38] and ensures the correct Gibbs equilibrium state for the free particle. The potential, the harmonic force and thermal fluctuations are symmetric with respect to time and space, so the only perturbation that breaks the symmetry of Eq. (3) and allows for the emergence of directed transport is the constant force F .
There are many physical systems that can be modeled in terms of the dynamics given by Eq. (3) including superionic conductors [39,40], dipoles rotating in external fields [41], charge density waves [42], Josephson junctions [43,44] and its variations like SQUIDs [45,46] as well as cold atoms dwelling in optical lattices [47,48], to name but a few.
In physics only the relations between characteristic scales of time, length and energy are relevant for the progress of observed phenomena, but not their absolute values. This fact suggests recasting of Eq. (3) into its dimensionless form in which all quantities are expressed as combinations of characteristic scales for the system. This procedure makes the analysis independent of the experimental setup provided that the mathematical structure of Eq. (3) is preserved. It also allows to reduce the number of free parameters appearing in the model. The scaling is based on choosing the length and time scales. The most obvious selection for the characteristic length is the potential period L. Depending on the definition of the time scale one arrives at different scalings. Below two most common variants are presented and the differences between them are elaborated.

First scaling with mass m = 1
The first time scale follows from the equation for frictionless movement of a particle in the periodic potential In such a case the time unit can be extracted as follows [49] It is related to the period of linearized oscillations within one potential well. We define the dimensionless particle coordinate and time aŝ Under such a choice the original Eq. (3) takes the following forṁ where v 1 = dx/dt 1 . The rescaled potential is defined aŝ and possesses the period equal to unity. The dimensionless thermal noise readŝ and still represents the δ-correlated Gaussian white noise of zero mean, c.f. Eq. (2). Its rescaled intensity D is a ratio of the thermal and potential barrier energies The external force parameters read The reader can note that in this scaling one may introduce also the dimensionless mass but it is fixed to m = 1. The rest of quantities appearing in Eq. (7) explicitly depend on the chosen time unit, namely The dimensionless friction coefficient γ can be expressed as a ratio of two characteristic time scales where τ 0 = M/Γ stands for the so-called Langevin time, i.e. the relaxation time for the velocity of a free Brownian particle. Nevertheless it is instructive to note that γ is proportional to the actual physical friction coefficient Γ.

Second scaling with friction coefficient γ = 1
Another time scale can be extracted from the equation of an overdamped motion of a particle in the periodic potential The time unit is then [49] It scales the characteristic time for the overdamped particle to move from the maximum to minimum of the potential U (x). In this case the original Eq. (3) transform as follows where now v 2 = dx/dt 2 . The parameters that do not depend on the time scale are the same as in the first scaling. The new quantities are Now the dimensionless mass m is proportional to the physical mass M and expressed as a ratio of two characteristic time scales τ 0 = M/Γ and τ 2 while the friction coefficient formally scales to γ = 1.

Differences between scalings
Since the two above presented variants of the rescaled dynamics refer to the same model given by Eq. (3) the parameter definitions appearing in both of them must be associated with each other. One can show that they obey the following relations Obviously, the inverse is also true It is important to note that the velocities v 1 and v 2 are calculated as derivatives ofx with respect to t 1 and t 2 , respectively. Thus, they cannot be compared directly, but a proper rescaling is needed As long as fixed values of γ and m are considered both scalings are equivalent. One can always recalculate all quantities from one scaling to the other. It might seem that thermal fluctuations terms are not equivalent since in the first scaling the thermal noise prefactor depends on both γ and D, whereas in the second one it is determined only by D. However, since the white noise is δ-correlated, one can show that which means that these terms are also equivalent. Nevertheless, there is one very important difference between these scalings which is the relation of the dimensionless time t 1 and t 2 to the actual physical time t in Eq. (3). In the first scaling the time unit τ 1 does not depend on the friction coefficient Γ, but on the mass M . It means that one can easily investigate the impact of damping Γ in this scaling without changing the time scale. Analyzing the influence of mass M is also possible, but for every value of M the time scale in the system would be different. Interpretation of the results in such an approach would be hardly possible. Similarly, the time unit τ 2 in the second scaling does not depend on the mass M , but on the damping Γ, which makes it unsuitable for studying the impact of friction Γ on the dynamics. This difference has the most profound consequences in the limiting situations when either γ or m approaches zero which is the case in the Hamiltonian and overdamped dynamical regimes, respectively. When γ → 0, then from the relation t 2 = t 1 /γ it follows that t 2 → ∞. This means that taking the limit γ → 0 in the first scaling is not equivalent to requiring m → ∞ in the second one. Similarly, the relation v 2 = γv 1 implies that the velocity in the second scaling would then tend to zero. For these reasons, the approach to the case of Hamiltonian dynamics in which formally γ = 0 have to be analyzed by taking the limit γ → 0 in the first scaling. The analogous discussion could be repeated for the overdamped regime when m = 0 that would lead us to the conclusion that such a scenario must be investigated by performing the limit m → 0 in the second scaling what is not equivalent to the case γ → ∞ in the first one.
In the remaining part of the article we will use only the dimensionless variables and therefore for simplicity we omit the hat-notation, i.e. we will write x instead ofx and so on.

Quantity of interest
The most fundamental quantity characterizing the directed transport is the average velocity defined as where · indicates averaging over all thermal noise realizations as well as initial conditions for the particle position x(0) and velocity v(0). The latter is obligatory especially in the limiting case of deterministic dynamics D = 0 when the ergodicity of the system may be broken and consequently the results are affected by those initial conditions [7,50]. The above definition is independent on the scaling selection. The time t can be either t 1 or t 2 and the velocity v may be v 1 or v 2 .

Absolute negative mobility effect
To make the paper self-contained in this part we briefly review the state of the art of the absolute negative mobility effect for the studied system.
We start with the observation that the underlying symmetries of the Langevin Eq. (3) imply that the average velocity v is an odd function of the constant bias f , [51]. Therefore the directed transport cannot emerge in absence of the static force f as then v (0) ≡ 0. We define the mobility µ [20] of the particle as v (f ) = µ(f )f (22) to describe its ability to move through the medium in response to the acting force. Typically the resultant particle displacement follows the direction of the applied bias µ(f ) > 0. The corresponding average velocity v renders a nonlinear function of the constant force f and it is expected that v increases for growing f . In the linear response regime the velocity v = µ 0 f is a linear function of the force f with the constant mobility coefficient µ 0 [20]. The term absolute negative mobility refers to the paradoxical case when the net particle movement is opposite to the direction of the static load around zero bias [31], i.e.
The counterintuitiveness of this phenomenon follows from the fact that in a linear system the influence of all forces can be analyzed separately and the collective effect is simply a sum of responses induced by all perturbations. Since the harmonic driving a cos (ωt) and thermal fluctuations ξ(t) have a vanishing mean, in a linear system their contribution to the net movement would be zero and the particle would follow the direction of the constant force f , what implies that the absolute negative mobility effect would not emerge. However, the potential U (x) is periodic and hence gives rise to a nonlinear system, where the superposition principle is no longer valid. The nonlinearity is another necessary condition for µ(f ) < 0 to occur [18,19].
One can argue that the net movement of the particle in a direction opposite to the constant force contradicts the Le Chatelier-Braun's principle [18,19]. This law is however no longer valid for systems out of equilibrium. Therefore the key requirement for the occurrence of the absolute negative mobility is that the system is driven far from thermal equilibrium into a nonequilibrium state [18,19]. In the considered model it is guaranteed by the presence of the external time periodic driving a cos (ωt).
Last but not least, it has been already shown in literature that the absolute negative mobility does not emerge in the limiting case of one-dimensional overdamped (m = 0) and Hamiltonian (γ = 0) regimes [27]. Omitting the dissipative term −γv is equivalent to the situation in which the system is coupled to infinitely hot bath for which the potential U (x) term becomes negligible. As it was discussed above the latter is essential for the emergence of the absolute negative mobility effect.
Three different mechanisms responsible for this counterintuitive phenomenon are currently known -deterministic chaotic, deterministic non-chaotic and thermal noise induced [18,19,27]. The deterministic dynamics given by Eq. (7) or Eq. (15) with D = 0 can be recasted into a set of three autonomous differential equations of the first order for which the corresponding phase space is three-dimensional {x, v, z = ωt} being the minimal requirement for the system to display chaotic evolution [1]. In such a case the absolute negative mobility effect emerges as a result of the subtle interplay between coexisting attractors and transient chaos [19]. Recently it has been demonstrated that this phenomenon can occur also in the deterministic system given by Eq. (15) which is in the non-chaotic dynamical regime [27], i.e. it exhibits regular attractors transporting the particle in the direction opposite v < 0 to the applied bias f > 0. Finally, the absolute negative mobility can be triggered solely by thermal equilibrium fluctuations [18]. In such a case in the deterministic dynamics the absolute mobility of the particle is positive µ(f ) > 0 but it takes negative values µ(f ) < 0 for certain temperature D window.

Methods
The Fokker-Planck equation corresponding to Eq. (7) or Eq. (15) is the second order parabolic partial differential equation with a nonlinear and time periodic drift coefficient due to the presence of potential U (x) and driving a cos (ωt), respectively. For this reason its solution is unattainable analytically and in order to analyze the transport behaviour of driven Brownian particle we carried out comprehensive numerical simulations. The studied system possesses a complex five dimensional parameter space {γ or m, a, ω, f, D}. Its systematic exploration was not possible until very recently due to limited computational capabilities of modern hardware and lack of innovative implementations of simulation methods. We performed numerical analysis by harvesting the GPU supercomputers [37] that allowed us to draw both qualitative and quantitative conclusions about the emergence of absolute negative mobility phenomenon in the parameter space. The latter were picked from a cuboid in {γ, a, ω} or {m, a, ω} subspace containing 400 × 400 × 400 values. This volume, with 14 combinations of f and D values, resulted in nearly 10 9 parameter regimes per each considered scaling.
We employed a weak second order predictor-corrector scheme [52] to simulate stochastic dynamics given by Eq. (7) or Eq. (15). The time step of integration was scaled as h = 10 −2 × T where T = 2π/ω is the fundamental period of the external driving a cos (ωt). The average velocity v was calculated over the ensemble of 1024 system trajectories each starting with different initial conditions. The initial positions x(0) and velocities v(0) were uniformly distributed over the intervals [0, 1] and [−2, 2], respectively. All trajectories lasted for 10 3 periods T of the external driving and spanned the interval used for calculating the time average in the definition of directed transport v .
Since the latter quantity is invariant under changes of the sign of a, we restrict our analysis only to positive values a ∈ [0, 25]. The limit can be put also on the frequency ω. For ω → 0 the adiabatic approximation may be employed and the velocity follows the external force. On the other hand, for ω → ∞ the average velocity can be expressed via the Bessel functions [43] and the absolute negative mobility does not emerge. Consequently, in our simulations we analyzed the interval ω ∈ [0.01, 20]. Moreover, since the directed transport v is an odd function of the static bias f , it is sufficient to consider only positive values f > 0 of the latter parameter. It is intuitive that when the constant force f is much larger than other perturbations it dominates the dynamics and the Brownian particle velocity follows its direction. Therefore the values of f were chosen from the interval f ∈ [0.01, 1.5]. Finally, we remind that the absolute negative mobility does not occur in the overdamped m = 0 and Hamiltonian γ = 0 regimes. Consequently, the simulations were performed for γ ∈ [0.1, 10] and m ∈ [0.01, 10]. In majority of the investigated parameter sets the thermal noise had a destructive influence on the occurrence of absolute negative mobility. For this reason most of the regimes corresponded to the deterministic system with D = 0, with a few runs for temperature up to D = 0.01. Overall, the analyzed parameter subspace allowed to cover almost entire range of values for which the net movement of the particle is in the direction opposite to the applied force and the absolute negative mobility phenomenon emerges.

Results
Although the simulations were performed for a wide range of parameters, the below presented results show only some subsets of the studied region. The presented areas were chosen to reflect the general dependence of the directed transport v on the model parameters.

First scaling with mass m = 1
We start our analysis with the investigation of the absolute negative mobility effect in the subspace of parameters characterizing the external harmonic driving. In Fig. 1 we depict the influence of dissipation γ on the absolute negative mobility phenomenon. The average velocity v 1 of the Brownian particle as a function of amplitude a and angular frequency ω 1 of the harmonic driving for four different values of γ with f = 0.1 and D = 0 is presented. One can observe an arc-like area of absolute negative mobility and its evolution with the change of γ. When γ increases, this region moves towards lower a and ω 1 . For γ = 1.1749 this area is maximal. When γ increases further the whole structure vanishes and adjacent regions of positive velocity become more intense. There is also a small area of relatively low negative velocity below the arc, accompanied by a similar region of positive velocity. Fig. 2 presents the impact of the amplitude a of the external harmonic driving on the absolute negative mobility effect. This panel depicts the average velocity v 1 of the Brownian particle as a function of friction coefficient γ and angular frequency ω 1 of the harmonic driving for different values of a with f = 0.1 and D = 0. There are two main areas of V-like shape, one with negative, the other with positive velocity. In the places where they overlap, they compensate each other and the velocity in these overlapping regions is much closer to zero than in the non-overlapping parts. When a grows the structure moves towards smaller γ and greater ω 1 . Moreover, for increasing a the velocity in both areas assumes more extreme values.
In Fig. 3 we illustrate the response of the absolute negative mobility effect induced by a change of the external driving frequency ω 1 . The directed transport v 1 of the Brownian particle versus the amplitude a and friction coefficient γ for selected values of ω 1 and fixed f = 0.1 and D = 0 is presented there. For small ω 1 there are ray-like structures with positive velocity with very thin stripes of the absolute negative mobility along them. It suggests that for small ω 1 there is almost linear relation between a and γ for which the transport occur. The ray-like pattern is deformed if ω 1 grows. First, the areas become larger and the absolute value of the average velocity increases. Second, they tend to orientate more horizontally indicating that if ω 1 is increased the absolute negative mobility is expected to emerge for greater a but smaller γ. At some point the ray-like structure is changed to V-like, similar to the one visible on Fig. 2. The areas of positive and negative velocity occur in pairs and partially overlap. Next, Fig. 4 shows the impact of the static bias f on the absolute negative mobility. The panel presents the average velocity v 1 of the Brownian particle as a function of the amplitude a and the frequency ω 1 of the harmonic driving for different f with γ = 1.1749 and D = 0. For small f the absolute negative mobility areas are barely visible, similarly to the regions of positive ones. This is consistent with the already mentioned requirement which states that the external constant force f breaks the symmetry of the system and is necessary to induce the transport. When f increases, the region of absolute negative mobility broadens, however for large f the absolute value of negative velocities decreases and the structure finally vanishes. The reason is that when the static bias f is much larger than other forces in the system, it dominates the dynamics and the absolute negative mobility effect cannot occur. The region of this anomalous transport behaviour is largest for f = 0.5, however for f = 0.1 there are many finer structures that vanish for greater f .
Finally, in Fig. 5 we discuss the influence of temperature D on the absolute negative mobility. The average velocity v 1 of the Brownian particle versus the amplitude a of the harmonic driving and dissipation γ is presented there for various values of D with γ = 1.1749 and f = 0.1. The reader can observe that temperature growth causes blurring of the structures of both positive and negative velocity. The fine details disappear first while the larger one are still present in the noisy case, however, for high enough temperature they all vanish. It is expected as eventually thermal noise dominates the dynamics and the particle behaves as the free one. Typically temperature has destructive impact on the emergence of absolute negative mobility in the parameter space.
Now we turn to the analysis of different parameter subspace associated with the propelling force f and dissipation γ. In Fig. 6 we present the average velocity v 1 of the Brownian particle versus the static bias f and damping γ for different values of the external driving amplitude a with ω 1 = 3.5 and D = 0. In all plots pairs of wedge-like areas with positive and negative velocity are visible. In places where they overlap they compensate each other. The borders of these regions are approximately linear, which suggest a simple relation between f and γ for which the directed transport occurs. The area of anomalous transport behaviour moves towards lower values of γ as f is increased.
In Fig. 7 we depict the directed transport v 1 for the same subspace, namely v 1 (f, γ) but for the fixed amplitude a = 24.4375 of the harmonic driving and different angular frequencies ω 1 of the latter perturbation. The reader can observe there the evolution of one wedge-like area of absolute negative mobility corresponding to alteration of ω 1 . When ω 1 increases, this structure is enlarged, with a maximum at ω 1 = 10.0, and then starts to disappear. Likewise, there is an optimal value of the static bias for which the absolute negative mobility emerges for the broadest interval of γ. If ω 1 grows the region of negative velocity is moved towards smaller dissipation γ. A closer look at Fig.  7 reveals barely visible areas of negative directed transport that look like copies of the main one, but shifted up or down and dimmed.
The above presented approach to analyze the dynamics of a driven Brownian particle in the multidimensional parameter space by dividing it onto several subspaces allowed us to draw a number of general qualitative conclusions regarding the emergence of the absolute negative mobility effect. Firstly, increasing the amplitude of the external perturbation applied to the particle, either static bias f or harmonic driving a, causes a shift of the absolute negative mobility regions towards lower values of dissipation γ. This fact can be seen explicitly in Fig. 2 or Fig. 6. Secondly, when the frequency ω 1 of the harmonic driving grows the regions of absolute negative mobility displace in the direction of smaller dissipation γ and greater constant bias f , c.f. Fig. 7. Thirdly, in most cases temperature D influences destructively the anomalous transport behavior which we exemplified in Fig. 5.
Exploiting the GPU supercomputers to investigate the absolute negative mobility of the driven Brownian particle let us to attain not only qualitative remarks about the emergence of this effect in the parameter space but also obtain the important quantitative results. We depict them in Fig. 8 where the fraction of the investigated space for which the absolute negative mobility phenomenon occurs (probability P ) versus different parameters of the studied model for several values of the static bias f is presented. Panel (a) illustrate the latter quantity as a function of dissipation γ. One can note that in the both extreme limits γ → 0 and γ → ∞ the anomalous transport behaviour completely disappears. It is consistent with the state of the art of this paradoxical effect. However, regardless of the magnitude of constant force f there is a common optimal value of γ ≈ 1.1749 for which P is maximal, i.e. the absolute negative mobility emerges most frequently in the parameter space. The reader can observe that the extremum of P is maximal for the bias f = 0.5. In panel (b) of the same Fig. we show P versus the amplitude a of the external harmonic driving. For a = 0 there is no transport in the negative direction which agrees with the statement that the harmonic driving is necessary for this anomalous behaviour to emerge. For all values of f the probability P initially rises when a increases from zero. Then it strongly depends on the value of f and can have both local maxima and minima. However, the general tendency is that for not too large a the absolute negative mobility occurs more frequently when the amplitude a is increased. In panel (c) the dependence of P on the frequency ω 1 is presented. For all values of the constant force f there are two intervals where P has pronounced extrema. Their locations vary for different f , but the first one is in the vicinity of ω 1 ≈ 3.5 and the other one is between ω 1 ∈ [7, 11]. In the former region the highest peak is reached for f = 0.1 while for the latter the optimal value of force is f = 0.5 for which the global maximum of P is attained when ω 1 ≈ 10.25. The higher the static bias f , the lower is the cutoff frequency ω 1 for which the absolute negative mobility ceases to exist. Regardless of the magnitude of the load f the probability P vanishes if ω 1 > 15. Finally, in panel (d) the impact of thermal fluctuations on the distribution P is illustrated for the fixed f = 0.1. For most of γ thermal noise have destructive influence on the occurrence of absolute negative mobility with the highest value of P in the deterministic case D = 0. However, careful inspection of the figure reveals that for small γ the probability P for D = 0.0001 is larger than for D = 0. This fact testifies that thermal fluctuations indeed can induce the absolute negative mobility of the driven Brownian particle.

Second scaling with friction coefficient γ = 1
For fixed values of dissipation γ and inertia m in the first and second scaling, respectively, one is able to transform the results obtained in each scaling to the other one by using the relations presented in Sec. 2.3. Therefore qualitative conclusions about the emergence  of absolute negative mobility in the parameter space which we made in the last section hold true provided that the relation m = 1/γ 2 is taken into account. It basically means that inertia m is inversely proportional to dissipation γ 2 . The exception from the above rule is when the directed transport v is studied in the parameter plane involving inertia m. Then the relation between the frequencies ω 1 = √ mω 2 tells that it is not possible to fix single ω 2 in such a way that it would correspond to ω 1 for all inertia m simultaneously and consequently there is no one-to-one correspondence between the analyzed parameter subspaces. For this reason it is not pointless to investigate the emergence of the absolute negative mobility effect also in the second scaling.
In Fig. 9 we show the probability P of emergence of this anomalous transport behaviour versus different quantity characterizing the studied model. These results are analogous to the corresponding ones depicted in Fig. 8 but now they are for the second scaling. In panel (a) P is plotted as a function of the dimensionless mass m. The reader can immediately realize that in both extreme limits m → 0 and m → ∞ the absolute negative mobility ceases to exist. Again it is consistent with the current state of knowledge. There are two values of m for which the probability P displays pronounced  maximum. The first one is m ≈ 0.16 and the second m ≈ 0.6. The corresponding optimal forces reads f = 0.1 and f = 0.5 for the former and latter inertia m, respectively. In panel (b) of the same Fig. we investigate the distribution P versus the amplitude a of the harmonic driving. For every static bias f there is a range of a for which the absolute negative mobility emerges most frequently, e.g. for f = 0.1 it is a ∈ [11,16]. It is expected that the probability P will eventually vanish for growing a. It is due to the fact that in such a case other perturbations in the dynamics which are necessary for the anomalous transport to arise will be negligible. The most significant difference between plots of P can be seen on Fig. 9 (c) and 8 (c), where the dependence of P on the angular frequency ω is captured. The similarity between them is that for each f there is the cutoff frequency ω above which the absolute negative mobility does not occur. Naturally it is because the dimensionless frequencies ω 1 and ω 2 are proportional to each other. Moreover, minor changes of ω in both cases can intensify or weaken the emergence of absolute negative mobility. However, in contrast to the panel Fig. 8 (c), there exists a common characteristic frequency ω 2 ≈ 11 for which the probability P displays the pronounced maximum for every studied force f . Finally, in panel (d) the influence of thermal noise intensity D on the occurrence of absolute negative mobility P is depicted. For the whole range of the investigated inertia m temperature D acts destructively on this anomalous transport feature.

Conclusions
In this work we revisited the problem of anomalous transport in driven periodic systems. Specifically, we considered a paradigmatic model of nonequilibrium statistical physics, namely, an inertial Brownian particle moving in a symmetric spatially periodic potential which in addition is exposed to both an external harmonic driving and a static load. We focused on the particular instance of anomalous transport behaviour in the form of the counter-intuitive absolute negative mobility effect, in which the net particle movement is opposite to the direction of the applied biasing force around zero load. It has been under investigation for many years, however, since the studied system possesses an extremely rich five-dimensional parameter space which for the same time has been too complex to explore systematically, research to date has been targeted mostly at the single case studies. In contrast, in the present work we exploited the state of the art computer simulations to analyze the emergence of absolute negative mobility phenomenon in the parameter space with unprecedented resolution of nearly 10 9 parameter regimes by harvesting the power of GPU supercomputers. This innovative method allowed us to draw a few general qualitative conclusions about the absolute negative mobility. Moreover, the astonishing number of considered regimes made it possible to obtain also a number of important quantitative results.
In particular, we determined how the occurrence of absolute negative mobility effect changes under the influence of the system parameters. The amplitude of the external perturbation, either static bias f or harmonic driving a, induce a shift of the anomalous transport regions towards lower values of dissipation γ. When the frequency ω of the harmonic driving grows the occurrence of absolute negative mobility displace in the direction of smaller damping γ and greater constant force f . Thermal fluctuations D typically have destructive impact on this anomalous transport behaviour. Surprisingly, it turned out that regardless of the magnitude of constant force f there is a common optimal value of dissipation γ ≈ 1.1749 for which the absolute negative mobility emerges most frequently in the parameter space. In contrast, for inertia m the dichotomous behaviour was detected, i.e. depending on the static bias f there are two global optima for this effect. When the constant force f ≈ 0.1 is significantly smaller than the potential barrier the optimal mass reads m ≈ 0.16, whereas for the load f ≈ 1.0 comparable to the potential barrier it reads m ≈ 0.6. The general tendency is that the absolute negative mobility occurs more frequently when the amplitude a of the external driving is increased. The effect of external harmonic driving frequency ω is most complex, however, for many studied loads f the frequency ω ≈ 11 seems to be the most optimal choice with respect to the emergence of absolute negative mobility.
Summarizing, our findings for the paradigmatic model of nonequilibrium statistical physics can be straightforwardly corroborated experimentally with a wealth of physical systems including, among others, the Josephson junctions and cold atoms dwelling in optical lattices which nowadays are intensively studied. We believe that the results presented in this work will serve as a platform for both experimenters and theorists to further investigate anomalous transport processes occurring in driven periodic systems.