Simulations of ion acceleration at nonrelativistic shocks:
i) Acceleration efficiency
Abstract
We use 2D and 3D hybrid (kinetic ions – fluid electrons) simulations to investigate particle acceleration and magnetic field amplification at nonrelativistic astrophysical shocks. We show that diffusive shock acceleration operates for quasiparallel configurations (i.e., when the background magnetic field is almost aligned with the shock normal) and, for large sonic and Alfvénic Mach numbers, produces universal powerlaw spectra , where is the particle momentum. The maximum energy of accelerated ions increases with time, and it is only limited by finite box size and run time. Acceleration is mainly efficient for parallel and quasiparallel strong shocks, where 10–20% of the bulk kinetic energy can be converted to energetic particles, and becomes ineffective for quasiperpendicular shocks. Also, the generation of magnetic turbulence correlates with efficient ion acceleration, and vanishes for quasiperpendicular configurations. At very oblique shocks, ions can be accelerated via shock drift acceleration, but they only gain a factor of a few in momentum, and their maximum energy does not increase with time. These findings are consistent with the degree of polarization and the morphology of the radio and Xray synchrotron emission observed, for instance, in the remnant of SN 1006. We also discuss the transition from thermal to nonthermal particles in the ion spectrum (suprathermal region), and we identify two dynamical signatures peculiar of efficient particle acceleration, namely the formation of an upstream precursor and the alteration of standard shock jump conditions.
Subject headings:
acceleration of particles — ISM: supernova remnants — magnetic fields — shock waves1. Introduction
The problem of finding the sources of the energetic extraterrestrial nuclei that we call cosmic rays (CRs) is a longlasting one. Baade & Zwicky (1934) proposed supernova remnants (SNRs) to be responsible for such nonthermal particles, requiring about 10–30% of the kinetic energy of the SN ejecta to be channelled into accelerated particles. This energetic argument, coupled with the acceleration mechanism put forward by Fermi (1949, 1954), was at the basis of the theory of diffusive shock acceleration (DSA) developed in the late ’70s by several authors (Krymskii, 1977; Axford et al., 1978; Bell, 1978a, b; Blandford & Ostriker, 1978). The most important feature of the DSA mechanism is that the spectrum of the accelerated particles is predicted to be a universal powerlaw, with a spectral index that depends only on the shock compression ratio. Since for strong (i.e., large Mach number) shocks the compression ratio tends to the asymptotic value of , particles are expected to be accelerated with a spectrum , with the modulus of the particle momentum (see, e.g., Blandford & Ostriker, 1978).
This universality is crucial to account for the spectrum of Galactic CRs observed at Earth, which extends as a powerlaw from a few GeV up to a few times GeV for protons, and up to an additional factor for heavier nuclei with charge . The discrepancy between the spectrum predicted by DSA ( in energy for relativistic particles) and the measured one can be explained by accounting for the CR transport in the Galaxy (the residence time is inferred to be a decreasing function of the energy) and for the differential escape from the sources (e.g., Ptuskin & Zirakashvili, 2005; Caprioli et al., 2011). Moreover, radio observations of SNRs suggest that the energy distribution of accelerated electrons (in the 1–10 GeV range) is consistent with the DSA prediction (e.g., Trushkin, 1998). Also, rays from SNRs are often interpreted to be of hadronic origin (see, e.g., Caprioli, 2011; Ackermann et al., 2013), suggesting that nuclei acceleration can be efficient at SNR forward shocks (about 10% in Tycho, see Morlino & Caprioli, 2012).
Understanding particle acceleration at nonrelativistic shocks, and the conditions that may favor it, is not a problem limited to SNRs, though. Most astrophysical shocks are collisionless, i.e., energy and momentum conversion is not mediated by interparticle collisions, the mean free path for Coulomb scattering being much larger the size of the system, but rather by collective electromagnetic processes. There are many examples of collisionless shocks besides SNR ones, in a very wide range of scales: in the Solar System (e.g., the Earth’s bow shock, the interplanetary shocks triggered by coronal mass ejections, and the solarwindrelated shocks), at the radio lobes of the jets of active galaxies, and also in clusters of galaxies. These shocks span a wide range of sonic Mach numbers, magnetization, and relative inclination between the shock velocity and the unperturbed magnetic field, but they are typically associated with some form of nonthermal emission, which attests to the presence of accelerated particles.
In many cases, astrophysical shocks are also associated with magnetic fields much larger than the ambient ones, up to levels that cannot be accounted for by simple compression. For instance, Xray observations of young SNRs lead to inferred magnetic fields as large as a few hundreds G, a factor of 50–100 larger than in the interstellar medium (see, e.g., Völk et al., 2005; Parizot et al., 2006, and references therein). Such amplified fields are likely produced by accelerated ions via different plasma instabilities (see, e.g., Bykov et al., 2013, for a review).
CR spectra, acceleration efficiency and magnetic field amplification have been calculated via a twofluid approach (O’C. Drury, 1983, for a review), via Monte Carlo simulations (e.g., Ellison et al., 1990; Jones & Ellison, 1991; Ellison et al., 1995, 1996; Niemiec & Ostrowski, 2004; Vladimirov et al., 2006), as well as by solving the CR transport equation either numerically (e.g., Berezhko & Völk, 1997, 2004; Kang & Jones, 1997; Kang et al., 2002; Kang & Jones, 2006; Zirakashvili & Aharonian, 2010), or analytically (Malkov, 1997; Blasi, 2002; Amato & Blasi, 2006; Caprioli et al., 2010a; Caprioli, 2012, and references therein). These methods return very consistent results (Caprioli et al., 2010b), and can tackle the large scale dynamics of the shock, also including the CR backreaction. Nevertheless, they need to be fed with microphysical prescriptions for particle scattering and injection, and for the excitation of the magnetic turbulence.
To overcome these limitations and to selfconsistently deal with the highly nonlinear interplay between particles and electromagnetic fields, numerical kinetic simulations are needed. Such simulations of nonrelativistic collisionless shocks have been carried out either in the full particleincell (PIC) approach (see, e.g, Amano & Hoshino, 2007, 2010; Riquelme & Spitkovsky, 2011; Niemiec et al., 2012), or in the hybrid (kinetic ions–fluid electrons) approach (see, e.g., Winske, 1985; Quest, 1988; Giacalone et al., 1993; Bennett & Ellison, 1995; Winske & Omidi, 1996; Giacalone et al., 1997; Giacalone & Ellison, 2000; Giacalone, 2004; Lipatov, 2002; Gargaté & Spitkovsky, 2012; Guo & Giacalone, 2013, and references therein). A novel approach to DSA exploiting Vlasov–Fokker–Plank simulations has also been recently proposed by Bell et al. (2013); such a method does not include a selfconsistent description of the injection physics, but it represents a powerful tool for studying accelerated particles on large scales.
In this paper we want to test, by means of kinetic hybrid simulations, the spectrum of ions accelerated at nonrelativistic collisionless shocks via DSA, also investigating under which conditions ion acceleration and magnetic field amplification are effective. Compared to PIC simulations, the hybrid approach does not resolve the small electron plasma scales, allowing us to simulate more macroscopic systems without losing any information about the shock dynamics, which is mostly driven by ions. Very importantly, the selfconsistent treatment of the interplay between accelerated particles and electromagnetic fields allows us to study the correlation between ion acceleration and magnetic field amplification.
We have performed two and threedimensional (2D and 3D) simulations of nonrelativistic shocks with different strengths (parametrized by how supersonic and superAlfvénic the upstream flow is), and with different orientations of the initial magnetic field. The presented simulations cover an unprecedentedly large range of parameters that includes both weak and strong shocks, and parallel, oblique and quasiperpendicular configurations. A preliminary investigation of such a parameter space has been performed, for instance, by Giacalone et al. (1997); Gargaté & Spitkovsky (2012) (hereafter, GS12), but the present one considerably improves on the previous analyses in several respects.
In particular, we improve on the work of GS12 by allowing for much larger computational boxes and much longer time evolution, for a much larger parameter space in shock strengths and inclinations. We have also checked the results obtained in 2D simulations against 3D ones, also performing an extensive study of the dependence of the results on time resolution, which is crucial to account for the kinematics of highenergy particles, especially at high Mach numbers. All of these ingredients are necessary to properly lead to the formation of a powerlaw tail, extending for orders of magnitude in the nonthermal regime with the slope predicted by the DSA theory, which has previously not been reported in the literature. Moreover, for the first time we present evidence for CRmodified shocks, i.e., shocks whose dynamics is relevantly affected by the pressure and energy density in the accelerated particles, which leads to very distinctive features like the formation of an upstream precursor and the modification of the standard jump conditions of gaseous shocks. The improved energy conservation and longer run times compared to GS12 work enables us to refine the measurements of the accelerating efficiency of shocks as a function of flow conditions.
The paper is structured as follows. In section 2 we outline the simulation technique and demonstrate that at strong parallel shocks the momentum spectrum of accelerated ions is perfectly consistent with the prediction of DSA theory (), also showing that its extent increases with time. In section 3 we characterize the transition between thermal and nonthermal particles, in particular discussing the properties of the peculiar suprathermal region of the ion spectrum, while in section 4 we illustrate the dependence of the CR acceleration efficiency on the shock Mach number and magnetic inclination. An important result is that, while at parallel shocks a fraction as large as 10–20% of the shock ram pressure is converted into accelerated ions, very oblique shocks are inefficient in producing energetic particles and magnetic turbulence (section 5). The observational implications of these findings are also outlined, with special attention to the case of SN1006. Sections 6 and 7 provide a description of the dynamics of the CRmodified shocks and of the role played by shock drift acceleration, respectively. The results of 3D simulations are shown in section 8. We conclude in section 9.
2. Diffusive shock acceleration
The simulations in this paper are carried out with the dHybrid code (Gargaté et al., 2007), a massively parallel, Newtonian (i.e., nonrelativistic), hybrid code, in both its 2D and 3D configurations. Ions are treated kinetically and electrons as a neutralizing fluid, with a prescribed polytropic equation of state. Lengths are measured in units of the ion skin depth , where is the light speed and is the ion plasma frequency, with and the ion mass, charge and number density, respectively; time is measured in inverse cyclotron times , with the strength of the initial magnetic field. Finally, velocities are normalized to the Alfvén speed . All simulations include the three spatial components of the particle momentum, and the three components of electric and magnetic fields.
Let us consider a 2D shock with Alfvénic Mach number , where we define to be the upstream fluid velocity in the downstream reference frame. The initial magnetic field is in the same direction of (parallel shock). Ions are initialized with thermal velocity , so that their temperature reads , with the Boltzmann constant. Electrons are initially in thermal equilibrium with ions, i.e., , and their adiabatic index is chosen in such a way to reproduce the expected jump conditions at the shock, as in GS12. The sound speed is and hence the sonic Mach number reads , with the ion adiabatic index. This regime is typical of the cold interstellar medium: with cm, G and K one has kms. Except when otherwise specified, throughout the paper we indicate the shock strength simply with . The reader may refer to, e.g., Giacalone et al. (1997) for a survey of hybrid simulations of nonrelativistic shocks with different sonic and Alfvénic Mach numbers.
The computational box measures , with two cells per ion skin depth and 4 (macro)particles per cell; the timestep is chosen as . The shock evolution is followed for many ion cyclotron times, until .
The shock is produced by sending a supersonic flow against a reflecting wall (left side in figures); the interaction between the initial stream and the reflected one produces a sharp discontinuity, which propagates to the right in the figures. As a consequence, in the simulation the downstream fluid is at rest, and the kinetic energy of the upstream flow is converted into thermal energy at the shock front. It is worth mentioning that — in the literature — the shock Mach number is often quoted as measured in the shock reference frame, here indicated as , which is related to through the implicit relation
(1) 
namely, for a strong shock with .
The downstream ion energy distribution is shown in figure 1, as a function of time; we can identify three main spectral regimes. At low energies, we have a thermal distribution, wellfitted with a Maxwellian (dashed line), the temperature of which — at late times — is lower than the temperature one would predict for a strong shock without CRs. The main reason for this reduced heating of the downstream plasma is that about 20% of the energy flux is channeled into nonthermal particles (also see section 6 for more details on how CRs modify the global shock dynamics).
Always at late times, we see that the ion spectrum goes as for , where we introduced
(2) 
Such a powerlaw is in remarkable agreement with the DSA prediction for strong shocks for nonrelativistic particles, as we now discuss. In a nutshell, the DSA mechanism relies on the fact that particles diffusing back and forth across the shock repeatedly gain energy because of firstorder Fermi acceleration (Fermi, 1954). The spectrum of the accelerated particles does not depend on the details of the scattering, but only on the density jump between upstream and downstream, . For the shock compression ratio is , and the spectrum of accelerated ions is predicted to be in momentum space, with .
The energy distribution can be calculated as
(3) 
In the nonrelativistic regime , so that and ; in the relativistic limit, instead, and .
In spite of the fact that simulations of nonrelativistic collisionless shocks have been performed for many years (see section 1), this is the first time —to our knowledge— that the DSA prediction for strong shocks has been convincingly recovered in selfconsistent simulations. Previous simulations, while indeed showing evidence of suprathermal ions and, occasionally, of powerlaw distributions, have never been run for long enough, and in sufficiently large computational boxes, to unequivocally see ions accelerated through DSA over almost 3 decades in energy (figure 1). Moreover, in this work we account for Mach numbers as large as 50, while most of the previous work has been done for shocks with ; in such a regime the magnetosonic Mach number is , implying and, in turn, .
dHybrid is a nonrelativistic code, and cannot directly test the regime. However, the dependence is common to both relativistic and nonrelativistic particles; therefore, it is reasonable to expect that the obtained momentum spectrum may really be universal.
By looking at the time evolution of the nonthermal ion distribution in figure 1, one notes that the spectral slope remains almost constant in time, while the highenergy cutoff, wellfitted by an exponential , with , moves to larger and larger energies at later and later times. We stress that, since our simulation box is very large in the direction, is not artificially limited by the finite size of the simulation until at least, but it is only determined by the acceleration time. We will discuss the properties of ion diffusion and the time evolution of in a more quantitative way in a forthcoming work.
3. Suprathermal particles
While far downstream the transition between the Maxwellian and the powerlaw tail is very sharp, immediately behind the shock the situation is quite different. As shown in figure 2 (red curve), the postshock ion distribution shows a “bridge” of about one order of magnitude in energy (from a few up to ), which can be fitted with a quite steep powerlaw . This spectral feature at suprathermal energies gradually disappears when moving downstream (figure 2). Far downstream ( behind the shock) the spectrum can be clearly separated into the thermal and the nonthermal distributions.
The suprathermal region is important because it contains information about the thermalization of the incoming plasma stream. It also provides a hint that behind the shock there is a pool of particles with mildly nonthermal energies, which is crucial for understanding the problem of particle injection, namely the determination of the conditions required for some particles to take part in the DSA process.
In the past decades, many efforts have been dedicated to the study of the effects of efficient acceleration at shocks (see references in section 1), and in particular to the investigation of the dynamical backreaction of CRs on the plasma flow and on the electromagnetic field. Excellent reviews on these CRmodified shocks are, e.g., O’C. Drury (1983); Jones & Ellison (1991); Malkov & O’C. Drury (2001).
Any nonlinear model of DSA, as well as any phenomenological model that aims to explain the nonthermal emission from astrophysical shocks, requires the knowledge of the fraction of particles injected in the powerlaw tail. Such a normalization can be worked out only within a selfconsistent description of the shock transition, both in terms of electromagnetic fields and particle distribution. Kinetic simulations can provide this information, which may then be fed into models that deal only with time and lengthscales much larger than the background plasma ones. An assumption common to all of the nonlinear approaches to DSA is that the shock is considered an infinitesimally thin transition where both isotropization and particle injection occur.
A popular way of dealing with injection is represented by the thermal leakage model, which assumes that particles living sufficiently far in the tail of the downstream Maxwellian, and sufficiently close to the shock, have gyroradii large enough to recross the shock in one orbit (see, e.g., Ellison et al., 1981; Malkov, 1997; Kang et al., 2002; Blasi et al., 2005). Conversely, some authors interpreted the output of their kinetic simulations as evidence that most of the particles are injected by being reflected at the shock surface (see Guo & Giacalone, 2013, and references therein), so that the idea of “leaking” may be misleading.
The asymptotic (far downstream) ion distribution found in our simulations shows that the fraction of injected particles can be effectively parametrized by defining a threshold energy, , which marks the boundary between the thermal and nonthermal distributions (see the bottom panel of figure 2). The number of nonthermal particles can be estimated by a simple continuity argument, the CR spectrum being steep enough to be dominated by number by its lowestenergy boundary. If is the Maxwellian distribution with the downstream temperature , the total number of injected ions must be calculated at .
In the context of the thermalleakage model, the injection momentum is expressed as a multiple of the downstream thermal momentum, namely
(4) 
For a strong shock, , with . Note that here is the velocity of the upstream fluid in measured in the downstream reference frame, while the jump conditions are usually calculated in the shock frame. From figure 2 we infer , and thereby
(5) 
A necessary caveat is that, even if we account for more than 2 orders of magnitude in the nonthermal tail, the energy spectrum of the CRs accelerated, e.g., in SNR shocks, is expected to extend into the relativistic regime for several decades. Since the “real” acceleration efficiency cannot be either larger than 1, or lower than what we find (15–20%), and since for a spectrum the CR energy per decade is almost constant, a simple energetic argument suggests the normalization of a “real” CR spectrum to be smaller by a factor of , depending on , on , and on the actual slope of the CR spectrum.
The dependence on of the fraction of injected particles, , is quite strong and reads:
(6) 
An increase of in the value of determined in selfconsistent kinetic simulations should effectively compensate for the necessarily limited extension of the obtained nonthermal spectra. All these effects considered, we can conclude that kinetic simulations suggest that for parallel nonrelativistic shocks a fraction of about of the particles crossing the shock is injected into the DSA process, and that the injection momentum is . The actual mechanisms leading to ion injection will be discussed in a forthcoming paper.
4. Acceleration efficiency
A crucial question is how ion injection and acceleration depend on shock strength () and obliquity, defined by the angle between the normal to the shock and the background magnetic field .
To address this question, we have run several hybrid simulations with box size , two cells for ion skin depth and 4 particles per cell; the timestep is chosen as , in order to allow for a constant Courant number () throughout all the runs. The shock evolution is followed until in all the cases. In order to capture the potential role of the filamentation instability, large computational boxes are needed (Caprioli & Spitkovsky, 2013, hereafter CS13), and large Mach number shocks require quite small time steps to enforce energy conservation, which is not guaranteed in explicit hybrid methods (see, e.g., Giacalone et al., 1993; Gargaté et al., 2007).
We consider several obliquities, corresponding to and several Mach numbers, , spanning a large parameter space meant to account for the weak shocks in the Solar System and in the intracluster medium, but also for the strong shocks relevant for SNRs.
We checked the convergence of our results against the number of particles per cell, the box size, and the time and space resolution. This work extends the investigation of GS12 to much larger boxes, crucial not to artificially limit the growth of (because of ions escaping from the upstream boundary), and to account for the effects of the filamentation instability. Also, the simulations presented here have been run at a fixed Courant number (), rather than at a fixed as in GS12: running a survey on at a fixed time step artificially suppresses ion acceleration at strong shocks, for which energy is systematically worse conserved. Finally, GS12 did not distinguish between suprathermal and nonthermal ions, thereby the acceleration efficiency they quote cannot be directly interpreted as proper DSA efficiency. For all these reasons, we argue that figure 3 provides a physicallyconsistent picture of CR acceleration efficiency at nonrelativistic shocks, while the results of GS12 are less conclusive. The present analysis also extends the range of up to 50, the range of up to 80, and supports the results obtained in the 2D configuration with unprecedentedlylarge 3D simulations.
In figure 3 the acceleration efficiency is plotted as a function of , for different Mach numbers, as in the legend. We define as the fraction of the postshock energy density in the shape of particles with energy larger than . This definition, independent of the shock Mach number, is consistent with the fact that these particles are observed to diffuse ahead of the shock (see returning particles in the top panel in figure 2), and to increase their energy with time. Alternative definitions or cuts in energy do not appreciably change the interpretation of the results, as long as the two physical conditions above are met.
The most striking result in figure 3 is that the acceleration efficiency drops above . At quasiparallel shocks (), is almost independent of the inclination, especially for low Mach numbers (); conversely, tends to zero at high obliquities (), at all the Mach numbers. Among the presented runs, the highest acceleration efficiency is achieved for , in the parallel configuration. For lowinclination shocks, the typical efficiency is always between 5 and 15% at the time considered (). Since the maximum energy is still increasing with time, these acceleration efficiencies may not have saturated, yet. Comparisons with longer runs show that the evolution of with time is not strong, and that the value at is a good proxy of the value at later times for all the shock strengths and inclinations.
Our simulations clearly attest to a strong dependence of the ion acceleration on the shock obliquity, favoring quasiparallel shocks over quasiperpendicular ones. These results are quite different from those by GS12, where the maximum acceleration efficiency was achieved i) for shocks, and ii) for . The reasons for these discrepancies are: i) the shock precursor is more extended for small , and is easily suppressed in simulations with small box in the longitudinal direction; ii) simulations in GS12 have fixed rather that fixed Courant numbers as those presented here; therefore, in GS12 the dynamics of highenergy particles is not properly accounted for at large .
5. Magnetic field amplification
Being driven by CRinduced instabilities, the magnetic field amplification naturally correlates with the abundance of nonthermal ions. At quasiperpendicular shocks, where acceleration is inefficient, the downstream magnetic field is ordered (along ) and only times larger than the initial one, consistent with the simple compression of the transverse component of the field. At quasiparallel shocks, instead, the field is first amplified in the precursor, then compressed at the shock, and eventually further enhanced by turbulent motions triggered by the shock corrugation driven by the CR precursor (CS13). These different phenomena are illustrated in figure 4, which shows the profile of the magnitude of the magnetic field, averaged over the transverse direction, for shocks with different inclinations and the same .
The magnetic field topology is illustrated for and different shock inclinations in figures 6–8. The total magnetic field map (figure 6) matches very well the density pattern shown in figure 5. In the parallel and quasiparallel configurations, the filamentation instability (in the upstream) and the Richtmeyer–Meshkov instability (at the shock transition) produce the characteristic pattern of cavities and filaments discussed by Reville & Bell (2012) and CS13. Filaments always develop in the direction of the initial magnetic field , being produced by fieldaligned currents. For high inclinations, no magnetic modes are excited to nonlinear levels in the upstream, and the shock discontinuity results much sharper than in the quasiparallel cases. The physical explanation of this phenomenon is outlined in section 7.
Figures 7 and 8 show the parallel and transverse (out of plane) components of the magnetic field; the ordered parallel component is subtracted from to facilitate the comparison among different obliquities. As discussed by CS13, when CR acceleration is efficient, the amplified field is coiled around the upstream cavities, and turbulent in the downstream. At quasiperpendicular shocks, instead, the downstream field is mainly oriented along , except in a thin layer (a few times the gyroradius of downstream thermal ions) behind the shock.
In the downstream, the spatial profile is quite different for parallel and oblique cases. Quasiperpendicular shocks show a spike of strong field immediately behind the shock (the magnetic overshoot associated with compression due to ion gyration, see, e.g., Alsop & Arons, 1988), after which the field relaxes to its asymptotic value, according to the RankineHugoniot condition . In parallel shocks, instead, the region with large field () is much more extended, because of the amplification granted by the turbulent dynamo mechanism triggered by the corrugation of the shock surface (CS13). Eventually, the magnetic field has to drop far downstream, in order to satisfy the asymptotic RankineHugoniot condition . The transverse magnetic field, which accelerated ions build up by winding and amplifying the initial field, must relax on some damping lengthscale, which is found to be larger than the typical diffusion length of the accelerated ions.
The amount of magnetic field amplification depends also on the shock strength, as illustrated by figure 9, which shows the profile of for shocks with different inclinations and Mach numbers. Low Mach number shocks show the same trend outlined for the case, but the factor is proportionally lower, both upstream and downstream. The fact that magnetic field amplification becomes more and more effective for stronger shocks is encouraging for large amplification factors inferred in young SNRs, whose fast shocks are characterized by Mach numbers as large as .
5.1. Observational consequences: SN1006
The strength and the degree of isotropy of the downstream magnetic field have observational consequences in the intensity and in the degree of polarization of the synchrotron emission detected from astrophysical shocks.
The Galactic field is inferred to be coherent on scales of few tens of pc, comparable with the diameter of young SNRs. It is thereby legitimate to wonder whether SNR blast waves might allow to test different field obliquities in different regions of the same object.
In spite of the fact that it is usually hard to determine the actual orientation of the Galactic magnetic field at SNR locations, Reynoso et al. (2013) have exploited stateoftheart radio observations to claim evidence for the bilateral shape of the nonthermal emission of SN 1006 to be consistent with a polar cap geometry. In other words, the two lobes showing more prominent (radio and Xray) synchrotron emission are inferred to be orthogonal to the direction of the local Galactic magnetic field. Very interestingly, these polar caps also show a low degree of polarization in the radio band. Conversely, in regions with low or null nonthermal Xray emission the degree of polarization is as high as 70%, a value consistent with an ordered magnetic field, inferred to be perpendicular to the shock normal. The observations by Reynoso et al. (2013) agree very well with our finding large, turbulent magnetic fields in the downstream of parallel shocks, and simplycompressed fields at quasiperpendicular shocks.
Turbulent fields are likely the result of an efficient ion acceleration, so that the azimuthal variation of the degree of radio polarization in SN 1006 seems to support our findings. Also, the observed variation of the synchrotron emissivity is consistent with a more effective magnetic field amplification in the portions of the shell where the shock is parallel. However, the emissivity might not depend on the magnetic file strength only, but also on how efficiently electrons are accelerated.
Electrons are indeed observed to be accelerated at oblique shocks, for instance at interplanetary shocks produced by solar coronal mass ejections (Wilson et al., 2012). This empirical evidence is supported by PIC simulations (e.g., Amano & Hoshino, 2010; Riquelme & Spitkovsky, 2011; Park et al., 2012), and may rely on the fact that whistler waves, which have been demonstrated to be important for electron acceleration (Riquelme & Spitkovsky, 2011), are preferentially excited at lowMachnumber oblique shocks. Whether quasiparallel shocks can accelerate electrons is currently under investigation (also see Masters et al., 2013, for some recent results about the Saturn’s bow shock).
As it is typically the case, the only smoking gun attesting to efficient ion acceleration at parallel shocks would be represented by rays from the decay of neutral pions produced in nuclear collision between CRs and the background gas. Interestingly, TeV rays have been detected from SN 1006 exactly in correspondence of the polar caps (Acero et al., 2010). This would be a strong hint of efficient ion acceleration, if such an emission were unequivocally proved to be hadronic. Detailed multiwavelength studies are needed to support or disprove this scenario.
6. Cosmicray–modified shocks
When acceleration is efficient, and pressure and energy density in accelerated particles become comparable with the respective bulk quantities, one enters the regime of CRmodified shocks. In this section we discuss the two main dynamical consequences of efficient ion acceleration, namely the formation of a shock precursor and the modification of the shock jump conditions that we observe in our simulations.
6.1. Upstream precursor
The most prominent signature of a CRmodified shock is the formation of a precursor, in which the upstream fluid is slowed down because of the pressure in CRs diffusing back and forth across the shocks. Indicating with the fluid velocity in the shock reference frame, in the stationary limit one has from mass flux conservation:
(7) 
so that deceleration also leads to compression, and hence to adiabatic heating. However, if some amplified magnetic modes were dissipated, for instance because of nonlinear Landau damping, there should be an additional source of heating, usually referred to as turbulent (or Alfvén) heating (see, e.g., McKenzie & Völk, 1982; Berezhko & Ellison, 1999; Amato & Blasi, 2006).
Let us consider a 2D simulation whose box size measures , with two cells per ion skin depth and 4 macroparticles per cell. The large transverse size is chosen in order to fully account for the effects of the filamentation instability, which can develop only if the box is larger than the gyroradius of the most energetic ions, at any time. The formation of rarefied cavities and dense filaments associated with the streaming of accelerated ions may significantly alter the results one would obtain with 1D simulations (CS13).
In figure 10 we show the ion distribution as a function of position for parallel shock. The flow is initialized at the upstream boundary with bulk velocity and thermal speed ; at the ion distribution retains basically the same properties. By looking at the bottom panel of figure 10, we see that the bulk flow becomes slower closer to the shock, and its thermal spread increases. In particular, immediately ahead of the shock the bulk velocity is reduced by about 12%, and the plasma temperature is more than 5 times larger than the initial one.
As expected from momentum flux conservation, the amount of fluid braking matches very well the energy channeled into nonthermal ions, which at the considered time is about 13%. The adiabatic compression in the precursor, which in filaments reaches , is not large enough to account for such a heating: a variation of a factor of in density would correspond to an increase of % only in . Very interestingly, the increase in the plasma pressure is comparable with the increase in the magnetic pressure, , which is due to the CRinduced instability. For the same run, the magnitude of the magnetic field, averaged over the transverse direction, is immediately ahead of the shock; the local field can be even larger (up to ) due to filamentary inhomogeneities.
The observed rough equipartition between thermal and magnetic pressures throughout the precursor strongly suggests that the free energy in the CR streaming is converted in both channels, likely due to the turbulent nature of the perturbation. The net result is that the preshock Mach number is reduced by about a factor of 3, as a consequence of both the decrease of the bulk flow speed (by about 15%), and the increase of (which is about 2.4 times larger than far upstream, see figure 10). In the presence of efficient particle acceleration, magnetic field amplification and turbulent heating of the background plasma lead to weaker subshocks, and in turn to a less efficient heating of the downstream plasma (see figure 1).
The ion temperature in collisionless shocks may be directly probed by measuring the optical H emission produced in SNRs expanding into partiallyneutral media (see, e.g., Ghavamian et al., 2007; Blasi et al., 2012). In particular, Morlino et al. (2013) showed how to extract consistent information about the lengthscale and the temperature of the upstream precursor by measuring the width of the narrow Balmer line produced by chargeexchange between neutrals and ions, even in the presence of efficient CR acceleration.
6.2. Modified jump conditions
Another important effect on the shock dynamics is the modification of the purelygaseous RankineHugoniot jump conditions in the presence of nonthermal particles. Such an effect can be outlined within a simple twofluid approach: ions are separated into a thermal population, which has an adiabatic index and feels the usual entropy variation at the discontinuity, and a CR population, constituted by particles with gyroradii larger than the shock thickness, the distribution of which is approximately constant throughout the discontinuity.
In twofluid models (see, e.g., O’C. Drury, 1983; Malkov & Völk, 1996), the CR component is usually assumed to have the adiabatic index of a relativistic gas, . When acceleration is efficient, the total compression ratio (, measured between upstream and downstream infinity) may become larger than 4, since it must be calculated with an effective adiabatic index (see eq. 1). Moreover, the possible escape of accelerated particles toward upstream infinity may make the shock behave as partially radiative, thereby contributing to increase the compressibility of the downstream plasma and, in turn, leading to .
Nevertheless, nonthermal particles alter the jump conditions in a similar way even in the fully nonrelativistic case: the very presence of energetic particles that do not obey the gaseous density (and entropy) jump at the shock is sufficient to lead to a total compression ratio larger than 4. The relativistic equation of state for the CRs adopted in twofluid approaches is actually nedeed to close the system of equations; conversely, in kinetic approaches to CRmodified shocks, closure is provided by the solution of the CR transport, through different techniques (see, e.g., Caprioli et al., 2010b).
Very generally, the modification in the shock velocity profile can be worked out as a function of the CR acceleration efficiency by using mass and momentum conservation only. We have already shown how the presence of a CRinduced shock precursor leads to a weaker subshock, the compression ratio of which reads
(8) 
Here is the sonic Mach number immediately upstream of the shock, with the fluid speed calculated in the shock reference frame. is related to through the implicit relation in eq.1. The stationary momentum conservation including gas and CR pressure reads
(9) 
where 0 and 1 correspond to quantities measured at upstream infinty and immediately in front of the subshock, and the subscripts and refer to thermal gas and CRs, respectively. For simplicity, we have neglected the magnetic field pressure in eqs. 8 and 9, but it is straigthforward to include such a contribution in the calculation of the actual jump conditions (see Caprioli et al., 2009).
We assume , and normalize all the quantities to the ram pressure , also introducing , so that eq. 9 can be rewritten as
(10) 
For strong shocks , hence the total compression ratio simply reads:
(11) 
For the shock shown in figure 10, one infers and , which corresponds to (eq.1). Plugging these values in eqs. 8 and 11 returns and finally , in good agreement with the output of the simulation.
Figure 11 shows the density profiles for shocks with different inclinations: total compression ratios are typically around 4.2–4.4 for quasiparallel shocks, where the acceleration efficiency is about 10% or larger, while they are systematically lower for inefficient, quasiperpendicular shocks. The determination of the actual value of is more complicated because of the shock broadening induced by the filamentation instability; in any case, the stationary calculation outlined above, while providing an estimate of the asymptotic shock dynamics, is not adequate to describe the subshock structure, which is intrinsically timedependent in the simulation reference frame (see also in figure 11 the density spike present at quasiperpendicular shocks).
CRinduced precursors and modified jumpconditions may also produce spectral features, since ions with different energies may in principle probe different compression ratios in their diffusive motion. However, in the presented simulations, it is difficult to quantify this effect for lowenergy particles, which are nominally sensitive to , because in the suprathermal region the spectrum is steeper for other reasons (see section 2). On the other hand, highenergy particles, which should probe the total compression ratio , are expected to have a spectrum , hardly distinguishable from the linear prediction of , especially considering the effects of the exponential cutoff close to .
The possibility of producing spectra significantly different from the standard DSA prediction has intrigued scientists for many years. However, concave spectra, flatter than at the highest energies, have never been convincingly observed in SNRs, even when CR acceleration is inferred to be efficient (Caprioli, 2012). An extension of kinetic simulations to the relativistic regime would be of primary importance to shed light on the very reason of this apparent discrepancy.
7. DSA vs SDA
The ability of a shock to accelerate ions in the nonthermal regime dramatically depends on the magnetic field topology, in the following sense. At quasiparallel shocks we invariably observe a stream of ions with energies larger than a few times propagating from the shock into the upstream plasma. The interaction between this stream and the incoming flow excites magnetic perturbations, which favor the diffusion of ions with larger and larger energies back and forth across the shock. This diffusion allows ions to undergo firstorder Fermi acceleration, and eventually produces the universal DSA powerlaw distribution.
At quasiperpendicular shocks, instead, particles do not seed the upstream plasma with selfgenerated waves beyond one ion gyroradius from the shock, and DSA cannot operate spontaneously. We cannot exclude that, in the presence of a suitable preexisting magnetic turbulence, DSA might operate at oblique shocks as well (see, e.g., Giacalone et al., 1994); nevertheless, in this case the process is not guaranteed to be selfsustaining, and the achievable may depend on the largest wavelength already present in the turbulence, rather than being a dynamicallyevolving quantity.
At quasiperpendicular shocks, ions penetrate into the upstream only while gyrating around the ordered magnetic field. The ions that can probe, during one gyration, the velocity jump between upstream and downstream, are accelerated via shock drift acceleration (SDA), which allows some particles to gain energy quite rapidly. The energy gain per cycle is proportional to , which means that suprathermal particles can almost double their velocity during the first shock crossing (Chen & Armstrong, 1975). However, after a few gyrations, particles experiencing SDA are advected downstream, and none of them can achieve energies larger than a few times .
Also, parallel shocks develop transverse magnetic fields in the upstream, but the obliquity is never expected to be larger than , since in the nonlinear () regime the parallel component is found to grow as well, in turn always returning . However, since is enhanced by compression at the shock, it is hard to envision a scenario in which the downstream magnetic field does not show a prominent component transverse to the bulk fluid motion. For this reason, some particles impinging on the shock surface are always reflected back by gyrating around this transverse magnetic field, gaining energy via SDA. This process, strictly related to the injection problem, will be discussed in greater detail in a forthcoming paper.
In general, SDA is expected to produce a population of suprathermal particles, which may look like a shifted Maxwellian with temperature , which means for a compression ratio (Park et al., 2012). Such a prediction is in good agreement with the “bump” around a few times observed at quasiperpendicular shocks in figure 12, which shows the postshock ion spectra for several shock strengths and inclinations. The normalization of such a bump seems to depend on the shock strength, being more prominent for larger , but this effect is only due to the fact that it takes more time to thermalize a more supersonic flow: when integrating in the downstream at a given time, the region of incomplete thermalization behind the shock (see figure 2) turns out to be more extended for high .
If we are interested in CR acceleration over several orders of magnitude, though, we cannot rely on SDA only, and we have to wonder under which conditions particles can be promoted into the DSA regime. By looking at figure 12, one can see that, at fixed , the total number of suprathermal particles does not depend dramatically on . The biggest difference is that at oblique shocks reflected ions only undergo SDA, and their maximum energy is always , while at quasiparallel shocks some ions contribute to form a powerlaw tail, the extent of which grows with time.
We argue that the difference may be due to the geometry of the field immediately upstream of the shock. An ordered transverse field prevents suprathermal ions from penetrating into the upstream for more than one gyroradius, which dramatically reduces the excitation of magnetic turbulence upstream, for two main reasons.

