Distribution Functions of Energetic Particles Experiencing Compound Subdiffusion

A comprehensive description of compound subdiffusion was presented in Webb et al. (2006). The latter authors employed an approach based on the so-called Chapman–Kolmogorov equation (see, e.g., Gardiner 1985)

$$\begin{eqnarray}&&{f}_{\perp }(x,y;t)={\int }_{-\infty }^{+\infty }{dz}\,{f}_{\mathrm{FL}}(x,y;z){f}_{\parallel }(z;t),\end{eqnarray} \tag{ 1 }$$

where the particle distribution in the perpendicular direction ${f}_{\perp }(x,y;t)$ is given as the convolution integral of the parallel distribution function ${f}_{\parallel }(z;t)$ and the field line distribution function ${f}_{\mathrm{FL}}(x,y;z)$. This means that particles are assumed to follow field lines and their statistics in the perpendicular direction is entirely controlled by parallel transport and the random walk of magnetic field lines.

All three distribution functions used in Equation (1) can be expressed by a Fourier representation so that

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\parallel }(z;t) & = & {\int }_{-\infty }^{+\infty }{{dk}}_{\parallel }\,{f}_{\parallel }({k}_{\parallel };t){e}^{{{ik}}_{\parallel }z},\\ {f}_{\mathrm{FL}}(x,y;z) & = & \int {d}^{2}k\,{f}_{\mathrm{FL}}({k}_{x},{k}_{y};z){e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }},\\ {f}_{\perp }(x,y;t) & = & \int {d}^{2}k\,{f}_{\perp }({k}_{x},{k}_{y};t){e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }},\end{array}\end{eqnarray} \tag{ 2 }$$

where the last two integrals cover the whole wave number space in x- and y-directions. Based on Equation (1), we find for the Fourier transform of the perpendicular distribution function

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }({k}_{x},{k}_{y};t) & = & {\int }_{-\infty }^{+\infty }{dz}\,{f}_{\mathrm{FL}}({k}_{x},{k}_{y};z){f}_{\parallel }(z;t)\\ & = & {\int }_{-\infty }^{+\infty }{{dk}}_{\parallel }\,{f}_{\parallel }({k}_{\parallel };t)\\ & & \times \ {\int }_{-\infty }^{+\infty }{dz}\,{f}_{\mathrm{FL}}({k}_{x},{k}_{y};z){e}^{{{ik}}_{\parallel }z}.\end{array}\end{eqnarray} \tag{ 3 }$$

This is evaluated for certain parallel transport and field line random walk models below.

It is also crucial to think about initial conditions. First we derive from the third line of Equation (2)

