Heterogeneous diffusion in comb and fractal grid structures

TL;DR

This paper provides exact analytical solutions for diffusion in comb and fractal grid structures with position-dependent diffusion coefficients, revealing various anomalous diffusion regimes depending on structural parameters.

Contribution

It introduces exact solutions for diffusion behavior in comb and fractal grid structures with power-law dependent diffusion coefficients, highlighting different anomalous regimes.

Findings

Mean square displacement varies as t^{1/(2-α)} for comb structures.

Fractal grid structures show mean square displacement scaling as t^{(1+ν)/(2-α)}.

Probability distributions are expressed using Fox H-functions.

Abstract

We give an exact analytical results for diffusion with a power-law position dependent diffusion coefficient along the main channel (backbone) on a comb and grid comb structures. For the mean square displacement along the backbone of the comb we obtain behavior $\langle x^2(t)\rangle\sim t^{1/(2-\alpha)}$, where $\alpha$ is the power-law exponent of the position dependent diffusion coefficient $D(x)\sim |x|^{\alpha}$. Depending on the value of $\alpha$ we observe different regimes, from anomalous subdiffusion, superdiffusion, and hyperdiffusion. For the case of the fractal grid we observe the mean square displacement, which depends on the fractal dimension of the structure of the backbones, i.e., $\langle x^2(t)\rangle\sim t^{(1+\nu)/(2-\alpha)}$, where $0<\nu<1$ is the fractal dimension of the backbones structure. The reduced probability distribution functions for both cases are obtained by help of the Fox $H$-functions.

Equations

$$\frac{d}{dt}x(t)=\sqrt{2\mathcal{D}(x)}\zeta (t),$$

$$\frac{\partial }{\partial t}P(x,t)=\frac{\partial }{\partial x}\left[\sqrt{\mathcal{D}(x)}\frac{\partial }{\partial x}\left(\sqrt{\mathcal{D}(x)}P(x,t)\right)\right].$$

$$\mathcal{D}(x)={\mathcal{D}}_x|x{|}^{\alpha },\alpha <2.$$

$$P(x,t)=\frac{|x{|}^{-\alpha /2}}{\sqrt{4\pi {\mathcal{D}}_{x}t}}\exp \left(-\frac{|x{|}^{2-\alpha }}{(2-\alpha {)}^2{\mathcal{D}}_{x}t}\right),$$

$$⟨x^2(t)⟩={\int }_{-\infty }^{\infty }dx{x}^{2}P(x,t)\simeq \frac{t^{\frac{2}{2-\alpha }}}{\Gamma \left(1+\frac{2}{2-\alpha }\right)}.$$

$$\frac{\partial }{\partial t}P(x,y,t)={\mathcal{D}}_x\delta (y)\frac{{\partial }^2}{\partial x^2}P(x,y,t)+{\mathcal{D}}_y\frac{{\partial }^2}{\partial y^2}P(x,y,t),$$

$$\frac{\partial }{\partial t}P(x,y,t)$$

$$P(x,y,t=0)=\delta (x)\delta (y),$$

$$\frac{\partial }{\partial t}P(x,y,t)={\mathcal{D}}_x\delta (y)\frac{\partial }{\partial x}\left[|x{|}^{\alpha /2}\frac{\partial }{\partial x}\left(|x{|}^{\alpha /2}P(x,y,t)\right)\right]+{\mathcal{D}}_y\frac{{\partial }^2}{\partial y^2}P(x,y,t).$$

$$sP(x,y,s)-P(x,y,t=0)$$

$$P(x,y,s)=g(x,s)\exp \left(-\sqrt{\frac{s}{{\mathcal{D}}_y}}|y|\right),$$

$$P(x,y=0,s)=g(x,s).$$

$$p_1(x,t)={\int }_{-\infty }^{\infty }dyP(x,y,t),$$

$$p_1(x,s)=2g(x,s)\sqrt{\frac{{\mathcal{D}}_y}{s}}.$$