First, the anisotropy of particles reflected at the shock is always in the direction of the flow, hence perpendicular to the magnetic field; in this configuration many instabilities, such as the firehose and the resonant streaming instability, are suppressed (see, e.g., Bykov et al., 2013, for a recent review).

Second, and probably most important, if accelerated ions can propagate into the upstream only for a gyroradius , the time available for any magnetic perturbation to grow (i.e., the advection time ) is drastically reduced with respect to the case of a parallel shock, in which energetic particles can propagate for a diffusion length .
Without CRinduced instabilities, the upstream magnetic field topology is left unperturbed, the shock is not corrugated, and particles cannot seed magnetic waves on larger and larger scales, that can enhance the diffusion of particles with larger and larger energies.
At parallel shocks, the upstream medium is seeded with longwavelength perturbations, which steepen in the shock layer and lead to downstream fields that are generally oblique and that can effectively reflect incoming ions. Nevertheless, reflected ions can fly away from the shock by spiraling around the upstream field, which may be quasiparallel at least for a fraction of the wave period. In other words, even if parallel shocks may often look oblique, their timedependent configuration, whose period is likely determined by the period of the upstream waves, is intrinsically different from the steady structure of quasiperpendicular shocks.
Initiallyparallel shocks can naturally realize both the downstream transverse field necessary for particle reflection and the upstream parallel field that allows ions to diffuse away from the shock. Perpendicular shocks, instead, do not trigger long wavelength modes, and cannot significantly alter their initial configuration.
In summary, ion injection into the DSA process is a natural phenomenon at parallel shocks, while at oblique shocks it may require an adhoc preexistent turbulence, with fluctuations of order , possibly on scales larger than the typical gyroradius of ions with .
8. 3D simulations
We also performed 3D simulations to ensure that the results outlined in the previous section are not an artifact of the reduced space dimensionality. This is especially important because it has been suggested that crossfield diffusion, which is possible only in 3D, may play a fundamental role in ion injection (e.g., Baring et al., 1995; Giacalone & Ellison, 2000).
We consider a shock with , and three inclinations . The box measures , with two cells for ion skin depth and 8 particles per cell; the timestep is .
The resulting downstream ion spectra are shown in figure 13. The acceleration efficiency as a function of inclination follows the same trend described for the 2D cases: it is , 3% and 1% for and , respectively.
Most interestingly, 3D simulations confirm the ineffectiveness of DSA at very oblique shocks, despite the exact account of the topology of the magnetic field lines. We do observe a large fraction of ions recrossing the shock and experiencing SDA, but no particles make it into the regime where DSA and magnetic field amplification are efficient; also, the maximum energy does not increase with time.
The 3D structure of the selfgenerated turbulence is shown in figure 14 for the three inclinations above. More precisely, the isovalue surfaces (10 levels in the interval ) and the color code (in the legend) show the value of the component of the magnetic field orthogonal to (which is initialized in the plane). The shock surface, found at , is always marked by an increase of , but the upstream and the downstream structures significantly depend on the shock obliquity.
In the parallel case, , there are several regions (also in the upstream) in which the field is amplified at levels of , and the turbulence is sustained up to large distances in the downstream; the magnetic structures have typical coherence lengths of a few tens to hundreds of .
In the quasiperpendicular case (), instead, is always less than in the upstream, and less than in the downstream. The actual magnetic field (not shown in the figure) is mainly along , and is consistently compressed in the shock transition, which is rather sharp, unlike in the parallel case. The case shows some intermediate features, exhibiting some magnetic structures in the upstream, even if less prominent than for the parallel configuration; the selfgenerated component of , however, tends to fade a few hundreds of ions skin depths behind the shock.
It is worth stressing that our results are not entirely inconsistent with the previous ones that claimed ion acceleration to be efficient at perpendicular shocks, possibly thanks to crossfield diffusion. The aforementioned results have been obtained either with Monte Carlo simulations (Baring et al., 1995), or by tracking testparticles on top of the output of hybrid simulations seeded —by hand— with largescale magnetic waves (Giacalone & Ellison, 2000). In these cases the ion scattering is either prescribed (through the specification of a meanfreepath, in Monte Carlo simulations) or artificially enhanced, because of the presence of some preexisting turbulence. Our simulations, instead, allow for a selfconsistent description of both the shock structure and the ion kinematics, and are long/large enough to potentially observe ions undergoing DSA in the selfgenerated magnetic turbulence. Throughout all of the parameter space considered in this paper, for both 2D and 3D setups, we never find quasiperpendicular shocks to be able to spontaneously generate longwavelength waves upstream of the shock, and, in turn, to accelerate particles into the DSA regime.
9. Conclusions
We used hybrid simulations of collisionless shocks to calculate under which conditions ion acceleration is efficient at nonrelativistic collisionless shocks. We can summarize our main results in the following points.

