**1. Introduction**

Research in nonequilibrium statistical physics provides a wealth of intriguing dynamics in which phenomena that are forbidden in equilibrium states may emerge. Prominent examples include anomalous diffusion [1–4], Brownian yet non-Gaussian diffusion [5–9], noise-assisted transport [10,11], and negative mobility [12–15], to name only a few. While the behaviour of low dimensional systems, where usually only one or two attractors rule the dynamics, has been studied intensively, much less is known for systems where several attractors coexist for a given set of the system parameters. This feature, called multistability is commonly found in different areas of science such as physics, chemistry, biology, economy, and in nature [16].

In this paper, we reinvestigate in this context the paradigmatic model of nonequilibrium statistical physics, namely, underdamped Brownian motion in a biased periodic potential. This nonlinear system enjoys never ending interest as its different aspects have already been studied for several decades [17–32]. The latter are mostly focused on the diffusive properties of the system. For instance, it may exhibit unusual phenomena such as the giant diffusion [19,20,25,32] or the non-monotonic temperature dependence of a diffusion coefficient [25,30,32,33]. Both these effects are related to a bistability observed in the velocity dynamics of the system. The later effect is well known due to the work by Risken et al. [34] who found that at low friction and appropriate bias values the velocity can be stable in a locked solution (the particle is trapped in a potential minimum) but also in a running solution (the motion is unbounded in space).

Here, we focus on multistability of the Brownian velocity dynamics in a tilted periodic potential. Despite so many years of intensive research on various aspects of this setup, the latter peculiar effect has been addressed only very recently [30] and later it was explained by

**Citation:** Spiechowicz, J.; Hänggi, P.; Łuczka, J. Velocity Multistability vs. Ergodicity Breaking in a Biased Periodic Potential. *Entropy* **2022**, *24*, 98. https://doi.org/10.3390/ e24010098

Academic Editor: Antonio Maria Scarfone

Received: 16 December 2021 Accepted: 5 January 2022 Published: 7 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

recoursing to the arcsine law [35], which is a cornerstone of extreme-value statistics. Specifically, we investigate the role of ergodicity breaking and its consequences on the velocity multistability. Ergodicity lies at the basis of statistical mechanics and implies that, over long enough observation times, the time averages of observables correspond to the equilibrium ensemble averages [1,36,37]. Equivalently, it states that a single trajectory is representative for the ensemble. An increasing number of systems exhibit nonergodic properties [1,36,37], in particular due to the ultra slow dynamics and non-exponential relaxation.

The paper is organized as follows. In Section 2, we recall the formulation of the model and introduce the dimensionless quantities. In Section 3, we discuss the results, in particular the effect of ergodicity breaking on the velocity multistability occurring in this paradigmatic system. Finally, Section 4 provides a discussion and concluding remarks.

## **2. Methods**

In this work, we study dynamics of a classical inertial Brownian particle of mass *M* moving in a spatially periodic and symmetric potential *U*(*x*) = *U*(*x* + *L*) of period *L* and subjected to a static bias *F*. This system can be described by the following Langevin equation

$$M\ddot{\mathbf{x}} + \Gamma \dot{\mathbf{x}} = -l\mathcal{U}(\mathbf{x}) + F + \sqrt{2\Gamma k\_B T} \mathcal{J}(t),\tag{1}$$

where the dot and prime denote differentiation with respect to the time *t* and the particle coordinate *x*, respectively. The parameter Γ is the friction coefficient, *T* is temperature, and *kB* denotes the Boltzmann constant. We consider the potential *U*(*x*) in the form

$$
\Delta I(\mathbf{x}) = -\Delta lL \sin\left(\frac{2\pi}{L}\mathbf{x}\right),
\tag{2}
$$

where Δ*U* denotes half of the potential barrier height. Thermal equilibrium fluctuations are modeled by the *δ*-correlated Gaussian white noise whose statistical characteristics read

$$
\langle \mathfrak{J}(t) \rangle = 0, \quad \langle \mathfrak{J}(t) \mathfrak{J}(s) \rangle = \delta(t - s). \tag{3}
$$

