Noise induced effects in multivariate systems driven by Lévy noise

(1)

Jagiellonian University

Faculty of Physics, Astronomy and Applied Computer Science

Noise induced effects in multivariate systems

driven by Lévy noises

Thesis for the degree of Doctor of Philosophy in Physics

Krzysztof Szczepaniec

supervisor:

(2)

(3)

Wydział Fizyki, Astronomii i Informatyki Stosowanej Uniwersytet Jagiello´nski

O´swiadczenie

Ja ni˙zej podpisany Krzysztof Szczepaniec (nr indeksu: 1014987) doktorant Wydziału Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiello´nskiego o´swiadczam, ˙ze przedło˙zona przeze mnie rozprawa doktorska pt. „Noise induced effects in multivariate systems driven by Lévy noises” jest oryginalna i przedstawia wyniki bada´n wykonanych przeze mnie osobi´scie, pod kierunkiem dr. hab. Bartłomieja Dybca. Prac˛e napisałem samodzielnie.

O´swiadczam, ˙ze moja rozprawa doktorska została opracowana zgodnie z Ustaw ˛a

o prawie autorskim i prawach pokrewnych z dnia 4 lutego 1994 r. (Dziennik Ustaw 1994 nr 24 poz. 83 wraz z pó´zniejszymi zmianami). Jestem ´swiadom, ˙ze niezgodno´s´c

ni-niejszego o´swiadczenia z prawd ˛a ujawniona w dowolnym czasie, niezale˙znie od skutków

prawnych wynikaj ˛acych z ww. ustawy, mo˙ze spowodowa´c uniewa˙znienie stopnia

naby-tego na podstawie tej rozprawy.

...

(4)

(5)

v Streszczenie

Niniejsza praca podejmuje problem opisu i numerycznego badania efektów

indu-kowanych wielowymiarowymi szumami α-stabilnym (szumami Lévy’ego) w układach

fizycznych. Na poziomie pojedynczej cz ˛astki, układy takie opisywane s ˛a przez

rów-nanie Langevina – stochastyczne rówrów-nanie ró˙zniczkowe – w którym wyst˛epujeα-stabilny

szum Lévy’ego. Na poziomie zespołu cz ˛astek równaniu Langevina odpowiada

(ułam-kowe) równanie Smoluchowskiego (Smoluchowskiego-Fokkera-Plancka), które opisuje

ewolucj˛e g˛esto´sci prawdopodobie´nstwa. Wyst˛epuj ˛ace w nim pochodne niecałkowitego

rz˛edu (ułamkowe) pojawiaj ˛a si˛e ze wzgl˛edu na szczególny rodzaj szumu w równaniu

Langevina. Rodzaj szumuα-stabilnego i jego parametry determinuj ˛a typ ułamkowego

operatora ró˙zniczkowego w równaniu Smoluchowskiego.

W pracy zostały przeanalizowane dwa podstawowe rodzaje wielowymiarowych

szu-mów α-stabilnych, zadanych poprzez dwie ró˙zne miary spektralne: jednorodn ˛a, ci ˛agł ˛a

miar˛e spektraln ˛a, okre´slon ˛a na okr˛egu jednostkowym, oraz dyskretn ˛a miar˛e spektraln ˛a,

zlokalizowan ˛a na przeci˛eciach okr˛egu jednostkowego z osiami układu współrz˛ednych.

Taki wybór miar spektralnych, w przypadku cz ˛astki swobodnej, prowadzi odpowiednio

do izotropowych i kartezja´nskich lotów Lévy’ego. Powy˙zsze miary spektralne stanowi ˛a

dwa najprostsze przypadki, posiadaj ˛ace zwart ˛a analityczn ˛a posta´c równania

Smoluchow-skiego, jednocze´snie s ˛a wystarczaj ˛ace do zademonstrowania wielu efektów wywołanych

szumamiα-stabilnymi. Przeprowadzona analiza zjawisk indukowanych szumami, oparta

na metodach Monte Carlo, obejmowała zagadnienia (i) ucieczki cz ˛astki swobodnej z

ogra-niczonych obszarów, (ii) istnienia stanów stacjonarnych w dwuwymiarowych studni-ach potencjału, (iii) aktywacji rezonansowej w dwu- i trójwymiarowych dychotomicznie

fluktuuj ˛acych potencjałach oraz (iv) rezonansu stochastycznego w dwuwymiarowej,

dwu-dołkowej studni potencjału.

Niniejsza praca pokazuje, ˙ze efekty wywołane szumami α-stabilnymi, wyst˛epuj ˛ace

w układach jednowymiarowych mog ˛a by´c obserwowane równie˙z w układach o wi˛ekszej

liczbie wymiarów. Pokazano równie˙z, ˙ze techniki i miary stosowane do opisu i

bada-nia efektów wywołanych szumami w układach jednowymiarowych mog ˛a by´c skutecznie

(6)

Abstract

The thesis addresses the problem of description and numerical study of effects

in-duced by multivariateα-stable Lévy type noise in physical systems. At the single particle

(trajectory) level, such systems are described by the Langevin equation (stochastic

differ-ential equation) with theα-stable noise component. At the ensemble level, the Langevin

equation corresponds to the (fractional) Smoluchowski (Smoluchowski-Fokker-Planck) equation, which describes the evolution of the probability density. Fractional derivatives in the Smoluchowski equation appear due to the special kind of noise in the Langevin

equation. Therefore, the α-stable noise type and its parameters define the type of the

fractional differential operator in the Smoluchowski equation.

In this work two basic types of multivariateα-stable noise, corresponding to two

spe-cial spectral measures, were analyzed. Namely, the uniform continuous spectral measure, defined on the unit circle, and the discrete spectral measure, concentrated on the intersec-tions of the unit circle with the axes. In particular, for a free particle, the choice of spectral measure leads to isotropic and Cartesian Lévy flights respectively. These spectral mea-sures are two simplest cases, resulting in compact analytical forms of the Smoluchowski equation. At the same time they are sufficient to demonstrate various effects induced by

theα-stable noise. The Monte Carlo analysis of the noise-induced phenomena focused on

(i) escape of a free particle from bounded domains, (ii) existence of stationary states for two-dimensional potential wells, (iii) resonant activation in two- and three-dimensional dichotomously fluctuating potentials and (iv) stochastic resonance in a two-dimensional, double-well potential.

The thesis shows that effects induced byα-stable noise, recorded in one-dimensional

systems can be also observed in multidimensional systems. Moreover, it is shown that techniques and measures used to describe and study noise-induced effects in one-dimen-sional systems can be effectively extended and applied in multidimenone-dimen-sional domains.

(7)

List of Figures

1.1 Univariate Gaussian distributions . . . 6

1.2 Sampleα-stable distributions . . . 13

1.3 2Dα-stable densities with α = 1 . . . 18

1.4 1D linear potential setup for the resonant activation . . . 24

1.5 Sample MFPT for resonant activation . . . 25

1.6 Sample 1D periodically modulated, double well potential. . . 26

1.7 Sample power spectrum and SNR for stochastic resonance . . . 27

2.1 MFPT for the 1D interval . . . 33

2.2 Trajectories of bivariate Lévy flights . . . 36

2.3 MFPT for the 2D disk (spherical and Cartesian LF) . . . 37

2.4 MFPT for hyperspheres (spherical LF) . . . 41

2.5 Ratio of MFPTs for 1D intervals and hyperspheres (spherical LF) . . . . 42

2.6 α maximizing MFPT for spherical LF in 1D and 2D . . . 43

2.7 MFPT for escape from hyperspheres (Cartesian LF) . . . 44

2.8 Ratio of MFPTs for intervals and hyperspheres (Cartesian LF) . . . 45

2.9 MFPT for the 2D disk with half of edge absorbing . . . 47

2.10 MFPT for the 2D disk with a small part of edge absorbing . . . 48

2.11 MFPT for hypercubes . . . 52

2.12 Rescaled MFPT for hypercubes . . . 52

2.13 MFPT for hypercubes with edge2L = 10 . . . 53

2.14 MFPT for hypercubes with edge2L = 20 . . . 53 ix

(10)

2.15 Survival probabilities (escape from the hypercube) . . . 57

2.16 Rescaled survival probabilities (escape from the hypercube) . . . 58

