EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 693 (January 2025) A&A, 693 (2025) A15 Full HTML
Open Access
Issue
A&A
Volume 693, January 2025
Article Number A15
Number of page(s) 8
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202451765
Published online 23 December 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

Cosmic rays (CRs) undergo random scattering due to magnetic field fluctuations present in astrophysical plasmas. When successive small-angle scattering leads to Brownian motion of the particles, this process can be described by Gaussian spatial diffusion along the magnetic field (Kulsrud & Pearce 1969; Skilling 1975). The time evolution of the particles’ distribution function may then be described by the transport equation (e.g., Parker 1965; Schlickeiser 2002), which has been successfully applied to model Galactic CR transport (see, e.g., Becker Tjus & Merten 2020; Mertsch 2020, for reviews).

Observations in the heliosphere, however, indicate that particle transport may not be Gaussian but instead anomalous (Perri et al. 2022). Such evidence comes from power-law profiles of ion and electron fluxes upstream of interplanetary shocks (Perri & Zimbardo 2007, 2008; Giacalone 2012; Perri et al. 2022) and at the solar wind termination shock (Perri & Zimbardo 2009, 2012) as well as solar energetic particle events (Trotta & Zimbardo 2011), which contradicts expectations for Gaussian diffusion. Anomalous diffusion may result from intermittent magnetic fields1, as particles scatter in magnetic field structures but move freely between them (Shukurov et al. 2017; Lübke et al. 2024).

Recently, we have shown how superdiffusive transport of energetic particles at shocks leads to the observed upstream power-law distributions and intensity peaks at the shock (Effenberger et al. 2024). A fractional diffusion and a Lévy flight model were introduced to describe superdiffusive transport. In this paper, we extend the Lévy flight model from Effenberger et al. (2024) to include the acceleration of particles at the shock. The expected energy spectrum and acceleration timescale give further insights into superdiffusion and its role in particle transport. We compare the energy spectra we obtain with the Lévy flight model to those obtained by standard diffusive shock acceleration and by a Lévy walk model as discussed in Perri & Zimbardo (2012).

In general, the term “anomalous diffusion” describes all dynamics for which the mean-square displacement is characterized by a non-linear dependence on time,

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

in contrast to normal diffusion (ζ = 1). Depending on the anomalous diffusion exponent ζ, the process is distinguished between subdiffusive (ζ < 1) and superdiffusive processes (ζ > 1). Usually, ζ ≤ 2 is the upper limit for superdiffusive transport, with ζ = 2 corresponding to ballistic transport or free streaming. It should be noted that the anomalous diffusion coefficient κ ∼ ζ $ \tilde{\kappa}_{\zeta} $ is in units of length2/timeζ.

Anomalous transport cannot be captured by the standard transport equation but can be described by fractional transport equations. Fractional transport equations can, for example, be solved by applying the generalized Itô lemma to obtain the corresponding stochastic differential equations (SDEs) that are driven by an alpha-stable Lévy distribution (Itô 1951; Magdziarz & Weron 2007). Section 2 gives an overview of the fractional transport equation and the Lévy flight model we employed to solve for the particle number density.

The Lévy flight model is more flexible than Fourier series approximations (Stern et al. 2014; Effenberger et al. 2024). The energy gain at the shock can be calculated along with the particle transport, the (anomalous) diffusion coefficient can be momentum dependent, and the 2D or 3D transport along magnetic field lines can be modeled with a similar simulation setup, which is not possible with the given Fourier series approximation. We present initial results for the particle acceleration at a 1D planar shock, and we contrast the constant and momentum-dependent anomalous diffusion coefficients.

For Gaussian diffusion, the energy spectrum of particles at a stationary non-relativistic planar shock with compression ratio q has the characteristic slope s = 2 − 3q/(q − 1) – for example, Drury (1983). Even for non-linear diffusion, this universal power law can be found (Walter et al. 2022). Considering Lévy flights, we find harder2 energy spectra at the shock, which is similar to what Perri & Zimbardo (2012) found for Lévy walks.

We compare the time-dependent spectra and number densities at the shock to those obtained with constant Gaussian diffusion as well as with a time-dependent Gaussian diffusion coefficient that mimics the mean-square displacement of the corresponding anomalous diffusion exponent ζ. This allows us to differentiate between effects coming from the larger mean-squared displacement of the Lévy flight model at a given time and those coming from the different underlying stochastic processes.

2. Fractional diffusion and Lévy flights

A common way to model superdiffusive processes is through the use of Lévy flights. As with Brownian motion, Lévy flights are Markov processes, but the jump length distribution follows an inverse power law (Metzler & Klafter 2000). From the jump length distribution and the assumption of finite characteristic waiting time between such jumps, a space-fractional diffusion equation can be obtained (see, e.g., Metzler & Klafter (2004), and references therein):

∂ f ( x , t ) ∂ t = κ α ∇ α f ( x , t ) , $$ \begin{aligned} \frac{\partial f(x,t)}{\partial t} = \kappa _{\alpha} \nabla ^\alpha f(x,t), \end{aligned} $$(2)

with the fractal dimension α and fractional diffusion tensor κ ̂ α $ \hat{\kappa}_{\alpha} $. The Riesz derivative (Gorenflo et al. 1999) is 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}{_x^{\alpha }} + *[_x]{D}{_{+\infty }^{\alpha }}) f(x), \end{aligned} $$(3)

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} *[_0]{D}{_t^{1-\beta }} 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} $$(4)

For α = 2, the normal diffusion equation is recovered. Figure 1 shows trajectories of Brownian motion α = 2 and α = 1.7 with the characteristic Lévy flights in 2D.

thumbnail Fig. 1.

