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[] 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 dynamics which sets down at a predefined target, in consistency with the associated kinetic (master) equation.
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 , 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 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 , 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 nonsymmetric 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 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” .
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 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 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 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 [4] disconnects the Langevin-modeled Fokker-Planck equation for any Lévy-stable noise 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 ): see, for example, [6, 7], where we assume that asymptotically sets down at a square root of a well-defined pdf . The semigroup potential follows from the compatibility condition:
In this particular context, while adopting a multiplicative decomposition of the time-dependent pdf , a novel fractional generalization of the Fokker-Planck equation governing the time evolution of 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: 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], and the Lévy measure in is given by It is the quantity which has an interpretation of the jump transition rate from the point to another point . The potential function can be chosen quite arbitrarily. However, we need to secure a normalization of the target pdf . We note that becomes a genuine stationary solution of (7) once we let and .
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 for a while) all jumps form a fixed point to any point in the set .
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 in , 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 . We mimic previously devised 1D steps [12], while adopting them to the 2D situation. Namely, (7) can be rewritten as follows: where 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 and .(iv)Using a random number generator draw from a uniform distribution.(v)Employing above and and the identities where and , find and corresponding to the “transition channel” .(vi)Draw a new number from a uniform distribution.(vii)Reset time label , 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 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, i , provided that one keeps in memory our convention not to reproduce the , index in 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 with , when .
3.1. Harmonic Confinement: Gaussian Target
Let us prescribe an asymptotic invariant pdf to be in a 2D Gaussian form: As an exemplary source of random noise we assume the Cauchy driver, that is, 2D Lévy-stable noise with the stability index . Accordingly, the jump transition rate reads
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 where Numerical routines were written in terms of C-codes [15].
In Figure 1 we have depicted the statistical data inferred from 100 000 trajectories, for three running time instants , , and , together with the asymptotic expression (14). The right-hand-side column depicts projections of those data upon the plane. A substantial increase of the analyzed ensemble data (like, e.g., , , or ) is merely a matter of the simulation time span and adds nothing inspiring to the obtained behavior.
(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 , in consistency with the ultimate target outcome, whose height directly follows from . 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 Figure 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 from the origin at the running time instant , while is a mean square distance from . In view of the dynamics should set down at , while this is of at . Figure 3 depicts the evolution of and , inferred from the simulated sample of jumping paths.
(a)
(b)
The observed convergence and 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: Like previously, we take the Cauchy driver, , as a reference Lévy stable noise. Accordingly, Proceeding like in the Gaussian case, for small and , (20) can be approximated by 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 Figure 4, where a surrogate pdf evolution is displayed. Figure 5 reproduces the 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.
(a)
(b)
(c)
(d)
3.3. Locally Periodic Confinement in
We consider target pdf whose Boltzmannian 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: The normalization constant actually reads .
Subsequently adopted numerical integration routines heavily rely on the experience gained during our previous case studies. In Figure 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 .
(a)
(b)
(c)
(d)
We are aware of the fact that the number of trajectories may be considered as a 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 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 Figure 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.
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 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 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 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.