2.17 Sample bivariate trajectories in a harmonic potential . . . 61

2.18 Stationary states in 2D systems driven by Gaussian white noise . . . 67

2.19 Stationary states in 2D systems driven by Cauchy noise . . . 68

2.20 Stationary states for bivariateα-stable noise with α = 1.9 . . . 69

2.25 CCD for marginal densities for various potentials . . . 72

2.26 2D linear potential setup for the resonant activation . . . 74

2.27 Sample trajectories of bivariate LF in a switching potential . . . 75

2.28 Various measures of RA in 2D (1 < α < 2) . . . 78

2.29 Various measures of RA in 2D (α < 1) . . . 79

2.30 Various measures of RA in 3D (1 < α < 2) . . . 81

2.31 Various measures of RA in 3D (α < 1) . . . 82

2.32 Strength of the RA in 2D and 3D . . . 82

2.33 Periodically modulated 2D double-well potential . . . 84

2.34 Sample trajectories of bivariate LF in a double-well potential . . . 85

2.35 Periodic response_{hx(t)i . . . 86}

2.36 AmplitudeAmaxof the periodic response . . . 87

2.37 Mean residence time in a double-well potential . . . 88

2.38 Residence time distributions in a double-well potential . . . 89

2.39 Area under the “first peak” of the residence time distribution . . . 90

(11)

List of Abbreviations

1D/2D/3D one/two/three dimension

CCD complementary cumulative distribution

CLT central limit theorem

FPT first passage time

GWN Gaussian white noise

LE Langevin equation

LF Lévy flight

MFPT mean first passage time

MSD mean square displacement

PDF probability density function

RA resonant activation

RTD residence time distribution

SFP Smoluchowski-Fokker-Planck (equation)

SNR signal-to-noise ratio

SODE stochastic ordinary differential equation

SPA spectral power amplification

SR stochastic resonance

(12)

(13)

Introduction

In gas and fluids, atoms (or molecules) are free to roam around, constantly bumping into each other and exercising local changes in density. If we now observe a much larger particle immersed in a fluid, such a particle experiences constant collisions with the fluid’s matter, resulting in a non-zero force acting upon it, and as a result, it seems to move in random directions.

Such a phenomenon was reported first by Dutch physiologist, biologist and chemist Jan Ingenhousz in 1785, who observed irregular motion of a coal dust on the surface of alcohol. The very same phenomenon was also reported by Robert Brown in 1827 [1]. Robert Brown was a Scottish botanist who studied motion of pollen grains immersed in water. The behaviour observed by Brown – the irregular, jittery motion of particles suspended in water is now called the Brownian motion. Analytical descriptions of the diffusion came independently from Albert Einstein (1905) [2] and Marian Smoluchowski (1906) [3]. Works of Einstein and Smoluchowski were verified experimentally, in partic-ular, by the French physicist Jean Baptiste Perrin in 1908, who received a Nobel prize in 1926 for confirmation of atomic nature of matter.

A straightforward mathematical description and modeling of the Brownian motion would be impossible due to the sheer number of particles. In one mole of substance one

have aroundNA ≈ 6.02 × 1023 (the Avogadro number) of particles. For a full,

Newto-nian description of motion for every particle, one would need6-dimensional phase space,

therefore, description of the entire physical system becomes a many-body problem, with

approximately1024 _{degrees of freedom. Fortunately, due to the time scale separation}

be-tween the observed Brownian particle (Brownian walker) and the molecules of the system 1

(14)

that interact with it, it is possible to use an effective statistical description. If we assume that interactions between the system and the random walker are independent (uncorre-lated) and identically distributed with the finite variance, the complex interactions can be simply modelled by a Gaussian white noise. This is possible due to the power of the central limit theorem, which states that the sum of independent, identically distributed random variables (characterized by a finite second moment) converges to a Gaussian dis-tribution [4]. Although the Gaussian noise is a very simple approximation, it provides an effective and useful description of many phenomena.

The concept of noise, as an effective approximation of complex interactions, is one of the crucial notions in statistical physics. Contrary to the everyday experience, noise in dynamical stochastic systems is more than an unwanted ingredient. Presence of noise is the key element of the so called noise induced effects, like stochastic resonance [5–7], res-onant activation [8–10] and noise-enhanced stability [11, 12] to name a few. Assumptions of the central limit theorem can be relaxed, therefore expanding the field of “normal” dif-fusion into the domain of “abnormal” – anomalous difdif-fusion. The central limit theorem relies on the independence and boundness (finite variance) of summands. If the distribu-tion of random variables is heavy-tailed, and therefore has an infinite variance, classical central limit theorem no longer holds. This give rise to a generalized central limit theo-rem. It has been proven that a sum of independent, identically distributed variables with infinite variance converges to a more general Lévy stable distribution [4], for which the Gaussian distribution is a special limiting case.

Heavy-tailed, non-Gaussian fluctuations have also been observed in multiple exper-iments and setups. Seminal example are the search and foraging patterns following the Lévy distribution [13–16] or anomalous transport effects [17, 18]. Stable distributions were also observed in laser cooling of atoms [19], diffusion of photons in hot vapour [20] and diffusion of supercold atoms in optical lattice [21]. Modeling with Lévy stable noises have also been used in economy [22, 23], epidemiology [24] or search optimization [25]. Lévy flights of photons have been observed in special optical materials [26].

(15)

INTRODUCTION 3 noise induced effects. The fact that the presence of noise can enhance the observed effect in the system is now a well-established statement. In the studies about the noise induced effects, a noise with heavy-tailed distribution is also often present. The best known exam-ples of noise induced effects, among others, are: stochastic resonance (SR) and resonant activation (RA). Stochastic resonance [5, 6, 27] is the amplification of a weak signal by noise, while the resonant activation [8–10] is the facilitation of noise assisted escape over a fluctuating potential barrier. Resonant effects in nature have been observed in various sce-narios. One of the most interesting examples comes from paleoclimatology [28, 29] and noise induced resonant effect in climate changes [30]. Resonant effects induced by the stable Lévy-type noise in various setups are the topic of earlier studies, and have been extensively analyzed in the recent years [31–38].

Modelling of the systems with noise is often done by the use of the Langevin equation, proposed by Paul Langevin [39, 40]. In the underdamped limit, the Langevin equation is a full Newtonian equation of motion with a stochastic term, i.e. it is a stochastic ordinary differential equation of the second order. In the overdamped limit, velocity can be adia-batically eliminated [5,41], resulting in the overdamped, first order equation for evolution of the position. The Langevin equation describes the dynamics of a single particle sub-jected to the deterministic and stochastic forces. Due to the presence of stochastic force, even for the fixed initial conditions, the trajectories of a single particle differ. Therefore, description on a single trajectory level is accomplished by the description on the ensemble level, which is provided by the Smoluchowski (or Smoluchowski-Fokker-Planck) equa-tion. The Smoluchowski equation describes the evolution of the probability density of trajectories associated with the Langevin equation. Within this thesis we mainly use the Langevin approach since it allows for a numerical evaluation and estimation of various statistical properties of the system [42–44] without facing complexity of the diffusion equation. Therefore, it is an alternative to the Smoluchowski equation, which is usually much more complicated to solve.

The main goal of this thesis is to extend the previous research on the noise induced effects into a multidimensional space and at the same time to prove that the noise induced

(16)

resonant effects, observed in 1D systems, can also be observed in multidimensional

sys-tems. All stochastic effects which are studied within this thesis are induced by anα-stable

Lévy-type noise, which gives rise to fractional derivatives in the Smoluchowski equation. Extensive use of the Monte Carlo simulations, along with analytical analysis have been used to gain insight into those complicated systems and determine, if the pronounced

resonant effects can be observed. The multivariateα-stable Lévy-type noise introduces

another difficulty, which is the need to define the spectral measure, replacing the scale and

asymmetry parameters of 1Dα-stable noises. Various spectral measures were considered

in the research, in order to provide a better understanding of the dynamics and illustrate how the choice of the spectral measure can influence properties of stochastic systems.

Basic concepts and preliminary information are included in Chapter 1. Original re-sults are presented in Chapter 2 which discusses the noise induced effects: escape from bounded domains (Chaps. 2.1 and 2.2), stationary states in 2D systems (Chap. 2.3), res-onant activation (Chap. 2.4) and stochastic resonance (Chap. 2.5). The thesis is closed with conclusions and outlook. The thesis has a form of a traditional dissertation, because it allows for extensive and comprehensive discussion of obtained results.