$$\begin{eqnarray}&&{f}_{\perp }(x,y;t=0)=\int {d}^{2}k\,{f}_{\perp }({k}_{x},{k}_{y};t=0){e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\end{eqnarray} \tag{ 4 }$$

and assume

$$\begin{eqnarray}&&{f}_{\perp }({k}_{x},{k}_{y};t=0)=\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\end{eqnarray} \tag{ 5 }$$

so that

$$\begin{eqnarray}&&{f}_{\perp }(x,y;t=0)=\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\int {d}^{2}k\,{e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}=\delta \left(x\right)\delta \left(y\right)\end{eqnarray} \tag{ 6 }$$

where we have used (see, e.g., Zwillinger 2012)

$$\begin{eqnarray}&&{\int }_{-\infty }^{+\infty }{dx}\,{e}^{{ikx}}=2\pi \delta \left(k\right)\end{eqnarray} \tag{ 7 }$$

twice. This means that any obtained distribution function needs to satisfy the constraint given by Equation (5). Only then the correct initial condition given by Equation (6) is satisfied.

The nth moment of a distribution is defined via

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle ={\int }_{-\infty }^{+\infty }{dx}\,{\int }_{-\infty }^{+\infty }{dy}\,{x}^{n}{f}_{\perp }(x,y;t).\end{eqnarray} \tag{ 8 }$$

Due to symmetry we find $\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle =0$ if n is odd. Thus, we focus on the case of even n. In order to rewrite this we follow Shalchi & Gammon (2019) and perform the following steps starting by combining Equation (8) with the third line of Equation (2)

$$\begin{eqnarray}&&\begin{array}{l}\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle \,=\,\int {d}^{2}k\,{f}_{\perp }({k}_{x},{k}_{y};t){\int }_{-\infty }^{+\infty }{dx}\,{\int }_{-\infty }^{+\infty }{dy}\,{x}^{n}{e}^{{{ik}}_{x}x+{{ik}}_{y}y}\\ \,=\,{\left(-i\right)}^{n}\int {d}^{2}k\,{f}_{\perp }({k}_{x},{k}_{y};t)\\ \,\times \ \displaystyle \frac{{d}^{n}}{{{dk}}_{x}^{n}}{\int }_{-\infty }^{+\infty }{dx}\,{\int }_{-\infty }^{+\infty }{dy}\,{e}^{{{ik}}_{x}x+{{ik}}_{y}y}\\ \,=\,{\left(2\pi \right)}^{2}{\left(-i\right)}^{n}\int {d}^{2}k\,{f}_{\perp }({k}_{x},{k}_{y};t)\delta \left({k}_{y}\right)\displaystyle \frac{{d}^{n}}{{{dk}}_{x}^{n}}\delta \left({k}_{x}\right)\\ \,=\,{\left(2\pi \right)}^{2}{i}^{n}{\left[\displaystyle \frac{{d}^{n}{f}_{\perp }({k}_{x},{k}_{y};t)}{{{dk}}_{x}^{n}}\right]}_{{k}_{x}={k}_{y}=0}\end{array}\end{eqnarray} \tag{ 9 }$$

where we have used Equation (7) and integration by parts n times. Equation (9) can be used for n = 0 to obtain the normalization condition

$$\begin{eqnarray}&&1={\left(2\pi \right)}^{2}{f}_{\perp }\left({k}_{x}={k}_{y}=0;t\right).\end{eqnarray} \tag{ 10 }$$

Any type of distribution function considered in the current article needs to satisfy this condition as well as the initial condition (5).

It is important to understand the relation between characteristic functions and Fourier transforms of distribution functions. First we note that the characteristic function is defined via

$$\begin{eqnarray}&&\left\langle {e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\right\rangle := \int {d}^{n}x\,f\left({\boldsymbol{x}};t\right){e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}},\end{eqnarray} \tag{ 11 }$$

where the integer number n denotes the dimensionality of the problem. Comparing this with the inverse Fourier transform

$$\begin{eqnarray}&&f\left({\boldsymbol{k}};t\right)=\displaystyle \frac{1}{{\left(2\pi \right)}^{n}}\int {d}^{n}x\,f\left({\boldsymbol{x}};t\right){e}^{-i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\end{eqnarray} \tag{ 12 }$$

yields

$$\begin{eqnarray}&&\left\langle {e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\right\rangle ={\left(2\pi \right)}^{n}f\left({\boldsymbol{k}};t\right).\end{eqnarray} \tag{ 13 }$$

Please note that here we assumed that all distribution functions are symmetric in the spectrum

$$\begin{eqnarray}&&f\left({\boldsymbol{k}};t\right)=f\left(-{\boldsymbol{k}};t\right)\equiv {f}^{* }\left({\boldsymbol{k}};t\right).\end{eqnarray} \tag{ 14 }$$

According to Equation (13), characteristic functions and Fourier transforms are directly proportional. The knowledge of characteristic functions is essential for the development of analytical theories for perpendicular diffusion (see, e.g., Lasuik & Shalchi 2017, 2018; Shalchi 2017).

For the axisymmetric case the perpendicular particle distribution function in configuration space can be written as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }(\rho ;t) & = & \int {d}^{2}k\,{f}_{\perp }({k}_{\perp };t){e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\\ & = & {\int }_{0}^{\infty }{dk}_{\perp }\,{k}_{\perp }{f}_{\perp }({k}_{\perp };t)\,{\int }_{0}^{2\pi }d{\rm{\Psi }}\,{e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\end{array}\end{eqnarray} \tag{ 15 }$$

where we switched from Cartesian to polar coordinates by using

$$\begin{eqnarray}\begin{array}{rcl}\rho & = & \sqrt{{x}^{2}+{y}^{2}},\quad \quad {\rm{\Phi }}=\arctan \left(\displaystyle \frac{y}{x}\right),\\ {k}_{\perp } & = & \sqrt{{k}_{x}^{2}+{k}_{y}^{2}},\quad \quad {\rm{\Psi }}=\arctan \left(\displaystyle \frac{{k}_{y}}{{k}_{x}}\right).\end{array}\end{eqnarray} \tag{ 16 }$$

In order to evaluate this further, we can employ

$$\begin{eqnarray}&&{e}^{i{{\boldsymbol{x}}}_{\perp }\cdot {{\boldsymbol{k}}}_{\perp }}={e}^{i\rho {k}_{\perp }\cos ({\rm{\Psi }}-{\rm{\Phi }})}.\end{eqnarray} \tag{ 17 }$$

However, the resulting distribution function in Equation (15) must not depend on the angle Φ due to symmetry. Therefore, we can set ${\rm{\Phi }}=0$ so that

$$\begin{eqnarray}\begin{array}{rcl}{e}^{i{{\boldsymbol{x}}}_{\perp }\cdot {{\boldsymbol{k}}}_{\perp }} & = & {e}^{i\rho {k}_{\perp }\cos ({\rm{\Psi }})}={e}^{i\rho {k}_{\perp }\sin ({\rm{\Psi }}+\pi /2)}\\ & = & \displaystyle \sum _{n=-\infty }^{+\infty }{J}_{n}(\rho {k}_{\perp }){e}^{{in}({\rm{\Psi }}+\pi /2)},\end{array}\end{eqnarray} \tag{ 18 }$$

where we have used Bessel functions of first kind J_n and the so-called Jacobi–Anger expansion (see, e.g., Cuyt et al. 2008). If we integrate the exponential factor given by Equation (18) over Ψ, we can employ

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }d{\rm{\Psi }}\,{e}^{{in}{\rm{\Psi }}}=2\pi {\delta }_{n0}\end{eqnarray} \tag{ 19 }$$

so that

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }d{\rm{\Psi }}\,{e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}=2\pi {J}_{0}(\rho {k}_{\perp }).\end{eqnarray} \tag{ 20 }$$

Using this in Equation (15) yields

$$\begin{eqnarray}&&{f}_{\perp }(\rho ;t)=2\pi {\int }_{0}^{\infty }{{dk}}_{\perp }\,{k}_{\perp }{f}_{\perp }({k}_{\perp };t){J}_{0}(\rho {k}_{\perp }),\end{eqnarray} \tag{ 21 }$$

which will be used later in this article when specific forms of the function ${f}_{\perp }({k}_{\perp };t)$ are considered. Furthermore, we can also compute the distribution function projected onto the x-axis defined via

$$\begin{eqnarray}&&{f}_{\perp }(x;t)={\int }_{-\infty }^{+\infty }{dy}\,{f}_{\perp }(x,y;t).\end{eqnarray} \tag{ 22 }$$

Using the third line of Equation (2) therein yields

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }(x;t) & = & \int {d}^{2}k\,{f}_{\perp }({k}_{x},{k}_{y};t){e}^{{{ik}}_{x}x}{\int }_{-\infty }^{+\infty }{dy}\,{e}^{{{ik}}_{y}y}\\ & = & 2\pi \int {d}^{2}k\,{f}_{\perp }({k}_{x},{k}_{y};t){e}^{{{ik}}_{x}x}\delta \left({k}_{y}\right)\\ & = & 2\pi {\int }_{-\infty }^{+\infty }{{dk}}_{x}\,{f}_{\perp }({k}_{x},{k}_{y}=0;t){e}^{{{ik}}_{x}x}\\ & = & 4\pi {\int }_{0}^{\infty }{{dk}}_{\perp }\,{f}_{\perp }({k}_{\perp };t)\cos ({k}_{\perp }x).\end{array}\end{eqnarray} \tag{ 23 }$$

The same result can be obtained if the projection on the y-axis is considered due to symmetry. So far we discussed some general relations. In order to obtain concrete results, we have to specify field line and parallel distribution functions. This is done in the next section.

If the field line random walk is diffusive and assuming a Gaussian distribution of field lines, the field line distribution function is given by

$$\begin{eqnarray}&&{f}_{\mathrm{FL}}(x,y;t)=\displaystyle \frac{1}{4\pi {\kappa }_{\mathrm{FL}}\left|z\right|}{e}^{-({x}^{2}+{y}^{2})/(4{\kappa }_{\mathrm{FL}}\left|z\right|)}.\end{eqnarray} \tag{ 24 }$$

In this case the Fourier transform of the field line distribution function is

$$\begin{eqnarray}&&{f}_{\mathrm{FL}}({k}_{\perp };z)=\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{e}^{-{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\left|z\right|}.\end{eqnarray} \tag{ 25 }$$

The characteristic function, on the other hand, is given by

$$\begin{eqnarray}&&{\left\langle {e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\right\rangle }_{\mathrm{FL}}={e}^{-{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\left|z\right|}.\end{eqnarray} \tag{ 26 }$$

Thus we find for the perpendicular particle distribution function given by Equation (3)

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }({k}_{\perp };t) & = & \displaystyle \frac{1}{2{\pi }^{2}}{\int }_{-\infty }^{+\infty }{{dk}}_{\parallel }\,{f}_{\parallel }({k}_{\parallel };t)\\ & & \times \ {\int }_{0}^{\infty }{dz}\,\cos \left({k}_{\parallel }z\right){e}^{-{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\left|z\right|}.\end{array}\end{eqnarray} \tag{ 27 }$$

Using

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{dz}\,\cos \left({az}\right){e}^{-{bz}}=\displaystyle \frac{b}{{a}^{2}+{b}^{2}}\end{eqnarray} \tag{ 28 }$$

yields

$$\begin{eqnarray}&&{f}_{\perp }({k}_{\perp };t)=\displaystyle \frac{1}{2{\pi }^{2}}{\int }_{-\infty }^{+\infty }{{dk}}_{\parallel }\,{f}_{\parallel }({k}_{\parallel };t)\displaystyle \frac{{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}}{{\left({\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\right)}^{2}+{k}_{\parallel }^{2}}.\end{eqnarray} \tag{ 29 }$$

If we also assume that parallel transport is diffusive and if we employ a Gaussian particle distribution function

$$\begin{eqnarray}&&{f}_{\parallel }(z;t)=\displaystyle \frac{1}{\sqrt{4\pi {\kappa }_{\parallel }t}}{e}^{-{z}^{2}/(4{\kappa }_{\parallel }t)}\end{eqnarray} \tag{ 30 }$$

we obtain

$$\begin{eqnarray}&&{f}_{\parallel }({k}_{\parallel };t)=\displaystyle \frac{1}{2\pi }{e}^{-{\kappa }_{\parallel }{k}_{\parallel }^{2}t}\end{eqnarray} \tag{ 31 }$$

yielding the characteristic function

$$\begin{eqnarray}&&\left\langle {e}^{{{izk}}_{\parallel }}\right\rangle ={e}^{-{\kappa }_{\parallel }{k}_{\parallel }^{2}t}.\end{eqnarray} \tag{ 32 }$$

Thus, Equation (29) turns into

$$\begin{eqnarray}&&{f}_{\perp }({k}_{\perp };t)=\displaystyle \frac{1}{2{\pi }^{3}}{\int }_{0}^{\infty }{{dk}}_{\parallel }\,\displaystyle \frac{{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}}{{\left({\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\right)}^{2}+{k}_{\parallel }^{2}}{e}^{-{\kappa }_{\parallel }{k}_{\parallel }^{2}t}.\end{eqnarray} \tag{ 33 }$$

The occurring integral can be solved by (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{dx}\,\displaystyle \frac{{e}^{-{{bx}}^{2}}}{{a}^{2}+{x}^{2}}=\displaystyle \frac{\pi }{2a}{e}^{{a}^{2}b}{\rm{erfc}}\left(a\sqrt{b}\right),\end{eqnarray} \tag{ 34 }$$

where we have used the complementary error function. Finally we obtain

$$\begin{eqnarray}&&{f}_{\perp }({k}_{\perp };t)=\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right)\end{eqnarray} \tag{ 35 }$$

where we have used

$$\begin{eqnarray}&&\xi ={\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\sqrt{{\kappa }_{\parallel }t}.\end{eqnarray} \tag{ 36 }$$

Equation (35) is in agreement with Equation (B5) of Webb et al. (2006). Furthermore, we can easily understand that the characteristic function in the perpendicular direction is now given by

$$\begin{eqnarray}&&\left\langle {e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\right\rangle ={e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right).\end{eqnarray} \tag{ 37 }$$

Typically, if particle transport is diffusive, we would expect a Gaussian distribution such as the characteristic functions given by Equations (26) and (32). Clearly this is not obtained for compound subdiffusive transport. Note that the scaled complementary error function is defined via

$$\begin{eqnarray}&&{\rm{erfcx}}\left(x\right)={\rm{erfc}}\left(x\right){e}^{{x}^{2}}.\end{eqnarray} \tag{ 38 }$$

Thus, we can simply write

$$\begin{eqnarray}&&{f}_{\perp }({k}_{\perp };t)=\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\rm{erfcx}}\left(\xi \right)\end{eqnarray} \tag{ 39 }$$

as well as

$$\begin{eqnarray}&&\left\langle {e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\right\rangle ={\rm{erfcx}}\left(\xi \right).\end{eqnarray} \tag{ 40 }$$

These forms are in particular useful for a numerical evaluation in order to avoid arithmetic underflow.

In order to determine the distribution function in configuration space we can use different approaches. One way of performing this task is to assume Gaussian distributions of field lines and particles in the parallel direction. Using Equations (24) and (30) in (1) yields

$$\begin{eqnarray}&&{f}_{\perp }(x,y;t)=\displaystyle \frac{1}{\sqrt{4\pi {\kappa }_{\parallel }t}}{\int }_{-\infty }^{+\infty }{dz}\,\displaystyle \frac{1}{4\pi {\kappa }_{\mathrm{FL}}\left|z\right|}{e}^{-\tfrac{{z}^{2}}{4{\kappa }_{\parallel }t}-\displaystyle \frac{{x}^{2}+{y}^{2}}{4{\kappa }_{\mathrm{FL}}\left|z\right|}}.\end{eqnarray} \tag{ 41 }$$

After the integral transformation λ = z/$\sqrt{4{\kappa }_{\parallel }t}$, this becomes

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }(\rho ;t) & = & \displaystyle \frac{2}{{\pi }^{3/2}}\displaystyle \frac{1}{4{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}}\\ & & \times \ {\int }_{0}^{\infty }\displaystyle \frac{d\lambda }{\lambda }\,{e}^{-{\lambda }^{2}-{\rho }^{2}/\left(\lambda 4{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}\right)}.\end{array}\end{eqnarray} \tag{ 42 }$$

Apart from the variable ρ itself, the distribution function depends only on one single parameter and this is $4{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}$. The remaining integral can be expressed via the Meijer G-function (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}&&{\int }_{0}^{\infty }\displaystyle \frac{{dx}}{x}\,{e}^{-{x}^{2}-a/x}=\displaystyle \frac{1}{2\sqrt{\pi }}{G}_{0,3}^{3,0}\left(\left.\begin{array}{c}-\\ 0,0,1/2\end{array}\right|\displaystyle \frac{{a}^{2}}{4}\right)\end{eqnarray} \tag{ 43 }$$

so that

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }(\rho ;t) & = & \displaystyle \frac{1}{{\pi }^{2}}\displaystyle \frac{1}{4{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}}\\ & & \times \ {G}_{0,3}^{3,0}\left(\left.\begin{array}{c}-\\ 0,0,1/2\end{array}\right|\displaystyle \frac{{\rho }^{4}}{256{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t}\right).\end{array}\end{eqnarray} \tag{ 44 }$$

For the projected function we can use

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }(x;t) & = & {\int }_{-\infty }^{+\infty }{dy}\,{f}_{\perp }(x,y;t)\\ & = & \displaystyle \frac{1}{\sqrt{4\pi {\kappa }_{\parallel }t}}{\int }_{-\infty }^{+\infty }{dz}\,\displaystyle \frac{1}{4\pi {\kappa }_{\mathrm{FL}}\left|z\right|}{e}^{-\tfrac{{z}^{2}}{4{\kappa }_{\parallel }t}-\tfrac{{x}^{2}}{4{\kappa }_{\mathrm{FL}}\left|z\right|}}\\ & & \times \ {\int }_{-\infty }^{+\infty }{dy}\,{e}^{-\tfrac{{y}^{2}}{4{\kappa }_{\mathrm{FL}}\left|z\right|}}\\ & = & \displaystyle \frac{1}{\sqrt{4\pi {\kappa }_{\parallel }t}}{\int }_{-\infty }^{+\infty }{dz}\,\displaystyle \frac{1}{\sqrt{4\pi {\kappa }_{\mathrm{FL}}\left|z\right|}}{e}^{-\tfrac{{z}^{2}}{4{\kappa }_{\parallel }t}-\tfrac{{x}^{2}}{4{\kappa }_{\mathrm{FL}}\left|z\right|}}.\end{array}\end{eqnarray} \tag{ 45 }$$

After using again the integral transformation $\lambda =z/\sqrt{4{\kappa }_{\parallel }t}$ this becomes

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }(x;t) & = & \displaystyle \frac{1}{\pi }\displaystyle \frac{1}{\sqrt{{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}}}\\ & & \times \ {\int }_{0}^{\infty }\displaystyle \frac{d\lambda }{\sqrt{\lambda }}\,{e}^{-{\lambda }^{2}-{x}^{2}/\left(\lambda 4{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}\right)}.\end{array}\end{eqnarray} \tag{ 46 }$$

The remaining integral can be expressed via generalized hypergeometric functions (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}\begin{array}{rcl}{\int }_{0}^{\infty }\displaystyle \frac{{dx}}{\sqrt{x}}\,{e}^{-{x}^{2}-a/x} & = & 2{\rm{\Gamma }}\left(\displaystyle \frac{5}{4}\right){}_{0}{F}_{2}\left(;\displaystyle \frac{1}{2},\displaystyle \frac{3}{4};-\displaystyle \frac{{a}^{2}}{4}\right)\\ & & -\ 2\sqrt{\pi a}{}_{0}{F}_{2}\left(;\displaystyle \frac{3}{4},\displaystyle \frac{5}{4};-\displaystyle \frac{{a}^{2}}{4}\right)\\ & & +\ 2a{\rm{\Gamma }}\left(\displaystyle \frac{3}{4}\right){}_{0}{F}_{2}\left(;\displaystyle \frac{5}{4},\displaystyle \frac{3}{2};-\displaystyle \frac{{a}^{2}}{4}\right).\end{array}\end{eqnarray} \tag{ 47 }$$

Using this in Equation (46) allows us to write

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\perp }(x;t)=\displaystyle \frac{2}{\pi }\displaystyle \frac{1}{\sqrt{{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}}}\\ \quad \times \ \left[{\rm{\Gamma }}\left(\displaystyle \frac{5}{4}\right){}_{0}{F}_{2}\left(;\displaystyle \frac{1}{2},\displaystyle \frac{3}{4};-\displaystyle \frac{{x}^{4}}{256{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t}\right)\right.\\ \quad -\ \sqrt{\displaystyle \frac{\pi {x}^{2}}{4{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}}}{}_{0}{F}_{2}\left(;\displaystyle \frac{3}{4},\displaystyle \frac{5}{4};-\displaystyle \frac{{x}^{4}}{256{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t}\right)\\ \quad +\ \left.\displaystyle \frac{{x}^{2}}{4{\kappa }_{\mathrm{FL}}\sqrt{4{\kappa }_{\parallel }t}}{\rm{\Gamma }}\left(\displaystyle \frac{3}{4}\right){}_{0}{F}_{2}\left(;\displaystyle \frac{5}{4},\displaystyle \frac{3}{2};-\displaystyle \frac{{x}^{4}}{256{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t}\right)\right].\end{array}\end{eqnarray} \tag{ 48 }$$

Problematic here is that we found solutions depending on more complicated special functions such as the Meijer G-function and generalized hypergeometric functions.

Alternatively, the perpendicular particle distribution function in configuration space is given by Equation (21). With Equation (35) this becomes

$$\begin{eqnarray}&&{f}_{\perp }(\rho ;t)=\displaystyle \frac{1}{2\pi }{\int }_{0}^{\infty }{{dk}}_{\perp }\,{k}_{\perp }{J}_{0}(\rho {k}_{\perp }){e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right),\end{eqnarray} \tag{ 49 }$$

which is in agreement with Equation (3.26) of Webb et al. (2006). The projected distribution function is given by Equation (23). For our case this becomes

$$\begin{eqnarray}&&{f}_{\perp }(x;t)=\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }{{dk}}_{\perp }\,\cos ({k}_{\perp }x){e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right).\end{eqnarray} \tag{ 50 }$$

However, Equations (49) and (50) are difficult to evaluate without employing approximations. Besides a numerical evaluation one can employ different approximations which is done in Section 3. At this point we only compute the values of these two functions at the origin. For x = 0, Equation (50) simplifies to (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }(x=0;t) & = & \displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }{{dk}}_{\perp }\,{e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right)\\ & = & \displaystyle \frac{\sqrt{2}}{\pi \sqrt{{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}{\rm{\Gamma }}\left(\displaystyle \frac{5}{4}\right).\end{array}\end{eqnarray} \tag{ 51 }$$

Employing (see, e.g., Abramowitz & Stegun 1974)

$$\begin{eqnarray}&&{\rm{\Gamma }}\left(z+1\right)=z{\rm{\Gamma }}\left(z\right)\end{eqnarray} \tag{ 52 }$$

allows us to write this as

$$\begin{eqnarray}&&{f}_{\perp }(x=0;t)=\displaystyle \frac{1}{\pi \sqrt{8{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}{\rm{\Gamma }}\left(\displaystyle \frac{1}{4}\right).\end{eqnarray} \tag{ 53 }$$

According to Abramowitz & Stegun (1974) we have ${\rm{\Gamma }}(1/4)\,\approx 3.62561$ so that

$$\begin{eqnarray}&&{f}_{\perp }(x=0;t)\approx \displaystyle \frac{0.408}{\sqrt{{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}.\end{eqnarray} \tag{ 54 }$$

Although this formula only provides the value at the origin, it is useful because it is exact. For $\rho =0$, on the other hand, Equation (49) turns into

$$\begin{eqnarray}&&{f}_{\perp }(\rho =0;t)=\displaystyle \frac{1}{2\pi }{\int }_{0}^{\infty }{{dk}}_{\perp }\,{k}_{\perp }{e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right),\end{eqnarray} \tag{ 55 }$$

which does not converge due to the behavior of the integrand at large k_⊥. Therefore, ${f}_{\perp }(\rho \to 0;t)\to \infty$.

One can easily show that solution (35) satisfies the normalization condition (10). Furthermore, by combining Equation (35) with (9) for n = 2, we can compute the second moment. The second derivative can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{f}_{\perp }}{\partial {k}_{\perp }^{2}}=\displaystyle \frac{{\partial }^{2}{f}_{\perp }}{\partial {\xi }^{2}}{\left(\displaystyle \frac{\partial \xi }{\partial {k}_{\perp }}\right)}^{2}+\displaystyle \frac{\partial {f}_{\perp }}{\partial \xi }\displaystyle \frac{{\partial }^{2}\xi }{\partial {k}_{\perp }^{2}}.\end{eqnarray} \tag{ 56 }$$

Since $\partial \xi /\partial {k}_{\perp }\propto {k}_{\perp }$ and we set ${k}_{x}={k}_{y}=0$ in Equation (9), the first term goes to zero. In the second term we can use Equation (36) to obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\xi }{\partial {k}_{\perp }^{2}}=2{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}\end{eqnarray} \tag{ 57 }$$

so that

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{2}\rangle =-{\left(2\pi \right)}^{2}2{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}{\left[\displaystyle \frac{{{df}}_{\perp }(\xi ;t)}{d\xi }\right]}_{\xi =0}.\end{eqnarray} \tag{ 58 }$$

Using (see, e.g., Abramowitz & Stegun 1974)

$$\begin{eqnarray}&&\displaystyle \frac{d}{d\xi }{\rm{erfc}}\left(\xi \right)=-\displaystyle \frac{2}{\sqrt{\pi }}{e}^{-{\xi }^{2}}\end{eqnarray} \tag{ 59 }$$

and Equation (35) yields

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{df}}_{\perp }(\xi )}{d\xi } & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\left[2\xi {e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right)+{e}^{{\xi }^{2}}\displaystyle \frac{d}{d\xi }{\rm{erfc}}\left(\xi \right)\right]\\ & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\left[2\xi {e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right)-\displaystyle \frac{2}{\sqrt{\pi }}\right].\end{array}\end{eqnarray} \tag{ 60 }$$

Using this in Equation (58) provides us with

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{2}\rangle =4{\kappa }_{\mathrm{FL}}\sqrt{\tfrac{{\kappa }_{\parallel }t}{\pi }}.\end{eqnarray} \tag{ 61 }$$

The latter formula is well-known and was derived in several previous papers (see, e.g., Kóta & Jokipii 2000). A time-dependent or running diffusion coefficient can be defined via

$$\begin{eqnarray}&&{d}_{\perp }(t):= \displaystyle \frac{1}{2}\displaystyle \frac{d}{{dt}}\langle {\left({\rm{\Delta }}x\right)}^{2}\rangle ,\end{eqnarray} \tag{ 62 }$$

which becomes for the mean square displacement given by Equation (61)

$$\begin{eqnarray}&&{d}_{\perp }(t)={\kappa }_{\mathrm{FL}}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}.\end{eqnarray} \tag{ 63 }$$

Very clearly we find a declining diffusion coefficient corresponding to subdiffusion. This type of behavior can also be observed in simulations (see, e.g., Qin et al. 2002a). We can also compute the fourth moment. In order to do this we need the following derivatives

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\partial }^{3}{f}_{\perp }}{\partial {k}_{\perp }^{3}} & = & \displaystyle \frac{{\partial }^{3}{f}_{\perp }}{\partial {\xi }^{3}}{\left(\displaystyle \frac{\partial \xi }{\partial {k}_{\perp }}\right)}^{3}+3\displaystyle \frac{{\partial }^{2}{f}_{\perp }}{\partial {\xi }^{2}}\displaystyle \frac{{\partial }^{2}\xi }{\partial {k}_{\perp }^{2}}\displaystyle \frac{\partial \xi }{\partial {k}_{\perp }}+\ \displaystyle \frac{\partial {f}_{\perp }}{\partial \xi }\displaystyle \frac{{\partial }^{3}\xi }{\partial {k}_{\perp }^{3}}\end{array}\end{eqnarray} \tag{ 64 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {k}_{\perp }^{4}} & = & \displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {\xi }^{4}}{\left(\displaystyle \frac{\partial \xi }{\partial {k}_{\perp }}\right)}^{4}+6\displaystyle \frac{{\partial }^{3}{f}_{\perp }}{\partial {\xi }^{3}}{\left(\displaystyle \frac{\partial \xi }{\partial {k}_{\perp }}\right)}^{2}\displaystyle \frac{{\partial }^{2}\xi }{\partial {k}_{\perp }^{2}}\\ & & +\ 4\displaystyle \frac{{\partial }^{2}{f}_{\perp }}{\partial {\xi }^{2}}\displaystyle \frac{{\partial }^{3}\xi }{\partial {k}_{\perp }^{3}}\displaystyle \frac{\partial \xi }{\partial {k}_{\perp }}+3\displaystyle \frac{{\partial }^{2}{f}_{\perp }}{\partial {\xi }^{2}}{\left(\displaystyle \frac{{\partial }^{2}\xi }{\partial {k}_{\perp }^{2}}\right)}^{2}\\ & & +\ \displaystyle \frac{\partial {f}_{\perp }}{\partial \xi }\displaystyle \frac{{\partial }^{4}\xi }{\partial {k}_{\perp }^{4}}.\end{array}\end{eqnarray} \tag{ 65 }$$

In the limit ${k}_{\perp }\to 0$ the only nonvanishing derivative of ξ with respect to k_⊥ is the one given by Equation (57). This means that all terms except the fourth one on the right-hand side of Equation (65) vanish. Thus, we find

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{d}^{4}{f}_{\perp }}{{{dk}}_{x}^{4}}\right]}_{{k}_{x}={k}_{y}=0}=12{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t{\left[\displaystyle \frac{{d}^{2}{f}_{\perp }(\xi ;t)}{d{\xi }^{2}}\right]}_{\xi =0}.\end{eqnarray} \tag{ 66 }$$

Furthermore, we find with the help of Equation (59)

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{f}_{\perp }(\xi )}{d{\xi }^{2}}=\displaystyle \frac{1}{2{\pi }^{2}}\left[\left(1+2{\xi }^{2}\right){e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right)-\displaystyle \frac{2}{\sqrt{\pi }}\xi \right]\end{eqnarray} \tag{ 67 }$$

