Royal Society Publishing

Two-particle anomalous diffusion: probability density functions and self-similar stochastic processes

Gianni Pagnini , Antonio Mura , Francesco Mainardi

Abstract

Two-particle dispersion is investigated in the context of anomalous diffusion. Two different modelling approaches related to time subordination are considered and unified in the framework of self-similar stochastic processes. By assuming a single-particle fractional Brownian motion and that the two-particle correlation function decreases in time with a power law, the particle relative separation density is computed for the cases with time sub-ordination directed by a unilateral M-Wright density and by an extremal Lévy stable density. Looking for advisable mathematical properties (for instance, the stationarity of the increments), the corresponding self-similar stochastic processes are represented in terms of fractional Brownian motions with stochastic variance, whose profile is modelled by using the M-Wright density or the Lévy stable density.

1. Introduction

Diffusion phenomena occur in many natural systems and they are investigated in several disciplines. In particular, to study diffusion it is important to compute the statistics of a tracer concentration field θ(x,t). In fact, let p(x;t | y,t0) be the probability density function (PDF) of finding a particle in x at time t knowing that it was in y at time t0, and let S(y,t0) be the initial tracer source, then it is well known that for single-particle diffusion the mean value of the tracer concentration field 〈θ〉 is Embedded Image 1.1 and for two-particle diffusion the two-point covariance is Embedded Image 1.2

Studying two-particle separation is of paramount importance in applied sciences because it is related to the local fluctuations of the concentration field θ′=θ−〈θ〉, whose variance, that is, Embedded Image 1.3 is strongly relevant for determining, for example, the safe concentration threshold in industrial accidents [1] or the segregation coefficient in non-premixed reacting mixtures [2].

In the formula of the concentration fluctuation variance (1.3), 〈θ2(x,t)〉 is obtained from (1.2) by setting x1=x2=x and t1=t2=t, then it emerges to be established by the backward diffusion of two particles with very close initial positions, in the limit x1=x2=x, at the instant t1=t2=t, and the separation is estimated at t01=t02=t0<t [3]. Hereinafter, the dependence on the initial condition of the PDFs will be dropped because initial particle positions are assumed to be coincident and stated equal to 0.

From (1.2) it follows that diffusion differences generate different quantitative results, which are important in the applied problems connected to concentration fluctuations. Diffusive processes are generally classified as normal when particle displacement variance grows linearly in time; otherwise they are classified as anomalous. Moreover, normal diffusion is also associated with the Gaussian PDF for particle displacement, and, by using the stochastic process terminology, normal diffusion is also called Brownian motion (BM). A Gaussian anomalous diffusion can be obtained from normal diffusion by assuming a diffusion coefficient that is dependent on time. This can lead, for instance, to a Gaussian process such as the fractional BM (fBM) that generalizes the standard BM.

Anomalous diffusion appears in several different fields, and fractional calculus emerged to be a powerful tool to manage memory effects and long-range dependence [4]. Fractional systems have been investigated in their dynamic evolution by Lyapunov exponent analysis [5] as well as in their ability to encode memory in nuclear magnetic resonance phenomena [6,7].

Anomalous diffusion is called fast diffusion when the variance of the particle spreading grows according to a power law with exponent greater than 1, and slow diffusion when such an exponent is less than 1.

In this paper, two-particle diffusion is analysed in the context of anomalous diffusion.

A physical idea useful to model anomalous diffusion is related to time subordination of the Gaussian process. The resulting PDF is no longer Gaussian, and the particle displacement variance no longer grows linearly in time. In this paper, two approaches for modelling anomalous diffusion based on time subordination are discussed and unified into the framework of self-similar stochastic processes by showing that the resulting one-time PDFs share the same integral representation. The two time-subordination mechanisms are the parent-directing process [8] and the product of two independent random variables [9,10].

The paper is organized as follows. In §2, the two time-subordination mechanisms, namely the parent-directing process and the product of random variables, are introduced and unified in the context of self-similar stochastic processes. In §3, two-particle anomalous diffusion is studied by assuming a fractional BM for a single-particle trajectory and a correlation function decreasing in time with a power law. In §4, the resulting integral representations of the one-time PDFs are evaluated, and the representations in terms of H-function and by series are given together with their asymptotic behaviours. In §5, a class of stochastic processes is chosen to model the two-particle anomalous diffusion corresponding to the derived PDFs. In §6, conclusions are discussed.

