© ESO, 2015

1. Introduction

For many decades now, astrophysicists have investigated the origin of cosmic rays, the high-energy charged particle radiation that permeates space. Various experiments, both earthbound and space-borne, measure a characteristic cosmic-ray energy spectrum, which follows a (more-or-less) universal power law. The individual particles have kinetic energies above several MeV, and even energies up to more than 10²⁰ eV have been detected experimentally, as shown by air-shower experiments (e.g., Abraham et al. 2008; Letessier-Selvon & Stanev 2011). Overall, the data come from a multitude of astrophysical experiments that have been conducted over past decades. The understanding of the origin of the distinct shape of the spectrum continues to be an active field of study, and new experiments like the Square Kilometer Array (SKA, e.g., Falcke et al. 2004) are being assembled.

While ineffective on its own, the original Fermi (1949) acceleration mechanism inspired later research, leading to the development of the so-called diffusive shock acceleration mechanism or of first-order Fermi acceleration (Axford et al. 1977; Krymsky 1977; Blandford & Ostriker 1978; Bell 1978). This theory outlines an effective acceleration mechanism for charged particles at magneto-hydrodynamic shock waves. Such shock waves of varying sizes are ubiquitous throughout space, from solar winds through supernova remnants (SNR) to relativistic shocks in exotic cosmic objects, such as active galactic nuclei or pulsars (see Balogh & Treumann 2013, for an introduction).

Thanks to the efficiency of the diffusive shock acceleration mechanism and the existence of plausible acceleration sites, it became accepted as the de-facto standard acceleration mechanism (cf. Abdo et al. 2010; Helder et al. 2012). There are, however, still many open questions that cannot be answered entirely yet. These include direct proof for acceleration at the shock, the injection problem, and the origin of ultra-high energy cosmic rays has not yet been understood (e.g., Abbasi et al. 2014). In general, the complexity of the physical problem at hand requires the application of multiple approximations to find results both analytically and numerically.

The work presented here will concentrate on non-relativistic shock waves in SNRs (for a review see, e.g., Reynolds 2008; Blasi 2013). Based on the numerical model used in the Padian code (Tautz 2010a), a test-particle simulation is extended to include shock waves. This approach is then applied in order to investigate the effect of diffusive shock acceleration by inspecting the obtained cosmic-ray energy spectra. This relatively simple test-particle model is compared to theory, on one hand, and the results of Caprioli & Spitkovsky (2014), on the other hand, who applied a sophisticated kinetic hybrid model to investigate particle acceleration at non-relativistic astrophysical shocks. Finally, the work presented here looks at the influence of the model of turbulence on the acceleration process. Most notably, a simple magnetostatic model is compared to an extended model based on Alfvén waves, following the research of Schlickeiser (2002) on the transmission of these waves through a parallel shock front.

This article is organized as follows. In Sect. 2, a brief outline of diffusive shock acceleration is given, followed by a summary of the numerical model in Sect. 3. The results are presented in Sect. 4. Section 5 provides a brief summary and a discussion of the results.

2. Diffusive shock acceleration

Shock waves are a special case of a general physical phenomenon, where the physical properties of a (magneto-)hydrodynamic environment change discontinuously across a surface area. They are characterized by being propagating structures; i.e., they have a non-zero mass flux across the discontinuity surface, which is commonly referred to as the shock front. In reality, a shock front is an irregular surface area and the discontinuity will occur over a certain small, but non-zero, distance. Any realistic description of such a system would be highly complex. Instead, one often approximates the shock front to be one-dimensional, planar, and to have infinitesimally small thickness. This simplification is justified on any scale far smaller than the total scale of the shock wave but much larger than the shock front’s thickness.