and, thus,

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{d}^{2}{f}_{\perp }(\xi ;t)}{d{\xi }^{2}}\right]}_{\xi =0}=\displaystyle \frac{1}{2{\pi }^{2}}\end{eqnarray} \tag{ 68 }$$

where we have used ${\rm{erfc}}(0)=1$. Combining these findings with Equation (9) yields

$$\begin{eqnarray}\begin{array}{rcl}\langle {\left({\rm{\Delta }}x\right)}^{4}\rangle & = & {\left(2\pi \right)}^{2}{i}^{4}{\left[\displaystyle \frac{{d}^{4}{f}_{\perp }({k}_{x},{k}_{y};t)}{{{dk}}_{x}^{4}}\right]}_{{k}_{x}={k}_{y}=0}\\ & = & 24{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t.\end{array}\end{eqnarray} \tag{ 69 }$$

Note that this result is in perfect agreement with Equation (2.24) of Webb et al. (2006). A general formula for the nth moment is derived in Appendix A. By using two different methods, it is shown that

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle =\displaystyle \frac{1}{\pi }{\rm{\Gamma }}\left(\displaystyle \frac{n+1}{2}\right){\rm{\Gamma }}\left(\displaystyle \frac{n+2}{4}\right){\left(8{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}\right)}^{n/2}\end{eqnarray} \tag{ 70 }$$

for even n. For odd n we have $\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle =0$ as mentioned before.

The exact relations derived so far were based on the Chapman–Kolmogorov approach as given by Equation (1). Naturally, the following question arises: what is the corresponding transport equation? In the Fourier space a usual diffusion equation for the axisymmetric case has the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{\perp }}{\partial t}=-{\kappa }_{\perp }{k}_{\perp }^{2}{f}_{\perp }=-{\kappa }_{\perp }\left({k}_{x}^{2}+{k}_{y}^{2}\right){f}_{\perp },\end{eqnarray} \tag{ 71 }$$

