Issue 
A&A
Volume 607, November 2017



Article Number  A7  
Number of page(s)  9  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201731179  
Published online  30 October 2017 
Fractional Parker equation for the transport of cosmic rays: steadystate solutions
^{1} Dipartimento di Fisica, Universitá della Calabria, Ponte P. Bucci, Cubo 31C, 87036 Rende CS, Italy
email: gaetano.zimbardo@fis.unical.it
^{2} International Space Science Institute, 3012 Bern, Switzerland
^{3} Bay Area Environmental Research Institute, Petaluma, CA 94952, USA
^{4} Department of Physics and KIPAC, Stanford University, Stanford, CA 94305, USA
^{5} Institut für Theoretische Physik IV, RuhrUniversität Bochum, Universitätsstraße 150, 44780 Bochum, Germany
Received: 16 May 2017
Accepted: 4 August 2017
Context. The acceleration and transport of energetic particles in astrophysical plasmas can be described by the socalled Parker equation, which is a kinetic equation comprising diffusion terms both in coordinate space and in momentum space. In the past years, it has been found that energetic particle transport in space can be anomalous, for instance, superdiffusive rather than normal diffusive. This requires a revision of the basic transport equation for such circumstances.
Aims. Here, we extend the Parker equation to the case of anomalous diffusion by means of fractional derivatives that generalize the usual secondorder spatial diffusion operator.
Methods. We introduce the left and right Caputo fractional derivatives in space. These derivatives are one of the tools used to describe anomalous transport. We consider the case of steadystate solutions upstream and downstream of a planar shock.
Results. We obtain an estimate of the particle acceleration time at shocks in the case of superdiffusion. An analytical solution of the steadystate fractional Parker equation is given by the MittagLeffler functions, which correspond to a powerlaw profile for the energetic particle intensity far upstream of the shock, in agreement with the results obtained from a probabilistic approach to superdiffusion. These functions also correspond to a stretched exponential close upstream of the shock.
Conclusions. These results can help to model more precisely the measured fluxes of energetic particles that are accelerated at both interplanetary shocks and supernova remnant shocks.
Key words: diffusion / acceleration of particles / shock waves / cosmic rays / methods: analytical / plasmas
© ESO, 2017
1. Introduction
The acceleration and transport of energetic particles in space plasmas and in astrophysical environments can be described by the socalled Parker equation (Parker 1965; Drury 1983; Zank et al. 2000; Moraal 2013), which is an equation for the evolution of the particle distribution function obtained under the assumption of nearisotropy in velocity space. In the case of a onedimensional system, for example in the presence of an infinite planar shock, ∂/∂y = ∂/∂z = 0, and the Parker equation reduces to (e.g., Drury 1983; Lagage & Cesarsky 1983a; Moraal 2013; Giacalone 2013) (1)where f(x,p,t) represents the omnidirectional distribution function of energetic particles with momentum p at position x and time t; x is the coordinate normal to the planar shock, U is the plasma (fluid) velocity in the x direction, k(x,p) is the spatial diffusion coefficient, D_{pp} the diffusion coefficient in momentum space, and Q is a source term. The various terms in the Parker equation represent, in the order from left to right, the time changes of f, the changes of f due to advection, spatial diffusion, the energy changes due to compressive flows, which represents firstorder Fermi acceleration, and the energy changes due to the motion of “magnetic irregularities”, as originally suggested by Fermi (1949), which is usually called secondorder Fermi acceleration.
The Parker equation is widely used for studying energetic particle transport, and the most popular mechanism of cosmic ray acceleration, that is, diffusive shock acceleration (Bell 1978; Lee & Fisk 1982; Drury 1983; Jones & Ellison 1991), can be obtained from it. The Parker equation is also used to study the acceleration of cosmic rays at supernova remnant shocks (Lagage & Cesarsky 1983a), the acceleration and transport of solar energetic particles at shocks driven by coronal mass ejections (Zank et al. 2000), and, with some variants, the acceleration and transport of energetic particles at the solar wind termination shock (Zank et al. 2015).
In recent years, however, it has become clear that anomalous transport regimes can be found that are different from normal diffusion. In these cases, the mean square deviation of particle positions grows as , with α> 1 corresponding to superdiffusion and α< 1 corresponding to subdiffusion (e.g., Metzler & Klafter 2000, 2004). For astrophysical plasmas in the presence of magnetic turbulence, transport both parallel and perpendicular to the average magnetic field B_{0} can be found to be either subdiffusive or superdiffusive. Anomalous transport corresponding to perpendicular subdiffusion was reported by Qin et al. (2002a) in the case of slab turbulence (when the magnetic turbulence wavevectors are parallel to B_{0}), while parallel transport was found to be normal. Such a perpendicular subdiffusion is the result of particles tracing back the field lines after pitchangle scattering, and is related to the socalled compound diffusion (Webb et al. 2006; Shalchi & Kourakis 2007). The latter depends on the rate of exponential separation of close field lines, which is lower for slab anisotropy and higher for 2D anisotropy (when the magnetic turbulence wavevectors are in the plane perpendicular to B_{0}, Zimbardo et al. 2009, 2012; Bitane et al. 2010). For a composite model of slab plus twodimensional (2D) turbulence with 80% of fluctuation energy in the 2D spectrum, Qin et al. (2002b) have indeed also recovered diffusion perpendicular to B_{0}. The effect of turbulence anisotropy on particle transport has been studied by Zimbardo et al. (2006) with a threedimensional (3D) anisotropic turbulence model, and they found that for quasislab turbulence, transport is anomalous and corresponds to parallel superdiffusion with α ≃ 1.2, while it corresponds to perpendicular subdiffusion with α ≃ 0.8. In the isotropic case, parallel transport is superdiffusive with α ≃ 1.4, while perpendicular transport is normal. In the quasi2D case, both parallel and perpendicular transport are found to be normal, confirming the results of Qin et al. (2002b). These results emphasize the importance of magnetic turbulence anisotropy. Parallel superdiffusion with α ≃ 1.3 and perpendicular subdiffusion with α ≃ 0.75 were also found by Shalchi & Kourakis (2007) by injecting test particles in a composite turbulence model with 20% slab turbulence and 80% 2D turbulence, thus expanding the range of cases when superdiffusion can be found. Furthermore, Pommois et al. (2007) have shown that parallel superdiffusion with α ≃ 1.2–2.0 and perpendicular subdiffusion with α ≃ 0.7–0.8 can be obtained, and it depends on parameters such as the turbulence level, the turbulence anisotropy, and the ratio between the particle gyroradius and the turbulence correlation lengths. An interesting study by Tautz (2010) shows that perpendicular subdiffusion is the case for magnetostatic slab turbulence, while the inclusion of timedependent electric and magnetic fields leads to parallel superdiffusion as well as to particle energization. A recent numerical study by Pucci et al. (2016), who used a very extended 3D isotropic spectrum, found perpendicular subdiffusion for very long times, while parallel superdiffusion was found for intermediate times corresponding to about 10^{3} gyroperiods, depending on the turbulence level and intermittency. On the other hand, numerical simulations by several authors (Beresnyak 2013; Xu & Yan 2013; Lazarian & Yan 2014; Servidio et al. 2016) have found perpendicular superdiffusion for intermediate times corresponding to particle displacements smaller than the turbulence correlation length, which is due to the effect of Richardson (1926) superdiffusion in the case of a magnetized plasma. This effect holds for scales up to the turbulence injection scale.
Thus, the transport properties strongly depend on the turbulence features and the adopted numerical representation. The implications of superdiffusion for particle acceleration and transport at shocks have recently been considered by Perri & Zimbardo (2012b); Zimbardo & Perri (2013); Lazarian & Yan (2014), and these works have attracted considerable theoretical interest (Uchaikin et al. 2015; Rocca et al. 2015, 2016; Saenko 2016).
Nondiffusive processes can be described by several tools, including Lévy walks, Lévy flights, and fractional derivatives (e.g., Perrone et al. 2013). Anomalous diffusion is due to the presence of a scalefree powerlaw distribution for the free path lengths and/or the trapping times, and the way these powerlaw distributions correspond to fractional derivatives is shown, for example, in Metzler & Klafter (2000, 2004), delCastilloNegrete et al. (2004), Bian & Browning (2008). A few studies have investigated the changes in acceleration and transport properties of energetic particles that are due to fractional derivatives. For instance, Milovanov & Zelenyi (2001) introduced a transport equation in velocity space where the stochastic acceleration term in Eq. (1) is written in terms of fractional derivatives; the solution of this equation leads to energy distributions with powerlaw tails. A similar approach was considered by Bian & Browning (2008), where fractional derivatives in velocity space were introduced to take the fragmented electric fields into account that are due to turbulent reconnection. These correspond to highly localized and intense electric fluctuations with a magnitude that is distributed according to a power law.
Several spacecraft observations have shown that the energetic particle profile upstream of interplanetary shocks is characterized by a powerlaw decay rather than an exponential decay, implying that transport can be superdiffusive (Perri & Zimbardo 2007, 2008, 2009a,b; Sugiyama & Shiota 2011). Recently, this powerlaw profile was also found for relativistic electrons upstream of supernova remnant shocks (Perri et al. 2016). The finding of superdiffusion in space suggests that the Parker equation could be extended by introducing fractional derivatives for the first term on the righthand side of Eq. (1) as well, representing spatial diffusion. This approach was considered for cosmic rays by Uchaikin (2010), but in a transport equation that did not include the advection term. Nonetheless, advection is crucial for shock physics, and an appropriate fractional diffusionadvection equation for energetic particle transport was written by Litvinenko & Effenberger (2014), who found analytical solutions in terms of Fourier transforms in the timedependent case.
Here we pursue the study of the Parker equation in the case of fractional spatial derivatives of order μ< 2, representing superdiffusion, and we concentrate on the steadystate case. In Sect. 2 we briefly summarize some wellknown predictions of the standard Parker equation, such as the upstream exponential particle decay and the acceleration time. In Sect. 3 the fractional Parker equation is used to estimate the acceleration time for particles accelerated at a planar shock in the case of superdiffusive propagation. In Sect. 4 we address the fractional Parker equation more quantitatively, introducing the left and right Caputo fractional derivatives and investigating the extent to which the advection term can be balanced by the fractional diffusion term. We show, for the first time in the astrophysical literature, that we can obtain a spatial profile for the steadystate particle distribution function in the form of MittagLeffler functions. These functions reduce to stretched exponentials close to the shock, and exhibit powerlaw tails in space far from the shock, rather than an exponential decay, thus recovering the original result of Perri & Zimbardo (2007, 2008). In Sect. 5 we give the conclusions.
2. Standard Parker equation in the shock frame
We consider the Parker equation (1) in the case of an infinite planar shock in the frame of reference at rest with the shock, and we assume that the particle source at the shock, that is, at x = 0, is modeled as Q = Φ_{0}δ(x). For shock acceleration, firstorder Fermi acceleration is thought to be much faster than secondorder Fermi acceleration, and the term p^{2}(∂/∂p)(p^{2}D_{pp}∂f/∂p) can be neglected (e.g., Drury 1983). The constant flow speed is U_{1} upstream of the shock, x< 0, and U_{2} downstream of the shock, x> 0, with both U_{1}, U_{2}> 0 (e.g., Jones & Ellison 1991). Then the divergence of the fluid velocity, ∂U/∂x, is different from zero only at the shock, in agreement with the fact that shock acceleration is due to shock compression (Lee & Fisk 1982). Therefore, we are left with (Bell 1978; Drury 1983; Jones & Ellison 1991) (2)that is, only advection and diffusion processes play a role. We search for stationary solutions of this equation: this means that the shock time evolution is slower than both the cosmic ray acceleration timescale and the energy loss processes. Setting ∂f/∂t = 0 yields (3)For simplicity, we consider a spatially independent diffusion coefficient k (e.g., Lee & Fisk 1982); in this context, we note that the data analysis in Perri & Zimbardo (2012b) has shown that the magnetic fluctuation level at the particle resonant scales, which determines the pitchangle scattering rate and the diffusion coefficient, is nearly constant both upstream and downstream of the shock. Integrating with respect to x, we obtain the flux of particles with fixed momentum p, and the solution of the obtained linear firstorder differential equation can be written as (4)where the subscript 1 (2) indicates upstream (downstream), and the boundedness of f requires that f_{2}(x,p) = const downstream (see also Giacalone 2012). Upstream of the shock, for x< 0, we have the usual exponential decay of the intensity of energetic particles, and the scale of the decay is Δx_{1} ~ k_{1}/U_{1}. The time needed for a particle of velocity v to explore the upstream region can be obtained by dividing the number N of particles of momentum p per unit surface in a sheath of thickness Δx_{1}, that is, f(p)Δx_{1}, by the density flux of particles going from upstream to downstream, dN/ dt = vf(p) / 4 for an isotropic distribution, so that we obtain the oneside cycle time as Δt_{1} = 4Δx_{1}/v. Considering a shockaccelerated particle that cycles back and forth, a cycle time (5)can be obtained (Drury 1983; Lagage & Cesarsky 1983a).
3. Fractional Parker equation: dimensional analysis
Now we consider superdiffusion: this is the result of a random walk that is composed of free paths whose probability distribution of the freepath length ℓ is a power law with diverging variance (which is obtained for μ ≤ 2); indeed, Lévy walks are characterized by a freepath probability distribution given by ψ(ℓ,t) = A  ℓ  ^{− (μ + 1)}δ(  ℓ  − vt) for  ℓ  >ℓ_{0}, while ψ(ℓ,t) has a smooth, nonsingular behavior for  ℓ  <ℓ_{0} (Zimbardo & Perri 2013). In Lévy walks, a coupling between the freepath length and the travel time is assumed, as expressed by the δ(  ℓ  − vt), while in Lévy flights, such a coupling is absent (Zumofen & Klafter 1993). Finitemass particles are more appropriately described by Lévy walks (Zimbardo & Perri 2013; Perrone et al. 2013), while Lévy flights can be considered as an approximate model for Lévy walks (Zumofen & Klafter 1993; Metzler & Klafter 2004; Perri et al. 2015). We can give a physical explanation of superdiffusive transport in terms of three physical effects: 1. the superdiffusive transport discovered by Richardson (1926), which is due to the presence in turbulence of a hierarchy of vortices such that largescale vortices move faster; this effect was considered for particle transport in magnetic turbulence by Beresnyak (2013), Xu & Yan (2013), Lazarian & Yan (2014), for instance, and is a rather general property of turbulent transport on scales smaller than the turbulence injection scale (e.g., Servidio et al. 2016). In this connection, we note that in astrophysical systems the turbulence injection scale can be very large, for instance, on the order of 50−100 parsec in the galactic disk (e.g., Amato 2014), or in other words, much larger than the typical supernova remnant shock size; 2. the finding by Perri & Zimbardo (2012b) that the magnetic variances of turbulence at the particle resonant scale, observed by Ulysses spacecraft, have a broad powerlaw distribution: these variances are the basic ingredient that enters the quasilinear pitchangle scattering rate (e.g., Amato 2014), and the observation of a broad distribution of magnetic variances implies that the scattering times also have a broad powerlaw distribution; 3. the finding in the test particle numerical simulations by Pucci et al. (2016) that the parallel velocity inversion times in the particle random walk in the presence of an isotropic threedimensional magnetic turbulence also have a long powerlaw distribution. The second and third findings are the ingredients of a particle Lévy walk, which leads to superdiffusion, because a powerlaw distribution of pitchangle scattering times gives rise to a powerlaw distribution of freepath lengths ℓ.
We furthermore consider that in a collisionless plasma (as is the case for cosmic rays during their transport and acceleration) many nonlocal effects are present: this is found in the long correlation scales of turbulence, in the measurements of solar energetic particles that propagate nearly scatterfree in the solar wind (Lin 1974; Reames 1999, and many others), and even in the measurements of the terminationshock energetic particles, which suggest mean free paths on the order of 30 AU (Giacalone 2013). These observations imply that the usual Fick law, J = − k∇f, with the particle flux J related to the local gradient ∇f of the particle density, should be extended to include the effects of nonlocal faraway particle density gradients.
When describing the transport process by a differential equation in phase space, the neglect of the derivative with respect to the pitch angle requires that the particle distribution function f(x,p,t) be nearly isotropic (e.g., Arthur & le Roux 2013; Giacalone 2013). This requires that enough scattering be present, and we point out that the Lévy walk model allows for a nonnegligible probability of very long free paths, but also a very high probability of very short free paths (Perri & Zimbardo 2015), which means frequent pitchangle scattering. We therefore consider the assumption of near isotropy still valid in the case of superdiffusion. For superdiffusive transport, the derivatives of the diffusion term in Eq. (2) are taken to be fractional with 1 <μ< 2, see below and the Appendix. Fractional derivatives are integrodifferential operators that reduce to ordinary derivatives when μ is integer, while because of the integral form, they are an appropriate tool to study nonlocal phenomena, longrange correlations, and scalefree transport (Podlubny 1999; Metzler & Klafter 2000; Zaslavsky 2002).
Then, the fractional Parker equation reads (Litvinenko & Effenberger 2014) (6)where k_{μ} has dimensions (length^{μ}/time) and is assumed to be independent of x, as for k(x,p) in Eq. (3). We note that this equation was written by Litvinenko & Effenberger (2014) using the Riesz fractional derivative, see the appendix; however, in this section we are mainly concerned with dimensional analysis, so that a discussion on the different forms of fractional derivatives is given later. The fractional derivative term corresponds to superdiffusion for 1 <μ< 2, with and α = 3 − μ. We note that in the case of superdiffusion, , which appears in the mean square displacement, is different from k_{μ}, which appears in the transport equation (see Perri et al. 2015, for details).
We search for stationary solutions of the fractional Parker equation. Either upstream or downstream of the shock (i.e., for x ≠ 0), we have (7)From dimensional analysis we have (8)which yields (9)that is, the spatial scale Δx of the region where advection and diffusion “balance” each other. This expression generalizes the normal diffusive scale of decay given by Eq. (4). The scale Δx is the same as the “diffusion length” considered in Litvinenko & Effenberger (2014). Furthermore, considering a probabilistic description based on Lévy walks, we can obtain (Perri et al. 2015; Zimbardo et al. 2015) (10)where Γ( − μ) is the Euler gamma function and ℓ_{0} is the powerlaw break length introduced above.
Now, the typical time for a particle of speed v, with the velocity distributed isotropically, to explore the thickness Δx is Δt ~ 4Δx/v, and the following expression is found for the cycle time of a particle going back and forth over the shock: (11)Therefore, the cycle time for superdiffusive shock acceleration scales with the anomalous diffusion constant k_{μ} and with the plasma velocities U_{1},U_{2} in a very different way than in the case of normal diffusion, see Eq. (5). The energeticparticle acceleration time is related to the cycle time by (Drury 1983; Lagage & Cesarsky 1983a), so we obtain (12)This new scaling of t_{acc} can have a direct influence on the maximum reachable cosmicray energy in supernova remnant shocks (Lagage & Cesarsky 1983b). As shown by Perri & Zimbardo (2015) for energetic particles accelerated at heliospheric shocks, considerably shorter acceleration times are indeed found in the case of superdiffusion than in the case of normal diffusion.
4. Steadystate solutions of the fractional Parker equation
We start from Eq. (6), still with the shock at x = 0 and with 1 <μ< 2 for superdiffusion. We search for stationary solutions, ∂f/∂t = 0, and we integrate from − ∞ to x, so that we have (13)for x< 0 (upstream) and (14)for x ≥ 0 (downstream). We note that on the lefthand side we have the sum of the advective flux and the diffusive flux, and that the second term on the lefthand side of the above equations corresponds to the fractional Fick law, J = − k_{μ}d^{μ − 1}f/ d  x  ^{μ − 1} , where J is the diffusive flux, considered for instance by Zanette (1998), Chaves (1998), Paradisi et al. (2001), Calvo et al. (2007). For the Caputo fractional derivative, adopted below, the derivative of a constant is zero (this does not hold for the more common RiemannLiouville fractional derivative). Therefore we can obtain the solution of the above nonhomogeneous equation for the downstream side as the sum of a constant function plus the solution of the associated homogeneous equation: (15)where we changed the symbol of the partial derivative to that of ordinary derivative, to emphasize that we are considering an ordinary differential equation (ODE; of course, the same applies to Eq. (13)).
We introduce the dimensionless variable z = x/ (k_{μ}/U)^{1 / (μ − 1)}, and we note that both k_{μ} and U are strictly positive constants. The change of variable gives (16)
4.1. Caputo fractional derivatives
We now specify the type of fractional derivative operator to be used. The completely symmetric Riesz derivative, given in the appendix, is appropriate for equations holding in a homogeneous infinite spatial domain. However, when considering spatial domains with given boundary conditions, the Caputo fractional derivative in space should be used (delCastilloNegrete et al. 2004). The presence of a planar shock at x = 0, which is the source of accelerated particles is precisely such a boundary. The left Caputo derivative on the order of β for 0 <β< 1 is defined as (Caputo 1967; Gorenflo & Mainardi 1997; delCastilloNegrete et al. 2004) (17)where df/ dx represents the first (integer) derivative of f(x). Integration intervals different from [ 0,x ] as well as β> 1 can also be used, but here we refrain from generality in favor of immediacy. The form of the Caputo derivative shows that it represents the average of the gradient df/ dx from x = 0 to the current position in x; each contribution to this average is multiplied by a weight inversely proportional to the distance, (x − x′)^{− β}, thus giving more weight to points nearby, but also taking into account distant points. This corresponds to the nonlocal character of fractional derivatives. This also embodies the properties of Lévy walks because the Caputo derivative allows the gradient of particle density far away at x′ to contribute to the flux in the given point x. In particular, the exponent β is related to the exponent of the freepath probability ψ(ℓ,t) by β = μ − 1, see Eq. (16).
When we apply the left Caputo derivative to power functions, f(x) = x^{n}, we easily find that for n ≥ 1 and x> 0(18)Of course, in the considered planar shock configuration, x> 0 corresponds to the downstream side.
We consider Eq. (13) upstream of shock, that is, for x< 0. This means that we have to introduce the nonlocal differential operator appropriate for x< 0: (19)We make the above choice for the following reason: being upstream of the shock, our form of the fractional derivative shows that it represents the average of the gradient df/ dx, weighted by (x′ − x)^{− β}, from the current position at x< 0 to the shock at x = 0. If we consider the fractional Fick law, J = − k_{μ}d^{μ − 1}f/ d  x  ^{μ − 1}, we have that if df/ dx is positive, the fractional derivative d^{β}f/d  x  ^{β} will also be positive; hence, the flux J will be negative, that is, toward upstream, as physical intuition suggests. We note that right Caputo derivatives are also defined, which are similar to Eq. (19), but with a different sign (e.g., Anastassiou 2009). Thus, we can write the above expression in terms of the right Caputo derivative as (20)Applying the right Caputo derivative to a power function of integer order n for x< 0, we have (21)Keeping in mind that the final power exponent will be n − β, that is, a real number, we express x in terms of its absolute value, x = −  x , and make a change of integration variable setting ξ = − x′/  x . Then we obtain (22)The last integral corresponds to the definition of Euler beta function, (23)As is known, the B function is related to the Γ function by B(z,w) = Γ(z)Γ(w) / Γ(z + w), so finally we obtain (24)This holds for x< 0, n ≥ 1 and 0 <β< 1, and completes Eq. (18), which holds for x> 0. We can also write the right Caputo derivative of the absolute value  x , for x< 0, as (25)This relation is used in the following subsection.
Solutions of normal and fractional transport equations.
4.2. MittagLeffler functions
The solution of fractional ODEs is given by the MittagLeffler functions, which are defined as (MittagLeffler 1903; Gorenflo & Mainardi Gorenflo & Mainardi 1997; Zaslavsky 2002; delCastilloNegrete et al. 2004; Mainardi 2014) (26)where β> 0 and z varies in the complex plane (although we 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 MittagLeffler functions depending on two and three indices are also used, but these functions are not needed here.
We note that in the case of normal transport, the steadystate condition implies a linear ODE, for instance, Eq. (3), and the exponential function has a central role in finding the solution. In the case of fractional ordinary differential equations (fODE), such a central role is played by MittagLeffler functions instead of the exponential function (Mainardi 2014). The relative roles of the exponential and the MittagLeffler functions are analogous to the role of the Gaussian propagator for normal diffusion and to that of Lévy distributions for anomalous diffusion in the timedependent case. These roles are summarized in Table 1.
Fig. 1 Plot of MittagLeffler functions E_{β}( −  z  ^{β}) for several β values (see legend) in loglin axes. Note that β = 1 corresponds to the standard exponential. A similar plot in linear axes is given by Mainardi (2014). 
4.3. Upstream behavior
For the region upstream of the shock we have x< 0, which implies z< 0, as well. Setting β = μ − 1 and using Eqs. (19), (20), we can write Eq. (16) as (27)Then we use  z  and for the negative real semiaxis we write the solution as (Gorenflo & Mainardi 1997; Mainardi 2014) (28)The fact that E_{β}( −  z  ^{β}) is the solution of Eq. (27) can be verified by direct substitution: (29)where we set m = n − 1. Plots of E_{β}( −  z  ^{β}) for a few β values are given in Figs. 1 and 2, where the comparison with the exponential shows that E_{β}( −  z  ^{β}) is steeper close to  z  = 0 and shallower for large  z . In particular, E_{β}( −  z  ^{β}) corresponds to a stretched exponential for a small argument  z  → 0, as can be seen from the first terms in the powerseries expansion: (30)Conversely, for a large argument  z  → ∞, we have a powerlaw decay to first order (Metzler & Klafter 2000; Mainardi 2014), (31)Returning to dimensional variables and recalling that β = μ − 1 < 1, we immediately recover the “longdistance” behavior originally derived by Perri & Zimbardo (2007, 2008) for the upstream energetic particle profile: (32)that is, a powerlaw decay rather than an exponential decay as given by Eq. (4), with the powerlaw exponent related to the superdiffusion exponent α by β = 2 − α. This behavior also corresponds to the behavior found by Litvinenko & Effenberger (2014) in the limit of  x  ≪ U_{1}t, as expressed by their Eq. (22). We point out that this powerlaw decay has been confirmed by spacecraft observations of energetic particles at several shock crossings in the heliosphere (e.g., Zimbardo et al. 2012). Conversely, for small  z  , the typical shape of E_{β}( −  z  ^{β}) is very similar to the energetic particle profiles reported by Giacalone (2012, 2015), Perri & Zimbardo (2015) close to the shock. Figures 1 and 2 also show rather clearly that the transition from the stretched exponential behavior, expressed by Eq. (30), to the powerlaw behavior, expressed by Eqs. (31) and (32), occurs at  z  ≃ 1, as expected. The break in the powerlaw behavior is therefore found at (33)in full agreement with Eq. (9) and with the methods developed by Perri & Zimbardo (2015) to determine the particle acceleration time at shocks.
Fig. 2 Plot of MittagLeffler functions E_{β}( −  z  ^{β}) in loglog axes for several β values (see legend). 
We furthermore note that Mainardi (2014) recently conjectured that two simple rational approximations (an upper and a lower bound) to the MittagLeffler funtion E_{β}( −  z  ^{β}) hold for every value of  z  and for 0 <β< 1: (34)This conjecture is supported by rather extensive numerical evaluations of the three functions above, and relevant graphs can be found in Mainardi (2014). It is worth noting that the above rational approximations for the MittagLeffler functions can be very useful for a quick comparison with spacecraft data.
4.4. Downstream behavior
For x> 0 we can write Eq. (16) with the left Caputo derivative, , and the solution can be conveniently expressed as (35)as can be verified by direct substitution of E_{β}(z^{β}) in Eq. (16) and using Eqs. (17), (18). Therefore, the complete solution to Eq. (14) is (36)The first term on the righthand side is a purely growing function of x, and it diverges as x goes far downstream. Since the solution for x → ∞ must be equal to Φ_{0}/U_{2}, we have to assume that the contribution from the diffusive term of Eq. (14) is zero, as for the case of normal diffusion. As argued by Drury (1983), the reason is that it is not possible to balance the diffusive flux against the advective flux in the downstream because both fluxes are in the same direction. Then for x> 0 we are left with f_{2}(x,p) = Φ_{0}/U_{2}.
4.5. Comparison between the Riesz and Caputo fractional equation approach to anomalous diffusion
Since several different forms of fractional derivatives are used in the literature, it is interesting to compare them. The relationship between the Caputo derivative, Eq. (17), and the RiemannLiouville derivative, Eq. (A.1), for β< 1 is easily found to be (see, e.g., Gorenflo & Mainardi 1997; Mainardi 2014) (37)so that they basically differ by an additional powerlaw term. In particular, if f(x) is a power function of positive exponent, they give the same result. These relationships may be the reason why here we can recover results obtained with the RiemannLiouville and the Riesz derivatives. For instance, Litvinenko & Effenberger (2014) solved the timedependent fractional Parker equation written with the Riesz derivative, and found an asymptotic solution for 0 <  x  ≪ Ut (their Eq. (22)) of the form (38)where A is a constant, corresponding to our Eq. (32). Furthermore, using the Riesz derivative, Stern et al. (2014) found solutions of a space fractional diffusion equation in terms of Fox Hfunctions and in terms of Fourier series representation, see their Eq. (2.20). In this connection, it is useful to recall that the MittagLeffler functions can be obtained from Fox Hfunction by a relatively simple choice of the parameters of these functions, see Eq. (B.19) in Metzler & Klafter (2000) and Eq. (A.9) in Qi & Jiang (2011). See also Eqs. (2.17) and (2.19) of Webb et al. (2006). Nevertheless, further study is needed to fully understand the relation between the timedependent and the steadystate solutions.
It is interesting to note, regarding the type of fractional derivative to be used in the fractional Fick law, that Bian et al. (2017) have very recently advocated the use of a nonlocal form for the particle diffusive flux in the context of electron transport in the solar corona (see their Eq. (13)): in this case, the flux is proportional to the spatial convolution of the density gradient, similar to our Eq. (17), but with a different spatial kernel. The form proposed by Bian et al. (2017) supports our assumptions that the Caputo derivative, with the density gradient df/ dx inside the integral, is most appropriate to express the fractional Fick law.
5. Discussion and conclusions
We have studied the extension of the Parker transport equation for cosmic rays to the case of superdiffusive transport, introducing fractional derivatives in the spatial diffusion term. We considered the case of steadystate solutions upstream and downstream of a planar shock, and obtained by dimensional analysis an estimate of the particle acceleration time at shocks in the case of superdiffusion. To study the fractional diffusionadvection equation in the presence of a shock, we proposed to use the left and right Caputo fractional derivatives in space. We have shown for the first time that an analytical solution of the steadystate fractional Parker equation is given by the MittagLeffler functions. These functions are a generalization of the standard exponential function. When considering the case of a planar shock with the shock as the source of energetic particles, the MittagLeffler functions correspond, for the energetic particle intensity, to a stretched exponential close to the shock, and to a powerlaw profile far upstream of the shock. The latter result is in agreement with the results obtained from a probabilistic approach to superdiffusion by Perri & Zimbardo (2007, 2008) and from a timeasymptotic solution obtained by Litvinenko & Effenberger (2014).
With these solutions, we have an analytical expression for the energetic particle density profile at every distance from the shock, which allows for a more detailed comparison with observations. This will be useful for the study of superdiffusive transport both at heliospheric shocks and at supernova remnant shocks (Perri et al. 2016); in particular, the possibility to fit the entire upstream profile with the MittagLeffler functions will allow a more precise determination of the break in the powerlaw profile, see Eqs. (9) and (33), and of the corresponding acceleration time for energetic particles (Perri et al. 2015; Perri & Zimbardo 2015). It is interesting to note that spacecraft data in the heliosphere show that the energetic particle profile close upstream of shocks is very steep (e.g., Giacalone 2012, 2015): while this corresponds to the stretched exponential limit of MittagLeffler functions, the data also show that in order to make a good fit, a much higher temporal resolution for energetic particles than currently available is needed. In this connection, the planned THOR mission will be able to resolve the energetic particle profile close upstream with a resolution of a few seconds, and it will be possible to fit the full MittagLeffler function with great accuracy.
A wellknown problem with anomalous transport is that fractional derivatives correspond to Lévy flights and not Lévy walks, the latter implying a spacetime delta coupling. This spacetime coupling translates into the fact that the Lévy walk propagator P(x,t) is similar to a Lévy distribution, but goes to zero for  x  >vt (see, e.g., Blumen et al. 1990; Zumofen & Klafter 1993; Perri et al. 2015). It is interesting, however, to note that in steadystate we can assume t → ∞, and thus the Lévy walk propagator fully coincides with a Lévy distribution. Hence, the solution of a fractional diffusionadvection equation in steadystate is entirely appropriate to describe the behavior of Lévy walks.
Acknowledgments
The work by G. Z and S. P. has been supported by the Agenzia Spaziale Italiana under the contract ASIINAF 2015039R.O “Missione M4 di ESA: Partecipazione Italiana alla fase di assessment della missione THOR”. F.E. acknowledges partial support from NASA grants NNX14AG03G and NNX17AK25G, and the International Space Science Institute, Bern.
References
 Amato, E. 2014, Int. J. Mod. Phys. D, 23, 1430013 [Google Scholar]
 Anastassiou, G. A. 2009, Chaos, Solitons and Fractals, 42, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Arthur, A. D., & le Roux, J. A. 2013, ApJ, 772, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 1978, MNRAS, 182, 14 [Google Scholar]
 Beresnyak, A. 2013, ApJ, 767, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Bian N. H., & Browning P. K. 2008, ApJ, 687, L111 [NASA ADS] [CrossRef] [Google Scholar]
 Bian, N. H., Emslie, A. G., & Kontar, E. P. 2017, ApJ, 835, 262 [NASA ADS] [CrossRef] [Google Scholar]
 Bitane, R., Zimbardo, G., Veltri, P. 2010, ApJ, 719, 1912 [NASA ADS] [CrossRef] [Google Scholar]
 Blumen, A., Klafter, J., & Zumofen, G. 1990, Europhys. Lett., 13, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Calvo, I., Sanchez, R., Carreras, B. A., & van Milligen, B. Ph. 2007, Phys. Rev. Lett., 99, 230603 [NASA ADS] [CrossRef] [Google Scholar]
 Caputo, M. 1967, Geophys. J., 13, 529 [Google Scholar]
 Chaves, A. S. 1998, Phys. Lett. A, 239, 13 [NASA ADS] [CrossRef] [Google Scholar]
 delCastilloNegrete, D., Carreras, B. A., & Lynch, V. E. 2004, Phys. Plasmas, 11, 3854 [NASA ADS] [CrossRef] [Google Scholar]
 Drury L. O’C. 1983, Rep. Progr. Phys., 46, 973 [Google Scholar]
 Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J. 2012, ApJ, 761, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J. 2013, Space Sci. Rev., 176, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J. 2015, ApJ, 799, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Gorenflo, R., & Mainardi, F. 1997, in Fractals and Fractional Calculus in Continuum Mechanics eds. A. Carpinteri and F. Mainardi, (Wien: Springer Verlag), 223 [Google Scholar]
 Gorenflo, R., De Fabritiis, G., & Mainardi, F. 1999, Physica A, 269, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, F. C., & Ellison, D. C. 1991, Space Sci. Rev., 58, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Lagage, P.O., & Cesarsky, C.J. 1983a, A&A, 118, 223 [NASA ADS] [Google Scholar]
 Lagage, P.O., & Cesarsky, C.J. 1983b, A&A, 125, 249 [NASA ADS] [Google Scholar]
 Lazarian, A., & Yan, H. 2014, ApJ, 784, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, M. A., & Fisk, L. A. 1982, Space Sci. Rev., 32, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, R. P. 1974, Space Sci. Rev., 16, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Litvinenko, Y., & Effenberger, F. 2014, ApJ, 796, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Mainardi, F. 2014, Discrete and Continuos Dynamical System – Series B, 19, 2267 [CrossRef] [Google Scholar]
 Metzler, R., & Klafter, J. 2000, Phys. Rep., 339,1 [NASA ADS] [CrossRef] [Google Scholar]
 Metzler, R., & Klafter, J. 2004, J. Phys. A: Math. Gen., 37, R161 [NASA ADS] [CrossRef] [Google Scholar]
 Milovanov, A. V., & Zelenyi, L. M. 2001, Phys. Rev. E, 64, 052101 [NASA ADS] [CrossRef] [Google Scholar]
 MittagLeffler, M. G., 1903, C. R. Acad. Sci. Paris, 137, 554 [Google Scholar]
 Moraal, H. 2013, Space Sci Rev., 176, 299 [Google Scholar]
 Paradisi, P., Cesari, R., Mainardi, F., & Tampieri, F. 2001, Phys. A: Stat. Mech. Applications, 293, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1965, Planet. Space Sci., 13, 9 [Google Scholar]
 Perri, S., & Zimbardo, G. 2007, ApJ, 671, L177 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, S., & Zimbardo, G. 2008, J. Geophys. Res., 113, A03107 [NASA ADS] [Google Scholar]
 Perri, S., & Zimbardo, G. 2009a, ApJ, 693, L118 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, S., & Zimbardo, G. 2009b, Adv. Space Res., 44, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, S., & Zimbardo, G. 2012a, ApJ, 750, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, S., & Zimbardo, G. 2012b, ApJ, 754, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, S., & Zimbardo, G. 2015, ApJ, 815, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Perri, S., Zimbardo, G., Effenberger, F., & Fichtner, H. 2015, A&A, 578, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perri, S., Amato, E., & Zimbardo, G. 2016, A&A, 596, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perrone, D., Dendy, R O., Furno, I., et al. 2013, Space Sci. Rev., 10.1007 [Google Scholar]
 Podlubny, I. 1999. Fractional Differential Equations, (San Diego: Academic Press) [Google Scholar]
 Pommois, P., Zimbardo, G., & Veltri, P. 2007, Phys. Plasmas, 14, 012311 [NASA ADS] [CrossRef] [Google Scholar]
 Pucci, F., Malara, F., Perri, S., et al. 2016, MNRAS, 459, 3395 [NASA ADS] [CrossRef] [Google Scholar]
 Qi, H. T., & Jiang, X. Y. 2011, Physica A, 390, 1876 [NASA ADS] [CrossRef] [Google Scholar]
 Qin, G., Matthaeus, W. H., & Bieber, J. W. 2002a, Geophys. Res. Lett., 29, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Qin, G., Matthaeus, W. H., & Bieber, J. W., 2002b, ApJ, 578, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Reames, D. V., 1999, Space Sci. Rev., 90, 413 [Google Scholar]
 Richardson, L. F. 1926, Proc. R. Soc. London A, 110, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Rocca, M. C., Plastino, A. R., Plastino, A., Ferri, G. L., & de Paoli, A. 2015, J. Stat. Phys., 161, 986 [NASA ADS] [CrossRef] [Google Scholar]
 Rocca, M. C., Plastino, A. R., Plastino, A., Ferri, G. L., & de Paoli, A. 2016, Phys. A, 447, 402 [CrossRef] [Google Scholar]
 Saenko, V. V. 2016, Physica A, 444, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Servidio, S., Haynes, C. T., Matthaeus, W. H. et al. 2016, Phys. Rev. Lett., 117, 095101 [NASA ADS] [CrossRef] [Google Scholar]
 Shalchi, A., & Kourakis I. 2007, A&A, 470, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stern, R., Effenberger, F., Fichtner, H., & Schäfer, T. 2014, Fract. Calc. Appl. Anal., 17, 171 [CrossRef] [Google Scholar]
 Sugiyama, T., & Shiota, D. 2011, ApJ, 731, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C. 2010, Plasma Phys. Control. Fusion, 52, 045016 [NASA ADS] [CrossRef] [Google Scholar]
 Uchaikin, V. V. 2010, JETP Lett., 91, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Uchaikin, V. V., Sibatov, R. T., & Byzykchi, A. N. 2015, Bull. Russian Acad. Sci. Physics, 79, 592 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, G. M., Zank, G. P., Kaghashvili, E. Kh., & le Roux, J. A. 2006, ApJ, 651, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, S., & Yan, H. 2013, ApJ, 779, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Zanette, D. 1998, Physica A: Stat. Mech. Applications, 252, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Zank, G. P., Rice, W. K. M., & Wu, C. C. 2000, J. Geophys. Res., 105, 25 [Google Scholar]
 Zank, G., Hunana, P., Mostafavi, P., et al. 2015, ApJ, 814, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Zaslavsky, G. M. 2002, Phys. Rep., 371, 461 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Zimbardo, G., & Perri, S. 2013, ApJ, 778, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Zimbardo, G., Pommois, P., & Veltri, P. 2006, ApJ, 639, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Zimbardo, G., Bitane, R., Pommois, P., & Veltri, P. 2009, Plasma Phys. Control. Fusion, 51, 015005 [NASA ADS] [CrossRef] [Google Scholar]
 Zimbardo, G., Perri, S., Pommois, P., & Veltri, P. 2012, Adv. Space Res. 49, 1633 [NASA ADS] [CrossRef] [Google Scholar]
 Zimbardo, G., Amato, E., Bovet, A., et al. 2015, J. Plasma Phys., 81, 495810601 [CrossRef] [Google Scholar]
 Zumofen, G., & Klafter, J. 1993, Phys. Rev. E, 47, 851 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Fractional derivatives
Fractional derivatives are integrodifferential operators that allow for noninteger order of differentiation μ and that reduce to normal derivatives when μ is integer (Podlubny 1999; Metzler & Klafter 2000; Zaslavsky 2002). Several definitions are possible. In particular, one popular definition is that of RiemannLiouville fractional derivatives, whose left and righthand expressions are defined as (A.1)where the integration domain is to the left of x, and (A.2)where the integration domain is to the right of x; for both expression, m − 1 <μ<m with m the smallest integer larger than μ. We can introduce the fractional Riesz derivative ∂^{μ}/∂  x  ^{μ} with the fractional index μ< 2 as a generalization of the
onedimensional Laplace operator. The symmetric Riesz derivative (Gorenflo et al. 1999; Metzler & Klafter 2000; Zaslavsky 2002) is then given by (A.3)The solution of the timedependent spacefractional diffusion equation can be expressed in terms of the Lévy distribution, whose Fourier transform is given by (A.4)with q the wavenumber (Metzler & Klafter 2000, 2004; Litvinenko & Effenberger 2014; Perri et al. 2015). This Fourier transform is easily obtained by writing the transport equation with the fractional Riesz derivative. The backFourier transform gives a Lévy distribution with powerlaw tails ~  x  ^{− μ − 1}, and this shows the relation between Lévy flights and fractional derivatives (see Zimbardo & Perri 2013; Perri et al. 2015, for more details).
All Tables
All Figures
Fig. 1 Plot of MittagLeffler functions E_{β}( −  z  ^{β}) for several β values (see legend) in loglin axes. Note that β = 1 corresponds to the standard exponential. A similar plot in linear axes is given by Mainardi (2014). 

In the text 
Fig. 2 Plot of MittagLeffler functions E_{β}( −  z  ^{β}) in loglog axes for several β values (see legend). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.