2. Two time-subordination mechanisms and their unification for self-similar stochastic processes

(a) The parent-directing subordination process

At the microscopic level, time subordination is defined as a process X(t)=Y (T(t)) that is obtained by randomizing the time clock of a stochastic process Y (τ) using a new clock τ=T(t), where T(t) is a random process with non-negative increments [8].

The resulting process Y (T(t)) is said to be subordinated to Y (τ), which is called the parent process, and is directed by T(t), which is called the directing process. The time variable τ is often referred to as the operational time. The process t=t(τ), which is the inverse of τ=T(t), is called the leading process and it plays a fundamental role in the parametric approach, introduced by Gorenflo & Mainardi [11,12], to simulate fractional diffusion processes. Even if the parent process Y (τ) is Markovian, the resulting subordinated process X(t) is in general non-Markovian, and the non-local memory effects are due to the random time process τ=T(t) and to its evolution, which is in general non-local in time.

At the macroscopic-level, i.e. in terms of the particle PDF, the time subordination is embodied by the following integral formula: Embedded Image 2.1 where Embedded Image is the PDF (of x evolving in τ) of the parent process, and φ(τ,t) is the PDF (of τ evolving in t) of the directing process.

(b) Independent random variables product and the time-subordination integral

It is well known that the PDF of the product of two independent random variables is given by the Mellin convolution1 of the two corresponding PDFs [9,10]. Let Z1 and Z2 be two real independent random variables whose PDFs are p1(z1) and p2(z2), respectively, with z1R and z2R+. Then, the joint PDF is p(z1,z2)=p1(z1)p2(z2). Let Z be the random variable obtained by the product of Z1 and Embedded Image, that is, Embedded Image 2.2 so that Embedded Image, then, carrying out the variable transformations z1=z/ξγ and z2=ξ, it follows that p(z,ξ)dzdξ=p1(z/ξγ)p2(ξ)Jdzdξ, where J=1/ξγ is the Jacobian of the transformation. Integrating in dξ, the PDF of Z emerges to be Embedded Image 2.3

By applying the changes of variable z=xtγω and ξ=τtω, the typical time-subordination integral (2.1) is recovered from (2.3) by setting Embedded Image 2.4

(c) Unification of time-subordination mechanisms for self-similar stochastic processes

A process W(t), t≥0, is a self-similar process with self-similarity exponent H (Hurst exponent) if, for all a>0, the processes W(at) and aHW(t) have the same finite-dimensional distributions.

If the parameter a is turned into a random variable, by setting Y (τ)=W(τ) and T(t)=at, according to definitions given in §2a, so that X(t)=Y (T(t))=W(at), and by setting Embedded Image, according to definitions given in §2b, a time subordination is recovered where the parent process is a self-similar process and the operational time is simply a line with stochastic slope.

Then, in the case described above, both mechanisms, i.e. explicit time subordination with a stochastic operational time and product of random variables, can be seen as equivalent in the context of subordinated processes.

(d) The M-Wright and Lévy directing processes

Two noteworthy examples of PDFs for anomalous diffusion are obtained, despite the effective microscopic stochastic formulation, by subordination-type integral (2.1) or (2.3) of the Markovian BM, with PDF Embedded Image 2.5 and variance 〈x2〉=2t, by selecting as the PDF of the directing process either a unilateral M-Wright function, as introduced in the so-called generalized grey BM (ggBM) [1315], also called Erdélyi–Kober fractional diffusion [16, 17], or an extremal Lévy stable density, as considered in Mainardi et al. [1820]. In fact, in the former case, for 0<β≤1 and 0<α<2, it holds [16] that Embedded Image 2.6 which includes as special cases the fBM for β=1, the grey BM, in the sense of Schneider [21,22], for 0<α=β<1 and the BM for α=β=1 [23]. In the last case, for 0<α≤2, we have [9,10] Embedded Image 2.7 The function Mν (with 0<ν≤1) is the so-called M-Wright function, also called the Mainardi function, and Embedded Image is the Lévy stable density of order 0<α≤2 and asymmetry parameter Embedded Image, which is symmetric for θ=0 and one-sided on the positive semi-axis for θ=−α<1. These functions can be defined, respectively, by their Mellin transform pair (N1) as (see eqn. (6.1) in Mainardi et al. [9]) Embedded Image 2.8 and (see eqn. (5.1) in Mainardi et al. [9]) Embedded Image 2.9 Remember here that the symmetric M-Wright function and the Lévy stable density provide the Green functions of the time-fractional and the space-fractional diffusion equations, respectively [18].

