# Thermalization of Lévy flights: Path-wise picture in 2D

## Abstract

We analyze two-dimensional (2D) random systems driven by a symmetric Lévy stable noise which, under the sole influence of external (force) potentials , asymptotically set down at Boltzmann-type thermal equilibria. Such behavior is excluded within standard ramifications of the Langevin approach to Lévy flights. There, the action of a conservative force field stands for an explicit reason for the emergence of an asymptotic invariant probability density function (pdf). However the latter cannot be represented in the Boltzmann form 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 . That is explicitly encoded in non-symmetric jump transition rates of the master equation for the pdf . We prescribe a priori the target pdf in the Boltzmann form 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 as time 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’s 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 pre-defined target pdf. We pay special attention to the response of the 2D Cauchy noise to an exemplary locally periodic ”potential landscape” .

## I 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 e.g. 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 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), (1) and (2).

As a necessary pre-requisite 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 field in the form , where stands for the normalized white noise: , and the mass parameter is scaled away. The corresponding Fokker-Planck equation for the probability density function reads

(1) |

and, in the confining regime, is known to enforce the existence of an asymptotic invariant pdf, as , in the explicit Boltzmann form , where .

By means of a standard substitution , (1), the Fokker-Planck equation can be transformed into a generalized diffusion equation for an auxiliary function . This Schrödinger-type equation (no imaginary unit ) reads

(2) |

where and .

By reintroducing a normalization constant (divide and multiply by a suitable number in the factorization formula for ), we can rewrite in the form , where while . Clearly, as goes to infinity. Moreover, we can rewrite the semigroup potential as follows: .

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, (3), disconnects the Langevin-modeled Fokker-Planck equation for any Lévy-stable noise

(3) |

from the very notion of the Boltzmann thermal equilibrium.

However, the thermal equilibrium notion remains a valid concept within an immediate Lévy transcript of the semigroup dynamics (2) (e.g. replace by ):

(4) |

see e.g. (4); (5), where we assume that asymptotically sets down at a square root of a well defined pdf . The semigroup potential follows from the compatibility condition:

(5) |

In this particular context, while adopting a multiplicative decomposition of the time-dependent pdf :

(6) |

a novel fractional generalization of the Fokker-Planck equation governing the time evolution of has been introduced in Refs. (6)-(8), see also (5); (9); (10), 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, (6)-(10), has the familiar master equation form, presently reproduced in the explicit 2D form:

(7) |

Here, anticipating the effectiveness of numerical routines to be described in below, from the start we impose cut-offs upon the size of jumps to be accounted for during the simulations (large jump sizes remain statistically significant if view of the involved Lévy distribution)

(8) |

and the Lévy measure in is given by

(9) |

It is the quantity which has an interpretation of the jump transition rate from the point to another point . The potential function can be chose quite arbitrarily. However, we need to secure a normalization of the target pdf . We note that becomes a genuine stationary solution of Eq. (7) once we let and .

## Ii Gillespie’s algorithm: Fine tuning in 2D.

The Gillespie’s algorithm has been originally constructed to give a dynamical picture of a finite chain of chemical reactions, (11); (12). There, switches between different chemical reaction channels can be re-interpreted 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 for a while) all jumps form a fixed point to any point in the set . Since we aim at numerical simulations, it is obvious that all admissible jumps cannot be computer-generated. Therefore, we must restrict considerations to a large number of suitably selected ”representative” jumps of a truncated jump process ruled by the truncated Lévy distribution of jumps.

Basic tenets of the modified Gillespie’s algorithm, fine tuned to account for Lévy jumps in have been described in detail elsewhere, (13). Presently, we shall give a brief outline of the algorithm that is capable to account for Lévy flights in . We mimic previously devised 1D steps, (13), while adopting them to the 2D situation.

Eq. (7) can be rewritten as follows

(10) |

where

(11) |

The algorithm outline reads as follows:

(i) Set time and the point of origin for jumps, .

(ii) Create the set of all admissible jumps from to , compatible with the transition rate .

(iii) Evaluate

(12) |

and .

(iv) Using a random number generator draw from a uniform distribution

(v) Employing above and and the identities

(13) |

where i , find and corresponding to the ”transition channel” .

(vi) Draw a new number from a uniform distribution.

(vii) Reset time label s where : .

(viii) Reset to a new location , considered as a reference point for the iterative procedure.

(ix) Return to step (ii) and repeat the procedure anew.

Remark: All precautions respected in 1D need to be respected in 2D as well (c.f. Comment 1 in (13)). The following jump size bounds (integration boundaries) were adopted for
exemplary numerical procedures to be described in below: i , provided one keeps in memory our convention not to
reproduce the index in anymore.

## Iii 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 with , when .

### iii.1 Harmonic confinement: Gaussian target

Let us prescribe an asymptotic invariant pdf to be in a 2D Gaussian form:

(14) |

As an exemplary source of random noise we assume the Cauchy driver, i.e. 2D Lévy-stable noise with the stability index . Accordingly, the jump transition rate reads:

(15) |

To generate sample paths of the process, we need first to evaluate integrals of (15) over suitable integration volumes. If the integration volume comprises which are close to the jump size boundaries , 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 are far away from integrals (15) are amenable to standard evaluation methods (like e.g. Simposon’s one). To be more explicit in this respect, let . If , then

(16) |

where

(17) |

Numerical routines were written in terms of C-codes, (14).

In Fig. 1 we have depicted the statistical data inferred from 100 000 trajectories, for three running time instants , , , together with the asymptotic expression (14). The right-hand-side column depicts projections of those data upon the the plane. A substanitial increase of the analyzed enesemble data (like e.g. , or ) is merely a matter of the simulation time span and adds nothing inspiring to the obtained behavior.

It is clear that the surrogate pdf evolution consistently goes towards an invariant asymptotic pdf (compare subfigures c) and d)). We note a lowering and flattening of the maximum around , in consistency with the ultimate target outcome, whose height is as follows from . Visually accessible in homogeneities of circular shapes in OXY projections b) i c), are a consequence of still relatively low number (100 000) of sample paths data and approximations involved in evaluating involved integrals. The shape of OXY projection a) is a consequence of the initial data choice. Our problem has a radial symmetry. Therefore in Fig. 2 we depict a projection of the trajectory induced data upon the plane. The projection shows as well a consistent convergence towards the target pdf.

An additional control method for the path-wise inferred pdf evolution, addresses the time evolution and an asymptotic behavior of the pdf moments and . Here is the mean distance of points of a trajectory form the origin at the running time instant , while is a mean square distance from . In view of

(18) |

the dynamics should set down at , while this of at . Fig. 3 depicts the evolution of and , inferred from the simulated sample of jumping paths.

The observed convergence and validates the number generator choice, we have used to arrive at sample jumping paths.

### iii.2 Logarithmic confinement: 2D Cauchy target.

We consider the target pdf in the 2D Cauchy form:

(19) |

Like previously, we take the Cauchy driver, , as a reference Lévy stable noise. Accordingly:

(20) |

Proceeding like in the Gaussian case, for small and , Eq. (20) can be approximated by

(21) |

where terms linear in i were preserved. The result can be analytically integrated term after term by employing Eqs. (17)..

Accumulated trajectory data have been analyzed to produce Fig. 4, where a surrogate pdf evolution is displayed. Fig. 5 reproduces the projection of the obtained pdf data. If compared with the 1D case analyzed in Ref. (13) a convergence to Boltzmannian equilibrium (target pdf) is substantially slowed down.

### iii.3 Locally periodic confinement in .

We consider target pdf whose Boltzmanian exponent is locally periodic (within a finite rectangle) and almost entirely localized within a finite spatial area due to harmonically confining tails of the potential:

(22) |

The normalization constant actually reads .

Subsequently adopted numerical integration routines heavily rely on the experience gained during our previous case studies. In Fig. 6 we report the surrogate pdf evolution, inferred from 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 .

We are aware of the fact that the number of trajectories may be considered as too small and not sufficiently representative sample. Our tentative 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 are very conspicuous if the number of involved trajectories gets significantly lowered by imposing constraints (like e.g. various spatial projections). All l 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 Fig. 7, where we have considered projections of the surrogate pdf data upon planes , and at time . We have set them in a direct comparison with respective target pdf (12) data.

## Iv Outlook

We have taken into consideration jump-type processes which can not 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 non-standard 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 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, (11); (12), 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 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 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 be random potentials of Ref. (7).

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 e.g. 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 re-analyze the path-wise data we have collected and stored in the trajectory generation process.

### References

- H. Risken, The Fokker-Plack Equation, (Springer-Verlag, Berlin, 1989)
- S. Jespersen, R. Metzler and H. C. Fogedby, Phys. Rev. E59, 2736, (1999)
- I. Eliazar and J. Klafter, J. Stat. Phys. 119, 165, (2005)
- P. Garbaczewski and R. Olkiewicz, J. Math. Phys. 40, 1057, (1999)
- P. Garbaczewski and V. A. Stephanovich, Phys. Rev. E 80, 031113, (2009)
- D. Brockmann and I. M. Sokolov, Chemical Physics, 284, 409, (2002)
- D. Brockmann and T. Geisel, Phys. Rev. Lett. 90, 170601, (2003)
- V. V. Belik and D. Brockmann, New. J. Phys. 9, 54, (2007)
- P. Garbaczewski and V. A. Stephanovich, Physica A 389, 4410, (2010)
- P. Garbaczewski and V. Stephanovich, Phys. Rev. E 84, 011142 (2011)
- T. D. Gillespie, J. Phys. Chemistry, 8 (25), 2340, (1977)
- T. D. Gillespie, J. Comput. Physics, 22 (4), 403, (1976)
- M. Żaba, P. Garbaczewski and V. Stephanovich, Lévy flights in confining environments: Random paths and their statistics, arXiv: 1209.5882
- C-codes available upon request