Multimodal stationary states under Cauchy noise

TL;DR

This paper investigates how Le9vy noise influences stationary states in double-well potentials, revealing conditions for bimodality and its impact on noise-induced escape phenomena.

Contribution

It demonstrates that high potential barriers can induce bimodal stationary states in each well under Le9vy noise, and analyzes how barrier height affects multimodality.

Findings

High barriers lead to bimodal stationary states.

Decreasing barrier height reduces multimodality.

Multimodality influences noise-induced escape dynamics.

Abstract

A L\'evy noise is an efficient description of out-of-equilibrium systems. The presence of L\'evy flights results in a plenitude of noise-induced phenomena. Among others, L\'evy flights can produce stationary states with more than one modal value in single-well potentials. Here, we explore stationary states in special double-well potentials demonstrating that a sufficiently high potential barrier separating potential wells can produce bimodal stationary states in each potential well. Furthermore, we explore how the decrease in the barrier height affects the multimodality of stationary states. Finally, we explore a role of the multimodality of stationary states on the noise induced escape over the static potential barrier.

Figures

Tables

Table 1. Table 1: The ratio of probabilities of escape in a single jump and two consecutive jumps, i.e. p ​ ( 1 ) / p ​ ( 2 ) 𝑝 1 𝑝 2 p(1)/p(2) .

$h$	$n=2$	$n=4$
$0.1$	0.75/0.14	0.73/0.15
$1$	0.82/0.12	0.79/0.13
$5$	0.86/0.09	0.80/0.11
$10$	0.87/0.09	0.81/0.11

Equations

$$\frac{dx}{dt}=-V^{\prime }(x)+\sigma {\zeta }_{\alpha }(t),$$

$$\Delta x$$

$$\varphi (k)=\exp \left[-|k{|}^{\alpha }\right].$$

$$\varphi (k)=\exp \left[-{\sigma }^{\alpha }|k{|}^{\alpha }\right],$$

$$f_2(x)=\frac{1}{\sqrt{4\pi {\sigma }^2}}\exp \left[-\frac{x^2}{4{\sigma }^2}\right].$$

$$f_1(x)=\frac{\sigma }{\pi (x^2+{\sigma }^2)},$$

$$\frac{\partial p(x,t)}{\partial t}=-\frac{\partial }{\partial x}{V}^{\prime }(x,t)p(x,t)+{\sigma }^{\alpha }\frac{{\partial }^{\alpha }p(x,t)}{\partial |x{|}^{\alpha }}.$$

$$V(x)=n{x}^n$$

$$n>2-\alpha .$$

$$p_{\alpha =1}(x)=\frac{4}{\pi (4\sigma {)}^{\text{nicefrac}13}}\times \frac{1}{{\left[\frac{4x}{(4\sigma {)}^{\text{nicefrac}13}}\right]}^4-{\left[\frac{4x}{(4\sigma {)}^{\text{nicefrac}13}}\right]}^2+1}.$$

$$p_{\infty }(x)=\frac{\Gamma (\alpha )(2L{)}^{1-\alpha }(L^2-x^2{)}^{\alpha /2-1}}{{\Gamma }^2(\alpha /2)}.$$

