## 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*,*t*_{0}) be the probability density function (PDF) of finding a particle in *x* at time *t* knowing that it was in *y* at time *t*_{0}, and let *S*(*y*,*t*_{0}) be the initial tracer source, then it is well known that for single-particle diffusion the mean value of the tracer concentration field 〈*θ*〉 is
1.1
and for two-particle diffusion the two-point covariance is
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,
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 *x*_{1}=*x*_{2}=*x* and *t*_{1}=*t*_{2}=*t*, then it emerges to be established by the *backward diffusion* of two particles with very close initial positions, in the limit *x*_{1}=*x*_{2}=*x*, at the instant *t*_{1}=*t*_{2}=*t*, and the separation is estimated at *t*_{01}=*t*_{02}=*t*_{0}<*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:
2.1
where 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 convolution^{1} of the two corresponding PDFs [9,10]. Let *Z*_{1} and *Z*_{2} be two real independent random variables whose PDFs are *p*_{1}(*z*_{1}) and *p*_{2}(*z*_{2}), respectively, with *z*_{1}∈*R* and *z*_{2}∈*R*^{+}. Then, the joint PDF is *p*(*z*_{1},*z*_{2})=*p*_{1}(*z*_{1})*p*_{2}(*z*_{2}). Let *Z* be the random variable obtained by the product of *Z*_{1} and , that is,
2.2
so that , then, carrying out the variable transformations *z*_{1}=*z*/*ξ*^{γ} and *z*_{2}=*ξ*, it follows that *p*(*z*,*ξ*)d*zdξ*=*p*_{1}(*z*/*ξ*^{γ})*p*_{2}(*ξ*)*J*d*z*d*ξ*, where *J*=1/*ξ*^{γ} is the Jacobian of the transformation. Integrating in d*ξ*, the PDF of *Z* emerges to be
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
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 *a*^{H}*W*(*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 §2*a*, so that *X*(*t*)=*Y* (*T*(*t*))=*W*(*at*), and by setting , according to definitions given in §2*b*, 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
2.5
and variance 〈*x*^{2}〉=2*t*, 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) [13–15], also called Erdélyi–Kober fractional diffusion [16, 17], or an extremal Lévy stable density, as considered in Mainardi *et al*. [18–20]. In fact, in the former case, for 0<*β*≤1 and 0<*α*<2, it holds [16] that
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]
2.7
The function *M*_{ν} (with 0<*ν*≤1) is the so-called *M*-Wright function, also called the Mainardi function, and is the Lévy stable density of order 0<*α*≤2 and asymmetry parameter , 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])
2.8
and (see eqn. (5.1) in Mainardi *et al.* [9])
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*(*x*_{1},*x*_{2};*t*) d*x*_{1} d*x*_{2}=*p*(*x*,*δr*;*t*) d*x* d*δr*, after introducing the variable *x*_{1}=*x* and *x*_{2}=*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*=*x*_{1}−*x*_{2} and of the centre-of-mass position *x*_{CM}=(*x*_{1}+*x*_{2})/2 that, for mathematical convenience (i.e. the Jacobian of the variable transformation is equal to 1), are transformed into: and .

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

Assuming that, before the Einstein regime, i.e. *σ*^{2}=2*t* for *t*≫1, the trajectories *x*_{1} and *x*_{2} 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
3.2

For the particle separation type variable *Δ*, it holds that
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 *σ*^{2}≃*t*^{2} with the exponential correlation function *ρ*=e^{−t}≃1−*t*, so that *σ*^{2}(1−*ρ*)≃*t*^{3}.

On the other hand, for the centre-of-mass type variable *Σ*, it emerges that in the considered context its variance is , when *t*≪1, and 〈*Σ*^{2}〉=2*t*, 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*=*x*_{1}, *Σ* and ,
3.4
3.5
3.6

However, because the centre-of-mass motion is a single-point trajectory then, in agreement with the assumed single-particle trajectory for *x*_{1} and *x*_{2}, it is governed by the fBM and its time-dependent diffusion coefficient is , but the two-point trajectory *Δ* can have long-range correlation and show anomalous diffusion. Then, hereinafter, the analysis is performed solely for . 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
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,
3.8
and
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 *P*_{M} and of the *P*_{L} probability density functions

### (a) The *P*_{M} probability density function

To determine the *P*_{M} 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
4.1
Then, by introducing the variable *τ*=*t*^{α}*z* and remembering the gamma function definition, i.e. , expression (4.1) turns into
4.2
To conclude, by the gamma function property 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 *P*_{M} PDF is
4.3
where 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 *P*_{M} is
4.4
and, by applying the residue theorem, the series representation of (4.3) is
4.5

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

### (b) The *P*_{L} probability density function

The determination of the *P*_{L} PDF defined in (3.9) is analogous to the determination of the *P*_{M} 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
4.7
Then, by introducing the variable *τ*=*t*^{2/α}*z* and using the gamma function definition, expression (4.7) becomes
4.8
To conclude, by the gamma function property 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 *P*_{L} PDF is
4.9
where is the contour path encircling the poles of *Γ*(*s*).

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

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

## 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 , *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
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 §2*b*. In fact, by comparing (2.3) and (2.4) with (2.6) it follows that *Z*=*Xt*^{−α/2} and, because , it holds that , where *γ*=1/2, the variable *Z*_{1} is Gaussian and *Z*_{2} is distributed according to *M*_{β}. Finally, representation (5.1) is recovered with *X*=*B*_{α,β}(*t*), *X*_{α}=*Z*_{1}*t*^{α/2} and *Λ*_{β}=*Z*_{2}.

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,
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
5.3
where the stochastic process *X*_{2/α}(*t*) is a standard fBM with Hurst exponent 1/*α* and ℓ_{α/2} is an independent non-negative random variable with PDF , *τ*≥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
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
5.6
where the stochastic process *X*_{2μ/α}(*t*) is a standard fBM with Hurst exponent *μ*/*α* and ℓ_{α/2} is an independent non-negative random variable with PDF , *τ*≥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]:
A 1
with
A 2
where an empty product is always interpreted as unity, {*m*,*n*,*p*,*q*}∈*N*_{0} with 1≤*m*≤*q* and 0≤*n*≤*p*, {*A*_{i},*B*_{j}}∈*R*^{+} and {*a*_{i},*b*_{j}}∈*R*, or *C*, with *i*=1,…,*p* and *j*=1,…,*q* such that *A*_{i}(*b*_{j}+*k*)≠*B*_{j}(*a*_{i}−ℓ−1) with *k* and ℓ∈*N*_{0}, *i*=1,…,*n* and *j*=1,…,*m*. The poles of the integrand in (A1) are assumed to be simple. The integration path encircles all the poles of *Γ*(*b*_{j}+*B*_{j}*s*) 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
For other existence conditions see Mathai *et al*. [26].

The asymptotic expansion for is obtained by integration around the poles of *Γ*(1−*a*_{i}−*A*_{i}*s*) with . 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 is of exponential type and determined for real *z* by the formula
A 3
where

## Footnotes

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

↵1 The Mellin transform pair , with

*r*∈*R*^{+}and*s*∈*C*, is defined by The transformed function*ψ**(*s*) exists if the integral 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*)|≤*Mr*^{−a}when*r*→0^{+}and |*ψ*(*r*)|≤*Mr*^{−b}when . The Mellin convolution integral corresponds to the pair

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.