Trajectory statistics of confined Lévy flights and Boltzmann-type equilibria

Trajectory statistics of confined Lévy flights and Boltzmann-type equilibria


We analyze a specific class of random systems that are driven by a symmetric Lévy stable noise, where Langevin representation is absent. In view of the Lévy noise sensitivity to environmental inhomogeneities, the pertinent random motion asymptotically sets down at the Boltzmann-type equilibrium, represented by a probability density function (pdf) . Here, we infer pdf based on numerical path-wise simulation of the underlying jump-type process. A priori given data are jump transition rates entering the master equation for and its target pdf . To simulate the above processes, we construct a suitable modification of the Gillespie algorithm, originally invented in the chemical kinetics context. We exemplified our algorithm simulating different jump-type processes and discuss the dynamics of real physical systems where it can be useful.

I Introduction

Despite many attempts to pin down the essential features of dynamics and relaxation in random systems, the problem is still far from its complete solution. It turns out that Lévy flight models are adequate for the description of different random systems ranging from the motion of defects in disordered solids to the dynamics of assets in stock markets, see, e.g. (1). They are especially useful to model the random systems on the semi-phenomenological, mesoscopic level, when the (often unknown) details of their microscopic random behavior are substituted by a properly tailored (e.g. based on experimental data) (Gaussian or Lévy) noise. Paradoxically, in disordered solids, the noise can promote order and organization, switching them between different equilibrium states. The latter situation emerges in disordered ferroelectrics, where the fluctuations of order parameter (spontaneous polarization) give rise to self-localization of charge carrier, generating a fluctuon, an analog of well-known polaron in disordered substance (2). These fluctuons make a substantial contribution into conductivity and optical properties of disordered dielectrics.

Many random processes admit a description based on stochastic differential equations. In such case there is a routine passage procedure from microscopic random variables to macroscopic (statistical ensemble) data. The latter are encoded in the time evolution of an associated probability density function (pdf) which is a solution of a deterministic transport equation. A paradigmatic example is the so-called Langevin modeling of diffusion-type and jump-type processes. The presumed microscopic model of the dynamics in external force fields is provided by the Langevin (stochastic) equation whose direct consequence is the Fokker-Planck equation, (3) and (4). We note that in case of jump-type processes the familiar Laplacian (Wiener noise generator) needs to be replaced by a suitable pseudo-differential operator (fractional Laplacian, in case of a symmetric Lévy-stable noise).

We pay a particular attention to jump-type processes which are omnipresent in Nature (see (5) and references therein). Their characterization is primarily provided by jump transition rates between different states of the system under consideration. In the present paper, our major focus is on a specific class of random systems which are plainly incompatible with a straightforward Langevin modeling of jump-type processes and, as such, are seldom addressed in the literature.

To this end we depart from the concept, coined in an isolated publication (6), of Lévy flights-driven models of disorder that, while at equilibrium, do obey detailed balance. The corresponding research line has been effectively initiated in Refs. (7)-(9). It has next been expanded in various directions, with a special emphasis on so-called Lévy-Schrödinger semigroup reformulation of the original probability density function (pdf) dynamics, (18); (19); (8) and (10)-(15), c.f. also (17); (18); (19). We note in passing that the familiar Fokker-Planck equation can be equally well reformulated in terms of the Schrödinger semigroup and this property is universally valid in the standard theory of Brownian motion, (3); (17). Its generalization to Lévy flights is neither immediate nor obvious. It is often considered in the prohibitive vein following (20); (23).

In fact, in relation to Lévy flights, a novel fractional generalization of the Fokker-Planck equation has been introduced in Refs. (7)-(9) to handle systems that are randomized by symmetric Lévy-stable drivers. In this case, contrary to the popular lore about properties of (Langevin-based) Lévy processes (c.f. Refs. (20)-(22) and (23)), the pertinent random systems are allowed to relax to (thermal) equilibrium states of a standard Boltzmann-Gibbs form.