$$s{p}_1(x,s)-p_1(x,t=0)={\mathcal{D}}_x\frac{\partial }{\partial x}\left[|x{|}^{\alpha /2}\frac{\partial }{\partial x}\left(|x{|}^{\alpha /2}g(x,s)\right)\right],$$

$$2\sqrt{{\mathcal{D}}_y}{s}^{1/2}g(x,s)-{\mathcal{D}}_x\frac{\partial }{\partial x}\left[|x{|}^{\alpha /2}\frac{\partial }{\partial x}\left(|x{|}^{\alpha /2}g(x,s)\right)\right]=\delta (x).$$

$$2\sqrt{{\mathcal{D}}_y}{s}^{1/2}|x{|}^{-\alpha /2}f(x,s)-{\mathcal{D}}_x\frac{\partial }{\partial x}\left[|x{|}^{\alpha /2}\frac{\partial }{\partial x}f(x,s)\right]=\delta (x).$$

$$2\sqrt{{\mathcal{D}}_y}{s}^{1/2}{z}^{-\alpha /2}f(z,s)-{\mathcal{D}}_x(\alpha /2)z^{\alpha /2-1}\frac{\partial }{\partial z}f(z,s)-{\mathcal{D}}_x{z}^{\alpha /2}\frac{{\partial }^2}{\partial z^2}f(z,s)$$

$$-2{\mathcal{D}}_x{z}^{\alpha /2}\frac{\partial }{\partial z}f(z,s)\delta (x)=\delta (x).$$

$$2\sqrt{{\mathcal{D}}_y}{s}^{1/2}{z}^{-\alpha /2}f(z,s)-{\mathcal{D}}_x(\alpha /2)z^{\alpha /2-1}\frac{\partial }{\partial z}f(z,s)-{\mathcal{D}}_x{z}^{\alpha /2}\frac{{\partial }^2}{\partial z^2}f(z,s)=0.$$

$$-2{\mathcal{D}}_x{z}^{\alpha /2}{\frac{\partial }{\partial z}f(z,s)|}_{z=0}=1.$$

$$f(x,s)$$

$$\mathcal{C}(s)=\frac{1}{\sqrt{4\pi }{\mathcal{D}}_x}\frac{2}{(2-\alpha {)}^{1/2}}{\left(\frac{{\mathcal{D}}_x}{2\sqrt{{\mathcal{D}}_y}}\right)}^{1/4}{s}^{-1/8}.$$

$$\displaystyle p_{1}(x,t)=\frac{|x|^{-\alpha/2}}{\sqrt{4\pi\frac{\mathcal{D}_{x}}{2\sqrt{\mathcal{D}_{y}}}t^{1/2}}}H_{1,2}^{2,0}\left[\left.\frac{1}{(2-\alpha)^{2}}\frac{2\sqrt{\mathcal{D}_{y}}}{\mathcal{D}_{x}}\frac{|x|^{2-\alpha}}{t^{1/2}}\right|\left.\begin{array}[]{l}(3/4,1/2)\\ (1/2,1),(0,1)\end{array}\right.\right].$$

$$⟨x^2(t)⟩\simeq \frac{t^{\frac{1}{2-\alpha }}}{\Gamma \left(1+\frac{1}{2-\alpha }\right)}\sim t^{\beta }.$$

$$\Delta x=\delta t\cdot \alpha /4\cdot |x{|}^{\alpha -1}+\sqrt{\delta t}\cdot \xi \cdot |x{|}^{\alpha /2},$$

$$\frac{\partial }{\partial t}P(x,y,t)$$

$$s{p}_1(x,s)-\delta (x)={\mathcal{D}}_x\frac{\partial }{\partial x}\left[|x{|}^{\alpha /2}\frac{\partial }{\partial x}\left(|x{|}^{\alpha /2}\sum _{l_j\in {\mathcal{S}}_{\nu }}P(x,y=l_j,s)\right)\right].$$

$$P(x,y,s)=g(x,s)\exp \left(-\sqrt{\frac{s}{{\mathcal{D}}_y}}|y|\right),$$