Shock waves with very high upstream Mach numbers M₁ ≫ 1 are called strong shock waves. In this limit, the compression ratio r, i.e. the ratio between the upstream density ρ₁ and the downstream density ρ₂, can be written as (e.g., Schlickeiser 2002, Sect. 16.2.2.1) $\mathit{r}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{1}}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{1}}}{{\mathit{v}}_{\mathrm{2}}}\mathrm{=}\frac{\mathit{\gamma }\mathrm{+}\mathrm{1}}{\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{\beta }\mathit{/}{\mathit{M}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{\approx }\frac{\mathit{\gamma }\mathrm{+}\mathrm{1}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{·}$(1)Here, the plasma β is the ratio of the plasma pressure to the magnetic pressure, which will be assumed to be negligible small compared to the square of the shock Mach number ${\mathit{M}}_{\mathrm{1}}^{\mathrm{2}}$. The right-hand side is written in terms of the ratio of specific heats of the gas, γ. The compression ratio reaches its minimum of r = 4 for a monoatomic gas. An example for a strong shock wave is the shell of a SNR, which expands with a velocity in the order of 10⁴ km s^-1. Compared to the speed of sound, the strong shock approximation is plausible.

Shocks with a velocity V_shock larger than the Alfvén speed V_A, which can be approximated to lie in the range of 30 km s^-1 to 50 km s^-1 in the interstellar medium (cf. Kang 2013), are said to be super-Alfvénic. These fast shocks are common elements of astrophysical environments. They can be found for example in s, the relativistic bulk outflow from active nuclei, pulsar wind nebulae, or gamma-ray bursts.

Averaging over particles crossing the shock at arbitrary angles and taking into account that, due to scattering processes both in the upstream and the downstream regions, the velocity vector will be quickly randomized, a particle’s energy gain for a complete round-trip is given by $⟨\frac{\mathit{\Delta E}}{\mathit{E}}⟩\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\frac{{\mathit{V}}_{\mathrm{shock}}}{\mathit{c}}\mathrm{·}$(2)Bell (1978) showed that only downstream are a small fraction of particles advected from the acceleration region of a non-relativistic strong shock. This leads to a probability for the particle to remain within the acceleration region of P = 1 − V_shock/c. The resulting energy spectrum power-law is then given by $\mathrm{d}\mathit{N}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{\propto }{\mathit{E}}^{-2}\mathrm{d}\mathit{E}\mathit{.}$(3)Looking at shells as potential sites for diffusive shock acceleration, Blandford & Ostriker (1978) estimate that at least 10% or even up to 50% of the energy of a in the typical standard Sedov solution with E = 10⁵¹erg can be used for particle acceleration. The limiting factor is rather the maximum size of the shock front, which it reaches at the end of the so-called adiabatic Sedov-Taylor expansion.

At the upper end of the cosmic-ray spectrum, particles are measured with such extreme energies that their gyroradii would exceed the dimensions of a . Additionally, the time spent by particles in the vicinity of the shock front plays a role: The faster a particle becomes, the faster it will be ejected from the acceleration zone. Hillas (1984) researched this topic extensively in the context of ultra-high-energy cosmic rays. By equating the particle’s gyroradius with the dimension L of an acceleration site, he finds a plausible upper bound to the maximum energy a particle can be accelerated to as ${\mathit{E}}_{\max }\mathrm{=}\mathit{Z}\mathrm{e}\mathit{B}{\mathit{V}}_{\mathrm{S}}\mathit{L}\mathit{.}$(4)This purely geometrical consideration is nowadays called the Hillas criterion (1984), and applies to any kind of acceleration site. For a single nucleon with Z = 1 and a supernova with a size L = 100 pc and velocity V_S = 10⁴ km s^-1 in a magnetic field of strength B = 10^-10 T, the above equation yields an upper limit to the particle energy of around 3 × 10¹⁵ eV. Cosmic rays consisting of heavier atoms with more nucleons, such as iron, could potentially be accelerated to even higher energies of around 10¹⁸ eV in the strong magnetic fields present in young (Longair 2011, Sect. 17.5.3). Beyond that regime, either different acceleration mechanisms must play a role, or other cosmic objects become the relevant acceleration sites. The Hillas plot (cf. Vukcevic & Schlickeiser 2007; Shalchi et al. 2009) shows a straight line for the solution of the Hillas criterion for a proton with E_max = 10²⁰ eV. Of the various potential accelerations sites shown in the plot, only those above the line could possibly accomplish such an acceleration.

Another open question is the origin of the initial particles. Both versions of the Fermi acceleration shown above assume that the induced particles already have relativistic velocities. Particles without a sufficiently high initial energy would not be accelerated due to ionization losses. This problem is called the injection problem. Various theories have been proposed to overcome this issue, which have been summarized in the review of Malkov & Diamond (2001). Of particular interest to the work presented here is the “thermal leakage” scenario, where the generation of Alfvén waves by accelerated particles is proposed. These waves could on one hand self-regulate the thermal particle injection and on the other hand lead to efficient pitch-angle scattering, confining accelerating particles to the shock wave and thereby increasing the efficiency of the diffusive acceleration mechanism (Malkov 1998; Winske & Quest 1988). Based on this idea, Caprioli & Spitkovsky (2014) used a Maxwell-Boltzmann distributed injection profile and could reproduce an efficient acceleration process, leading to a power spectrum with a spectral index of −1.5 in conformance with the theory under the assumption of non-relativistic particle speeds.

3. Numerical model

In the Padian code, a Monte-Carlo method is applied to calculate the transport parameters of cosmic-ray test-particles. The particle equation of motion – the Newton-Lorentz equation – is solved numerically over a defined period of time, e.g., several thousand gyro periods. This procedure is repeated for a sufficiently large number of particles, until, the individual results are combined and statistically evaluated to reach an overall result.

3.1. Magnetic field

Padian splits the magnetic field into two components: the turbulent component δB and the mean magnetic field B₀, which models the large-scale magnetic background field. In the work presented here, a uniform constant magnetic background field directed along the z-axis is applied.

The turbulent component δB is calculated for every spatial coordinate (x,y,z) and is not discretized. A detailed description of the model of magnetic turbulence applied in Padian is given in Tautz (2010a) and Tautz & Dosch (2013). The following will reiterate the most important aspects of it.

Based on the work of Shalchi & Weinhorst (2009), the power spectrum of the turbulence is modeled by a modified kappa distribution $\mathit{G}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}{\mathit{l}}_{\mathrm{0}}\mathit{k}{\mathrm{)}}^{\mathit{q}}}{{\left[\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{l}}_{\mathrm{0}}\mathit{k}{\mathrm{)}}^{\mathrm{2}}\right]}^{\mathrm{(}\mathit{s}\mathrm{+}\mathit{q}\mathrm{)}\mathit{/}\mathrm{2}}}\mathit{,}$(5)where q and s define the energy range and inertial range spectral indices, respectively. The transition between these two ranges occurs at the so-called bend-over scale l₀, which will be used for normalization purposes as discussed below. For q = 0 and s = 5/3 one obtains the Kolmogorov spectrum of turbulence with a flat energy range. An important alternative to the above spectrum has been given by Iroshnikov (1964) and Kraichnan (2003), who extended the original theory to fluids carrying a magnetic field: In their model, which is nowadays called Iroshnikov-Kraichnan turbulence, they assume locality of interactions, weak interactions, and isotropy. For turbulence in a plasma influenced by a uniform magnetic field they then arrive at an energy spectrum of the form $\mathit{E}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\propto }\mathrm{(}\mathit{ϵ}{\mathit{V}}_{\mathrm{A}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}$(6)where ${\mathit{V}}_{\mathrm{A}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$ denotes the Alfvén velocity and where “epsilon” is the rate of energy dissipation. Numerous experiments successfully reproduced either form of the energy spectrum when investigating turbulent systems (Bruno & Carbone 2013). But both of these theories, and many others, assume an isotropic turbulence which is contradicted by many observational results, in particular in the solar wind. More recently, Goldreich & Sridhar (1997) presented another model that takes the anisotropy in the turbulence of astrophysical plasmas into account. However, no model has been found which can explain all experimental data nor predict the results of numerical simulations. A recent review on this matter is given by Schekochihin & Cowley (2007).

Padian furthermore follows the ideas of Giacalone & Jokipii (1999) and generates the turbulent magnetic field by superposing multiple plane waves with random phase angles and directions of propagation. According to Batchelor (1953), one can model isotropic, spatially homogeneous turbulence this way, when using a large number of waves modes. For an individual mode, the magnetic field can then be calculated with (Tautz & Dosch 2013)$\mathit{\delta }{{B}}_{\mathit{n}}\mathrm{=}\mathit{\zeta ̂}\mathrm{\perp }\mathit{A}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\hspace{0.17em}\cos \left({\mathit{k}}_{\mathit{n}}\mathit{\zeta }\mathrm{+}{\mathit{\psi }}_{\mathit{n}}\right)\mathit{,}$(7)where $\mathit{A}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\mathrm{\propto }\sqrt{\mathit{G}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\hspace{0.17em}\mathit{\Delta }{\mathit{k}}_{\mathit{n}}}$ is the amplitude function, which is defined by the energy spectrum (cf. Eq. (5)). The unit vector is perpendicular to the direction of ζ, which in turn is obtained by applying a random rotation matrix to the spatial coordinate (x,y,z) for each wave mode k_n. Additionally, each wave mode has a random phase given by ψ_n.

The real part of a summation over a number of plane wave modes N then yields the total turbulent magnetic field as $\mathit{\delta }{B}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathit{\delta }{{B}}_{\mathit{n}}\mathit{.}$(8)

3.1.1. Alfvénic turbulence

The turbulence model introduced above can be extended to support propagating Alfvén waves (e.g., Michałek & Ostrowski 1996; Tautz 2010a). This is achieved by modifying the turbulent magnetic field components of individual plane wave modes δB_n in Eq. (8)with an oscillation frequency, yielding $\mathit{\delta }{{B}}_{\mathit{n}}\mathrm{=}\mathit{\zeta ̂}\mathrm{\perp }\mathit{A}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\hspace{0.17em}\cos \left[{\mathit{k}}_{\mathit{n}}\mathit{\zeta }\mathrm{-}\mathit{\omega }\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\mathit{t}\mathrm{+}{\mathit{\psi }}_{\mathit{n}}\right]\mathit{.}$(9)Here, ω(k_n) denotes the dispersion relation of the wave, which, for Alfvén waves, is given by ω(k_n) = ± V_Ak_∥. It has to be noted that this model is only valid in absence of any shock effects, i.e. in the upstream region. Section 3.2.1 below discusses the downstream region.

Optionally, Padian can also include the effect of the electric field components of the Alfvén waves. These are perpendicular to the magnetic field components and obtained from applying the Faraday induction equation ${B}\mathrm{=}\frac{\mathrm{c}}{\mathit{\omega }}{k}\mathrm{\times }{E}$(10)by taking the polarization properties of Alfvén waves into account.

3.2. Shock wave

In order to investigate the effects of diffusive shock acceleration, Padian was extended to model a one-dimensional shock front: The ambient plasma will flow through the shock rest frame with a defined velocity toward the shock front. At the position of the shock front, a discontinuity in the magnetic field and ambient gas velocity is added according to the Rankine-Hugoniot jump conditions (see, e.g., Schlickeiser 2002, Chap. 16 and Longair 2011, Chap. 11). Note that the back reaction of the cosmic-ray particles on the shock wave will be neglected (cf. Riquelme & Spitkovsky 2010).

Padian internally works with dimensionless parameters, which simplifies the support for relativistic particles and improves the generality of the results. The time t is normalized with the relativistic gyrofrequency to τ = Ω_relt = qB₀t/ (γmc) with γ the relativistic Lorentz factor. In place of the momentum, Padian instead uses the particle’s rigidity R normalized to the turbulence bend-over scale l₀, as given by R = R_L/ℓ₀ = v/ (Ω_relℓ₀.

Using these parameters, the equation of motion is solved in the gas rest frame instead of in the observer frame. According to Faraday’s law, the flow of the ambient plasma will also induce an electric field ${{E}}_{\mathrm{ind}\mathrm{.}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{c}}{{R}}_{\mathrm{flow}}\mathrm{\times }{B}\mathit{,}$(11)where R_flow represents the flow rigidity, and B the total magnetic field. In the shock frame of reference, which is used in the simulations presented in Sect. 4, this will influence the particle rigidity even without them crossing the shock front. The effect of the flow of the ambient plasma gas on the test particles is taken into account by adapting the Newton-Lorentz equation to include the inductive electric field from Eq. (11), yielding $\frac{\mathrm{d}}{\mathrm{d}\mathit{\tau }}{R}\mathrm{=}\mathrm{(}{R}\mathrm{-}{{R}}_{\mathrm{flow}}\mathrm{)}\mathrm{\times }{B}\mathit{.}$(12)The relative rigidity R_flow follows from the Lorentz transformation into the gas rest frame. It is always parallel to the shock normal, and scales analogously to the gas velocity across the shock front. For a given value of the shocks rigidity R_S and compression ratio r, one can calculate the relative rigidity upstream and downstream of the shock front using $\begin{array}{ccc}{{R}}_{\mathrm{rel,}\mathrm{1}}& \mathrm{=}& {{R}}_{\mathrm{S}}\mathit{,}\\ {{R}}_{\mathrm{rel,}\mathrm{2}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{r}}{{R}}_{\mathrm{S}}\mathit{.}\end{array}$Upstream, the magnetic field B₁ is calculated directly according to Eq. (8). Downstream, the tangential magnetic field component B_2,t is scaled with the compression ratio as defined by the Rankine-Hugoniot equations, while the normal component B_n is continuous, leading to $\begin{array}{ccc}{\mathit{B}}_{\mathrm{2}\mathrm{,n}}& \mathrm{=}& {\mathit{B}}_{\mathrm{n}}\\ {\mathit{B}}_{\mathrm{2}\mathrm{,t}}\hspace{0.17em}& \mathrm{=}& \end{array}$This simple model allows the investigation of parallel, oblique, and perpendicular fast, strong shock waves.

3.2.1. Interaction with Alfvén waves

Alfvén waves interact with parallel shock fronts, which process was described in detail by Campeanu & Schlickeiser (1992), Schlickeiser et al. (1993), Vainio & Schlickeiser (1998, 1999), and Vainio et al. (2003). A review of this combined work is given in Schlickeiser (2002, Chap. 16.3).

The important result for the work presented here are the transmission coefficients at constant wavenumber k. They are found to be (Schlickeiser 2002, Eq. (16.3.26f)) $\begin{array}{ccc}{\mathit{T}}_{\mathit{k}}& \mathrm{=}& \frac{\sqrt{\mathit{r}}\mathrm{+}\mathrm{1}}{\mathrm{2}\sqrt{\mathit{r}}}{\left(\mathit{r}\frac{\mathit{M}\mathrm{+}{\mathit{H}}_{\mathrm{c}\mathit{,}\mathrm{1}}}{\mathit{M}\mathrm{+}\sqrt{\mathit{r}}{\mathit{H}}_{\mathrm{c}\mathit{,}\mathrm{1}}}\right)}^{\mathrm{(}\mathit{s}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}\\ {\mathit{R}}_{\mathit{k}}& \mathrm{=}& \frac{\sqrt{\mathit{r}}\mathrm{-}\mathrm{1}}{\mathrm{2}\sqrt{\mathit{r}}}{\left(\mathit{r}\frac{\mathit{M}\mathrm{+}{\mathit{H}}_{\mathrm{c}\mathit{,}\mathrm{1}}}{\mathit{M}\mathrm{-}\sqrt{\mathit{r}}{\mathit{H}}_{\mathrm{c}\mathit{,}\mathrm{1}}}\right)}^{\mathrm{(}\mathit{s}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}\mathrm{·}\end{array}$Here, r is the shock compression ratio and s is the inertial range spectral index with values 1 <s< 2. For the work presented here, it is set to 5/3 (cf. Sect. 3.1). Finally, the equation above includes the normalized cross helicity state H_c of the waves, which is given by ${\mathit{H}}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{(}\mathit{\delta }{\mathit{B}}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}\mathit{\delta }{\mathit{B}}_{\mathrm{b}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}\mathit{\delta }{\mathit{B}}_{\mathrm{f}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{\delta }{\mathit{B}}_{\mathrm{b}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(16)In Padian, only forward moving waves are injected upstream so that H_c,1 = + 1. There, the turbulent magnetic field is calculated using Eq. (9). Downstream of the shock front, the formula is adapted to take both the forward and backward traveling Alfvén waves with the corresponding transmission coefficient into account, leading to $\begin{array}{ccc}\mathit{\delta }{{B}}_{\mathrm{2}}& \mathrm{=}& \sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{m}}}\mathit{\zeta ̂}\mathrm{\perp }\mathit{A}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\mathrm{\left[}{\mathit{T}}_{\mathit{k}}\cos \mathrm{\right(}{\mathit{k}}_{\mathit{n}}\mathit{\zeta }\mathrm{-}\mathit{\omega }\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\mathit{t}\mathrm{+}{\mathit{\psi }}_{\mathit{n}}\mathrm{)}\\ & & \end{array}$(17)This model for the downstream turbulent magnetic field already leads to a significant increase in the magnetic field strength, compared to the upstream values. As such, contrary to the magnetostatic model of turbulence of Eq. (8), no artificial compression following Eq. (14b)must be applied.

Finally, it is important to state that this model is currently limited to parallel shock configurations. Vainio & Schlickeiser (1999) explicitly mention in their work that it might not be easily extensible to oblique shock configurations: “Note that our treatment holds only for unidirectional circularly polarized upstream wave fields, since variations in the magnitude of the perpendicular magnetic field component may induce motion of the shock front itself from its equilibrium position (e.g., Achterberg & Blandford 1986). The situation gets even more complicated, if the upstream wave vectors are not aligned with the shock normal, since the assumption of a planar shock is then probably violated, also. This is why we do not expect our model to be generalized to oblique shocks very easily”.

3.2.2. Spatial diffusion

In presence of a shock, the typical way of computing the spatial diffusion coefficients, κ_x,y,z, must be altered to take the ambient gas flow velocity into account. Originally, these coefficients are obtained from the particles’ mean square displacements via ${\mathit{\kappa }}_{\mathit{i,}\mathrm{orig}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{t}}_{\max }}\mathrm{⟨}({\mathit{x}}_{\mathit{i}}\mathrm{\left(}{\mathit{t}}_{\max }\mathrm{\right)}\mathrm{-}{\mathit{x}}_{\mathit{i}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{)}^{\mathrm{2}}\mathrm{⟩}\mathit{,} \mathit{i}\mathrm{\in }\mathrm{\left\{}\mathit{x,y,z}\mathrm{\right\}}\mathit{.}$(18)In the shock rest frame, the gas flow upstream and downstream of the shock front would also be included in the equation above. This is undesired, as one is rather interested in the diffusion perpendicular or parallel to the magnetic field lines, which, in the flux-freezing approximation, flow together with the ambient plasma. To account for that, time-dependent “running” diffusion coefficients are used (Tautz & Shalchi 2010), which are flow-corrected and read ${\mathit{\kappa }}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{t}}\mathrm{⟨}({\mathit{x}}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{\mathit{x}}_{\mathit{i}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{-}\mathit{\Delta }{\mathit{x}}_{\mathrm{flow}\mathit{,i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}{)}^{\mathrm{2}}\mathrm{⟩}\mathit{.}$(19)Here, the correction factor Δx_flow(t) is introduced, which corresponds to the displacement of the initial particle location x₀ due to the ambient gas flow. It is defined as $\mathit{\Delta }{{x}}_{\mathrm{flow}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}{{R}}_{\mathrm{flow}}\mathrm{\left(}{{x}}_{\mathrm{0}}\mathrm{\right(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathit{,}$(20)where R_flow(x₀(t)) is either the upstream or downstream flow velocity, depending on the relative position of x₀(t) to the shock front. In the Padian simulation, this factor can only be approximated by a discrete summation over a number N of time steps$\mathit{\Delta }{{x}}_{\mathrm{flow}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathit{N}}{{R}}_{\mathrm{flow}}({\mathit{x}}_{\mathrm{0}}\mathrm{(}\mathit{n}\hspace{0.17em}\mathit{\Delta t}\mathrm{)})\mathit{\Delta t}\mathit{,}$(21)where NΔt = t. It is expected that the impact on the final values of the diffusion coefficient is negligible because the flow velocity is small compared to the particle velocity.

3.3. Injection

The work presented here investigates two methods of particle injection. In the simplest case, all particles of the simulation are injected with a constant absolute rigidity value. The second injection profile is inspired by the work of Caprioli & Spitkovsky (2014) and follows the thermal leakage model outlined by Malkov & Diamond (2001). There, the particles are initially assigned an absolute non-relativistic rigidity value following an offsetted Maxwell-Boltzmann distribution of the form $\mathit{R}\mathrm{=}{\mathit{R}}_{\mathrm{offs}}\mathrm{+}{\mathit{f}}_{\mathrm{Gamma}}\left(\mathit{R}\mathrm{;}\mathrm{3}\mathit{/}\mathrm{2}\mathit{,\beta }\right)\mathit{.}$(22)Here, a C++11 library function can be used for the gamma distribution. Its probability density function in terms of the shape factor α and rate factor β is given by ${\mathit{f}}_{\mathrm{Gamma}}\mathrm{\left(}\mathit{x}\mathrm{;}\mathit{\alpha ,\beta }\mathrm{\right)}\mathrm{=}\frac{{\mathit{\beta }}^{\mathit{\alpha }}}{\mathit{\Gamma }\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}}{\mathit{x}}^{\mathit{\alpha }\mathrm{-}\mathrm{1}}{\mathrm{e}}^{\mathrm{-}\mathit{\beta x}}\mathit{,}$(23)where Γ denotes the gamma function. For a shape factor α = 3/2 this yields $\mathit{\Gamma }\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\sqrt{\mathit{\pi }}\mathit{/}\mathrm{2}$. Substituting β = 1 /k_BT, one then obtains the Maxwell-Boltzmann distribution function of energy ${\mathit{f}}_{\mathrm{MW}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\mathrm{2}\sqrt{\frac{\mathit{E}}{\mathit{\pi }}}{\left(\frac{\mathrm{1}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathit{E}\mathit{/}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathit{.}$(24)The rate parameter β is configurable in the simulations and defines the width of the distribution. The offset value R_offs introduced in Eq. (22)moves the whole distribution into a user defined area.

By default, this distribution would also produce very low rigidity values, which are uninteresting for the shock acceleration process, as their energy will not change considerably. Thus, to improve the runtime performance of the simulation, the distribution function is additionally modified to include a cutoff at R_cutoff. Rigidity values below this value are discarded and the distribution function queried again. Combined, this scheme produces an injection profile which is shown in Fig. 1.

Fig. 1 Rigidity injection profile following a Maxwell distribution according to (22)with β = (Rmax − Rmin)/2 and Roffs = Rmin, with Rmin = 10-4 and Rmax = 10-2.1. The red line corresponds to 104 particles used in a simulation, obtained from the original Maxwell distribution shown by the black dotted line with a cutoff value of Rcutoff = 10-3.

In both injection profiles, the initial rigidity direction of the individual particles is randomized. Additionally, a distinct turbulent magnetic field is created for each particle by varying the seeds of the random number generators used in the calculation of Eq. (8).

3.4. Acceleration spectrum

An energy spectrum comparable to the measured cosmic-ray energy spectrum can be obtained from the numerical results as follows. Following the ideas of Giacalone & Jokipii (1996), the number of particles with a momentum in the range [ p,p + Δp ] is given by $\mathit{N}\mathrm{\left(}\mathit{p,\Delta p}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{p}}^{\mathrm{2}}\mathit{f}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathit{\Delta p,}$(25)where f(p) represents the distribution function for the momentum p. In terms of the flux dj/ dE this can then be rewritten to yield$\mathrm{4}\mathit{\pi \Delta E}\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{E}}\mathrm{=}\mathit{N}\mathrm{\left(}\mathit{E,\Delta E}\mathrm{\right)}\mathit{v,}$(26)where v is the particle’s velocity. The above equation can be transformed to use the normalized particle rigidity R applied by Padian, instead of the energy. Starting with the energy$\mathit{E}\mathrm{=}\mathit{\gamma m}{\mathrm{c}}^{\mathrm{2}}\mathit{,}$(27)one can substitute the Lorentz factor $\mathit{\gamma }\mathrm{=}\sqrt{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{p}\mathit{/}\mathit{m}\mathrm{c}{\mathrm{)}}^{\mathrm{2}}}$ yielding$\mathit{E}\mathrm{=}\sqrt{{\left(\mathit{m}{\mathrm{c}}^{\mathrm{2}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\mathit{p}\mathrm{c}\right)}^{\mathrm{2}}}\mathit{.}$(28)Replacing the momentum with p = RΩml₀ and simplifying the expression leads to $\mathit{E}\mathrm{=}\mathit{\Omega m}{\mathit{l}}_{\mathrm{0}}\mathrm{c}\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{R}}_{\mathrm{c}}^{\mathrm{2}}}\mathit{,}$(29)with R_c = c/Ωl₀. Under the assumption of R ≫ R_c, which is valid as long as $\mathit{v}\mathit{/}\mathit{c}\mathrm{\gg }\mathrm{1}\mathit{/}\sqrt{\mathrm{2}}$, one arrives at E ∝ R. The equation for the power spectrum can thus be written as $\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{E}}\mathrm{\propto }\frac{\mathit{N}\mathrm{\left(}\mathit{R,\Delta R}\mathrm{\right)}}{\mathit{R}\hspace{0.17em}\mathit{\Delta R}}\mathrm{·}$(30)In the non-relativistic limit R ≪ R_c on the other hand, E_kin ∝ R², leading to $\frac{\mathrm{d}\mathit{j}}{\mathrm{d}{\mathit{E}}_{\mathrm{kin}}}\mathrm{\propto }\frac{\mathit{N}\mathrm{\left(}{\mathit{R}}^{\mathrm{2}}\mathit{,\Delta }{\mathit{R}}^{\mathrm{2}}\mathrm{\right)}}{{\mathit{R}}^{\mathrm{2}}\mathit{\Delta }{\mathit{R}}^{\mathrm{2}}}\mathrm{·}$(31)In the spectra plots below, logarithmic spacing is chosen by sampling the above equation with a constant ratio ΔR/R.

4. Simulation results

The following chapter presents various results obtained from running Padian under different parameter configurations. To reduce the parameter space, the following values where set in all simulations:

• The compression ratio r is set to 4.

• The magnetic mean field B₀ is set to be uniform and points along the z-axis.

• Two types of turbulent magnetic field δB generation are studied. On one hand, the magnetostatic case according to Eq. (8), or via Alfvén waves and Eq. (9). In both cases, N_m = 128 wave modes are superposed in slab geometry. A Kolmogorov energy spectrum is chosen for the magnetic turbulence, i.e. s = 5/3, with an additional energy range proportional to k^q, here with q = 0. The minimum wavenumber exponent in units of the bend-over scale is set to −5 and the maximum one to 3. Finally, the turbulent magnetic field is normalized to fulfill the condition δB/B₀ = 1, corresponding to a model of strong turbulence.

• For every particle injected into the simulation, a separate turbulent field is created with different random seed values. This ensures the independence of individual test particles when analyzing all as a statistical ensemble.

For all simulations, the particle energy spectrum at different times is obtained. This allows one to estimate: (i) the influence of the initial particle energy spectrum; (ii) the method of magnetic field turbulence generation; and (iii) the direction of the shock has on the final spectrum. Within the resulting particle energy spectra, a range is selected and its double-logarithmic data fitted linearly to find the spectral index. Note that the errors given for the spectral index come from this analysis alone. While all results are qualitatively reproducible, differently seeded Padian simulations with otherwise equal configuration parameters would yield slightly different numerical results in the order of the fit errors.

Fig. 2 Parallel magnetostatic shock with particles injected at a fixed rigidity of R(0) = 10-2. The time evolution of the power spectrum, i.e., the flux dj/ dR over total rigidity gain R/R(0), over the runtime of a Padian simulation of a strong parallel shock with 104 test particles is shown. The feature on the left side indicates that some particles lose energy. To the right, a smooth constant power law establishes itself over time.

4.1. Benchmark results

Initially, different configurations of a system with magnetostatic turbulence are investigated. The turbulent magnetic field component is generated using Eq. (8)and can easily be adapted to oblique shocks.

First, the effects of shocks on an ensemble of 10⁴ particles injected at the origin with an equal rigidity value of R = 10^-2 is investigated. All particles are assigned random initial rigidity directions. The shock rigidity is set to ${\mathit{R}}_{\mathrm{S}}\mathrm{=}\mathit{R}\mathit{/}\sqrt{\mathrm{3}}\mathrm{\approx }\mathrm{5.8}\mathrm{\times }{\mathrm{10}}^{-3}$, following the previous work by Giacalone & Jokipii (1999). The simulation runs in the shock rest frame, where the shock front lies in a plane through the origin, such that all particles are injected directly at the shock. This system is then simulated until the normalized time τ = Ω_relt = 10⁶ is reached with Ω_rel the relativistic gyrofrequency. At 100 equidistant time points, the power spectrum of all particles is computed according to Eq. (30).

For a parallel shock, where the shock front lies in the x-y-plane, the results are shown in Fig. 2. One can observe multiple effects: First, note how a fraction of the particles lose energy. This effect is a result of the relative motion of the ambient gas on the rigidity of the injected particles and unrelated to the shock acceleration. Second, all other particles are efficiently accelerated by the shock. Already after the first snapshot, a power law is exhibited for parts of the energy spectrum. Over time, this range grows as particles continue to be accelerated. A linear fit of the double-logarithmic power spectrum data in this range, as is shown in Fig. 3, yields a spectral index of about x ≈ −2.01 ± 0.03. This value is compatible with the theoretically predicated value of −2.

Fig. 3 Same as Fig. 2 but only for the last time step. A linear fit of the double logarithmic data of the final power spectrum for the last time step at τ = 106 is shown. For the data range spanned by the blue fit line, a spectral index x ≈ −2.01 ± 0.03 is obtained.

Inspection of the individual particles’ relative rigidity boost over the shock distance shows some interesting insights. First, some particles are still close to the shock front and thus within the acceleration zone, which explains the continuous rise of the total energy of the particle ensemble. Second, it appears that the bulk of particles gets advected downstream of the shock. These particles have a distance close to −R_Sτ/r, indicating that, under the flux-freezing assumption, they follow the magnetic field lines which flow with the downstream ambient plasma velocity of R_S/r.

The total kinetic energy of the system increases by a factor of close to 7 over the simulation time. The total aggregated rigidity of all particles still rises toward the end of the simulation. This is an indication that the simulation did not yet reach a steady state. By increasing the simulation time further, even higher maximum particle energies are to be expected. However, no qualitative changes to the power law are expected, just the high-energy tail will become smoother.

4.2. Maxwellian injection profile

Fig. 4 Parallel magnetostatic shock with a Maxwellian injection profile as indicated by the black dotted line. This figure shows how the initial power spectrum is modified by the interaction with the shock front. Particles get accelerated and, over time, form a smooth constant power law, close to parallel to the abscissa.

Inspired by the work of Caprioli & Spitkovsky (2014), the effect of the shock acceleration on a particle ensemble with an initially Maxwell-distributed rigidity spectrum is also investigated. The simulation was adapted such that each of the 10⁴ particles is injected with an absolute rigidity value given by Eq. (22). For the simulations below, the distribution parameters were set to β = (R_max − R_min)/2 and R_offs = R_min, with R_min = 10^-4 and R_max = 10^-2.1. The rigidity cutoff was set to R_cutoff = 10^-3. The initial direction of each particle is, as before, still random. The shock rigidity is set to R_S = 10^-2 and only a limited time period up to τ = 10⁵ is simulated. Increasing this value improves the smoothness of the power spectrum in the high-rigidity regime, at the cost of significantly longer simulation run times, but does not influence the results for the power spectra qualitatively.

For comparison purposes with the original work of Caprioli & Spitkovsky (2014), the power spectra are now obtained by plotting the flux to the power of 1.5 over the square of the rigidity. This relation corresponds to a particles’ kinetic energy dependence on the rigidity in the non-relativistic limit (cf. Sect. 3.4). This is necessary because the theory of the Maxwell-Boltzman distribution is only applicable to non-relativistic energies. Contrary to the previous injection model with relativistic constant rigidity, one now expects spectral indices of 1.5 from theory (Caprioli & Spitkovsky 2014).

Furthermore, one cannot directly calculate diffusion coefficients for this model of injection. The formula in Eq. (19)is only meaningful for particles of similar starting energy. Here, one would need to group the injected particles by energy and calculate individual coefficients, which requires a significant increase in the number of injected particles. Due to the higher computational effort this would require, the spatial diffusion coefficients are omitted below.

4.2.1. Parallel shock

For a parallel shock configuration a smooth power spectrum with a spectral index of about −1.5 quickly establishes itself, as can be seen in Fig. 4. The high-rigidity tail of the power spectrum is still in flux and expected to smooth out for longer simulation times.

It is notable how the injection profile is influenced also in the low-rigidity region below the shock rigidity. This effect might be due to the cutoff applied to the injection profile. However, while removing this may influence the shape of the power spectra, no qualitative change to the spectral index in the constant power law region is expected. Rather, one can expect that the direct injection at the shock front leads to this result. If significantly more particles are injected, in a larger space around the shock front, the low-energy tail of the power spectrum should stay rather unaffected. But, since this requires significantly higher computational effort, this assumption was not verified in the work presented here.

While Caprioli & Spitkovsky (2014, do not list a value for the spectral index they obtain for a parallel shock, the power spectrum also exhibits a region close to parallel to the abscissa, in agreement with the results presented here. Compared to the previous results obtained from the simulation of a constant injection profile discussed in Sect. 4.1, this result only differs due to the non-relativistic calculation of the power spectrum.

Fig. 5 Same as Fig. 4 but only for the last time step. A linear fit of the double logarithmic data of the spectrum for the last time step at τ = 105 yields a spectral index of x ≈ −1.49 ± 0.02.

Fig. 6 Perpendicular magnetostatic shock with a Maxwellian injection profile as indicated by the black dotted line. Again, the injection profile is modified by the shock acceleration, but the slope is not parallel to the abscissa.

Fig. 7 Same as Fig. 6 but only for the last time step. The double logarithmic data of the spectrum for the last time step at τ = 105 is scattered more strongly around the linear fit and yields a spectral index of x ≈ −1.98 ± 0.04.

4.2.2. Perpendicular shock

Fig. 8 Flow-corrected spatial diffusion coefficients κ according to Eq. (19)for the z-direction parallel to the magnetic background field B0, and the perpendicular components x and y. For a parallel shock as shown in the upper panel, the diffusion in z-direction, parallel to the shock normal, is the highest. The perpendicular diffusion coefficients in x and y direction are approximately equivalent, as expected. Throughout the shock simulation, the values of all diffusion coefficients continue to rise. A perpendicular shock configuration on the other hand produces significantly different results, as shown in the lower panel. Here, all three components differ significantly from another, with the perpendicular y-component, parallel to the shock normal, having the lowest diffusion coefficient. The parallel z-component, here perpendicular to the shock normal, is again the highest. The x-component lies in between these two. Furthermore, the diffusion coefficients κx and κy show a peak at around τ ≈ 250. The parallel diffusion coefficient on the other hand continues to rise after this time, but significantly slower than before.

For a perpendicular shock, the simulation was set up as before in Sect. 4.2.1, but the shock front now lies in the x–z-plane. The background magnetic field B₀ still points in the z-direction, and is thus transverse to the shock normal. Downstream of the shock, the transverse components of the combined magnetic field B are amplified by the shock compression ratio. In this case, the results differ strongly from those of a parallel shock. Most notably, the particles are only accelerated during a short time period, as seen in Fig. 6. This is in agreement with previous studies (Giacalone 2005; Caprioli & Spitkovsky 2014, and references therein), who showed that perpendicular shocks accelerate particles at a higher rate than parallel shocks. In particular, Giacalone (2005) found that the flux at the highest energies is dominated by particles accelerated at a perpendicular shock.

The particles diffusion coefficient in these directions is significantly larger than in the normal, y-direction, as shown in Fig. 8 for a constant injection profile. Once the particles cross the shock front, they are advected from the acceleration zone with the ambient downstream plasma flow. Because of that, a much shorter simulation time of only τ = 3 × 10³ is sufficient to capture the dynamics of the acceleration process, indicated also by the plateau reached in the total rigidity. Within this time, the diffusion coefficients leave the initial ballistic regime (see Fig. 8). After reaching a peak at around τ = 250, the transverse diffusion coefficients κ_x,y decrease over time. Due to ⟨ B_x,y ⟩ = 0 one expects these values to approach zero for even higher times. The diffusion along z, i.e. parallel to the magnetic mean field but transverse to the shock normal, also decreases its growth at this time, but still increases slightly. With ⟨B_z⟩ = B₀, one can expect κ_z to approach a steady state for higher simulation times. The strong initial increase of the diffusion coefficients, especially in the y-direction, correlates to the acceleration process: with increased energy, the Larmor radii of the particles increases and thus the diffusion coefficients increase.

The resulting power spectrum for this shock configuration is shown in Fig. 7. Fitting the more strongly scattered numerical data yields a spectral index with a value of −1.98 ± 0.04 corresponding to a constant power law. Note that this result is equivalent to that of a perpendicular shock with a constant injection profile. One can still observe efficient acceleration during the short acceleration period, with individual particle rigidities being boosted by about two orders of magnitude. Overall, the combined system increases its kinetic energy by a factor of close to 2.5.

4.2.3. Oblique shocks

Fig. 9 Oblique magnetostatic shock with a Maxwellian injection profile. The spectral indices obtained for Padian simulations of strong shocks with varying obliqueness and a Maxwell rigidity injection profile.

Fig. 10 Spatial diffusion coefficients for various shock obliqueness angles, showing a variation by multiple orders of magnitude for larger shock angles. They seem to be correlated to the spectral indices (cf. Fig. 9, with higher diffusion coefficients also yielding higher spectral indices. Note that for the calculation of diffusion coefficients, a constant injection profile has been assumed.

Between the two extreme cases of a parallel or perpendicular shock, the so called oblique shocks can be set up (Sironi & Spitkovsky 2009). For simulations using a varying level of the shock obliqueness, the resulting spectral indices are plotted in Fig. 9, showing a clear angular dependency. Again, the behavior is qualitatively equivalent to the results obtained from the simulations on a particle system with a constant rigidity injection profile.

The quasi-parallel shocks produce a power spectrum with a spectral index slightly below −1.5. For quasi-perpendicular shocks at 45° and higher, the acceleration process becomes less efficient. The spectral index decreases significantly, which is correlated to the less efficient acceleration at such shocks. For the perpendicular shocks then, a spectral index of about −2.0 is measured.

Our results are, to some extend, corroborated by the findings of Caprioli & Spitkovsky (2014, their Fig. 12), although their explanation relies on the central fact that the back-reaction of the particles on the shock structure is included. Accordingly, the structure of the self-generated turbulent magnetic fields in the upstream regime depends strongly on the shock angle. In contrast, Giacalone (2005) found that, in the high-energy end of the particle spectrum, the spectral index tends to be compatible with the predictions of diffusion shock acceleration model, independently of the shock normal angle. However, it should be noted that they obtained the high-energy tail of the spectrum with the help of a special steady-state method. At intermediate energies, the spectral index becomes steeper with an increasing shock angle, similar to our findings.

A possible explanation for our results can be found by modifying the core assumptions of diffusive shock acceleration. In particular, the escape probability might be related to the diffusive behavior of the particles and, thus, might be energy-dependent. To illustrate this further, note that, to lowest order, particles can be assumed to follow the magnetic field lines, which, on average, are given by the magnetic mean field B₀. The diffusion coefficient parallel to the direction of B₀, κ_∥ is significantly higher than that in the perpendicular direction, κ_⊥, as shown in Fig. 8. In addition, the diffusion coefficients also have a different energy dependence. In quasi-perpendicular shocks, the magnetic field lines are advected with the downstream plasma flow, and the particles will then follow this flow, with κ_y being the smallest of the three diffusion coefficient components. Additionally, the increased magnetic field strength downstream reduces the particles’ Larmor radius. Thus the acceleration process only occurs for a short time period, and just a few particles reach high energies, leading to the steep power law. The results for the diffusion coefficients for varying shock obliqueness, as shown in Fig. 10, also seem to validate this hypothesis. There, a correlation of the spectral index and the diffusion coefficients is visible: the higher diffusion coefficients in the direction parallel to the shock front enable more shock front crossings, leading to higher particle acceleration and thus a flatter power spectrum. For further discussion, the reader is referred to the investigations of shock drift acceleration (e.g., Ball & Melrose 2001; Park et al. 2013, and references therein), which underline the complex interactions of particles with quasi-perpendicular shocks.

Fig. 11 Parallel Alfvénic shock with particles injected at a fixed rigidity of R(0) = 10-2. The time evolution of the power spectrum shows shows how efficient the particles get accelerated, and that a constant power law steady state is approached.

Fig. 12 Same as Fig. 11 but only for the last time step. A (partially) linear fit of the double logarithmic data at the last time step is shown, which yields a spectral index of x ≈ −1.95 ± 0.03.

4.3. Alfvénic turbulent magnetic field

The previous results were all obtained from simulations with magnetostatic turbulent magnetic fields in slab geometry. In this section, the turbulence is generated by the Alfvén-wave model outlined in Sects. 3.1.1 and 3.2.1, and its influence on the acceleration process investigated. This model can only be applied to parallel shock configurations. In the first simulation, only the magnetic component of Alfvén plasma waves is enabled, followed by a second simulation which includes the electric field component to investigate the effect of stochastic acceleration. All other parameters for the simulations presented in the following are copied from Sect. 4.1. Most notably, particles are injected with a constant rigidity of 10^-2, and the shock rigidity is set to ${\mathrm{10}}^{-2}\mathit{/}\sqrt{\mathrm{3}}$.

Fig. 13 Parallel Alfvénic shock with particles injected at a fixed rigidity of R(0) = 10-2. Shown are the cases without (left panels) and with (right panels) turbulent electric fields. As shown in the upper panels, a strong shock efficiently accelerates the particle ensemble. The acceleration is strongest at the beginning, when all particles are close to the shock front. Due to parallel diffusion, many particles will not encounter the shock again after the initial boost, whereas others continue to get accelerated. This is apparent from the lower panels, which also show that some particles initially lose energy or never cross the shock to get accelerated.

Running Padian with such a configuration yields power spectra as shown in Fig. 11, which again indicate a clear acceleration process occurring at the shock front. Compared to the results for magnetostatic turbulence, a slightly lower value of x ≈ −1.95 ± 0.03 is measured for the spectral index of the power spectrum at the last time step (Fig. 12). Otherwise, no fundamental differences are apparent. Likewise, the growth of the total kinetic energy of the particle system, as shown in Fig. 13, resembles the previous data for magnetostatic turbulence that was presented in Sect. 4.2.1.

4.3.1. Influence of the electric field

Fig. 14 Parallel Alfvénic shock with particles injected at a fixed rigidity of R(0) = 10-2. In contrast to Figs. 11 and 12, the turbulent electric field components for an Alfvén rigidity of RA = 10-6 is also included. The figure shows how initially the acceleration process yields similar results to the simulation without electric fields. However, over time, the stochastic acceleration process influences the form of the spectra heavily, with a flux peak forming at around R/R0 ≈ 15.

When the electric field component of the Alfvén plasma waves are included in the simulation, the results are significantly modified. Figure 14 shows the power spectra obtained from such a configuration, and Fig. 13 depicts the acceleration process in terms of the kinetic energy boost of the total particle system as well as for a selection of individual particles. Compared to the simulation without electric fields, one sees a larger fraction of particles losing energy, indicated by the relatively large flux of particles with relative rigidities below 1 in the power spectrum. Furthermore, no constant power law is visible, but rather a second peak at around R/R₀ ≈ 15 grows over time. A fit of the final power spectrum (Fig. 15) to find a spectral index for comparison reasons yields a value of x ≈ −1.9 ± 0.1. The total rigidity of the particle system (cf. Fig. 13) keeps growing after the characteristic strong initial boost due to shock acceleration.

Interestingly, the total rigidity boost is of the order of one magnitude, whereas individual particle boost lie in the range of up to three orders of magnitude. The total particle ensemble gains energy by a factor of about 7, whereas individual particles encounter rigidity boosts of over more than two orders of magnitude. When the electric field component of the Alfvén model of turbulence is enabled, the total rigidity boost for short times resembles the pure shock acceleration. Over time though, the influence of the stochastic acceleration due to the electric field becomes apparent, with the total rigidity approaching a constant linear growth instead of converging toward a maximum value when the particles are advected away from the shock front (cf. Tautz et al. 2013). In the absence of a shock, the stochastic acceleration due to turbulent electric fields is well known and has been investigated theoretically in terms of momentum diffusion (e.g., Skilling 1975; Schlickeiser 1989, 1994) as well as numerically (e.g., Michałek et al. 1999; Tautz 2010b). Time-dependent turbulence in the form of magnetohydrodynamic plasma waves such as Alfvén waves can be interpreted as an acceleration mechanism comparable to a first-order Fermi process (Schlickeiser 2009).

In combination, the rigidity boosts of individual particles show the energy gain in steps, characteristic for shock acceleration, but also constant growths due to stochastic acceleration. Compared to the previous results which neglected the electric field, the effect of stochastic acceleration are thus clearly apparent.

Fig. 15 Same as Fig. 14 but only for the last time step. For comparison reasons, a region of the final power spectrum is fitted linearly. This yields a qualitative result of x ≈ −1.9 ± 0.1.

Note that Padian’s runtime performance significantly degrades for high-energy particles because the numeric integration algorithm applies a constant error tolerance, leading to an increase in integration steps with growing energy. Due to that, only a reduced simulation time up to τ = 2 × 10⁵ was simulated here, as the combined effect of both, stochastic and diffusive particle acceleration create too many high-energy particles.

5. Discussion and conclusion

In this article, the process of diffusive shock acceleration was investigated using an ensemble of independent test particles. In comparison to particle-in-cell simulations, the main advantage is that, aside from the model being easily adjustable, it is possible to correlate the acceleration process to the spatial diffusivity of the injected particles.

When the energy of the total particle ensemble is investigated, one finds power spectra being formed over time in all simulations. Independent of which model turbulence was applied, i.e. for both the magnetostatic and the Alfvén-wave model of turbulence, a spectral index of close to −2 was found for parallel shock configurations. (For non-relativistic particles, the spectral index is accordingly reduced to −1.5.) This value equals the one predicted by the theory of first-order Fermi acceleration. In both simulations, particles remain within the acceleration region over long times and the total energy of the particle ensemble keeps growing. However, the power spectrum establishes itself quickly and only the high-energy tail changes toward the end of the simulation. The diffusivity of the particles over time indicates that the ballistic regime is not left within the maximum simulation time. This is due to the continued acceleration of some particles. To exclude this effect one either needs to apply even higher simulation times, or manually exclude particles from the acceleration region.

The simulation with a full model of Alfvén waves, including their electric field components, illuminates the interplay of diffusive shock acceleration and stochastic acceleration. This has some minor influences on the resulting energy spectrum, leading to higher scattering of the data around the numerical fit corresponding to a spectral index of, in this case, −1.9 ± 0.1. Qualitatively however, no large differences are otherwise apparent.

The constructed magnetic fields of the two applied models of turbulence differ in the magnitude of the magnetic field downstream of the shock. For the magnetostatic model, the perpendicular components are scaled manually by the compression ratio, whereas the theory of Schlickeiser (2002) for Alfvén wave transmission through parallel shocks leads to a more than twice as strong downstream magnetic field. As such, with both models leading to comparable results for the power spectral indices, one might argue that, despite the increased complexity of the Alfvén-wave model of turbulence is unnecessary. However, the general behavior is different with regard to the evolution of the coupled wave-particle system. In the presence of self-generated Alfvén waves, the relative contribution to the overall acceleration of the shock and the wave field can vary and lead to a wrong estimation of shock parameters if one were to neglect the influence of these waves. This effect may therefore complicate the interpretation of astrophysical gamma radiation in terms of shock acceleration (Aharonian et al. 2004). This subject is currently under active investigation, and results will be presented in due time. In addition, note that the model used here is strongly restricted by the underlying theoretical understanding of Alfvén-wave transmission, which is so far only valid for parallel shock configurations.

In addition, the work presented here investigated the shock acceleration process on a non-relativistic particle ensemble with a Maxwell-Boltzmann distributed rigidity injection profile. In conformance with the findings of Caprioli & Spitkovsky (2014), power spectra with spectral indices of close to −1.5 are obtained, and no injection problem is apparent. Similar to the relativistic simulations, which used a constant rigidity injection profile, the acceleration process becomes less efficient for oblique shocks. For a perpendicular shock then, a spectral index of ca. −2 is found. Caprioli & Spitkovsky (2014), too, reported a significant decrease of the acceleration efficiency for quasi-perpendicular shocks with angles above 45°. Their results also show that parallel shocks are most efficient at accelerating particles, compatible to our findings. But, for quasi-parallel shocks, they also report a drop in acceleration efficiency, which cannot be seen in the results presented here, where no significant change in the spectral index is observable for oblique shocks below 45°.

A limitation of the results presented here is the slab geometry that was used in the magnetostatic model of turbulence applied throughout the simulations presented in this work. This was chosen to ensure both the applicability of the Alfvén wave transmission and reflection and the comparability with the magnetostatic model. However, a more realistic model might have a potentially large influence on the quasi-perpendicular shock simulations, because this geometry is symmetric in the directions perpendicular to the magnetic background field. Thus, in the limit of a perpendicular shock, the total turbulent magnetic field is static, even when it flows with the ambient plasma upstream or downstream of the shock. As such, particles only see the increase of the magnetic field strength, but are otherwise unaffected by the shock front. Future work certainly has to conduct the simulations presented here for different turbulence geometries, in order to investigate the effect this has on the acceleration process. Potentially, a more sophisticated geometry, such as the one proposed by Goldreich & Sridhar (1995), will increase the diffusivity along the shock normal for quasi-perpendicular shocks and thereby improve the acceleration efficiency, increasing the number of high-energy particles and yielding higher spectral indices in the energy spectrum of the particles.

In conclusion, it has been shown that, even with the limited configurations used by the test-particle simulation presented in this work, reliable results can be obtained. This lays the necessary foundation for future studies. Most notably, the effect of the turbulence geometry on the acceleration efficiency of quasi-perpendicular shocks should be investigated, as mentioned earlier. Additionally, the influence of many other simulation parameters has to be checked. These parameters include, but are not limited to, the compression ratio, the shock speed, as well as the energy spectrum of the magnetic turbulence and the relative turbulence strength. Furthermore, if an extension of the Alfvén wave transmission theory of Schlickeiser (2002) to oblique shocks can be derived, then it should be tested here as well. Also, a more sophisticated model of the background magnetic field, e.g., following the research by Uyaniker et al. (2002), Eriksen et al. (2011), could be included. Finally, one should lessen the simplifying assumptions on the shock front itself. Instead, a three-dimensional model of an expanding shock wave of finite thickness, potentially including precursor instabilities and/or Rayleigh-Taylor instabilities, could be used. In such a model, the shock front will have a varying obliqueness toward the magnetic background field, with so-far unexplored effects on the final upstream particle energy spectrum.

Acknowledgments

M.W. thanks Dieter Breitschwerdt and the ZAA for support. R.C.T. acknowledges helpful discussions with Bram Achterberg and Fathallah Alouani Bibi.

References

All Figures

	Fig. 1 Rigidity injection profile following a Maxwell distribution according to (22)with β = (Rmax − Rmin)/2 and Roffs = Rmin, with Rmin = 10-4 and Rmax = 10-2.1. The red line corresponds to 104 particles used in a simulation, obtained from the original Maxwell distribution shown by the black dotted line with a cutoff value of Rcutoff = 10-3.
In the text

Fig. 2 Parallel magnetostatic shock with particles injected at a fixed rigidity of R(0) = 10-2. The time evolution of the power spectrum, i.e., the flux dj/ dR over total rigidity gain R/R(0), over the runtime of a Padian simulation of a strong parallel shock with 104 test particles is shown. The feature on the left side indicates that some particles lose energy. To the right, a smooth constant power law establishes itself over time.

In the text

	Fig. 3 Same as Fig. 2 but only for the last time step. A linear fit of the double logarithmic data of the final power spectrum for the last time step at τ = 106 is shown. For the data range spanned by the blue fit line, a spectral index x ≈ −2.01 ± 0.03 is obtained.
In the text

	Fig. 4 Parallel magnetostatic shock with a Maxwellian injection profile as indicated by the black dotted line. This figure shows how the initial power spectrum is modified by the interaction with the shock front. Particles get accelerated and, over time, form a smooth constant power law, close to parallel to the abscissa.
In the text

	Fig. 5 Same as Fig. 4 but only for the last time step. A linear fit of the double logarithmic data of the spectrum for the last time step at τ = 105 yields a spectral index of x ≈ −1.49 ± 0.02.
In the text

	Fig. 6 Perpendicular magnetostatic shock with a Maxwellian injection profile as indicated by the black dotted line. Again, the injection profile is modified by the shock acceleration, but the slope is not parallel to the abscissa.
In the text

	Fig. 7 Same as Fig. 6 but only for the last time step. The double logarithmic data of the spectrum for the last time step at τ = 105 is scattered more strongly around the linear fit and yields a spectral index of x ≈ −1.98 ± 0.04.
In the text

Fig. 8 Flow-corrected spatial diffusion coefficients κ according to Eq. (19)for the z-direction parallel to the magnetic background field B0, and the perpendicular components x and y. For a parallel shock as shown in the upper panel, the diffusion in z-direction, parallel to the shock normal, is the highest. The perpendicular diffusion coefficients in x and y direction are approximately equivalent, as expected. Throughout the shock simulation, the values of all diffusion coefficients continue to rise. A perpendicular shock configuration on the other hand produces significantly different results, as shown in the lower panel. Here, all three components differ significantly from another, with the perpendicular y-component, parallel to the shock normal, having the lowest diffusion coefficient. The parallel z-component, here perpendicular to the shock normal, is again the highest. The x-component lies in between these two. Furthermore, the diffusion coefficients κx and κy show a peak at around τ ≈ 250. The parallel diffusion coefficient on the other hand continues to rise after this time, but significantly slower than before.

In the text

	Fig. 9 Oblique magnetostatic shock with a Maxwellian injection profile. The spectral indices obtained for Padian simulations of strong shocks with varying obliqueness and a Maxwell rigidity injection profile.
In the text

	Fig. 10 Spatial diffusion coefficients for various shock obliqueness angles, showing a variation by multiple orders of magnitude for larger shock angles. They seem to be correlated to the spectral indices (cf. Fig. 9, with higher diffusion coefficients also yielding higher spectral indices. Note that for the calculation of diffusion coefficients, a constant injection profile has been assumed.
In the text

	Fig. 11 Parallel Alfvénic shock with particles injected at a fixed rigidity of R(0) = 10-2. The time evolution of the power spectrum shows shows how efficient the particles get accelerated, and that a constant power law steady state is approached.
In the text

	Fig. 12 Same as Fig. 11 but only for the last time step. A (partially) linear fit of the double logarithmic data at the last time step is shown, which yields a spectral index of x ≈ −1.95 ± 0.03.
In the text

Fig. 13 Parallel Alfvénic shock with particles injected at a fixed rigidity of R(0) = 10-2. Shown are the cases without (left panels) and with (right panels) turbulent electric fields. As shown in the upper panels, a strong shock efficiently accelerates the particle ensemble. The acceleration is strongest at the beginning, when all particles are close to the shock front. Due to parallel diffusion, many particles will not encounter the shock again after the initial boost, whereas others continue to get accelerated. This is apparent from the lower panels, which also show that some particles initially lose energy or never cross the shock to get accelerated.

In the text

Fig. 14 Parallel Alfvénic shock with particles injected at a fixed rigidity of R(0) = 10-2. In contrast to Figs. 11 and 12, the turbulent electric field components for an Alfvén rigidity of RA = 10-6 is also included. The figure shows how initially the acceleration process yields similar results to the simulation without electric fields. However, over time, the stochastic acceleration process influences the form of the spectra heavily, with a flux peak forming at around R/R0 ≈ 15.

In the text

	Fig. 15 Same as Fig. 14 but only for the last time step. For comparison reasons, a region of the final power spectrum is fitted linearly. This yields a qualitative result of x ≈ −1.9 ± 0.1.
In the text