Issue 
A&A
Volume 580, August 2015



Article Number  A58  
Number of page(s)  12  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201525907  
Published online  31 July 2015 
Cosmicray acceleration at collisionless astrophysical shocks using MonteCarlo simulations
Zentrum für Astronomie und Astrophysik, Technische Universität
Berlin,
Hardenbergstraße 36,
10623
Berlin,
Germany
email:
robert.c.tautz@gmail.com
Received: 16 February 2015
Accepted: 30 May 2015
Context. The diffusive shock acceleration mechanism has been widely accepted as the acceleration mechanism for galactic cosmic rays. While selfconsistent hybrid simulations have shown how powerlaw spectra are produced, detailed information on the interplay of diffusive particle motion and the turbulent electromagnetic fields responsible for repeated shock crossings are still elusive.
Aims. The framework of testparticle theory is applied to investigate the effect of diffusive shock acceleration by inspecting the obtained cosmicray energy spectra. The resulting energy spectra can be obtained this way from the particle motion and, depending on the prescribed turbulence model, the influence of stochastic acceleration through plasma waves can be studied.
Methods. A numerical MonteCarlo simulation code is extended to include collisionless shock waves. This allows one to trace the trajectories of test particle while they are being accelerated. In addition, the diffusion coefficients can be obtained directly from the particle motion, which allows for a detailed understanding of the acceleration process.
Results. The classic result of an energy spectrum with E^{2} is only reproduced for parallel shocks, while, for all other cases, the energy spectral index is reduced depending on the shock obliqueness. Qualitatively, this can be explained in terms of the diffusion coefficients in the directions that are parallel and perpendicular to the shock front.
Key words: cosmic rays / plasmas / shock waves / turbulence / acceleration of particles / diffusion
© ESO, 2015
1. Introduction
For many decades now, astrophysicists have investigated the origin of cosmic rays, the highenergy charged particle radiation that permeates space. Various experiments, both earthbound and spaceborne, measure a characteristic cosmicray energy spectrum, which follows a (moreorless) universal power law. The individual particles have kinetic energies above several MeV, and even energies up to more than 10^{20} eV have been detected experimentally, as shown by airshower experiments (e.g., Abraham et al. 2008; LetessierSelvon & 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 socalled diffusive shock acceleration mechanism or of firstorder Fermi acceleration (Axford et al. 1977; Krymsky 1977; Blandford & Ostriker 1978; Bell 1978). This theory outlines an effective acceleration mechanism for charged particles at magnetohydrodynamic 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 defacto 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 ultrahigh 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 nonrelativistic 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 testparticle 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 cosmicray energy spectra. This relatively simple testparticle 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 nonrelativistic 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 nonzero 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 nonzero, distance. Any realistic description of such a system would be highly complex. Instead, one often approximates the shock front to be onedimensional, 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} ≫ 1 are called strong shock waves. In this limit, the compression ratio r, i.e. the ratio between the upstream density ρ_{1} and the downstream density ρ_{2}, can be written as (e.g., Schlickeiser 2002, Sect. 16.2.2.1) (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 . The righthand 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^{4} 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 superAlfvé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 gammaray 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 roundtrip is given by (2)Bell (1978) showed that only downstream are a small fraction of particles advected from the acceleration region of a nonrelativistic 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 powerlaw is then given by (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^{51}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 socalled adiabatic SedovTaylor expansion.
At the upper end of the cosmicray 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 ultrahighenergy 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 (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^{4} 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^{15} eV. Cosmic rays consisting of heavier atoms with more nucleons, such as iron, could potentially be accelerated to even higher energies of around 10^{18} 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^{20} 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 selfregulate the thermal particle injection and on the other hand lead to efficient pitchangle 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 MaxwellBoltzmann 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 nonrelativistic particle speeds.
3. Numerical model
In the Padian code, a MonteCarlo method is applied to calculate the transport parameters of cosmicray testparticles. The particle equation of motion – the NewtonLorentz 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_{0}, which models the largescale magnetic background field. In the work presented here, a uniform constant magnetic background field directed along the zaxis 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 (5)where q and s define the energy range and inertial range spectral indices, respectively. The transition between these two ranges occurs at the socalled bendover scale l_{0}, 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 IroshnikovKraichnan 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 (6)where 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)(7)where 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 (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 (9)Here, ω(k_{n}) denotes the dispersion relation of the wave, which, for Alfvén waves, is given by ω(k_{n}) = ± V_{A}k_{∥}. 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 (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 onedimensional 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 RankineHugoniot jump conditions (see, e.g., Schlickeiser 2002, Chap. 16 and Longair 2011, Chap. 11). Note that the back reaction of the cosmicray 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 τ = Ω_{rel}t = qB_{0}t/ (γmc) with γ the relativistic Lorentz factor. In place of the momentum, Padian instead uses the particle’s rigidity R normalized to the turbulence bendover scale l_{0}, as given by R = R_{L}/ℓ_{0} = v/ (Ω_{rel}ℓ_{0}.
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 (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 NewtonLorentz equation to include the inductive electric field from Eq. (11), yielding (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 Upstream, the magnetic field B_{1} 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 RankineHugoniot equations, while the normal component B_{n} is continuous, leading to 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)) 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 (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 (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 (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 fluxfreezing approximation, flow together with the ambient plasma. To account for that, timedependent “running” diffusion coefficients are used (Tautz & Shalchi 2010), which are flowcorrected and read (19)Here, the correction factor Δx_{flow}(t) is introduced, which corresponds to the displacement of the initial particle location x_{0} due to the ambient gas flow. It is defined as (20)where R_{flow}(x_{0}(t)) is either the upstream or downstream flow velocity, depending on the relative position of x_{0}(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(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 nonrelativistic rigidity value following an offsetted MaxwellBoltzmann distribution of the form (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 (23)where Γ denotes the gamma function. For a shape factor α = 3/2 this yields . Substituting β = 1 /k_{B}T, one then obtains the MaxwellBoltzmann distribution function of energy (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 β = (R_{max} − R_{min})/2 and R_{offs} = R_{min}, with R_{min} = 10^{4} and R_{max} = 10^{2.1}. The red line corresponds to 10^{4} particles used in a simulation, obtained from the original Maxwell distribution shown by the black dotted line with a cutoff value of R_{cutoff} = 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 cosmicray 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 (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(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(27)one can substitute the Lorentz factor yielding(28)Replacing the momentum with p = RΩml_{0} and simplifying the expression leads to (29)with R_{c} = c/Ωl_{0}. Under the assumption of R ≫ R_{c}, which is valid as long as , one arrives at E ∝ R. The equation for the power spectrum can thus be written as (30)In the nonrelativistic limit R ≪ R_{c} on the other hand, E_{kin} ∝ R^{2}, leading to (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_{0} is set to be uniform and points along the zaxis.

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 bendover scale is set to −5 and the maximum one to 3. Finally, the turbulent magnetic field is normalized to fulfill the condition δB/B_{0} = 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 doublelogarithmic 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 10^{4} 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^{4} 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 , 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 τ = Ω_{rel}t = 10^{6} 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 xyplane, 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 doublelogarithmic 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 τ = 10^{6} 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 fluxfreezing 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 highenergy 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 Maxwelldistributed rigidity spectrum is also investigated. The simulation was adapted such that each of the 10^{4} 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^{5} is simulated. Increasing this value improves the smoothness of the power spectrum in the highrigidity 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 nonrelativistic limit (cf. Sect. 3.4). This is necessary because the theory of the MaxwellBoltzman distribution is only applicable to nonrelativistic 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 highrigidity 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 lowrigidity 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 lowenergy 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 nonrelativistic 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 τ = 10^{5} 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 τ = 10^{5} 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 Flowcorrected spatial diffusion coefficients κ according to Eq. (19)for the zdirection parallel to the magnetic background field B_{0}, and the perpendicular components x and y. For a parallel shock as shown in the upper panel, the diffusion in zdirection, 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 ycomponent, parallel to the shock normal, having the lowest diffusion coefficient. The parallel zcomponent, here perpendicular to the shock normal, is again the highest. The xcomponent 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–zplane. The background magnetic field B_{0} still points in the zdirection, 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, ydirection, 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^{3} 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_{0}, one can expect κ_{z} to approach a steady state for higher simulation times. The strong initial increase of the diffusion coefficients, especially in the ydirection, 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 quasiparallel shocks produce a power spectrum with a spectral index slightly below −1.5. For quasiperpendicular 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 backreaction of the particles on the shock structure is included. Accordingly, the structure of the selfgenerated turbulent magnetic fields in the upstream regime depends strongly on the shock angle. In contrast, Giacalone (2005) found that, in the highenergy 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 highenergy tail of the spectrum with the help of a special steadystate 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 energydependent. 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_{0}. The diffusion coefficient parallel to the direction of B_{0}, κ_{∥} 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 quasiperpendicular 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 quasiperpendicular 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énwave 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 .
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 R_{A} = 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/R_{0} ≈ 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_{0} ≈ 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). Timedependent turbulence in the form of magnetohydrodynamic plasma waves such as Alfvén waves can be interpreted as an acceleration mechanism comparable to a firstorder 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 highenergy 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^{5} was simulated here, as the combined effect of both, stochastic and diffusive particle acceleration create too many highenergy 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 particleincell 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énwave model of turbulence, a spectral index of close to −2 was found for parallel shock configurations. (For nonrelativistic particles, the spectral index is accordingly reduced to −1.5.) This value equals the one predicted by the theory of firstorder 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 highenergy 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énwave model of turbulence is unnecessary. However, the general behavior is different with regard to the evolution of the coupled waveparticle system. In the presence of selfgenerated 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énwave transmission, which is so far only valid for parallel shock configurations.
In addition, the work presented here investigated the shock acceleration process on a nonrelativistic particle ensemble with a MaxwellBoltzmann 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 quasiperpendicular shocks with angles above 45°. Their results also show that parallel shocks are most efficient at accelerating particles, compatible to our findings. But, for quasiparallel 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 quasiperpendicular 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 quasiperpendicular shocks and thereby improve the acceleration efficiency, increasing the number of highenergy 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 testparticle 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 quasiperpendicular 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 threedimensional model of an expanding shock wave of finite thickness, potentially including precursor instabilities and/or RayleighTaylor instabilities, could be used. In such a model, the shock front will have a varying obliqueness toward the magnetic background field, with sofar 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
 Abbasi, R., Abe, M., AbuZayyad, T., Allen, M., et al. 2014, ApJ, 790, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Science, 327, 1103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abraham, J., Abreu, P., Aglietta, M., et al. 2008, Astropart. Phys., 29, 243 [Google Scholar]
 Achterberg, A., & Blandford, R. D. 1986, MNRAS, 218, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Aharonian, F. A., Akhperjanian, A. G., Aye, K.M., et al. 2004, Nature, 432, 75 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Axford, W., Leer, E., & Skadron, G. 1977, 15th Internat. CosmicRay Conf. [Google Scholar]
 Ball, L., & Melrose, D. B. 2001, PASA, 18, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Balogh, A., & Treumann, R. A. 2013, Physics of Collisionless Shocks (Heidelberg: Springer) [Google Scholar]
 Batchelor, G. K. 1953, The Theory of Homogeneous Turbulence (Cambridge University Press) [Google Scholar]
 Bell, A. 1978, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29 [Google Scholar]
 Blasi, P. 2013, A&ARv, 21, 70 [Google Scholar]
 Bruno, R., & Carbone, V. 2013, Liv. Rev. Sol. Phys., 10, 2 [Google Scholar]
 Campeanu, A., & Schlickeiser, R. 1992, A&A, 263, 413 [NASA ADS] [Google Scholar]
 Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, K. A., Hughes, J. P., Badenes, C., et al. 2011, ApJ, 728, L28 [Google Scholar]
 Falcke, H., Gorham, P., & Protheroe, R. J. 2004, New Astron. Rev., 48, 1487 [NASA ADS] [CrossRef] [Google Scholar]
 Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J. 2005, ApJ, 624, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J., & Jokipii, J. 1996, J. Geophys. Res.: Space Phys. (1978–2012), 101, 11095 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J., & Jokipii, J. 1999, ApJ, 520, 204 [Google Scholar]
 Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Sridhar, S. 1997, ApJ, 485, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Helder, E. A., Vink, J., Bykov, A. M., et al. 2012, Space Sci. Rev., 173, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Hillas, A. M. 1984, ARA&A, 22, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Iroshnikov, P. 1964, Sov. Astron., 7, 566 [NASA ADS] [Google Scholar]
 Kang, H. 2013, J. Astron. Space Sci., 30, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Kraichnan, R. H. 2003, Phys. Fluids (1958–1988), 11, 945 [Google Scholar]
 Krymsky, G. F. 1977, Doklady Akademiya Nauk SSSR, 234, 130608 [Google Scholar]
 LetessierSelvon, A., & Stanev, T. 2011, Rev. Mod. Phys., 83, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Longair, M. S. 2011, High Energy Astrophysics (Cambridge University Press), 3 [Google Scholar]
 Malkov, M. 1998, Phys. Rev. E, 58, 4911 [NASA ADS] [CrossRef] [Google Scholar]
 Malkov, M., & Diamond, P. 2001, Phys. Plasmas, 8, 2401 [NASA ADS] [CrossRef] [Google Scholar]
 Michałek, G., & Ostrowski, M. 1996, Nonlin. Processes Geophys., 3, 66 [Google Scholar]
 Michałek, G., Ostrowski, M., & Schlickeiser, R. 1999, Sol. Phys., 184, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Park, J., Ren, C., Workman, J. C., & Blackman, E. G. 2013, ApJ, 765, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Reynolds, S. P. 2008, ARA&A, 46, 89 [Google Scholar]
 Riquelme, M. A., & Spitkovsky, A. 2010, ApJ, 717, 1054 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., & Cowley, S. C. 2007, in Magnetohydrodynamics (Springer), 85 [Google Scholar]
 Schlickeiser, R. 1989, ApJ, 336, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 1994, ApJS, 20, 929 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Springer) [Google Scholar]
 Schlickeiser, R. 2009, Mod. Phys. Lett. A, 24, 1461 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R., Campeanu, A., & Lerche, L. 1993, A&A, 276, 614 [NASA ADS] [Google Scholar]
 Shalchi, A., & Weinhorst, B. 2009, Adv. Space Res., 43, 1429 [NASA ADS] [CrossRef] [Google Scholar]
 Shalchi, A., Škoda, T., Tautz, R. C., & Schlickeiser, R. 2009, Phys. Rev. D, 80, 023012 [Google Scholar]
 Sironi, L., & Spitkovsky, A. 2009, ApJ, 698, 1523 [NASA ADS] [CrossRef] [Google Scholar]
 Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C. 2010a, Comput. Phys. Comm., 181, 71 [Google Scholar]
 Tautz, R. C. 2010b, Plasma Phys. Contr. Fusion, 52, 045016 [Google Scholar]
 Tautz, R. C., & Dosch, A. 2013, Phys. Plasmas, 20, 022302 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C., & Shalchi, A. 2010, J. Geophys. Res., 115, A03104 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C., Lerche, I., & Kruse, F. 2013, A&A, 555, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Uyaniker, B., Kothes, R., & Brunt, C. M. 2002, ApJ, 565, 1022 [NASA ADS] [CrossRef] [Google Scholar]
 Vainio, R., & Schlickeiser, R. 1998, A&A, 331, 793 [NASA ADS] [Google Scholar]
 Vainio, R., & Schlickeiser, R. 1999, A&A, 343, 303 [NASA ADS] [Google Scholar]
 Vainio, R., Virtanen, J., & Schlickeiser, R. 2003, A&A, 409, 821 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vukcevic, M., & Schlickeiser, R. 2007, A&A, 467, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Winske, D., & Quest, K. 1988, J. Geophys. Res.: Space Physics (1978–2012), 93, 9681 [CrossRef] [Google Scholar]
All Figures
Fig. 1 Rigidity injection profile following a Maxwell distribution according to (22)with β = (R_{max} − R_{min})/2 and R_{offs} = R_{min}, with R_{min} = 10^{4} and R_{max} = 10^{2.1}. The red line corresponds to 10^{4} particles used in a simulation, obtained from the original Maxwell distribution shown by the black dotted line with a cutoff value of R_{cutoff} = 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 10^{4} 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 τ = 10^{6} 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 τ = 10^{5} 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 τ = 10^{5} is scattered more strongly around the linear fit and yields a spectral index of x ≈ −1.98 ± 0.04. 

In the text 
Fig. 8 Flowcorrected spatial diffusion coefficients κ according to Eq. (19)for the zdirection parallel to the magnetic background field B_{0}, and the perpendicular components x and y. For a parallel shock as shown in the upper panel, the diffusion in zdirection, 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 ycomponent, parallel to the shock normal, having the lowest diffusion coefficient. The parallel zcomponent, here perpendicular to the shock normal, is again the highest. The xcomponent 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 R_{A} = 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/R_{0} ≈ 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.