Energetic Particle Superdiffusion in Solar System Plasmas: Which Fractional Transport Equation?

Abstract

Superdiffusive transport of energetic particles in the solar system and in other plasma environments is often inferred; while this can be described in terms of Lévy walks, a corresponding transport differential equation still calls for investigation. Here, we propose that superdiffusive transport can be described by means of a transport equation for pitch-angle scattering where the time derivative is fractional rather than integer. We show that this simply leads to superdiffusion in the direction parallel to the magnetic field, and we discuss some advantages with respect to approaches based on transport equations with symmetric spatial fractional derivates.

1. Introduction

The random motion of particles in homogeneous media can often be described by diffusive propagation, such that the mean square displacement of the particle position grows linearly in time, $\langle \Delta x^2\rangle \propto t$, where $\Delta x=x\left(t\right)-x_0$, and where, for the moment, we consider one dimension only. On the other hand, many physical systems exhibit anomalous transport, meaning that the mean square displacement grows nonlinearly in time, i.e., $\langle \Delta x^2\rangle \propto t^{\gamma }$, with $\gamma <1$ for subdiffusion and $\gamma >1$ for superdiffusion [1,2,3,4]. In astrophysics and space physics there are several indications of superdiffusive transport of energetic particles, which come from the observations of nonrelativistic solar electrons [5], from the analysis of energetic particle profiles observed by spacecraft at shock crossings [6,7,8,9,10], from the analysis of electron transport at the Coma cluster of galaxies [11], from the analysis of the extended precursor of supernova remnants [12,13], and from the radio-derived energy spectra of galaxy cluster merger shocks [14,15]. Furthermore, we notice that many observations of solar energetic particles are consistent with very weak pitch-angle scattering and long (even larger than one astronomical unit) mean free paths, so that “scatter free” propagation is also stemming out from measurements. Typically, pitch-angle scattering is due to the interaction of charged particles with low-frequency magnetohydrodynamic turbulence in the case of ions [16,17], and with whistler mode waves in the case of electrons [18,19], and the variability of the statistical properties of such fluctuations (i.e., intensity, spectral extension, spectral anisotropy) can lead to several different transport regimes.

Superdiffusive propagation can be appropriately described by the continuous time random walk model [1,4,20,21]; the concept of Lévy walks, which involves a coupling between the free path lengths and the time needed to cover the free paths (resulting in a time cost for long displacements), turns out to be pivotal to avoid some divergences associated with Lévy flights (see below). Based on the Lévy walk model, Perri and Zimbardo [22,23] have extended the theory of diffusive shock acceleration (DSA) to the case of superdiffusive shock acceleration (SSA) in a consistent way, and SSA has been applied to the interpretation of energetic particle acceleration in the heliosphere and at galaxy merger shocks [8,15,24].

On the other hand, a complete description of superdiffusion also requires the availability of a transport equation, such as, e.g., a Fokker–Planck equation. Indeed, the normal diffusion equation, here in one dimension,

$$\frac{\partial f}{\partial t}=\frac{\partial }{\partial x}{D}_{xx}\frac{\partial f}{\partial x}$$

(1)

readily implies normal diffusion, with $\langle \Delta x^2\rangle =2D_{xx}t$. In order to obtain a transport equation describing anomalous diffusion, fractional time and/or spatial derivatives are to be used [1,2,3,25,26,27,28], and the transport equation is sometimes termed fractional Fokker–Planck equation. Fractional derivatives are integrodifferential operators which take into account the non-Markovian, nonlocal character of anomalous transport. Fractional transport equations have been introduced by several authors, with fractional time derivatives usually describing subdiffusion [1,3,25], and fractional spatial derivatives usually describing superdiffusion [2,3,25]. The relation with persistent and antipersistent processes, as described by the Hurst exponent, has been considered, too [29].

In this paper, we address the problem, never discussed so far, of whether particle superdiffusion (i.e., the superdiffusive spreading in coordinate space) should be described by fractional derivatives in space (as commonly assumed) or fractional derivatives in time. We discuss the problem of the diverging moments of the solutions that stem from the fractional derivative in space, i.e., those corresponding to Lévy flights. We show how this problem is not encountered if a time fractional equation in velocity space is introduced and the Lévy walk properties are implemented by considering pitch-angle scattering only. Thus, the latter approach allows us to describe superdiffusion in coordinate space. This may be surprising, as time fractional derivatives are usually associated with subdiffusion.

