# Thermalization of Lévy Flights: Path-Wise Picture in 2D

Authors: Mariusz Żaba; Piotr Garbaczewski
Journal: International Journal of Statistical Mechanics
 (2013)

---

## Abstract

We analyze two-dimensional (2D) random systems driven by a symmetric Lévy stable noise which in the presence of confining potentials may asymptotically set down at Boltzmann-type thermal equilibria. In view of the Eliazar-Klafter no-go statement, such dynamical behavior is plainly incompatible with the standard Langevin modeling of Lévy flights. No explicit path-wise description has been so far devised for the thermally equilibrating random motion we address, and its formulation is the principal goal of the present work. To this end we prescribe a priori the target pdfρ∗ in the Boltzmann form ~exp[−Φ(x)] and next select the Lévy noise (e.g., its Lévy measure) of interest. To reconstruct random paths of the underlying stochastic process we resort to numerical methods. We create a suitably modified version of the time honored Gillespie algorithm, originally invented
in the chemical kinetics context. A statistical analysis of generated sample trajectories allows us to infer a surrogate pdf ρ(x,t) dynamics which sets down at a predefined target, in consistency with the associated kinetic (master) equation.


---

## Body

## 1. Introduction

Various random processes in real physical systems admit a simplified description based on stochastic differential equations. Then, there is a routine passage procedure from microscopic random variables to macroscopic (statistical ensemble, mean field) data, like, for example, the time evolution of an associated probability density function (pdf) which is a solution of a deterministic transport equation. A paradigm example is so-called the Langevin modeling of diffusion-type and jump-type processes. The presumed microscopic model of random dynamics is provided by the Langevin (stochastic) equation, which additively decomposes into a (Newtonian by origin) drift and purely random (perturbing noise) term. Its direct consequence is the Fokker-Planck equation for an associated probability density function (pdf); confer [1] for a discussion of the Brownian motion and [2, 3] for that of Lévy flights in external forces. We note that the Lévy-Langevin formulation results in the space-fractional Fokker-Planck equation.The subject of our further discussion is two-dimensional (2D) random systems driven by a symmetric Lévy stable noise which, under the sole influence of external (force) potentialsΦ(x), asymptotically set down at Boltzmann-type thermal equilibria. Such behavior is excluded within standard ramifications of the Langevin approach to Lévy flights, where the action of a conservative force field ~-∇Φ(x) stands for an explicit reason for the emergence of an asymptotic invariant probability density function (pdf). The latter cannot be represented in the Boltzmann form ρ*(x)~exp[-Φ(x)], and the thermal equilibrium concept appears to be alien to Langevin-modeled Lévy flights.In the present paper we address the response of Lévy noise not to an external conservative force field, but directly to its potentialΦ(x). That is explicitly encoded in nonsymmetric jump transition rates of the master equation for the pdf ρ(x,t).We prescribe a priori the target pdfρ* in the Boltzmann form ~exp[-Φ(x)] and next select the Lévy noise of interest. Given suitable initial data, this allows to infer a reliable path-wise approximation to a true (albeit analytically beyond the reach) solution of the pertinent master equation, with the property ρ(x,t)→ρ*(x) as time t goes to infinity.No explicit path-wise description has been so far devised for such thermally equilibrating random motion. To reconstruct random paths of the underlying stochastic process we resort to numerical methods, where long jumps of the Lévy stable processes are statistically significant but are truncated to become amenable to simulation procedures. We create a suitably modified version of the time honored Gillespie algorithm, originally invented in the chemical kinetics context.A statistical analysis of generated sample trajectories allows us to infer a surrogate pdf dynamics which consistently sets down at a predefined target pdf. We pay special attention to the response of the 2D Cauchy noise to an exemplary locally periodic “potential landscape”Φ(x),x∈R2.As a necessary prerequisite for our further discussion, let us discuss a transformation of the Fokker-Planck equation into the Schrödinger-type (generalized diffusion) equation, often employed in the theoretical framework of the Brownian motion. Here, the Langevin equation, the induced Fokker-Planck equation, and its Schrödinger-type image are dynamically equivalent and describe the same diffusion-type process. This is not the case if one turns over to jump-type processes.For clarity of arguments, let us consider the Langevin equation for a one-dimensional diffusion process in an external conservative force fieldF(x)=-dV(x)/dx in the form dx/dt=F(x)+2νb(t), where b(t) stands for the normalized white noise: 〈b(t)〉=0, 〈b(t′)b(t)〉=δ(t-t′) and the mass parameter is scaled away. The corresponding Fokker-Planck equation for the probability density function ρ(x,t) reads
(1)∂tρ=νΔρ-∇(Fρ)
and, in the confining regime, is known to enforce the existence of an asymptotic invariant pdf, ρ(x,t)→ρ*(x) as t→∞, in the explicit Boltzmann form exp(-Φ/2), where Φ(x)=V(x)/ν.By means of a standard substitutionρ(x,t)=ψ(x,t)exp[-V(x)/2ν] [1], the Fokker-Planck equation can be transformed into a generalized diffusion equation for an auxiliary function ψ(x,t). This Schrödinger-type equation (no imaginary unit i) reads
(2)∂tψ=νΔψ-𝒱(x)ψ,
where 𝒱(x)=(1/2)((F2/2ν)+∇F) and F=F(x).By reintroducing a normalization constant (divide and multiply by a suitable numberZ1/2 in the factorization formula for ρ(x,t)), we can rewrite ρ(x,t) in the form ρ(x,t)=Ψ(x,t)ρ*1/2(x), where ρ*1/2=Z-1/2exp(-Φ/2) while Ψ=Z1/2ψ. Clearly, Ψ(x,t)→ρ*1/2(x) as t goes to infinity. Moreover, we can rewrite the semigroup potential as follows: 𝒱=(νΔρ*1/2)/ρ*1/2.The transformation of (1) into (2) cannot be adopted to Lévy jump-type processes, where the Langevin and Schrödinger-type (semigroup) modeling are known to be incompatible. Moreover, the Eliazar-Klafter no go statement [4] disconnects the Langevin-modeled Fokker-Planck equation for any Lévy-stable noise
(3)x˙=b(x)+Aμ(t)⟹∂tρ=-∇(b·ρ)-λ|Δ|μ/2ρ.
From the very notion of the Boltzmann thermal equilibrium. (We note in passing that an ample literature is available on various aspects of the Lévy-Langevin random motion (3), specifically on the associated asymptotic invariant pdfs whose inverse polynomial decay may be much steeper than that of any Lévy-stable pdf; see, for example, [5].)However, the thermal equilibrium notion remains a valid concept within an immediate Lévy transcript of the semigroup dynamics (2) (e.g., replace νΔ by -λ|Δ|μ/2):
(4)∂tΨ=-λ|Δ|μ/2Ψ-𝒱Ψ;
see, for example, [6, 7], where we assume that Ψ(x,t) asymptotically sets down at a square root ρ*1/2(x) of a well-defined pdf ρ*. The semigroup potential 𝒱(x) follows from the compatibility condition:
(5)𝒱=-λ|Δ|μ/2ρ*1/2ρ*1/2.In this particular context, while adopting a multiplicative decomposition of the time-dependent pdfρ(x,t)→ρ*(x),
(6)ρ(x,t)=Ψ(x,t)ρ*1/2(x),
a novel fractional generalization of the Fokker-Planck equation governing the time evolution of ρ(x,t) has been introduced in [8, 9, 16] (see also [7, 10, 11]) to handle systems that are randomized by symmetric Lévy-stable drivers and may asymptotically set down at Boltzmann-type equilibria under the influence of external potentials (thus not Newtonian forces anymore).The pertinent Fokker-Planck type equation, whose origin has been discussed before in a number of papers [8–11, 16], has the familiar master equation form, presently reproduced in the explicit 2D form:
(7)∂tρ(x,y)=∬A[wϕ(x,y∣u,v)ρ(u,v)-wϕ(u,v∣x,y)ρ(x,y)]νμ(x,y,dx,dy).
Here, anticipating the effectiveness of numerical routines to be described below, from the start we impose cut-offs upon the size of jumps to be accounted for during the simulations. A contribution of larger jump sizes nonetheless remains statistically significant. This  happens because the   unconstrained   Lévy distribution is used in sampling procedures. We have verified this property before in the 1D considerations of [12],
(8)A={ɛ2y(x,y)∈ℝ2;ɛ1x⩽|x-u|⩽ɛ2x∧ɛ1y⩽|y-v|⩽ɛ2yℝ2},wϕ(x,y∣u,v)=exp[(Φ(u,v)-Φ(x,y))2],
and the Lévy measure νμ(x,y,dx,dy) in ℝ2 is given by
(9)νμ(x,y,dx,dy)=2μΓ((2+μ)/2)π|Γ(-μ/2)|1(x2+y2)(2+μ)/2dxdy=ν~μ(x,y)dxdy,μ∈(0,2).
It is the quantity wϕ(x,y∣u,v)ν~μ(x,y) which has an interpretation of the jump transition rate from the point (u,v)∈R2 to another point (x,y)∈R2. The potential function Φ(x,y) can be chosen quite arbitrarily. However, we need to secure a L1(R2) normalization of the target pdf ρ*(x,y)~exp(-Φ(x,y)). We note that exp(-Φ(x,y)) becomes a genuine stationary solution of (7) once we let ɛ1x,y→0 and ɛ2x,y→∞.