$$V(x)=n(|x|-1)^{n}=\left\{\begin{array}[]{ccl}n(x+1)^{n}&\mbox{for}&x<0\\ n(x-1)^{n}&\mbox{for}&x\geqslant 0\\ \end{array}\right.,$$

$$p(x)=\left\{\begin{array}[]{ccl}p_{\alpha=1}(x+1)/2&\mbox{for}&x<0\\ p_{\alpha=1}(x-1)/2&\mbox{for}&x>0\\ \end{array}\right.,$$

$$V(x)=(|x|-h^{1/n}{)}^n,$$

$$p(x)=\frac{2}{\pi }\times \frac{1}{(2x{)}^2+1}.$$

$$p(x)=\frac{2}{\pi }\times \frac{1}{(2x{)}^4-(2x{)}^2+1}.$$

$$V(x)=h(|x|-1{)}^n,$$

$$\tau =\min \{\tau :x(0)=1\text{mbox}andx(\tau )⩽0\}.$$

$$\frac{dx}{dt}=-x^3+{\zeta }_{\alpha }(t).$$

$$f_{\alpha =1}(x)=\frac{1}{\pi (x^4-x^2+1)}.$$

$$\frac{dx}{dt}=-x^3+\sigma {\zeta }_1(t)$$

$$\tilde{x}=\frac{x}{x_0}$$

$$\tilde{t}=\frac{t}{t_0}.$$

$$\sigma {\zeta }_{\alpha }(t)$$

$$\frac{x_0}{t_0}\frac{d\tilde{x}}{d\tilde{t}}=-x_0^3{\tilde{x}}^3+\sigma t_0^{\frac{1}{\alpha }-1}{\zeta }_{\alpha }(\tilde{t}).$$

$$\left\{\begin{array}[]{lcl}t_{0}x_{0}^{2}&=&1\\ \sigma\frac{t_{0}^{\frac{1}{\alpha}}}{x_{0}}&=&1\\ \end{array}.\right.$$

$$\left\{\begin{array}[]{lcl}x_{0}&=&\sigma^{\frac{1}{3}}\\ t_{0}&=&\sigma^{-\frac{2}{3}}\\ \end{array}.\right.$$

$$f_{\alpha =1}(x)$$

$$\frac{dx}{dt}$$

$$a\frac{dy}{dt}=-\lambda a^{n-1}{n}^2{y}^{n-1}+\sigma {\zeta }_{\alpha }$$

Full text

I Introduction

Noise is one of the fundamental concepts in the theory of stochastic systems Horsthemke and Lefever (1984). It is used to effectively approximate complex interactions of an observed particle with its environment. After a long lasting dominance of a Gaussian paradigm, it has been documented and accepted that many situations significantly departs from the Gaussian approximation. In the non-equilibrium regime, noise can be of the $\alpha$-stable type Dubkov et al. (2008) resulting in unexpected properties of noise-driven systems.

In the recent two decades, significant progress in the development of the theory of systems driven by the Lévy noise has been achieved Metzler and Klafter (2000); Barkai (2001); Anh and McVinish (2003); Brockmann and Sokolov (2002); Chechkin et al. (2006); Jespersen et al. (1999); Yanovsky et al. (2000); Schertzer et al. (2001). Lévy flights have been studied in various situations Klages et al. (2008). Their applications cover, among others, economy and finance Bouchaud and Georges (1990), superdiffusion of micellar systems Bouchaud et al. (1991), studies of turbulence Shlesinger et al. (1986), description of photons in hot atomic vapours Mercadier et al. (2009), laser cooling Bardou et al. (2002); Barkai et al. (2014) and therapeutic aspects Cabrera and Milton (2004). These studies cover experimental Solomon et al. (1993); Barthelemy et al. (2008) and various theoretical aspects Eliazar and Klafter (2003); Sokolov and Metzler (2004); Garbaczewski and Stephanovich (2009); Garbaczewski (2010) of Lévy flights.

In the context of current research, the problem of stationary states in systems driven by the Lévy noises is especially important. Stationary states in a single-well potential can be bimodal what is well known and documented Chechkin et al. (2002, 2003a, 2004, 2006, 2008). Moreover, conditions for a steepness of single-well potentials that can bound Lévy flights have been developed Dybiec et al. (2010). Furthermore, in Capała and Dybiec (2019) we have presented sample, fine tailored, single-well potentials resulting in stationary states with an arbitrary number of modal values. In Capała and Dybiec (2019) and Chechkin et al. (2002, 2004) the number of modal values has been attributed to the number of maxima of the potential curvature.

The multimodality of stationary states can be also dynamically induced Dybiec and Schimansky-Geier (2007); Calisto et al. (2017). For example in single-well potentials perturbed by the Gaussian white noise and the Markovian dichotomous noise the stationary state can be bimodal. For such a system the maxima of stationary densities are placed in the vicinity of minima of the altered potential Dybiec and Schimansky-Geier (2007) indicating differences between mechanisms producing multimodal stationary states by solely action of the Lévy noise and the combined action of the Gaussian white noise and the Markovian dichotomous noise.

Within the current manuscript we further inspect the problem of multimodal stationary states in systems driven by the Lévy noises, however, we focus on double-well potentials. Therefore, we explore a situation when the number of modal values in the stationary state is larger than the number of minima of the potential. We demonstrate that, if the potential minima are deep enough, they are sufficient to produce multiple maxima of a stationary state in any of potential wells.

The studied model is presented in the next section (Sec. II). The main results are included in Secs. II.2 – II.3 and Sec. III, while Sec. II.1 plays an introductory role. The paper is closed with Summary and Conclusions (Sec. IV) and supplemented with Appendices A – C presenting relevant technicalities.

II Model and Results

We start with the basic information regarding overdamped stochastic dynamics in single-well potentials (Sec. II.1). Afterwards, in Sec. II.2, using this information, we study various types of double-well potentials which are able to produce stationary states characterized by a larger number of modal values than the number of minima of the potential. Finally, in Sec. II.3, we investigate a model with a varying barrier height which allows to increase understanding of the role played by the potential barrier separating minima of double-well potentials.

The overdamped motion in an external potential $V(x)$ is described by the following Langevin equation

[TABLE]

where $\zeta(t)$ is the symmetric $\alpha$-stable Lévy type noise, i.e. the formal time derivative of the symmetric $\alpha$-stable motion Janicki and Weron (1994a). The parameter $\sigma$, in Eq. (1), might be interpreted as a noise strength. Eq. (1) is supplemented with the initial condition $x(0)=x_{0}$. Accordingly, the stochastic process $\{X(t),t\geqslant 0\}$ governed by Eq. (1) has increments

[TABLE]

In Eq. (II) $\zeta_{t}$ represents independent, identically distributed random variables Chambers et al. (1976); Weron and Weron (1995); Weron (1996) following the symmetric $\alpha$-stable density Janicki and Weron (1994b); Janicki (1996) with the unity scale parameter and the characteristic function $\phi(k)$

[TABLE]

Within simulations we use $\sigma=1$, nevertheless, due to various types of considered potentials, the comparison with analytical results requires reintroduction of the scale parameter $\sigma$, see Appendix A and B. The non-unity scale parameter is reintroduced by the multiplication of $\zeta_{t}$ by $\sigma$. Such a multiplied (rescaled) random variable has the characteristic function

[TABLE]

which is the typical form of the characteristic function of symmetric $\alpha$-stable densities. The stability index $\alpha$ ($0<\alpha\leqslant 2$) describes an asymptotic power-law decay of $\alpha$-stable densities, which for $\alpha<2$ is of a $|x|^{-(\alpha+1)}$ type. The scale parameter $\sigma$ controls the distribution width. For $\alpha<2$, the variance of an $\alpha$-stable density is infinite, thus the distribution width can be defined by the interquantile width or fractional moments. For $\alpha=2$, the characteristic function (4) reduces to the characteristic function of the normal (Gaussian) distribution with the probability density

[TABLE]

The case of $\alpha=1$ corresponds to the Cauchy distribution

[TABLE]

which is extensively used here. For clarity and practical reasons, the scale parameter is extracted from the noise, see Eq. (1). From Appendices A and B it can be deducted how the scale parameter can be reintroduced.

The evolution of the probability density generated by Eq. (1) is described by the fractional Smoluchowski-Fokker-Planck equation Samorodnitsky and Taqqu (1994); Podlubny (1999); Yanovsky et al. (2000)

[TABLE]

The operator $\partial^{\alpha}/\partial|x|^{\alpha}$ is the fractional Riesz-Weil derivative Podlubny (1999); Samko et al. (1993) which can be defined via the Fourier transform $\mathcal{F}_{k}\left(\frac{\partial^{\alpha}f(x)}{\partial|x|^{\alpha}}\right)=-|k|^{\alpha}\mathcal{F}_{k}\left(f(x)\right).$

Results included in following subsections have been constructed numerically by methods of stochastic dynamics. Eq. (1) was integrated by the Euler-Maryuama method with the time step of integration $\Delta t=10^{-5}$ and $10^{6}$ – $10^{8}$ repetitions. More details can be found in Appendix C.

II.1 Single-well potentials

First, we study the motion of a particle driven by the Lévy noise in single-well potentials of

[TABLE]

type, where $n$ is even and greater than 0. Stationary states for potentials given by Eq. (8) exist for the sufficiently large steepness (exponent) of the potential $n$ Dybiec et al. (2010), i.e.

[TABLE]

The steepness of the potential, which is sufficient to produce stationary states, depends on the stability index $\alpha$ characterizing asymptotics of random pulses. Therefore, the existence of a stationary state does not depend solely on the potential, but it is also sensitive to the stability index $\alpha$ characterizing the noise. Furthermore, if a stationary state exists, it is not of the Boltzmann-Gibbs type Eliazar and Klafter (2003).

In Appendices A and B, starting from known formulas Chechkin et al. (2002, 2003a, 2004, 2006, 2008), it is demonstrated that for $n=4$ with $\alpha=1$ the stationary solution of Eq. (1) is

[TABLE]

The stationary density (10) depends on the scale parameter $\sigma$. Please note, that for a finite $n$, the stationary density is sensitive to the scale parameter $\sigma$ characterizing strength of noise pulses, see Dubkov and Spagnolo (2007) and Appendix B.

The infinite rectangular potential well is recovered in the $n\to\infty$ limit of Eq. (8), see Kharcheva et al. (2016). For the infinite rectangular potential well the stationary state is Denisov et al. (2008)

[TABLE]

In contrast to the finite $n$, the stationary density (11) does not depend on the scale parameter $\sigma$.

In order to verify the correctness of implemented numerical methods, see Appendix C, we have performed a computer simulation for $V(x)=4x^{4}$. Fig. 1 demonstrates the perfect agreement between results of a stochastic simulation (points) and the exact formula (dashed line) given by Eq. (10).

Further numerical studies are devoted to the investigation of changes in the stationary densities in single-well potentials of $V(x)=nx^{n}$ type induced by the increasing exponent (steepness) $n$. Analytical results for such class of potentials can be obtained by a transformation of (Dubkov and Spagnolo, 2007, Eqs. (38) and (39)), see Appendix B. In Fig. 2 these formulas are depicted by dashed lines and compared with results of stochastic simulations marked by points. Finally, the solid line presents the stationary density for the infinite rectangular potential well, see Eq. (11). From Fig. 2 it is clearly visible that results of stochastic simulations perfectly agree with analytical formulas. Moreover, analogously like in Kharcheva et al. (2016), it is demonstrated that the stationary density approaches the stationary state for the infinite rectangular potential well in the limit of $n\to\infty$, see Eq. (11) and Denisov et al. (2008). Due to symmetry of the potential (8) stationary states are also symmetric. With the increasing $n$, modal values of the stationary density are shifted towards $x=\pm 1$. Finally, in the limit of $n\to\infty$, modal values are located at impenetrable boundaries. This test is crucial for the considerations performed in Sec. II.2, where double-well potentials based on Eq. (8) are considered.

Multimodal stationary states are produced because of the interplay between the deterministic and random forces. The deterministic force $-V^{\prime}(x)$ always acts towards the minimum of the potential $V(x)$. The random force is the only factor which causes excursions towards large $|x|$. Therefore, the emerging stationary density is the distribution that balances random excursions and deterministic sliding, which is interrupted by long jumps. For the large exponent $n$, the potential $V(x)$ given by Eq. (8) is very close to the infinite rectangular potential well. A random jump to the right or to the left moves the particle to the point where the restoring force is very large. Consequently, the particle immediately slides to $|x|\approx 1$ and waits there for the next random pulse. With the decreasing exponent $n$ the point where particle can slide down moves towards $x=0$. At the same time, for $n\geqslant 2$, the time needed to reach the origin by a deterministic sliding is infinite. The impossibility of reaching $x=0$ in a deterministic way in a finite time makes the stationary state bimodal.

II.2 Double-well potentials

As it is well known, and also demonstrated in Sec. II.1, it is possible to create a bimodal stationary state in a single well potential Chechkin et al. (2002, 2003a, 2004, 2006, 2008). As the next step, it is intriguing to verify conditions necessary to produce stationary states with higher modality than the number of potential minima. Moreover, it is interesting to inspect how the shape of the stationary state depends on the height of the potential barrier separating minima of a double-well potential. In order to explore these issues, a special potential needs to be selected. Firstly, it has to behave like one of single-well potentials with $n\geqslant 2$ in the vicinity of its minima. Moreover, potential wells must be wide enough to avoid interference of structures (e.g. peaks), emerging in different potential wells. In addition, the potential barrier needs to be steep enough in order to produce internal minima of the stationary probability density. Finally, in the limit of $n\rightarrow\infty$, the height of the potential barrier should be infinite. The simplest choice fulfilling all required properties is a double-well potential based on the single-well potential given by Eq. (8). Therefore, we use the following double-well potential

[TABLE]

where $n$ is even and positive. Sample double-well potentials given by Eq. (12) with various $n$ are depicted in Fig. 3. With the increasing exponent $n$ the potential $V(x)=n(|x|-1)^{n}$ approaches the sum of two infinite rectangular potential wells. The height of the potential barrier is equal to $n$, i.e. $V(0)=n$ for every $n$. Importantly, for a finite $n$, minima of the potential (12) are located at $\pm 1$ regardless of the steepness exponent $n$.

Figure 4 shows results of a stochastic simulation for the double-well potential (12) with $n=4$. The dashed line presents

[TABLE]

where $p_{\alpha=1}(x)$ is given by Eq. (10). The exact density (points) differs from $p(x)$ given by Eq. (13), nevertheless these differences are not so pronounced. These differences originate in the finite height of the barrier separating the left and the right minimum of the double-well potential. Therefore, internal maxima of the stationary density are lower than the outer one.

Subsequently, we have explored the behaviour of stationary states for the increasing values of the exponent $n$ in Eq. (12). With the increasing $n$, positions of internal modes shift toward $x=0$. Finally, in the limit of $n\to\infty$, internal maxima of the stationary density merge, resulting in the amplification of the central mode. This phenomenon produces a stationary state with three modal values. At the same time, external modal values shift towards $x=\pm 1$ and their height increase. For a large enough steepness (exponent) $n$, the stationary state $p_{\infty}(x)$ is a sum of stationary states in two rectangular potential wells, see Eq. (11). Analogously like in Fig. 4, dashed lines in Fig. 5 present exact solutions given by transformed formulas (Dubkov and Spagnolo, 2007, Eqs. (38) and (39)), see Eq. (13) for the reference. A comparison of computer simulations (points) with dashed lines demonstrates how the asymptotic density is reached. It indicates how the interference of internal modes is responsible for significant deviations from theoretical densities. Finally, this comparison is especially useful for observing the behaviour of outer modes, since it nicely shows how, with the increasing $n$, outer maxima of probability densities approach impenetrable boundaries arising at $\pm 1$. Here, the qualitative explanation of the observed effect is the same as for the single-well potential. The only difference is that the random pulses can induce transition over the finite potential barrier separating left and right potential wells.

II.3 Barrier height

In Sec. II.2, we have verified how the shape of stationary states changes with the increasing height of the potential barrier separating potential wells. We have also explored how the limit of double-well infinite rectangular potential well is reached. Here, we would like to examine the opposite situation. Namely, we want to investigate how the shape of stationary states changes with the decreasing height of the potential barrier separating both potential wells. Since the double-well potential given by Eq. (12) is not very suitable for such tests, we use the potential

[TABLE]

with $n=2$ and $n=4$, see Fig. 6. In Eq. (14), $h$ controls the height of the potential barrier separating minima of the double-well potential. The parameter $h$ controls also a position of potential minima. In the limit of $h=0$, the double-well potential reduces to a single-well potential with a minimum at the origin. Therefore, for a decreasing $h$, it is possible to see when the barrier separating potential wells is not sufficient to produce maxima of the stationary state in left and right potential wells ($n=2$) or when bimodality within a potential well disappears ($n=4$), see Fig. 7.

The top panel of Fig. 7 presents results for $V(x)=(|x|-\sqrt{h})^{2}$. For $h$ large enough, the potential is built by two parabolas separated by a high potential barrier. Therefore, within each potential well, there is a stationary state which is given by the stationary state for the harmonic potential, see Fig. 8. With the decreasing $h$, both peaks of the stationary density are approaching each other and start to interfere. In the limit of $h=0$, the potential is of a single-well, parabolic type. Consequently, the stationary state is the same as for the stochastic (overdamped) harmonic Cauchy oscillator. From computer simulations, one may observe that for $h=10^{-4}$ the stationary state is very close to the stationary state for $V(x)=x^{2}$. A perfect agreement is obtained for $h=0$. However, due to the cusp in the potential, even for a very small $h$, $p(x=0)$ is reduced in comparison to the stationary density for the stochastic harmonic oscillator. The cusp at $x=0$ moves from the origin some of particles, which want to accumulate there. Finally, for $h=0$, the single-well parabolic potential is recovered and the unimodal stationary state is given by the appropriately rescaled Cauchy density, which for $V(x)=x^{2}$ and $\sigma=1$ is given by Eq. (B) with $\lambda=1/2$, i.e.

[TABLE]

Similar situation is observed for $n=4$, see Eq. (14) and the bottom panel of Fig. 7. For $n=4$ and large enough $h$, the potential is built by two quartic, single-well potentials. Therefore, within each potential well, the stationary state is bimodal, see Eq. (10) and Figs. 1 and 8. With the decreasing $h$, internal maxima of the stationary state decrease. Simultaneously, outer modes shift towards the origin and their height increase. Finally, for $h=0$, the single-well quartic potential $V(x)=x^{4}$ is recovered. For $\sigma=1$ in the Langevin equation, the bimodal stationary state is given by the equation derived from Eq. (36) with $n=4$ and $\lambda=1/4$, i.e.

[TABLE]

The potential given by Eq. (14) is one of many possible potentials that allow to investigate the role of a height of a barrier separating potential’s minima on the shape of stationary states. The potential (14) have been used due to its properties. In Eq. (14) the $h$ parameter controls the barrier height. On the one hand, single well potentials, e.g. harmonic ($n=2$) or quartic ($n=4$), are recovered in the limit of $h\to 0$. On the other hand, for very large $h$, there are no transitions between potential minima and the potential is practically built from two independent single-well potentials. One of the drawbacks of the potential (14) is the fact that the increase in $h$ shifts also positions of potential minima into the direction of larger $|x|$, see Fig. 6. Fig. 8 presents stationary states for large $h$: $h=10^{4}$ (quartic potential) and $h=10^{2}$ (parabolic potential). For such values of $h$, potential minima are located at $x=\pm 10$, where central parts of within the well probability densities are located. In every potential well, probability densities perfectly follow shifted and renormalized densities for single well potentials, see Eqs. (15) and (16). The perfect agreement between theoretical single-well densities and numerically estimated histograms originates in the height and width of the potential barrier, which makes the transition from one well to another very unlikely. Alternatively, one can generalize potential (12) to non-integer $n$ or use $V(x)=h(|x|-1)^{n}$. Unfortunately, these generalizations also possess some disadvantages. For the generalized Eq. (12), with the change in $n$ not only the barrier height is changed, but also the steepness (exponent) and the “slope” (prefactor) of the potential are modified. In the case of $V(x)=h(|x|-1)^{n}$, the steepness is fixed but the “slope” changes. Noticeably, in the limit of $h\to 0$, this potential does not reduce to the single-well potential, because it disappears. Nevertheless, we use this potential to inspect the role of multimodality on noise induced effects, see Sec. III.

Naturally, it is still possible to fine-tailor other, more complicated, potentials bearing required properties. Nevertheless, we have used potential (14) as a simple potential that allows to study the role of the barrier height on the shape of stationary states. Already, from the potential (14) it is possible to draw many general conclusions, see Sec. IV.

III Escape from a potential well

In Sec. II we have inspected conditions for emergence of multiple modes within a single potential well. In the current section we compare the noise induced escape from the potential well in order to explore the role of a multimodality of stationary states.

We consider a motion of the particle driven by the Cauchy noise in the double-well potential

[TABLE]

with $n\in\{2,4\}$. We are studying the escape from a potential well Kramers (1940); Dybiec et al. (2007); Chechkin et al. (2007); Ditlevsen (1999), due to symmetry of the problem, we can assume $x(0)=1$ and consider $x>0$ only. $x=0$ is the point separating left and right potential wells. Therefore, the stochastic dynamics is continued as long as $x\geqslant 0$, see Appendix C. Every first escape from the right potential well is associated with the first passage time

[TABLE]

The mean first passage time (MFPT) is the average of first passage times. In Fig. 9 the dependence of the MFPT $\langle\tau\rangle$ on the barrier height $h$ is presented. For the parabolic potential the mean first passage time is larger than for the quartic potential. With the increasing barrier height ($h>10$) the MFPTs become comparable. The MFPT is influenced also by size of the domain in which the random walk can be performed, i.e. by the outer ($|x|\gg 1$) parts of the potential. For the quartic potential, the particle is effectively bounded in a smaller fraction of space than for the parabolic potential, see top panels of Figs. 10 and 11. The effective domain size depends also on the $h$ parameter. For instance, stationary states for $V(x)=hx^{4}$ and $V(x)=hx^{2}$ can be found from Eq. (B) with $\lambda=h/4$ and from Eq. (B) with $\lambda=h/2$, respectively. The decrease of the $h$ parameter makes excursions to the large $|x|$ more probable. In turn these excursion are responsible for the increase in the MFPT. In the limit of $h\to 0$, a motion in the potential given by Eq. (17) reduces to the problem of the escape from the positive half-line which is characterized by the diverging mean first passage time Sparre Andersen (1953, 1954); Redner (2001); Chechkin et al. (2003b); Dybiec et al. (2016). Therefore, in the limit of $h\to 0$ the MFPT $\langle\tau\rangle\to\infty$. For large $|x|$, the parabolic potential is less steep, consequently, increase in the MFPT with the decrease of the barrier height for the parabolic potential is more rapid than for the quartic potential. Therefore, the role of a multimodality on escape kinetics can be assessed for a high enough potential barrier. The top panel of Figs. 10 and 11 suggest that $h=1$ is sufficient because, in the stationary state, the majority of the probability mass is localized in comparable domains.

Changes in the shape of stationary states, see top panels of Fig. 10 and 11, are also reflected in the dependence of the MFPT on the barrier height $h$, see Fig. 9. For $n=2$, the MFPT monotonically decays with the increase of $h$. At the same time the maximum of the stationary density does not move and it is located at $x\approx 1$, see top panel of Fig. 10. For $n=4$, the non-monotonous dependence of the MFPT, with a shallow minimum at $h\approx 0.5$, is due to changes in the shape of the stationary states, see top panel of Fig. 11.

The significant differences between the parabolic and the quartic potential are visible on the level of stationary states, see top panel of Figs. 10 and 11. Analogously, like for the potentials considered in Secs. II.2 and II.3, if the barrier is high enough, the stationary state in the quartic potential has four modes. At the same time, the stationary state for the parabolic potential is characterized by two modes located in minima of the potential, i.e. $x=\pm 1$. In order to elucidate further differences between both potentials we have studied the distribution of the last hitting points.

The last hitting point $x_{\mathrm{last}}$ is the last point visited before escape from the domain of motion, i.e. it is $x(\tau_{i}-\Delta t)$, where $\Delta t$ is the integration time step and $\tau_{i}$ is the first passage time associated with this particular realization of the trajectory. For the selected setup, it is the last point such that $x>0$. Due to presence of noise in Eq. (1) trajectories are randomized. Consequently, the first passage time $\tau$ and the last hitting point $x_{\mathrm{last}}$ are random variables. The distribution of the last hitting points for both setups is very different, compare middle panels of Figs. 10 and 11. The distribution of last hitting points is related to the stationary state and the initial condition. A particle is more likely to escape from points where a probability of finding the particle is larger. Therefore, the last hitting point density follow very different pattern for $n=2$ and $n=4$. The most probable escape scenario, even for finite $\sigma$, is the escape via a single long jump Ditlevsen (1999); Imkeller and Pavlyukevich (2006a, b) because in the weak noise limit the Lévy noise can be decomposed into the Gaussian white noise and the Poisson compound process. This does not mean that the particle escapes immediately, because it can be wandering due to the Gaussian component of the noise. For the parabolic potential, if the barrier has an appropriate height, the escape from the potential well typically is not immediate ($P(\tau>0.1)>96\%$) and most probably finishes from the minimum of the potential. For the quartic potential the situation is more complicated due to two modes of the stationary state within a potential well. Thus, for the high enough potential barrier, the last hitting point density has three modes: two modes in points where the stationary state has modes and the additional mode in the minimum of the potential, due to immediate escapes from the starting point (initial condition). Nevertheless, a immediate escape is not very probable as $P(\tau<0.1)<4\%$, while the integration time step is $10^{-5}$. The peak arising due to the initial condition can be diminished or removed by the change in the initial condition. For instance for $x(0)$ uniformly sampled from $(0.5,1.5)$, i.e. $x(0)\sim\mathcal{U}(0.5,1.5)$, there is no third peak, compare middle and bottom panels of Fig. 11. The height of peaks in the last hitting density is reversed in comparison to the height of peaks in the stationary state, because it is easier to escape from the point which is closer to the boundary. Moreover, for the decreasing barrier height, the height of peaks located at $|x|\neq 1$ decrease and finally they disappear due to the larger role played by trajectories approaching the boundary separating states. For a very low potential barrier, regardless of a potential type, the particle can also approach the vicinity of the point separating states, due to “weak” noise pulses ruled by the central part of the Lévy distribution. Therefore, in addition to the peak in the minimum of the potential, the peak at the boundary emerges.

Lastly, we have compared escape scenarios in both potentials. We have estimated probabilities of surmounting the potential barrier via a sequence of given number of jumps into the direction of the boundary. Tab. 1 presents the ratio of probabilities of the passage over the barrier in a single jump $p(1)$ and two consecutive jumps $p(2)$, i.e. $p(1)/p(2)$. From numerical simulations, it looks that the results included in Tab. 1 are robust with respect to the integration time step as long as the integration time step is small enough. Tab. 1 confirms the hypothesis that the main scenario of the Lévy noise induced escape is an escape via a single long jump Ditlevsen (1999); Imkeller and Pavlyukevich (2006a, b); Vezzani et al. (2018). With the increasing barrier height, the escape via a single long jump becomes dominating. With the increasing value of the stability index $\alpha$ longer sequences are observed. Finally, in the limits of $\alpha=2$, i.e. for the Gaussian white noise, the escape takes place in a sequence of jumps. This effect is fully coherent with the decomposition of the Lévy noise into the Gaussian white noise and the Poisson compound process Imkeller and Pavlyukevich (2006a, b); Imkeller et al. (2010). In the weak noise limits (small $\sigma$), transitions over the potential barrier are due to Poissonian component of the noise.

IV Summary and Conclusions

We have studied stationary states in single-well and double-well potentials. It has been demonstrated that, for a sufficiently high potential barrier separating minima of double-well potentials, stationary states for systems perturbed by the Lévy noises can be at least bimodal, due to the possibility of producing more than one maximum of the stationary state within each potential well. In order to produce two maxima of the stationary density within a potential well, the potential well needs to be steeper than parabolic. For instance, for the double-well potential build from two single-well potentials of $x^{4}$ type, the stationary density can have four modal values. If single-well potentials are of $x^{2}$ type, the stationary state can be bimodal, with modal values located in potential wells. To conclude, the number of modal values is sensitive to the height of the potential barrier separating potential wells and the leading term describing the potential shape in the vicinity of its minima.

In general, it also possible to construct double-well potentials producing more than four modal values of the stationary state. Such potentials can be built from two single-well potentials of the type considered in Capała and Dybiec (2019), i.e. from single-well potentials with non-monotonous dependence of the potential curvature. The multimodality of stationary states can be also produced by the combined action of the Gaussian white noise and the Markovian dichotomous noise Dybiec and Schimansky-Geier (2007); Calisto et al. (2017). The dynamically induced multimodalidy emerges due to sliding to minima of the altered potential Dybiec and Schimansky-Geier (2007) while the Lévy noise induced multimodality is produced due to the competition between random long excursion and the deterministic sliding towards minima of the static potential, which is interrupted by long jumps.

Finally, we have studied the role of the multimodality on an escape from a potential well. We have observed that the escape process from a quartic potential well is slightly faster than the escape from the parabolic potential well, because particles are more likely to concentrate closer to the boundary, due to inner modes of the stationary state. The concentration of the probability mass closer to the boundary allow for faster escape because a shorter jump is more likely to occur than a longer one. This effect is the most pronounced for moderate barrier heights. In addition to stationary states, the most striking differences are visible on the level of last hitting point densities, which follow very different patterns due to different shapes of stationary states.

We have performed our considerations for the Cauchy noise but they can be extended to other Lévy noises. On the one hand, such an extension is straight forward. On the other hand, we are not expecting to observe different results (except weakening of effects for $1<\alpha<2$) that one described within the current manuscript. We have used the Cauchy density due to analytically known stationary states, which can be compared with numerical simulations.

Appendix A Stationary state for $V(x)={x^{n}}/{n}$ with $n=4$ and $\sigma\neq 1$

The quartic oscillator $V(x)=\frac{x^{4}}{4}$ is described by the following Langevin equation is

[TABLE]

As proved in Chechkin et al. (2002, 2003a, 2004, 2006, 2008), for $\alpha=1$, the stationary state is

[TABLE]

The stationary state for

[TABLE]

can be re-constructed by rescaling the space variable $x$ and time $t$. Let us introduce

[TABLE]

and

[TABLE]

The non-trivial part is the noise transformation

[TABLE]

where $L(t)$ is the $\alpha$-stable Lévy motion, which is $\frac{1}{\alpha}$ self-similar process, whose derivative is the $\alpha$-stable (Lévy) noise $\zeta_{\alpha}(t)$. Eq. (21) transforms into

[TABLE]

First, we multiply both sides of Eq. (25) by $\frac{t_{0}}{x_{0}}$. Next, we require that prefactors in the deterministic force and in the random force (noise) are equal to unity

[TABLE]

For $\alpha=1$, the solution of the above equation is

[TABLE]

After rescaling of variables Eq. (21) takes the form of Eq. (19) with $x$ and $t$ interchanged with $\tilde{x}$ and $\tilde{t}$. Consequently, the stationary state is given by Eq. (20) with $x$ replaced by $\tilde{x}$ because the rescaling of time does not change the stationary state. The stationary state for $x$ can be reconstructed by the change of variables

[TABLE]

Appendix B Stationary states for $V(x)=\lambda nx^{n}$ with $n=4$ and $n=2$

For $V(x)=\lambda nx^{n}$, the Langevin equation takes the following form

[TABLE]

Please note that the scale parameter of noise $\zeta_{\alpha}$ is set to 1. Thus $\sigma\zeta_{\alpha}$ acts as the $\alpha$-stable noise with the scale parameter $\sigma$. By the appropriate change of variables, the prefactor in the deterministic force can be incorporated into the $\sigma$ parameter allowing us to use results of Appendix A. Let us substitute $x=ay$ in Eq. (B)

[TABLE]

and divide it by $a$

[TABLE]

For $n\neq 2$, the value of $a$ is determined by the condition

[TABLE]

resulting in

[TABLE]

Therefore, Eq. (31) transforms into

[TABLE]

where

[TABLE]

Consequently, from the stationary solution $f(y)$ of Eq. (34) the stationary state $f(x)$ of Eq. (B) can be obtained by the transformation of variables

[TABLE]

In particular for $n=4$, $\lambda=1$ and $\alpha=1$ we have $a=\frac{1}{4}$ and $\sigma^{\prime}=4\sigma$. In the special case of the quartic Cauchy oscillator the stationary solution is

[TABLE]

Similar calculations can be performed for any $n>0$, however for $n=2$ one needs to transform both space and time in the full Langevin equation. Following calculations from Appendix A, for $n=2$ with $\alpha=1$, one obtains $t_{0}=1/4\lambda$ and $x_{0}=\sigma/4\lambda$. Consequently, the stationary state is the Cauchy density

[TABLE]

Appendix C Numerical methods

The Langevin equation, see Eq. (1)

[TABLE]

is integrated with the Euler-Maryuama scheme

[TABLE]

where $\zeta_{t}$ represents independent, identically distributed random variables Chambers et al. (1976); Weron and Weron (1995); Weron (1996) following the symmetric $\alpha$-stable density Janicki and Weron (1994b); Janicki (1996). The Euler-Maryuama scheme assures proper interpretation of stochastic integrals Burrage et al. (2006), which become more complex due to discontinuity of Lévy flights. Eq. (C) is accompanied with the initial condition on $x(0)$ which can be deterministic, i.e. $x(0)=x_{0}$, or random, e.g. $x(0)\sim\mathcal{U}(a,b)$. In the absence of impenetrable boundaries, stationary states are independent to the selection of $x(0)$, while for the first escape problems results are usually sensitive to the choice of $x(0)$.

Scheme (C) is used both for the construction of stationary states and examination of first escape problems. Crucial parameters in the approximation (C) are the integration time step $\Delta t$ and the number of repetitions $N$. The number of repetitions controls the statistical error of the approximation and fluctuations of the histogram. The integration time step controls the systematic error. These parameters are either selected by comparison of results of computer simulations with theoretical predictions or by self-consistency tests. Therefore, we have performed a series of simulations with the decreasing integration time step. This sequence was stopped when results were coherent with theoretical predictions or statistically equivalent to results with a larger integration time step. After fixing $\Delta t$ we have increased the number of repetitions in order to decrease statistical errors. In order to asses whether a stationary state was reached we have measured the interquantile widths. Simulations were performed till the time $t_{\mathrm{max}}$ which is long enough to assure that the interquantile width is not increasing. The interquantile width is more robust than the standard deviation because it can be defined also for distributions with the diverging variance.

The problem of first escape is also studied with the help of Eq. (C). After fixing $x(0)$ the simulation is performed as long as the particle has not crossed the prescribed boundary. Here, as long as $x>0$ . From many simulations of escape events we have recorded a series of first passage times $\tau_{i}$ and last hitting points $x_{\mathrm{last}}$, i.e. positions of the last visited point before leaving the domain of motion. In the next step from these series we can calculate the mean first passage time $\langle\tau\rangle$ along with its error and estimate statistics of last hitting points, e.g. $p(x_{\mathrm{last}})$.

Within simulations the integration time step varied between $10^{-3}$ – $10^{-5}$, while the number of repetitions was adjusted to $10^{5}$ – $10^{8}$. Stationary states were constructed with a smaller integration time step and a larger number of repetitions because they were constructed using CUDA on GPU (graphics cards). The MFPT was estimated on conventional CPUs, consequently statistics was poorer than for stationary states.

Acknowledgements.

This project was supported by the National Science Center grant (2014/13/B/ST2/02014). This research was supported in part by PLGrid Infrastructure. Computer simulations have been performed at the Academic Computer Center Cyfronet, AGH University of Science and Technology (Kraków, Poland) under CPU grant DynStoch.