The noise prefactor 2Γ*kBT* satisfies the fluctuation–dissipation theorem that ensures the canonical Gibbs statistics when the system is at the equilibrium state.

The above Langevin Equation (1) can be transformed into the dimensionless form

$$
\ddot{\mathfrak{X}} + \gamma \dot{\mathfrak{X}} = -\mathcal{U}'(\mathfrak{X}) + \sqrt{2\gamma\theta} \,\, \hat{\mathfrak{Y}}(\hat{\mathfrak{Y}}) \tag{4}
$$

by introducing the rescaled coordinate *x*ˆ and time ˆ*t*,

$$\mathbf{f} = \frac{2\pi}{L}\mathbf{x}, \quad \mathbf{f} = \frac{\mathbf{f}}{\tau\_0}, \quad \tau\_0 = \frac{L}{2\pi} \sqrt{\frac{M}{\Delta U}} \tag{5}$$

where the characteristic time *τ*<sup>0</sup> is proportional to the inverse of frequency *ω*<sup>0</sup> of small oscillations in the potential well of *U*(*x*). The effective dimensionless potential is

$$\mathcal{U}(\mathfrak{k}) = -\sin \mathfrak{k} - f \mathfrak{k}.\tag{6}$$

The dimensionless friction coefficient *γ* and bias *f* read

$$\gamma = \frac{1}{2\pi} \frac{L}{\sqrt{M\Delta l I}} \Gamma, \quad f = \frac{1}{2\pi} \frac{L}{\Delta l I} F. \tag{7}$$

The rescaled temperature *θ* is the ratio of thermal energy *kBT* to half of the barrier height the particle needs to overcome the original potential well, namely,

$$
\theta = \frac{k\_B T}{\Delta U}.\tag{8}
$$

The dimensionless thermal noise ˆ *ξ*(ˆ*t*) is statistically equivalent to *ξ*(*t*), meaning that it is a stationary Gaussian stochastic process with vanishing mean. Later, we use only the rescaled quantities, and therefore in order to improve the readability of notation from now on we omit the hat appearing in Equation (4).

The model of a Brownian particle moving in a washboard potential formulated in terms of the Langevin Equation (4) served for decades as a tool for the investigation of transport effects occurring in both classical and quantum systems. For instance, it has been employed for understanding the dynamics of phase across the Josephson junction [38], rotating dipoles in external fields [39], superionic conductors [40], charge density waves [41], and cold atoms dwelling in optical lattices [42–44]. Further systems are mentioned in Ref. [17].

The analytical methods of solution for the Fokker–Planck equation corresponding to Equation (4) are not yet elaborated, therefore in doing so, we rely solely on precise numerical simulations. All calculations have been completed using a Compute Unified Device Architecture (CUDA) environment implemented on a modern desktop graphics processing unit (GPU). This method allowed to speed up necessary calculations by a factor of the order 103 as compared to the traditional methods. We refer interested reader to Refs. [45,46] where more details on this scheme can be found. Here, we only mention that all quantities of interest were averaged over the ensemble of 219 = 524,288 system trajectories.

## **3. Results**

In this paper, we investigate various aspects of multistability in the velocity dynamics of a Brownian particle dwelling in a tilted periodic potential. This interesting phenomenon has been addressed for this setup only very recently [35], however, it has been also reported in systems driven by other types of noise. Examples include fractional Gaussian noise [29], Ornstein–Uhlenbeck, and harmonic Levy noise [47,48], to name but a few.

In Figure 1 we exemplify the velocity multistability phenomenon occurring in this system. The probability distribution *p*(*v*) for the instantaneous long time velocity *v* obtained from the histogram of the latter quantity is depicted for fixed time *t* = 104 (but is time-invariant) and for different dimensionless temperatures *θ*. The issue of measurement of the instantaneous velocity of a Brownian particle is presented in Ref. [49]. In Figure 1, one can observe three well pronounced maxima. One of them corresponds to the velocity *v* = 0 (the locked state) and the other two with *v* = 0 are related to running solutions. This means that these values occur significantly more frequently than the others and therefore are more stable. This observation matches the common definition of multistability for stochastic systems [50]. We note that, as temperature *θ* increases, the difference between each maximum becomes less pronounced and eventually disappears.