We showed, for the first time within a PIC/hybrid approach, that DSA can be very efficient (10–20% in energy) at quasiparallel shocks, producing nonthermal ion spectra with the expected universal distribution, a powerlaw in momentum, i.e., for nonrelativistic ions.

The extent of the nonthermal powerlaw tail increases with time: energetic ions are able to selfgenerate magnetic turbulence to larger and larger scales, enhancing their own scattering and the probability of recrossing the shock. The maximum energy achievable is only limited by the available computational resources (finite box size and running time).

The ion spectrum immediately behind the shock shows a peculiar shape at energies , which deviates from both a Maxwellian and the DSA powerlaw. This suprathermal region contains information about the postshock thermalization, and the processes responsible for injection of particles into the acceleration mechanism.

At quasiparallel shocks, the fraction of ions that undergo acceleration is of the total. A good parametrization of the thermal–nonthermal transition is obtained by assuming that the DSA powerlaw tail is attached to the Maxwellian at the injection momentum , where is the downstream thermal momentum.

DSA occurs only at parallel and quasiparallel shocks. Conversely, shocks show evidence of SDA, but they are ineffective in accelerating particles above . Such a limitation is intrinsic and does not depend either on the computational box, or on the reduced dimensionality of the simulations. 3D runs turn out to be very consistent with 2D simulations in this respect. These results apply only to initially “laminar” flows and ordered magnetic fields: the addition of preexisting turbulence may modify the acceleration properties of the shocks.

