EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 686 (June 2024) A&A, 686 (2024) A219 Full HTML
Open Access
Issue
A&A
Volume 686, June 2024
Article Number A219
Number of page(s) 7
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202449334
Published online 13 June 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

There is growing evidence from in situ observations of energetic particles near heliospheric shocks that the transport behavior of these particles can deviate from standard Gaussian diffusion. This includes observed power law profiles of energetic particles upstream of interplanetary shocks (Perri & Zimbardo 2007, 2008; Giacalone 2012; Perri et al. 2022), in solar energetic particle events (Trotta & Zimbardo 2011), and at the solar wind termination shock (Perri & Zimbardo 2009, 2012). This is different from Galactic cosmic ray propagation, where anomalous transport has not been explored systematically (see, e.g., Becker Tjus & Merten 2020, for a review).

Full-orbit particle simulations in weak, isotropic, random-phase Kolmogorov turbulence added to a uniform background field show different propagation regimes depending on the particles’ reduced gyro radius (Reichherzer et al. 2020, 2022a,b). All of them show either ballistic or Gaussian diffusive behavior with varying energy dependencies. However, from particle tracing in dynamo-generated fields (Shukurov et al. 2017) or synthetic turbulence models (Pucci et al. 2016) extended phases of non-diffusive behavior are observed. Zimbardo et al. (2006) found that in a generalized turbulence geometry, in a quasi-slab scenario, parallel particle transport can show prolonged superdiffusive transport, while the perpendicular propagation remains subdiffusive. Most studies to date only consider random-phase synthetic turbulence without accurate structure formation due to a lack of phase correlations (e.g., Dundovic et al. 2020; Mertsch 2020). Similarly, Magnetohydrodynamic (MHD) simulations (e.g., Beresnyak 2019), while capturing in principle the full dynamics of turbulence, exhibit only a very short and often spectrally distorted inertial range due to a lack of resolution. Next to such fundamental studies, it is thus important to also develop effective models that go beyond the paradigm of diffusive transport and that can be compared to simulation results and observations.

Anomalous diffusion encompasses both so-called subdiffusion and superdiffusion, which, to the first order, are characterized by a nonlinear temporal evolution of the particle mean-square displacement; that is, the second moment of the particle distribution. This is expressed quantitatively as

⟨ ( Δ x ) 2 ⟩ = κ ζ t ζ , $$ \begin{aligned} \langle (\Delta x)^2\rangle = \kappa _\zeta t^\zeta ,\end{aligned} $$(1)

with ζ ≠ 1 being the anomalous diffusion index and κζ a generalized diffusion constant, where one can distinguish the three diffusion regimes depending on ζ, namely subdiffusion (0 < ζ < 1), normal or Brownian or Gaussian diffusion (ζ = 1), and superdiffusion (1 < ζ < 2).

For a comprehensive overview of anomalous diffusion processes and their description, including the necessity of a generalized central limit theorem (e.g., Bouchaud & Georges 1990), see Metzler & Klafter (2000) and Metzler & Klafter (2004). More recent reviews of nonclassical transport processes in laboratory, heliophysical, and astrophysical plasmas have been presented by Perrone et al. (2013), Zimbardo et al. (2015), and Perri et al. (2022).

In the context of energetic particle transport, a number of theoretical studies have been performed, which have investigated the physical foundations and consequences of anomalous transport. Duffy et al. (1995) and Kirk et al. (1996) studied subdiffusive shock acceleration, assuming long-range coherent fields at perpendicular shocks. This way, a more rapid acceleration may be possible. Ragot & Kirk (1997) investigated the anomalous diffusion of electrons in a galaxy cluster. Kóta & Jokipii (2000) analyzed the classical concept of effective subdiffusion due to compound perpendicular diffusion. The general result is a non-Markovian perpendicular transport process with an anomalous diffusion index, ζ = 1/2, resulting from the superposition of field-line wandering and parallel transport. Webb et al. (2006) introduced (amongst other approaches) a fractional Fokker–Planck model for the compound perpendicular diffusion of cosmic rays. Shalchi & Kourakis (2007), Shalchi et al. (2007) and le Roux et al. (2010) extended the compound diffusion model to nonlinear theories. The relation between intermittent heliospheric turbulence caused by magnetic flux-tube structures and anomalous transport has been explored more systematically in recent years (e.g., Alouani-Bibi & le Roux 2014; Malandraki et al. 2019; le Roux & Zank 2021; le Roux 2022). Anomalous transport has also been studied in connection with: high-energy cosmic rays in galactic superbubbles (Barghouty & Schnee 2012); superdiffusive acceleration in the heliospheric termination shock, supernova shocks, and galaxy cluster shocks (Perri & Zimbardo 2012; Perri et al. 2016; Zimbardo & Perri 2017), yielding harder energy spectra than predicted from standard theory; and galactic propagation (Buonocore & Sen 2021; Hu et al. 2022). The acceleration model has been augmented by the consideration of Lévy walks in the studies by Zimbardo & Perri (2013) and Prete et al. (2019).

