2. A simple derivation of fractional Fick’s law

The diffusion equation can be obtained by combining the continuity equation,

(2.1) \begin{equation} \frac {\partial n}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = 0 \end{equation}

with the standard Fick’s law

(2.2) \begin{equation} \boldsymbol{J} = -\kappa \boldsymbol{\nabla }n, \end{equation}

where $\kappa$ is the diffusion coefficient; assuming that $\kappa$ is independent of coordinates, in the case of one spatial dimension, we obtain

(2.3) \begin{equation} \frac {\partial n}{\partial t} = \kappa \frac {\partial ^2 n}{\partial x^2} \, . \end{equation}

We now try to generalise Fick’s law to take into account the non-local character of the particle flux $\boldsymbol{J}$ and the scale-free properties of the particle random walk. For this purpose, we follow the logical steps adopted by Richard Feynman in his famous lectures, since we believe that Feynman’s approach is clearer from the physical point of view than using (1.2). Indeed, the assumption that the contribution of far away points to the non-local flux is proportional to the density derivative in $x'$ , as in (1.2), is not always clearly justified.

Consider the net flow in the $x$ direction that crosses a plane surface perpendicular to the $x$ axis. We have to count the number of particles that cross the plane in the direction of positive $x$ , which is $n_{-}v_{x}$ , and subtract from this value the number of particles crossing the surface in the negative $x$ direction, which is $n_{+}v_{x}$ . Here, $n_{-}$ is the number density of particles to the left of the point $x$ at which we intend to evaluate the net flux, $n_{+}$ is the number density of particles on the right side and $v_{x}$ is the mean particle velocity in absolute value along $x$ . The net flux density $J$ can be expressed as

(2.4) \begin{equation} J = (n_{-}-n_{+})v_{x} \equiv -(n_{+}-n_{-})v_{x}, \end{equation}

which is indeed (43.22) of Feynman, Leighton & Sands (Reference Feynman, Leighton and Sands1989). At this point, Feynman et al. (Reference Feynman, Leighton and Sands1989, vol. I, pp. 43–47) say

In the standard case of a local flux, for small distances $\Delta x$ , we can writeFootnote ¹

(2.5) \begin{align} (n_{+}-n_{-}) = \frac {\partial n}{\partial x} \ \Delta x = \frac {\partial n}{\partial x} \ 2\lambda . \end{align}

Thus, we can rewrite (2.4) as

(2.6) \begin{equation} J = -2\lambda v_{x} \ \frac {\partial n}{\partial x}. \end{equation}

Therefore, the flow of particles along $x$ is proportional to the derivative of the density, and the diffusion coefficient is given by $\kappa \simeq \lambda v_{x}$ .

Here, we would like to obtain a non-local flux which also includes the contribution to the flux of particles coming from ‘far away’ regions: how far do we mean? We mean from $x'=x+\xi$ and $x'=x-\xi$ , that is, from a range of distances $\pm \xi$ from the current position $x$ . Therefore, $n_{+}$ and $n_{-}$ in (2.4) are to be replaced by $n(x+\xi )$ and $n(x-\xi )$ . We introduce the small flux $\mathrm{d}J$ coming from a small extent $\text{d}\xi$ of positions as

(2.7) \begin{equation} \mathrm{d}J = - [{n(x+\xi ,t) - n(x-\xi ,t)}] \ v_x \ W(\xi ) \ \frac {\text{d}\xi }{\xi _0}, \end{equation}