Ten pseudoparticle trajectories starting at (0, 0) for a Gaussian diffusion process (α = 2.0, with κ2 = 1 top) and Lévy process (α = 1.7, κ1.7 = 1, bottom). We discuss how to obtain pseudoparticle trajectories and their meaning in Sect. 2.2.

The power-law distribution of the jump length leads to a diverging variance since at a given time t, an arbitrarily large distance ℓ may be traveled. However, Metzler & Klafter (2000) define a "pseudo" mean-square displacement [x2] ∝ t2/α to relate the fractional dimension α to the anomalous diffusion exponent ζ = 2/α. The fractional diffusion coefficient κα is in units of lengthα/time.

Still, the diverging mean-square displacement is an issue for the spatial transport of massive particles. This can be solved by using Lévy walks instead. By coupling the jump probability of a distance ℓ to the time t that is needed to travel such a distance, a well-defined mean-square displacement is recovered (Metzler & Klafter 2000, 2004; Zaburdaev et al. 2015). Lévy walks have already been used to model particle transport and acceleration at shocks (Perri & Zimbardo 2012). The spatio-temporal coupling leads to a different dynamical process that cannot be easily modeled with our approach. However, Lévy flights are a good approximation (we refer to Perri et al. (2015) for a comparison of a Lévy flight and Lévy walk model). Thus, we investigate how superdiffusive acceleration is different when applying a Lévy flight model compared to Lévy walks.

To investigate superdiffusive shock acceleration with Lévy flights, we solved the one-dimensional fractional transport equation while assuming isotropic particle distributions in momentum space. The time evolution of the differential number density 𝒩 = p2f(x, p, t) in space x and momentum p, with f being the particle distribution function, is determined by

∂ N ∂ t = κ α ∇ α N − u ( x ) · ∂ N ∂ x − p 3 ∂ u ∂ x ∂ N ∂ p + S ( x , p , t ) , $$ \begin{aligned} \frac{\partial \mathcal{N} }{\partial t}&= \kappa _{\alpha} \nabla ^{\alpha} \mathcal{N} - u(x) \cdot \frac{\partial \mathcal{N} }{\partial x} - \frac{p}{3} \frac{\partial u}{\partial x} \frac{\partial \mathcal{N} }{\partial p} + S(x, p, t), \end{aligned} $$(5)

where the spatial diffusion is generalized to superdiffusion, and we considered a spatially constant anomalous diffusion coefficient κα. In this macroscopic picture (e.g., Krymskii 1977; Kirk et al. 1988; Drury 1983) diffusive or superdiffusive shock acceleration results from adiabatic heating due to the divergence of the background flow u(x). Sources and sinks are given by S(x, p, t).

2.1. Solving the fractional Fokker-Planck equation with stochastic differential equations

The fractional transport equation, Eq. (5), is difficult to solve analytically due to the non-local Riesz derivative. When energy changes are neglected, the solution of the fractional diffusion equation or fractional diffusion-advection equation can, for example, be approximated by Fourier series (Stern et al. 2014; Effenberger et al. 2024). A more flexible method is to sample the solution with a Monte Carlo approach. For that, the transport equation, Eq. (5), is written into a fractional Fokker-Planck equation (FFPE), which corresponds to a set of SDEs according to Itô calculus (Itô 1951). In order to model one-dimensional superdiffusive shock acceleration, diffusion in momentum is neglected, and the diffusion coefficient κα is considered to be constant in space. The spatial displacement is then described by the SDE:

d x ( t ) = u d t + 2 κ α 1 / α d L α ( t ) . $$ \begin{aligned} \mathrm{d} x(t) = u \,\mathrm{d} t + \sqrt{2}\kappa _\alpha ^{1/\alpha }\,\mathrm{d} L_\alpha (t). \end{aligned} $$(6)

Here, dLα(t) is an α-stable Lévy process, for α = 2 is identical to a Wiener process and describes Brownian motion. The first term describes the deterministic motion due to advection, and the second term describes the stochastic motion, characterized by the fractional diffusion coefficient κα and the Lévy process dLα(t).

The change in momentum is given by an ordinary differential equation since diffusion in momentum space is neglected,

d p ( t ) = − p 3 ∂ u ∂ x d t . $$ \begin{aligned} \mathrm{d} p(t) = - \frac{p}{3} \frac{\partial u}{\partial x} \,\mathrm{d} t \;. \end{aligned} $$(7)

Equations (6) and (7) do not describe the time evolution of actual particles but sample the differential number density 𝒩.

Stochastic differential equations can be approximated numerically in their integral form (see, e.g., Kloeden & Platen 1992; Gardiner 2009). A modified version of the CR propagation framework CRPropa 3.2 (Alves Batista et al. 2022) is used to solve the stochastic differential equation in space, which uses the Euler-Maruyama scheme for integrating the SDE (6):

x n + 1 − x n = u Δ t + 2 κ α 1 / α Δ t 1 / α η α , t . $$ \begin{aligned} x_{n+1}-x_n = u \Delta t + \sqrt{2} \kappa _{\alpha }^{1/\alpha} \Delta t^{1/\alpha} \eta _{\alpha , t}. \end{aligned} $$(8)

In each simulation time step, Δt, a random number ηα, t is drawn from the α-stable Lévy distribution, and the position of the pseudoparticles is updated to xn + 1. The random number generation for ηα, t is based on the Chambers-Mallows-Stuck algorithm (Chambers et al. 1976). Lévy flights result from random numbers drawn from the heavy power-law tails of the distribution. For α = 2, the Euler-Maruyama scheme with a Wiener process is recovered.