In the present paper, we aim to explicitly formulate and solve a fractional transport equation for particles near a plasma shock wave that can describe superdiffusive transport. In the following section, we introduce the transport equation based on fractional derivatives and describe the solution method based on the equivalent stochastic differential equation (SDE) with a Lévy symmetric process for the superdiffusive term. In Sect. 3, we explore the solutions to the equation. Finally, in Sect. 4 we discuss the results and future prospects for our approach.

2. Anomalous transport with fractional equations and Lévy flights

In general, subdiffusive behavior corresponds to a time-fractional diffusion equation, while superdiffusion can be described by a space-fractional diffusion operator (e.g., Chukbar 1995; Zimbardo et al. 2017). In the present work, we concentrate on superdiffusion, and thus consider a symmetric, fractional diffusion operator. This approach is complementary to a propagator approach employing asymptotic properties (e.g., Perri & Zimbardo 2009). Solving a fractional transport equation allows for the detailed computation of the distribution function of energetic particles in dependence of the phase-space coordinates and a subsequent extension to higher dimensions in future studies.

A general form of the fractional Fokker-Planck equation (FFPE) for the distribution function, f(x, t), at position x and time t, with a prescribed potential, V(x), and a generalized, spatially constant diffusion coefficient, κα, is given by (Metzler & Klafter 2000; Magdziarz & Weron 2007)

∂ f ( x , t ) ∂ t = 0 D t 1 − β ( ∂ ∂ x V ′ ( x ) + κ α ∇ α ) f ( x , t ) , $$ \begin{aligned} \frac{\partial f(x,t)}{\partial t} = {_0D}^{1-\beta }_{t}\left(\frac{\partial }{\partial x} V^\prime (x) +\kappa _\alpha \nabla ^\alpha \right) f(x,t) , \end{aligned} $$(2)

with the Riemann–Liouville fractional derivative defined as

0 D t 1 − β f ( t ) = 1 Γ ( β ) d d t ∫ 0 t ( t − s ) β − 1 f ( s ) d s , $$ \begin{aligned} {_0D}^{1-\beta }_{t} f(t) = \frac{1}{\Gamma (\beta )}\frac{\mathrm{d}}{\mathrm{d}t}\int _0^t(t-s)^{\beta -1}f(s)\,\mathrm{d}s, \end{aligned} $$(3)

and the Riesz derivative (Podlubny 1998; Metzler & Klafter 2004) given by

∇ α f ( x ) = − 1 2 cos ( α π / 2 ) ( − ∞ D x α + x D + ∞ α ) f ( x ) . $$ \begin{aligned} \nabla ^{\alpha } f(x) = -\frac{1}{2\cos (\alpha \pi /2)} ({_{-\infty }D}^{\alpha }_{x} + {_xD}_{+\infty }^{\alpha }) f(x) . \end{aligned} $$(4)

From dimensional analysis, the diffusion coefficient in Eq. (2) is related to the one in Eq. (1) via κ α = κ ζ 1/ζ $ \kappa_\alpha = \kappa_\zeta^{1/\zeta} $ with α = 2/ζ, so that the dimension of κα is (length)α/(time) and that of κζ is (length)2/(timeζ) (Metzler & Klafter 2000; Fichtner et al. 2014). In the special case of Gaussian diffusion – that is, ζ = 1 and α = 2 – one has κα = κζ ≡ κ.

The FFPE can also be derived from a generalized master equation (Metzler et al. 1999). In general, the Riesz derivative can be regarded as the proper fractional generalization of the Laplace operator, due to its symmetry properties and its Fourier transform characterization, F { ∇ α f ( x ) } = − | k | α f ∼ ( k ) $ \mathcal{F}\{\nabla^\alpha f(x)\} = -|k|^\alpha \tilde{f}(k) $, where the Fourier transform is defined in the usual way:

f ∼ ( k , t ) = F { f ( x , t ) } = 1 2 π ∫ − ∞ ∞ f ( x , t ) exp ( − i k x ) d x . $$ \begin{aligned} \tilde{f}(k,t) = \mathcal{F} \{f(x,t)\} = \frac{1}{2\pi }\int _{-\infty }^{\infty }f(x,t)\exp (-ikx)\,\mathrm{d}x. \end{aligned} $$(5)

2.1. Shock transport model

We employed a simplified model of an interplanetary or astrophysical nonrelativistic shock and the transport of energetic particles in its vicinity. The shock is assumed to be one-dimensional and planar. The considerable complication of an oblique or spherical shock and the introduction of further dimensions was left to future studies, because here we were primarily interested in the effect of superdiffusive particle propagation close to the shock. Thus, also ignoring all momentum-dependent and subdiffusive effects, the time-fractional Riemann–Liouville derivative in the FFPE gives just the identity operator (setting β = 1). For the space-fractional index, α, we allowed for values between 1 and 2. Using further a potential function giving the usual advection term and introducing a delta-functional, time-homogeneous, mono-energetic source of the particles at the origin, the FFPE (Eq. (2)) reduces to the following fractional diffusion-advection equation:

∂ ∂ t f ( x , t ) = a ∂ ∂ x f ( x , t ) + κ α ∂ α ∂ | x | α f ( x , t ) + δ ( x ) , $$ \begin{aligned} \frac{\partial }{\partial t} f(x,t) = a \frac{\partial }{\partial x} f(x,t) + \kappa _\alpha \frac{\partial ^{\alpha }}{\partial |x|^{\alpha }} f(x,t) + \delta (x) ,\end{aligned} $$(6)