**Figure 1.** The probability distribution *p*(*v*) for the instantaneous long time velocity *v* = *v*(*t*) of the Brownian particle is illustrated for *t* = 10<sup>4</sup> and selected values of temperatures *θ* of the system. The used parameters read *γ* = 0.66 and *f* = 0.91.

In Ref. [35], the origin of the multistability effect is explained in terms of the arcsine law for the velocity dynamics at the zero temperature limit *θ* = 0, i.e., as the trace of deterministic dynamics perturbed by thermal noise. In such a case in the long time regime, the velocity *v*(*t*) of the particle is a time-periodic function. Moreover, the ergodicity of the setup is strongly broken, which means that its phase space can be divided into two non-intersecting invariant sets corresponding to the locked and running state [51]. We visualize this in Figure 2, where the time averaged particle velocity

$$\mathbf{v} = \lim\_{t \to \infty} \frac{1}{t} \int\_0^t ds \, \dot{\mathbf{x}}(s) \tag{9}$$

is depicted as a function of the initial conditions for the coordinate *x*(0) = *x*<sup>0</sup> and velocity *v*(0) = *v*0. The black region corresponds to the locked state with **v** = 0 whereas the grey one indicates the regime of a running solution for which **v** = 0. Therefore, different initial conditions {*x*0, *v*0} can lead to a distinct average velocity **v**. It is a disturbing situation, as typically in experiments the initial conditions are not known a priori or can be settled only with a finite resolution. To get rid of the dependence of the obtained results on the initial conditions, one needs to average over them. In Ref. [35], the authors distributed *x*<sup>0</sup> and *v*<sup>0</sup> uniformly over the intervals [0, 2*π*] and [−2, 2], respectively. Moreover, they found that in such a case the initial conditions induce an almost uniformly distributed phase shift *ϕ* in the time-periodic dependence of the velocity *v*(*t*) in the long time regime. This in turn results in the arcsine law for the velocity probability density *p*(*v*) which constitutes the backbone of multistability in this system.

**Figure 2.** The basins of attraction for the time averaged velocity **v** of the particle. The black colour codes the locked state **v** = 0 whereas the grey part indicates the regime with running solutions **v** = 0. Parameters read *γ* = 0.66, *f* = 0.91 and *θ* = 0.

In this work, we present a complementary study. Namely, we investigate in detail the influence of various distributions of initial conditions {*x*0, *v*0} on the velocity multistability phenomenon. In Figure 3, we show the probability distribution *p*(*v*) for the instantaneous long time velocity *v* of the Brownian particle for the deterministic system *θ* = 0 and different choice of the initial conditions. In simulations, the moment of time is fixed *ti* = 104. In the inset we depict the corresponding probability distribution *P*(**v**) for the time averaged velocity **v**. In panel (a) the initial position and velocity are fixed, *x*<sup>0</sup> = 0, *v*<sup>0</sup> = 0. The corresponding probability densities are represented by the Dirac-delta *px*<sup>0</sup> (*x*) = *δ*(*x*) and *pv*<sup>0</sup> (*v*) = *δ*(*v*), respectively. Consequently, as the system is noiseless *θ* = 0, the resulting probability distributions *p*(*v*) and *P*(**v**) for the instantaneous long time *v* and time averaged velocity **v**, respectively, take the Dirac-delta forms. All phase space trajectories follow the same route, and the multistability effect is absent. The situation changes drastically already if the initial position of the particle is distributed uniformly over the period *L* = 2*π* of the potential *U*(*x*), i.e., *px*<sup>0</sup> (*x*) = U(0, 2*π*) (see panel (b)). Here, U(*a*, *b*) indicates the uniform distribution over the interval [*a*, *b*]. The starting velocity of the particle can be fixed *pv*<sup>0</sup> (*v*) = *δ*(*v*) but still the systems display multimodality in the probability density *p*(*v*). In fact, in such a case, even four distinct maxima are visible there. In the inset, we note that both locked **v** = 0 and running **v** = 0 states are represented in the ensemble of system trajectories. If we permute the initial conditions, i.e., the starting coordinate *px*<sup>0</sup> (*x*) = *δ*(*x*) but *pv*<sup>0</sup> (*v*) = U(−2, 2) (see panel (c)), the multistability emerges but the locked state is not sampled at all. This situation can be modified depending on the choice of the initial coordinate as it is demonstrated in panel (d) where, in contrast, *px*<sup>0</sup> (*x*) = *δ*(*x* − *π*).