In general, the Euler-Maruyama scheme has no constraints on the time step. However, the choice of time step can be crucial to obtaining correct results, especially when simulating diffusive shock acceleration. This is discussed in detail in such works as Kruells & Achterberg (1994), Strauss & Effenberger (2017), and Aerdker et al. (2024). Furthermore, to enhance statistics at high energies that only a fraction of pseudoparticles ever reach, we used the CandidateSplitting module of CRPropa. Pseudoparticles that cross specified boundaries in energy are split into two copies. In the later analysis they are weighted accordingly to obtain the correct spectra. Boundaries are chosen depending on the expected spectral slope to balance the number of pseudoparticles in each energy bin (see Aerdker et al. 2024, for details of this method).

The SDE approach is quite flexible, and it is easy to extend to higher dimensions and other geometries compared to Fourier series approximations. In CRPropa, Eq. (8) is defined in the lab frame. However, the diffusive step can be calculated in the orthonormal base of the magnetic field by integrating along the magnetic field line for parallel diffusion, calculating the perpendicular diffusive step, and transforming back to the lab frame in each time step (see Merten et al. 2017, for details of the transformation).

The ordinary differential equation, Eq. (7), describing the time evolution of the pseudoparticles’ momentum is integrated with a Euler scheme. CRPropa 3.2 does not allow for continuous injection of pseudoparticles during the simulation time. Instead, all pseudoparticle positions in phase-space are stored at times Ti = iΔT during simulation. After simulation, the differential number density n(x, p, T) at time T is obtained by summing over all contributions n(x, p, Ti) weighted by the time interval ΔT. The time interval ΔT does not necessarily need to be the same as the simulation time step Δt. We refer to Merten et al. (2018), Aerdker et al. (2024) for more details. This method has already been applied to model superdiffusive transport and compared to a Fourier series approximation in a previous work (Effenberger et al. 2024).

2.2. Solving the fractional diffusion equation: Two-dimensional pseudoparticle trajectories

As a first application, we obtained the pseudoparticle trajectories for Gaussian diffusion and superdiffusion in 2D shown in Fig. 1 by setting u = 0. To model the process in 2D, we made the following changes to the Euler-Maruyama scheme (8): The scattering direction was randomly chosen within [0, 2π), and the diffusive step length Δ r diff = Δ x diff 2 + Δ y diff 2 $ \Delta r_{\mathrm{diff}} = \sqrt{\Delta x_{\mathrm{diff}}^2 +\Delta y_{\mathrm{diff}}^2} $ was calculated at each time step. For this, two random numbers were drawn from the Lévy distribution with

Δ r diff = 2 κ α 1 / α η x 2 + η y 2 t 1 / α $$ \begin{aligned} \Delta r_{\mathrm{diff}} = \sqrt{2} \kappa _{\alpha }^{1/\alpha} \sqrt{\eta _x^2 + \eta _y^2} \, t^{1/\alpha} \end{aligned} $$(9)

while assuming that the diffusion coefficient κα is the same in the x and y direction and with Lévy random numbers ηα, x, ηα, y. This approach ensured a positive step length, and for α = 2, it is equal to drawing a random number from the two-dimensional chi distribution. This approach can easily be extended to three dimensions.

3. Modeling superdiffusive shock acceleration

For the superdiffusive shock acceleration (SSA) with Lévy walks, Perri & Zimbardo (2012) found that the energy spectrum is slightly harder than predicted for the normal diffusive shock acceleration (DSA). For the Lévy flights, we thus expected the spectral slope to be flatter as well but not necessarily the same as for the Lévy walks.

Walter et al. (2022) found that the spectral slope does not differ from Gaussian diffusive shock acceleration for non-linear diffusion where the diffusion coefficient depends on the distribution function. Also, for time-dependent diffusion, we did not expect a different spectral slope. However, similar to findings by Walter et al. (2022), the acceleration timescale is affected by the increasing diffusion coefficient over time.

In the macroscopic picture of the transport equation, acceleration at the shock arises from the interplay between (anomalous) diffusion and advection. In the following, we briefly explain the constraints on modeling DSA or SSA with SDEs. We then compare the results for SSA with Lévy flights to time-dependent diffusion coefficient modeling, a process with the same (pseudo) mean-square displacement.

3.1. One-dimensional planar shock

In the diffusive picture, momentum gain at shocks is described by Eq. (7) so that the advective speed u(x) must be continuously differentiable. Thus, a finite shock width is assumed instead of a discrete shock transition. The one-dimensional planar shock at x = 0 is described by

u ( x ) = u 1 + u 2 2 − u 1 − u 2 2 tanh ( x L sh ) , $$ \begin{aligned} u(x) = \frac{u_1 + u_2}{2} - \frac{u_1 - u_2}{2} \tanh \left( \frac{x}{L_{\mathrm{sh} }}\right), \end{aligned} $$(10)

with upstream and downstream velocity, u1 and u2 = u1/q, compression ratio q, and shock width Lsh. The velocity profile u(x) has also been used in other studies of DSA (Kruells & Achterberg 1994; Achterberg & Schure 2011; Walter et al. 2022; Aerdker et al. 2024).

In the following, all units have been normalized so that

x ∼ = x x 0 , u ∼ = u u 0 , t ∼ = t t 0 , κ ∼ = κ κ 0 , p ∼ = p p 0 , $$ \begin{aligned} \tilde{x} = \frac{x}{x_0},\quad \tilde{u} = \frac{u}{u_0}, \quad \tilde{t} = \frac{t}{t_0}, \quad \tilde{\kappa} = \frac{\kappa }{\kappa _0},\quad \tilde{p} = \frac{p}{p_0}, \end{aligned} $$(11)