The underlying jump-type processes, in the stationary (equilibrium) regime, respect the principle of detailed balance by construction (15); (16). Their distinctive feature, if compared with the standard Langevin modeling of Lévy flights, is that they have a built-in response not to external forces but rather to external force potentials. These potentials are interpreted to form confining ”potential landscapes” that are specific to the environment. Lévy jump-type processes appear to be particularly sensitive to environmental inhomogeneities, (7); (14).

We note in passing that confining Lévy flights resolve themselves to the problem of truncated distributions, which can be addressed in many ways. beside a simple cut-off, by considering fast falling tails. The random walk theory predominantly employs to the master equation as the major tool to quantify random motion. Inhomogeneous transition rates have been used in the literature before and the emergence of non-trivial target distributions has been demonstrated, c.f. (24).

Lévy flights are pure jump (jump-type) processes. Therefore, it seems useful to indicate that various model realizations of standard jump processes (jump size is bounded from below and above) can be thermalized by means of a specific scenario of an energy exchange with the thermostat. It is based on the principle of detailed balance. We have discussed this issue in some detail before (15) along with an extension of this conceptual framework to Lévy-stable processes. Not to reproduce easily available arguments of past publications, we shall be very rudimentary in our motivations, see however (16).

We quantify a probability density evolution, compatible with a jump-type process on (this limitation may in principle be lifted in favor of ), in terms of the master equation:


where and are, respectively, the lower and upper bounds of jump size and


is the jump transition rate from to . We stress that is a non-symmetric function of and .

An implicit Boltzmann-type weighting involves a square root of a target pdf and accounts for the a priori prescribed ”potential landscape” whose confining features affect the jump-type process. What matters is a relative impact of a confinement strength of (level of attraction, see Ref. (9)) upon jumps of the size , both at the point of origin and that of destination . In principle, may be an arbitrary function that secures a normalization of . In this case, the resultant pdf is a stationary solution of the transport equation (1) with unbounded jump length, e.g. and .

We note that the presence of lower and upper bounds of the jump size , that are necessary for an implementation of numerical algorithms, enforce a truncation of the jump-type process (without any cutoffs) to a standard jump process. The transition rates of the latter, however, are ruled by Lévy measures of symmetric Lévy stable noises with . A lower bound for the jump size is usually removed while evaluating the corresponding integrals in the sense of their Cauchy principal values. An upper bound is less innocent and its effects need to be controlled by long tailed pdfs which stands for a distinctive feature of Lévy flights, see a discussion of Lévy stable limits of step processes in Ref. (10). There is also pertinent discussion of a long time behavior of (unconfined, e.g. free) truncated Lévy flights in Ref. (28).

In contrast to procedures based on the Langevin modeling of Lévy flights in external force fields, (4); (21); (22), there is no known path-wise approach underlying the transport equation (1). With no direct access to sample trajectories of the stochastic process in question, a method must be devised to generate random paths directly from jump transition rates (2). The additional requirement here is that we set a priori a ”potential landscape” for a chosen jump-type (symmetric Lévy stable) noise driver.

The outline of the paper is as follows. First we describe our modification of the Gillespie algorithm which entails a numerical generation of random paths for the dynamics determined by Eqs (1) and (2). Next the statistics of random paths is addressed and various accumulated data are analyzed with a focus on inherent compatibility issues.

We analyze generic (Cauchy, quadratic Cauchy) and non-generic (Gaussian and locally periodic) examples of target pdfs for the jumping dynamics. Random paths are generated in conjunction with representative Lévy stable drivers, like e.g. those indexed by . Their qualitative typicality is emphasized.

Statistical data, acquired from our modification of Gillespie algorithm, have been employed to generate the dynamical patterns of behavior . Effectively, that entails a path-wise reconstruction of the solution of the master equation (1), (2), or - in case this solution it is available prior to trajectory simulations - to verify a compatibility of the transport (master) equation (1), (2) and its underlying path-wise representation. All that comes from the predefined knowledge of the target pdf and non-symmetric (biased) jump transition rates.

Ii Modified Gillespie algorithm.