In addition, this work tries to build a bridge between the solar and heliospheric communities, interested in energetic particle pitch-angle scattering, and the mathematical physics community, interested in fractional differential equations. We believe that this approach sheds some light on the microscopic dynamics of superdiffusion in space and astrophysical plasmas.

2. Spatial Fractional Derivatives and Lévy Flights

The diffusion Equation (1) consistently leads to normal diffusion, with a diffusion coefficient approximately given by

$$D_{xx}\simeq \frac{\langle {\ell }^2\rangle }{\langle \tau \rangle }$$

(2)

where $\langle {\ell }^2\rangle$ is the mean square free path length, and $\langle \tau \rangle$ is the mean scattering time (or travel time). On the other hand, anomalous diffusion is due to the fact that either $\langle {\ell }^2\rangle$ and/or $\langle \tau \rangle$ are not finite, so that the central limit theorem does not hold. In particular, this means that both ℓ and/or $\tau$ do not have a typical value. This happens when the probability distribution of the values of ℓ and $\tau$ has power-law tails with small exponents such that $\langle {\ell }^2\rangle$ and/or $\langle \tau \rangle$ diverge, so that long displacements and long travel times are possible. These possibilities are taken into account by fractional derivatives [1,3,25,26,27,30]; in the time domain, they can be defined as

$$\frac{{\partial }^{\beta }f\left(t\right)}{\partial t^{\beta }}\equiv {}_0^C{D}_t^{\beta }f\left(t\right)=\frac{1}{\Gamma (1-\beta )}{\int }_0^t\frac{\partial f\left(t^{\prime }\right)}{\partial t^{\prime }}\frac{d{t}^{\prime }}{{(t-t^{\prime })}^{\beta }}.$$

(3)

with $0<\beta <1$ for the fractional time derivative in the Caputo form [25,31], which exploits the non-Markovian property by considering $\partial f/\partial t^{\prime }$ also in times earlier than the current time t. Here, $\Gamma$ is Euler gamma function. In the spatial domain, fractional derivatives can be defined by

$$\begin{array}{c}\hfill \frac{{\partial }^{\alpha }f\left(x\right)}{{\partial \left|x\right|}^{\alpha }}=\frac{1}{\pi }sin\left(\frac{\pi }{2}\alpha \right)\Gamma (1+\alpha ){\int }_0^{\infty }\frac{f(x+\xi )-2f\left(x\right)+f(x-\xi )}{{\xi }^{1+\alpha }}d\xi ,\end{array}$$

(4)

with $0<\alpha <2$ for the symmetric fractional derivative in Riesz’s form (see also [30]), which clearly shows the nonlocal character of transport, since the distribution function at positions $x+\xi$ and $x-\xi$ is involved in the integral. It should be noted that several other forms of fractional derivatives can be defined [32], and that the Caputo form is a favored choice for time derivatives since it allows us to take into account the initial conditions [25,33]. Then, a fully fractional transport equation can be written as [3,25,28]

$${}_0^C{D}_t^{\beta }f(x,t)={\mathcal{D}}_{\alpha ,\beta }\frac{{\partial }^{\alpha }f(x,t)}{{\partial \left|x\right|}^{\alpha }}.$$

(5)

Dimensional analysis shows that the mean square displacement should scale as

$$\langle \Delta x^2\rangle \propto t^{2\beta /\alpha }$$

(6)

so that by varying $\alpha$ and $\beta$ we can obtain a variety of anomalous transport regimes [1,2,3,25]. In particular, if the spatial derivatives are integer, i.e., if $\alpha =2$, subdiffusion with

$$\langle \Delta x^2\rangle \propto t^{\beta }$$

(7)

and $\beta <1$ is found [1,3]. Therefore, subdiffusion is readily described by fractional time derivatives. On the other hand, the scenario is somewhat more involved for superdiffusion: for integer time derivative, $\beta =1$, and fractional spatial derivative, we obtain