$$\sum _{l_j\in {\mathcal{S}}_{\nu }}P(x,y=l_j,s)=g(x,s)\frac{1}{\Gamma (\nu )}{\int }_0^{\infty }dl{l}^{\nu -1}\exp \left(-\sqrt{\frac{s}{{\mathcal{D}}_y}}l\right)$$

$$=g(x,s){\left(\frac{{\mathcal{D}}_y}{s}\right)}^{\nu /2}=\frac{s^{(1-\nu )/2}}{2{\mathcal{D}}_y^{(1-\nu )/2}}{p}_1(x,s).$$

Full text

††journal: Chaos Solitons & Fractals

1 Introduction

In many physical systems, such as transport in inhomogeneous media and plasmas [1], and diffusion on random fractals [2], the diffusion coefficient is not a constant but depends on the particle position, like in turbulent diffusion [3, 4], including turbulent two-particle diffusion [5]. The heterogeneous diffusion equation has been investigated within the continuous time random walk (CTRW) theory in [6, 7, 8], and the mean first passage time of such systems was analyzed in [9]. Lévy processes in inhomogeneous media [10] and ergodicity breaking in heterogeneous diffusion processes [11] have been investigated, as well as the influence of external potentials on heterogeneous diffusion processes was recently considered in [12]. Time and power-law position dependent diffusion coefficient were also considered in the literature [13] in analysis of $N$-dimensional diffusion equation.

The displacement $x(t)$ of a particle in a heterogeneous medium with space dependent diffusivity $\mathcal{D}(x)$ is described by the Langevin equation

[TABLE]

where $\zeta(t)$ is a white Gaussian noise with $\left\langle\zeta(t)\zeta(t^{\prime})\right\rangle=\delta(t-t^{\prime})$ and zero mean $\left\langle\zeta(t)\right\rangle=0$. In the Stratonovich interpretation this Langevin equation corresponds to the diffusion equation for the probability distribution function (PDF) [11]

[TABLE]

It is supplemented with the initial condition $P(x,t=0)=\delta(x)$, and the boundary conditions are set to zero at infinities. The diffusion coefficient has the power-law position dependent form

[TABLE]

The solution of Eq. (2) is obtained in the stretched exponential form [11]

[TABLE]

and the mean square displacement (MSD) has the power-law dependence on time

[TABLE]

This expression describes different diffusiove regimes, where for $\alpha<0$ one observes subdiffusion, normal diffusion for $\alpha=0$, superdiffusion for $0<\alpha<1$, ballistic motion for $\alpha=1$ and hyperdiffusion for $1<\alpha<2$. The case with $\alpha=2$ leads to exponentially fast spreading [6, 14]. The case with $\alpha>2$ yields localization111In Ref. [7], where inhomogeneous advection in a comb was considered, this regime has been named by negative superdiffusion. with the decay MSD $~{}t^{-\frac{2}{\alpha-2}}$.

A random walk in a simple comb structure consisting of main diffusion channel (backbone) and trapping fingers leads to anomalous diffusion with a transport exponent equal to $1/2$ [15]. It can be described by the two-dimensional diffusion equation [16]

[TABLE]

where $P(x,y,t)$ is the probability distribution function (PDF), $\mathcal{D}_{x}\delta(y)$ is the diffusion coefficient in $x$ direction with dimension $[\mathcal{D}_{x}]=\mathrm{m}^{3}/\mathrm{s}$, and $\mathcal{D}_{y}$ is the diffusion coefficient in $y$ direction with dimension $[\mathcal{D}_{y}]=\mathrm{m}^{2}/\mathrm{s}$. The $\delta$-function in Eq. (6) means that the diffusion along the $x$ direction occurs only at $y=0$ (the backbone) and the fingers play the role of traps. The comb model (6) is used to describe diffusion in low-dimensional percolation clusters [16, 17]. Comb models can further be generalized to grid and fractal grid structures [18] in which the diffusion along the $x$ direction may appear in many backbones, even infinite number of backbones which positions belong to a fractal set $\mathcal{S}_{\nu}$ with fractal dimension $0<\nu<1$. In this case anomalous diffusion is observed and the transport exponent depends on the fractal dimension $\nu$. In this paper we consider heterogeneous diffusion on such comb and fractal grid structures, where the diffusivity is position dependent with power-law diffusion coefficient of Eq. (3).