Here we adopt (29) (and properly adjust to handle Lévy flights) basic tenets of so-called Gillespie’s algorithm (25); (26). Originally, this algorithm had been devised to simulate random properties of coupled chemical reactions. The advantage of the algorithm is that it permits to generate random trajectories of the corresponding stochastic process directly from its (jump) transition rates, with no need for any stochastic differential equation and/or its explicit solution. We emphasize that this feature of Gillespie’s algorithm is vitally important, since Langevin modeling is not operational in our framework.

We rewrite Eq. (1) in the form ()


To construct a reliable path generating algorithm consistent with Eq. (3) we first note that chemical reaction channels in the original Gillespie’s algorithm may be re-interpreted as jumps from one spatial point to another, like transition channels in the spatial jump process. An obvious provision is that the set of possible chemical reaction channels is discrete, while we are interested in continuous set of all admissible jumps from a chosen point of origin . With a genuine computer simulation in mind, we must respect standard numerical assistance limitations. Surely we cannot admit all conceivable jump sizes. As well, the number of destination points, even if potentially enormous, must remain finite for any fixed point of origin.

Our modified version of the Gillespie’s algorithm, appropriate for handling of spatial jumps is as follows (30):

  1. Set time and the point of origin .

  2. Create the set of all admissible jumps from to that is compatible with the transition rate .

  3. Evaluate


    and .

  4. Using a random number generator draw from a uniform distribution.

  5. Using above and identities


    find corresponding to the ”transition channel” .

  6. Draw a new number from a uniform distribution.

  7. Reset time label where .

  8. Reset to a new value .

  9. Return to the step 2 and repeat the procedure anew.

Iii Statistics of random paths for different jump-type processes

In this section, we select the jump-type processes, which best demonstrate the capabilities of our algorithm of random paths generation. Once suitable path ensemble data are collected, we ca n reconstruct the pdf dynamics and next verify whether statistical (ensemble) features of generated random trajectories are compatible with the asymptotic a priori imposed upon solutions of the master equation (1). That refers to a control of an asymptotic behavior when .

iii.1 Harmonic confinement

Let us first consider an asymptotic invariant (target) pdf in the Gaussian form:


The corresponding -family of transition rates reads


The C-codes for trajectory generating algorithm, (30), were employed to get the trajectory statistics data for Lévy drivers with .

Figure 1: Gaussian target: Time evolution of the pdf second moment for 75 000 trajectories. Inset visualize the oscillations smoothing in the asymptotic regime for ; figures near curves correspond to values.

Figure 2: Gaussian target: Time evolution of inferred from 75 000 trajectories: (left panel), (middle panel) and (right panel). All trajectories originate from , i.e. refer to the -type initial probability distribution.

The results of our numerical simulations are reported in Figs. 1 and 2. We note, that in Fig. 1 the second moment oscillates near its equilibrium value . A numerical convergence to is consistent with an analytic equilibrium value of the second moment of the chosen 6. The rate of this convergence is higher for larger . Clearly, for small the big jumps are frequent which enlarges the inferred time intervals in the Gillespie’s algorithm. Thus, the relaxation to equilibrium is slow. It gets faster for larger , when big jumps are rare and time intervals are generically very small.

Fig. 2 displays a probability density evolution, inferred from the ensemble statistics of trajectories. All of them have started form the same point . Although the data fidelity grows with the number of contributing paths, we have not found significant qualitative differences to justify a presentation of data for 100 000, 200 000, 250 000 and more trajectories. The relaxation time rate dependence on is clearly visible as well. It suffices to analyze differences between three curves for and/or . We observe a conspicuous lowering of their maxima with the growth of (take care of different scales on the vertical axes on Fig.2 panels). The simulated pdfs at are practically indistinguishable from an exact analytical asymptotic pdf 6. The convergence of towards appears to be relatively fast irrespective of the chosen -driver.

Although our reasoning is definitely path-wise and all data have been extracted from simulated trajectory ensembles, our focus was on the inferred dynamics of . We do not reproduce figures visualizing generic sample paths. However, we indicate their (paths) distinctive qualitative typicalities, if one would compare e.g. motion scenarios corresponding to stability indices , and , The structural impact (probability of occurrence) of larger against smaller jumps, even on the visual level, conforms with standard simulations of Lévy stable sample paths (with no forces or potentials involved), c.f. (31).