where $\xi _0$ is a scale length introduced for dimensional reasons. We note that, in the case of non-local flux contributions, the use of (2.6) with the density derivative is not appropriate. We also introduced a dimensionless weight $W(\xi )$ , which takes into account the fact that, physically, the contribution to the flux of nearby regions should be larger than that of distant regions. This weight could be taken as either a Gaussian function or a bi-exponential decrease. For instance, in the context of laser-produced plasmas, for electron heat flux, Luciani et al. (Reference Luciani, Mora and Virmont1983) assume a weight proportional to the exponential decrease of an integral of density over distance; the same weight is adopted by Karpen & DeVore (Reference Karpen and DeVore1987) in the astrophysical context. For electron transport in the solar corona, Bian et al. (Reference Bian, Emslie and Kontar2017) derived from a kinetic equation in velocity space, assuming a finite mean free path, a weight $W(\xi )\sim \exp ({-{\sqrt {45}} \xi /\lambda })$ . In all those cases, the contribution to the flux becomes exponentially small for $\xi \gt \lambda$ . However, superdiffusion is basically related to the presence of power-law probability distributions for the free path lengths $\ell$ , which implies a diverging mean free path and a broader range of distances $\xi$ for the contribution to the flux than that corresponding to a bi-exponential weight. Therefore, to have a more decidedly non-local flux, we choose an inverse power law $W(\xi )\sim (\xi /\xi _0)^{-(1+\beta )}$ , so that $\mathrm{d}J$ can be written as

(2.8) \begin{equation} \mathrm{d}J = - [{n(x+\xi ,t) - n(x-\xi ,t)}] \ v_x \ \left ( \frac {\xi }{\xi _0} \right )^{-(1+\beta )} \ \frac {\text{d}\xi }{\xi _0}, \end{equation}

that is,

(2.9) \begin{equation} \mathrm{d}J = - v_x \xi _0^\beta \ \frac {n(x+\xi ,t) - n(x-\xi ,t)}{\xi ^{1+\beta }}\, \text{d}\xi \ . \end{equation}

We can see that the length $\xi _0$ separates the regions where the weight $W(\xi )\gt 1$ from those where $W(\xi )\lt 1$ . Assuming that $n(x)$ remains finite for $x\to \pm \infty$ and that $\beta \gt 0$ , we can extend the integration over $\xi$ from $0$ to $\infty$ to obtain

(2.10) \begin{equation} J(x,t) \simeq - v_x \xi _0^{\beta } \int _0^{\infty } \frac {n(x+\xi ,t) - n(x-\xi ,t)}{\xi ^{1+\beta }}\, \text{d}\xi \ . \end{equation}

This expression of the non-local flux depends on two new parameters, i.e. $\beta$ and $\xi _0$ . In a way similar to the power-law distributions of Lévy processes, $\beta$ determines the degree of non-locality, that is, how heavy are the tails of the non-local flux. Due to the singularity for $\xi \to 0$ , we require that $\beta \lt 1$ to avoid the divergence, considering that the density $n(x)$ is a smooth function so that the numerator in (2.10) will be proportional to $\xi$ . However, the value of $\xi _0$ contributes to determine the value of the fractional diffusion coefficient $\kappa _{\beta }$ , which is given as

(2.11) \begin{equation} \kappa _{\beta } = v_x \xi _0^{\beta } \, . \end{equation}

Despite the partial similarity to the standard diffusion coefficient, we note that the scale length $\xi _0$ does not correspond to the mean free path $\lambda$ (the latter being divergent in the case of superdiffusion). As we show later, $\xi _0$ has the same physical meaning as the scale parameter $\ell _0$ of Perri & Zimbardo (Reference Perri and Zimbardo2015), namely it is related to the minimum free path length such that for $\ell \gt \ell _0$ , the distribution of the free path lengths $\psi (\ell )$ decays as a power-law. Thus, we can write the fractional Fick’s law as

(2.12) \begin{equation} J(x,t) \simeq - \kappa _{\beta } \int _0^{\infty } \frac {n(x+\xi ,t) - n(x-\xi ,t)}{\xi ^{1+\beta }} \,\text{d}\xi \ , \end{equation}

with $0\lt \beta \lt 1$ .

We now proceed to show that the expression of the flux, (2.12), leads to a fractional diffusion equation with the Riesz derivative when $J(x,t)$ is inserted into the continuity equation. Let us consider

