On the Fractional Diffusion-Advection Equation for Fluids and Plasmas

: The problem of studying anomalous superdiffusive transport by means of fractional transport equations is considered. We concentrate on the case when an advection ﬂow is present (since this corresponds to many actual plasma conﬁgurations), as well as on the case when a boundary is also present. We propose that the presence of a boundary can be taken into account by adopting the Caputo fractional derivatives for the side of the boundary (here, the left side), while the Riemann-Liouville derivative is used for the unbounded side (here, the right side). These derivatives are used to write the fractional diffusion–advection equation. We look for solutions in the steady-state case, as such solutions are of practical interest for comparison with observations both in laboratory and astrophysical plasmas. It is shown that the solutions in the completely asymmetric cases have the form of Mittag-Lefﬂer functions in the case of the left fractional contribution, and the form of an exponential decay in the case of the right fractional contribution. Possible applications to space plasmas are discussed.


Introduction
A fundamental goal of fluid and plasma physics is to understand transport properties, given that transport governs, for example, the dispersion of pollutants in the atmosphere and the oceans, the confinement of plasmas in laboratory devices, and the acceleration and propagation of energetic particles in astrophysics [1][2][3]. Non-diffusive anomalous transport, that is, different from standard Brownian diffusion, is observed in a large variety of physical systems [4][5][6], including both laboratory and astrophysical plasmas [1,3,[7][8][9][10][11][12][13]. This is characterized by a nonlinear growth with time t of the mean square displacement, ∆x 2 ∝ t α , and both superdiffusive (α > 1) and subdiffusive (α < 1) regimes are observed in many systems [4,5,14]. In fluid and plasma physics, the presence of turbulence-with its broad range of scales involved in the eddy size distribution-is one of the main factors causing anomalous transport [11,[15][16][17][18][19][20]. Further, in magnetized plasmas, an intriguing role that is yet to be fully understood is due to the scale-free distribution of pitch-angle scattering times [21][22][23].
In laboratory plasmas, both numerical [8,9,24] and experimental studies [7,25] point to superdiffusive and non-local transport. The uphill transport observed in laboratory plasmas [26], the nearly scatter-free propagation of solar energetic particles [27,28], as well as the 30 AU mean free paths deduced for particles accelerated at the solar wind termination shock [29] imply that the local Fick's law, J = −k∇n, is not appropriate for all situations. Here, J is the diffusive flux, n is the number density, and k is a suitable diffusion coefficient. The relevance of non-locality phenomena and non-local transport in laboratory plasmas has recently been reviewed by [30,31], and an integral non-local relation between the heat flux and the temperature gradient is given by Equation (1) of [30], while a non-local relation between the particle flux and the density gradient is given by Equation (42) of [31]. In astrophysical plasmas, numerical simulations of particle transport in the presence of magnetohydrodynamic turbulence, like those of Ref. [18,32], and the analysis of energetic particle fluxes upstream of shock waves [12,13,33] show that energetic particle transport can be superdiffusive.
From a microscopic point of view, superdiffusive transport is due to the scale-free non-local character of the particle random walk process, which typically involves long-range correlations and a Lévy walk. The Lévy walk is composed of displacements which, rather than having a typical mean free path, have a power-law distribution of lengths with a diverging second-order moment [17,19,34]. From a macroscopic (i.e., fluid) point of view, the description of superdiffusion in terms of fractional diffusion equations is often adopted [4,35,36]. For instance, a common approach is to change the second-order derivative in the standard diffusion equation with a fractional derivative of order µ, with 1 < µ < 2 [4,5,9,[35][36][37][38][39], where n(x, t) is the particle density (or the tracer density), and k µ is an anomalous diffusion constant with dimensions L µ /T. Several different definitions of fractional derivatives are found in the literature (see [40]). For the fractional diffusion equation, a popular choice is the completely symmetric Riesz fractional derivative [1,4,36,39]. The latter can be expressed as which clearly shows the non-local character of transport, since the density of particles at positions x + ξ and x − ξ is involved in the transport equation. Here, Γ is the Euler gamma function. Assuming a delta function as an initial distribution, upon Fourier transforming, Equation (1) yields the characteristic function exp(−k µ t|q| µ ) (where q is the wavenumber), which is the Fourier representation of the Lévy distribution [4,36]. As is known, Lévy distributions having power law tails are one of the prime tools to study superdiffusion. On the other hand, other definitions of fractional derivatives are frequently used, such as the Riemann-Liouville and Caputo derivatives (see below). A feature of the fractional diffusion Equations (1) and (2) is that the Riesz derivative requires the integration of the particle density over space from −∞ to +∞-that is, it is appropriate for an infinite domain. However, in many cases we are interested in understanding anomalous transport in a bounded system (including boundaries). For instance, this is the case of wall impurities diffusing into the scrape-off layer of a laboratory plasma device [25,41], of hot plasma diffusing away from the axis of a toroidal device along the minor radial direction r [42], of shock-accelerated particles propagating upstream of an astrophysical shock, and of solar energetic particles propagating from the Sun to the Earth. This means that because of the boundaries which are limiting the integration domain, fractional derivatives other than Riesz's have to be used. Ref. [42] has shown that this problem can be treated by using the Caputo fractional derivatives.
In this paper, we give the form of the fractional diffusion-advection equation in terms of the Caputo derivative for the left side, where the boundary is found, and of the Riemann-Liouville derivative for the right side, which corresponds to the unbounded domain. In Section 2, we show that in the steady-state case the solutions can be expressed in terms of the Mittag-Leffler functions, and of exponential functions, respectively, in Section 3. The conclusions and possible applications to space plasmas are given in Section 4.

The Fractional Diffusion-Advection Equation
The problem of transport in the presence of both anomalous diffusion and advection is of great interest for many applications. This includes the case of energetic particle transport in front of a shock wave, considered by [39,43], where the plasma turbulence is generated by backstreaming particles upstream of the shock [44,45]; the case of plasma transport in the presence of turbulence and of an external electric field, as in the Earth's magnetotail; and the case of plasma and impurity transport in the presence of turbulence and an imposed velocity field, as in recent plasma medical applications [46]. Fractional diffusion-advection equations have already been considered in the literature (e.g., [4]), but mostly considering the time fractional derivative appropriate for the subdiffusive case. Conversely, here we are interested in the superdiffusive case, and consider space fractional derivatives which explicitly consider the non-local character of transport [5,9,36,42]. Space fractional derivatives accurately correspond to Lévy flights, and in the case of particles having mass (i.e., particles which cannot propagate at arbitrarily large velocities), the distinction between Levy flights and Levy walks is to be made-see Refs. [6,34,47] for a discussion. In this paper we study the transport of plasma particles with a simplified analytical model, that is, by means of a diffusion-advection equation in one dimension. This is clearly an oversimplification for many plasma systems, nevertheless, this simple modeling can show the potentiality of anomalous superdiffusive transport and of fractional equations when trying to reproduce the observations.
In order to introduce the advection term, let us start from the flux-gradient relationship which is expressed by Fick's law In the cases of non-local transport the corresponding fractional Fick's law in one dimension can be written as [37,42,48,49] where β = µ − 1 (see Ref. [37] for a discussion on the form needed for D β n in order to recover Equation (1)). We note that here D β represents a generic fractional derivative operator, whose form is specified below. Since 1 < µ < 2 for superdiffusion, we have that 0 < β < 1. From the continuity equation we obtain the transport equation in the absence of an external velocity field as ∂n ∂t Considering the terms corresponding to advection and fractional transport means to add the advection flux nV to Equation (5): where V is the plasma velocity. This is the fractional diffusion-advection equation that we would like to study. Time-dependent solutions to this equation, also including a delta-function source term at x = 0, have been given by [39]. We now express the non-local diffusive flux due to particles coming both from the region 0 < x < x to the left of the current position x and from the unbounded domain Indeed, in order to treat domains with a boundary at x = 0, we adopt the left Caputo derivative where d f /dx represents the first (integer) derivative of f (x ). Also, to treat domains unbounded on the right, we adopt the right Riemann-Liouville derivative for x > x: (see Ref. [40] for a review of the most common forms of fractional derivatives). Here, we consider the order of differentiation β < 1, since this is appropriate for superdiffusion described by Fick's law, but higher-order fractional derivatives are easily defined [4,36,40,42]. We point out that if, on the right side, the right Caputo derivative (i.e., with the integer derivative inside the integral) on an open domain would be used, this would give the same result as the right Riemann-Liouville derivative [50].
In other words, if one wishes, it is possible to write both the left and the right contributions in terms of Caputo derivatives; this would give a flux-gradient relationship in the same form as Equation (1) of Ref. [30] and Equation (42) of Ref. [31], that is, with the gradient of the temperature/density inside the integral. Note that Equation (7) above encompasses a single diffusion constant k β for both the left and right contributions, but expressions with different diffusion constants, that is, different weights for the two sides, can be introduced [31,42,48,49]. Here we keep a single constant since below we consider the completely asymmetric fractional Fick's law [48,49], that is, the weights for the two sides might be represented as (1, 0) and (0, 1). Basically, we are considering the same kind of particle transport in the two domains.

Steady-State Solutions of Fractional Diffusion-Advection Equations
The solution of the above fractional diffusion equation gives information on the time evolution of the particle density n(x, t); for instance, if the initial distribution is a delta function, at later times the density becomes a Lévy distribution of time-increasing width, shifted towards the direction of the flow [5,37,39,51]. Less attention has been given to steady-state solutions of the transport problem, corresponding to studying the flux-gradient relationship which is expressed by the fractional Fick's law. The time independence, ∂n/∂t = 0, allows a first integration to be made in Equation (6) to obtain J = const.
To fix the ideas, we consider the direction x pointing to the right, and we assume that charged particles are subject to advection due to a fluid motion to the left, which we indicate by V < 0 (since the x-axis is pointing to the right). In reality, V will change with the position in most cases, but for the purpose of illustration we keep it constant. Further, charged particles are influenced by magnetic and plasma turbulence, which create a diffusive or a superdiffusive motion. Let us first consider the normal diffusive case, so that the total particle flux J in the x direction can be written as We further assume that in steady state the advection and diffusion balance each other, that is, which for V < 0 and x > 0 yields an exponential decay, n(x) = n 0 exp(−|V|x/k). Let us now consider the superdiffusive case: in the one-dimensional case, let us write the density flux in terms of fractional derivatives as in Equation (7), where k β is an anomalous transport coefficient, with dimensions L 1+β /T, different from those of the normal diffusion coefficient. For the total flux, we proceed as above in the diffusive case, and we have We point out the following property of this equation: for both the Caputo derivative and the Riemann-Liouville derivative with terminals going from x to ∞, the derivative of a constant is zero. This means that we can find solutions of Equation (12) in the non-homogeneous case when a net flow is present, J = 0, in terms of a constant functionn = J/V to be added to the solutions of the associated homogeneous equation.
Balancing again the advective flux with the diffusive flux, we have J = 0. In order to illustrate the effects of either the left contribution or the right contribution to the fractional flux, we now look for solutions in the completely asymmetric cases, that is, considering either the left Caputo derivative only, or the right Riemann-Liouville derivative only. Asymmetric fractional transport has been considered by several authors [31,42,48,49], and asymmetry gives the possibility "to tune the strength of transport coming from one direction with respect to the other" [31].

Solutions with the Left Contribution Only
Starting from Equation (12) with J = 0, we keep only the left Caputo fractional derivative and we can write The solution of this linear fractional differential equation for x ≥ 0 is given by the Mittag-Leffler functions E β (−ζ β ) [4,36,50,52], with the dimensionless variable ζ > 0 defined via ζ β = |V|x β /k β , that is The fact that the Mittag-Leffler functions are solutions of the above equation is shown in, for example, [50], and also by [43] for the case when the x-axis is pointing in the same direction as the flow velocity V. In addition, it can be easily verified by recalling the effect of fractional derivatives on power functions f (x) = x γ for γ > −1, γ = 0, and x > 0, that is, and differentiating Equation (14) term-by-term. Here we can recall the limiting forms for |V|x β /k β 1 corresponding to a stretched exponential, and for |V|x β /k β 1 corresponding to a power law decay. Sample solutions are shown in Figures 1-3 for various β values. Notably, the power-law decay of Equation (17) corresponds to the results already found by [12,13] for the energetic particle profiles upstream of shock waves in the heliosphere, by considering the propagator for Lévy walks [34]. We recall that the Mittag-Leffler functions are the general solutions of fractional ordinary differential equations; they are defined as [9,36,50,53,54]: where β > 0 and z varies in the complex plane (although we can restrict ourselves to the real axis).
Clearly, E β (z) is a generalization of the exponential function and it reduces to the exponential for β = 1, in which case we have Γ(βn + 1) = n! We also note that E β (0) = 1 for any β. Generalized expressions of the Mittag-Leffler functions depending on two and three indices are also introduced.  Compared to the case of normal diffusion, which would give an exponential profile, we can see that in the case of superdiffusion and advection the density profile obtained by considering only the left contribution to the flux is steeper close to the source of particles at x = 0, and more shallow than an exponential farther away. This is similar to the behavior found for energetic particles propagating upstream of interplanetary shock waves [12,13,45,51,55] and upstream of supernova remnant shocks [56,57].

Solutions with the Right Contribution Only
We now look for solutions of Equation (12) in the completely asymmetric case when only the contribution of the domain to the right of x is taken into account. We have Remembering the minus sign implicit in the right Riemann-Liouville derivative of order β < 1, we can re-write 1 whose solution is as can be verified by direct substitution and remembering that with a > 0 a constant. We note that if the right Caputo derivative were used instead of the right Riemann-Liouville derivative, an identical result would be obtained, in agreement with the fact that the two derivatives only differ in the presence of a boundary (typically at x = 0) [50].
The argument of the exponential decay above can be written in terms of the dimensionless variable ζ = (|V|/k β ) 1/β x defined above, n(x) = n(0) exp(−ζ) . (23) Even if the dimensionless variable ζ is the same, these solutions are different from those obtained considering the left contribution only, that is, n(x) = n(0)E β (−ζ β ). Indeed, the Mittag-Leffler functions are steeper close to the boundary at x = 0, and are flatter than an exponential decay for large ζ-see Equations (16) and (17) and Figures 2 and 3 (of course, those differences vanish when β = 1 and transport is normal). The differences in the solutions obtained considering the left and right contributions separately confirm the richness of behaviors that can be described by asymmetric fractional derivatives [31,42,48,49]. It is instructive to compare our results to those reported by del-Castillo-Negrete [42]. He considered anomalous heat transport for a plasma column with the magnetic axis (the high-density, high-temperature side) at x = 0 and the outer boundary at x = L (the low-density, low-temperature side). Transport was described by means of a normal diffusive term, a left fractional operator, and a right fractional operator, and a space-dependent heat conductivity was also used. In his Figure 4 he shows the steady-state temperature profiles for various cases, including the completely asymmetric cases: it can be seen there that the temperature profile obtained by the right fractional operator was much steeper than that obtained from the left fractional operator. This is similar to our results for the density for ζ > 1, as can be seen from Figures 1-3. Conversely, for ζ → 0 we find a steeper profile when only the left fractional operator is taken into account: such a feature is not reported by [42], probably because in that paper the heat conductivity is adjusted to go to zero for x → 0.
We further notice that both the solutions for the left and right fractional operators depend on the dimensionless variable ζ = (|V|/k β ) 1/β x. A typical lengthL(β) for the density variation can be obtained by setting ζ = 1, which givesL This length is of the order of the distance from x = 0, at which the density decreases to 1/e of its value at x = 0. As can be seen,L(β) sensitively depends on the order of differentiation 0 < β < 1, which is related to the anomalous diffusion exponent α as β = 2 − α [12,13,33].

Discussion and Conclusions
In this paper we addressed the problem of studying anomalous superdiffusive transport by means of fractional transport operators of order β < 1, which generalizes the standard Fick's law. We concentrate on the case when an advection flow is present, since this corresponds to many actual plasma configurations, and on the case when a boundary is also present. In agreement with previous studies, we propose that the presence of a boundary can be taken into account by using the Caputo fractional derivatives. Therefore, we write the diffusive flux as the sum of a left contribution (left being the side of the boundary at x = 0) expressed by the left Caputo derivative and a right contribution (right being the side of the unbounded domain) expressed by the right Riemann-Liouville derivative. To obtain a simple set of solutions, we considered the completely asymmetric fractional Fick's law, which considers either only the left contribution to the flux or only the right contribution to the flux [31,42,48,49]. Within these limitations, we looked for solutions in the steady-state cases, as such solutions are of practical interest for comparison with observations both in laboratory and in astrophysical plasmas. It is shown that the solutions for the left contribution only had the form of Mittag-Leffler functions, which can be seen as a generalization of the exponential function. Since these functions can have many applications in fluid and plasma physics, some typical properties on the real axis are described. Instead, the solutions for the right fractional operator were only given by an exponential decay, with the typical scale of decay sensitively depending on β.
We recall that several definitions of fractional derivatives are available [36,40,50], and that this variety reflects the enhanced possibility offered by fractional derivatives to adapt the non-local flux operator, corresponding to the fractional Fick's law, to the diverse physical systems under consideration. In particular, in space plasmas the density of high-energy particles upstream of a shock wave is found to exhibit a sharp decrease close to the shock, and a slow power-law decay far upstream [12,13,33,43,47,51]. Such energetic particle density profiles are very similar to the properties of the Mittag-Leffler functions shown above. This suggests that the left Caputo derivative could be giving the main contribution to the flux upstream of the shock, with the contribution to the flux from the right domain being negligible. We can tentatively justify this interpretation by saying that we considered only the non-local flux coming from the high-density side-that is, the left side in the considered case of negative (leftward) advection velocity.
Finally, we suggest for future work the possibility to interpret the observations on the spatial distribution of multiply charged oxygen ions in the Earth's magnetotail [58] as the result of the advective Earthward flow due to the cross-tail electric field (giving rise to an approximately constant Earthward velocity), and of the turbulent superdiffusion due to the presence of a hierarchy of vortices and magnetic fluctuations typical of turbulence, as regularly observed in the magnetotail [59][60][61].
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.