$$\frac{\partial f(x,t)}{\partial t}={\mathcal{D}}_{\alpha }\frac{{\partial }^{\alpha }f(x,t)}{{\partial \left|x\right|}^{\alpha }}.$$

(8)

Let us search for solutions by Fourier transforming the x coordinate. Then, in Fourier space we obtain

$$\frac{\partial \widehat{f}(k,t)}{\partial t}=-{\mathcal{D}}_{\alpha }{\left|k\right|}^{\alpha }\widehat{f}(k,t),$$

(9)

in agreement with the general rule for the Fourier transform of Riesz derivatives [25,27,29]. We can solve with respect to time to find

$$\widehat{f}(k,t)=exp\left[-{\mathcal{D}}_{\alpha }{\left|k\right|}^{\alpha }t\right],$$

(10)

which, for $\alpha <2$, is the Fourier transform of a symmetric Lévy distribution. Backtransforming in the limit $\left|x\right|\gg {\left({\mathcal{D}}_{\alpha }t\right)}^{1/\alpha }$ yields [3,9,23].

$$f(x,t)\approx \frac{\Gamma (\alpha +1)}{\pi }sin\left(\frac{\pi }{2}\alpha \right)\frac{{\mathcal{D}}_{\alpha }t}{{\left|x\right|}^{\alpha +1}},$$

(11)

So we can see that the probability density $f(x,t)$, also called the propagator, has power-law tails, $f(x,t)\sim t/{\left|x\right|}^{\alpha +1}$, typical of Lévy distributions. In particular, a Lévy distribution has tails much longer than a Gaussian distribution. However, if we would like to compute the mean square displacement of particles starting at $x_0=0$, so that $\Delta x=x\left(t\right)$,

$$\langle \Delta x^2\rangle \left(t\right)={\int }_{-\infty }^{\infty }{x}^{2}f(x,t)dx\propto {\int }_{-\infty }^{\infty }\frac{x^2}{{\left|x\right|}^{\alpha +1}}dx,$$

(12)

we would find, for $\alpha \le 2$, that $\langle \Delta x^2\rangle$ diverges for any time t. This is the well-known property of Lévy distribution, i.e., that moments of order $\delta \ge \alpha$ are diverging [1,20,21], so that superdiffusion, i.e., $\langle \Delta x^2\rangle \propto t^{2/\alpha }$, is not recovered. Processes described by Equations (8)–(11) are called Lévy flights. They are appropriate to describe a number of systems [2,3], but cannot describe the propagation of particles with finite mass, since Lévy flights may imply an infinite velocity: ”For such massive particles, a finite velocity of propagation exists, making long instantaneus jumps impossible” [1]. Therefore, in spite of the formal appeal of Equation (8), this spatial fractional equation has some limits in describing superdiffusion, and one has to use a coupled scheme called Lévy walks [4,20,34].

3. Time Fractional Derivatives and Lévy Walks

Superdiffusion of material particles, i.e., of particles having mass, can be described by a statistical process called Lévy walk [4,20,34]. In this process, the random walk of a particle is described by a symmetric probability density $\Psi (\ell ,\tau )$ of making a free path of length ℓ (positive or negative) in a time $\tau$ given by $\Psi (\ell ,\tau )=\frac{1}{2}\psi \left(\tau \right)\delta \left(\right|\ell |-v\tau )$, where $\psi \left(\tau \right)$ is a power-law tailed distribution of times, $\psi \left(\tau \right)\sim 1/{\tau }^{\alpha +1}$. The long tails of $\psi \left(\tau \right)$ make for the nonlocal, non-Markovian character of Lévy walks. Here, the essential point is that the free path length, $|\ell |$, is related to the travel time, $\tau$, by the delta coupling, so that $|\ell |=v\tau$. This means that longer free paths require longer times; that is, particles move at constant velocity (further models where $|\ell |=v{\tau }_0{(\tau /{\tau }_0)}^{\nu }$ can be introduced [21], but here, we restrict ourselves to the more common case).