where we have used the perpendicular diffusion coefficient ${\kappa }_{\perp }$. The corresponding transport equation in configuration space can easily be obtained via the formal replacement

$$\begin{eqnarray}&&{k}_{n}\to -i\displaystyle \frac{\partial }{\partial {x}_{n}}\end{eqnarray} \tag{ 72 }$$

so that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{\perp }}{\partial t}={\kappa }_{\perp }{{\rm{\Delta }}}_{\perp }{f}_{\perp }\end{eqnarray} \tag{ 73 }$$

where we have used the Laplace operator

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\perp }=\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {y}^{2}}.\end{eqnarray} \tag{ 74 }$$

However, in the case considered in the current paper, the transport is subdiffusive and, therefore, the transport equation needs to have a different form. One way of approaching this problem is to assume that the transport equation still contains only a first order time-derivative. This can be justified because the distribution function has to be determined in the sense that if the initial distribution is known, one has to be able to compute the distribution for all later times. Considering the first time-derivative of Equation (35) yields

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {f}_{\perp }}{\partial t} & = & \displaystyle \frac{\partial f}{\partial \xi }\displaystyle \frac{\partial \xi }{\partial t}=\displaystyle \frac{1}{8{\pi }^{2}}{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{t}}\left[2\xi {e}^{{\xi }^{2}}{\rm{erfc}}\left(\xi \right)-\displaystyle \frac{2}{\sqrt{\pi }}\right],\end{array}\end{eqnarray} \tag{ 75 }$$

