SIMULATIONS OF ION ACCELERATION AT NON-RELATIVISTIC SHOCKS. I. ACCELERATION EFFICIENCY

The problem of finding the sources of the energetic extraterrestrial nuclei that we call cosmic rays (CRs) is a long-lasting one. Baade & Zwicky (1934) proposed supernova remnants (SNRs) to be responsible for such non-thermal particles, requiring about 10%–30% of the kinetic energy of the SN ejecta to be channeled 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 1970s by several authors (Krymskii 1977; Axford et al. 1977; Bell 1978a, 1978b; 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 power-law, 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 toward the asymptotic value of r = 4, particles are expected to be accelerated with a spectrum f(p)∝p^{−3r/(r − 1)}∝p⁻⁴, with p being 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 power-law ∝E^−2.7 from a few GeV up to a few times 10⁶ GeV for protons, and up to an additional factor Z for heavier nuclei with charge Z. The discrepancy between the spectrum predicted by DSA (∝E⁻² 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 non-relativistic shocks, and the conditions that may favor it, is not a problem limited to SNRs, however. 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 than 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 solar-wind-related 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 non-thermal 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, X-ray 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 two-fluid approach (see O'C. Drury 1983 for a review), via Monte Carlo simulations (e.g., Ellison et al. 1990, 1995, 1996; Jones & Ellison 1991; 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, 2006; Kang et al. 2002; 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 back-reaction. 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 self-consistently deal with the highly non-linear interplay between particles and electromagnetic fields, numerical kinetic simulations are needed. Such simulations of non-relativistic collisionless shocks have been carried out either in the full particle-in-cell (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 self-consistent 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 non-relativistic 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 self-consistent 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 three-dimensional (2D and 3D) simulations of non-relativistic shocks with different strengths (parameterized by how supersonic and super-Alfvé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 quasi-perpendicular configurations. A preliminary investigation of such a parameter space has been performed, for instance, by Giacalone et al. (1997) and 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 high-energy particles, especially at high Mach numbers. All of these ingredients are necessary to properly lead to the formation of a power-law tail, extending for orders of magnitude in the non-thermal 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 CR-modified 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 (∝p⁻⁴), also showing that its extent increases with time. In Section 3 we characterize the transition between thermal and non-thermal particles, in particular discussing the properties of the peculiar supra-thermal 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 CR-modified shocks and of the role played by shock drift acceleration (SDA), respectively. The results of 3D simulations are shown in Section 8. We conclude in Section 9.

The simulations in this paper are carried out with the dHybrid code (Gargaté et al. 2007), a massively parallel, Newtonian (i.e., non-relativistic), 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 c/ω_p, where c is the light speed and $\omega _p=\sqrt{4\pi n e^2/m}$ is the ion plasma frequency, with m, e, and n the ion mass, charge, and number density, respectively; time is measured in inverse cyclotron times $\omega _c^{-1}=mc/eB_0$, with B₀ the strength of the initial magnetic field. Finally, velocities are normalized to the Alfvén speed $v_A=B/\sqrt{4\pi m n}=c\omega _p/\omega _c$. 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 M_A = v_sh/v_A = 20, where we define v_sh to be the upstream fluid velocity in the downstream reference frame. The initial magnetic field B₀ = B₀x is in the same direction of v_sh = −v_shx (parallel shock). Ions are initialized with thermal velocity v_th = v_A, so that their temperature reads $T_0 = ({1}/{2})mv_{A}^2 / k_B$, with k_B the Boltzmann constant. Electrons are initially in thermal equilibrium with ions, i.e., T_e = T_i = T₀, 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 $c_s=\sqrt{2\gamma k_BT_0/m}$ and hence the sonic Mach number reads $M_s=M_A\sqrt{\gamma }$, with γ = 5/3 the ion adiabatic index. This regime is typical of the cold interstellar medium: with n = 0.1 cm⁻³, B₀ = 3 μG and T = 10⁴ K one has v_A ≈ c_s ≈ 15 km s⁻¹. Except when otherwise specified, throughout the paper we indicate the shock strength simply with M = M_A ≈ M_s. The reader may refer to, e.g., Giacalone et al. (1997) for a survey of hybrid simulations of non-relativistic shocks with different sonic and Alfvénic Mach numbers.

The computational box measures (L_x, L_y) = (10⁵, 10²)[c/ω_p]², with two cells per ion skin depth and 4 (macro)particles per cell; the time-step is chosen as $\Delta t=5\times 10^{-3}\omega _c^{-1}$. The shock evolution is followed for many ion cyclotron times, until $t=2500\omega _c^{-1}$.

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 $\tilde{M}$, which is related to M through the implicit relation

$$\begin{equation} \tilde{M}=M\left[ 1+\frac{1}{r(\tilde{M})}\right],\quad r=\frac{(\gamma +1)\tilde{M}^2}{(\gamma -1)\tilde{M}^2+2} \end{equation} \tag{ 1 }$$

namely, $\tilde{M}=5/4 M$ for a strong shock with r = 4.

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, well-fitted with a Maxwellian (dashed line), the temperature of which—at late times—is ∼20% 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 non-thermal particles (also see Section 6 for more details on how CRs modify the global shock dynamics).

Zoom In Zoom Out Reset image size

Always at late times, we see that the ion spectrum goes as E^−1.5 for E ≳ 3E_sh, where we introduced

$$\begin{equation} E_{{\rm sh}}=\frac{1}{2}mv_{{\rm sh}}^2=\frac{1}{2}mM^2v_A^2. \end{equation} \tag{ 2 }$$

Such a power-law is in remarkable agreement with the DSA prediction for strong shocks for non-relativistic 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 first-order 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, r. For M ≫ 1 the shock compression ratio is r ≃ 4, and the spectrum of accelerated ions is predicted to be ∝p^−q in momentum space, with q = 3r/(r − 1) ≃ 4.

The energy distribution f(E) can be calculated as

$$\begin{equation} 4\pi p^2 f(p){dp}= f(E){dE} \rightarrow f(E)=4\pi p^2 f(p)\frac{{dp}}{{dE}}. \end{equation} \tag{ 3 }$$

In the non-relativistic regime E = p²/2m, so that (dp/dE)∝1/p∝E^−1/2 and f(E)∝E^−1.5; in the relativistic limit, instead, E∝p and f(E)∝E⁻².

In spite of the fact that simulations of non-relativistic 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 self-consistent simulations. Previous simulations, while indeed showing evidence of supra-thermal ions and, occasionally, of power-law distributions, have never been run for long enough, and in sufficiently large computational boxes, to unequivocally see ions accelerated through DSA over almost three 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 M ≲ 10; in such a regime the magnetosonic Mach number is $v_{{\rm sh}}/\sqrt{\vphantom{A^A}\smash{{c_s^2+v_A^2}}}\lesssim 7$, implying r < 4 and, in turn, q > 4.

dHybrid is a non-relativistic code, and cannot directly test the E⁻² regime. However, the p⁻⁴ dependence is common to both relativistic and non-relativistic particles; therefore, it is reasonable to expect that the obtained momentum spectrum may really be universal.

By looking at the time evolution of the non-thermal ion distribution in Figure 1, one notes that the spectral slope remains almost constant in time, while the high-energy cut-off, well-fitted by an exponential ∝exp (− E/E_max)^τ, with τ ∼ 1.5, moves to larger and larger energies at later and later times. We stress that, since our simulation box is very large in the x-direction, E_max(t) is not artificially limited by the finite size of the simulation until $t\simeq 2000\omega _c^{-1}$ at least, but it is only determined by the acceleration time. We will discuss the properties of ion diffusion and the time evolution of E_max in a more quantitative way in a forthcoming work.

While far downstream the transition between the Maxwellian and the power-law tail is very sharp, immediately behind the shock the situation is quite different. As shown in Figure 2 (red curve), the post-shock ion distribution shows a "bridge" of about one order of magnitude in energy (from a few E_sh up to ∼10E_sh), which can be fitted with a quite steep power-law ∝E⁻³. This spectral feature at supra-thermal energies gradually disappears when moving downstream (Figure 2). Far downstream (≳ 3000c/ω_p behind the shock) the spectrum can be clearly separated into the thermal and the non-thermal distributions.

Zoom In Zoom Out Reset image size

The supra-thermal 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 non-thermal 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 back-reaction of CRs on the plasma flow and on the electromagnetic field. Excellent reviews on these CR-modified shocks are, e.g., O'C. Drury (1983), Jones & Ellison (1991), and Malkov & O'C. Drury (2001).

Any non-linear model of DSA, as well as any phenomenological model that aims to explain the non-thermal emission from astrophysical shocks, requires the knowledge of the fraction of particles injected in the power-law tail. Such a normalization can be worked out only within a self-consistent 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 length-scales much larger than the background plasma ones. An assumption common to all of the non-linear 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 parameterized by defining a threshold energy, E_inj, which marks the boundary between the thermal and non-thermal distributions (see the bottom panel of Figure 2). The number of non-thermal particles can be estimated by a simple continuity argument, the CR spectrum being steep enough to be dominated by number by its lowest-energy boundary. If f_{th, d}(E) is the Maxwellian distribution with the downstream temperature T_d, the total number of injected ions must be calculated at ≈f_th(E_inj).

In the context of the thermal-leakage model, the injection momentum $p_{{\rm inj}}=\sqrt{2 m E_{{\rm inj}}}$ is expressed as a multiple of the downstream thermal momentum, namely,

$$\begin{equation} p_{{\rm inj}}=\xi _{{\rm inj}}p_{{\rm th}}; \quad p_{{\rm th}}=\sqrt{2mk_B T_d}. \end{equation} \tag{ 4 }$$

For a strong shock, $p_{{\rm th}}=({4\gamma \sqrt{\gamma -1}}/{(\gamma +1)^2})m v_{{\rm sh}}\simeq 0.77 m v_{{\rm sh}}$, with γ = 5/3. Note that here v_sh 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 E_inj ≃ 4–5E_sh, and thereby

$$\begin{equation} \xi _{{\rm inj}}\simeq 3\hbox{--} 3.5 \,. \end{equation} \tag{ 5 }$$

A necessary caveat is that, even if we account for more than two orders of magnitude in the non-thermal 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 p⁻⁴ 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 ∼log₁₀(p_max/p_inj) ≃ 5–10, depending on v_sh/c, on p_max, and on the actual slope of the CR spectrum.

The dependence on ξ_inj of the fraction of injected particles, η, is quite strong and reads:

$$\begin{equation} \eta \simeq \frac{4\pi p_{{\rm inj}}^3 f_{{\rm th}}(p_{{\rm inj}})}{n} \propto \xi _{{\rm inj}}^3\exp {\left(-\xi _{{\rm inj}}^2 \right)}. \end{equation} \tag{ 6 }$$

An increase of ∼0.2–0.5 in the value of ξ_inj determined in self-consistent kinetic simulations should effectively compensate for the necessarily limited extension of the obtained non-thermal spectra. All these effects considered, we can conclude that kinetic simulations suggest that for parallel non-relativistic shocks a fraction of about 10⁻³–10⁻⁴ of the particles crossing the shock is injected into the DSA process, and that the injection momentum is p_inj ≃ 3–4 p_th. The actual mechanisms leading to ion injection will be discussed in a forthcoming paper.

A crucial question is how ion injection and acceleration depend on shock strength (M) and obliquity, defined by the angle ϑ between the normal to the shock and the background magnetic field B₀.

To address this question, we have run several hybrid simulations with box size (L_x, L_y) = (40000, 500)[c/ω_p]², two cells for ion skin depth and four particles per cell; the time-step is chosen as $\Delta t=(0.01/M) \omega _c^{-1}$, in order to allow for a constant Courant number (=v_shΔt/Δx) throughout all the runs. The shock evolution is followed until $t=200\omega _c^{-1}$ 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 ϑ = 0°, 20°, 30°, 45°, 50°, 60°, 80° and several Mach numbers, M = 5, 10, 30, 50, spanning a large parameter space meant to account for the weak shocks in the solar system and in the intra-cluster 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 E_max (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 ($\mathcal {C}\sim v_{{\rm sh}}\Delta t/\Delta x$), rather than at a fixed Δt as in GS12: running a survey on M 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 supra-thermal and non-thermal 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 physically consistent picture of CR acceleration efficiency at non-relativistic shocks, while the results of GS12 are less conclusive. The present analysis also extends the range of M up to 50, the range of ϑ up to 80°, and supports the results obtained in the 2D configuration with unprecedentedly large 3D simulations.

Zoom In Zoom Out Reset image size

In Figure 3 the acceleration efficiency ε_cr is plotted as a function of ϑ, for different Mach numbers, as in the legend. We define ε_cr as the fraction of the post-shock energy density in the shape of particles with energy larger than 10E_sh. 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 ϑ ∼ 45°. At quasi-parallel shocks (ϑ ≲ 40°), ε_cr is almost independent of the inclination, especially for low Mach numbers (M = 5, 10); conversely, ε_cr tends to zero at high obliquities (ϑ ≳ 60°), at all the Mach numbers. Among the presented runs, the highest acceleration efficiency is achieved for M = 50, in the parallel configuration. For low-inclination shocks, the typical efficiency is always between 5% and 15% at the time considered ($t=200\omega _c^{-1}$). 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 ε_cr with time is not strong, and that the value at $t=200\omega _c^{-1}$ 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 quasi-parallel shocks over quasi-perpendicular ones. These results are quite different from those by GS12, where the maximum acceleration efficiency was achieved (1) for 15° shocks, and (2) for M = 6. The reasons for these discrepancies are (1) the shock precursor is more extended for small ϑ, and is easily suppressed in simulations with small box in the longitudinal direction; (2) simulations in GS12 have fixed Δt rather that fixed Courant numbers as those presented here; therefore, in GS12 the dynamics of high-energy particles is not properly accounted for at large M.

Being driven by CR-induced instabilities, the magnetic field amplification naturally correlates with the abundance of non-thermal ions. At quasi-perpendicular shocks, where acceleration is inefficient, the downstream magnetic field is ordered (along $\bf y$) and only ∼4 times larger than the initial one, consistent with the simple compression of the transverse component of the field. At quasi-parallel 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 M = 30.

Zoom In Zoom Out Reset image size

The magnetic field topology is illustrated for M = 50 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 quasi-parallel configurations, the filamentation instability (in the upstream) and the Richtmyer–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 B₀, being produced by field-aligned currents. For high inclinations, no magnetic modes are excited to non-linear levels in the upstream, and the shock discontinuity results much sharper than in the quasi-parallel cases. The physical explanation of this phenomenon is outlined in Section 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 7 and 8 show the parallel and transverse (out of plane) components of the magnetic field; the ordered parallel component B₀cos ϑ is subtracted from B_x 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 quasi-perpendicular shocks, instead, the downstream field is mainly oriented along y, except in a thin layer (a few times the gyroradius of downstream thermal ions) behind the shock.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the downstream, the spatial profile is quite different for parallel and oblique cases. Quasi-perpendicular 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 Rankine–Hugoniot condition B_{⊥, d} = rB_{⊥, 0}. In parallel shocks, instead, the region with large field (B_d ∼ 10B₀) 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 Rankine–Hugoniot condition B_{⊥, d} = 0. The transverse magnetic field, which accelerated ions build up by winding and amplifying the initial field, must relax on some damping length-scale, 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 B_tot for shocks with different inclinations and Mach numbers. Low Mach number shocks show the same trend outlined for the M = 50 case, but the δB/B 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 ∼50–500.

Zoom In Zoom Out Reset image size

5.1. Observational Consequences: SN1006

When acceleration is efficient, and pressure and energy density in accelerated particles become comparable with the respective bulk quantities, one enters the regime of CR-modified 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 CR-modified 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 $\tilde{u}$ the fluid velocity in the shock reference frame, in the stationary limit one has from mass flux conservation:

$$\begin{equation} n(x) \tilde{u}(x) = {\rm const}, \end{equation} \tag{ 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 non-linear 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 (L_x, L_y) = (4 × 10⁴, 10³)[c/ω_p]², with two cells per ion skin depth and four 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 p_x distribution as a function of position for M = 20 parallel shock. The flow is initialized at the upstream boundary with bulk velocity −20v_A and thermal speed v_th = v_A; at x ≃ 10⁴c/ω_p 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 $T\propto v_{{\rm th}}^2$ is more than five times larger than the initial one.

Zoom In Zoom Out Reset image size

As expected from momentum flux conservation, the amount of fluid braking matches very well the energy channeled into non-thermal ions, which at the considered time is about 13%. The adiabatic compression in the precursor, which in filaments reaches δn/n ∼ 2, is not large enough to account for such a heating: a variation of a factor of ∼2 in density would correspond to an increase of ∼60% only in T∝P/n∝n^2/3. Very interestingly, the increase in the plasma pressure is comparable with the increase in the magnetic pressure, P_B = B²/8π, which is due to the CR-induced instability. For the same run, the magnitude of the magnetic field, averaged over the transverse direction, is B_tot ≈ 2–3B₀ immediately ahead of the shock; the local field can be even larger (up to ∼6B₀) 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 pre-shock Mach number is reduced by about a factor of three, as a consequence of both the decrease of the bulk flow speed (by about 15%), and the increase of v_th (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 partially neutral 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 length-scale and the temperature of the upstream precursor by measuring the width of the narrow Balmer line produced by charge-exchange 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 purely gaseous Rankine–Hugoniot jump conditions in the presence of non-thermal particles. Such an effect can be outlined within a simple two-fluid 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 two-fluid 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, γ = 4/3. When acceleration is efficient, the total compression ratio (r_tot, measured between upstream and downstream infinity) may become larger than 4, since it must be calculated with an effective adiabatic index 4/3 ⩽ γ ⩽ 5/3 (see Equation (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 r_tot ≫ 4.

Nevertheless, non-thermal particles alter the jump conditions in a similar way even in the fully non-relativistic 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 two-fluid approaches is actually needed to close the system of equations; conversely, in kinetic approaches to CR-modified 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 CR-induced shock precursor leads to a weaker subshock, the compression ratio of which reads

$$\begin{equation} r_{{\rm sub}}=\frac{(\gamma +1)\tilde{M}_1^2}{(\gamma -1)\tilde{M}_1^2+2}. \end{equation} \tag{ 8 }$$

Here $\tilde{M_1}$ is the sonic Mach number immediately upstream of the shock, with the fluid speed calculated in the shock reference frame. $\tilde{M}$ is related to M_s through the implicit relation in Equation (1). The stationary momentum conservation including gas and CR pressure reads

$$\begin{equation} \rho _0\tilde{u}_0^2+P_{g,0}+P_{{\rm cr},0}=\rho _1\tilde{u}_1^2+P_{g,1}+P_{{\rm cr},1}, \end{equation} \tag{ 9 }$$

where 0 and 1 correspond to quantities measured at upstream infinity and immediately in front of the subshock, and the subscripts g and cr refer to thermal gas and CRs, respectively. For simplicity, we have neglected the magnetic field pressure in Equations (8) and (9), but it is straightforward to include such a contribution in the calculation of the actual jump conditions (see Caprioli et al. 2009).

We assume P_{cr, 0} = 0, and normalize all the quantities to the ram pressure $\rho _0\tilde{u}_0^2$, also introducing $\Xi =P_{{\rm cr},1}/\rho _0\tilde{u}_0^2$, so that Equation (9) can be rewritten as

$$\begin{equation} 1+\frac{1}{\gamma \tilde{M}_0^2}=\frac{r_{{\rm sub}}}{r_{{\rm tot}}} \left[1+\frac{1}{\gamma \tilde{M}_1^2} \right]+\Xi \,. \end{equation} \tag{ 10 }$$

For strong shocks $\tilde{M}_0^2\gg 1$, hence the total compression ratio simply reads

$$\begin{equation} r_{{\rm tot}}\simeq r_{{\rm sub}}\frac{1+1/\left(\gamma \tilde{M}_1^2 \right)}{1-\Xi }. \end{equation} \tag{ 11 }$$

For the shock shown in Figure 10, one infers Ξ ≈ 0.12 and M_{s, 1} ≈ 5.7, which corresponds to $\tilde{M}_1\approx 5.6$ (Equation (1)). Plugging these values in Equations (8) and (11) returns r_sub ≈ 3.65 and finally r_tot ≈ 4.23, in good agreement with the output of the simulation.

Figure 11 shows the density profiles for M = 30 shocks with different inclinations: total compression ratios are typically around 4.2–4.4 for quasi-parallel shocks, where the acceleration efficiency is about 10% or larger, while they are systematically lower for inefficient, quasi-perpendicular shocks. The determination of the actual value of r_sub 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 time-dependent in the simulation reference frame (see also in Figure 11 the density spike present at quasi-perpendicular shocks).

Zoom In Zoom Out Reset image size

CR-induced precursors and modified jump-conditions 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 low-energy particles, which are nominally sensitive to r_sub, because in the supra-thermal region the spectrum is steeper for other reasons (see Section 2). On the other hand, high-energy particles, which should probe the total compression ratio r_tot ≈ 4.3, are expected to have a spectrum ∝E^−1.43, hardly distinguishable from the linear prediction of E^−1.5, especially considering the effects of the exponential cut-off close to E_max.

The possibility of producing spectra significantly different from the standard DSA prediction has intrigued scientists for many years. However, concave spectra, flatter than p⁻⁴ 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.

The ability of a shock to accelerate ions in the non-thermal regime dramatically depends on the magnetic field topology, in the following sense. At quasi-parallel shocks we invariably observe a stream of ions with energies larger than a few times E_sh 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 first-order Fermi acceleration, and eventually produces the universal DSA power-law distribution.

At quasi-perpendicular shocks, instead, particles do not seed the upstream plasma with self-generated waves beyond one ion gyroradius from the shock, and DSA cannot operate spontaneously. We cannot exclude that, in the presence of a suitable pre-existing 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 self-sustaining, and the achievable E_max may depend on the largest wavelength already present in the turbulence, rather than being a dynamically evolving quantity.

At quasi-perpendicular 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 SDA, which allows some particles to gain energy quite rapidly. The energy gain per cycle is proportional to ≈v/v_sh, which means that supra-thermal 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 E_sh.

Also, parallel shocks develop transverse magnetic fields in the upstream, but the obliquity is never expected to be larger than ∼45°, since in the non-linear (δB/B₀ ≳ 1) regime the parallel component is found to grow as well, in turn always returning B_⊥/B_∥ ∼ 1. However, since B_⊥ 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 supra-thermal particles, which may look like a shifted Maxwellian with temperature T_SDA = T_dexp [8r(r − 1)/(r + 1)³], which means T_SDA ≈ 2.2T_d for a compression ratio r = 4 (Park et al. 2012). Such a prediction is in good agreement with the "bump" around a few times E_sh observed at quasi-perpendicular shocks in Figure 12, which shows the post-shock 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 M, 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 M.

Zoom In Zoom Out Reset image size

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 M, the total number of supra-thermal 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 ≲ 10E_sh, while at quasi-parallel shocks some ions contribute to form a power-law 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 supra-thermal ions from penetrating into the upstream for more than one gyroradius, which dramatically reduces the excitation of magnetic turbulence upstream, for two main reasons.

• 1.  
  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 fire-hose and the resonant streaming instability, are suppressed (see, e.g., Bykov et al. 2013 for a recent review).
• 2.  
  Second, and probably most important, if accelerated ions can propagate into the upstream only for a gyroradius ∼r₀, the time available for any magnetic perturbation to grow (i.e., the advection time ${\sim} r_0/v_{{\rm sh}}\approx \omega _c^{-1}$) is drastically reduced with respect to the case of a parallel shock, in which energetic particles can propagate for a diffusion length ∼v/v_shr₀ ≫ r₀.

Without CR-induced 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 long-wavelength 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 quasi-parallel at least for a fraction of the wave period. In other words, even if parallel shocks may often look oblique, their time-dependent configuration, whose period is likely determined by the period of the upstream waves, is intrinsically different from the steady structure of quasi-perpendicular shocks.

Initially parallel 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 ad-hoc pre-existent turbulence, with fluctuations of order δB/B ∼ 1, possibly on scales larger than the typical gyroradius of ions with v ∼ v_sh.

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 cross-field 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 M = 6, and three inclinations ϑ = 0°, 45°, 80°. The box measures (L_x, L_y, L_z) = (2000, 200, 200)[c/ω_p]³, with two cells for ion skin depth and 8 particles per cell; the time-step is $\Delta t=0.01 \omega _c^{-1}$.

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 ε_cr ≈ 12%, 3% and 1% for ϑ = 0°, 45° and 80°, respectively.

Zoom In Zoom Out Reset image size

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 self-generated turbulence is shown in Figure 14 for the three inclinations above. More precisely, the iso-value surfaces (10 levels in the interval B_z ∈ [ − 1, 1]) and the color code (in the legend) show the value of the component of the magnetic field orthogonal to B₀ (which is initialized in the xy-plane). The shock surface, found at x ≈ 600c/ω_p, is always marked by an increase of B_z, but the upstream and the downstream structures significantly depend on the shock obliquity.

Zoom In Zoom Out Reset image size

In the parallel case, ϑ = 0°, there are several regions (also in the upstream) in which the field is amplified at levels of δB/B₀ ≈ 1, 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 c/ω_p.

In the quasi-perpendicular case (ϑ = 80°), instead, B_z is always less than ∼0.1B₀ in the upstream, and less than ∼0.3B₀ in the downstream. The actual magnetic field (not shown in the figure) is mainly along y, and is consistently compressed in the shock transition, which is rather sharp, unlike in the parallel case. The ϑ = 45° case shows some intermediate features, exhibiting some magnetic structures in the upstream, even if less prominent than for the parallel configuration; the self-generated component of B_z, 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 cross-field diffusion. The aforementioned results have been obtained either with Monte Carlo simulations (Baring et al. 1995), or by tracking test-particles on top of the output of hybrid simulations seeded—by hand—with large-scale magnetic waves (Giacalone & Ellison 2000). In these cases the ion scattering is either prescribed (through the specification of a mean-free-path, in Monte Carlo simulations) or artificially enhanced, because of the presence of some pre-existing turbulence. Our simulations, instead, allow for a self-consistent description of both the shock structure and the ion kinematics, and are long/large enough to potentially observe ions undergoing DSA in the self-generated magnetic turbulence. Throughout all of the parameter space considered in this paper, for both 2D and 3D setups, we never find quasi-perpendicular shocks to be able to spontaneously generate long-wavelength waves upstream of the shock, and, in turn, to accelerate particles into the DSA regime.

We used hybrid simulations of collisionless shocks to calculate under which conditions ion acceleration is efficient at non-relativistic collisionless shocks. We can summarize our main results in the following points.

• 1.  
  We showed, for the first time within a PIC/hybrid approach, that DSA can be very efficient (10%–20% in energy) at quasi-parallel shocks, producing non-thermal ion spectra with the expected universal distribution, a power-law ∝p⁻⁴ in momentum, i.e., ∝E^−1.5 for non-relativistic ions.
• 2.  
  The extent of the non-thermal power-law tail increases with time: energetic ions are able to self-generate 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).
• 3.  
  The ion spectrum immediately behind the shock shows a peculiar shape at energies 2 ≲ E/E_sh ≲ 10, which deviates from both a Maxwellian and the DSA power-law. This supra-thermal region contains information about the post-shock thermalization, and the processes responsible for injection of particles into the acceleration mechanism.
• 4.  
  At quasi-parallel shocks, the fraction of ions that undergo acceleration is ∼10⁻³ of the total. A good parameterization of the thermal–non-thermal transition is obtained by assuming that the DSA power-law tail is attached to the Maxwellian at the injection momentum p_inj ≈ 3–4p_th, where p_th is the downstream thermal momentum.
• 5.  
  DSA occurs only at parallel and quasi-parallel shocks. Conversely, ϑ ≳ 45° shocks show evidence of SDA, but they are ineffective in accelerating particles above ∼10E_sh. 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 pre-existing turbulence may modify the acceleration properties of the shocks.
• 6.  
  Magnetic field amplification goes along with efficient ion acceleration: the self-generated turbulence is present only at quasi-parallel shocks, and the total amplification achievable in the simulation is larger for larger Mach numbers (up to δB/B₀ ≈ 6 upstream, for M = 50).
• 7.  
  The presence of accelerated ions produces an upstream precursor, in which the fluid is slowed down and compressed. The magnetic field amplification induced by CR-induced instabilities also leads to a non-adiabatic heating of the background plasma (turbulent heating), in such a way that plasma and magnetic field pressure are almost in equipartition throughout the precursor.
• 8.  
  When acceleration is efficient, the global shock dynamics is significantly modified by the presence of non-thermal particles: in this regime of CR-modified 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 non-relativistic 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 self-generated magnetic turbulence (also see CS13), the details of particle diffusion, and the mechanisms leading to ion injection.

We thank L. Gargaté for providing a version of dHybrid, P. Blasi, E. Amato, and M. Kunz for stimulating discussions, and an anonymous referee for thorough comments. This research was supported by NSF grant AST-0807381 and NASA grant NNX12AD01G, and facilitated by the Max-Planck/Princeton Center for Plasma Physics. This work was also partially supported by a grant from the Simons Foundation (grant 267233 to A.S.), and by the NSF under grant No. PHYS-1066293 and the hospitality of the Aspen Center for Physics. Simulations were performed on the computational resources supported by the PICSciE-OIT 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. DE-AC02-05CH11231, and Teragrid/XSEDE's Ranger and Stampede under allocation No. TG-AST100035.