with x0/v0 = t0 and p0 being the momentum of pseudoparticles injected at the shock. For modeling transport and acceleration in the heliosphere, the normalization can be set, for example, to x0 = 1 AU and v0 = 400 km/s. The compression ratio is assumed to be q = 4.

3.1.1. Constraints on the simulation time step

The finite shock width in Eq. (10) leads to constraints on the simulation time step and diffusion coefficient when modeling an ideal shock (for discussion on normal diffusion, see, e.g., Kruells & Achterberg 1994; Achterberg & Schure 2011; Aerdker et al. 2024). From a numerical perspective, pseudoparticles need to encounter the changing advection to be accelerated, which implies a sufficiently small time step. However, with small time steps, the diffusive step Δ x diff = 2 κ α 1 / α Δ t 1 / α $ \Delta x_{\mathrm{diff}} = \sqrt{2} \kappa_{\alpha}^{1/\alpha} \Delta t^{1/\alpha} $ – a measure for the stochastic step – may become smaller than the shock width. In that case, pseudoparticles would not make it back to the shock to be repeatedly accelerated. This constraint essentially depends on the chosen diffusion coefficient and shock width.

To model an ideal shock, the parameter ϵ = u1Lsh/κ2 should be lower than one for Gaussian diffusion. Where the diffusion coefficient is in units of lengthα/time when Lévy flights are considered. We found, however, that the constraints in the time step are less restrictive for superdiffusive transport since occasionally Lévy flights make it possible to cross the shock front again, even if the mean diffusive step is small.

3.1.2. Simulation results

We compared the time-dependent number density and momentum spectrum at the shock with superdiffusion, α = 1.7, to normal diffusion, α = 2, and to normal diffusion with a time-dependent diffusion coefficient that mimics the same dependency of the mean-square displacement on time, κ2(t) = t2/α − 1 with α = 1.7. Figure 2 shows the time evolution of the spectra fp2 at the shock (x = [0, 1]) weighted by p2 for the different diffusion processes and diffusion coefficients. The dotted lines show the fitted slope of the spectra when the stationary solution for p <  = 102 p0 is reached. The spectra were fit in the range p = [10, 102] p0 in order to exclude effects from pseudoparticle injection at p = p0.

thumbnail Fig. 2.

Time evolution of the weighted spectrum f(p, t)p4 at the shock ( x ∼ = [ 0 , 1 ] $ \tilde{x} = [0,1] $) for normal diffusion, κ ∼ 2 = 1 $ \tilde{\kappa}_2 = 1 $ (left); time-dependent diffusion, κ ∼ 2 ( t ) ∝ t ∼ 2 / 1.7 − 1 $ \tilde{\kappa}_2(t) \propto \tilde{t}^{2/1.7-1} $ (middle); and Lévy flights with α = 1.7 and κ ∼ 1.7 = 1 $ \tilde{\kappa}_{1.7} = 1 $ (right). Time-dependent diffusion and Lévy flights with α = 1.7 have the same (pseudo) mean-square displacement. The dashed lines indicate the slope of the fitted spectra at time t ∼ = 800 $ \tilde{t} = 800 $.

In the case of a constant Gaussian diffusion coefficient (left), the stationary spectrum has the expected slope −2.017 ± 0.0033, which is already reached at t ∼ = 200 $ \tilde{t} = 200 $. Acceleration at the shock with a constant diffusion coefficient is quicker compared to that with a diffusion coefficient that increases over time (middle). The time-dependent diffusion coefficient slows down the acceleration at the shock; thus, at t ∼ = 200 $ \tilde{t} = 200 $ the stationary solution for p < 102 is not yet reached, as expected for an increased diffusion coefficient in the estimate for the acceleration time τacc (see Eq. (14)). The stationary spectral slope is the same as for a constant Gaussian diffusion coefficient. The dip at low momentum (approx. p < 3) comes from the vanishing diffusion coefficient early in time (κ2(t = 0) = 0) and constraints on the chosen shock width and simulation time step.

The right panel of Fig. 2 shows the time evolution resulting from superdiffusive shock acceleration. The stationary spectrum is harder compared to normal diffusion, s1.7 = 1.750 ± 0.007 for p < 102. Also, the stationary solution at t ∼ = 200 $ \tilde{t} = 200 $ for p < 102 p0 is not yet reached. However, acceleration is more efficient than in the case of time-dependent diffusion. The results are consistent with the work of Perri & Zimbardo (2012), who also found flatter spectra for superdiffusive shock acceleration with Lévy walks.

Figure 3 compares the time evolution of the differential number density integrated over momentum. With increasing Gaussian diffusion over time, more particles reach the upstream region against the background flow compared to constant Gaussian diffusion. With Lévy flights, the characteristic power-law distributions emerge upstream. For a detailed discussion, see the prpevious work by Effenberger et al. (2024). There, the number densities differ since here the advection speed drops at the shock.

thumbnail Fig. 3.

Time evolution of the number density integrated over momentum n = ∫𝒩 dp for normal diffusion, κ ∼ 2 = 1 $ \tilde{\kappa}_2 = 1 $ (left); time-dependent diffusion, κ ∼ 2 ( t ) ∝ t ∼ 2 / 1.7 − 1 $ \tilde{\kappa}_2(t) \propto \tilde{t}^{2/1.7-1} $ (middle); and Lévy flights with α = 1.7 and κ ∼ 1.7 = 1 $ \tilde{\kappa}_{1.7} = 1 $ (right). The dashed line shows the number density at t ∼ = 800 $ \tilde{t} = 800 $.

The dashed lines show the number density profiles at t ∼ = 800 $ \tilde{t} = 800 $. Since there are no significant changes in the range x ∼ = [ − 20 , 50 ] $ \tilde{x} = [-20, 50] $, the steady-state solution is already reached for Gaussian diffusion and almost for time-dependent Gaussian diffusion and Lévy flights.