iii.2 Logarithmic confinement

Quadratic Cauchy target

Let us consider a long-tailed asymptotic pdf which is a special case of the one-parameter -family of equilibrium (Boltzmann-type) states, associated with a logarithmic potential , , see (15); (13); (12); (14) :


The transition rate 2 for any takes the form


Figure 3: Quadratic Cauchy target: Time evolution of the pdf second moment for 200 000 trajectories.

Figure 4: Quadratic Cauchy target: Time evolution of inferred from trajectories for (left panel), (middle panel) and (right panel). All trajectories are started from . Note scale differences on vertical axes.

Simulation results for this case are displayed in Figs. 3 and 4. If we compare Fig. 3 with Fig. 1 we see the existence of small oscillations in the asymptotic regime about the value . Those from Fig.1 were relatively small, while those on Fig. 3 are more noticeable. This is related to much slower decay of transition rates (9) (determined by slow-decaying asymptotic pdf (8)), as compared to those for Gaussian case (7).

The second moment of the present , (8), equals 1 and the convergence towards this value is clearly seen in Fig. 3. This convergence is much slower than in the Gaussian (harmonic confinement) case which is not a surprise: (7) and (9) indicate that the present rate of convergence should be logarithmically slower. Fig. 5, quite alike Fig. 2, convincingly demonstrates a convergence of to the asymptotic . For definitely large times around , and become practically indistinguishable. Similarly to the Gaussian case, the rate of convergence becomes larger with the growth of .

Figure 5: Cauchy target: Time evolution of pdf for and trajectories, starting from .

Cauchy target

Now we consider an asymptotic pdf of the form :


In this case, the transition rate from to reads


We consider Cauchy driver corresponding to . In Fig. 5, we report the time evolution of the statistically inferred for different time instants. An approach to the asymptotic pdf (10) is clearly seen, together with a convergence of a half-width to its asymptotic value 1. It can be shown that the same convergence pattern is valid for cumulative probability distribution which approaches the asymptotic function .

iii.3 Locally periodic confinement

To set firm grounds for future research it is instructive to study our model for more complicated forms of confining potentials. In view of their physical relevance, it is appealing to address an issue of confining (trapping) environments with a periodic spatial structure. This problem may be relevant to the random motion of defects in doped semiconductors and to phenomena like hopping conductivity. Here, we encounter a major difficulty with integrability of the Boltzmann-type weighting function as for periodic the corresponding integral is divergent. Periodicity and integrability can here be reconciled by means of locally periodic potentials that take a definite confining form (harmonic or polynomial) for larger values of . Let us consider the following asymptotic pdf


where is a normalization constant. The transition rate from to reads


where the potential has the following locally confining form

Figure 6: Time evolution of inferred from trajectories at . The data for and trajectories (not displayed) do not show qualitative differences.

Time evolution of the inferred pdf is reported in Fig. 6 for . All sample trajectories were started form which corresponds to the -type initial distribution. The probability density spreads out with time in conformity with the trapping (confining) properties of the locally periodic enclosure (environment or ”potential landscape”). For large running times t=400 the trajectory statistics produces data that are indistinguishable from those for the asymptotic pdf. We have checked that beginning from about 100 000 trajectories, further accumulation of the trajectories number like e.g. 200 000 (displayed) and 300 000 (not displayed) for the data statistics is inessential. In such cases the curves are almost the same, we merely improve a fidelity of the statistics.

Iv Conclusions.

If a random process does not admit the description in terms of a stochastic differential equation (e.g. Langevin modeling), its direct numerical simulation becomes impossible by means of existing popular algorithms. In the present paper, for the first time in the literature, we propose a working method to generate stochastic trajectories (sample paths) of a random jump-type process without resorting to any explicit (or numerical) solution of a stochastic differential equation.

To this end we have modified the Gillespie algorithm (25); (26), normally devised for sample paths generation, if transition rates occur between 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 procedure for stochastic trajectories generation has been changed accordingly.