3. Two-particle anomalous diffusion

The problem of two-particle diffusion can be restated taking into account the single-particle motion and the relative separation between the particles. In terms of the PDF, this means that p(x1,x2;t) dx1 dx2=p(x,δr;t) dx dδr, after introducing the variable x1=x and x2=xδr.

However, following the turbulent dispersion literature [24,25], the problem of two-particle diffusion can be analysed in terms of the particle separation δr=x1x2 and of the centre-of-mass position xCM=(x1+x2)/2 that, for mathematical convenience (i.e. the Jacobian of the variable transformation is equal to 1), are transformed into: Embedded Image and Embedded Image.

Because the two particles are identical, it holds that Embedded Image and in general their motion is correlated 〈x1x2〉=σ2ρ, then the joint Gaussian PDF is Embedded Image 3.1 so that Embedded Image and the following three normalization conditions hold: Embedded Image, Embedded Image and Embedded Image. From (3.1), it follows that 〈Δ2〉=σ2(1−ρ) and 〈Σ2〉=σ2(1+ρ).

Assuming that, before the Einstein regime, i.e. σ2=2t for t≫1, the trajectories x1 and x2 follow a fBM, then the same holds for x, and assuming also that the correlation function decreases in time according to a power law then Embedded Image 3.2

For the particle separation type variable Δ, it holds that Embedded Image 3.3 It is worth noting that the regime established for t≪1 in (3.3) is a generalization of that in the most common Ornstein–Uhlenbeck process, i.e. the combination of the ballistic regime σ2t2 with the exponential correlation function ρ=et≃1−t, so that σ2(1−ρ)≃t3.

On the other hand, for the centre-of-mass type variable Σ, it emerges that in the considered context its variance is Embedded Image, when t≪1, and 〈Σ2〉=2t, when t≫1, so that the simple power law scaling when t≪1 is lost.

Hence, in the range t≪1, it holds that, according to the notation introduced at the beginning of the section, i.e. x=x1, Σ and Embedded Image, Embedded Image 3.4 Embedded Image 3.5 Embedded Image 3.6

However, because the centre-of-mass motion is a single-point trajectory then, in agreement with the assumed single-particle trajectory for x1 and x2, it is governed by the fBM and its time-dependent diffusion coefficient is Embedded Image, but the two-point trajectory Δ can have long-range correlation and show anomalous diffusion. Then, hereinafter, the analysis is performed solely for Embedded Image. Such non-local memory effects can be modelled by the subordination integral (2.1) and, in this framework, the PDF of particle relative separation emerges to be determined by Embedded Image 3.7

Finally, adopting the same directing processes of formulae (2.6) and (2.7), the resulting PDFs for the two-particle anomalous diffusion emerge to be, with 0<β≤1 and 0<α<2, Embedded Image 3.8 and Embedded Image 3.9

Clearly, if anomalous diffusion is assumed also for the single-particle trajectory x, the analysis follows from the present one by setting μ=q. Moreover, it is important to highlight here also that, if the single-particle trajectory is assumed to be Gaussian, i.e. q=1, then in the anomalous diffusion framework the PDF of x is given by (2.6) or (2.7) depending on the selected directing process.

4. Evaluation, representations and asymptotic behaviour of the PM and of the PL probability density functions

(a) The PM probability density function

To determine the PM PDF defined in (3.8), by applying the Mellin transform (N1) to the first line of (3.8) and introducing the variable Δ2=4τμξ, the right-hand side becomes Embedded Image 4.1 Then, by introducing the variable τ=tαz and remembering the gamma function definition, i.e. Embedded Image, expression (4.1) turns into Embedded Image 4.2 To conclude, by the gamma function property Embedded Image and the Mellin transform pair for the M-Wright function (2.8), after the Mellin inverse transformation (N1), the Mellin–Barnes integral representation of the PM PDF is Embedded Image 4.3 where Embedded Image is the contour path encircling the poles of Γ(s) and those of Γ(1−μ/2+μs/2).