## 2. Gillespie’s Algorithm: Fine Tuning in 2D

Gillespie’s algorithm has been originally constructed to give a dynamical picture of a finite chain of chemical reactions [13, 14]. There, switches between different chemical reaction channels can be reinterpreted as jumps between points in a finite state space. Since the number of allowed chemical channels is finite, a serious modification of an original algorithm must be created to account for very large (virtually infinite) state space we need to consider in connection with Lévy flights. As an example, a jump process analog of chemical reaction channels could comprise (take R1 for a while) all jumps form a fixed point x0 to any point in the set [x0-ϵ2,x0-ϵ1]∪[x0+ϵ1,x0+ϵ2].Since numerical simulations impose intrinsic lower and upper bounds upon the jump size, it is obvious that what we actually implement is a truncated jump process, with the “representative” truncated Lévy distribution of jumps.Basic tenets of the modified Gillespie’s algorithm, fine-tuned to account for Lévy jumps inR1, have been described in detail elsewhere [12]. Presently, we shall give a brief outline of the algorithm that is capable of accounting for Lévy flights in R2. We mimic previously devised 1D steps [12], while adopting them to the 2D situation. Namely, (7) can be rewritten as follows:
(10)∂tρ(x,y)=∬A′Cμ[wϕ(x,y∣ξ+x,η+y)ρ(ξ+x,η+y)-wϕ(ξ+x,η+y∣x,y)ρ(x,y)]dξdη,
where
(11)A′={(ξ,η)∈ℝ2;ɛ1x⩽|ξ|⩽ɛ2x∧ɛ1y⩽|η|⩽ɛ2y},Cμ=2μsin(πμ/2)[Γ((μ+2)/2)]2π2(ξ2+η2)(μ+2)/2.
The algorithm outline reads as follows.(i)
Set timet=0 and the point of origin for jumps, (x0,y0)∈R2.(ii)
Create the set of all admissible jumps from(x0,y0) to (x0+ξ,y0+η), compatible with the transition rate wϕ(x,y∣x+ξ,y+η).(iii)
Evaluate(12)W1(x0,y0)=∫-ɛ2x-ɛ1x∫-ɛ2y-ɛ1yCμwϕ(ξ+x0,η+y0∣x0,y0)dξdη,W2(x0,y0)=∫-ɛ2x-ɛ1x∫ɛ1yɛ2yCμwϕ(ξ+x0,η+y0∣x0,y0)dξdη,W3(x0,y0)=∫ɛ1xɛ2x∫-ɛ2y-ɛ1yCμwϕ(ξ+x0,η+y0∣x0,y0)dξdη,W4(x0,y0)=∫ɛ1xɛ2x∫ɛ1yɛ2yCμwϕ(ξ+x0,η+y0∣x0,y0)dξdη,
and W(x0,y0)=W1(x0,y0)+W2(x0,y0)+W3(x0,y0)+W4(x0,y0).(iv)
Using a random number generator drawpx,py∈[0,1] from a uniform distribution.(v)
Employing abovepx and py and the identities
(13)∫-ɛ2xbx(∫-ɛ2y-ɛ1yCμwϕ(ξ+x0,η+y0∣x0,y0)dη+∫ɛ1yɛ2yCμwϕ(ξ+x0,η+y0∣x0,y0)dη)dξ=pxW(x0,y0),px<W12(x0,y0)W(x0,y0),W12(x0,y0)+∫ɛ1xbx(∫-ɛ2y-ɛ1yCμwϕ(ξ+x0,η+y0∣x0,y0)dη+∫ɛ1yɛ2yCμwϕ(ξ+x0,η+y0∣x0,y0)dη)dξ=pxW(x0,y0),px⩾W12(x0,y0)W(x0,y0),∫-ɛ2yby(∫-ɛ2x-ɛ1xCμwϕ(ξ+x0,η+y0∣x0,y0)dξ+∫ɛ1xɛ2xCμwϕ(ξ+x0,η+y0∣x0,y0)dξ)dη=pyW(x0,y0),py<W13(x0,y0)W(x0,y0),W13(x0,y0)+∫ɛ1yby(∫-ɛ2x-ɛ1xCμwϕ(ξ+x0,η+y0∣x0,y0)dξ+∫ɛ1xɛ2xCμwϕ(ξ+x0,η+y0∣x0,y0)dξ)dη=pyW(x0,y0),py⩾W13(x0,y0)W(x0,y0),
where W12(x0,y0)=W1(x0,y0)+W2(x0,y0) and W13(x0,y0)=W1(x0,y0)+W3(x0,y0), find bx and by corresponding to the “transition channel” (x0,y0)→(bx,by).(vi)
Draw a new numberq∈(0,1) from a uniform distribution.(vii)
Reset time labelt=t+Δt, where Δt=-lnq/W(x0,y0).(viii)
Reset(x0,y0) to a new location (x0≐x0+bx,y0≐y0+by), considered as a reference point for the iterative procedure.(ix)
Return to step (ii), and repeat the procedure anew.Remark 1.
All precautions respected in 1D need to be respected in 2D as well (cf. comment  1 in [12]). The following jump size bounds (integration boundaries) were adopted for exemplary numerical procedures to be described, ɛ1x=ɛ1y=ɛ1=0.001 i ɛ2x=ɛ2y=ɛ2=1, provided that one keeps in memory our convention not to reproduce the x, y index in ɛ1,2x,y anymore.