3.2. Comparison to Lévy walks

The question of why fractional diffusion changes the universal shock spectrum in contrast to time-dependent or even non-linear diffusion coefficients remains. At a one-dimensional planar shock, the spectral slope only depends on the shock compression ratio and the escape probability (see, e.g., Drury 1983). Since the compression ratio remains the same, superdiffusion modeled with Lévy flights must result in a lower escape probability, leading to more efficient acceleration.

For Lévy walks, Perri & Zimbardo (2012) using the propagator approach obtained the escape probability Pesc and with that the power law for superdiffusive shock acceleration:

γ = P esc E Δ E = 6 μ − 2 μ − 1 1 q − 1 + 1 . $$ \begin{aligned} \gamma = P_{\mathrm{esc}} \frac{E}{\Delta E} = 6 \frac{\mu - 2}{\mu - 1} \frac{1}{q - 1} + 1. \end{aligned} $$(12)

Here, ΔE/E is the relative energy gain of relativistic particles at a planar shock with compression ratio q. Their parameter μ can be related to the anomalous diffusion exponent ζ = 4 − μ. Thus, the fractional dimension α = 2/(4 − μ) has the same (pseudo) mean-square displacement. For the relation between the escape probability and spectral index, we also refer to Kirk et al. (1996), as they found softer spectra for subdiffusive shock acceleration.

Figure 4 compares the spectra of superdiffusive shock acceleration obtained with Lévy flights to those with Lévy walks with the same anomalous diffusion exponent ζ. We did not find the same spectral slopes for the Lévy flights but found even harder spectra. Thus, the jump length distribution of Lévy flights seems to lower the escape probability even more, which leads to more efficient acceleration.

thumbnail Fig. 4.

Stationary spectra for p <  = 100p0 for different fractional dimension α considering Lévy flights. On the left, the resulting energy spectra are shown. The dots are simulation results, and lines indicate the slopes fitted to the range p = (10, 100)p0. The fitted slopes are shown in comparison to the slopes for Lévy walks in the right panel, depending on the fractional dimension α and (pseudo) mean-square displacement. The fitted values for the Lévy flights are (with decreasing fractional dimension): s = −2.012 + −0.004,  s = −1.934 ± 0.003,  s = −1.844 ± 0.005 s = −1.764 ± 0.005,  s = −1.659 ± 0.013,  s = −1.576 ± 0.012.

3.3. Comparison of particle number densities

Another approach to analyzing the escape probability is to have a look at the differential number densities for the different diffusion processes. Figure 5 compares the differential number density of all considered diffusion processes in the range x ∼ = [ − 10 , 10 ] $ \tilde{x} = [-10, 10] $.

thumbnail Fig. 5.

Difference in differential number densities comparing Gaussian diffusion and Lévy flights (orange dots), constant and time-dependent Gaussian diffusion (blue squares), time-dependent Gaussian diffusion and Lévy flights (light blue triangles).

Differences between the processes are subtle. For both time-dependent Gaussian diffusion and Lévy flights, more particles are upstream and close to the shock, as indicated by the negative ratio in the upper and middle panel of Fig. 5. For Lévy flights, the upstream power laws are visible as more particles make it far upstream, even compared to time-dependent Gaussian diffusion – at least at t ∼ = 800 $ \tilde{t} = 800 $, without an upper boundary for the time-dependent diffusion coefficient κ2(t). Also, the peak structure that forms at the shock for Lévy flights is visible at x = 0 in the upper panel.

More interesting is the ratio between the diffusion processes that have the same (pseudo) mean-square displacement in the bottom panel. For Lévy flights, the particle number density far upstream is higher, indicating a more efficient scattering from the downstream region back over the shock. Also, the number density right at the shock is higher. In the macroscopic approach, energy gain at the shock comes from the divergence of the velocity profile (see Eq. (7)). With more particles close to the shock, they may be accelerated more efficiently.

3.4. Momentum-dependent diffusion

A more realistic description of the diffusion coefficient includes a dependence on the particles’ momentum. With the SDE approach, momentum-dependent diffusion in 1D can be modeled within the same framework. In the following, we show the resulting spectra at a shock while considering momentum-dependent (fractional) diffusion coefficients

κ α ( p ) = κ 0 ( p p 0 ) δ . $$ \begin{aligned} \kappa _{\alpha }(p) = \kappa _0 \left(\frac{p}{p_0}\right)^\delta . \end{aligned} $$(13)

The fractional transport equation is solved analogously to the previous section but with the diffusion coefficient in the SDE (6) now being momentum dependent.

For Gaussian diffusion, the mean time to reach momentum p depends on the momentum dependence of the diffusion coefficient (Drury 1983). For α = 1, the mean acceleration time is given by

t ¯ ( p ) ∝ τ acc ( p p 0 − 1 ) , $$ \begin{aligned} \bar{t}(p) \propto \tau _{\mathrm{acc}} \left( \frac{p}{p_0} - 1\right), \end{aligned} $$(14)

with τacc = 3/(u1 − u2)(κ1/u1 + κ2/u2). For Gaussian diffusion, the mean acceleration time we obtained for the simulation agrees with the one given by Eq. (14).

With δ > 0, momentum-dependent diffusion slows down the acceleration at the shock over time since the particles reach a higher momentum, and with that the diffusion coefficient increases. Figure 6 shows the spectrum of momentum-dependent Gaussian diffusion and Lévy flights with κ0 = 1 and δ = 1 at the same times as for constant diffusion coefficients in Fig. 2. The cutoff due to the finite acceleration time is visible in both cases.

thumbnail Fig. 6.