(17)

Chapter 1 Preliminaries

1.1 Normal diffusion

Diffusion is a fundamental transport process in fluids and gases. It is random in nature and arises as a consequence of the microscopic nature of matter. Diffusive transport is caused by continuous collisions between molecules of the fluid with the molecules of a diffusing substance. Fundamental quantity of interest is the concentration, which is an amount of given substance within a volume of mixture. For example, when a drop of ink is poured into water, the ink is concentrated in a small volume. With time, the ink diffuses in water, until the stationary state is reached, i.e. when the concentration of ink is equal within entire volume of water.

The diffusion process can be described by the diffusion equation. In the absence of deterministic forces, the one-dimensional diffusion equation has the form

∂u(x, t)

∂t = D

∂2_{u(x, t)}

∂x2 , (1.1)

whereD is the diffusion coefficient and u(x, t) is the concentration of substance at

posi-tionx at time t. u(x, t) can also be interpreted as a probability distribution of diffusing

particles (probability of finding a particle at position x at time t). The diffusion

coef-ficient is a proportionality constant, depending on size of molecules, temperature, and 5

(18)

other characteristics of the material. Mathematically, the Eq. (1.1) is a partial differential equation, which can be solved explicitly, if provided with adequate initial and boundary conditions. In the unrestricted space, the solution of Eq. (1.1) with the initial condition u(x, 0) = δ(x) is

u(x, t) = √ 1

4πDte

−x2

4Dt. (1.2)

The Eq. (1.2) is the Gaussian distribution with the variance σ2 _{= 2Dt and the mean}

µ = 0. The variance is time-dependent because the concentration of substance changes in time, see Fig. 1.1.

0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4 u (x, t) x t = 0.1 t = 0.2 t = 0.3 t = 0.6

Figure 1.1: Solution of Eq. (1.1) given by Eq. (1.2). The concentrationu(x, t) is the

Gaussian distribution with the time-dependent varianceσ2 _{= 2Dt and the mean µ = 0.}

Various lines presents concentrations at various times. The diffusion coefficient isD = 1.

Diffusion equation, Eq. (1.1), can also be generalized to a multidimensional domain

in a straightforward manner. Assuming that the diffusion coefficientD is constant, the

diffusion equation can be written as ∂u(x, t)

∂t = D∆u(x, t), (1.3)

(19)

1.1. NORMAL DIFFUSION 7

1.1.1 Gaussian distribution and the central limit theorem

The Normal distribution (Gaussian distribution), usually denoted asN (µ, σ2_{), is arguably}

one of the most important distribution in statistics. The probability density function of the (one-dimensional) normal distribution is

pµ,σ(x) = 1 √ 2πσ2e −(x−µ)2 2σ2 , (1.4)

whereµ is the mean value (first moment) and σ2 _{is the variance (second central moment),}

whileσ is the standard deviation. The characteristic function

φ(k) =

Z +∞

−∞

eikx_p

µ,σ(x)dx (1.5)

of the normal distribution is

φ(k) = eiµk−σ2k2₂ _.

(1.6) The distribution of the sum of two independent Gaussian variables can be calculated from the characteristic function given by Eq. (1.6). From Eq. (1.5) it implies that the char-acteristic function of the sum of independent, identically distributed random variables is the product of characteristic functions of every variable. Therefore, the sum of two

in-dependent Gaussian variables is also normally distributed. Namely, ifx1 ∼ N(µ1, σ21)

andx2 ∼ N(µ2, σ22), then x1 + x2 ∼ N(µ1+ µ2, σ12+ σ22). Analogously a1x1+ a2x2 ∼

N (a1µ1+ a2µ2, a12σ12+ a22σ22).

The central limit theorem states that the sum of large number of independent and identically distributed random variables (with finite variance) converges to the Gaussian distribution [4, 41, 45]. Therefore, the Gaussian distribution is a very good approximation in statistical analysis of complex processes. Seminal example is the Brownian motion, where a test particle experiences stochastic force due to the random interactions with the particles (molecules) of the environment. Assuming that collisions are independent and bounded, i.e leading to displacement with a finite variance, the Brownian motion can be

(20)

incrementsW (t)_{− W (t}0₎

∼ N(0, t − t0_{) [46].}

1.1.2 Langevin equation

Alternatively to the diffusion equation, the process of diffusion can also be described on a single particle trajectory level. The Langevin equation, proposed by Paul Langevin [39, 40] describes a single particle subjected to random and deterministic forces. The

Langevin equation for a particle with massm immersed in a fluid is

md

2_r(t)

dt2 =−γ

dr(t)

dt + f (t), (1.7)

whereγ is the friction coefficient and f (t) is the stochastic term representing the random

force originating from collisions with other particles. Eq. (1.7) is stochastic in nature and two consecutive realizations of r(t) with the same initial conditions are not the same. Therefore, a valuable information can be extracted by calculating averages for an ensem-ble of possiensem-ble realizations of r(t).

In the simplest case (the one-dimensional Brownian motion) the stochastic termf (t)

can be approximated by the Gaussian white noise. Every particle in the ensemble starts at x = 0. This corresponds to the initial condition u(x, 0) = δ(x) for the diffusion equation.

In the overdamped limit (with largeγ), when the velocity can be disregarded [41, 47], the

diffusion equation (1.1) with such an initial condition describes the time evolution of the probability distribution for the ensemble of particles described by the Langevin equation (1.7).

For the one-dimensional Brownian motion, the Gaussian white noiseξ(t) fulfills

hξ(t)i = 0 (1.8)

hξ(t)ξ(t0

)i = 2Dδ(t − t0

), (1.9)

whereD is the diffusion coefficient (strength of the fluctuations, see Chap. 1.1). The first

(21)

1.1. NORMAL DIFFUSION 9

likely to move in any direction. Multiplying both sides of Eq. (1.7) byx and taking the

ensemble average we get

m xd 2_x dt2 =_−γ xdx dt +_{hxξ(t)i .} (1.10)

The last term in the Eq. (1.10) dissappears, because_{hxξ(t)i = hxi hξ(t)i = 0, see [47]}

and Eq. (1.8). Eq. (1.10) can be rewritten as

mdhx ˙xi dt = kT − γ hx ˙xi , (1.11) where xd2_x dt2 = m _{d(x ˙}_x) dt − ˙x

2 _{and by the equipartition theorem} _{h ˙x}2_{i =} kT

m, see [48].

Solution of the Eq. (1.11) for_{hx ˙xi is}

hx ˙xi = Ce−t

+kT

γ , (1.12)

whereC is the integration constant and = _mγ. For the initial conditionx = 0 at t = 0

we getC = −kT

γ and Eq. (1.12) can be rewritten as

hx ˙xi = 1 − e−t kT

γ . (1.13)

Furthermore, because_{hx ˙xi =} 1₂_dtd_hx2_{i we get the MSD (mean square displacement) hx}2_i

of such process hx2 i = 2 t₋ 1 (1− e −t ) kT γ . (1.14)

The−1 _{is the characteristic time scale of the system, and in the limiting case of a long}

timet −1 _{the Eq. (1.14) reduces to}

hx2

i ' 2kT

γ t, (1.15)

which is the characteristic property of a random walk. Eq. (1.15) demonstrates that linear scaling of the mean square displacement is a typical property of many stochastic systems.

(22)

In the overdamped limit, where the damping of the particle movement is much larger

than the inertial termγdr(t)_dt

m

d2_r(t)

dt2

, the Langevin equation (1.7) can be simplified

[41, 49] to

dr(t)

dt = f (t), (1.16)

where f(t) is rescaled by 1

γ in comparison to Eq. (1.7) . The overdamped Langevin

equation is a very useful tool for numerical evaluation. If f(t) is the Gaussian white noise, the mean square displacement for Eq. (1.16) scales as predicted by Eq. (1.15), what is a direct consequence of Eqs. (1.8) and (1.9). Moreover, the Langevin equation (1.16) driven by the Gaussian white noise is associated with the diffusion equation (1.1), namely, the diffusion equation (1.1) describes the time evolution of the probability distribution of the ensemble of particles described by the Eq. (1.16).

Furthermore, the Eq. (1.16) can also be generalized to the situation when a particle experiences both stochastic and deterministic forces. When a particle moves in a potential V (r) the Langevin equation (1.16) attains the form

dr(t)

dt =−∇V (r) + f(t), (1.17)

where _{−∇V (r) is the deterministic force acting on a particle. Other generalizations}

are also possible, for example for time-dependent potentialsV (r, t) [35, 50, 51] or

non-potential forces.

1.2 Anomalous diffusion

Normal diffusion is characterized by the mean square displacement being a linear func-tion of time. It was eventually discovered that this is not an universal relafunc-tion. Various experiments [52,53] showed that in more general cases the mean square displacement can be a power function of time, i.e

(23)

1.2. ANOMALOUS DIFFUSION 11

where theδ is a parameter used to distinguish the type of diffusion. In general, for δ < 1

the diffusive process is called subdiffusion, while forδ > 1 we have superdiffusion. For

δ = 1 normal diffusion is recovered. Nevertheless, the discrimination of the diffusion type solely on Eq. (1.18) can be misleading [54, 55] because, in general realms, the com-petition between long jumps and long waiting times can hide anomalous diffusion. There-fore, discrimination of the diffusion type should take into account microscopic system’s dynamics.

Brownian motion is a natural consequence of the central limit theorem, which requires the random increments to be independent and ruled by a distribution with a finite variance. These assumptions can be weakened, resulting in an anomalous diffusion. Therefore, the generalized central limit theorem states that a sum of independent and identically distributed random variables drawn from a heavy-tailed distribution with power-law tails

of_|x|−α−1 type, with0 < α < 2, will tend to an α-stable distribution with the increasing

number of summands.

1.2.1 α-stable Lévy type noise and the generalized CLT

Random variable X is said to be stable if a linear combination of two (or more)

inde-pendent copies of the variable X has the same distribution as X (up to location and

scale) [56]. Formally, for any positive numbersa and b, there is a positive number c and

a real numberd such that

aX1+ bX2

d

= cX + d, (1.19)

whereX1 and X2 are independent copies of X and

d

= denotes equality in distributions.

Ifd = 0 random variables are said to be strictly stable. For example, condition (1.19) is

fulfilled by independent Gaussian random variables, see Eqs. (1.4) and (1.6).

Using characteristic functions, see Eq. (1.5), and their properties it is possible to rewrite Eq. (1.19) in the Fourier space

(24)

Eq. (1.20) can be solved [56–58], resulting in

φα(θ) =

 



exp_−σα_|θ|α₍₁_{− iβsign(θ) tan}πα

2 ) + iµθ