## 3. Statistics of Random Paths in 2D: Case Studies of the PDF Evolution

Our main purpose in the present section is to analyze the response of free Lévy noise to generic external potentials, that has been previously found to be encoded in jump transition rates. Now we address the same problem in a path-wise fashion. We shall faithfully follow the outlined random path generation procedure. Once suitable path ensemble data are collected, we shall verify whether statistical (ensemble) features of generated random trajectories are compatible with those predicted by the master equation (7). This includes a control of an asymptotic behavior ρ(x,t)→ρ*(x) with x∈R2, when t→∞.


### 3.1. Harmonic Confinement: Gaussian Target

Let us prescribe an asymptotic invariant pdfρ* to be in a 2D Gaussian form:
(14)ρ*(x,y)=1πe-x2-y2.
As an exemplary source of random noise we assume the Cauchy driver, that is, 2D Lévy-stable noise with the stability index μ=1. Accordingly, the jump transition rate reads
(15)C1wϕ(ξ+x,η+y∣x,y)=12πe-ξ2/2-xξ-η2-yη(ξ2+η2)3/2.To generate sample paths of the process, we need first to evaluate integrals of (15) over suitable integration volumes. If the integration volume comprises x, y which are close to the jump size boundaries ±ɛ1, one develops a numerator of an expression into Taylor series and keeps terms up to the quadratic one. Errors induced by such approximation procedure are marginal. On the other hand, if x,y are far away from ±ɛ1 integrals (15) are amenable to standard evaluation methods (like, e.g., Simposon’s one). To be more explicit in this respect, let ɛ12=0.05. If |a|,|b|,|c|,|d|⩽ɛ12; then
(16)∫ab∫cde-ξ2/2-xξ-η2-yη(ξ2+η2)3/2dξdη≈∫ab∫cd((ξ2+η2)-3/2(1-xξ+((x2-1)/2)ξ2)×(1-yη+((y2-1)/2)η2)×(ξ2+η2)-3/2)dξdη≈∑i=16‍𝕀i,
where
(17)𝕀1=∫ab∫cddξdη(ξ2+η2)3/2=∫ab(dξ2ξ2+d2-cξ2ξ2+c2)dξ=-b2+d2bd+a2+d2ad+b2+c2bc-a2+c2ac,𝕀2=-∫ab∫cdxξdξdη(ξ2+η2)3/2=xln((d+b2+d2)(c+a2+c2)(c+b2+c2)(d+a2+d2)),𝕀3=-∫ab∫cdyηdξdη(ξ2+η2)3/2=yln((b+b2+d2)(a+a2+c2)(a+a2+d2)(b+b2+c2)),𝕀4=∫ab∫cdxyξηdξdη(ξ2+η2)3/2=xy(a2+d2+b2+c2-b2+d2-a2+c2),𝕀5=∫ab∫cdx2-12ξ2dξdη(ξ2+η2)3/2=x2-12(dln(b+b2+d2a+a2+d2)-cln(b+b2+c2a+a2+c2)),𝕀6=∫ab∫cdy2-12η2dξdη(ξ2+η2)3/2=y2-12(bln(d+b2+d2c+b2+c2)-aln(d+a2+d2c+a2+c2)).
Numerical routines were written in terms of C-codes [15].In Figure1 we have depicted the statistical data inferred from 100 000 trajectories, for three running time instants t=0.2, t=0.7, and t=15, together with the asymptotic expression (14). The right-hand-side column depicts projections of those data upon the OXY plane. A substantial increase of the analyzed ensemble data (like, e.g., 300000, 500000, or 106) is merely a matter of the simulation time span and adds nothing inspiring to the obtained behavior.Gaussian target; statistics of100000 random paths: surrogate pdf evolution and its OXY projection, displayed in the vertical order for t=0.2, t=0.7, and t=15. The bottom subfigure refer to the asymptotic pdf (14). All trajectories were started form the origin (0,0).
(a)
(b)
(c)
(d)It is clear that the surrogate pdf evolution consistently goes towards an invariant asymptotic pdf. We note a lowering and flattening of the maximum around(0,0), in consistency with the ultimate target outcome, whose height directly follows from ρ*(0,0)=1/π. Visually accessible inhomogeneities of circular shapes in OXY projections, are a consequence of still relatively low number (100 000) of sample paths data and approximations involved in evaluating involved integrals. The shape of the top (initial) OXY projection is a consequence of the initial data choice.Our problem has a radial symmetry. Therefore in Figure2 we depict a projection of the trajectory induced data upon the y=0 plane. The projection shows as well a consistent convergence towards the target pdf.Figure 2
We have displayed they=0 projection of the previous pdf data (100000 trajectories) for running time instants (a) t=0.2, (b) t=0.7, and (c) t=15; (d) a projection of the target pdf is depicted.An additional control method for the path-wise inferred pdf evolution addresses the time evolution and an asymptotic behavior of the pdf moments〈R(t)〉 and 〈R2(t)〉. Here 〈R(t)〉 is the mean distance of (x,y) points of a trajectory from the origin (0,0) at the running time instant t, while 〈R2(t)〉 is a mean square distance from (0,0). In view of
(18)〈R〉as=1π∬ℝ2x2+y2e-x2-y2=π2≈0.886,〈R2〉as=1π∬ℝ2(x2+y2)e-x2-y2=1,
the 〈R(t)〉 dynamics should set down at π/2, while this is of 〈R2(t)〉 at 1. Figure 3 depicts the evolution of 〈R(t)〉 and 〈R2(t)〉, inferred from the simulated sample of 100000 jumping paths.Time evolution of (a) the first moment (averager=x2+y2 value), (b) second moment (mean value of r2). All 100000 trajectories originate from (0,0).
(a)
(b)The observed convergence〈R(t)〉→1/π and 〈R2(t)〉→1 validates the number generator choice we have used to arrive at sample jumping paths.