where we followed the same convention for the advective term as in Litvinenko & Effenberger (2014) and Stern et al. (2014). This equation was written explicitly for only one dimension using the common partial derivative notation for the Riesz derivative. The background velocity, a, and the (fractional) diffusion coefficient, κα, were assumed to be constant.

In the following, we consider this equation as the basis of our modeling and we introduce nondimensional units. In a heliophysical scenario, this means for example that the length is given in AU, L0 = 1 AU = 1.49 × 1013 cm, and the velocity in solar wind speed units, a0 = 400 km s−1 (for our form of the diffusion advection equation, this velocity has to be negative to allow for a solar wind flowing from negative to positive x values). This gives the normalization time as the time it takes the solar wind to cross 1 AU, as t0 = 372 500 s ≈ 4.3 days. The diffusion coefficient is thus measured in units of κ0 = 5.96 × 1020 cmα/s.

In the limit of Gaussian diffusion (α = 2), Eq. (6) reduces to the usual diffusion-advection equation

∂ ∂ t f ( x , t ) = a ∂ ∂ x f ( x , t ) + κ ∂ 2 ∂ x 2 f ( x , t ) + C 0 δ ( x ) , $$ \begin{aligned} \frac{\partial }{\partial t} f(x,t) = a \frac{\partial }{\partial x} f(x,t) + \kappa \frac{\partial ^2}{\partial x^2} f(x,t) + C_0\delta (x) , \end{aligned} $$(7)

of which the steady state solution is readily given as

f ( x , ∞ ) = C 0 1 | a | exp ( − a κ x ) + C 1 , $$ \begin{aligned} f(x,\infty ) = C_0 \frac{1}{\vert a \vert }\exp \left(-\frac{a}{\kappa } x \right) + C_1, \end{aligned} $$(8)

for x < 0, and f(x, ∞) = C0|a|−1 + C1 for x > 0. Here, C0 and C1 are arbitrary integration constants that can be utilized to model the source strength and a constant particle background. This choice implies that a < 0 is equivalent to an advection in a positive x direction (i.e., opposite to the choice made in Litvinenko & Effenberger 2014 and Stern et al. 2014).

2.2. Fourier series solution

Stern et al. (2014) have found a Fourier series representation of the fundamental solution to the fractional diffusion equation (i.e., Eq. (6) without a source and advection term) on a finite periodic domain of length L:

f ( x , t ) = ∑ n = 1 ∞ ( 1 + ( − 1 ) n + 1 ) L cos ( n π x L ) · exp ( − ( n π L ) α κ α t ) . $$ \begin{aligned} f(x,t)&= \sum _{n=1}^{\infty }\frac{(1+(-1)^{n+1})}{L}\cos \left(\frac{n\pi x}{L} \right)\nonumber \\&\quad \cdot \exp \left(-\left(\frac{n \pi }{L} \right)^{\alpha } \kappa _\alpha t \right) . \end{aligned} $$(9)

This series representation coincides with the exact solution given in terms of Fox’s H-function in the limit L → ∞ (Mainardi et al. 2005). To include the advection contribution, this Fourier series solution was modified as

f ( x , t ) = ∑ n = 1 ∞ ( 1 + ( − 1 ) n + 1 ) L cos ( n π L ( x + a t ) ) · exp ( − ( n π L ) α κ α t ) . $$ \begin{aligned} f(x,t)&= \sum _{n=1}^{\infty } \frac{(1+(-1)^{n+1})}{L} \cos \left(\frac{n \pi }{L} (x+at) \right)\nonumber \\&\quad \cdot \exp \left(- \left(\frac{n \pi }{L} \right)^{\alpha } \kappa _\alpha t \right) . \end{aligned} $$(10)

Finally, to obtain the full solution to Eq. (6) including the source, the previous expression can be integrated with respect to time, which yields

f ( x , t ) = ∫ 0 t ∑ n = 1 ∞ ( 1 + ( − 1 ) n + 1 ) L cos ( n π L ( x + a t ′ ) ) · exp ( − ( n π L ) α κ α t ′ ) d t ′ = ∑ n = 1 ∞ { ( 1 + ( − 1 ) n + 1 ) [ ( n π L ) α κ α L cos ( n π x L ) − n π a sin ( n π L x ) − ( n π L ) α κ α L cos ( n π L ( x + a t ) ) exp ( − ( n π L ) α κ α t ) + n π a sin ( n π L ( x + a t ) ) exp ( − ( n π L ) α κ α t ) ] / ( ( n π L ) 2 α κ α 2 L 2 + n 2 π 2 a 2 ) } . $$ \begin{aligned} f(x,t)&= \int _{0}^{t} \sum _{n=1}^{\infty } \frac{(1+(-1)^{n+1})}{L} \cos \left(\frac{n \pi }{L} (x+at^\prime ) \right)\nonumber \\&\cdot \exp \left(- \left(\frac{n \pi }{L} \right)^{\alpha } \kappa _\alpha t^\prime \right) \mathrm{d}t^\prime \nonumber \\&=\sum _{n=1}^{\infty }\Biggl \{(1+(-1)^{n+1}) \biggl [ \left(\frac{n \pi }{L}\right)^{\alpha } \kappa _\alpha L \cos \left(\frac{n \pi x}{L} \right)\nonumber \\&-n \pi a \sin \left(\frac{n \pi }{L}x \right) \nonumber \\&-\left(\frac{n \pi }{L} \right)^{\alpha } \kappa _\alpha L\cos \left(\frac{n \pi }{L}(x+at) \right) \exp \left(-\left(\frac{n \pi }{L} \right)^{\alpha } \kappa _\alpha t \right) \nonumber \\&+ n \pi a \sin \left( \frac{n \pi }{L} (x+at) \right) \exp \left(-\left(\frac{n \pi }{L} \right)^{\alpha } \kappa _\alpha t \right) \biggr ] {\Bigg /} \nonumber \\&\left( \left(\frac{n \pi }{L}\right)^{2 \alpha }\kappa _\alpha ^2L^2 + n^2 \pi ^2a^2 \right) \Biggr \} . \end{aligned} $$(11)