Time evolution of the weighted spectrum f(p, t)p4 at the shock (x = [0, 1]) for energy-dependent normal diffusion, κ2 = 1p (left) and Lévy flights with α = 1.7 and κ1.7 = 1p (right).

Considering superdiffusion, Eq. (14) is hard to define. However, Fig. (6) shows that it takes longer to reach the same momentum in case of energy-dependent Lévy flights compared to Gaussian diffusion. The Gaussian diffusion and Lévy flight process have the same value, κ0 = 1 here, but different physical units. At a given time, the Lévy flight model has a larger mean-square displacement, which slows down the acceleration process. It is, thus, difficult to compare the different processes.

One attempt by Perri & Zimbardo (2012) scales the diffusion coefficient for Lévy walks with the fractional dimension μ. This leads to a smaller anomalous diffusion coefficient, and they concluded that superdiffusive acceleration is faster than Gaussian diffusion, which is similar to our results in the previous section.

Over time, more particles make it to the upstream region due to the increasing diffusion coefficient. Figure 7 shows the number density upstream of the shock over time. The upstream number density breaks from a Gaussian core to power-law tails. While the slope of the tails is determined by the fractional dimension, the break position depends on the diffusion coefficient (see Effenberger et al. 2024). The effective diffusion coefficient grows over time when particles are accelerated to a higher momentum.

thumbnail Fig. 7.

Time evolution of the distribution f(x, t) integrated over momentum p for energy-dependent normal diffusion, κ ∼ 2 = 1 p / p 0 $ \tilde{\kappa}_2 = 1 p/p_0 $ (left) and Lévy flights with α = 1.7 and κ ∼ 1.7 = 1 p / p 0 $ \tilde{\kappa}_{1.7} = 1 p/p_0 $ (right).

4. Discussion and conclusions

We have presented a Lévy flight model of superdiffusive acceleration based on stochastic differential equations that is an extension of the Lévy flight approach described in Effenberger et al. (2024) for superdiffusive transport. The energy gain of CRs at the shock is given by Eq. (7), which was solved along with Eq. (6) to describe the Lévy flight motion of CRs. We used a modified version of CRPropa3.2 to solve the system of (stochastic) differential equations with a Euler-Maruyama scheme.

We find that the energy spectra at the shock are harder when CRs are subject to Lévy flights than for Gaussian diffusion. For Lévy flights, the spectral slope depends on the fractional dimension α but not on the anomalous diffusion coefficient. This is analogous to Gaussian diffusion, where the diffusion coefficient impacts the acceleration timescale but not the stationary spectrum.

The energy spectra are even harder than previously found by Perri & Zimbardo (2012) for Lévy walks. The spectral slope was determined by the compression ratio of the shock and the escape probability of the particles (Drury 1983). With the compression ratio being the same, the probability to escape the shock must be changed due to the power-law jump length distribution of the Lévy flights. Compared to Gaussian diffusion, more particles are upstream, and thus the transport back to the shock must be more efficient. Since Lévy flights do not have a cutoff at large jumps similar to Lévy walks, it is plausible that the resulting energy spectra are even harder for the same anomalous diffusion exponent ζ.

Comparison of such fundamentally different processes as Gaussian diffusion, Lévy flights, and Lévy walks is not trivial. To eliminate the effect of the increasingly higher mean-square displacement of the Lévy flight processes, we further compared our results to a Gaussian process with a time-dependent diffusion coefficient. The mean-square displacement of the time-dependent Gaussian and Lévy process are always the same, but the underlying scattering process is different. Compared to the time-dependent Gaussian diffusion, Lévy flights have a shorter acceleration timescale. This is different from the approach by Perri & Zimbardo (2012), who scaled down the anomalous diffusion coefficient depending on its fractional dimension for Lévy walks in order to compare the acceleration time to that of Gaussian diffusion. For future work in this context, a study of the effect of so-called tempered Lévy motion (Baeumer & Meerschaert 2010) may be of interest (i.e., the effect of exponentially truncated Lévy distributions on particle transport (see, e.g., le Roux 2024)).

This work has investigated anomalous diffusion at the level of particle distributions and energy spectra. However, it remains unclear what physical process is responsible for the non-Brownian motion. This side of the problem can be approached by, for example, studying test particle motion in magnetohydrodynamic or synthetic turbulence with intermittent coherent structures (Lübke et al. 2024). Recently, Lemoine (2023) and Kempski et al. (2023) have found that CRs can be scattered on bent magnetic field lines when their gyroradius exceeds the curvature radius of the magnetic field. Such localized, strong scattering events can lead to non-Brownian diffusion on small scales (Lemoine 2023).

To make meaningful predictions on superdiffusive particle transport in astrophysical scenarios, both the fractional dimension and anomalous diffusion coefficient must be known. For Gaussian diffusion, the diffusion coefficient may be calculated from theory (e.g., Shalchi 2020, 2021) or obtained in full-orbit test particle simulations by deriving running diffusion coefficients ⟨Δx2⟩/Δt (e.g., Mertsch 2020; Reichherzer et al. 2022). For Lévy flights, the running anomalous diffusion coefficient would converge over ⟨Δx2⟩/Δtζ. Thus, when the fractional dimension α = 2/ζ (or α = 4 − ζ for Lévy walks) is known, the anomalous diffusion coefficient can be determined from test particle simulations. For a similar study investigating subdiffusion of test particles depending on the considered turbulence, we refer to Tautz & Shalchi (2010). The anomalous diffusion coefficient and fractional dimension can also be fit to observations by taking into account the slope of the upstream power-law distribution (determined by the fractional dimension) and the distance from the shock at which the distribution turns into power laws (determined by the anomalous diffusion coefficient) (Effenberger et al. 2024).

