ANALYTICAL SOLUTIONS OF A FRACTIONAL DIFFUSION-ADVECTION EQUATION FOR SOLAR COSMIC-RAY TRANSPORT

When energetic particles, accelerated in the solar corona or the solar wind, propagate in a turbulent heliospheric medium, the particle transport is often diffusive (Parker 1965). The diffusion approximation is a standard tool in the description of evolving cosmic-ray distributions (e.g., Schlickeiser & Shalchi 2008; Artmann et al. 2011, and references therein).

A generalization of the diffusion equation is obtained by replacing the usual derivatives by fractional derivatives. Formally, standard diffusion, characterized by a linear growth of the variance of a particle displacement, separates the processes of faster superdiffusion and slower subdiffusion (e.g., Saichev & Zaslavsky 1997). These processes are governed by partial differential equations with fractional operators (Samko et al. 1993). Physically, fractional differential equations conveniently describe stochastic transport when an effective mean free path is comparable with a macroscopic length scale. The corresponding diffusive process is non-local: it is described by an integral equation that can be rewritten in terms of fractional derivatives (for a clear discussion, see Chukbar 1995).

In sharp contrast to classical diffusion, solutions to fractional diffusion equations typically are not Gaussian but rather have power-law tails. The key question in concrete applications is whether the postulated anomalous diffusion is based on physical processes that can be described by a fractional evolution equation (see, e.g., Metzler & Klafter 2000; Perrone et al. 2013, for numerous potential applications).