The convergence of this series depends on the period, L, with f(x, t) = f(x + 2L, t) and consequently it does not converge on the exact solution for t → ∞, since the periodic domains begin to overlap. In practice, L has to be chosen to be sufficiently large to prevent these influences. We used this solution to validate our Lévy flight model.

2.3. Numerical method based on a stochastic differential equation

We used a modified version of CRPropa3.2 (Alves Batista et al. 2022) based on the equivalence between a certain class of SDEs and the Fokker–Planck equation (Gardiner 2009). The standard SDE methods have become increasingly popular in the energetic particle modeling context due to their simplicity and scalability with modern computer architecture (see, e.g., Kopp et al. 2012 for an account of numerical techniques). They have been applied successfully to, for example, solar energetic particle transport (Dröge et al. 2010) and cosmic ray modulation problems (Strauss et al. 2011), as well as high-energy cosmic ray transport. including acceleration in the case of Gaussian diffusion (Achterberg & Schure 2011; Aerdker et al. 2024). Strauss & Effenberger (2017) give an extended review of the relevant methods and applications.

For the full FFPE (Eq. (2)), an SDE-type method has been developed by Magdziarz & Weron (2007) that resembles a competition between subdiffusion and Lévy flights. In this approach, the subdiffusive aspect is represented by a time rescaling of the stochastic process, called subordination, while the space-fractional aspect of Lévy flights is modeled by an α-stable random distribution, which deviates from the Gaussian exponential behavior. Some more details on particle transport modeling with this approach are also described in Effenberger (2014). Since in our model we are only considering superdiffusion, the main aspect is the generalization to α-stable distributions. This gives the following SDE,

d X ( t ) = − a d t + 2 κ α 1 / α d L α ( t ) , $$ \begin{aligned} \mathrm{d}X(t) = - a \,\mathrm{d}t + \sqrt{2}\kappa _\alpha ^{1/\alpha }\,\mathrm{d}L_\alpha (t) , \end{aligned} $$(12)

where dLα(t) is a Lévy α-stable distribution that fulfils the Fourier transform characteristic

F { e i k L α ( t ) } = e − t | k | α . $$ \begin{aligned} \mathcal{F} \{e^{ikL_{\alpha }(t)}\} = e^{-t|k|^{\alpha }} . \end{aligned} $$(13)

In each simulation time step, Δt, a random number, ηα, is drawn from the Lévy α-stable distribution and the pseudo-particle position is updated based on an advective step, −aΔt, and a diffusive step, 2 κ α 1 / α η α Δ t 1 / α $ \sqrt{2}\kappa_\alpha^{1/\alpha}\eta_{\alpha} \Delta t^{1/\alpha} $. The random number generation is based on the Chambers–Mallows–Stuck method (Chambers et al. 1976), which is explained in more detail in Appendix A.

Figure 1 shows a comparison between pseudo-particle orbits driven by a Wiener process (α = 2) and a Lévy motion (α = 1.7) for a simple SDE with a = 0 and κα = 1 obtained with CRPropa3.2. All pseudo-particles are injected at t = 0 and followed for a unit time interval. The trajectories show the prominent Lévy flight behavior of the pseudo-particles, giving occasional large spatial jumps that result from the heavy tails in the random distribution (compare to Fig. A.1). These jumps become more pronounced and frequent for smaller values of α. It has to be emphasized at this point that these trajectories represent only samples of the stochastic process and not real particle trajectories. Their jumps thus do not violate any finite particle speed effects. Still, the mean squared displacement of the Lévy distribution diverges (Metzler & Klafter 2000). This is a potential issue for the spatial transport of massive particles, which could be solved by replacing Lévy flights with Lévy walks (Prete et al. 2019, 2021). Lévy walks are spatiotemporal-coupled and long jumps get a time penalty, which recovers a well-defined mean squared displacement (see Metzler & Klafter 2000, 2004 and references therein).

thumbnail Fig. 1.

Five pseudo-particle orbits starting at x = 0 for a Gaussian process (α = 2.0, top panel) and a Lévy process (α = 1.7, bottom panel) for a unit diffusion coefficient, κα, and vanishing background velocity, a.