**Figure 3.** The probability distribution *p*(*v*) for the instantaneous long time velocity *v* of the Brownian particle is depicted in the deterministic regime *θ* = 0 for *t* = 10<sup>4</sup> and different choice of the initial conditions for the system. Panel (**a**): *px*<sup>0</sup> (*x*) = *δ*(*x*), *pv*<sup>0</sup> (*v*) = *δ*(*v*); (**b**): *px*<sup>0</sup> (*x*) = U(0, 2*π*), *pv*<sup>0</sup> (*v*) = *δ*(*v*); (**c**): *px*<sup>0</sup> (*x*) = *δ*(*x*), *pv*<sup>0</sup> (*v*) = U(−2, 2); (**d**): *px*<sup>0</sup> (*x*) = *δ*(*x* − *π*), *pv*<sup>0</sup> (*v*) = U(−2, 2); (**e**): *px*<sup>0</sup> (*x*) = U(0, 2*π*), *pv*<sup>0</sup> (*v*) = U(−2, 2); and (**f**): *px*<sup>0</sup> (*x*) = N(0, 1), *pv*<sup>0</sup> (*v*) = N(0, 1), where U(*a*, *b*) indicates the uniform distribution over the interval [*a*, *b*]. Likewise, N(*μ*, *σ*2) is the Gaussian distribution with the mean *μ* and the variance *σ*2. In the inset, the corresponding probability distribution *P*(**v**) for the time averaged velocity **v** is shown. Parameters read *γ* = 0.66, *f* = 0.91, and *θ* = 0.

Overall, if the ergodicity of the system is broken, the initial conditions are never forgotten and therefore crucially impact the results. Depending on the circumstances, this behaviour may be seen as a feature, not a bug. Nevertheless, the only way to cure it is to properly average over the initial conditions. In doing so, each of them must be taken into

account equally; none can be preferred. This requirement translates to the fact that initial conditions must be equally probable and therefore uniformly distributed over the whole phase space. As the system under consideration is spatially periodic *U*(*x*) = *U*(*x* + *L*), the periodic boundary condition can be employed to yield *px*<sup>0</sup> (*x*) = U(0, 2*π*). It is not the case for the starting velocity *v*<sup>0</sup> of the particle which in principle is unbounded. However, naturally such a situation cannot be implemented in numerical simulations, and therefore one needs to carefully check the impact of the initial velocity subspace volume on the obtained results. As we demonstrated, if this is not performed thoroughly one can significantly spoil the outcomes and, e.g., break the inherent symmetries of the system [52]. We checked that, in the considered regime, the condition *pv*<sup>0</sup> (*v*) = U(−2, 2) is sufficient and that further increase in the initially chosen velocity subspace volume would not alter the outcomes. In Figure 3e, we reproduce the result from Ref. [35] obtained for *px*<sup>0</sup> (*x*) = U(0, 2*π*) and *pv*<sup>0</sup> (*v*) = U(−2, 2). The characteristic U-shape part which portrays the arcsine law corresponding to the running state is visible in the probability density *p*(*v*). Consequently, the velocity dynamics is multistable.

One can claim that the initial conditions, especially the velocity, should be distributed according to the Gaussian probability density, as then it obeys the canonical Gibbs statistics (Maxwell–Boltzmann distribution) valid for equilibrium systems. Obviously, such a choice does not satisfy the above discussed condition of equal probability. In panel (f), we show that, as a consequence of the non-uniformity for *px*<sup>0</sup> (*x*) = N(0, 1) and *pv*<sup>0</sup> (*v*) = N(0, 1), where N(*μ*, *σ*2) is the Gaussian distribution with the mean *μ* and the variance *σ*2, the results are deformed and the arcsine law is not properly recovered. There is one more argument that the condition of equal probability is the only one to be correct and consistent with the case of non-zero temperature. In the running state, the long time velocity trajectory *v*(*t*) is a periodic function of time and can be well approximated by the simple periodic function [35]