The investigation of anomalous diffusion processes in complex systems leads to appearance of fractional differintegration in the corresponding stochastic and kinetic equation representing the memory effect in the system. Therefore, the mathematical background of the theory of fractional differential and integral equations [19, 20, 21], and associated Mittag-Leffler and Fox $H$-functions [22, 23] for analysis of such processes are of the primary importance. From the other side, diffusion on fractal structures, and the connection between the fractal dimension and fractional differintegration, as well as description of fractal processes by fractional calculus have been discussed in the scientific community [24].

The paper is organized as follows. In Sec. 2 we consider a two-dimensional diffusion equation for a comb with the position dependent (power-law) diffusion coefficient along the backbone. Exact results for the PDF and MSD are obtained and various diffusion regimes are observed, such as anomalous subdiffusion, superdiffusion and hyperdiffusion. The case of heterogeneous diffusion on a fractal grid structure is considered in Sec. 3, and exact results for the PDF and MSD are derived. The summary is given in Sec. 4.

2 Heterogeneous diffusion on a comb

We consider the following heterogeneous two dimensional diffusion equation on a comb for the PDF $P(x,y,t)$

[TABLE]

where $\mathcal{D}(x)$ is the position dependent diffusion coefficient along the backbone, $\mathcal{D}_{y}$ is the diffusion coefficient along the fingers. This equation is a generalization of the one-dimensional heterogeneous diffusion equation (2) to a two-dimensional comb structure. The initial condition is

[TABLE]

and the boundary conditions for $P(x,y,t)$ and $\frac{\partial}{\partial q}P(x,y,t)$, $q=\{x,y\}$ are set to zero at infinities, $x=\pm\infty$, $y=\pm\infty$. The position dependent diffusion coefficient has power-law form (3) with $\alpha<2$, therefore the physical dimension of the diffusion coefficient along the backbone $\mathcal{D}_{x}\delta(y)$ is $\left[\mathcal{D}_{x}\delta(y)\right]=\mathrm{m}^{2-\alpha}\mathrm{s}^{-1}$, and the physical dimension of $\mathcal{D}_{y}$ is $\left[\mathcal{D}_{y}\right]=\mathrm{m}^{2}\mathrm{s}^{-1}$.

Inserting the diffusion coefficient (3) in Eq. (7) one obtains

[TABLE]

From the Laplace transform222The Laplace transform of a given function $f(t)$ is defined by $f(s)=\mathcal{L}[f(t)]=\int_{0}^{\infty}dt\,e^{-st}f(t)$., it follows

[TABLE]

We present the solution of the Eq. (10) in the form of the ansatz

[TABLE]

from where it follows that

[TABLE]

We also introduce the reduced PDF, which describes the transport along the backbones only

[TABLE]

and yields

[TABLE]

Integrating Eq. (9) over $y$, one finds

[TABLE]

where the initial condition $p_{1}(x,t=0)=\delta(x)$. Therefore, from Eqs. (14) and (13) we obtain the differential equation

[TABLE]

After the substitution $f(x,s)=|x|^{\alpha/2}g(x,s)$, from Eq. (15) we obtain

[TABLE]

We take into account symmetrical property of the equation, which is invariant with respect to inversion $x\rightarrow-x$. Therefore, in order to solve this equation, we use $z=|x|$, from where by partial differentiation with respect to $x$ we find333We also use here the following property $x=|x|{\rm sign}(x)$, and ${\rm sign}(x)\partial_{z}=\partial_{x}$.

[TABLE]

This equation splits into the system of equations

[TABLE]

and

[TABLE]

Eq. (18) is the Lommel equation, with the solution

[TABLE]