The solution to the fractional diffusion-advection (Eq. (6)) is given by the conditional probability of the random process described by the SDE (Eq. (12)). In practice, the distribution function, f(x, t), is obtained from an ensemble average over a large number of pseudo-particle trajectories sampling the stochastic process. We only consider the time-forward approach here, although in principle one can also integrate trajectories backward in time from the phase-space point of interest (see again Kopp et al. 2012 for details on the difference between forward and backward methods in the classical SDE context).

To model the chosen source term in the transport equation, pseudo-particles need to be continuously injected at the origin. The structure of CRPropa3.2, however, does not allow for an intrinsic description of such continuous injection. Instead, the pseudo-particles position in phase-space are stored at times Ti = iΔT1 during the simulation. In the later analysis, the distribution function, f(x, T), at time T can be constructed by summing over all distributions, f(x, Ti), weighted by the time interval, ΔT. The time resolution, ΔT, is sufficient, if the distribution function does not change too much during the time interval, ΔT. For more details, we refer to Merten et al. (2017, 2018), Aerdker et al. (2024).

3. Results for the transport of shock particles

To validate our solution methods presented in the previous chapter, we compared them in the framework of our shock model. We considered the time evolution of the spatial structure of the distribution function for Gaussian diffusion and for a case of superdiffusion with α = 1.7, which is similar to the values estimated in the data analysis of Perri & Zimbardo (2009) for the termination shock.

The result at three different times for Gaussian diffusion is given in the left panel of Fig. 2 in a semi-logarithmic plot. The solid lines giving the Fourier series solution (as is presented in Sect. 2.2) coincide very well with the symbols, which show the numerical solution to the SDE (Eq. (12)) calculated with a total number of 104 pseudo-particles released at the origin. In case of Gaussian diffusion, the steady state is reached at t = 100 for x < 30.

thumbnail Fig. 2.

Fourier series (solid lines) and SDE (symbols) solution of the distribution function, f, for the shock setup with Gaussian diffusion (α = 2, left panel) and superdiffusion (α = 1.7, right panel). We took the parameters a = −1 (corresponding to an advection in positive x direction) and κα = 1, and f is given at times t1 = 1, t2 = 10, and t3 = 100. The dashed green line indicates the time-asymptotic steady state solution for Gaussian diffusion (Eq. (8)).

The resulting distributions in the superdiffusive case (Fig. 2, right panel) show an equally good agreement between both solution methods. In particular, the upstream solution shows some prominent deviation from the exponential behavior of the classical result. A peaked structure at the origin is formed as well. This is visible more clearly in Fig. 3, where the peak is shown for different values of α.

thumbnail Fig. 3.

Fourier series (lines) and SDE (symbols) solution at t = 100 close to the shock at x = 0 for different values of α. The peak at the shock is more distinct for lower values of α. The solution for Gaussian diffusion is shown by the dotted line.

We find that for different values of κα (keeping the velocity at unity) the upstream Fourier series solution has different turning distances, as is illustrated in Fig. 4. It can be seen that close to the shock the solution resembles more the behavior of a Gaussian diffusion model before turning into a power law for larger distances. This enables us to estimate possible values for the superdiffusion coefficient, κα, from the turnover into the power law in comparison with data (Perri et al. 2015). In Fig. 4, the upstream solution of the distribution is compared for different values of the diffusion coefficient, κα, for Gaussian diffusion and superdiffusive transport. It can clearly be seen that the turnover into the power laws depends on the (super)diffusion coefficient.

thumbnail Fig. 4.

Upstream Fourier series solution (Eq. (11)) (line) compared to the corresponding SDE solution (symbols) at t = 100 with different values of κα, namely κα = 0.5 (dot), κα = 1 (square), and κα = 2 (triangle). Superdiffusion with α = 1.7 (solid line) is compared to Gaussian diffusion (dotted lines).

This behavior is also clearly visible in Fig. 5, where we investigate the transition to the power law as dependent on α. It can be seen that the break scale is dependent on alpha, while the Gaussian diffusion case of α = 2 remains an exponential rollover that is matched closely by the other cases in the vicinity of the shock.

thumbnail Fig. 5.

Upstream Fourier series solution (Eq. (11)) (line) compared to the corresponding SDE solution (symbols) at t = 100 with different values of α and κα = 1. Superdiffusion is compared to Gaussian diffusion (dotted line).

4. Conclusions

We have presented a stochastic solution method based on Lévy flights for the fractional diffusion-advection equation for cosmic ray transport. With the modified CRPropa3.2 code2, we are able to reproduce time-dependent spatial profiles of the full solution of the space-fractional equation as ensemble averages of pseudo-particles. Our numerical solutions are validated by semi-analytic Fourier series solutions of the space-fractional equation.

A particular feature of these solutions is a pronounced peak in the spatial profile near the shock position. Figure 3 illustrates this and the dependence of the feature on the superdiffusivity parameter, α. This may provide a further observational constraint on superdiffusion models and can also be of relevance to the effective acceleration spectra in a superdiffusion scenario. While the acceleration of particles itself was beyond the scope of our present study, we aim to investigate the superdiffusive acceleration of particles in a comprehensive SDE model in future work.