In terms of the H-function (see appendix A), when |Δ/tαμ/2|→0, the density PM is Embedded Image 4.4 and, by applying the residue theorem, the series representation of (4.3) is Embedded Image 4.5

The asymptotic behaviour for Embedded Image can be computed by formula (A3) and it emerges to be of the exponential type as follows: Embedded Image 4.6 where Embedded Image

(b) The PL probability density function

The determination of the PL PDF defined in (3.9) is analogous to the determination of the PM PDF. Actually, by applying the Mellin transform (N1) to the first line of (3.9) and after introducing the variable Δ2=4τμξ, the right-hand side becomes Embedded Image 4.7 Then, by introducing the variable τ=t2/αz and using the gamma function definition, expression (4.7) becomes Embedded Image 4.8 To conclude, by the gamma function property Embedded Image and the Mellin transform pair for the Lévy stable density (2.9), after the Mellin inverse transformation (N1), the Mellin–Barnes integral representation of the PL PDF is Embedded Image 4.9 where Embedded Image is the contour path encircling the poles of Γ(s).

In terms of the H-function (see appendix A), it holds that Embedded Image 4.10 and Embedded Image 4.11

Then, by applying the residue theorem, the corresponding series representations of (4.9) are Embedded Image 4.12 when |Δ/tμ/α|→0, and Embedded Image 4.13 when Embedded Image.

5. H-sssi processes for two-particle anomalous diffusion

In order to choose a class of stochastic processes to model two-particle anomalous diffusion with PDFs according to (3.8) and (3.9), the same approach introduced by Mura [13] (see also [14,15]) to characterize the so-called ggBM is adopted.

The ggBM is a class of H-sssi processes, where sssi means self-similar with stationary increments, that generalize Gaussian processes (which are recovered as a special case) and defined by only the autocovariance structure. This property can be easily deduced by noting that the ggBM can be represented by the process Embedded Image, t≥0, where Xα(t) is a Gaussian stochastic process and Λβ is a suitable chosen independent non-negative random variable.

In fact, let Bα,β(t), t≥0, be a ggBM, then Embedded Image 5.1 where the symbol =d denotes the equality of the finite-dimensional distribution, the stochastic process Xα(t) is a standard fBM with Hurst exponent α/2 and Λβ is an independent non-negative random variable with PDF Mβ(τ), τ≥0.

Representation (5.1) can be backwardly obtained from the subordination definition given in §2b. In fact, by comparing (2.3) and (2.4) with (2.6) it follows that Z=Xtα/2 and, because Embedded Image, it holds that Embedded Image, where γ=1/2, the variable Z1 is Gaussian and Z2 is distributed according to Mβ. Finally, representation (5.1) is recovered with X=Bα,β(t), Xα=Z1tα/2 and Λβ=Z2.

It is worth noting here that representation (5.1) permits a number of questions to be solved, in particular those related to the distribution properties of Bα,β(t), because they can be reduced to questions concerning the fBM Xα(t) that, because Xα(t) is a Gaussian process, has been largely studied in the literature.

For instance, the Hölder continuity of the Bα,β(t) trajectories follows immediately from those of Xα(t), that is, Embedded Image 5.2 Moreover, representation (5.1) highlights the stationary increment property of the ggBm, and it emerges to be suitable for path simulations.

With the same backward-derivation method described to obtain representation (5.1) from (2.3) plus (2.4) and (2.6), it follows that the symmetric Lévy process governed by (2.7) is represented by Embedded Image 5.3 where the stochastic process X2/α(t) is a standard fBM with Hurst exponent 1/α and ℓα/2 is an independent non-negative random variable with PDF Embedded Imageτ≥0.

Finally, for the two-particle anomalous diffusion, the self-similar stochastic process representations are straightforwardly obtained. In fact, comparing (2.3) plus (2.4) and (3.8), it follows that Embedded Image 5.4 where the stochastic process Xμα(t) is a standard fBM with Hurst exponent μα/2 and Λβ is an independent non-negative random variable with PDF Mβ(τ), τ≥0.