Our work on superdiffusive transport and acceleration at 1D planar shocks sets the basis for more elaborate studies of anomalous transport. With the flexibility of the SDE approach, 3D superdiffusive transport can also be modeled analogously to Eq. (9). We have already shown in Sect. 2.2 how pseudoparticle trajectories in two dimensions can be obtained, and this can easily be extended to three dimensions. Also, anisotropic superdiffusion parallel and perpendicular to the magnetic field lines can be studied using the field line integration of CRPropa 3.2. Furthermore, the impact of a spherical shock geometry and the corresponding cooling due to the expanding wind can be taken into account in a future analysis.

1

This refers to non-Gaussian random magnetic fields governed by strong coherent structures, filaments, and voids, which are areas of weaker magnetic field strength.

2

These correspond to smaller values for the power-law slope parameter p−s, which is s > −2 here. More particles have high energies.

3

The spectra were fit for p < 102 p0 at t ∼ = 800 $ \tilde{t} = 800 $ by a least-squares method.

Acknowledgments

We acknowledge support from the DFG within the Collaborative Research Center SFB1491 “Cosmic Interacting Matters - From Source to Signal” (project number 445990517) and via the DFG grants EF 98/4-1 and FI 706/26-1.

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, JCAP, 2024, 068 [CrossRef] [Google Scholar]
  3. Alves Batista, R., Becker Tjus, J., Dörner, J., et al. 2022, JCAP, 2022, 035 [CrossRef] [Google Scholar]
  4. Baeumer, B., & Meerschaert, M. M. 2010, J. Comput. Appl. Math., 233, 2438 [NASA ADS] [CrossRef] [Google Scholar]
  5. Becker Tjus, J., & Merten, L. 2020, Phys. Rep., 872, 1 [NASA ADS] [CrossRef] [Google Scholar]
  6. Chambers, J. M., Mallows, C. L., & Stuck, B. W. 1976, J. Am. Stat. Assoc., 71, 340 [CrossRef] [Google Scholar]
  7. Drury, L. O. 1983, Rep. Prog. Phys., 46, 973 [Google Scholar]
  8. Effenberger, F., Aerdker, S., Merten, L., & Fichtner, H. 2024, A&A, 686, A219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Gardiner, C. W. 2009, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 4th edn. (Berlin: Springer) [Google Scholar]
  10. Giacalone, J. 2012, ApJ, 761, 28 [NASA ADS] [CrossRef] [Google Scholar]
  11. Gorenflo, R., De Fabritiis, G., & Mainardi, F. 1999, Phys. A: Stat. Mech. Appl., 269, 79 [Google Scholar]
  12. Itô, K. 1951, On Stochastic Differential Equations No. 4 (American Mathematical Soc.) [Google Scholar]
  13. Kempski, P., Fielding, D. B., Quataert, E., et al. 2023, MNRAS, 525, 4985 [NASA ADS] [CrossRef] [Google Scholar]
  14. Kirk, J. G., Schlickeiser, R., & Schneider, P. 1988, ApJ, 328, 269 [NASA ADS] [CrossRef] [Google Scholar]
  15. Kirk, J. G., Duffy, P., & Gallant, Y. A. 1996, A&A, 314, 1010 [NASA ADS] [Google Scholar]
  16. Kloeden, P. E., & Platen, E. 1992, Stochastic Differential Equations (Springer), 23 [Google Scholar]
  17. Kruells, W. M., & Achterberg, A. 1994, A&A, 286, 314 [NASA ADS] [Google Scholar]
  18. Krymskii, G. F. 1977, Dokl. Akad. Nauk SSSR, 234, 1306 [NASA ADS] [Google Scholar]
  19. Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 [NASA ADS] [CrossRef] [Google Scholar]
  20. le Roux, J. A. 2024, ApJ, 968, 112 [CrossRef] [Google Scholar]
  21. Lemoine, M. 2023, J. Plasma Phys., 89, 175890501 [CrossRef] [Google Scholar]
  22. Lübke, J., Effenberger, F., Wilbert, M., Fichtner, H., & Grauer, R. 2024, EPL, 146, 43001 [CrossRef] [EDP Sciences] [Google Scholar]
  23. Magdziarz, M., & Weron, A. 2007, Phys. Rev. E, 75, 056702 [Google Scholar]
  24. Merten, L., Becker Tjus, J., Fichtner, H., Eichmann, B., & Sigl, G. 2017, JCAP, 2017, 046 [CrossRef] [Google Scholar]
  25. Merten, L., Bustard, C., Zweibel, E. G., & Becker Tjus, J. 2018, ApJ, 859, 63 [CrossRef] [Google Scholar]
  26. Mertsch, P. 2020, Ap&SS, 365, 135 [NASA ADS] [CrossRef] [Google Scholar]
  27. Metzler, R., & Klafter, J. 2000, Phys. Rep., 339, 1 [CrossRef] [Google Scholar]
  28. Metzler, R., & Klafter, J. 2004, J. Phys. A Math. Gen., 37, 161 [Google Scholar]
  29. Parker, E. N. 1965, Planet. Space Sci., 13, 9 [Google Scholar]
  30. Perri, S., & Zimbardo, G. 2007, ApJ, 671, L177 [NASA ADS] [CrossRef] [Google Scholar]
  31. Perri, S., & Zimbardo, G. 2008, J. Geophys. Res.: Space Phys., 113, 3107 [CrossRef] [Google Scholar]
  32. Perri, S., & Zimbardo, G. 2009, ApJ, 693, L118 [NASA ADS] [CrossRef] [Google Scholar]
  33. Perri, S., & Zimbardo, G. 2012, ApJ, 750, 87 [NASA ADS] [CrossRef] [Google Scholar]
  34. Perri, S., Zimbardo, G., Effenberger, F., & Fichtner, H. 2015, A&A, 578, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Perri, S., Bykov, A., Fahr, H., Fichtner, H., & Giacalone, J. 2022, Space Sci. Rev., 218, 26 [NASA ADS] [CrossRef] [Google Scholar]
  36. Reichherzer, P., Merten, L., Dörner, J., et al. 2022, SN Appl. Sci., 4, 15 [NASA ADS] [CrossRef] [Google Scholar]
  37. Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Springer Science& Business Media) [CrossRef] [Google Scholar]
  38. Shalchi, A. 2020, Space Sci. Rev., 216, 23 [NASA ADS] [CrossRef] [Google Scholar]
  39. Shalchi, A. 2021, ApJ, 923, 209 [NASA ADS] [CrossRef] [Google Scholar]
  40. Shukurov, A., Snodin, A. P., Seta, A., Bushby, P. J., & Wood, T. S. 2017, ApJ, 839, L16 [NASA ADS] [CrossRef] [Google Scholar]
  41. Skilling, J. 1975, MNRAS, 172, 557 [CrossRef] [Google Scholar]
  42. Stern, R., Effenberger, F., Fichtner, H., & Schäfer, T. 2014, Fract. Calc. Appl. Anal., 17, 171 [Google Scholar]
  43. Strauss, R. D. T., & Effenberger, F. 2017, Space Sci. Rev., 212, 151 [CrossRef] [Google Scholar]
  44. Tautz, R. C., & Shalchi, A. 2010, J. Geophys. Res.: Space Phys., 115, A03104 [CrossRef] [Google Scholar]
  45. Trotta, E. M., & Zimbardo, G. 2011, A&A, 530, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Walter, D., Effenberger, F., Fichtner, H., & Litvinenko, Y. 2022, Phys. Plasmas, 29, 072302 [NASA ADS] [CrossRef] [Google Scholar]
  47. Zaburdaev, V., Denisov, S., & Klafter, J. 2015, Rev. Mod. Phys., 87, 483 [Google Scholar]