Magnetic field amplification goes along with efficient ion acceleration: the selfgenerated turbulence is present only at quasiparallel shocks, and the total amplification achievable in the simulation is larger for larger Mach numbers (up to upstream, for ).

The presence of accelerated ions produces an upstream precursor, in which the fluid is slowed down and compressed. The magnetic field amplification induced by CRinduced instabilities also leads to a nonadiabatic heating of the background plasma (turbulent heating), in such a way that plasma and magnetic field pressure are almost in equipartition throughout the precursor.

When acceleration is efficient, the global shock dynamics is significantly modified by the presence of nonthermal particles: in this regime of CRmodified shocks, the total compression ratio may become larger than 4, and the downstream plasma is heated up less effectively.
This paper is conceived as the first of a series aimed to investigate several aspects of ion acceleration at nonrelativistic shocks through hybrid simulations. Important points that are not addressed in the present work, but that will be covered in forthcoming publications, are the nature and the properties of the selfgenerated magnetic turbulence (also see CS13), the details of particle diffusion, and the mechanisms leading to ion injection.
We wish to thank L. Gargaté for providing a version of dHybrid, P. Blasi, E. Amato and M. Kunz for stimulating discussions, and an anonymous referee for his/her thorough comments. This research was supported by NSF grant AST0807381 and NASA grant NNX12AD01G, and facilitated by the MaxPlanck/Princeton Center for Plasma Physics. This work was also partially supported by a grant from the Simons Foundation (grant #267233 to AS), and by the NSF under Grant No. PHYS1066293 and the hospitality of the Aspen Center for Physics. Simulations were performed on the computational resources supported by the PICSciEOIT TIGRESS High Performance Computing Center and Visualization Laboratory. This research also used the resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DEAC0205CH11231, and Teragrid/XSEDE’s Ranger and Stampede under allocation No. TGAST100035.
References
 Acero et al. (2010) Acero et al., F. 2010, A&A, 516, A62
 Ackermann et al. (2013) Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807
 Alsop & Arons (1988) Alsop, D., & Arons, J. 1988, Physics of Fluids, 31, 839
 Amano & Hoshino (2007) Amano, T., & Hoshino, M. 2007, ApJ, 661, 190
 Amano & Hoshino (2010) —. 2010, Physical Review Letters, 104, 181102
 Amato & Blasi (2006) Amato, E., & Blasi, P. 2006, MNRAS, 371, 1251
 Axford et al. (1978) Axford, W. I., Leer, E., & Skadron, G. 1978, in International Cosmic Ray Conference, Vol. 11, ICRC, 132
 Baade & Zwicky (1934) Baade, W., & Zwicky, F. 1934, Proceedings of the National Academy of Science, 20, 259
 Baring et al. (1995) Baring, M. G., Ellison, D. C., & Jones, F. C. 1995, Advances in Space Research, 15, 397
 Bell (1978a) Bell, A. R. 1978a, MNRAS, 182, 147
 Bell (1978b) —. 1978b, MNRAS, 182, 443
 Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415
 Bennett & Ellison (1995) Bennett, L., & Ellison, D. C. 1995, J. Geophys. Res., 100, 3439
 Berezhko & Ellison (1999) Berezhko, E. G., & Ellison, D. C. 1999, ApJ, 526, 385
 Berezhko & Völk (1997) Berezhko, E. G., & Völk, H. J. 1997, APh, 7, 183
 Berezhko & Völk (2004) —. 2004, A&A, 427, 525
 Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJL, 221, L29
 Blasi (2002) Blasi, P. 2002, APh, 16, 429
 Blasi et al. (2005) Blasi, P., Gabici, S., & Vannoni, G. 2005, MNRAS, 361, 907
 Blasi et al. (2012) Blasi, P., Morlino, G., Bandiera, R., Amato, E., & Caprioli, D. 2012, ApJ, 755, 121
 Bykov et al. (2013) Bykov, A. M., Brandenburg, A., Malkov, M. A., & Osipov, S. M. 2013, Space Sci. Rev., arXiv:1304.7081 [astroph.HE]
 Caprioli (2011) Caprioli, D. 2011, J. Cosmology Astropart. Phys, 5, 26
 Caprioli (2012) —. 2012, J. Cosmology Astropart. Phys, 7, 38
 Caprioli et al. (2010a) Caprioli, D., Amato, E., & Blasi, P. 2010a, APh, 33, 307
 Caprioli et al. (2011) Caprioli, D., Blasi, P., & Amato, E. 2011, APh, 34, 447
 Caprioli et al. (2009) Caprioli, D., Blasi, P., Amato, E., & Vietri, M. 2009, MNRAS, 395, 895
 Caprioli et al. (2010b) Caprioli, D., Kang, H., Vladimirov, A. E., & Jones, T. W. 2010b, MNRAS, 407, 1773
 Caprioli & Spitkovsky (2013) Caprioli, D., & Spitkovsky, A. 2013, ApJ, 765, L20
 Chen & Armstrong (1975) Chen, G., & Armstrong, T. P. 1975, in International Cosmic Ray Conference, Vol. 5, International Cosmic Ray Conference, 1814
 Ellison et al. (1995) Ellison, D. C., Baring, M. G., & Jones, F. C. 1995, Ap. J., 453, 873
 Ellison et al. (1996) —. 1996, Ap. J., 473, 1029
 Ellison et al. (1981) Ellison, D. C., Jones, F. C., & Eichler, D. 1981, Journal of Geophysics Zeitschrift Geophysik, 50, 110
 Ellison et al. (1990) Ellison, D. C., Moebius, E., & Paschmann, G. 1990, Ap. J., 352, 376
 Fermi (1949) Fermi, E. 1949, Physical Review, 75, 1169
 Fermi (1954) —. 1954, Ap. J., 119, 1
 Gargaté et al. (2007) Gargaté, L., Bingham, R., Fonseca, R. A., & Silva, L. O. 2007, Comp. Phys. Commun., 176, 419
 Gargaté & Spitkovsky (2012) Gargaté, L., & Spitkovsky, A. 2012, ApJ, 744, 67
 Ghavamian et al. (2007) Ghavamian, P., Laming, J. M., & Rakowski, C. E. 2007, ApJ, 654, L69
 Giacalone (2004) Giacalone, J. 2004, ApJ, 609, 452
 Giacalone et al. (1993) Giacalone, J., Burgess, D., Schwartz, S. J., & Ellison, D. C. 1993, ApJ, 402, 550
 Giacalone et al. (1997) Giacalone, J., Burgess, D., Schwartz, S. J., Ellison, D. C., & Bennett, L. 1997, J. Geophys. Res., 102, 19789
 Giacalone & Ellison (2000) Giacalone, J., & Ellison, D. C. 2000, J. Geophys. Res., 105, 12541
 Giacalone et al. (1994) Giacalone, J., Jokipii, J. R., & Kota, J. 1994, J. Geophys. Res., 99, 19351
 Guo & Giacalone (2013) Guo, F., & Giacalone, J. 2013, ArXiv eprints, arXiv:1303.5174 [astroph.HE]
 Jones & Ellison (1991) Jones, F. C., & Ellison, D. C. 1991, Space Science Reviews, 58, 259
 Kang & Jones (1997) Kang, H., & Jones, T. W. 1997, Ap. J., 476, 875
 Kang & Jones (2006) —. 2006, APh, 25, 246
 Kang et al. (2002) Kang, H., Jones, T. W., & Gieseler, U. D. J. 2002, Ap. J., 579, 337
 Krymskii (1977) Krymskii, G. F. 1977, Akademiia Nauk SSSR Doklady, 234, 1306
 Lipatov (2002) Lipatov, A. S. 2002, The hybrid multiscale simulation technology: an introduction with application to astrophysical and laboratory plasmas (Berlin; New York: Springer, Scientific computation)
 Malkov (1997) Malkov, M. A. 1997, Ap. J., 485, 638
 Malkov & O’C. Drury (2001) Malkov, M. A., & O’C. Drury, L. 2001, Rep. Prog. Phys., 64, 429
 Malkov & Völk (1996) Malkov, M. A., & Völk, H. J. 1996, ApJ, 473, 347
 Masters et al. (2013) Masters, A., Stawarz, L., Fujimoto, M., et al. 2013, Nature Physics, 9, 164
 McKenzie & Völk (1982) McKenzie, J. F., & Völk, H. J. 1982, A&A, 116, 191
 Morlino et al. (2013) Morlino, G., Blasi, P., Bandiera, R., Amato, E., & Caprioli, D. 2013, ApJ, 768, 148
 Morlino & Caprioli (2012) Morlino, G., & Caprioli, D. 2012, A&A, 538, A81
 Niemiec & Ostrowski (2004) Niemiec, J., & Ostrowski, M. 2004, ApJ, 610, 851
 Niemiec et al. (2012) Niemiec, J., Pohl, M., Bret, A., & Wieland, V. 2012, ApJ, 759, 73
 O’C. Drury (1983) O’C. Drury, L. 1983, Reports of Progress in Physics, 46, 973
 Parizot et al. (2006) Parizot, E., Marcowith, A., Ballet, J., & Gallant, Y. A. 2006, A&A, 453, 387
 Park et al. (2012) Park, J., Workman, J. C., Blackman, E. G., Ren, C., & Siller, R. 2012, Physics of Plasmas, 19, 062904
 Ptuskin & Zirakashvili (2005) Ptuskin, V. S., & Zirakashvili, V. N. 2005, A&A, 429, 755
 Quest (1988) Quest, K. B. 1988, J. Geophys. Res., 93, 9649
 Reville & Bell (2012) Reville, B., & Bell, A. R. 2012, MNRAS, 419, 2433
 Reynoso et al. (2013) Reynoso, E. M., Hughes, J. P., & Moffett, D. A. 2013, AJ, 145, 104
 Riquelme & Spitkovsky (2011) Riquelme, M. A., & Spitkovsky, A. 2011, ApJ, 733, 63
 Trushkin (1998) Trushkin, S. A. 1998, Bulletin of the Special Astrophysics Observatory, 46, 62
 Vladimirov et al. (2006) Vladimirov, A., Ellison, D. C., & Bykov, A. 2006, ApJ, 652, 1246
 Völk et al. (2005) Völk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2005, A&A, 433, 229
 Wilson et al. (2012) Wilson, III, L. B., Koval, A., Szabo, A., et al. 2012, Geophys. Res. Lett., 39, 8109
 Winske (1985) Winske, D. 1985, Space Sci. Rev., 42, 53
 Winske & Omidi (1996) Winske, D., & Omidi, N. 1996, J. Geophys. Res., 101, 17287
 Zirakashvili & Aharonian (2010) Zirakashvili, V. N., & Aharonian, F. A. 2010, ApJ, 708, 965