Comparing (2.3) plus (2.4) and (3.9), it follows that Embedded Image 5.6 where the stochastic process X2μ/α(t) is a standard fBM with Hurst exponent μ/α and ℓα/2 is an independent non-negative random variable with PDF Embedded Image, τ≥0.

6. Conclusion

In the present paper, two-particle dispersion is investigated in the context of anomalous diffusion. Two-particle diffusion is important for the statistical analysis of fluctuations of a tracer concentration field. Because anomalous diffusion can be modelled by time subordination of Gaussian stochastic processes, here it is shown that, in the framework of self-similar stochastic processes, two different subordination modelling approaches, namely the parent-directing process and the product of random variables, can be unified because they give the same integral formula for the resulting PDF.

By assuming a single-particle fBM and that the particle correlation function decreases in time with a power law, it follows that the centre-of-mass trajectory is an fBM, too, whose variance no longer has a simple power law scaling because of the two-particle correlation. Furthermore, the two-particle separation PDF is computed for the cases with time subordination directed by the M-Wright density and by the Lévy stable density. The resulting PDFs are evaluated and represented by the Mellin–Barnes integral, in terms of the H-function and by series. Moreover, the asymptotic behaviours are also shown. In particular, in the case with operational time driven by the M-Wright density, the PDF asymptotically decreases with a stretched exponential law, whereas in the case driven by the Lévy stable density, the PDF asymptotically decreases with a power law.

In order to choose a class of stochastic processes with PDFs according to those previously established, the same method introduced for the derivation of the ggBM is adopted. In analogy with the ggBM, this class of stochastic processes emerges to be determined by the product of a Gaussian process with an opportune Hurst exponent, i.e. an opportune fBM, and an independent random variable distributed according to the M-Wright density or the Lévy stable density. Then the processes belonging to such a class are referred to as H-sssi and are defined by only the autocovariance structure.

The determination of the master equation of this class of H-sssi processes that is satisfied by the evaluated PDFs will be the basis of future developments of the present research.

Acknowledgements

G.P. was funded by Regione Autonoma della Sardegna (PO Sardegna FSE 2007–2013 sulla L.R. 7/2007 ‘Promozione della ricerca scientifica e dell’innovazione tecnologica in Sardegna’).

Appendix A

The H-function is defined by means of a Mellin–Barnes type integral as follows [26]: Embedded Image A 1 with Embedded Image A 2 where an empty product is always interpreted as unity, {m,n,p,q}∈N0 with 1≤mq and 0≤np, {Ai,Bj}∈R+ and {ai,bj}∈R, or C, with i=1,…,p and j=1,…,q such that Ai(bj+k)≠Bj(ai−ℓ−1) with k and ℓ∈N0, i=1,…,n and j=1,…,m. The poles of the integrand in (A1) are assumed to be simple. The integration path Embedded Image encircles all the poles of Γ(bj+Bjs) with j=1,…,m.

The H-function is an analytic function of z and exists for all z≠0 when q≥1 and μ>0 or for 0<|z|<Δ when q≥1 and μ=0 where Embedded Image For other existence conditions see Mathai et al. [26].

The asymptotic expansion for Embedded Image is obtained by integration around the poles of Γ(1−aiAis) with Embedded Image. Actually this is similar to the exchange s→−s and then passing from the series of powers of z to a series of powers of 1/z, from which the asymptotic expansion follows.

In the particular case with n=0, the asymptotic behaviour for Embedded Image is of exponential type and determined for real z by the formula Embedded Image A 3 where Embedded Image

Footnotes

  • One contribution of 14 to a Theme Issue ‘Fractional calculus and its applications’.

  • 1 The Mellin transform pair Embedded Image, with rR+ and sC, is defined by Embedded Image The transformed function ψ*(s) exists if the integral Embedded Image is bounded and this constraint is met in the vertical strip a<σ=Re(s)<b, where the boundary values a and b follow from the analytic structure of ψ(r) provided that |ψ(r)|≤Mra when r→0+ and |ψ(r)|≤Mrb when Embedded Image. The Mellin convolution integral corresponds to the pair Embedded Image

References

View Abstract