We emphasize the advantages of the SDE approach in such a context. For superdiffusive shock acceleration, a second ordinary differential equation can be solved along with Eq. (12), describing the energy gain of pseudo-particles due to the divergence of the advective background flow, a(x), at the shock (for more details on simulating diffusive shock acceleration with SDEs, see e.g., Achterberg & Schure 2011 and Aerdker et al. 2024). An additional advantage of the SDE approach is that three-dimensional superdiffusive transport could also be modeled by a system of SDEs, analogously to Eq. (12). This would enable the study of spherical shocks and, potentially, galactic transport.

In summary, a quantitative model of superdiffusive energetic particle transport, such as the one presented in this paper, enables us to bridge our theoretical understanding of anomalous diffusion with specific questions related to the modeling of observed energetic particle profiles. Combined with new insights from full-orbit simulations in turbulence, extensions to higher dimensions, and the inclusion of particle acceleration, future work should give us a more complete picture of anomalous diffusion in space plasmas.

1

Note, that Δt and ΔT do not necessarily need to be the same.

2

We used the CRPropa3.2 version with the git hash 83a85e54, which will become publicly available in the future.

Acknowledgments

We acknowledge support from the German Science Foundation (DFG; SFB1491, and grants EF98/4-1 and FI706/26-1) and the International Space Science Institute (ISSI). Preliminary work was already published in the PhD thesis by FE. We thank the participants of a workshop on the topics of anomalous diffusion held in Bochum in May 2023 for helpful discussions.