### 3.2. Logarithmic Confinement: 2D Cauchy Target

We consider the target pdfρ* in the 2D Cauchy form:
(19)ρ*(x,y)=12π1(1+x2+y2)3/2.
Like previously, we take the Cauchy driver, μ=1, as a reference Lévy stable noise. Accordingly,
(20)C1wϕ(ξ+x,η+y∣x,y)=12π1(ξ2+η2)3/2(1+x2+y21+(x+ξ)2+(y+η)2)3/4.
Proceeding like in the Gaussian case, for small ξ and η, (20) can be approximated by
(21)12π1(ξ2+η2)3/2(1-3x2(1+x2+y2)ξ-3y2(1+x2+y2)η),
where terms linear in ξ i η were preserved. The result can be analytically integrated term after term by employing (17).Accumulated trajectory data have been analyzed to produce Figure4, where a surrogate pdf evolution is displayed. Figure 5 reproduces the y=0 projection of the obtained pdf data. If compared with the 1D case analyzed in [12] a convergence to Boltzmannian equilibrium (target pdf) is substantially slowed down.Cauchy target; surrogate pdf evolution inferred from100000 sample trajectories and the OXY projection for running times (vertical order) t=0.2, t=3, and t=500. The bottom subfigure refers to the target Cauchy pdf (19). In the course of simulations all trajectories were started from (0,0).
(a)
(b)
(c)
(d)Figure 5
The projection of the previous pdf data on they=0 plane, at time instants t=0.2, t=3, and t=500. The latter curve actually appears to coincide with the target pdf (asymptote) projection, at least within adopted scales.