(2.13) \begin{align} \frac {\partial }{\partial x} \int _{0}^{\infty } \frac {n(x+\xi ,t)-n(x-\xi ,t)}{\xi ^{1+\beta }}\,\text{d}\xi &=\int _{0}^{\infty } \frac {\partial n(x+\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi\nonumber \\[5pt] &\quad - \int _{0}^{\infty } \frac {\partial n(x-\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi , \end{align}

where

(2.14) \begin{align} \int _{0}^{\infty } \frac {\partial n(x+\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi \ = \int _{0}^{\infty } \frac {\partial n(x+\xi ,t) }{\partial \xi } \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi . \end{align}

If we integrate by parts, we obtain for the ‘right’ contribution

(2.15) \begin{align} \int _{0}^{\infty } \frac {\partial n(x+\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi \ = \frac {n(x+\xi ,t)}{\xi ^{1+\beta }}\Big |_{0}^{\infty }+(1+\beta ) \int _{0}^{\infty } \frac {n(x+\xi ,t)}{\xi ^{2+\beta }}\,\text{d}\xi. \end{align}

In an analogous way, for the ‘left’ contribution, we have

(2.16) \begin{align} \int _{0}^{\infty } \frac {\partial n(x-\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi \ = -\int _{0}^{\infty } \frac {\partial n(x-\xi ,t) }{\partial \xi } \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi , \end{align}

that is,

(2.17) \begin{align} \int _{0}^{\infty } \frac {\partial n(x-\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi \ = - \left [ \frac {n(x-\xi ,t)}{\xi ^{1+\beta }}\Big |_{0}^{\infty }+(1+\beta ) \int _{0}^{\infty } \frac {n(x-\xi ,t)}{\xi ^{2+\beta }}\,\text{d}\xi \right ]. \end{align}

Using the relation (Bayın Reference Bayın2016)

(2.18) \begin{equation} \frac {f(x-\xi )}{\xi ^{1+\beta }}\Big |_{0}^{\infty } = -\displaystyle \lim _{\xi \to 0} \frac {f(x)}{\xi ^{1+\beta }} = -(1+\beta ) \int _{0}^{\infty } \frac {f(x)}{\xi ^{2+\beta }}\,\text{d}\xi , \end{equation}

we obtain

(2.19) \begin{align} \int _{0}^{\infty } \frac {\partial n(x+\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi \ = (1+\beta ) \int _{0}^{\infty } \frac {n(x+\xi ,t)-n(x,t)}{\xi ^{2+\beta }}\,\text{d}\xi \end{align}

and

(2.20) \begin{align} \int _{0}^{\infty } \frac {\partial n(x-\xi ,t) }{\partial x} \frac {1}{\xi ^{1+\beta }}\,\text{d}\xi \ = -(1+\beta ) \int _{0}^{\infty } \frac {n(x-\xi ,t)-n(x,t)}{\xi ^{2+\beta }}\,\text{d}\xi, \end{align}

which lead to

(2.21) \begin{align} \frac {\partial }{\partial x} J(x,t) \simeq -(1+\beta ) \kappa _{\beta } \int _{0}^{\infty } \frac {n(x+\xi ,t)-2n(x,t)+n(x-\xi ,t)}{\xi ^{2+\beta }}\,\text{d}\xi. \end{align}

We now recall that the Riesz derivative of order $\mu$ can be written as (e.g. Litvinenko & Effenberger Reference Litvinenko and Effenberger2014)

(2.22) \begin{align} \frac {\partial ^{\mu } n(x,t)}{\partial |x|^{\mu }} = \frac {\varGamma (1+\mu )}{\pi }\sin \left ( \frac {\pi \mu }{2} \right ) \int _{0}^{\infty } \frac {n(x+\xi ,t)-2n(x,t)+n(x-\xi ,t)}{\xi ^{1+\mu }}\,\text{d}\xi. \end{align}

Setting $\mu =1+\beta$ and using the properties of the Euler gamma function, we have $ \displaystyle \varGamma (2+\beta ) = (1+\beta ) \varGamma (1+\beta )$ . By comparing with (2.22), we can obtain from $\partial J{\kern-1.5pt/}\partial x$ the Riesz derivative by multiplying $J$ by a constant $A$ given as

(2.23) \begin{equation} A = \frac {\varGamma (1+\beta )}{\pi }\sin \left [ \frac {\pi }{2}(1+\beta ) \right ], \end{equation}

which is of order 0.3 for $0\lt \beta \lt 1$ . Finally, we have

(2.24) \begin{equation} J(x,t) = - \kappa _{\beta } \ \frac {\varGamma (1+\beta )}{\pi }\sin \left [ \frac {\pi }{2}(1+\beta ) \right ] \int _{0}^{\infty } \frac {n(x+\xi ,t)-n(x-\xi ,t)}{\xi ^{1+\beta }}\,\text{d}\xi . \end{equation}

We point out that (2.24) is equivalent to that proposed by Zanette (Reference Zanette1998), and coincides to that derived by Calvo et al. (Reference Calvo, Sánchez, Carreras and van Milligen2007) in terms of Riemann–Liouville fractional derivatives, except that those authors do not give an expression for $\kappa _{\beta }$ .

We can conclude that with the flux defined as in (2.24), we see that $ {\partial J}/{\partial x}$ is exactly equal to $ -\kappa _\beta {\partial ^{\mu } n(x,t)}/{\partial |x|^{\mu }}$ . Then, from the continuity equation, we obtain

(2.25) \begin{align} \frac {\partial n(x,t) }{\partial t} = \kappa _{\beta } \frac {\partial ^{1+\beta }}{\partial |x|^{1+\beta }} n(x,t), \end{align}

which is exactly the fractional diffusion equation for the number density of particles, having as solution a Lévy distribution of index $\mu = 1+\beta$ .

We notice that if our plasma domain is finite, we can stop the integration domain in (2.10) and (2.24) to the distances corresponding to our domain, including the case of asymmetric configurations, and still be able to compute a non-local flux. Another advantage of having the flux as in (2.24) is that the Riesz derivative is originally defined in Fourier space and, therefore, implies that the density $n(x,t)$ should admit a Fourier transform. However, (2.10) and (2.24) allow one to compute the flux even if $n(x=\pm \infty )$ does not go to zero, provided that the difference $n(x+\xi ) - n(x-\xi )$ remains finite.

2.1. Fractional diffusion coefficient $\kappa _\beta$ and scale length $\xi _0$

It is worth recalling that in the case of anomalous diffusion, the fractional diffusion coefficient $\kappa _\beta = v_x\xi _0^\beta$ does not coincide with the anomalous diffusion coefficient $\mathcal{D}_\alpha$ , which enters the relation defining anomalous diffusion, that is, $\langle \Delta x^2 \rangle =2 {\mathcal{D}}_\alpha t^{\alpha }$ . In fact, $\kappa _\beta$ and $\mathcal{D}_\alpha$ even have different physical dimensions. In the case of superdiffusion, the relation between these two constants has been given by Perri et al. (Reference Perri, Zimbardo, Effenberger and Fichtner2015), and can be written as

(2.26) \begin{equation} {\mathcal{D}}_\alpha = \frac {\varGamma (4-\alpha )}{\pi (\alpha -1)} \sin {\left [\frac {\pi }{2}(3-\alpha )\right ]} {v_x}^{\alpha - 1} \kappa _\beta, \end{equation}

which holds for superdiffusion.Footnote ² Furthermore, Perri & Zimbardo (Reference Perri and Zimbardo2012) and Zimbardo & Perri (Reference Zimbardo and Perri2013) have shown that ${\mathcal{D}}_\alpha$ can be written as

(2.27) \begin{equation} {\mathcal{D}}_\alpha = \frac {2(2-\alpha )}{(3-\alpha )(4-\alpha )} \varGamma (\alpha -1) \ \ell _0^{2-\alpha } v_x^{\alpha }, \end{equation}

where $\ell _0$ is the shortest free path length for which the free path distribution $\varPsi (\ell ,\tau )$ has a power law shape. Since $\alpha =2-\beta$ , a comparison of (2.26) and (2.27) shows that

(2.28) \begin{equation} \xi _0 = \left [ \frac {2\pi (2-\alpha ) \varGamma (\alpha )}{(3-\alpha ) \sin {\left [({\pi }/{2})(3-\alpha )\right ]} \varGamma (5-\alpha ) } \right ]^{1/\beta } \ell _0 \end{equation}

with the numerical factor being of order one for $1.2 \lt \alpha \lt 1.6$ . These values for $\alpha$ are typically found for the transport of energetic particles upstream of interplanetary shocks (Zimbardo, Prete & Perri Reference Zimbardo, Prete and Perri2020). Therefore, we can see that the scale length $\xi _0$ represents the distance which, in the non-local flux given by (2.8), separates regions with larger weight $W(\xi )\gt 1$ from regions with smaller weight $W(\xi )\lt 1$ , and is of the same order of the shortest free path length for which the free path distribution has a power law shape. Notably, Perri & Zimbardo (Reference Perri and Zimbardo2012), Zimbardo & Perri (Reference Zimbardo and Perri2013), Perri et al. (Reference Perri, Zimbardo, Effenberger and Fichtner2015) and Perri & Zimbardo (Reference Perri and Zimbardo2015) have given a method to extract $\ell _0$ from the measured profiles of energetic particles at shock waves.

3. Numerical computation of the non-local flux

To compute the value of the diffusive flux given by the fractional Fick’s law, we need to specify the number density of particles $n(x,t)$ . Indeed, the determination of $J(x,t)$ , apart from a constant, reduces to the calculation of the following integral:

(3.1) \begin{equation} I(x,t) = \int _{0}^{\infty } \frac {n(x+\xi ,t)-n(x-\xi ,t)}{\xi ^{1+\beta }}\,\text{d}\xi. \end{equation}

For practical applications, we aim to also address the cases when the analytic form of $n(x,t)$ is not known in advance; instead, its values are known at a finite number of points of coordinates $x_{j}$ . In such cases, the integral expression (3.1) can be evaluated using the trapezoidal formula, which is convenient since it is not very restrictive in terms of the regularity of the functions we intend to integrate. Thus, it can be effectively used when dealing with experimental data and spacecraft data, where the behaviour of the physical quantities of interest can be quite irregular. Obviously, the unbounded integration interval of (3.1) has to be truncated and an upper bound shall be introduced. Hence, instead of $[0,\infty ]$ , we will have the bounded interval $[0,\xi _{ \mathrm{max}}]$ . We note that in correspondence of $ \xi = 0$ , the integrand function diverges. To deal with such singularity, it is convenient to write the integral (3.1) as the sum of the following two contributions:

(3.2) \begin{align} I(x,t) &= I_{1} + I_{2} = \int _{0}^{\varDelta } \frac {n(x+\xi ,t)-n(x-\xi ,t)}{\xi ^{1+\beta }}\,\text{d}\xi \nonumber \\&\quad +\int _{\varDelta }^{\xi _{ \mathrm{max}}} \frac {n(x+\xi ,t)-n(x-\xi ,t)}{\xi ^{1+\beta }}\,\text{d}\xi \ \end{align}

and to treat them separately. Assuming that the distance $\varDelta$ is much smaller than the typical length scales of the density $n(x)$ under consideration, we can adopt the following approximation inside the singular term $I_{1}$ :

(3.3) \begin{align} \frac {n(x+\xi ,t)-n(x-\xi ,t)}{\xi } \approx 2n'(x,t), \end{align}

where the first derivative of the density is

(3.4) \begin{align} n'(x,t) = \frac {\partial n}{\partial x} = \frac {n(x+\varDelta ,t)-n(x-\varDelta ,t)}{2\varDelta }+o(\varDelta ^{2}) \, . \end{align}

Therefore, the first integral can be written as

(3.5) \begin{align} I_{1} = \int _{0}^{\varDelta } 2n'(x,t)\frac {1}{\xi ^{\beta }}\,\text{d}\xi , \end{align}

which yields, after trivial calculations,

(3.6) \begin{equation} I_{1} = \frac {1}{1-\beta } \frac {n(x+\varDelta ,t)-n(x-\varDelta ,t)}{\varDelta ^{\beta }} \, . \end{equation}

The term $I_{2}$ is easily evaluated by means of the trapezoidal formula, which is implemented by running a programme written in Fortran 90.

As a validation of the above-mentioned method, we compute the diffusive current $J(x)$ associated with the non-normalised Gaussian density profile $ \displaystyle n(x) = \exp {(-x^{2})}$ , which corresponds to the case considered by Zanette (Reference Zanette1998). In our numerical set-up, the $x$ axis is discretised as follows:

(3.7) \begin{align} x_{j} = j\varDelta , \quad j = -N_{x},\ldots ,N_{x}, \end{align}

so that a total number of $ 2N_{x}+1$ equally spaced nodal points is considered. Here, $J(x)$ is evaluated in correspondence of $ 2N_{1} + 1$ nodal points, with $N_1\lt N_x$ . For the present case, $\varDelta = 0.1$ . As mentioned previously, the integration interval has to be reduced to a bounded domain $[0,\xi _{N_{2}}]$ , which means that the number of $ \xi$ values involved in the computation of the integral (3.1) is $N_{2}+1$ , including the singularity contribution, with $N_{1} + N_{2} \leqslant N_{x}$ .

Figure 1. Fractional flux, normalised to $\kappa _\beta$ , for a Gaussian density profile, for several values of $\beta$ (see legend). The normal first derivative of the Gaussian is also shown by the dashed blue line. The inset shows the same values in log–log axes, for $1\lt x\lt 400$ .

Figure 2. Values of the fractional flux, normalised to $\kappa _\beta$ , for a Gaussian function of variance $\sigma ^2$ and for $\beta =0.5$ . Several values of $\sigma$ are used (see legend); as expected for normal derivatives, also for the fractional flux, the larger $\sigma$ , the wider and less peaked the flux.

Figure 1 shows the results of the numerical calculation of $ J(x)/\kappa _\beta$ for various values of the parameter $\beta$ , which is related to the degree of non-locality, as explained before. We can notice that, as $ \beta$ decreases, the tails of the current profile become higher, with slower decays (namely the non-locality is increasing and the process becomes increasingly more superdiffusive). For $|x|\gg 1$ , the tails follow a power law whose slope depends on $ \beta$ , $J\propto |x|^{-(1+\beta )}$ , as clearly shown in the log–log scale inset in figure 1. Such a feature, noted by Zanette (Reference Zanette1998), is in agreement with the fact that diminishing $\beta$ is equivalent to increasing the value of the weight $ \displaystyle {\xi ^{-(1+\beta )}}$ inside (3.1) for large $\xi$ . Figure 1 shows that as $ \beta$ approaches 1, the flux $J(x)$ tends to the first derivative of $ n(x)$ (apart from the sign), which means that the numerical method we are using is correct and consistent. Moreover, the inset shows the very fast decay of the first derivative of $n(x)$ (dashed curve in figure 1) in contrast to the fractional flux, which decays slowly as a power law. The numerical calculation is also carried out assuming the normalised Gaussian profile $ \displaystyle n(x) = ({1}/{\sqrt {2\pi } \sigma }) \exp { (-({x^{2}}/{2 \sigma ^{2}} ))}$ to evaluate the influence of the length parameter $ \sigma$ on the values of $J(x)$ for a fixed $ \beta =0.5$ . The results are presented in figure 2, which shows that a larger $\sigma$ results in a smaller peak value of $|J|$ and a broader range of $|x|$ , where $J$ is appreciably different from zero, as expected for the first derivative of a Gaussian function.

4. Uphill transport of energetic particles at collisionless shock waves

As an application, now we compute the flux associated with the density profile of energetic particles around a collisionless shock wave, which is a powerful accelerator of particles (e.g. Bell Reference Bell1978; Drury Reference Drury1983; Giacalone Reference Giacalone2012). Investigating the flux of energetic particles around a shock wave is relevant because the most popular theory of cosmic ray acceleration at shocks, based on a diffusive–advection equation and called diffusive shock acceleration (DSA), foresees a density profile, for energetic particles, which is an exponential growth in the region upstream of the shock and a constant, flat profile in the region downstream of the shock (Giacalone Reference Giacalone2012). However, this is at variance with what is observed around shocks in interplanetary space and, in some cases, around shocks of supernova remnants (e.g. Perri & Zimbardo Reference Perri and Zimbardo2007, Reference Perri and Zimbardo2009; Perri, Amato & Zimbardo Reference Perri, Amato and Zimbardo2016). In practice, a long power-law-like density profile is observed upstream and a decreasing, non-flat density profile downstream (Perri & Zimbardo Reference Perri and Zimbardo2012; Perrone et al. Reference Perrone2013). This creates a tension with the solutions of the standard Fick’s law (Drury Reference Drury1983; Effenberger et al. Reference Effenberger2025).

If we assume a fractional diffusion process based on (2.24), the non-local flux $J(x)$ is an integral operator depending on the global properties of the density profile $n(x)$ , which can give rise to uphill transport (e.g. del Castillo-Negrete Reference del Castillo-Negrete2006). In this case, we can also get a negative fractional diffusive flux in the presence of a negative density gradient.

Let us assume to be in the shock frame, with the shock at $x=0$ and the plasma velocity $V_{\mathrm{bulk}}$ going in the positive $x$ direction, so the upstream side is found for $x\lt 0$ and the downstream side for $x\gt 0$ . The source of energetic particles is at the shock, and, in addition to performing a random walk, particles are convected downstream by the plasma flow. The density profile of energetic particles is obtained from test particle numerical simulations (Prete et al. Reference Prete, Perri and Zimbardo2019, Reference Prete, Perri and Zimbardo2021). In this simulation, particles are injected at $x=0$ and make a one-dimensional motion according to

(4.1) \begin{equation} \text{d}x = V_{\mathrm{bulk}} \text{d}t + v_{\mathrm{ran}} \,\text{d}t. \end{equation}

Here, $V_{\mathrm{bulk}}=V_1$ upstream and $V_2$ downstream, with $V_1=600$ km s $^{-1}$ and $V_2=200$ km s $^{-1}$ . The random velocity $v_{\mathrm{ran}}$ can be modelled to reproduce both normal diffusion and superdiffusion, by a suitable choice of the scattering times $\tau$ ; in particular, for superdiffusion, $\tau$ is distributed according to a power-law, $\psi (\tau )\sim \tau ^{-(1+\mu )}$ (the reader is referred to Prete, Perri & Zimbardo (Reference Prete, Perri and Zimbardo2019, Reference Prete, Perri and Zimbardo2021) for a full description of the code. See also Effenberger et al. (Reference Effenberger, Aerdker, Merten and Fichtner2024)). We recall that the superdiffusion exponent of the mean square displacement $\alpha$ is related to $\mu$ and $\beta$ by $\alpha =3-\mu =2-\beta$ (Perri & Zimbardo Reference Perri and Zimbardo2015). Typically, $10^6$ particles are injected and their motion is integrated within a very large box that allows us to reach a steady-state density shape, namely the simulation is only marginally affected by the presence of lost particles. After some trials, the simulation box size has been fixed to $\pm 5\times 10^9$ m (almost one-tenth of an astronomical unit), while the vertical axis is in arbitrary units.

Figure 3. (a) Density profile $n(x)$ of energetic particles accelerated at a shock in $x=0$ , with the upstream region on the left and the downstream region on the right, in the case of normal diffusion. (b) Corresponding fractional flux obtained for $\beta =0.25$ . It can be seen that the flux is negative both upstream and downstream. If we used the local Fick’s law, the downstream flux would be zero.

4.2. Superdiffusion

The density as a function of the distance from the shock $x$ is shown in figure 4(a) for the case of superdiffusion with an exponent $\alpha =1.75$ , which corresponds to $\beta =0.25$ . The upstream density profile corresponds to a power-law decrease, that is, the upstream density is larger than in the case of normal diffusion, while the downstream density is decreasing. Both these features are in agreement with the studies of superdiffusion of shock-accelerated particles by Perri & Zimbardo (Reference Perri and Zimbardo2007, Reference Perri and Zimbardo2008, Reference Perri and Zimbardo2012); Zimbardo & Perri (Reference Zimbardo and Perri2013), Litvinenko & Effenberger (Reference Litvinenko and Effenberger2014), Aerdker et al. (Reference Aerdker, Merten, Effenberger, Fichtner and Tjus2025). Figure 4(b) shows the corresponding non-local flux: the upstream values are negative but smaller (in absolute value) than those found in the case of normal diffusion, see figure 3. This can be interpreted as the result of the larger upstream densities due to superdiffusion, so that the upstream gradients are smaller. In contrast, in the downstream region, the flux is negative (i.e. towards the left) up to approximately $2.8 \times 10^9$ m, that is, in a region where the density is decreasing, as shown in panel (a). This gives evidence of uphill transport in the downstream region, and can be understood considering again the integrand in (2.24): when $x$ corresponds to downstream and $n(x-\xi )$ is to be evaluated in the upstream region, its contribution is smaller than that of $n(x+\xi )$ , so a negative contribution to the flux is obtained. Interestingly, for $x$ larger than approximately $2.8 \times 10^9$ m, the flux becomes positive, that is, downhill: being far downstream, the effects of the upstream region become of secondary importance and the decreasing density in the downstream region prevails. From the microscopic point of view, i.e. particle motion, we consider that uphill transport is due to the fact that, starting from downstream, particles performing a Lévy walk can make long displacements also in the direction upstream, that is, can move upwind more easily than particles performing a normal random walk.

Figure 4. (a) Density profile $n(x)$ of energetic particles accelerated at a shock in $x=0$ , with the upstream region on the left and the downstream region on the right, in the case of superdiffusion. (b) Corresponding fractional flux obtained for $\beta =0.25$ . It can be seen that the flux is negative both upstream and downstream, so that uphill transport is obtained in the downstream region.

The finding of an uphill transport in the downstream region, in connection with a decreasing density profile $n(x)$ , is consistent with the continuity equation (2.1) under steady-state conditions when the advection flux $nV_{\mathrm{bulk}}$ is also considered,

(4.2) \begin{equation} \frac {\partial }{\partial x} (J + n V_{\mathrm{bulk}}) = 0 \, . \end{equation}

For the total flux to be constant, an increasing flux $J$ (that is negative for uphill transport), similar to that in figure 4, is needed when $n(x)$ is decreasing. However, the ordinary Fick’s law is not consistent with the downstream increasing diffusive flux, in the presence of a decreasing number density. In other words, the case of a constant downstream density profile is not consistent with the continuity (4.2) and non-local flux: a constant density profile would require a constant flux $J$ downstream, but this is not found in the case of non-local flux, as shown in figure 3. Therefore, spacecraft observations around a shock wave of a decreasing downstream density profile of energetic particles, as shown in figure 4, are an indication of non-local transport, that is, of particle superdiffusion.

These results show that uphill transport can be obtained in a simple and natural way in a collisionless plasma, and are useful for understanding energetic particle transport and acceleration at astrophysical shocks.