Bimodality and hysteresis in systems driven by confined L\'evy flights

We demonstrate occurrence of bimodality and dynamical hysteresis in a system describing an overdamped quartic oscillator perturbed by additive white and asymmetric L\'evy noise. Investigated estimators of the stationary probability density profiles display not only a turnover from unimodal to bimodal character but also a change in a relative stability of stationary states that depends on the asymmetry parameter of the underlying noise term. When varying the asymmetry parameter cyclically, the system exhibits a hysteresis in the occupation of a chosen stationary state.


INTRODUCTION
The Langevin description of an overdamped Brownian motion in a potential V (x) constitutes a basic paradigm to study the effects of fluctuations at the mesoscopic scale [1,2,3]. Here prime stands for differentiation over x and ζ(t) is a commonly assumed white Gaussian noise representing close-to-equilibrium fluctuations of the dynamic variable x. In contrast, in farfrom-equilibrium situations, the Gaussianity of the noise term may be questionable, due to e.g. strong interaction with the surrounding "bath" [4,5]. A natural generalizations to the Brownian motion are then δ-correlated Lévy stable processes which can be interpreted as fluctuations resulting from strong collisions between the test particle and the environment. In particular, the scale-free, self similar feature of Lévy distributions gives rise to the occurrences of large increments of the position coordinates ∆x during small time intervals causing a non-local character of the motion. Within the paper we address the problem of kinetics as described by Eq. (1) under the action of ζ(t) representing a stationary white Lévy noise [6,7]. Accordingly, the position of the Brownian particle subjected to additive white Lévy noise is calculated by direct integration of Eq. (1) with respect to the α-stable measure [6,7,8,9,10,11] L α,β (s), i.e., where ζ is distributed according to the α-stable Lévy type distribution L α,β (ζ; σ, µ = 0) whose representation [6,8,12] is given by the characteristic function φ(k) defined in the Fourier and for α = 1 The stability index α, determining tails of the probability density function (PDF) takes values α ∈ (0, 2], the skewness of the distribution is modeled by the asymmetry parameter β ∈ [−1, 1]. Indices α and β classify the type of stable distributions up to translations and dilations [6,8]. Two other parameters of scaling σ ∈ (0, ∞) and location µ ∈ (−∞, ∞) can vary, although replacing ζ − µ and σζ in the original coordinates, shifts the origin and rescales the abscissa without altering function L α,β (ζ). For simplicity, we will restrict here to strictly stable Lévy noises with µ = 0 [6,8]. Generally, for β = µ = 0 PDFs are symmetric while for β = ±1 and α ∈ (0, 1) they are totally skewed, i.e., ζ is always positive or negative only, depending on the sign of asymmetry parameter β (cf. Fig. 2). Asymptotically, for ζ → ∞ with α < 2, stable PDFs behave as p(ζ) ∝ |ζ| −(α+1) causing divergence of moments ζ ν = ∞ for ν > α. The asymmetry is reflected in a biased distribution lim ζ→∞ P (Z>ζ) P (|Z|>ζ) = 1+β 2 . Equivalent to the stochastic ordinary differential Eq. (1) is a fractional Fokker-Planck equation [13,14,15,16,17] (FFPE) for the probability distribution function with p(x, t) representing the probability density functions for finding a particle at time t in the vicinity of x, φ(k, t) = exp[ikx(t)] = exp ik t 0 ζ(s)ds standing for the characteristic function of the stable process and −∆ α f (x) denoting fractional Laplacian [18,19] −(∆) α/2 f (x) = F −1 |k| αf (k) , with α = 2 corresponding to the standard Brownian diffusion case. The addition of the potential force −V ′ (x) to Eq. (1) adds the classical drift term ∂ ∂x [V ′ (x)p(x, t)] to Eq. (4). In the approach presented herein, instead of solving Eq. (4), information on the system is drawn from the statistics of numerically [6,8,20] generated trajectories satisfying the generalized Langevin equation (1). At a single trajectory level sampled from the stochastic dynamic study of the problem, the initial condition for Eq. (1) has been set to x(0) d = U[−1, 1], i.e., initial position of the particle is drawn from the uniform distribution over the interval [−1, 1]. For simplicity, Eq. (1) has been studied in dimensionless variables [21,22] with additionally setting σ = 1. The time independent potential V (x) is assumed to be of the form V (x) = x 4 /4 which guarantees the confinement of the trajectory x(t) within the potential well [21,23,24] leading to a finite variance of the stationary PDF. For the general α-stable driving and the quartic potential the stationary PDF fulfills whereP (k) stands for a Fourier transform of the stationary PDF, P (x) = lim t→∞ P (x, t). Analytical solutions of Eq. (5) can be readily obtained for a Gaussian case α = 2: P (x) ∝ exp(−x 4 /4) and for a Cauchy additive noise α = 1: [21,23,24]. They display a perfect agreement with profiles of PDFs obtained by numerical simulations of Eq. (1) performed with the Janicki-Weron algorithm [6,8,25,26].