where we have also used Equation (60). Using the distribution function therein allows us to derive

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{\perp }}{\partial t}={\kappa }_{\parallel }{\kappa }_{\mathrm{FL}}^{2}{k}_{\perp }^{4}{f}_{\perp }+g({k}_{\perp },t),\end{eqnarray} \tag{ 76 }$$

where we have used the inhomogeneity

$$\begin{eqnarray}&&g({k}_{\perp },t)=-\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}.\end{eqnarray} \tag{ 77 }$$

Transport Equation (76) is clearly in disagreement with Equation (71). The next step is to determine the function $g(x,y,t)$ in configuration space. We derive

$$\begin{eqnarray}\begin{array}{rcl}g(x,y,t) & = & \int {d}^{2}k\,g({k}_{\perp },t){e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\\ & = & -\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\kappa }_{\mathrm{FL}}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}\int {d}^{2}k\,\left({k}_{x}^{2}+{k}_{y}^{2}\right){e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\\ & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\kappa }_{\mathrm{FL}}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}\int {d}^{2}k\,\left(\displaystyle \frac{{d}^{2}}{{{dx}}^{2}}+\displaystyle \frac{{d}^{2}}{{{dy}}^{2}}\right){e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\\ & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\kappa }_{\mathrm{FL}}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}\left(\displaystyle \frac{{d}^{2}}{{{dx}}^{2}}+\displaystyle \frac{{d}^{2}}{{{dy}}^{2}}\right)\int {d}^{2}k\,{e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\\ & = & {\kappa }_{\mathrm{FL}}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}\left[\delta (x)\displaystyle \frac{{d}^{2}}{{{dy}}^{2}}\delta (y)+\delta (y)\displaystyle \frac{{d}^{2}}{{{dx}}^{2}}\delta (x)\right],\end{array}\end{eqnarray} \tag{ 78 }$$