Motion at constant speed avoids the divergence of $\langle \Delta x^2\rangle$, and indeed superdiffusion with $\langle \Delta x^2\rangle \propto t^{3-\alpha }$ is found for Lévy walks [4,20,21,23,34]. Nevertheless, it is not immediate to write a fractional transport equation for Lévy walks in coordinate space (see, however, [2,35]).

In the current paper, we propose an alternative scheme, where a transport equation in velocity space is introduced, with the time derivative being fractional. We start from the cosmic ray transport equation, which can be written as a Fokker–Planck equation in phase space [27,36,37,38,39,40,41,42]; averaging over position, we have (e.g., [43], Equation (2))

$$\frac{\partial f}{\partial t}=\frac{\partial }{\partial p_i}{D}_{p_i{p}_i}\frac{\partial f}{\partial p_i}$$

(13)

We introduce spherical coordinates in momentum space $\mathbf{p}=(p,\theta ,\varphi )$. In magnetized plasmas such as those of the solar system, the polar angle $\theta$ represents the angle between the magnetic field and the momentum vector, that is, the pitch angle. Following common usage, we also introduce the pitch-angle cosine $\mu =cos\theta$, obtaining

$$\frac{\partial f}{\partial t}=\frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2{D}_{pp}\frac{\partial f}{\partial p}\right)+\frac{\partial }{\partial \mu }{D}_{\mu \mu }\frac{\partial f}{\partial \mu }+\frac{\partial }{\partial \varphi }{D}_{\varphi \varphi }\frac{\partial f}{\partial \varphi }.$$

(14)

Further, we assume that $f(p,\mu ,\varphi )$ is gyrotropic, so that $\partial f/\partial \varphi =0$. Then, we can write

$$\frac{\partial f}{\partial t}=\frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2{D}_{pp}\frac{\partial f}{\partial p}\right)+\frac{\partial }{\partial \mu }{D}_{\mu \mu }\frac{\partial f}{\partial \mu }.$$

(15)

This equation describes the diffusive evolution of f in momentum (and, as a consequence, velocity) space. The first term on the right-hand side (r.h.s.) implies diffusion in energy, and corresponds to second-order Fermi acceleration, while the second term on the r.h.s. represents pitch-angle scattering. We now impose the condition characterizing Lévy walks versus Lévy flights, i.e., that $v=\left|\mathbf{v}\right|=$ const for each particle displacement. Of course, this also means $\left|\mathbf{p}\right|=$ const (except for ultrarelativistic particles, which, however, are less relevant to solar system plasmas). Then, Fermi acceleration is not to be expected and the first term on the r.h.s., i.e., $D_{pp}$, is to be set to zero. We may notice that if Fermi acceleration would be active, with the increase of speed in time it would be possible to have transport regimes that are not only superdiffusive but even superballistic [44,45].

Therefore, with the constraint $v=$ const on particle motion, we are left with

$$\frac{\partial f}{\partial t}=\frac{\partial }{\partial \mu }{D}_{\mu \mu }\frac{\partial f}{\partial \mu },$$

(16)

corresponding to the well-known pitch-angle diffusion equation [38,46,47]. We emphasize that this equation is consistent with motion at constant speed, so that the particle velocity diffuses on a constant velocity shell, i.e., a constant energy shell. This embodies a particle random walk at constant speed where only the velocity direction is changing—this is just one of the properties of Lévy walks, but long-range correlations are also needed.

With proper boundary conditions, Equation (16) describes normal “diffusion” in $\mu$ space. To obtain anomalous pitch-angle diffusion, we make the fundamental assumption of the present paper, i.e., we replace the integer time derivative by a fractional time derivative [31]:

$$\frac{{\partial }^{\beta }f}{\partial t^{\beta }}=\frac{\partial }{\partial \mu }{\mathcal{D}}_{\mu \mu }\frac{\partial f}{\partial \mu }.$$

(17)