for α_{6= 1,}

exp−σ|θ| 1 + iβ2

πsign(θ) ln|θ| + iµθ for α = 1.

(1.21)

For the stability indexα = 2, the characteristic function, Eq. (1.21), reduces to the

char-acteristic function of the Gaussian distribution (see. Eq. (1.6)).

The α-stable distributions are characterized by 4 parameters, namely, the stability

indexα, skewness parameter β, scale σ and location parameter (shift) µ. Multiple

alter-native parametrization of the stable distribution are used in the literature. In this thesis,

the parametrizationSα(σ, β, µ) with the characteristic function (1.21) [56] is used. The

parameters of theα-stable distribution are constrained in the following way:

0 < α 6 2, −1 6 β 6 1, σ > 0,

µ_{∈ R.}

The stability indexα controls distribution tails, which for α < 2 are of|x|−α−1 _{type. The}

asymmetry parameterβ defines which distribution tail is heavier. For β = 0 distributions

are symmetric with respect toµ. µ is the location parameter and σ scales the overall

distri-bution width. Forα > 1 the α-stable densities are characterized by the finite mean, while

their variance diverges for everyα < 2. If a random variable X is distributed according to

the probability density with the characteristic function (1.21), it is schematically denoted as

X ∼ Sα(σ, β, µ). (1.22)

Effects of a various parameter on the shape of an α-stable distribution can be seen in

Fig. 1.2.