where we have used relation (7) several times.

Using the operator replacement given by Equation (72) in Equation (76), yields

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{\perp }}{\partial t}={\kappa }_{\parallel }{\kappa }_{\mathrm{FL}}^{2}\left(\displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {x}^{4}}+2\displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {x}^{2}\partial {y}^{2}}+\displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {y}^{4}}\right)+g(x,y,t)\end{eqnarray} \tag{ 79 }$$

where the function $g(x,y,t)$ is given by the last line of Equation (78). In Appendix C a more detailed discussion of this differential equation and its solution can be found.

Alternatively, we can use the relation

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dx}}\delta (x)=-\displaystyle \frac{1}{x}\delta (x)\end{eqnarray} \tag{ 80 }$$

to derive

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}}{{{dx}}^{2}}\delta (x)=\displaystyle \frac{2}{{x}^{2}}\delta (x).\end{eqnarray} \tag{ 81 }$$

Therewith we can derive the following transport equation

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {f}_{\perp }}{\partial t} & = & {\kappa }_{\parallel }{\kappa }_{\mathrm{FL}}^{2}\left(\displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {x}^{4}}+2\displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {x}^{2}\partial {y}^{2}}+\displaystyle \frac{{\partial }^{4}{f}_{\perp }}{\partial {y}^{4}}\right)\\ & & +\ 2{\kappa }_{\mathrm{FL}}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}\left[\displaystyle \frac{1}{{x}^{2}}+\displaystyle \frac{1}{{y}^{2}}\right]\delta (x)\delta (y)\end{array}\end{eqnarray} \tag{ 82 }$$

where ${f}_{\perp }\equiv {f}_{\perp }(x,y;t)$.

We note that transport Equation (79) corresponds to the one obtained in Urch (1977) if we did not have the additional function $g(x,y,t)$. However, it is questionable whether Equation (79) corresponds to the correct transport equation. Usually one assumes that compound subdiffusion incorporates a memory of past history but there is no obvious hint of past history in Equation (82). Therefore, Webb et al. (2006) discussed using fractional Fokker–Planck equations which retain memory of the diffusive flux at earlier times due to the convolution term therein. To work with fractional transport equations became more popular during recent years (see, e.g., Strauss & Effenberger 2017 and Zimbardo et al. 2017).

First we consider an approximation based on asymptotic properties. According to Abramowitz & Stegun (1974) we have

$$\begin{eqnarray}&&{\rm{erfc}}\left(\xi \right)\approx 1\quad {\rm{for}}\quad \xi \ll 1\end{eqnarray} \tag{ 83 }$$

and

$$\begin{eqnarray}&&{\rm{erfc}}\left(\xi \right)\approx \displaystyle \frac{1}{\sqrt{\pi }\xi }{e}^{-{\xi }^{2}}\quad {\rm{for}}\quad \xi \gg 1.\end{eqnarray} \tag{ 84 }$$

