EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 607 (November 2017) A&A, 607 (2017) A7 Full HTML
Free Access
Issue
A&A
Volume 607, November 2017
Article Number A7
Number of page(s) 9
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201731179
Published online 30 October 2017

© ESO, 2017

1. Introduction

The acceleration and transport of energetic particles in space plasmas and in astrophysical environments can be described by the so-called 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 near-isotropy in velocity space. In the case of a one-dimensional 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 omni-directional 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, Dpp 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 first-order Fermi acceleration, and the energy changes due to the motion of “magnetic irregularities”, as originally suggested by Fermi (1949), which is usually called second-order 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 B0 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 B0), while parallel transport was found to be normal. Such a perpendicular subdiffusion is the result of particles tracing back the field lines after pitch-angle scattering, and is related to the so-called 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 B0, Zimbardo et al. 2009, 2012; Bitane et al. 2010). For a composite model of slab plus two-dimensional (2D) turbulence with 80% of fluctuation energy in the 2D spectrum, Qin et al. (2002b) have indeed also recovered diffusion perpendicular to B0. The effect of turbulence anisotropy on particle transport has been studied by Zimbardo et al. (2006) with a three-dimensional (3D) anisotropic turbulence model, and they found that for quasi-slab 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 quasi-2D 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 time-dependent 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 103 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 scale-free power-law distribution for the free path lengths and/or the trapping times, and the way these power-law distributions correspond to fractional derivatives is shown, for example, in Metzler & Klafter (2000, 2004), del-Castillo-Negrete 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 power-law 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 power-law decay rather than an exponential decay, implying that transport can be superdiffusive (Perri & Zimbardo 2007, 2008, 2009a,b; Sugiyama & Shiota 2011). Recently, this power-law 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 right-hand 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 diffusion-advection equation for energetic particle transport was written by Litvinenko & Effenberger (2014), who found analytical solutions in terms of Fourier transforms in the time-dependent 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 steady-state case. In Sect. 2 we briefly summarize some well-known 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 steady-state particle distribution function in the form of Mittag-Leffler functions. These functions reduce to stretched exponentials close to the shock, and exhibit power-law 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, first-order Fermi acceleration is thought to be much faster than second-order Fermi acceleration, and the term p-2(∂/∂p)(p2Dpp∂f/∂p) can be neglected (e.g., Drury 1983). The constant flow speed is U1 upstream of the shock, x< 0, and U2 downstream of the shock, x> 0, with both U1, U2> 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 pitch-angle 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 first-order differential equation can be written as (4)where the subscript 1 (2) indicates upstream (downstream), and the boundedness of f requires that f2(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 Δx1 ~ k1/U1. 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 Δx1, that is, f(p)Δx1, 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 one-side cycle time as Δt1 = 4Δx1/v. Considering a shock-accelerated 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 free-path length ℓ is a power law with diverging variance (which is obtained for μ ≤ 2); indeed, Lévy walks are characterized by a free-path 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 free-path 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). Finite-mass 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 large-scale 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 power-law distribution: these variances are the basic ingredient that enters the quasi-linear pitch-angle 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 power-law 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 three-dimensional magnetic turbulence also have a long power-law distribution. The second and third findings are the ingredients of a particle Lévy walk, which leads to superdiffusion, because a power-law distribution of pitch-angle scattering times gives rise to a power-law distribution of free-path 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 scatter-free in the solar wind (Lin 1974; Reames 1999, and many others), and even in the measurements of the termination-shock 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 far-away 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 non-negligible probability of very long free paths, but also a very high probability of very short free paths (Perri & Zimbardo 2015), which means frequent pitch-angle 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 integro-differential operators that reduce to ordinary derivatives when μ is integer, while because of the integral form, they are an appropriate tool to study nonlocal phenomena, long-range correlations, and scale-free 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 power-law 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 U1,U2 in a very different way than in the case of normal diffusion, see Eq. (5). The energetic-particle acceleration time is related to the cycle time by (Drury 1983; Lagage & Cesarsky 1983a), so we obtain (12)This new scaling of tacc can have a direct influence on the maximum reachable cosmic-ray 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. Steady-state 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 left-hand side we have the sum of the advective flux and the diffusive flux, and that the second term on the left-hand side of the above equations corresponds to the fractional Fick law, J = − kμdμ − 1f/ 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 Riemann-Liouville fractional derivative). Therefore we can obtain the solution of the above non-homogeneous 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 (del-Castillo-Negrete 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; del-Castillo-Negrete 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 free-path probability ψ(ℓ,t) by β = μ − 1, see Eq. (16).

When we apply the left Caputo derivative to power functions, f(x) = xn, 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μ − 1f/ 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.

Table 1

Solutions of normal and fractional transport equations.

4.2. Mittag-Leffler functions

The solution of fractional ODEs is given by the Mittag-Leffler functions, which are defined as (Mittag-Leffler 1903; Gorenflo & Mainardi Gorenflo & Mainardi 1997; Zaslavsky 2002; del-Castillo-Negrete 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 Mittag-Leffler 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 steady-state 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 Mittag-Leffler functions instead of the exponential function (Mainardi 2014). The relative roles of the exponential and the Mittag-Leffler 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 time-dependent case. These roles are summarized in Table 1.

thumbnail Fig. 1

Plot of Mittag-Leffler functions Eβ( − | z | β) for several β values (see legend) in log-lin 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 semi-axis 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 power-series expansion: (30)Conversely, for a large argument | z | → ∞, we have a power-law decay to first order (Metzler & Klafter 2000; Mainardi 2014), (31)Returning to dimensional variables and recalling that β = μ − 1 < 1, we immediately recover the “long-distance” behavior originally derived by Perri & Zimbardo (2007, 2008) for the upstream energetic particle profile: (32)that is, a power-law decay rather than an exponential decay as given by Eq. (4), with the power-law 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 | ≪ U1t, as expressed by their Eq. (22). We point out that this power-law 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 power-law behavior, expressed by Eqs. (31) and (32), occurs at | z | ≃ 1, as expected. The break in the power-law 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.

thumbnail Fig. 2

Plot of Mittag-Leffler functions Eβ( − | z | β) in log-log 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 Mittag-Leffler 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 Mittag-Leffler 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 right-hand 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/U2, 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 f2(x,p) = Φ0/U2.

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 Riemann-Liouville 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 power-law 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 Riemann-Liouville and the Riesz derivatives. For instance, Litvinenko & Effenberger (2014) solved the time-dependent 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 H-functions and in terms of Fourier series representation, see their Eq. (2.20). In this connection, it is useful to recall that the Mittag-Leffler functions can be obtained from Fox H-function 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 time-dependent and the steady-state 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 steady-state 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 diffusion-advection 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 steady-state fractional Parker equation is given by the Mittag-Leffler 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 Mittag-Leffler functions correspond, for the energetic particle intensity, to a stretched exponential close to the shock, and to a power-law 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 time-asymptotic 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 Mittag-Leffler functions will allow a more precise determination of the break in the power-law 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 Mittag-Leffler 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 Mittag-Leffler function with great accuracy.

A well-known problem with anomalous transport is that fractional derivatives correspond to Lévy flights and not Lévy walks, the latter implying a space-time delta coupling. This space-time 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 steady-state we can assume t → ∞, and thus the Lévy walk propagator fully coincides with a Lévy distribution. Hence, the solution of a fractional diffusion-advection equation in steady-state 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 ASI-INAF 2015-039-R.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

  1. Amato, E. 2014, Int. J. Mod. Phys. D, 23, 1430013 [Google Scholar]
  2. Anastassiou, G. A. 2009, Chaos, Solitons and Fractals, 42, 365 [NASA ADS] [CrossRef] [Google Scholar]
  3. Arthur, A. D., & le Roux, J. A. 2013, ApJ, 772, L26 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bell, A. R. 1978, MNRAS, 182, 14 [Google Scholar]
  5. Beresnyak, A. 2013, ApJ, 767, L39 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bian N. H., & Browning P. K. 2008, ApJ, 687, L111 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bian, N. H., Emslie, A. G., & Kontar, E. P. 2017, ApJ, 835, 262 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bitane, R., Zimbardo, G., Veltri, P. 2010, ApJ, 719, 1912 [NASA ADS] [CrossRef] [Google Scholar]
  9. Blumen, A., Klafter, J., & Zumofen, G. 1990, Europhys. Lett., 13, 223 [NASA ADS] [CrossRef] [Google Scholar]
  10. Calvo, I., Sanchez, R., Carreras, B. A., & van Milligen, B. Ph. 2007, Phys. Rev. Lett., 99, 230603 [NASA ADS] [CrossRef] [Google Scholar]
  11. Caputo, M. 1967, Geophys. J., 13, 529 [Google Scholar]
  12. Chaves, A. S. 1998, Phys. Lett. A, 239, 13 [NASA ADS] [CrossRef] [Google Scholar]
  13. del-Castillo-Negrete, D., Carreras, B. A., & Lynch, V. E. 2004, Phys. Plasmas, 11, 3854 [NASA ADS] [CrossRef] [Google Scholar]
  14. Drury L. O’C. 1983, Rep. Progr. Phys., 46, 973 [Google Scholar]
  15. Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] [Google Scholar]
  16. Giacalone, J. 2012, ApJ, 761, 28 [NASA ADS] [CrossRef] [Google Scholar]
  17. Giacalone, J. 2013, Space Sci. Rev., 176, 73 [NASA ADS] [CrossRef] [Google Scholar]
  18. Giacalone, J. 2015, ApJ, 799, 80 [NASA ADS] [CrossRef] [Google Scholar]
  19. 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]
  20. Gorenflo, R., De Fabritiis, G., & Mainardi, F. 1999, Physica A, 269, 79 [NASA ADS] [CrossRef] [Google Scholar]
  21. Jones, F. C., & Ellison, D. C. 1991, Space Sci. Rev., 58, 259 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lagage, P.-O., & Cesarsky, C.-J. 1983a, A&A, 118, 223 [NASA ADS] [Google Scholar]
  23. Lagage, P.-O., & Cesarsky, C.-J. 1983b, A&A, 125, 249 [NASA ADS] [Google Scholar]
  24. Lazarian, A., & Yan, H. 2014, ApJ, 784, 38 [NASA ADS] [CrossRef] [Google Scholar]
  25. Lee, M. A., & Fisk, L. A. 1982, Space Sci. Rev., 32, 205 [NASA ADS] [CrossRef] [Google Scholar]
  26. Lin, R. P. 1974, Space Sci. Rev., 16, 189 [NASA ADS] [CrossRef] [Google Scholar]
  27. Litvinenko, Y., & Effenberger, F. 2014, ApJ, 796, 125 [NASA ADS] [CrossRef] [Google Scholar]
  28. Mainardi, F. 2014, Discrete and Continuos Dynamical System – Series B, 19, 2267 [CrossRef] [Google Scholar]
  29. Metzler, R., & Klafter, J. 2000, Phys. Rep., 339,1 [NASA ADS] [CrossRef] [Google Scholar]
  30. Metzler, R., & Klafter, J. 2004, J. Phys. A: Math. Gen., 37, R161 [NASA ADS] [CrossRef] [Google Scholar]
  31. Milovanov, A. V., & Zelenyi, L. M. 2001, Phys. Rev. E, 64, 052101 [NASA ADS] [CrossRef] [Google Scholar]
  32. Mittag-Leffler, M. G., 1903, C. R. Acad. Sci. Paris, 137, 554 [Google Scholar]
  33. Moraal, H. 2013, Space Sci Rev., 176, 299 [Google Scholar]
  34. Paradisi, P., Cesari, R., Mainardi, F., & Tampieri, F. 2001, Phys. A: Stat. Mech. Applications, 293, 130 [NASA ADS] [CrossRef] [Google Scholar]
  35. Parker, E. N. 1965, Planet. Space Sci., 13, 9 [Google Scholar]
  36. Perri, S., & Zimbardo, G. 2007, ApJ, 671, L177 [NASA ADS] [CrossRef] [Google Scholar]
  37. Perri, S., & Zimbardo, G. 2008, J. Geophys. Res., 113, A03107 [NASA ADS] [Google Scholar]
  38. Perri, S., & Zimbardo, G. 2009a, ApJ, 693, L118 [NASA ADS] [CrossRef] [Google Scholar]
  39. Perri, S., & Zimbardo, G. 2009b, Adv. Space Res., 44, 465 [NASA ADS] [CrossRef] [Google Scholar]
  40. Perri, S., & Zimbardo, G. 2012a, ApJ, 750, 87 [NASA ADS] [CrossRef] [Google Scholar]
  41. Perri, S., & Zimbardo, G. 2012b, ApJ, 754, 8 [NASA ADS] [CrossRef] [Google Scholar]
  42. Perri, S., & Zimbardo, G. 2015, ApJ, 815, 75 [NASA ADS] [CrossRef] [Google Scholar]
  43. Perri, S., Zimbardo, G., Effenberger, F., & Fichtner, H. 2015, A&A, 578, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Perri, S., Amato, E., & Zimbardo, G. 2016, A&A, 596, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Perrone, D., Dendy, R O., Furno, I., et al. 2013, Space Sci. Rev., 10.1007 [Google Scholar]
  46. Podlubny, I. 1999. Fractional Differential Equations, (San Diego: Academic Press) [Google Scholar]
  47. Pommois, P., Zimbardo, G., & Veltri, P. 2007, Phys. Plasmas, 14, 012311 [NASA ADS] [CrossRef] [Google Scholar]
  48. Pucci, F., Malara, F., Perri, S., et al. 2016, MNRAS, 459, 3395 [NASA ADS] [CrossRef] [Google Scholar]
  49. Qi, H. T., & Jiang, X. Y. 2011, Physica A, 390, 1876 [NASA ADS] [CrossRef] [Google Scholar]
  50. Qin, G., Matthaeus, W. H., & Bieber, J. W. 2002a, Geophys. Res. Lett., 29, 1048 [NASA ADS] [CrossRef] [Google Scholar]
  51. Qin, G., Matthaeus, W. H., & Bieber, J. W., 2002b, ApJ, 578, L117 [NASA ADS] [CrossRef] [Google Scholar]
  52. Reames, D. V., 1999, Space Sci. Rev., 90, 413 [Google Scholar]
  53. Richardson, L. F. 1926, Proc. R. Soc. London A, 110, 709 [NASA ADS] [CrossRef] [Google Scholar]
  54. 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]
  55. Rocca, M. C., Plastino, A. R., Plastino, A., Ferri, G. L., & de Paoli, A. 2016, Phys. A, 447, 402 [CrossRef] [Google Scholar]
  56. Saenko, V. V. 2016, Physica A, 444, 765 [NASA ADS] [CrossRef] [Google Scholar]
  57. Servidio, S., Haynes, C. T., Matthaeus, W. H. et al. 2016, Phys. Rev. Lett., 117, 095101 [NASA ADS] [CrossRef] [Google Scholar]
  58. Shalchi, A., & Kourakis I. 2007, A&A, 470, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Stern, R., Effenberger, F., Fichtner, H., & Schäfer, T. 2014, Fract. Calc. Appl. Anal., 17, 171 [CrossRef] [Google Scholar]
  60. Sugiyama, T., & Shiota, D. 2011, ApJ, 731, L34 [NASA ADS] [CrossRef] [Google Scholar]
  61. Tautz, R. C. 2010, Plasma Phys. Control. Fusion, 52, 045016 [NASA ADS] [CrossRef] [Google Scholar]
  62. Uchaikin, V. V. 2010, JETP Lett., 91, 105 [NASA ADS] [CrossRef] [Google Scholar]
  63. Uchaikin, V. V., Sibatov, R. T., & Byzykchi, A. N. 2015, Bull. Russian Acad. Sci. Physics, 79, 592 [NASA ADS] [CrossRef] [Google Scholar]
  64. Webb, G. M., Zank, G. P., Kaghashvili, E. Kh., & le Roux, J. A. 2006, ApJ, 651, 211 [NASA ADS] [CrossRef] [Google Scholar]
  65. Xu, S., & Yan, H. 2013, ApJ, 779, 140 [NASA ADS] [CrossRef] [Google Scholar]
  66. Zanette, D. 1998, Physica A: Stat. Mech. Applications, 252, 159 [NASA ADS] [CrossRef] [Google Scholar]
  67. Zank, G. P., Rice, W. K. M., & Wu, C. C. 2000, J. Geophys. Res., 105, 25 [Google Scholar]
  68. Zank, G., Hunana, P., Mostafavi, P., et al. 2015, ApJ, 814, 137 [NASA ADS] [CrossRef] [Google Scholar]
  69. Zaslavsky, G. M. 2002, Phys. Rep., 371, 461 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  70. Zimbardo, G., & Perri, S. 2013, ApJ, 778, 35 [NASA ADS] [CrossRef] [Google Scholar]
  71. Zimbardo, G., Pommois, P., & Veltri, P. 2006, ApJ, 639, L91 [NASA ADS] [CrossRef] [Google Scholar]
  72. Zimbardo, G., Bitane, R., Pommois, P., & Veltri, P. 2009, Plasma Phys. Control. Fusion, 51, 015005 [NASA ADS] [CrossRef] [Google Scholar]
  73. Zimbardo, G., Perri, S., Pommois, P., & Veltri, P. 2012, Adv. Space Res. 49, 1633 [NASA ADS] [CrossRef] [Google Scholar]
  74. Zimbardo, G., Amato, E., Bovet, A., et al. 2015, J. Plasma Phys., 81, 495810601 [CrossRef] [Google Scholar]
  75. Zumofen, G., & Klafter, J. 1993, Phys. Rev. E, 47, 851 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Fractional derivatives

Fractional derivatives are integro-differential operators that allow for non-integer 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 Riemann-Liouville fractional derivatives, whose left- and right-hand 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

one-dimensional 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 time-dependent space-fractional 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 back-Fourier transform gives a Lévy distribution with power-law 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

Table 1

Solutions of normal and fractional transport equations.

In the text

All Figures

thumbnail Fig. 1

Plot of Mittag-Leffler functions Eβ( − | z | β) for several β values (see legend) in log-lin axes. Note that β = 1 corresponds to the standard exponential. A similar plot in linear axes is given by Mainardi (2014).
In the text
thumbnail Fig. 2

Plot of Mittag-Leffler functions Eβ( − | z | β) in log-log axes for several β values (see legend).
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.