The Gaussian distribution N (µ, σ2_{) with the mean value µ and the variance σ}2 _(see

(25)

1.2. ANOMALOUS DIFFUSION 13 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −4 −2 0 2 4 Sα (σ = 1 ,β = 0 ,µ = 0) x α = 0.5 α = 1.0 α = 1.5 α = 2.0 0 0.05 0.1 0.15 0.2 0.25 0.3 −4 −2 0 2 4 Sα =1 .5 (σ = 1 ,β ,µ = 0) x β = −1 β = 0 β = 1

Figure 1.2: Sampleα-stable distributions with various values of the (a) stability index α

and (b) skewness parameterβ.

skewness parameterβ has no effect, but traditionally it is set to β = 0

N (µ, σ2)≡ S2(

σ √

2, 0, µ). (1.23)

The Gaussian distribution is the onlyα-stable density having all moments finite.

Proba-bility density function is given by Eq. (1.4).

For the stability indexα = 1 and the skewness parameter β = 0 an α-stable

distribu-tion becomes a Cauchy distribudistribu-tionC(µ, γ) with the location parameter µ and the scale

parameterγ > 0

C(µ, γ)≡ S1(γ, 0, µ). (1.24)

The probability density function of the Cauchy distribution is

cµ,γ(x) = 1 πγ 1 +x−µ γ 2 . (1.25)

It has power law tails of the _|x|12 type and it is symmetric with respect toµ.

Furthermore, forα = 1

2 and maximal positive skewness β = 1, an α-stable

distribu-tion reduces to a Lévy-Smirnov distribudistribu-tionF (µ, c) with the location parameter µ and the

scale parameterc > 0

F (µ, c)≡ S1

(26)

The probability density function of the Lévy-Smirnov distribution is (_{x > µ)} fµ,c(x) = r c 2π e−2(x−µ)c (x− µ)3/2. (1.27)

The Gaussian, Cauchy and Lévy-Smirnov distributions, are the only known cases, where the probability distribution function can be expressed in a simple analytical form.

The central limit theorem (Chap. 1.1.1) states that the sum of independent, identically distributed random variables with well-defined (finite) expected value and variance (first two moments), converges to a normal distribution. Assumption about finite variance can be relaxed [4], giving raise to the generalized central limit theorem, which states that the sum of independent, identically distributed random variables, with a power-law tail

distribution, decaying as _|x|α+11 where0 < α < 2 will tend to a stable distribution as the

number of variables grows. For_{α > 2 the sum tends to a Gaussian distribution, and the}

central limit theorem is recovered. Formally, following the definition of stable variable

Y , for a sequence of independent and identically distributed random variables X1, X2, . . .,

and sequence of positive numbersdnand real numbersanwe have

X1+ X2+· · · + Xn

dn

+ an

d

→ Y, (1.28)

where_{→ denotes a convergence in the distribution [4].}d

Arithmetic properties of theα-stable densities can be deduced from the characteristic

function (1.21), see also [56, 59, 60]:

• if Xi ∼ Sα(σi, βi, µi) then X1+ X2 ∼ Sα(σ, β, µ), (1.29) withσ = (σα 1 + σ2α)1/α,β = β1σα1+β2σα2 σα 1+σ2α andµ = µ1+ µ2. • if X ∼ Sα(σ, β, µ) then X + d∼ Sα(σ, β, µ + d); (1.30)

(27)

1.2. ANOMALOUS DIFFUSION 15

aX _∼

 



Sα(|a|σ, sign(a)β, aµ) for α6= 1,

S1(|a|σ, sign(a)β, aµ − 2_πaσβ ln|a|) for α = 1,

, (1.31) − X ∼ Sα(σ,−β, −µ); (1.32) • if Xj ∼ Sα(σ, β, µ) then aX1+bX2+c∼    Sα((aα+ bα)1/ασ, β, (a + b)µ + c) for α6= 1,

S1((a + b)σ, β, (a + b)µ− _π2σβ(a ln(a) + b ln(b)) + c) for α = 1,

, (1.33)

wherea, b∈ R+_and_c_{∈ R.}

1.2.2 Multivariate

α-stable variables

For some fundamental probability distributions, a generalization into multidimensional domain is very simple. For example the Gaussian distribution, Eq. (1.4), can be extended

tod dimensions in a straightforward manner

p(d)_µ,Σ(x) = 1 p(2π)d_{det Σ}e −1 2(x−µ) T_Σ−1_(x−µ) , (1.34)

where x is ad-dimensional random vector with the mean value µ and a covariance matrix

Σ.

Very different situation arises for theα-stable distribution (see Chap. 1.2.1). The

gen-eralization to a multivariate distribution requires the introduction of a spectral measure

Λ(·) [56], which encapsulates the information about the symmetry and scale of the

prob-ability density function. Therefore, the spectral measure replaces the skewness parameter β and the scale parameter σ of the univariate distribution, see Eq. (1.21).

Definition of a multivariate stable random variable is a natural extension of a

one-dimensional case. We define a random vector X = (X1, X2, . . . , Xd) in Rd. Following a

(28)

positive numbersa and b there is a positive number c and a vector D _{∈ R}d_{such that}

aX(1)_{+ bX}(2) d_{= cX + D,}

(1.35)

where X(1) and X(2) are independent copies of X, see [56].

Spectral measure

Let X = (X1, X2, . . . , Xd) be an α-stable random vector in Rd. Characteristic function

Φα(θ) of such random vector is then given by

Φα(θ) = Eeihθ,Xi = E h eiPd k=1θkXk i , (1.36)

which can be re written as

Φα(θ) =      expn−R Sd|hθ, si|

α_{1 − isign(hθ, si) tan}πα

2 Λ(ds) + ihθ, µ

0_io _for _α_{6= 1,}

expn−R

Sd|hθ, si|

α_{1 + i}2

πsign(hθ, si) ln(hθ, si) Λ(ds) + ihθ, µ

0_io _for _{α = 1,}

where_{hθ, si is a scalar product, and Λ(·) is called the spectral measure (compare Eq. (1.21)).}

The spectral measure is defined on a unit sphereSdin Rd[56]. ForSd=1 unit sphere, the

spectral measure is defined on a set of two points_{{−1, 1} and for S}d=2it is defined on a

unit circle. The spectral measure is the defining object of a multidimensional distribution, encapsulating skewness (β) and scale (σ) parameters that described the univariate stable distribution (see Chap. 1.2.1). The choice of a spectral measure defines the final distri-bution. If the spectral measure is symmetric, the stable distribution with such a spectral measure is said to be symmetric. Moreover, the spectral measure does not need to be

continuous. It can be discrete, i.e. concentrated on a discrete set of points onSd.

In general, components of a random vector X are not independent. Independency of the vector components is a special case, occurring if and only if the spectral measure

(29)

1.2. ANOMALOUS DIFFUSION 17 example, for a two-dimensional case that would be a set of four points:

{(1, 0), (−1, 0), (0, 1), (0, −1)}. (1.37)

Clear example of how the choice of the spectral measure influence the shape of the

prob-ability density is a 2D Cauchy distribution, i.e. α-stable density with the stability index

α = 1 (see Chap. 1.2.1). Discrete, symmetric spectral measure, concentrated on the

in-tersections of the axes with the unit circleS2 reduces to a two-dimensional distribution,

which is a product of two one-dimensional Cauchy distributions

p(x, y) = 1 π σ x2_{+ σ}2 × 1 π σ y2_{+ σ}2. (1.38)

While for the uniform continuous spectral measure, uniformly distributed on a unit circle

S2, theα-stable distribution is a two-dimensional Cauchy density

p(x, y) = 1

2π

σ

(x2_{+ y}2_{+ σ}2₎3/2. (1.39)

Fig. 1.3 illustrates the difference in shape of both 2D Cauchy densities. The left picture

depicts the α-stable density with a discrete spectral measure (see Eq. (1.38)), while the

right picture shows theα-stable density with a uniform continous spectral measure,

con-centrated on a unit circleS2 (see Eq. (1.39)). Both densities have the same value of the

stability indexα = 1.

Langevin equation with theα-stable Lévy-type noise

The motion of a free particle subjected to the action of anα-stable noise can be described

using the overdamped Langevin equation (see Chap. 1.1.2) dx

(30)

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

Figure 1.3: Effect of various spectral measures on a shape of a 2Dα-stable densities with

α = 1 (Cauchy densities). The left picture depicts the α-stable density with a discrete

symmetric spectral measure (see Eq. (1.38)), while the right picture shows theα-stable

density with a uniform continuous spectral measure, concentrated on a unit circleS2(see

Eq. (1.39)).

wherex(t) is a particle position and ζα(t) is an α-stable white noise, which is a formal

time-derivative of the α-stable motion ζα(t) =

dLα(t)

dt . Lα(t) is a generalization of the

Wiener process [61], i.e. a stochastic process with independent increments, distributed

according to anα-stable density. Eq. (1.40) can be easily generalized to higher number

of dimensions and situations with the deterministic force alongside the random force. If

the deterministic force, arising from the potentialV (x) is present, the Langevin equation

(1.40) attains the form

dx

dt =−V

0

(x) + σζα(t). (1.41)

And a multivariate generalization of Eq. (1.41) with theα-stable noise can be written as

dr(t)

dt =−∇V (r) + σζα(t), (1.42)

which is the straightforward generalization of the Langevin equation (1.17).

1.3 Smoluchowski-Fokker-Planck equation

The diffusion equation (see Chap. 1.1) describes the time evolution of the probability

(31)

prob-1.3. SMOLUCHOWSKI-FOKKER-PLANCK EQUATION 19

ability distributionp(x, t) describes the probability of finding a particle at position x at

time t. The classical diffusion equation has its limitation, i.e. it describes only

nor-mal diffusion, or correspondingly, spreading of an ensemble of particles following the Langevin equation with the Gaussian white noise term, see Eq. (1.16). Therefore, for

systems with theα-stable Lévy type noise, a more general equation must be used. The

Smoluchowski-Fokker-Planck [48] (SFP, sometimes called Fokker-Planck, or the Smolu-chowski equation) describes the time evolution of the probability density function of an ensemble of particles subjected to a deterministic and a random force.

For the symmetric,α-stable noise, the Eq. (1.40) can be associated with the fractional

Smoluchowski-Fokker-Planck [62–72] equation

∂p(x, t_|x0, 0)

∂t = σ

α∂αp(x, t|x0, 0)

∂|x|α , (1.43)

describing the evolution of the probability densityp(x, t|x0, 0) of finding a particle at

posi-tionx at time t, assuming the beginning of motion at position x0at time0. The term ∂

α

∂|x|α

is the fractional Riesz-Weil derivative, defined by its Fourier transform _F

h

∂α_p(x,t)

∂|x|α i

=

F−(−∆)α/2_{p(x, t)}_{= −|k|}α_{F [p(x, t)] [73, 74]. Analogously, Eq. (1.41) with the}

de-terministic force term is associated with the following fractional Smoluchowski-Fokker-Planck equation ∂p(x, t|x0, 0) ∂t = d dx[V 0 (x)p(x, t|x0, 0)] + σα ∂α_{p(x, t}_|x 0, 0) ∂|x|α . (1.44)

In case of a more general asymmetricα-stable noise the fractional

Smoluchowski-Fokker-Planck equation can be found in [74, 75].

The Smoluchowski-Fokker-Planck equation can also be generalized into a

multivari-ate domains. If the random term in the corresponding Langevin equation is anα-stable

Lévy type noise (see Eq. (1.42)), a special attention is required. The multivariateα-stable

noise requires the choice of a spectral measure (see Chap. 1.2.2), which is reflected in the corresponding 2D Smoluchowski-Fokker-Planck equation. The two-dimensional SFP

(32)

written as

∂p(r, t)

∂t =∇ · [∇V (r)p(r, t)] + σ

α

Ξp(r, t), (1.45)

where Ξ is a fractional operator arising due to two-dimensional α-stable noise [75].

The exact form of the fractional operator depends on the chosen spectral measure, see Chaps. 1.2.2 and 2.1. For the uniform continuous spectral measure, the fractional opera-tor has the form

Ξ=−(−∆)α/2

, (1.46)

and it is defined by its Fourier transform_F−(−∆)α/2p(r, t) = −|k|α_{F [p(r, t)]. For}

a discrete, symmetric spectral measure, concentrated on the intersections of the axes with

the unit circleS2, the fractional operator has the form

Ξ = ∂ α ∂|x|α + ∂α ∂|y|α, (1.47) where _∂|x|∂αα and ∂α

∂|y|α are the fractional Riesz-Weil derivatives, defined by the Fourier

transforms_F h ∂α_p(r,t) ∂|x|α i =−|k|α_{F [p(r, t)] and F}h∂α_p(r,t) ∂|y|α i =−|l|α_{F [p(r, t)]. The drift}

term _{∇ · [∇V (r)p(r, t)] in the Eq. (1.45) arises due to a deterministic force F (r) =}

−∇V (r), acting on a test particle, originating from the potential V (r). There is no known general solution for the Smoluchowski-Fokker-Planck equation, except for a few simples cases. Therefore, in many applications, the Monte Carlo simulations (see Chap. 1.4) can provide an alternative for estimating the time evolution of the probability density.

1.4 Monte Carlo methods for SODE

Numerical simulations of systems driven by noise can be conducted in terms of the Langevin equation (see Chap. 1.1.2), corresponding to the Smoluchowski-Fokker-Planck

equation. The Langevin equation is used to simulate an ensemble of trajectories x(t),

which can be used to estimate the statistical properties of the system. The univariate over-damped Langevin equation driven by the Gaussian white noise (see Eq. (1.16)) for a free particle driven by the Gaussian white noise can be numerically integrated [76–78] with

(33)

1.4. MONTE CARLO METHODS FOR SODE 21 the following Euler-Maruyama algorithm [79, 80]

xt+∆t = xt+

√

2D∆tξt, (1.48)

where ∆t is the time step of the numerical integration and ξt is the sequence of

inde-pendent and identically distributed random Gaussian variables N (0, 1). The numerical

scheme (Eq. (1.48)) can be extended for case with the deterministic force, originating

from the potentialV (x) in the straightforward manner:

xt+∆t= xt− V0(xt)∆t +

√

2D∆tξt. (1.49)

The equation (1.49) can be extended for a more general,α-stable noise. For the stability

index α < 2, the numerical recipe can be deduced from the algorithm for the

Gaus-sian white noise (1.49) and the arithmetic properties of theα-stable variables [81], see

Chap. 1.2.1. Numerical estimation of trajectories determined by Eq. (1.41) can be per-formed with the following scheme [42, 43, 60]

xt+∆t = xt− V0(xt)∆t + σ(∆t)1/αθα,t, (1.50)

whereθα,t is the sequence of independent and identically distributedα-stable Lévy-type

random variables Sα(β, 1, 0). For α = 2 the recipe (1.50) reproduces the numerical

scheme for the Gaussian white noise with√2 times larger σ, see Eq. (1.49).

Extension of Eq. (1.50) for the multivariate system is also straightforward, and reads

xt+∆t = xt− ∇V (xt)∆t + (∆t)1/αθα,t, (1.51)

where θα,tis the sequence of independent and identically distributed multivariateα-stable

random variables.d-dimensional multivariate α-stable random vector X with the discrete

spectral measure Λ(·) = d X j=1 γjδsj(·), (1.52)

(34)

where γj is the weight of the j component, and δsj(·) is the point mass at sj, can be numerically generated with the following recipe [82, 83]:

X =    Pd j=1(γj)1/αZjsj α6= 1 Pd j=1γj Zj+ 2_πln γj sj α = 1, (1.53)

whereZj, j ∈ (1, . . . , d) are the independent and identically distributed, totally skewed

one-dimensionalα-stable random variables, namely Zj ∼ Sα(1, 1, 0), see Chap. 1.2.1.

Multivariateα-stable random variables corresponding to the continuous uniform

spec-tral measures can be generated using the sub-Gaussian random vectors [56]. By

generat-ing anα-stable random variable A distributed as

A∼ Sα/2 (cos πα 4 ) 2/α , 1, 0, (1.54)

forα < 2 and a d-dimensional zero-mean Gaussian random vector

G= (G1, G2, . . . , Gd) , (1.55)

we get ad-dimensional α-stable random vector X as

X =√AG1, √ AG2, . . . , √ AGd . (1.56)

Vector (1.56) is anα-stable random vector distributed according to a multivariate α-stable

Lévy distribution with the stability indexα and the continuous uniform spectral measures.

1.5 Noise induced effects

Presence of noise can be beneficial for the systems’ dynamics. It can result in the so-called noise induced effects, i.e. effects manifesting a constructive role of noises. Such effects cannot be recorded without noise. Among various effects, the noise induced escape plays a special role, underlying two famous phenomena: the stochastic resonance and the

(35)

1.5. NOISE INDUCED EFFECTS 23 resonant activation. Those effects can be modelled by a system described by the Langevin equation (see Chap. 1.1.2) with the deterministic force, originating from the potential and a stochastic force, see Eq. (1.17). For a static potential, in the absence of noise, the particle deterministically slides towards the closest local minimum of the potential. The presence of the noise allows the particle to surmount the potential barrier, and in particular the particle can jump over the barrier due to the stochastic force.

Noise induced escape

Noise-induced escape can be observed already in the absence of any deterministic forces. The simplest model demonstrating the noise induced escape is a motion of a free particle in a bounded domain. In a one-dimensional case, it can be an interval restricted by two

absorbing boundaries located at_{±L. The particle’s motion is described by the Langevin}

equation (1.16). A particle starts atx0 (−L < x0 < L) and can move freely, until it

crosses one of the absorbing boundaries, when it is removed from the system and the first

passage timeτ is recorded

τ = min_{{t > 0 : x(0) = x}0 and |x(t)| > L} . (1.57)

The first passage time provides an information how much time is needed to leave the domain of motion for the first time. Another commonly investigated quantity is the

mean first passage time (MFPT) _{hτi, which is the average of first passage times,}

cal-culated for an ensemble of independent particles. The mean first passage time for the one-dimensional system described by Eq. (1.16) can be calculated analytically for the

interval[−L, L] where both boundaries are absorbing. For x(0) = 0 and escape induced

by anα-stable Lévy-type noise with the stability index α, see Eq. (1.40), the formula for

the MFPT [84] is hτi = L σ α 1 Γ(1 + α). (1.58)

The noise induced escape can be also observed in the presence of deterministic force, originating from the potential. In that case, the system can be described by the Langevin

(36)

equation (1.41). In such a case, forα = 2, the formula for the MFPT is hτ(x0 → L)i = 1 σ2 Z L x0 dx0exp V (x 0₎ σ2 Z x0 −L exp −V (x 00₎ σ2 dx00, (1.59)

see [41]. Moreover the left boundary can be shifted to_{−∞ or eliminated by a potential}

acting as a reflecting boundary. Examination of the resonant activation (RA) and stochas-tic resonance (SR) is related to studies of the MFPT in time dependent potentials.

Equa-tion (1.59) can be simplified using the saddle point approximaEqua-tion to_{hτi ≈ exp}∆V_σ2 ,

where∆V is the depth of the potential well.

Resonant activation

The simplest model of the resonant activation [8, 9] is a one-dimensional interval with the potential barrier, dichotomously switching between two linear configurations, char-acterized by various heights, see Fig. 1.4. The barrier switching process is Markovian

0 L H− H+ V (x ) x

Figure 1.4: The basic setup for the resonant activation. Potential dichotomously switches

between two configurations, characterized by various heightsH−andH+.

dichotomous [46] and the switching rate isγ. An analogous effect is observed for the

deterministic modulation of the partial barrier [85]. Similarly like in the escape from a bounded domain, the efficiency of the resonant activation is measured by the mean first passage time (MFPT). If the barrier is dichotomously modulated, the MFPT becomes a

(37)

1.5. NOISE INDUCED EFFECTS 25

non-monotonic function of the switching rateγ. Therefore, the resonant activation is the

effect resulting in the optimal escape kinetics over a time-dependent potential barrier, as measured by the MFPT, see Fig. 1.5 and [8]. Its efficiency originates from the interplay of the barrier modulation and noise.

−0.6 −0.3 0 0.3 0.6 0.9 1.2 1.5 −6 −4 −2 0 2 4 6 log 10 hτ i log10γ H−= 4, H+= 8 H−= 0, H+= 8 H−=−8, H+ = 8

Figure 1.5: Example of the optimal escape kinetics in a one-dimensional escape problem over modulated potential barrier [86]. The resonant activation can be captured as an

optimal (minimal) mean first passage time_hτi.

Stochastic resonance

Stochastic resonance [5, 87] is the celebrated noise induced effect, resulting in the op-timal input/output synchronization. The archetypal model of the SR is a motion of the overdamped particle in a periodically modulated double well potential. The Langevin equation describing the system is given by Eq. (1.17) with the time-dependent potential

V (x, t) =−a 2x 2₊ b 4x 4_{+ Ax sin(Ωt),} _(1.60) see Fig. 1.6.

In order to quantify the SR, there are various measures indicating the presence of the resonant effect. First of all, one can rely on the intuitive definition of the SR. The best input/output synchronization is recorder when during one period of the potential modulation, a particle can jump forth and back. Therefore, the SR is expected to take

(38)

−5 0 5 10 15 20 25 30 −3 −2 −1 0 1 2 3 V (x, t) x t = 0 t = π/2 t = 2π

Figure 1.6: The generic setup for the stochastic resonance. A double well potential (see Eq. (1.60)) is periodically modulated, changing the relative depths of the local minima.

place when

TΩ ≈ 2Te, (1.61)

whereTΩis the modulation period of the potential barrier, andTeis the escape time from

one of the potential’s minima. In classical examples, the modulation periodTΩ is fixed,

while the escape time Te is tuned by changing the noise intensity. Other measures of

stochastic resonance can be derived from power spectra, for example the spectral power amplification and signal-to-noise ratio [5, 7, 35]. The power spectrum is defined as a Fourier transform of the autocorrelation function

S(ω) =

Z ∞

−∞

e−iωτ_{hx(t + τ)x(t)idτ,} (1.62)

while the signal-to-noise ratio (SNR), following the definition in [5] is

SN R = 2 h lim∆ω→0 RΩ+∆ω Ω−∆ω S(ω)dω i SN(Ω) , (1.63)

where theSN(Ω) is a background at the modulation frequency Ω, while the numerator

gives the power carried by the signal. Fig. 1.7 presents sample power spectrum and SNR for Eq. (1.41) with the time-dependent potential (1.60), see [86].

(39)

1.5. NOISE INDUCED EFFECTS 27 100 101 102 103 104 105 106 107 10−2 ₁₀−1 ₁₀0 ₁₀1 ₁₀2 ₁₀3 P (ω ) ω −2 0 2 4 6 8 10 12 14 0 1 2 3 4 5 6 7 8 9 10 10 log 10 S N R [dB] σ

Figure 1.7: Sample power-spectra based measures of the stochastic resonance [86]. The left figure shows the power spectrum (see Eq. (1.62)), while the right column presents the signal-to-noise ratio (see Eq. (1.63)).

To detect the stochastic resonance, one can also use the probability of a given number of transitions [88,89] or the area under “first peak” of the residence time distribution [90]. These measures are further discussed in Chap. 2.5.

(40)

(41)

Chapter 2 Noise induced effects

The Langevin equation (see Chap. 1.1.2) can be generalized by replacing Gaussian white

noise with a more general α-stable Lévy type noise, containing the Gaussian noise as

a special case. A lot of such generalizations have been explored in great details. The most common form of the Langevin equation used in inspection of noise induced effects is a one-dimensional equation. Therefore, another generalization that can be done is to consider the Langevin equation in a multidimensional domain. The presence of a

mul-tidimensionalα-stable noise in the Langevin equation requires special care, due to the

appearance of the spectral measure (see Chap. 1.2.2).

The following sections discuss the noise induced effects observed in a multidimen-sional systems, where both the state variable and noise are multivariate. Throughout the

following chapters, the noise term in the Langevin equation is an α-stable Lévy type

noise, unless stated otherwise. We begin (Chap. 2.1) with a noise induced motion of a free particle in bounded domains, namely a 2D disk and a 3D sphere, with further gener-alizations into hyperspheres and hypercubes (Chap. 2.2). An analysis of a noise induced motion of a free particle in the bounded domains is performed mainly in terms of the mean first passage time, revealing its non-monotonic behaviour. Chapter 2.3 investigates

the behaviour of a test particle that is subjected to a two-dimensionalα-stable noise in a

bivariate single-well potential. It is shown that if the potential is steep enough, it produces stationary states. Next, we consider a motion of a random walker in time-dependent

(42)

tentials. First, we study a system with a dichotomously switching potential. Cooperative action of the noise and barrier modulation can result in optimal escape kinetics, namely a resonant activation (Chap. 2.4). Finally, we comment on a possibility of observing a stochastic resonance in a bivariate double-well potential (Chap. 2.5).

(43)

2.1. ESCAPE FROM HYPERSPHERES 31

2.1 Escape from hyperspheres

Escape in one dimension

A noise induced escape from the bounded domains is one of seminal problems studied in statistical physics. Canonical problem is the one-dimensional diffusion, induced by the Gaussian white noise [41, 91], eventually extended into Gaussian [92] and

non-Markovian regimes [93, 94], including α-stable noises. A bounded domain introduces

constraints on the particle movement. In the one dimension, this can be, for example,

a finite interval[−L, L] with absorbing boundaries. Such a setup introduces boundary

conditions that require special attention. Since trajectories ofα-stable motions are

dis-continuous forα < 2, the boundary conditions cannot be local [94–96].

The motion of a free particle subjected to the action of anα-stable noise can be

de-scribed using the overdamped Langevin equation dx

dt = σζα(t), (2.1)

wherex(t) is a particle position and ζα(t) is an α-stable white noise, which is a formal

time-derivative of theα-stable motion ζα(t) =

dLα(t)

dt . Lα(t) is a stochastic process with

independent increments, distributed according to anα-stable density, see Chap. 1.2.1.

Equation (2.1) can be associated with the Smoluchowski equation

∂p(x, t_|x0, 0)

∂t = σ

α∂αp(x, t|x0, 0)

∂|x|α , (2.2)

describing the evolution of the probability density p(x, t|x0, 0) of finding a particle at

positionx at time t, assuming the beginning of motion at position x0 at time0. ∂

α

∂|x|α is the

fractional Riesz-Weil derivative, see Chap. 1.3 and [73, 74]. For theα-stable noise with

the stability indexα < 2, and motion confined to a finite interval [−L, L], the equation

(2.2) must be associated with non-local boundary conditions [96, 97], namely

(44)

This is caused by a discontinuity of the random walker trajectory, which is responsible for the already mentioned non-local boundary conditions and leapovers [98, 99]. The mean

first passage time (MFPT) [100] can be calculated analytically for the interval [_{−L, L]}

where both boundaries are absorbing [84, 96]. Forx(0) = 0 the formula for the MFPT

[84] is hτi = L σ α 1 Γ(1 + α), (2.4)

which forα = 2 yields

hτi = L

2

2σ2. (2.5)

The MFPT might have a non-monotonic dependence on the stability index α. This is

caused by an interplay of the α-stable noise and non-local boundary conditions, which

results in a complex behaviour of the system. Analytical analysis reveals the non-triviality

of MFPT as a function of the stability indexα, see below.

Equation (2.4) can be differentiated with respect to the stability indexα

∂_hτi ∂α = L σ α _log L σ − ψ(1 + α) Γ(1 + α) , (2.6)

whereψ(z) is the polygamma function. Equation (2.6) can be used to calculate value of

the ratio L_σ for which the maximum of the mean first passage time _{hτi is reached (for a}

fixedα). Small values of the ratio L

σ result in maximum value of MFPT for low values of

the stability indexα ≈ 0. With the increasing value of L

σ, the maximum moves towards

α = 2, reaching it for L

σ ≈ 2.51. Intermediate values of

L

σ result in a non-monotonic

dependence of MFPT on the stability indexα. The interval length is fixed, therefore, for a

given interval half-widthL, one can always find such a noise intensity σ, which produces

the non-monotonic dependence of MFPT as a function of the stability indexα. Mean first

passage times for various values of the ratio L_σ are depicted in Fig. 2.1. Non-monotonic

(45)

2.1. ESCAPE FROM HYPERSPHERES 33 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 hτ i α L/σ = 0.2 L/σ = 0.7 L/σ = 1.3 L/σ = 1.7 L/σ = 2.1 L/σ = 2.5 L/σ = 2.8

Figure 2.1: Mean first passage times for the noise induced escapes from the

one-dimensional interval[_{−L, L] as a function of the stability index α for various values of the}

L

σ ratio. Low values of

L

σ produce maximum value of MFPT forα ≈ 0. With increasing

value of L_σ, the maximum moves towardsα = 2, reaching it for L

σ ≈ 2.51.

Escape in two dimensions

Analogously to the 1D case, the random walker position is described by the Langevin equation

dx

dt = σζα(t), (2.7)

where ζ_α(t) is the multidimensional α-stable noise. Eq. (2.7) is a straightforward

gener-alization of the one-dimensional Eq. (2.1). In two (and more) dimensions the multivariate α-stable noise is defined by the chosen spectral measure (see Chap. 1.2.2). Various spec-tral measures alter the escape kinetics in the system [101]. Thus, additional attention is required when moving to a multidimensional domain. The choice of a spectral measure implies an implicit choice of stochastic kinetics. Consequently, the equation that describes the system is also altered. The most common choices of spectral measures are the uniform continuous spectral measure and uniform discrete spectral measure, concentrated on the intersections of the unit circle with the axes. The uniform continuous spectral measure

(continuous onS2) leads to a sphericalα-stable motion, so called spherical Lévy flights,

while the discrete spectral measure leads to Cartesian Lévy flights [102]. Uniform con-tinuous spectral measure results in the following Smoluchowski equation (compare Eq.

(46)

(2.2)) ∂p(x, t|x0, 0) ∂t =−σ α (−∆)α/2 p(x, t|x0, 0), (2.8)

where_−(−∆)α/2 _{is the fractional Riesz-Weil derivative, see Chap. 1.3 and [73, 74]. The}

symmetric discrete spectral measure, concentrated on intersections of the unit circle with the axes, leads to the Smoluchowski equation of the form

∂p(x, t_|x0, 0) ∂t = σ α ∂α ∂|x|α + ∂α ∂|y|α p(x, t|x0, 0), (2.9)

where fractional derivatives are defined in Chap. 1.3.

Similar to 1D, trajectories of bivariateα-stable motion are discontinuous, therefore,

the boundary conditions must be non-local. In the simplest case, boundary conditions can by deducted directly from the one-dimensional case

p(x, t_|x0, 0) = 0 for |x| > R. (2.10)

Boundary conditions (2.10) clearly impose a domain of motion to be a two-dimensional

disk of radius R. This is the simplest choice, matching the symmetry of the uniform

continuous spectral measure. Because no direction is favored, the symmetry of the domain of motion can be later utilized to simplify the analysis of the first passage times.

Equations (2.8) and (2.9) are in general case fractional differential equations with not known analytical solutions. Solutions for the mentioned equations can be estimated by the Monte Carlo methods [43, 60, 103, 104]. Estimation is done by simulating the random walker trajectories, following the Langevin dynamics (2.7). From a large number

of trajectories x(t) the probability density p(x, t_|x0, 0) following Eqs (2.8) or (2.9) can

be estimated. Moreover, from the ensemble of trajectories, other properties of the process can be directly estimated, for example the mean first passage time. Within the model, a

random walker starts at the center of the disk x(0) = x0 = (0, 0) and performs a

two-dimensional Lévy flights (noise inducedα-stable motion). When the trajectory reaches or

(47)

2.1. ESCAPE FROM HYPERSPHERES 35 passage time is recorded. The first passage time is defined as the shortest time after which the random walker leaves the system, namely

τ = min_{{t > 0 : x(0) = 0 and |x(t)| > R} .} (2.11)

The whole edge of the disk and its exterior are absorbing, as implied by boundary condi-tions, see Eq. (2.10). Sample trajectories of the described process x(t) are presented in Fig. 2.2. The left column presents spherical Lévy flights (see Eq. (2.8)), while the right column depicts the Cartesian Lévy flights, see Eq. (2.9) and [102]. The top row presents

trajectories induced by noise with the stability indexα = 1.0, while the bottom row shows

sample trajectories for the stability indexα = 1.9. For α = 1.0 (top row) the pronounced

long jumps of the random walker are clearly visible, for both spherical and the Cartesian Lévy flights. For the Cartesian Lévy flights, the jump length density is a product of two 1D Cauchy densities, while for the spherical Lévy flights, jumps are distributed according to bivariate isotropic Cauchy density, see also Chapter 1.2.2.

The Monte Carlo simulations reveal that the mean first passage time in two dimensions

can also be a non-monotonic function of the stability indexα. Therefore, this property of

theα-stable noise induced escape transfers directly from the one-dimensional system to

higher number of dimensions. It is possible to find such a scale parameterσ for which the

mean first passage time is a non-monotonic function of the stability indexα. Similar to the

one-dimensional setup, the position of the maximum of the MFPT shifts with the change

of the scale parameter σ. Estimated values of the MFPT for the uniform continuous

spectral measure are presented in the left column of Fig. 2.3. For large enough value of

the scale parameterσ the mean first passage is a monotonically decreasing function of the

stability indexα. With the decrease of the scale parameter σ, non-monotonic behaviour

can be observed. Finally, for low values of the scale parameter, the MFPT is an increasing

function of the stability indexα, see the left column of Fig. 2.3.

For the uniform continuous spectral measure in the limit ofα = 2, a spherically

(48)

(a) Spherical LF, α= 1.0 (b) Cartesian LF, α= 1.0

(c) Spherical LF, α= 1.9 (d) Cartesian LF, α= 1.9

Figure 2.2: Sample trajectories of bivariate Lévy flights confined in a circular domain. The left column demonstrates spherical Lévy flights, while the right column shows the Cartesian Lévy flights. The top row presents trajectories induced by noise with the

sta-bility index α = 1.0, while the bottom row depicts sample trajectories for the stability

indexα = 1.9. Once the trajectory reaches the edge of the domain of motion (or crosses

it), it is removed from the system. Long jumps of Lévy flights are clearly visible for the

stability indexα = 1.0 (top row) for both spectral measures. Every picture contains three

independent trajectories.

time _{hτi can be calculated directly from the bivariate backward Smoluchowski}

equa-tion [41, 100]

∆_{hτi = −} 1

(49)

2.1. ESCAPE FROM HYPERSPHERES 37 0 0.5 1 1.5 2 0 0.5 1 1.5 2 hτ i α σ = 2.0 σ = 2.5 σ = 3.0 σ = 4.0 σ = 5.0 0 0.5 1 1.5 2 0 0.5 1 1.5 2 hτ i α σ = 2.0 σ = 2.5 σ = 3.0 σ = 4.0 σ = 5.0 σ = 6.0

Figure 2.3: Mean first passage time for the 2D disk for spherical Lévy flights (left column)

and the Cartesian Lévy flights (right column), induced by various type ofα-stable noises.

The whole edge of the disc is absorbing. Points represent values estimated by the Monte Carlo methods. Solid lines depict exact values of MFPT calculated from Eqs. (2.29) and

(2.32). The disk radius isR = 4.

with two-dimensional Laplace operator∆ = ∂2

∂x2 + ∂ 2

∂y2 and boundary conditionhτi(R) =

0. Eq. (2.12) can be rewritten in the polar coordinates

∂2_{hτi(r, ϕ)} ∂r2 + 1 r ∂_{hτi(r, ϕ)} ∂r + 1 r2 ∂2_{hτi(r, ϕ)} ∂ϕ2 =− 1 σ2. (2.13)

Because of the spherical symmetry of the system, the MFPT _{hτi does not depend on}

the angleϕ. Eq. (2.13) transforms into a univariate, second order, ordinary differential

equation d2_hτi(r) dr2 + 1 r dhτi(r) dr =− 1 σ2 (2.14)

with the solution

hτi(r) = 1

4σ2 R

2

− r2_,

(2.15)

wherer is the initial distance from the center of the disk. For r = 0 the solution (2.15)

reduces to

hτi = R

2

4σ2. (2.16)

A direct comparison between results for 1D (Eq. (2.5)) and 2D (Eq. (2.16)) reveals that

the mean first passage time for the disk of radius R is two times smaller, than for the

Noise induced effects in multivariate systems driven by Lévy noise

Jagiellonian University

Faculty of Physics, Astronomy and Applied Computer Science