### 3.3. Locally Periodic Confinement inR2

We consider target pdf whose Boltzmannian exponent is locally periodic (within a finite rectangle) andalmostentirely localized within a finite spatial area due to harmonically confining tails of the potential:
(22)ρ*(x,y)=1C{e-3sin2(2πx)sin2(2πy),|x|⩽2i|y|⩽2,e-x2-y2+8,|x|>2i|y|>2,e-x2+4,|x|>2i|y|⩽2,e-y2+4,|x|⩽2i|y|>2.
The normalization constant C actually reads C=13.5921.Subsequently adopted numerical integration routines heavily rely on the experience gained during our previous case studies. In Figure6 we report the surrogate pdf evolution, inferred from 100000 sample trajectories. A convergence rate to the asymptotic (target) pdf is satisfactory, although a reasonable agreement with the target data has been achieved for relatively large running time values, here t=300.Locally periodic target; surrogate pdf evolution and itsOXY projection inferred from 100000 trajectories at running time instants (vertical order) t=0.2, t=3 and t=300. The bottom subfigure depicts the asymptotic pdf (22). All trajectories have been started from (0,0).
(a)
(b)
(c)
(d)We are aware of the fact that the number of100000 trajectories may be considered as a too small and not sufficiently representative sample. Our tentative 300000 paths data do not show significant qualitative changes in the obtained evolution picture.We should mention that there are significant statistical fluctuations to be kept under control. They become very conspicuous if the number of involved trajectories gets significantly lowered by imposing constraints (like, e.g., various spatial projections). All trajectory data, after being gathered, are safely stored in the computer memory. Therefore we can get access to any conceivable and more detailed statistical picture of what is going on, even if the outcome is hampered by significant random deviations from the reference (target) pdf data.A sample of such fluctuating data is provided in Figure7, where we have considered projections of the surrogate pdf data upon planes y=0, y=1/4, and y=x at time t=300. We have set them in a direct comparison with respective target pdf (12) data.Figure 7
Projections of the surrogate pdf data att=300 upon planes (vertical order) y=0, y=1/4 and y=x. An “asymptote” refers to the target pdf projection.


## 3.1. Harmonic Confinement: Gaussian Target

Let us prescribe an asymptotic invariant pdfρ* to be in a 2D Gaussian form:
(14)ρ*(x,y)=1πe-x2-y2.
As an exemplary source of random noise we assume the Cauchy driver, that is, 2D Lévy-stable noise with the stability index μ=1. Accordingly, the jump transition rate reads
(15)C1wϕ(ξ+x,η+y∣x,y)=12πe-ξ2/2-xξ-η2-yη(ξ2+η2)3/2.To generate sample paths of the process, we need first to evaluate integrals of (15) over suitable integration volumes. If the integration volume comprises x, y which are close to the jump size boundaries ±ɛ1, one develops a numerator of an expression into Taylor series and keeps terms up to the quadratic one. Errors induced by such approximation procedure are marginal. On the other hand, if x,y are far away from ±ɛ1 integrals (15) are amenable to standard evaluation methods (like, e.g., Simposon’s one). To be more explicit in this respect, let ɛ12=0.05. If |a|,|b|,|c|,|d|⩽ɛ12; then
(16)∫ab∫cde-ξ2/2-xξ-η2-yη(ξ2+η2)3/2dξdη≈∫ab∫cd((ξ2+η2)-3/2(1-xξ+((x2-1)/2)ξ2)×(1-yη+((y2-1)/2)η2)×(ξ2+η2)-3/2)dξdη≈∑i=16‍𝕀i,
where
(17)𝕀1=∫ab∫cddξdη(ξ2+η2)3/2=∫ab(dξ2ξ2+d2-cξ2ξ2+c2)dξ=-b2+d2bd+a2+d2ad+b2+c2bc-a2+c2ac,𝕀2=-∫ab∫cdxξdξdη(ξ2+η2)3/2=xln((d+b2+d2)(c+a2+c2)(c+b2+c2)(d+a2+d2)),𝕀3=-∫ab∫cdyηdξdη(ξ2+η2)3/2=yln((b+b2+d2)(a+a2+c2)(a+a2+d2)(b+b2+c2)),𝕀4=∫ab∫cdxyξηdξdη(ξ2+η2)3/2=xy(a2+d2+b2+c2-b2+d2-a2+c2),𝕀5=∫ab∫cdx2-12ξ2dξdη(ξ2+η2)3/2=x2-12(dln(b+b2+d2a+a2+d2)-cln(b+b2+c2a+a2+c2)),𝕀6=∫ab∫cdy2-12η2dξdη(ξ2+η2)3/2=y2-12(bln(d+b2+d2c+b2+c2)-aln(d+a2+d2c+a2+c2)).
Numerical routines were written in terms of C-codes [15].In Figure1 we have depicted the statistical data inferred from 100 000 trajectories, for three running time instants t=0.2, t=0.7, and t=15, together with the asymptotic expression (14). The right-hand-side column depicts projections of those data upon the OXY plane. A substantial increase of the analyzed ensemble data (like, e.g., 300000, 500000, or 106) is merely a matter of the simulation time span and adds nothing inspiring to the obtained behavior.Gaussian target; statistics of100000 random paths: surrogate pdf evolution and its OXY projection, displayed in the vertical order for t=0.2, t=0.7, and t=15. The bottom subfigure refer to the asymptotic pdf (14). All trajectories were started form the origin (0,0).
(a)
(b)
(c)
(d)It is clear that the surrogate pdf evolution consistently goes towards an invariant asymptotic pdf. We note a lowering and flattening of the maximum around(0,0), in consistency with the ultimate target outcome, whose height directly follows from ρ*(0,0)=1/π. Visually accessible inhomogeneities of circular shapes in OXY projections, are a consequence of still relatively low number (100 000) of sample paths data and approximations involved in evaluating involved integrals. The shape of the top (initial) OXY projection is a consequence of the initial data choice.Our problem has a radial symmetry. Therefore in Figure2 we depict a projection of the trajectory induced data upon the y=0 plane. The projection shows as well a consistent convergence towards the target pdf.Figure 2
We have displayed they=0 projection of the previous pdf data (100000 trajectories) for running time instants (a) t=0.2, (b) t=0.7, and (c) t=15; (d) a projection of the target pdf is depicted.An additional control method for the path-wise inferred pdf evolution addresses the time evolution and an asymptotic behavior of the pdf moments〈R(t)〉 and 〈R2(t)〉. Here 〈R(t)〉 is the mean distance of (x,y) points of a trajectory from the origin (0,0) at the running time instant t, while 〈R2(t)〉 is a mean square distance from (0,0). In view of
(18)〈R〉as=1π∬ℝ2x2+y2e-x2-y2=π2≈0.886,〈R2〉as=1π∬ℝ2(x2+y2)e-x2-y2=1,
the 〈R(t)〉 dynamics should set down at π/2, while this is of 〈R2(t)〉 at 1. Figure 3 depicts the evolution of 〈R(t)〉 and 〈R2(t)〉, inferred from the simulated sample of 100000 jumping paths.Time evolution of (a) the first moment (averager=x2+y2 value), (b) second moment (mean value of r2). All 100000 trajectories originate from (0,0).
(a)
(b)The observed convergence〈R(t)〉→1/π and 〈R2(t)〉→1 validates the number generator choice we have used to arrive at sample jumping paths.


## 3.2. Logarithmic Confinement: 2D Cauchy Target

We consider the target pdfρ* in the 2D Cauchy form:
(19)ρ*(x,y)=12π1(1+x2+y2)3/2.
Like previously, we take the Cauchy driver, μ=1, as a reference Lévy stable noise. Accordingly,
(20)C1wϕ(ξ+x,η+y∣x,y)=12π1(ξ2+η2)3/2(1+x2+y21+(x+ξ)2+(y+η)2)3/4.
Proceeding like in the Gaussian case, for small ξ and η, (20) can be approximated by
(21)12π1(ξ2+η2)3/2(1-3x2(1+x2+y2)ξ-3y2(1+x2+y2)η),
where terms linear in ξ i η were preserved. The result can be analytically integrated term after term by employing (17).Accumulated trajectory data have been analyzed to produce Figure4, where a surrogate pdf evolution is displayed. Figure 5 reproduces the y=0 projection of the obtained pdf data. If compared with the 1D case analyzed in [12] a convergence to Boltzmannian equilibrium (target pdf) is substantially slowed down.Cauchy target; surrogate pdf evolution inferred from100000 sample trajectories and the OXY projection for running times (vertical order) t=0.2, t=3, and t=500. The bottom subfigure refers to the target Cauchy pdf (19). In the course of simulations all trajectories were started from (0,0).
(a)
(b)
(c)
(d)Figure 5
The projection of the previous pdf data on they=0 plane, at time instants t=0.2, t=3, and t=500. The latter curve actually appears to coincide with the target pdf (asymptote) projection, at least within adopted scales.


## 3.3. Locally Periodic Confinement inR2

We consider target pdf whose Boltzmannian exponent is locally periodic (within a finite rectangle) andalmostentirely localized within a finite spatial area due to harmonically confining tails of the potential:
(22)ρ*(x,y)=1C{e-3sin2(2πx)sin2(2πy),|x|⩽2i|y|⩽2,e-x2-y2+8,|x|>2i|y|>2,e-x2+4,|x|>2i|y|⩽2,e-y2+4,|x|⩽2i|y|>2.
The normalization constant C actually reads C=13.5921.Subsequently adopted numerical integration routines heavily rely on the experience gained during our previous case studies. In Figure6 we report the surrogate pdf evolution, inferred from 100000 sample trajectories. A convergence rate to the asymptotic (target) pdf is satisfactory, although a reasonable agreement with the target data has been achieved for relatively large running time values, here t=300.Locally periodic target; surrogate pdf evolution and itsOXY projection inferred from 100000 trajectories at running time instants (vertical order) t=0.2, t=3 and t=300. The bottom subfigure depicts the asymptotic pdf (22). All trajectories have been started from (0,0).
(a)
(b)
(c)
(d)We are aware of the fact that the number of100000 trajectories may be considered as a too small and not sufficiently representative sample. Our tentative 300000 paths data do not show significant qualitative changes in the obtained evolution picture.We should mention that there are significant statistical fluctuations to be kept under control. They become very conspicuous if the number of involved trajectories gets significantly lowered by imposing constraints (like, e.g., various spatial projections). All trajectory data, after being gathered, are safely stored in the computer memory. Therefore we can get access to any conceivable and more detailed statistical picture of what is going on, even if the outcome is hampered by significant random deviations from the reference (target) pdf data.A sample of such fluctuating data is provided in Figure7, where we have considered projections of the surrogate pdf data upon planes y=0, y=1/4, and y=x at time t=300. We have set them in a direct comparison with respective target pdf (12) data.Figure 7
Projections of the surrogate pdf data att=300 upon planes (vertical order) y=0, y=1/4 and y=x. An “asymptote” refers to the target pdf projection.


## 4. Outlook

We have taken into consideration jump-type processes which cannot be handled by standard stochastic differential equation methods (e.g., the Langevin modeling, where a conspicuous motion “tendency” quantified by an additive drift term can be unambiguously isolated from the noise contribution). Existing popular algorithms cannot provide a direct numerical simulation of sample paths of such nonstandard processes. In the present paper, we have proposed a working method to generate stochastic trajectories (sample paths) of a random jump-type process that avoids any reference to a stochastic differential equation. An additional gain of that procedure is that we are in fact capable of reliably approximating the time evolutionρ(x,t) of a true (typically not available in a closed analytic form) solution of the master equation.To this end we have modified the Gillespie algorithm, [13, 14], normally devised for sample paths generation if the transition rates refer to a finite number of states of a system. The essence of our modification is that we take into account the continuum of possible transition rates, thereby changing the finite sums in the original Gillespie algorithm into integrals. The corresponding procedures for stochastic trajectories generation have been changed accordingly.In other words, here we are able (i) to extract the background sample paths of a jump process and (ii) to infer a reliable approximation of an actual (analytically unavailable) solution of the master equation (7)-(8). We emphasize once more here that we have focused on those jump-type processes that cannot be modeled by any stochastic differential equation of the Langevin type.Although heavy-tailed Lévy stable drivers were involved in the present considerations, we have clearly confirmed that a large variety of stationary target distributions is dynamically accessible for each particularμ∈(0,2) Lévy driver choice. That variety comprises not only the standard Gaussian pdf, casually discussed in relation to the Brownian motion (e.g., the Wiener process), but the whole non-Gaussian family, associated with the Lévy stable conceptual imagery.Among heavy-tailed distributions, we have paid attention to the Cauchy pdf which can stand for an asymptotic target for anyμ∈(0,2) driver, provided a steering environment (e.g., “potential landscape”) is properly devised. In turn, the Cauchy driver, while excited in a proper environment, may lead to an asymptotic pdf with an arbitrarily large number of moments, the previously mentioned Gaussian case being included.An example of the locally periodic environment has been considered as a toy model for more realistic physical systems. Our major hunch are strongly inhomogeneous “potential landscapes,” modeled by relatively smooth potentials. We note that a radically extreme variant could comprise random potentials of [16].In connection with the master equation which was our departure point let us stress that, even if various mean field data are available in experimentally realizable systems, it is of vital interest to gain some knowledge about the microscopic dynamics (random paths) realized by the random system under consideration. The detailed analysis of sample path data with a focus on their specific features, like, for example, ergodicity, mixing, or lack of those properties, deserves a separate analysis. These goals can be achieved as well within the present simulation framework. It suffices to reanalyze the path-wise data we have collected and stored in the trajectory generation process.



---

*Source: 738345-2013-10-03.xml*