This choice is motivated by the fact that the particle pitch-angle scattering times, in the presence of magnetic turbulence or magnetic rotational discontinuities, are found to be characterized by power-law probability distributions [48,49,50,51,52]. Those probabilities exhibit a power-law behaviour over more than three decades both in the experimental data measured in the solar wind [48,50] and in the numerical simulations of particle transport either in a three-dimensional spectrum of magnetic fluctuations [49] or in the numerical study of the effects of magnetic rotational discontinuities [51]. Such power-law distributions of scattering times correspond to a fractional time derivative in the Caputo form [1,25,27,31], and the presence of very long scattering times substantiates the non-Markovian character of transport.

For stationary, homogeneous turbulence, i.e., ${\mathcal{D}}_{\mu \mu }$ independent of time, Equation (17) can easily be solved by separation of variables. We can write the solutions in the form

$$f(\mu ,t)=g\left(t\right)h\left(\mu \right)$$

(18)

and we find for the time-depending function $g\left(t\right)$

$${}_0^C{D}_t^{\beta }g\left(t\right)=-Ag\left(t\right)$$

(19)

where A is the separation constant. The solution of this time-fractional, linear, differential equation is given by the Mittag-Leffler functions $E_{\beta }$ of index $\beta$:

$$g\left(t\right)=g\left(0\right)E_{\beta }[-A{t}^{\beta }]=g\left(0\right)\sum _{n=0}^{\infty }\frac{{\left(-A{t}^{\beta }\right)}^n}{\Gamma (\beta n+1)}$$

(20)

see [1,31,53]. The Mittag-Leffler functions reduce to an exponential for $\beta \to 1$, and have a power-law decay for large t otherwise (see Figure 1). The diffusion coefficient ${\mathcal{D}}_{\mu \mu }$ can depend on $\mu$ in several ways [47,54]. For sufficiently strong turbulence, an isotropic scattering process in pitch angle can be assumed [54], i.e., ${\mathcal{D}}_{\mu \mu }=(1-{\mu }^2){\mathcal{D}}_0$, and the full solution can be expressed in terms of Legendre polynomials $P_l\left(\mu \right)$ as [1,31,47]

$$f(\mu ,t)=\frac{1}{2}+\sum _{l=1}^{\infty }\frac{2l+1}{2}{P}_l\left({\mu }_0\right)P_l\left(\mu \right)E_{\beta }[-{\mathcal{D}}_{0}l(l+1)t^{\beta }],$$

(21)

where ${\mu }_0$ represents the initial pitch angle. For long times, the leading terms are

$$f(\mu ,t)=\frac{1}{2}+\frac{3}{2}{P}_1\left({\mu }_0\right)P_1\left(\mu \right)E_{\beta }[-2{\mathcal{D}}_0{t}^{\beta }]$$

(22)

which decay as $1/2{\mathcal{D}}_0{t}^{\beta }$, i.e., much more slowly than in the case of an exponential decay, towards an isotropic distribution.

We now consider the motion along the magnetic field, that is, along the z direction, considering that $v_z=v\mu$. Transport along this direction can be expressed in terms of the velocity autocorrelation function $\langle v_z\left(t\right)v_z\left(0\right)\rangle =v^2\langle \mu \left(t\right){\mu }_0\rangle$. The average can be obtained from the probability density $f(\mu ,t)$, and, remembering that $P_1\left(\mu \right)=\mu$ and the orthogonality condition for Legendre polynomials, one readily finds [31,47]

$$\langle v_z\left(t\right)v_z\left(0\right)\rangle =v^2{\mu }_0^2{E}_{\beta }[-2{\mathcal{D}}_0{t}^{\beta }].$$

(23)

Equation (23) corresponds to a basic result, i.e., that the velocity autocorrelation function does not decay exponentially, but, for long times, it decays as a power law in time,

$$\langle v_z\left(t\right)v_z\left(0\right)\rangle \simeq \frac{v^2{\mu }_0^2}{\Gamma (1-\beta )}\frac{1}{2{\mathcal{D}}_0{t}^{\beta }}$$

(24)

with $0\le \beta \le 1$. The latter result can be derived by considering the Laplace transform of Mittag-Leffler functions of index $\beta$, i.e., $s^{\beta -1}/(s^{\beta }+2{\mathcal{D}}_0)$, then taking the limit for small s, which corresponds to long times t, and finally taking the Laplace back-transform [33,53]. Equation (24) implies that the memory of the initial velocity, $v{\mu }_0$, is retained for long times. Then, the mean square displacement along the z direction can be obtained by means of the generalized Taylor–Green–Kubo formula [55]