In other words, here we ”extract” the background sample paths of a jump process, whose pdf is to obey the generalized Fokker-Planck dynamics 1, 2, . We emphasize once more here, that we have focused our attention on those background jump-type processes that cannot be modeled by any stochastic differential equation of the Langevin type. However, our ultimate result was a reconstruction of the pdf dynamics from the path-wise data. Thus, our modfication of the Gillespie’s algorithm serves two purposes: (i) getting access to path-wise data and (ii) reconstructing the pdf dynamics, compatible with the master equation in question, without resorting to other methods of its solution

Although heavy-tailed Lévy stable drivers were involved in the present considerations, we have clearly confirmed that an enormous variety of stationary target distributions is dynamically accessible in each particular case. That comprises not only a standard Gaussian pdf, discussed usually in relation to the Brownian motion like Wiener process. We can reach Cauchy pdf in terms of any driver as well, provided a steering environment is properly devised. In turn, the Cauchy driver in a proper environment may lead to an asymptotic pdf with a finite number of moments, the Gaussian case included ((15)).

An example of the locally periodic environment has been considered as a toy model for more realistic physical systems. We mention possible generalizations of our method to the Brownian motor concept (see, e.g., Ref. (32) for recent review) to include a non-Gaussian jumping component. In those systems it is the properly tailored periodic ”potential landscape” which enforces a conversion of a homogenous stochastic process (Brownian motion for reference) into the directed motion of particles at nanometer scales. It is closely related to the problem of so-called sorting in periodic potentials (33). Other problems to be addresses concern ultracold atoms in optical lattices subject to random potentials (34), which are promising not only from purely scientific point of view, but also for many technological applications.

We note that the theoretical description of the above mentioned topics relies essentially on the Langevin-like equation input. Our approach offers prospects for generalizations, where systems with non-Langevin response to external potentials may come into consideration, along with more traditional ones. What we actually need to implement our version of Gillespie’s algorithm is the knowledge of jump transition rates of those random systems only.

Our preliminary work (in progress) shows that an extension of our algorithm to higher dimensions is operational. In particular, the planar case is worth further analysis, possibly with more realistic ”potential landscapes”. The direct numerical modeling in three-dimensional case can also be important, for instance, in the investigations of the dynamics of amorphous materials. Theoretical studies of their physical properties have a longer history, starting in the eighties and even earlier (35); (36). These materials may comprise conventional glasses (see, e.g., (37)), so-called spin (35); (37) and orientational (for example, dipole) glasses (36) as well as virtually any disordered solid like doped semiconductors. In these substances, the random motion of defects can be regarded as real jump-type process to be modeled numerically. Moreover, if defects have only rotational degrees of freedom (like spins or electric/elastic dipoles), their dynamics can be well represented by some effective jump-type process. Such representation can be complementary to existing Monte-Carlo modeling techniques which for disordered systems are very computer intensive. For instance, we hope to tackle the issue of interplay between exponential and long-time (logarithmic or stretched-exponential) relaxation in these substances. Also, the 3D modification of our algorithm can be used in the simulation of switching processes in disordered ferroelectrics (38). Similar fictitious jump-type process can be assigned to the dynamics of inhomogeneous broadening of resonant lines (39). This broadening occurs in condensed matter and/or biological species, in a number of spectroscopic manifestations like the electron paramagnetic resonance, nuclear magnetic resonance, optical and neutron scattering methods. It arises due to random electric and magnetic fields, strains and other perturbations from defects in a substance containing the centers whose resonant transitions between energy levels are studied. Here also, the standard approaches cannot describe all the details of experiment and we hope that our algorithm can improve the situation.

Let us finally note that in the 3D problems admitting full spherical symmetry (like spin glasses with random exchange interactions), the problem becomes reducible to 1D case so that our algorithm might work in its present form.