RESULTS
Numerical results were constructed for a time step of integration ∆t = 10 −3 , simulation length T max = 10 with an overall statistics of N = 5 × 10 4 realizations. To check whether results are influenced by the length of simulations, results for various T max (T max = 10, T max = 15) were compared showing consistency of estimated PDFs for both values. Further details on values of parameters are included in the text underlying the figures. Numerical examination of Eq. (1) allows for construction of PDF estimators for the whole range of parameters α and β. Likewise, by direct integration of Eq. (1) it is also possible to investigate time-dependent PDFs [20] and noise-induced bimodality of the probability distribution [20,21,23,24].
The change of α (for β = 0) from values greater than 1 to values smaller than 1 results in change of location of a modal value of stable densities, see left vs. right panel of Fig. 2. The shift of modal values is reflected in properties of stationary states. For α > 1 with β < 0 the modal value is located for x > 0 (right panel of Fig. 3) while for α < 1 with β < 0 it shifts to the negative interval x < 0 (left panel of Fig. 3). For β < 0 modal values are located on the opposite side of the origin than for β > 0. Therefore, the position of the modal value can be moved from the positive to negative real lines by change of α while β is kept constant (left vs. right panel of Fig. 3) or by change of β to −β with a preset value of α (cf. different rows in Fig. 3).
The fact that changes in β can change the location of the modal value suggests that periodic changes in this parameter can lead to a phenomenon resembling dynamical hysteresis [27,28,29]. In order to register a hysteretic behavior of the system, we have performed an analysis of trajectories of Eq. (1) based on a two-state approximation. For that purpose, we have defined an occupation probability of being in left/right state according to Transition between the states is induced by time dependent asymmetry parameter β which is periodically modulated over time, i.e., β = β max cos Ωt = β max cos 2π TΩ . In the Fig. 4 values of p(left) for various T Ω (with β max = 1) are presented. Stable random variables for α < 1 with |β| = 1 take only positive/negative values, depending on the sign of β. Therefore higher level of saturation is observed for α = 0.9 (left panel) than for α = 1.1 (right panel), i.e., when |β| = 1 probability of being in the left/right state for α = 0.9 is higher than for α = 1.1 (right panel). The direction of the hysteresis loop is a direct consequence of the fact that changes in β move modal values from positive/negative real line to the negative/positive real line. Due to the initial condition imposed on x(0), the starting point for each hysteresis loop is (0,0.5) and the first part of the curve describes approaching to the proper hysteresis loop. With decreasing Ω (increasing T Ω ) area of hysteresis loop decreases because the system has more time for relaxation and response to the changes of β. For large Ω (small T Ω ) the response of the system is more delayed and consequently area of the hysteresis loop is larger.
Finally the influence of β max on the shape of hysteresis loop has been examined. In Fig. 5 hysteresis loops for various β max (with T Ω = 9) are presented. Smaller β max makes stable distribution less skewed and as a consequence, less probability mass becomes located in the left state and the effect of saturation is less visible. Furthermore, with decreasing β max , loops become more oval and finally for β max = 0 the hysteresis phenomenon disappears.

SUMMARY AND CONCLUSIONS
The modulation of stable noise parameters modify shape of stationary densities corresponding to Eq. (1). In particular skewed noise characterized by nonzero asymmetry parameter β can induce asymmetry of stationary states in symmetric potentials. Furthermore, totally skewed stable noises (α < 1 with |β| = 1) could move the probability mass to one side of the x-axis making stationary states totally skewed. In such situations, the whole probability mass can be located on the left hand side or on the right hand side of the origin x = 0 depending on the sign of the asymmetry parameter β. For α > 1 with a nonzero β, stable noises still are asymmetric. Nevertheless, even |β| = 1 is not sufficient to induce totally skewed stationary states. Consequently, dynamical hysteresis loops induced by cyclic variation of β for α < 1 and α > 1 are characterized by various level of saturation, see left vs. right panel of Fig. 4.
The stationary densities depicted in Fig. 3 are recorded in the system described by Eq. (1) which is subjected to the action of stable noises. The very same stationary densities can be also observed in the equilibrium system perturbed by the white Gaussian noise, i.e., . The dynamical hysteresis detected in the model described by Eq. (1) emerges as a consequence of periodical modulation of the asymmetry parameter β. In the effective potential model, periodical modulation of the asymmetry parameter corresponds to the periodical modulation [30] of the effective potential V eff (x). Nevertheless both models are not fully equivalent. They have the same one-point densities p(x, t), while other characteristic of the process {x(t)} are distinct [22]. Therefore, examination of stationary densities itself is not a fully conclusive method for discrimination of underlying types of noises.
The dynamical hysteresis can be also observed in the generic double well potential model subject to the joint action of the deterministic periodic modulation and stochastic α-stable fluctuations [20]. In such a case, however, the shape of the dynamical hysteresis loop is affected both by the character of the noise and by the type of periodic modulation [20].
The research has been supported by the Marie Curie TOK COCOS grant (6th EU Framework Program under Contract No. MTKD-CT-2004-517186). Computer simulations have been performed at the Academic Computer Center CYFRONET AGH, Kraków. Additionally, BD acknowledges the support from the Foundation for Polish Science and the hospitality of the Humboldt University of Berlin and the Niels Bohr Institute (Copenhagen).