$$\langle \Delta z^2\rangle =2t{\int }_0^{t}d{t}^{\prime }\langle v_z\left(t^{\prime }\right)v_z\left(0\right)\rangle =2v^2{\mu }_0^{2}t{\int }_0^{t}d{t}^{\prime }{E}_{\beta }[-2{\mathcal{D}}_0{t}^{\prime \beta }]$$

(25)

whose long time limit is

$$\langle \Delta z^2\rangle \simeq \frac{v^2{\mu }_0^2}{\Gamma (1-\beta )}\frac{1}{{\mathcal{D}}_0}\frac{t^{2-\beta }}{1-\beta }.$$

(26)

Finally, for $\beta <1$ we have superdiffusion along z with anomalous diffusion exponent given by

$$\gamma =2-\beta .$$

(27)

Thus, we can see that the assumption of a fractional time derivative for the pitch-angle diffusion equation leads rather directly to superdiffusion along the magnetic field. We notice that the above results can also be obtained by means of a fractional Langevin equation with a power-law tailed memory kernel [56].

4. Discussion and Conclusions

In this paper, we have discussed the topic of superdiffusive propagation of energetic particles, which is frequently observed in solar system plasmas and also inferred for remote astrophysical systems such as supernova remnant shocks and galaxy cluster shocks. Anomalous transport can be described by fractional differential equations, and we derive the well-known result that a diffusion equation with fractional spatial derivatives has, as solutions, the Lévy distribution. These functions have long power-law tails, but the lack of a coupling between free path length and travel time leads to the divergence of the mean square displacement at any given time. Therefore, a different approach, based on Lévy walks, is needed. Here, we propose to use a transport equation in velocity space, as commonly used for the description of cosmic ray transport. The condition of constant velocity, necessary for Lévy walks, is implemented by setting to zero the diffusion coefficient $D_{pp}$ in momentum modulus p. This yields the well-known pitch-angle diffusion equation. The anomalous character of transport is implemented by changing the normal time derivative to a fractional time derivative of order $\beta <1$, which leads to subdiffusion in pitch-angle, but superdiffusion for motion along the z direction, i.e., along the magnetic field. Indeed, the time fractional differential equation can be solved by separation of variables, and the solution of the temporal part is given in terms of the Mittag-Leffler functions. Therefore, an exact solution is obtained both for short and for long times; then, for long times, a power-law decay of the velocity autocorrelation function is obtained. When this is inserted in the Taylor–Green–Kubo formula, superdiffusion along z is readily obtained.

Another approach to superdiffusive transport equations is based on the fact that the delta coupling $\delta \left(\right|\ell |-v\tau )$ of Lévy walks implies a definite relation between the probability distributions of free path lengths and of travel times, that is, they should have the same exponent, and this leads to fractional derivatives of the same order in space and time. Indeed, a fractional material derivative was proposed by Sokolov and Metzler [35], where a single fractional operator encompasses both spatial and time derivatives. Use of the Fourier–Laplace transforms allows us to recover the Lévy walk propagator and the ensuing properties [2,35]. Still, we believe that the approach presented here for studying superdiffusion clearly identifies the assumptions needed to describe Lévy walks by means of fractional differential equations, and is somewhat simpler than the approach based on the fractional material derivative.

It is important to point out that other approaches to fractional transport equations have been used for solar energetic particles and cosmic rays. For instance, Refs. [26,27] use a Fokker–Planck equation in energy space where the time derivative is integer, while the energy equation is fractional, and this allows us to obtain a power-law energy spectrum. This is not possible with the present approach since we neglect changes in momentum amplitude (that is, in energy). In a recent study, Ref. [57] used a continuous time random walk model in momentum space to study particle acceleration, and power-law energy distributions were also derived.

To study the acceleration of particles, future research will consider how to derive superdiffusive shock acceleration from a fractional transport equation considering derivatives in both coordinate and momentum space.