Acknowledgement: P. G. would like to thank Professor I. M. Sokolov for focusing our attention on the Gillespie’s algorithm. All numerical simulations were completed by means of the facilities of the Platon - Science Services Platform of the Polish Pionier Network.


  1. H. S. Wio, R.R. Deza, J.M. Lopez, An Introduction to Stochastic Processes and Nonequilibrium Statistical Physics, (World Scientific Publishing, Singapore, 2012).
  2. V.A.Stephanovich, JETP Letters 65, 446 (1997).
  3. H. Risken, The Fokker-Plack Equation, (Springer-Verlag, Berlin, 1989)
  4. S. Jespersen, R. Metzler and H. C. Fogedby, Phys. Rev. E 59, 2736, (1999)
  5. R. Metzler, J. Klafter, J. Phys. A:Math. Gen. 37, R161, (2004)
  6. L. Chen and M. W. Deem, Phys. Rev. E 65, 011109, (2001)
  7. D. Brockmann and I. M. Sokolov, Chemical Physics, 284, 409, (2002)
  8. D. Brockmann and T. Geisel, Phys. Rev. Lett. 90, 170601, (2003)
  9. V. V. Belik and D. Brockmann, New. J. Phys. 9, 54, (2007)
  10. P. Garbaczewski and R. Olkiewicz, J. Math. Phys. 40, 1057, (1999)
  11. P. Garbaczewski and R. Olkiewicz, J. Math. Phys. 41, 6843, (2000)
  12. P. Garbaczewski and V. A. Stephanovich, Phys. Rev. E 80, 031113, (2009)
  13. P. Garbaczewski and V. Stephanovich, Open Systems & Information Dynamics, 17, 287, (2010)
  14. P. Garbaczewski and V. A. Stephanovich, Physica A 389, 4410, (2010)
  15. P. Garbaczewski and V. Stephanovich, Phys. Rev. E 84, 011142 (2011)
  16. P. Garbaczewski and V. A. Stephanovich, Acta Phys. Pol. B 43, 977, (2012)
  17. J. Lőrinczi, F. Hiroshima and V. Betz, Feyman-Kac-Type Theorems and Gibbs Measures on Path Space, (Studies in Mathematics 34, Walter de Gruyter, Berlin, 2011)
  18. R. Vilela Mendes, J. Math. Phys. 27,178, (1986)
  19. S. Eleutério and R. Vilela Mendes, Phys. Rev. B 50, 5035, (1993).
  20. I. Eliazar and J. Klafter, J. Stat. Phys. 111, 739, (2003)
  21. V. V. Janovsky et al., Physica A 282, 13, (2000)
  22. A. Dubkov and B. Spagnolo, Fluct. Noise Lett. 5, L267, (2005)
  23. H. Touchette and E. G. D. Cohen, Phys. Rev. E 76, 020101, (2007)
  24. T. Srokowski, Physica A 388, 1057, (2009)
  25. T. D. Gillespie, J. Phys. Chemistry, 81 (25): 2340, (1977)
  26. T. D. Gillespie, J. Comput. Physics, 22 (4): 403, (1976)
  27. M. Matsumoto and T. Nishimura, ACM Trans. Model. Comput. Simul. 8, 3 (1998)
  28. R. N. Mantegna and H. E. Stanley, Phys. Rev. Lett. 73, 2946, (1994)
  29. I. M. Sokolov, private communication
  30. C-code is available upon request
  31. A. Janicki and A. Weron, Simulation and chaotic behavior of -stable stochastic processes, (M. Dekker, New York, 1994)
  32. P. Hänggi, F. Marchesoni, Rev. Mod. Phys. 81, 387, (2009)
  33. S. I. Denisov, T. V. Lyutyy et al, Phys. Rev. E 79, 051102, (2009)
  34. H. Gimperlein, S. Wessel, J. Schmiedmayer, and L. Santos, Phys. Rev. Lett. 95, 170401, (2005)
  35. K. Binder, A.P. Young, Rev. Mod. Phys. 58, 801 (1986).
  36. U.T. Höchli, K. Knorr, A. Loidl, Adv. Phys. 39, 405 (1990).
  37. G. Parisi, Field Theory, Disorder and Simulations, (World Scientific Publishing, Singapore, 1992).
  38. D. Kedzierski, E.V. Kirichenko, V.A. Stephanovich, Phys. Lett. A 375, 685 (2011).
  39. A.M. Stoneham, Rev. Mod. Phys. 41, 1 (1969).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description