$$V(t) = A\sin\left(\omega t + \phi\right) + c.\tag{10}$$

For a fixed set of the system parameters, the constants (*A*, *ω*, *c*) are the same for all initial conditions {*x*0, *v*0}. However, the distribution of the phase shift *φ* depends on the distribution of initial conditions {*x*0, *v*0}. This fact is reflected in different probability densities *p*(*v*) for the instantaneous velocity depicted in Figure 3. As the ergodicity of the system is broken, the distributions *p*(*v*) generally depend on the measurement time *t* = *ti*. The exception is the uniform distribution for the phase *φ* corresponding to the panel (e) in Figure 3 for which *p*(*v*) is time-invariant [35]. The latter feature is characteristic for ergodic systems and is crucial from the experimental point of view.

As we just reported, the ergodicity of the deterministic system with *θ* = 0 is broken for the parameter regimes in which it exhibits the multistability phenomenon [32]. One may argue that the case *θ* = 0 is only an idealization and, in practice, there exists no realistic situation with zero temperature. However, the ergodicity breaking in a deterministic system often also carries prominent consequences for non-zero temperature. In particular, for any positive temperature *θ* > 0 the system described by Equation (4) is always ergodic, although it is not a trivial fact, as it is driven by noise [53]. At non-zero temperatures, the whole phase space is accessible due to thermally activated escape events connecting the coexisting deterministic disjoint attractors. However, if the temperature tends to zero *θ* → 0, the time *τ* after it is fully sampled becomes extremely long and goes to infinity as *τ* → ∞. From an experimental point of view, due to finite observation time, the system seemingly behaves as being non-ergodic although in fact it is ergodic. Such a situation is often termed as weak ergodicity breaking [36,37,51] and can be identified with an unusually slow relaxation of the system towards its steady state which manifests itself as the nonequivalence of time and ensemble averages. In the latter case, the initial conditions do not fade, but in fact modify the results. We exemplify this feature in Figure 4a, where we depict the probability distribution *p*(*v*) for the instantaneous long time velocity *v* of the Brownian particle for different initial conditions and low temperature *θ* = 0.0001. Clearly, when the particle starts from *x*<sup>0</sup> = 0 and *v*<sup>0</sup> = 0 (see the red solid line), even in the long time limit, there are only running solutions. On the other hand, if the initial position *x*<sup>0</sup> and velocity *v*<sup>0</sup> are either uniformly or normal distributed (see the blue or green line, respectively), the multistability emerges, but one can still note quantitative difference between these two initial conditions. In contrast, in panel (b) we depict the same characteristics but for higher temperature *θ* = 0.05. Then, thermal fluctuations are strong enough to recover the ergodicity of the system and there are no longer differences between different initial conditions. Even when the particle trajectories start from the same point in the phase space *x*<sup>0</sup> = 0, *v*<sup>0</sup> = 0 (see the red solid line), the whole density is obtained.

**Figure 4.** The probability distribution *p*(*v*) for the instantaneous long time velocity *v* of the Brownian particle is depicted for *t* = 10<sup>4</sup> and different initial conditions of the system. The red solid line indicates *px*<sup>0</sup> (*x*) = *δ*(*x*), *pv*<sup>0</sup> (*v*) = *δ*(*v*). The blue dotted line corresponds to *px*<sup>0</sup> (*x*) = U(0, 2*π*), *pv*<sup>0</sup> (*v*) = U(−2, 2). The green dashed line denotes *px*<sup>0</sup> (*x*) = N(0, 1), *pv*<sup>0</sup> (*v*) = N(0, 1). In panel (**a**) temperature is *θ* = 0.0001 while in (**b**) *θ* = 0.05. Other parameters read *γ* = 0.66, *f* = 0.91.