where $\mathcal{C}(s)$ is a function which depends on $s$, $K_{\nu}(z)$ is the modified Bessel function (of the third kind) and $H_{p,q}^{m,n}(z)$ is the Fox $H$-function. The function $\mathcal{C}(s)$ is obtained by Eqs. (19) and (22), by using series representation of the modified Bessel function when $z\rightarrow 0$ (80). Therefore, we find

[TABLE]

Thereafter, by using the inverse Laplace transform formula (62) for the Fox $H$-function, and the properties (68) and (73), for the reduced PDF, we finally obtain the solution as follows

[TABLE]

Note that for $\alpha=0$ we recover the result for the classical comb (6). Graphical representation of solution (26) is plotted in Fig. 1.

From the reduced PDF we calculate the MSD

[TABLE]

This corresponds to subdiffusion with the transport exponent $0<\beta<1/2$ for $\alpha<0$. For $0<\alpha<1$ subdiffusion is observed as well. Superdiffusion takes place for $1<\alpha<3/2$ since the transprt exponent is $1<\beta<2$. For the case with $3/2<\alpha<2$ one observes hyperdiffusion since $\beta>2$.

2.1 Numerical analysis of MSD

Furthermore, we give the numerically calculated ensemble averaged MSD $\langle x^{2}(t)\rangle$ for the CTRW model of Eq. (7). The waiting time in between successive jumps has PDF $\psi(t)\simeq t^{-3/2}$, which is the waiting time PDF in the fingers of the comb, with a cut-off at $t=1$. Each jump (the increment $\Delta x$ of the actual position $x$) is a sum of a deterministic part (drift) and a random part:

[TABLE]

where, $\delta t$ is some fictitious time step of an integrator and determines the balance between drift and diffusion (in the simulations it is fixed to $\delta t=0.01$), $\xi$ is a gaussian random number with variance = 1. The drift is a consequence of the conversion of the Stratonovich stochastic differential equation Eq. (1) into an Ito stochastic differential equation which is then integrated using the Euler Maruyama scheme. Then the system is run where time is given by the sum of all waiting times. Graphical representation of the obtained MSD is given in Figure 2. There exists good agreement in the long time limit with the analytical results (28).

3 Heterogeneous diffusion on a fractal grid

Considering heterogeneous diffusion on a fractal grid comb, we introduce an infinite number of backbones at positions $y$ which belong to the fractal set $\mathcal{S}_{\nu}$ with fractal dimension $0<\nu<1$ [18]. This situation is described by the diffusion equation

[TABLE]

Such kind of fractal structures is an idealization of more complex comb-like fractal networks in anisotropic porous media [25, 26].

Integration of Eq. (29) over $y$, and then the Laplace transform, yield

[TABLE]

We present the PDF $P(x,y,t)$ in the Laplace space by means the same ansatz (11)

[TABLE]

from where it follows $p_{1}(x,s)=2g(x,s)\sqrt{\frac{\mathcal{D}_{y}}{s}}$. The summation in Eq. (29) is over the fractal set $\mathcal{S}_{\nu}$, which corresponds to integration over the fractal measure $\mu_{\nu}\sim l^{\nu}$, such that $\sum_{l_{j}\in\mathcal{S}_{\nu}}\delta(l-l_{j})\rightarrow\frac{1}{\Gamma(\nu)}t^{\nu-1}$ is the fractal density, and $d\mu_{\nu}=\frac{1}{\Gamma(\nu)}l^{\nu-1}dl$ [27], from where it follows

[TABLE]

From Eqs. (30) and (3), we obtain the equation for the reduced PDF in the closed form

[TABLE]

Here we also take into account symmetrical property of the equation. Therefore we are considering the symmetrical PDF $p_{1}(x,s)=p_{1}(|x|,s)$ and substitution $f(x,s)=|x|^{\alpha/2}p_{1}(x,s)$, and by exchanging $|x|=z$, we find

[TABLE]

from where one find the following system of differential equations:

[TABLE]

[TABLE]

Therefore, for the final results we obtain

[TABLE]

from where the MSD becomes

[TABLE]

Therefore, the diffusion is enhanced in comparison to the heterogeneous diffusion in a comb, as it is expected.