References

  1. Achterberg, A., & Schure, K. M. 2011, MNRAS, 411, 2628 [CrossRef] [Google Scholar]
  2. Aerdker, S., Merten, L., Becker Tjus, J., et al. 2024, J. Cosmol. Astropart. Phys., 2024, 068 [CrossRef] [Google Scholar]
  3. Alouani-Bibi, F., & le Roux, J. A. 2014, ApJ, 781, 93 [NASA ADS] [CrossRef] [Google Scholar]
  4. Alves Batista, R., Becker Tjus, J., Dörner, J., et al. 2022, J. Cosmol. Astropart. Phys., 2022, 035 [CrossRef] [Google Scholar]
  5. Barghouty, A. F., & Schnee, D. A. 2012, ApJ, 749, 178 [NASA ADS] [CrossRef] [Google Scholar]
  6. Becker Tjus, J., & Merten, L. 2020, Phys. Rep., 872, 1 [NASA ADS] [CrossRef] [Google Scholar]
  7. Beresnyak, A. 2019, Liv. Rev. Comput. Astrophys., 5, 2 [CrossRef] [Google Scholar]
  8. Bouchaud, J., & Georges, A. 1990, Phys. Rep., 195, 127 [CrossRef] [Google Scholar]
  9. Buonocore, S., & Sen, M. 2021, AIP Adv., 11, 055221 [NASA ADS] [CrossRef] [Google Scholar]
  10. Chambers, J. M., Mallows, C. L., & Stuck, B. W. 1976, J. Am. Stat. Assoc., 71, 340 [CrossRef] [Google Scholar]
  11. Chukbar, K. V. 1995, Sov. J. Exp. Theoret. Phys., 81, 1025 [NASA ADS] [Google Scholar]
  12. Dröge, W., Kartavykh, Y. Y., Klecker, B., & Kovaltsov, G. A. 2010, ApJ, 709, 912 [CrossRef] [Google Scholar]
  13. Duffy, P., Kirk, J. G., Gallant, Y. A., & Dendy, R. O. 1995, A&A, 302, L21 [NASA ADS] [Google Scholar]
  14. Dundovic, A., Pezzi, O., Blasi, P., Evoli, C., & Matthaeus, W. H. 2020, Phys. Rev. D, 102, 103016 [NASA ADS] [CrossRef] [Google Scholar]
  15. Effenberger, F. 2014, in Numerical Modeling of Space Plasma Flows, Astronum-2013, ASP Conf. Ser., 488, 195 [Google Scholar]
  16. Fichtner, H., Stern, R., & Effenberger, F. 2014, in 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013), eds. N. V. Pogorelov, E. Audit, & G. P. Zank, ASP Conf. Ser., 488, 17 [Google Scholar]
  17. Gardiner, C. W. 2009, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (Berlin: Springer) [Google Scholar]
  18. Giacalone, J. 2012, ApJ, 761, 28 [NASA ADS] [CrossRef] [Google Scholar]
  19. Hu, Y., Lazarian, A., & Xu, S. 2022, MNRAS, 512, 2111 [NASA ADS] [CrossRef] [Google Scholar]
  20. Kirk, J. G., Duffy, P., & Gallant, Y. A. 1996, A&A, 314, 1010 [NASA ADS] [Google Scholar]
  21. Kopp, A., Büsching, I., Strauss, R. D., & Potgieter, M. S. 2012, Comput. Phys. Commun., 183, 530 [NASA ADS] [CrossRef] [Google Scholar]
  22. Kóta, J., & Jokipii, J. R. 2000, ApJ, 531, 1067 [Google Scholar]
  23. le Roux, J. A. 2022, ApJ, 930, 125 [CrossRef] [Google Scholar]
  24. le Roux, J. A., & Zank, G. P. 2021, ApJ, 913, 84 [CrossRef] [Google Scholar]
  25. le Roux, J. A., Webb, G. M., Shalchi, A., & Zank, G. P. 2010, ApJ, 716, 671 [NASA ADS] [CrossRef] [Google Scholar]
  26. Litvinenko, Y. E., & Effenberger, F. 2014, ApJ, 796, 125 [NASA ADS] [CrossRef] [Google Scholar]
  27. Magdziarz, M., & Weron, A. 2007, Phys. Rev. E, 75, 056702P [CrossRef] [Google Scholar]
  28. Mainardi, F., Pagnini, G., & Saxena, R. K. 2005, J. Comput. Appl. Math., 178, 321 [NASA ADS] [CrossRef] [Google Scholar]
  29. Malandraki, O., Khabarova, O., Bruno, R., et al. 2019, ApJ, 881, 116 [Google Scholar]
  30. Merten, L., Becker Tjus, J., Fichtner, H., Eichmann, B., & Sigl, G. 2017, J. Cosmol. Astropart. Phys., 2017, 046 [CrossRef] [Google Scholar]
  31. Merten, L., Bustard, C., Zweibel, E. G., & Becker Tjus, J. 2018, ApJ, 859, 63 [CrossRef] [Google Scholar]
  32. Mertsch, P. 2020, Ap&SS, 365, 135 [NASA ADS] [CrossRef] [Google Scholar]
  33. Metzler, R., & Klafter, J. 2000, Phys. Rep., 339, 1 [CrossRef] [Google Scholar]
  34. Metzler, R., & Klafter, J. 2004, J. Phys. A Math. Gen., 37, 161 [Google Scholar]
  35. Metzler, R., Barkai, E., & Klafter, J. 1999, EPL (Europhysics Letters), 46, 431 [NASA ADS] [CrossRef] [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. (Space Phys.), 113, 3107 [NASA ADS] [CrossRef] [Google Scholar]
  38. Perri, S., & Zimbardo, G. 2009, ApJ, 693, L118 [NASA ADS] [CrossRef] [Google Scholar]
  39. Perri, S., & Zimbardo, G. 2012, ApJ, 750, 87 [NASA ADS] [CrossRef] [Google Scholar]
  40. Perri, S., Zimbardo, G., Effenberger, F., & Fichtner, H. 2015, A&A, 578, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Perri, S., Amato, E., & Zimbardo, G. 2016, A&A, 596, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Perri, S., Bykov, A., Fahr, H., Fichtner, H., & Giacalone, J. 2022, Space Sci. Rev., 218, 26 [NASA ADS] [CrossRef] [Google Scholar]
  43. Perrone, D., Dendy, R. O., Furno, I., et al. 2013, Space Sci. Rev., 178, 233 [NASA ADS] [CrossRef] [Google Scholar]
  44. Podlubny, I. 1998, Fractional Differential Equations: Mathematics in Science and Engineering (Elsevier Science) [Google Scholar]
  45. Prete, G., Perri, S., & Zimbardo, G. 2019, AdSpR, 63, 2659 [NASA ADS] [Google Scholar]
  46. Prete, G., Perri, S., & Zimbardo, G. 2021, Nat. Astron., 87, 101605 [Google Scholar]
  47. Pucci, F., Malara, F., Perri, S., et al. 2016, MNRAS, 459, 3395 [NASA ADS] [CrossRef] [Google Scholar]
  48. Ragot, B. R., & Kirk, J. G. 1997, A&A, 327, 432 [NASA ADS] [Google Scholar]
  49. Reichherzer, P., Becker Tjus, J., Zweibel, E. G., Merten, L., & Pueschel, M. J. 2020, MNRAS, 498, 5051 [NASA ADS] [CrossRef] [Google Scholar]
  50. Reichherzer, P., Becker Tjus, J., Zweibel, E. G., Merten, L., & Pueschel, M. J. 2022a, MNRAS, 514, 2658 [NASA ADS] [CrossRef] [Google Scholar]
  51. Reichherzer, P., Merten, L., Dörner, J., et al. 2022b, SN Appl. Sci., 4, 15 [NASA ADS] [CrossRef] [Google Scholar]
  52. Shalchi, A., & Kourakis, I. 2007, A&A, 470, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Shalchi, A., Kourakis, I., & Dosch, A. 2007, J. Phys. A Math. Gen., 40, 11191 [NASA ADS] [CrossRef] [Google Scholar]
  54. Shukurov, A., Snodin, A. P., Seta, A., Bushby, P. J., & Wood, T. S. 2017, ApJ, 839, L16 [NASA ADS] [CrossRef] [Google Scholar]
  55. Stern, R., Effenberger, F., Fichtner, H., & Schäfer, T. 2014, Fract. Calcul. Appl. Anal., 17, 171 [CrossRef] [Google Scholar]
  56. Strauss, R. D. T., & Effenberger, F. 2017, Space Sci. Rev., 212, 151 [CrossRef] [Google Scholar]
  57. Strauss, R. D., Potgieter, M. S., Büsching, I., & Kopp, A. 2011, ApJ, 735, 83 [NASA ADS] [CrossRef] [Google Scholar]
  58. Trotta, E. M., & Zimbardo, G. 2011, A&A, 530, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Webb, G. M., Zank, G. P., Kaghashvili, E. K., & le Roux, J. A. 2006, ApJ, 651, 211 [Google Scholar]
  60. Zimbardo, G., & Perri, S. 2013, ApJ, 778, 35 [CrossRef] [Google Scholar]
  61. Zimbardo, G., Pommois, P., & Veltri, P. 2006, ApJ, 639, L91 [NASA ADS] [CrossRef] [Google Scholar]
  62. Zimbardo, G., Amato, E., Bovet, A., et al. 2015, J. Plasma Phys., 81, 495810601 [NASA ADS] [CrossRef] [Google Scholar]
  63. Zimbardo, G., & Perri, S. 2017, Nat. Astron., 1, 0163 [Google Scholar]
  64. Zimbardo, G., Perri, S., Effenberger, F., & Fichtner, H. 2017, A&A, 607, A7 [Google Scholar]

Appendix A: Random number generation

In the modified version of CRPropa3.2, a random number, ηα, is drawn from the α-stable Lévy distribution for each time step, Δt, and for every pseudo-particle to get the diffusive step Δ x diff = 2 ( κ α Δ t ) 1 / α η α $ \Delta x_{\mathrm{diff}} = \sqrt{2}\left(\kappa_\alpha \Delta t \right)^{1/\alpha} \eta_{\alpha} $. In Figure A.1 pseudo-random numbers (PRNs) drawn from an α-stable Lévy distribution with α = 1.5 and α = 2 are shown. Compared to the Gaussian distribution, the α-stable Lévy distribution with α ≠ 2 has enhanced tails and an enhanced center.

thumbnail Fig. A.1.

Histogram (200 bins) of 107 PRNs drawn from a symmetric α-stable distribution with α = 1.5 and α = 2.

The algorithm from Chambers et al. (1976) (see also Magdziarz & Weron 2007) was used to generate the PRNs:

η α = 1 2 sin ( α U ) [ cos ( U ) ] 1 / α [ cos ( ( 1 − α ) U ) V ] ( 1 − α ) / α . $$ \begin{aligned} \eta _{\alpha } = \frac{1}{\sqrt{2}} \frac{\sin (\alpha U)}{ \left[\cos (U)\right]^{1/\alpha }} \left[ \frac{\cos ((1 - \alpha ) U)}{ V} \right]^{(1-\alpha )/\alpha }. \end{aligned} $$(A.1)

Here, U = π(ξ − 0.5) and V = −log(ξ), and ξ is a random number uniformly distributed on (0, 1).

Appendix B: Delta injection

The power law asymptotics of the superdiffusive process can already be studied in the simpler case of pure fractional diffusion,

∂ f ∂ t = κ α ∇ α f , f ( x , 0 ) = δ ( x ) . $$ \begin{aligned} \frac{\partial f}{\partial t} = \kappa _\alpha \nabla ^{\alpha } f,\quad f(x, 0) = \delta (x). \end{aligned} $$(B.1)

Figure B.1 shows the solution of the distribution function, f, decaying over time as compared to the Fourier series approximation (Equation 9) with excellent agreement, giving further validation of our method.

thumbnail Fig. B.1.

Fourier series (dash-dotted line) and SDE (dots) solution of the distribution function for the delta injection of pseudo-particles at t = 0.

All Figures

thumbnail Fig. 1.

Five pseudo-particle orbits starting at x = 0 for a Gaussian process (α = 2.0, top panel) and a Lévy process (α = 1.7, bottom panel) for a unit diffusion coefficient, κα, and vanishing background velocity, a.
In the text
thumbnail Fig. 2.

Fourier series (solid lines) and SDE (symbols) solution of the distribution function, f, for the shock setup with Gaussian diffusion (α = 2, left panel) and superdiffusion (α = 1.7, right panel). We took the parameters a = −1 (corresponding to an advection in positive x direction) and κα = 1, and f is given at times t1 = 1, t2 = 10, and t3 = 100. The dashed green line indicates the time-asymptotic steady state solution for Gaussian diffusion (Eq. (8)).
In the text
thumbnail Fig. 3.

Fourier series (lines) and SDE (symbols) solution at t = 100 close to the shock at x = 0 for different values of α. The peak at the shock is more distinct for lower values of α. The solution for Gaussian diffusion is shown by the dotted line.
In the text
thumbnail Fig. 4.

Upstream Fourier series solution (Eq. (11)) (line) compared to the corresponding SDE solution (symbols) at t = 100 with different values of κα, namely κα = 0.5 (dot), κα = 1 (square), and κα = 2 (triangle). Superdiffusion with α = 1.7 (solid line) is compared to Gaussian diffusion (dotted lines).
In the text
thumbnail Fig. 5.

Upstream Fourier series solution (Eq. (11)) (line) compared to the corresponding SDE solution (symbols) at t = 100 with different values of α and κα = 1. Superdiffusion is compared to Gaussian diffusion (dotted line).
In the text
thumbnail Fig. A.1.

Histogram (200 bins) of 107 PRNs drawn from a symmetric α-stable distribution with α = 1.5 and α = 2.
In the text
thumbnail Fig. B.1.

Fourier series (dash-dotted line) and SDE (dots) solution of the distribution function for the delta injection of pseudo-particles at t = 0.
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.