Observed solar cosmic-ray particle distributions often appear to exhibit power-law tails, suggesting an interpretation in terms of superdiffusion. Distributions of electrons and protons, accelerated both in the solar corona and at interplanetary shocks, have recently been analyzed in terms of superdiffusive transport models (Perri & Zimbardo 2007, 2009; Sugiyama & Shiota 2011; Trotta & Zimbardo 2011; Zimbardo & Perri 2013). Notably, those studies relied on an asymptotic expression for a non-Gaussian propagator (equivalent to the Green's function of a fractional diffusion equation), which has a limited validity range.

In this paper, we develop more accurate analytical solutions and argue that they should be used to put the superdiffusive particle transport in the heliosphere on a firmer footing. As a concrete example, we examine a prototypical transport problem, described by a one-dimensional fractional diffusion-advection equation (Section 2). We derive an exact solution by the Fourier transform (Section 3), and we obtain an approximate solution in terms of elementary functions in a weak diffusion approximation (Section 4). We demonstrate the accuracy of the approximation by comparing the new solution with a semi-numerical solution based on a formally exact Fourier series expansion (Section 5), and we discuss the implication of our results on the interpretation of the observed particle distributions as evidence of superdiffusive transport in the heliosphere (Section 6).

Consider first the usual diffusion-advection equation for a one-dimensional particle distribution function f (x, t) that depends on position x and time t > 0:

$$\begin{equation} \frac{\partial f}{\partial t} = \kappa \frac{\partial ^2 f}{\partial x^2} + a \frac{\partial f}{\partial x} + \delta (x). \end{equation} \tag{ 1 }$$

Here a is a constant advection speed and κ is a constant diffusion coefficient. This transport equation follows from the Fokker–Planck equation for energetic particles if their energy losses are neglected and pitch-angle scattering is strong (e.g., Schlickeiser & Shalchi 2008; Litvinenko & Schlickeiser 2013, and references therein). A delta-functional source term on the right may correspond to energetic particles injected at a shock, and a can be interpreted as the background solar wind speed.

Standard methods (e.g., Carslaw & Jaeger 1959) give the solution of the initial value problem f(x, 0) = 0 on the interval −∞ < x < ∞:

$$\begin{equation} f(x, t) = \int _0^t (4 \pi \kappa t^{\prime })^{-1/2} \exp \left[ -\frac{(x+at^{\prime })^2}{4 \kappa t^{\prime }} \right] dt^{\prime }. \end{equation} \tag{ 2 }$$

A steady distribution is established in the limit t → ∞:

$$\begin{equation} f(x, \infty) = \frac{1}{a} \exp \left(-\frac{a}{\kappa } x \right), \quad x>0, \end{equation} \tag{ 3 }$$

$$\begin{equation} f(x, \infty) = a^{-1}, \quad x<0. \end{equation} \tag{ 4 }$$

In the remainder of the paper, we investigate the following fractional differential equation that generalizes the usual diffusion-advection equation:

$$\begin{equation} \frac{\partial f}{\partial t} = \kappa \frac{\partial ^{\alpha } f}{\partial |x|^{\alpha }} + a \frac{\partial f}{\partial x} + \delta (x) \end{equation} \tag{ 5 }$$

(e.g., Stern et al. 2014). The equation governs the evolution of a distribution function f(x, t) for t > 0 and −∞ < x < ∞. For simplicity, we assume the initial condition

$$\begin{equation} f(x, 0) = 0. \end{equation} \tag{ 6 }$$

The advection speed a and diffusion coefficient κ (now with dimensions length^α/time) are positive constants. We use the Riesz derivative to define a fractional spatial derivative:

$$\begin{eqnarray} \frac{\partial ^{\alpha } f(x,t)}{\partial |x|^{\alpha }} &=& \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (1+\alpha) \nonumber \\ &&\times \int _0^{\infty } \frac{f(x+\xi)-2f(x)+f(x-\xi)}{\xi ^{1+\alpha }} d\xi \end{eqnarray} \tag{ 7 }$$

(Samko et al. 1993; Saichev & Zaslavsky 1997). Because the derivative corresponds to a fractional Laplacian operator in higher dimensions, an alternative notation −(− Δ)^α/2f is also used (Mainardi et al. 2001). Although the regularized form above is defined for 0 < α < 2, in the following we are mainly interested in the superdiffusive case 1 < α < 2.

The Fourier transform gives a convenient method of solving the fractional diffusion-advection equation on the interval −∞ < x < ∞. Taking the Fourier transform of Equation (5) yields

$$\begin{equation} \frac{\partial \tilde{f}}{\partial t} = -\kappa |k|^{\alpha } \tilde{f} + iak \tilde{f} + \frac{1}{2 \pi }, \end{equation} \tag{ 8 }$$

where

$$\begin{equation} \tilde{f} (k, t) = \frac{1}{2 \pi } \int _{-\infty }^{\infty } f (x, t) \exp (-ikx) dx. \end{equation} \tag{ 9 }$$

Integration of the first-order equation for $\tilde{f}$ yields

$$\begin{equation} \tilde{f} (k, t) = \frac{1}{2 \pi } \frac{1 - \exp [(iak - \kappa |k|^{\alpha }) t]}{\kappa |k|^{\alpha } - iak}, \end{equation} \tag{ 10 }$$

where an integration constant is specified by the initial condition $\tilde{f} (k, 0) = 0$. Now the inverse Fourier transform gives the solution

$$\begin{equation} f(x, t) = \frac{1}{2 \pi } \int _{-\infty }^{\infty } \frac{1 - \exp [(iak - \kappa |k|^{\alpha }) t]}{\kappa |k|^{\alpha } - iak} \,\exp\, (ikx)\, dk. \end{equation} \tag{ 11 }$$

In the non-diffusive case, κ = 0 and the integral is evaluated to give an expanding "top-hat" solution

$$\begin{equation} \left. f (x, t) \right|_{\kappa =0} = f_0 (x, t) = \frac{1}{2a} \left[ {\rm sgn}(x+at) - {\rm sgn}(x) \right]. \end{equation} \tag{ 12 }$$

The solution of Equation (5) can also be expressed as

$$\begin{equation} f(x, t) = \int _0^t G(x+at^{\prime }, t^{\prime }) dt^{\prime }, \end{equation} \tag{ 13 }$$

where the Green's function G(x, t) satisfies the fractional diffusion equation

$$\begin{equation} \frac{\partial G}{\partial t} = \kappa \frac{\partial ^{\alpha } G}{\partial |x|^{\alpha }} + \delta (x) \delta (t). \end{equation} \tag{ 14 }$$

Its Fourier transform is given by

$$\begin{equation} \frac{\partial \tilde{G}}{\partial t} = -\kappa |k|^{\alpha } \tilde{G} + \frac{1}{2 \pi } \delta (t). \end{equation} \tag{ 15 }$$

It follows that

$$\begin{equation} \tilde{G} (k, t) = \frac{1}{2 \pi } \exp (-\kappa |k|^{\alpha } t), \end{equation} \tag{ 16 }$$

and the inverse Fourier transform yields

$$\begin{equation} G(x, t) = \frac{1}{\pi } \int _{0}^{\infty } \exp (-\kappa k^{\alpha } t) \cos (kx) dk \end{equation} \tag{ 17 }$$

(e.g., Chukbar 1995).

An asymptotic expression for G(x, t) is obtained by integrating Equation (17) by parts and applying an analogue of Watson's lemma for Fourier type integrals (e.g., Ablowitz & Fokas 1997). For x ≫ (κt)^1/α, the result is

$$\begin{equation} G(x, t) \approx \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (1+\alpha) \frac{\kappa t}{|x|^{1+\alpha }}. \end{equation} \tag{ 18 }$$

For x > 0, substitution into Equation (13) yields

$$\begin{equation} f(x, t) \approx \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (1+\alpha) \int _0^t \frac{\kappa t^{\prime }}{(x+at^{\prime })^{1+\alpha }} dt^{\prime }, \end{equation} \tag{ 19 }$$

and so the distribution function f(x, t) is given by

$$\begin{equation} f(x, t) \approx \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (\alpha -1) \frac{\kappa }{a^2} \left[ x^{1-\alpha } - \frac{x+\alpha a t}{(x+at)^{\alpha }} \right] \end{equation} \tag{ 20 }$$

as long as x + at ≫ (κt)^1/α and 1 < α < 2. Two limiting cases are as follows:

$$\begin{eqnarray} f(x, t) &\approx \frac{1}{2 \pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (1+\alpha) \frac{\kappa t^2}{x^{1+\alpha }}, \quad x \gg at, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} f(x, t) \approx \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (\alpha -1) \frac{\kappa }{a^2} x^{1-\alpha },\quad 0 &<& x \ll at.\nonumber\\ \end{eqnarray} \tag{ 22 }$$

Equation (22) is essentially the asymptotic power law that was previously used in the data analysis of energetic particles, accelerated in the solar corona (Trotta & Zimbardo 2011) and at interplanetary shocks (Perri & Zimbardo 2007, 2009; Sugiyama & Shiota 2011). For instance, Perri & Zimbardo (2007, 2009) used a formula for the particle density at a fixed point due to a shock-associated source moving with a speed V_sh. The formula follows from our analysis by changing the reference frame. Suppose a shock, initially located at x₀ = −V_sht₀, moves at speed V_sh and arrives at x = 0 when t = 0. On changing variables, Equation (20) yields

$$\begin{eqnarray} f(x, t) &\approx& \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (\alpha -1) \frac{\kappa }{V_{\rm sh}^2} \\ &&\times \left[ (x - V_{\rm sh} t)^{1-\alpha } - \frac{x + (\alpha - 1) V_{\rm sh} t + \alpha V_{\rm sh} t_0}{(x+ V_{\rm sh} t_0)^{\alpha }} \right] \nonumber \,, \end{eqnarray} \tag{ 23 }$$

which reduces to Equation (4) in Perri & Zimbardo (2007) on assuming that the shock is coming from a very large distance, t₀ → ∞. Note for clarity that Perri & Zimbardo (2007) do not specify a normalization constant b in their Equation (1), and that their notation is different: their μ is our 1 + α, and their α is our 3 − α. The latter expression appears in the dependence of the variance of a particle displacement on time ∼t^{3 − α} when a finite particle speed is taken into account.

A limitation of the analysis in the previous section is that we used a short-time asymptotic (18) for the Green's function G(x,t) to derive Equation (20) for the particle distribution function f(x, t), and so it is not clear whether the results are valid for at ≫ x. In addition, the results are only valid for x > 0 because the integral in Equation (13) diverges for x < 0 if Equation (18) is used to evaluate the integral. To remove these limitations of the analysis, we solve for f(x,t) in a weak diffusion approximation that allows us to evaluate the integral in Equation (11) for both x > 0 and x < 0.

The superdiffusive term in Equation (5) can be treated as a perturbation sufficiently far from the locations where the advective solution (12) predicts jumps in f (x, t). Suppose l_d is the distance from a jump where the diffusive and advective terms in Equation (5) are comparable. To an order of magnitude, $\kappa \partial ^{\alpha } f / \partial |x|^{\alpha } \sim \kappa f / l_d^{\alpha }$ and a∂f/∂x ∼ af/l_d. The diffusion length is thus defined as

$$\begin{equation} l_d =(\kappa / a)^{1/(\alpha -1)}. \end{equation} \tag{ 24 }$$

Now assuming that

$$\begin{equation} l_d \ll |x|, |x+at|, at, \end{equation} \tag{ 25 }$$

we can formally treat κ as a small parameter, and so

$$\begin{equation} f(x, t) \approx f_0 (x, t) + \kappa \left. \frac{df(x, t)}{d\kappa } \right|_{\kappa =0}, \end{equation} \tag{ 26 }$$

where f₀(x, t) is given by Equation (12). Differentiation of Equation (11) with respect to κ yields

$$\begin{eqnarray} f(x, t) &\approx& f_0 (x, t) \nonumber \\ &&+\frac{\kappa }{\pi a^2} \int _{0}^{\infty } \left [ \frac{\cos kx - \cos k(x+at)}{k^{2-\alpha }}\right. \nonumber \\ &&\left. - \frac{(\alpha -1)at}{(x+at)} \frac{\cos k(x+at)}{k^{2-\alpha }} \right ] dk, \end{eqnarray} \tag{ 27 }$$

where the last term in the integrand is obtained by integrating by parts and neglecting a rapidly varying term containing exp [ik(x + at)] at the upper integration limit. On simplifying and using Watson's lemma, we get

$$\begin{eqnarray} f(x, t) &\approx& f_0 (x, t) + \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (\alpha -1) \frac{\kappa }{a^2}\nonumber \\ &&\times \left[ |x|^{1-\alpha } - \frac{x+\alpha a t}{(x+at) |x+at|^{\alpha -1}} \right]. \end{eqnarray} \tag{ 28 }$$

For x > 0, we recover Equation (20) and its limiting cases, confirming the analysis of the previous section. For x < 0, we have

$$\begin{equation} f(x, t) \approx \frac{1}{2 \pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (1+\alpha) \frac{\kappa t^2}{|x|^{1+\alpha }}, \quad |x| \gg at, \end{equation} \tag{ 29 }$$

$$\begin{equation} f(x, t) \approx \frac{1}{a} + \frac{1}{\pi } \sin \left(\frac{\pi }{2} \alpha \right) \Gamma (\alpha -1) \frac{\kappa }{a^2} |x|^{1-\alpha }, \quad |x| \ll at. \end{equation} \tag{ 30 }$$

The discontinuity at x = −at broadens into a smoother transition. The solution is inapplicable in the vicinity of the jump at |x + at| ≈ 0, because the weak diffusion approximation is valid only as long as |x|, |x + at|, and at are large in comparison with the diffusion length l_d.

Stern et al. (2014) gave the following exact Fourier series solution to Equation (5) on a domain of length L:

$$\begin{eqnarray} f(x,t) &=& \sum _{n=1}^{\infty }\left \lbrace (1+(-1)^{n+1}) \left [ \left(\frac{n \pi }{L}\right)^{\alpha } \kappa L \cos \left(\frac{n \pi x}{L} \right)\right.\right.\nonumber \\ &&\left.\left. -n \pi a \sin \left(\frac{n \pi }{L}x \right)\right.\right. \nonumber \\ &&\left.\left. -\left(\frac{n \pi }{L} \right)^{\alpha } \kappa L \cos \left(\frac{n \pi }{L}(x+at) \right) \exp \left(-\left(\frac{n \pi }{L} \right)^{\alpha } \kappa t \right)\right.\right. \nonumber \\ &&\left.\left. + n \pi a \sin \left(\frac{n \pi }{L} (x+at) \right) \exp \left(-\left(\frac{n \pi }{L} \right)^{\alpha } \kappa t \right) \right ]{{\left/\right. }}\right. \nonumber \\ &&\left. \left(\left(\frac{n \pi }{L}\right)^{2 \alpha }\kappa ^2L^2 + n^2 \pi ^2a^2 \right) \right \rbrace \,. \end{eqnarray} \tag{ 31 }$$

Because f(x, t) = f(x + 2L, t), the series does not represent the solution of an initial value problem on an infinite interval for t → ∞. For a localized source and finite t, however, the series solution accurately represents f (x, t) on an infinite interval if L is sufficiently large. In practice, we achieve accuracy by choosing L ≫ at.

We compare the new analytical solutions of the previous sections with a semi-numerical solution based on the Fourier series expansion. We use Equation (31) and sum up N = 10⁶ terms with L = 1000 to achieve high accuracy. We set a = 1, which simply means that in the following we measure speeds in units of the solar wind speed. We choose α = 1.5 in agreement with the range of values inferred from the heliospheric particle data (Perri & Zimbardo 2007, 2009; Sugiyama & Shiota 2011). The superdiffusion coefficient κ is a key parameter of the theory. Work is underway to estimate κ from the data (S. Perri et al., in preparation). We adopt κ = 0.5 as an illustration and investigate how the particle distribution evolves over a few hundred advection times.

Figure 1 shows the results in a semi-logarithmic plot at time t = 10. We find a good agreement between the analytical (solid black line) and semi-numerical (black symbols) solutions, except for the region around x = −at where the weak diffusion approximation is not valid. The red box in the figure illustrates the non-diffusive solution given by Equation (12). We have truncated the analytical solution in the vicinity of x = −at over a length l of 10 times the diffusion length (l = 10 l_d = 2.5). The green and blue lines give the steady-state solution for the Gaussian diffusion case, given by Equation (3), and the approximate steady-state solution, given by Equation (30).

Zoom In Zoom Out Reset image size

Figure 2 gives the same solutions as in Figure 1 at a later time t = 200. The weak diffusion solution and the Fourier series remain in good agreement as they slowly approach the steady state. The downstream region in Figure 2 is already completely filled since |x| ≪ at. An interesting feature of the solution is a peak at the injection site x = 0, which is not present for Gaussian diffusion.

Zoom In Zoom Out Reset image size

We used the Fourier transform to analytically solve a fractional diffusion-advection equation for cosmic-ray transport, and we applied the solution to the problem of describing the transport of energetic particles, accelerated at a traveling heliospheric shock. We also developed a weak diffusion approximation based on the exact Fourier transform solution. We confirmed the validity of the approximation for both early and late times by comparing it with an exact Fourier series solution. Our analysis is motivated by recent applications of superdiffusive transport models to the observed shock-accelerated particle distributions (Perri & Zimbardo 2007, 2009; Sugiyama & Shiota 2011).

Our new solution quantifies the limited validity of the asymptotic expressions, used previously to interpret the particle data. Specifically, the formula used by Perri & Zimbardo (2007, 2009) and Sugiyama & Shiota (2011) is basically our Equation (23) in the limit t₀ → ∞, corresponding to a shock approaching an observer from a very large distance V_sht₀. As our results show, however, it may take a very long time for the asymptotic expression to become accurate. The ratio of the second term in Equation (23) to the first one is (− t/t₀)^{α − 1}(α + (α − 1)t/t₀) at x = 0, and so our more accurate solution differs from the t₀ = ∞ asymptotic expression by about a factor of two when the distance between the observer and the shock is as short as one tenth of the initial distance V_sht₀ between them.

To sum up, solar cosmic-ray data in various settings appear to be consistent with asymptotic propagator solutions to a fractional diffusion equation or more general continuous-time random-walk models (Zimbardo & Perri 2013). We argued, however, that more accurate solutions of an appropriate transport equation should be used for validating the superdiffusive transport of energetic particles in the heliosphere. In the context of the transport of particles, accelerated at a traveling heliospheric shock, our analysis strongly suggests that we should not assume the initial distance V_sht₀ of the shock from the observer to be infinite. The shock travel time t₀ should be a parameter of the superdiffusive transport model.

Y.L. thanks Horst Fichtner for drawing his attention to the topic of this paper. F.E. appreciates support from the International Space Science Institute (ISSI) during a team meeting in Bern on superdiffusive transport and thanks the team members for many insightful discussions. The authors acknowledge the anonymous referee for several useful suggestions.