4 Summary

The study is concerned with an inhomogeneous diffusivity of media in the comb-like geometry. We presented exact analytical results for various realizations of anomalous diffusion with a power-law-range-dependent diffusion coefficient $D(x)\sim|x|^{\alpha}$ with $\alpha<2$ along the main channels (backbones) on a comb (with one backbone) and a grid comb structures with a fractal backbone structure. For both cases of the comb and the fractal grid, the exact analytical solutions for the PDFs are obtained in the form of the Fox $H$-function and the mean square displacements are rigorously estimated as well. For the comb, the analytical form of the MSD $\langle x^{2}(t)\rangle\sim t^{1/(2-\alpha)}=t^{\beta}$, in Eq. (28) reflects different realizations of anomalous diffusion, depending on the values of $\alpha$. These regimes are (i) subdiffusion with the transport exponent $0<\beta<1$ for $\alpha<1$; (ii) for $1<\alpha<3/2$ superdiffusion takes place with $1<\beta<2$; (iii) for the case with $3/2<\alpha<2$ one observes hyperdiffusion with the transport exponent $\beta>2$. Another important result is chaotic localization with the negative transport exponent $\beta<0$, when $\alpha>2$. In the case of the fractal grid, the fractal dimension of the backbone structure increases the transport exponent, which depends on the fractal dimension of the structure of the backbones, i.e., $\langle x^{2}(t)\rangle\sim t^{(1+\nu)/(2-\alpha)}$, where $0<\nu<1$ is the fractal dimension of the backbones structure. Here also, all regimes from subdiffusion to hyperdiffusion and localization take place as well.

Acknowledgments

TS acknowledges support within DFG – Deutsche Forschungsgemeinschaft project “Random search processes, Lévy flights, and random walks on complex networks”. TS and AI thank the hospitality at the Max-Planck Institute for the Physics of Complex Systems in Dresden, Germany. AI was also supported by the Israel Science Foundation (ISF).

Appendix A Fox $H$-function

The Fox $H$-function is defined as the inverse Mellin transform for a set of gamma functions [23]

[TABLE]

where

[TABLE]

with $0\leq n\leq p$, $1\leq m\leq q$, $a_{i},b_{j}\in C$, $A_{i},B_{j}\in R^{+}$, $i=1,\ldots,p$, and $j=1,\ldots,q$. The contour $\Omega$ starting at $c-i\infty$ and ending at $c+i\infty$ separates the poles of the function $\Gamma(b_{j}+B_{j}s)$, $j=1,...,m$ from those of the function $\Gamma(1-a_{i}-A_{i}s)$, $i=1,...,n$.

The Mellin-cosine transform of Fox $H$-function is given by [23]

[TABLE]

The Mellin transform of the Fox $H$-function yieds

[TABLE]

where $\theta(\xi)$ is defined in Eqs. (47) and (48).

The inverse Laplace transform of the Fox $H$-function reads [23]

[TABLE]

Fox $H$-function has the following properties [23]

[TABLE]

[TABLE]

Appendix B Lommel equation

The solution of the Lommel differential equation

[TABLE]

is given in terms of the Bessel functions [28]

[TABLE]

The Bessel function $Z_{\frac{1}{2\zeta}}(x)$ is given by $Z_{\frac{1}{2\zeta}}(x)=C_{1}J_{\frac{1}{2\zeta}}(x)+C_{2}N_{\frac{1}{2\zeta}}(x)$, where $J_{\frac{1}{2\zeta}}(x)$ is the Bessel function of the first kind and $N_{\frac{1}{2\zeta}}(x)$ is the Bessel function of the second kind (Neumann function). In our case, the Bessel function is with imaginary argument, therefore the solution of the Lommel equation (74) is given in terms of the modified Bessel function (of the third kind) [28]

[TABLE]

which satisfies the zero boundary conditions at infinity.

The modified Bessel function (of the third kind) $K_{\nu}(z)$ is a special case of the Fox $H$-function [23]

[TABLE]

Its series representation for $z\rightarrow 0$ is given by

[TABLE]