Thus, one can use the following approximation for the Fourier transform of the distribution function as given by Equation (35)

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\perp }({k}_{\perp };t) & \approx & \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\displaystyle \frac{1}{1+C\xi }\\ & \equiv & \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\displaystyle \frac{1}{1+C{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\sqrt{{\kappa }_{\parallel }t}}\end{array}\end{eqnarray} \tag{ 85 }$$

where we find the correct asymptotic limits for $C\,=\sqrt{\pi }\approx 1.77$. Based on approximation (85), the characteristic function becomes

$$\begin{eqnarray}&&\left\langle {e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\right\rangle \approx \displaystyle \frac{1}{1+C\xi }\equiv \displaystyle \frac{1}{1+C{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\sqrt{{\kappa }_{\parallel }t}}.\end{eqnarray} \tag{ 86 }$$

The characteristic function found here is visualized together with the exact result and another approximation in Figures 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The particle distribution function in configuration space can be derived by combining Equations (85) and (21) to obtain

$$\begin{eqnarray}&&{f}_{\perp }(\rho ;t)=\displaystyle \frac{1}{2\pi }{\int }_{0}^{\infty }{{dk}}_{\perp }\,\displaystyle \frac{{k}_{\perp }}{1+C{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\sqrt{{\kappa }_{\parallel }t}}{J}_{0}\left(\rho {k}_{\perp }\right).\end{eqnarray} \tag{ 87 }$$

This last integral can be solved via (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{{dk}}_{\perp }\,\displaystyle \frac{{k}_{\perp }}{1+\alpha {k}_{\perp }^{2}}{J}_{0}\left(\rho {k}_{\perp }\right)=\displaystyle \frac{1}{\alpha }{K}_{0}\left(\displaystyle \frac{\rho }{\sqrt{\alpha }}\right),\end{eqnarray} \tag{ 88 }$$

where we have used the modified Bessel function of the second kind ${K}_{0}(x)$. Thus, we finally find for the distribution function the following approximation

$$\begin{eqnarray}&&{f}_{\perp }(\rho ;t)\approx \displaystyle \frac{1}{2\pi C{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}{K}_{0}\left(\displaystyle \frac{\rho }{\sqrt{C{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}\right).\end{eqnarray} \tag{ 89 }$$

The latter function is visualized in Figure 3 where we can clearly see that the function is sharply peaked at the origin. We can explore the properties of this function at the origin by using (see, e.g., Abramowitz & Stegun 1974)

$$\begin{eqnarray}&&{K}_{0}\left(z\right)=-\mathrm{ln}\left(\displaystyle \frac{z}{2}\right)-\gamma \quad {\rm{for}}\quad \left|z\right|\to 0,\end{eqnarray} \tag{ 90 }$$

where we have used the Euler–Mascheroni constant γ. Clearly we find ${f}_{\perp }(\rho \to 0;t)\to \infty$. One can show by using the integral (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}&&{\int }_{0}^{\infty }d\rho \,\rho {K}_{0}\left(a\rho \right)=\displaystyle \frac{1}{{a}^{2}}\end{eqnarray} \tag{ 91 }$$

that the distribution function given by Equation (89) is correctly normalized. This can also be obtained simply by showing that approximation (85) satisfies condition (10). It is a strength of this approximation that normalization and initial conditions are satisfied exactly.

Zoom In Zoom Out Reset image size

Furthermore, we can compute the distribution function projected on the x-axis. In order to do this we combine Equations (85) and (23) to obtain

$$\begin{eqnarray}&&{f}_{\perp }(x;t)=\displaystyle \frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{{dk}}_{\perp }\,\displaystyle \frac{\cos ({k}_{\perp }x)}{1+C{\kappa }_{\mathrm{FL}}{k}_{\perp }^{2}\sqrt{{\kappa }_{\parallel }t}}.\end{eqnarray} \tag{ 92 }$$

The k_⊥-integral can be evaluated by employing Cauchy's residue theorem. After some algebra we find

$$\begin{eqnarray}&&{\int }_{-\infty }^{+\infty }{{dk}}_{\perp }\displaystyle \frac{\cos ({k}_{\perp }x)}{1+{{ak}}_{\perp }^{2}}=\displaystyle \frac{\pi }{\sqrt{a}}{e}^{-\left|x\right|/\sqrt{a}},\end{eqnarray} \tag{ 93 }$$

where we have used

$$\begin{eqnarray}&&a=C{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}.\end{eqnarray} \tag{ 94 }$$

Therewith we derive for the projected distribution function

$$\begin{eqnarray}&&{f}_{\perp }(x;t)=\displaystyle \frac{1}{2\sqrt{C{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}{e}^{-| x| /\sqrt{C{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}.\end{eqnarray} \tag{ 95 }$$

The latter function is shown in Figure 4 for different times. Again the distribution is sharply peaked but now we find a finite result at the origin, namely

$$\begin{eqnarray}&&{f}_{\perp }(x=0;t)=\displaystyle \frac{1}{2\sqrt{C{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}.\end{eqnarray} \tag{ 96 }$$

Below it will be shown that we find the exact second moment as given by Equation (61) if we set $C=2/\sqrt{\pi }$. For this value we obtain

$$\begin{eqnarray}&&{f}_{\perp }(x=0;t)\approx \displaystyle \frac{0.471}{\sqrt{{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}},\end{eqnarray} \tag{ 97 }$$

which is close to the exact result given by Equation (54).

Zoom In Zoom Out Reset image size

We can use Equation (95) to compute the nth moment via

$$\begin{eqnarray}\begin{array}{rcl}\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle & = & {\int }_{-\infty }^{+\infty }{dx}\,{x}^{n}{f}_{\perp }(x;t)\\ & = & \displaystyle \frac{1}{\sqrt{a}}{\int }_{0}^{\infty }{dx}\,{x}^{n}{e}^{-x/\sqrt{a}}\\ & = & n!{a}^{n/2}={\rm{\Gamma }}\left(n+1\right){C}^{n/2}{\left({\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}\right)}^{n/2},\end{array}\end{eqnarray} \tag{ 98 }$$

where we assumed that n is even. For odd n the moments are, of course, zero due to symmetry. For the second moment we have n = 2 and Equation (98) becomes

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{2}\rangle =2C{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}\end{eqnarray} \tag{ 99 }$$

and the running diffusion coefficient turns into

$$\begin{eqnarray}&&{d}_{\perp }(t)=\displaystyle \frac{1}{2}C{\kappa }_{\mathrm{FL}}\sqrt{\displaystyle \frac{{\kappa }_{\parallel }}{\pi t}}.\end{eqnarray} \tag{ 100 }$$

If we want to ensure that approximation (85) provides the correct second moment, we need to set $C=2/\sqrt{\pi }\approx 1.13$. The 4th moment can also easily be obtained from Equation (98) by setting n = 4 therein. We derive

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{4}\rangle =24{C}^{2}{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t.\end{eqnarray} \tag{ 101 }$$

With $C=2/\sqrt{\pi }$ this becomes

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{4}\rangle =24\displaystyle \frac{4}{\pi }{\kappa }_{\mathrm{FL}}^{2}{\kappa }_{\parallel }t.\end{eqnarray} \tag{ 102 }$$

Since $4/\pi \approx 1.2732$ this is very close to the exact result given by Equation (69). Problematic in Equation (98) is that the result sensitively depends on C. Therefore, for higher moments, the used approximation will be less accurate. Figure 5 shows a comparison of the moments up to n = 10 based on the exact formula given by Equation (70) and different approximations. Clearly we can observe a discrepancy for larger values of n.

Zoom In Zoom Out Reset image size

In the current subsection we employ an alternative approximation. Karagiannidis & Lioumpas (2007) have stated that the complementary error function can be well approximated by

$$\begin{eqnarray}&&{\rm{erfc}}\left(x\right)\approx \displaystyle \frac{1}{B\sqrt{\pi }x}\left(1-{e}^{-{Ax}}\right){e}^{-{x}^{2}},\end{eqnarray} \tag{ 103 }$$

where A = 1.98 and B = 1.135. Therewith the function given by Equation (35) can be approximated by

$$\begin{eqnarray}&&{f}_{\perp }({k}_{\perp };t)\approx \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\displaystyle \frac{1}{B\sqrt{\pi }\xi }\left(1-{e}^{-A\xi }\right)\end{eqnarray} \tag{ 104 }$$

and the characteristic function can be written as

$$\begin{eqnarray}&&\left\langle {e}^{i{{\boldsymbol{k}}}_{\perp }\cdot {{\boldsymbol{x}}}_{\perp }}\right\rangle \approx \displaystyle \frac{1}{B\sqrt{\pi }\xi }\left(1-{e}^{-A\xi }\right).\end{eqnarray} \tag{ 105 }$$

Problematic here is that the initial condition (5) is not exactly satisfied. According to Equation (36) we have $\xi =0$ for t = 0. In this case Equation (104) becomes

$$\begin{eqnarray}&&{f}_{\perp }({k}_{\perp };t=0)\approx \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\displaystyle \frac{A}{B\sqrt{\pi }}.\end{eqnarray} \tag{ 106 }$$

However, $A/(B\sqrt{\pi })\approx 0.9842$ which yields a condition very close to Equation (5). The same can be said about the normalization condition (10). Using Equation (104) in (21) yields

$$\begin{eqnarray}&&{f}_{\perp }(\rho ;t)=\displaystyle \frac{1}{2\pi B\sqrt{\pi }}{\int }_{0}^{\infty }{{dk}}_{\perp }\,{k}_{\perp }{J}_{0}(\rho {k}_{\perp })\displaystyle \frac{1-{e}^{-A\alpha {k}_{\perp }^{2}}}{\alpha {k}_{\perp }^{2}},\end{eqnarray} \tag{ 107 }$$

where we have used

$$\begin{eqnarray}&&\alpha ={\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}\end{eqnarray} \tag{ 108 }$$

and Equation (36). Employing approximation (104) in the fourth line of Equation (23), on the other hand, yields

$$\begin{eqnarray}&&{f}_{\perp }(x;t)=\displaystyle \frac{1}{B{\pi }^{3/2}}{\int }_{0}^{\infty }{{dk}}_{\perp }\,\cos \left({{xk}}_{\perp }\right)\displaystyle \frac{1-{e}^{-A\alpha {k}_{\perp }^{2}}}{\alpha {k}_{\perp }^{2}}.\end{eqnarray} \tag{ 109 }$$

Using (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}&&\begin{array}{l}{\int }_{0}^{\infty }{{dk}}_{\perp }\,\cos \left({{xk}}_{\perp }\right)\displaystyle \frac{1-{e}^{-A\alpha {k}_{\perp }^{2}}}{\alpha {k}_{\perp }^{2}}\\ \quad =\,\displaystyle \frac{2\sqrt{\alpha A\pi }{e}^{-{x}^{2}/(4\alpha A)}-\pi \left|x\right|{\rm{erfc}}\left(\left|x\right|/(2\sqrt{\alpha A}\right)}{2\alpha },\end{array}\end{eqnarray} \tag{ 110 }$$

the projected distribution function can be written as

$$\begin{eqnarray}&&{f}_{\perp }(x;t)=\displaystyle \frac{2\sqrt{\alpha A\pi }{e}^{-{x}^{2}/(4\alpha A)}-\pi \left|x\right|{\rm{erfc}}\left(\left|x\right|/(2\sqrt{\alpha A}\right)}{2B{\pi }^{3/2}\alpha }.\end{eqnarray} \tag{ 111 }$$

We can compute the value of this distribution function at the origin. We can easily derive

$$\begin{eqnarray}&&{f}_{\perp }(x=0;t)=\displaystyle \frac{\sqrt{A}}{\pi B\sqrt{\alpha }}=\displaystyle \frac{\sqrt{A}}{\pi B\sqrt{{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}}\approx \displaystyle \frac{0.3946}{\sqrt{{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}}},\end{eqnarray} \tag{ 112 }$$

which is similar but not identical compared to the exact value given by Equation (54) and the previous approximation given by Equation (97).

Calculation of the nth moment requires us to solve the integral

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle =2{\int }_{0}^{\infty }{dx}\,{x}^{n}{f}_{\perp }(x;t)\end{eqnarray} \tag{ 113 }$$

with the distribution function given by Equation (111). Therefore, in order to obtain the moments we need to evaluate the following two integrals (see, e.g., Gradshteyn & Ryzhik 2000)

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{dx}\,{x}^{n}{e}^{-{{cx}}^{2}}=\displaystyle \frac{1}{2}{c}^{-(n+1)/2}{\rm{\Gamma }}\left(\displaystyle \frac{n+1}{2}\right)\end{eqnarray} \tag{ 114 }$$

and

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{dx}\,{x}^{n+1}{\rm{erfc}}\left({cx}\right)=\displaystyle \frac{1}{(n+2)\sqrt{\pi }}{c}^{-(n+2)}{\rm{\Gamma }}\left(\displaystyle \frac{n+3}{2}\right),\end{eqnarray} \tag{ 115 }$$

where, of course, n is assumed to be an even integer number. Using these integrals in Equation (113) yields after some straightforward calculations

$$\begin{eqnarray}&&\begin{array}{l}\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle \\ \quad =\ \displaystyle \frac{1}{2\pi \alpha B}\left[{\rm{\Gamma }}\left(\displaystyle \frac{n+1}{2}\right)-\displaystyle \frac{2}{n+2}{\rm{\Gamma }}\left(\displaystyle \frac{n+3}{2}\right)\right]{\left(2\sqrt{\alpha A}\right)}^{n+2}.\end{array}\end{eqnarray} \tag{ 116 }$$

We can further simplify this by employing Equation (52). We derive

$$\begin{eqnarray}\begin{array}{rcl}\langle {\left({\rm{\Delta }}x\right)}^{n}\rangle & = & \displaystyle \frac{2}{\pi (n+2)B}{A}^{\left(n+2)/2\right.}\\ & & \times \,{\rm{\Gamma }}\left(\displaystyle \frac{n+1}{2}\right){\left(4{\kappa }_{\mathrm{FL}}\sqrt{{\kappa }_{\parallel }t}\right)}^{n/2}.\end{array}\end{eqnarray} \tag{ 117 }$$

The moments obtained here are visualized in Figure 5. We observe that the moments computed based on the approximated distribution function (111) yields a result that is too small. Whereas the lower moments still agree well with the exact values, for higher moments there is clearly disagreement.