All Figures

thumbnail Fig. 1.

Ten pseudoparticle trajectories starting at (0, 0) for a Gaussian diffusion process (α = 2.0, with κ2 = 1 top) and Lévy process (α = 1.7, κ1.7 = 1, bottom). We discuss how to obtain pseudoparticle trajectories and their meaning in Sect. 2.2.
In the text
thumbnail Fig. 2.

Time evolution of the weighted spectrum f(p, t)p4 at the shock ( x ∼ = [ 0 , 1 ] $ \tilde{x} = [0,1] $) for normal diffusion, κ ∼ 2 = 1 $ \tilde{\kappa}_2 = 1 $ (left); time-dependent diffusion, κ ∼ 2 ( t ) ∝ t ∼ 2 / 1.7 − 1 $ \tilde{\kappa}_2(t) \propto \tilde{t}^{2/1.7-1} $ (middle); and Lévy flights with α = 1.7 and κ ∼ 1.7 = 1 $ \tilde{\kappa}_{1.7} = 1 $ (right). Time-dependent diffusion and Lévy flights with α = 1.7 have the same (pseudo) mean-square displacement. The dashed lines indicate the slope of the fitted spectra at time t ∼ = 800 $ \tilde{t} = 800 $.
In the text
thumbnail Fig. 3.

Time evolution of the number density integrated over momentum n = ∫𝒩 dp for normal diffusion, κ ∼ 2 = 1 $ \tilde{\kappa}_2 = 1 $ (left); time-dependent diffusion, κ ∼ 2 ( t ) ∝ t ∼ 2 / 1.7 − 1 $ \tilde{\kappa}_2(t) \propto \tilde{t}^{2/1.7-1} $ (middle); and Lévy flights with α = 1.7 and κ ∼ 1.7 = 1 $ \tilde{\kappa}_{1.7} = 1 $ (right). The dashed line shows the number density at t ∼ = 800 $ \tilde{t} = 800 $.
In the text
thumbnail Fig. 4.

Stationary spectra for p <  = 100p0 for different fractional dimension α considering Lévy flights. On the left, the resulting energy spectra are shown. The dots are simulation results, and lines indicate the slopes fitted to the range p = (10, 100)p0. The fitted slopes are shown in comparison to the slopes for Lévy walks in the right panel, depending on the fractional dimension α and (pseudo) mean-square displacement. The fitted values for the Lévy flights are (with decreasing fractional dimension): s = −2.012 + −0.004,  s = −1.934 ± 0.003,  s = −1.844 ± 0.005 s = −1.764 ± 0.005,  s = −1.659 ± 0.013,  s = −1.576 ± 0.012.
In the text
thumbnail Fig. 5.

Difference in differential number densities comparing Gaussian diffusion and Lévy flights (orange dots), constant and time-dependent Gaussian diffusion (blue squares), time-dependent Gaussian diffusion and Lévy flights (light blue triangles).
In the text
thumbnail Fig. 6.

Time evolution of the weighted spectrum f(p, t)p4 at the shock (x = [0, 1]) for energy-dependent normal diffusion, κ2 = 1p (left) and Lévy flights with α = 1.7 and κ1.7 = 1p (right).
In the text
thumbnail Fig. 7.

Time evolution of the distribution f(x, t) integrated over momentum p for energy-dependent normal diffusion, κ ∼ 2 = 1 p / p 0 $ \tilde{\kappa}_2 = 1 p/p_0 $ (left) and Lévy flights with α = 1.7 and κ ∼ 1.7 = 1 p / p 0 $ \tilde{\kappa}_{1.7} = 1 p/p_0 $ (right).
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.