ION DYNAMICS AND ACCELERATION IN RELATIVISTIC SHOCKS

The dynamics of relativistic particle acceleration in collisionless shocks is of great importance to several astrophysical scenarios. The acceleration of electrons, positrons, or ions in various structures such as active galactic nuclei, gamma-ray bursts, pulsar wind nebulae, and supernova remnants results in energetic particles that can then scatter in magnetic and electric fields and emit synchrotron radiation (see, for instance, Jones & Ellison 1991). Despite their relevance to understanding the radiation collected in astronomical observations, the underlying processes inherent to the acceleration are not yet fully understood. In this context, full kinetic simulations play an important role in the assessment of particular physical mechanisms relevant for astrophysical shocks, namely in the study of the nonlinear growth of the Weibel instability (Weibel 1959; Medvedev & Loeb 1999) with magnetic field generation (Silva et al. 2003; Fonseca et al. 2003; Frederiksen et al. 2004; Nishikawa et al. 2005), and particle acceleration in the precursor region (Silva 2006). Nevertheless, the self-consistent modeling of a relativistic collisionless shock, from first principles, is computationally very demanding, as large temporal and spatial scales push the numerical techniques to conditions still unexplored. Hence, full kinetic simulations require massive computational resources, optimized algorithms, methods for improved energy conservation, and also advanced visualization diagnostics.

Recently, progress was made in the understanding of electron–ion shock formation and electron–positron acceleration in relativistic unmagnetized shocks (Spitkovsky 2008a, 2008b). These studies with particle in cell (PIC) simulations confirmed the capability of these structures to effectively accelerate electrons, which is identified by the development of a nonthermal tail in the energy spectrum. Similar approaches with large-scale self-consistent modeling can provide valuable input to improve Monte Carlo methods (e.g., Bednarz & Ostrowski 1998; Ellison & Double 2004), and to support the development of analytical models (e.g., Kirk et al. 2000; Achterberg et al. 2001; Keshet 2006).

Here, we examine a relativistic electron–ion unmagnetized shock with ab initio relativistic PIC simulations, and apply a full particle tracking diagnostic to better understand the particle dynamics and the acceleration mechanism. We then show that nonthermal particle acceleration occurs through a small number of scatterings in the shock front, in a Fermi-like process (ΔE = αE) with increasing time between each energy gain. Our simulations and initial data analysis follow and confirm previous results obtained with a different PIC framework (Spitkovsky 2008a, 2008b), which reveals the robustness of the numerical results with distinct algorithms and implementations. In Section 2, we present the simulation results of the electron–ion plasma shock with a reduced ion-to-electron mass ratio of 32. In Section 3, we leverage on the OSIRIS (Fonseca et al. 2002) particle tracking and data processing tools to analyze the particle dynamics and their acceleration process, by focusing on the time evolution of the main physical quantities of the most energetic particles. A discussion of the acceleration mechanism and the conclusions are presented in Section 4.

Numerical simulations were performed with OSIRIS, a fully relativistic, electromagnetic, and massively parallel PIC code, which has been used in many different physical scenarios, such as astrophysics (e.g., Silva et al. 2003), laser/plasma accelerators (e.g., Mangles et al. 2004; Tsung et al. 2004), nanoplasma dynamics (e.g., Peano et al. 2006), and fast-ignition (e.g., Ren et al. 2004).

We simulate a two-dimensional system of a cold unmagnetized electron–ion plasma with mass ratio m_i/m_e = 32 (m_i and m_e the ion and electron mass, respectively), and evolve it until evidence of a nonthermal acceleration tail in the downstream particle spectrum is achieved. To generate the shock, the plasma stream is launched from the right wall with proper velocity u = γβ = −20, and minimal thermal dispersion from randomized particle injection. This neutral plasma stream is reflected from a rigid boundary at the left wall (this is one of the most direct methods of generating shocks in simulations; see, for instance, Forslund & Freidberg 1970; Jones & Ellison 1991, or Spitkovsky 2008a, 2008b). The computational domain is 50c/ω_p in the transverse direction and 280c/ω_p in the longitudinal direction, with c/ω_p = (4πe²n/γm_ic²)^−1/2 the ion skin depth for a number density n and relativistic ion mass γm_i; e is the elementary charge and c the speed of light in vacuum. A time step of 0.012/ω_p is used. The system is numerically resolved with 10 cells per electron skin depth in both directions, thus ensuring that the dynamics of the lighter species is accurately modeled. Taking advantage of second-order particle shapes (Esirkepov 2001) and current smoothing compensation, we use two particles per cell (ppc) for each species, which is equivalent to 6 ppc in the shocked gas, for which the energy conservation is equivalent to 16 ppc with linear particle shapes (Fonseca et al. 2008). Lower mass ratios, higher grid resolutions, and more ppc were tested, showing an overall qualitative and quantitative result convergence.

The main physical processes that generate the shock are observed when the reflected particle stream interacts with the incoming plasma, which leads to the growth of the Weibel instability, particle thermalization, and the generation of electric and magnetic turbulence. The turbulence slows down the flow, which generates the shock as a density compression that propagates in the positive x₁ direction (Figure 1(a)). Simulation results agree with the hydrodynamical jump conditions (Blandford & McKee 1976; Spitkovsky 2008a); the steady state velocity of the shock, β_shock ≃ 0.48, and the corresponding density compression obtained directly by particle number conservation, n₂/n₁ = 3.1. We emphasize that, given the simulation configuration, all quantities discussed are measured in the downstream frame.

Zoom In Zoom Out Reset image size

Figure 1 also shows other relevant physical quantities after the shock has achieved a steady state. Of relevance to the acceleration process, transverse electric fields (Figure 1(b)) arise in the linear stage of the Weibel instability, associated with space charge effects, as the two counter-propagating plasma streams pinch/filament with different rates because of their different temperatures (Tzoufras et al. 2006). The spatial symmetries of this field have direct impact on the overall transverse momentum acquired by the particles (compared to Section 3). In addition, the energy deposited in the magnetic field reaches _B ≡ B²/4πγnmc²≃ 15%–20%, similar to pair shocks (Spitkovsky 2008b), and to the m_i/m_e = 16 case. This observation, coupled with the average value of $\mathbf v_{\mathrm{drift}}= \mathbf E \times \mathbf {B}/B^2$ across the shock front associated with the structure of the self-consistent fields, also suggests the origin of the particle trapping mechanism in the shock front for both positrons/ions and electrons (S. F. Martins et al. 2009, in preparation), similar to what is observed in Earth bow shocks (Burgess 2007).

In accordance with results obtained by Spitkovsky (2008b), the energy spectrum of the ions (Figure 1(c)–(f)) is significantly different across the longitudinal direction. The upstream region is dominated by the quasi-monoenergetic negative flow of particles, and contains a lower density returning stream of heated particles that already escaped the shock region (or were never trapped). Despite the strong thermalization of the shocked gas, evidence of nonthermal particle acceleration can be observed in the downstream ion spectrum (Figure 1(d)), where a fit to a pure relativistic Maxwellian does not account for the high-energy tail. The nonthermal spectrum of both electrons and ions can be fitted with a power law (γ^−p, with p = 2.3 − 2.6) and an exponential cutoff defined by exp [ − (γ − γ_cut)/Δγ_cut], with γ_cut = 80 and Δγ_cut = 15 for the ion species. Electrons reach energies that are higher by the mass ratio m_i/m_e, thus spanning more than three orders of magnitude in energy. The high-energy tail (for γ>40) has ∼2% of the total number of ions in the downstream slice analyzed, and accounts for ∼10% of the total ion energy in that spatial region. This confirms previous results obtained for a pair plasma configuration by Spitkovsky (2008b), leading to the conclusion that the generated spectra and acceleration efficiency are not very sensitive to the mass ratio of the species, at least in two-dimensional simulations. Finally, additional propagation time leads to a linear increase in the nonthermal tail span to higher energies, and results indicate a fit with the same power-law index.

For a detailed analysis of the particle dynamics and acceleration mechanism, the OSIRIS particle tracking features (Fonseca et al. 2008) were used to follow the trajectories of the most energetic particles, selected in a first scanning simulation.

Acceleration occurs for both longitudinal and transverse momentum, but with different dynamics in each direction (Figure 2). On one hand, the approximately null space average of the transverse electric field (Figure 1(b)) leads to a symmetric acceleration of the particles in the x₂ direction (Figure 2(c)–(d)). On the other hand, the longitudinal momentum shows an average increase over time, and the particles reach farther into the upstream region, until escaping (Figure 2(a)) mostly to the downstream region. The propagation of the particles encompasses several magnetic rotations in which momentum is transferred between directions, leading to the oscillations observed in u₁ × t and u₂ × t. The ions are able to remain inside the accelerating structure by performing long drifts in the transverse direction (u₂ ≫ u₁). For the most energetic particle (colored trajectory), the final energy gain occurs with high angle variation from the transverse to the longitudinal direction (u₂ nearly constant in Figure 2(c)–(d)); the particle is then reflected from the upstream, and escapes the shock region into the downstream plasma (final u₁ < 0). After escaping, the particle sustains a large transverse momentum, and performs long drifts in the transverse direction with constant energy. Also, it should be emphasized that the behavior observed for this particle is representative of the accelerated particles.

Zoom In Zoom Out Reset image size

The energy gain and the interaction with the shock region are depicted in Figure 3(a), which shows the time evolution of the longitudinal position relative to the shock front for the 80 most energetic ions. Identical to the particle motion in pair plasmas observed in Spitkovsky 2008b, ions gain energy after being trapped as they perform multiple oscillations in the shock region until they finally escape, mainly to the downstream region. The wall reflections observed in Figure 3(a) are close to the injection point only at the beginning of the simulation, and do not affect the overall process. After the shock is formed, particles from the unshocked gas can be directly trapped, without reaching the downstream. The particle gains energy from the electric fields of the upstream, and then crosses the shock region until being reflected in the downstream (Figure 3(b)). We emphasize that, since the simulation is performed in the downstream frame, no significant acceleration occurs when particles are reflected on the downstream shocked gas. Accelerations occur rapidly and ΔE ≃ E is typically observed (Figure 3(b)), as expected in a Fermi process (Fermi 1949). The maximum energy reached is γ_final ≃ 170. Since for a single bounce ΔE ≃ E, and the initial ion energy is γ_initial ≃ 18, that final energy can only be achieved after several shock crossings. It is important to note that, unlike a Fermi mechanism, the final energy is usually obtained after only 3–5 effective collisions in a continuously evolving shock.

Zoom In Zoom Out Reset image size

In the standard formalism of Fermi acceleration, the energy evolution is written as E(N) = E₀exp(N), with E₀ the initial energy and N the number of energy gains (assumed very large: N ≫ 1), and its time dependence can therefore be estimated by relating N with time t. The initial estimate by Fermi assumed a constant time between energy gains τ_coll = τ₀, thereby leading to an exponential growth of the energy with time: E(t) = E₀exp(t/τ₀). To account for the small number of discrete energy gains observed in the simulation, we now write E(N) = E₀(1 + α)^N, where we include the constant fractional energy gain α ⩽ 1 in the form ΔE = αE, already considered by Fermi. The value α was estimated with an individual analysis of the energy gains of each tracked particle. Fitting the data of 52 of the 80 trajectories to (1 + α)^N yields α ≃ 0.81 (particular energy fluctuations of the remaining particles did not allow for a clear identification of all scatterings, defined as x₁ propagation direction inversion with greater than 50% energy gain—see arrows in Figure 3(b)).

Further analysis of several particle trajectories in the simulation shows an effect that decreases the energy growth with time, namely the increase of time between energy gains as the particle accelerates and the shock structure evolves. A simple model can be obtained using τ_coll = τ₀ + st, with τ₀ the initial collision time, and where s is the rate of change of the time between collisions. We thus get an approximate fit N(t) = t/(τ₀ + st), and the energy evolution in time becomes $E(t) \simeq E_0 (1+\alpha)^{t/(\tau _0 + st)}$. We emphasize that this expression is only valid for a limited time t, and thus implicitly assumes that particles escape the accelerating region with a maximum of N(t → ∞) = 1/s collisions. Figure 3(c) presents a fit to E(t) with s = 0.11 and τ₀ = 50.6/ω_p. Alternatively, the collision time can be written as a function of energy with τ_coll = τ₀[E(t)/E₀]^d (an energy dependence of τ_coll is also observed in other scenarios, as in the Earth's bow shock acceleration, Kis et al. 2004). For this case, a numerical fit yields τ_coll ∝ E(t)^0.24, for the same τ₀, which assumes no time domain restrictions, as opposed to τ_coll ∝ t. The parameters of the model, particularly the fractional energy gain α, and the growth rate s, are very similar to those obtained with a mass ratio of 16. Nevertheless, the parameter study required to explore these dependencies, and to completely understand the microphysics underlying the parameter values, is beyond the scope of the present paper and will be tackled in future work.

Ab initio full PIC simulations have been presented for a two-dimensional relativistic collisionless shock propagating in an initially unmagnetized electron–ion plasma (mass ratio of 32). The shock structure and jump conditions confirm previous results obtained with a different PIC framework (Spitkovsky 2008a, 2008b). Nonthermal particle acceleration is observed and occurs as particles are trapped and oscillate in the shock front, similar to a Fermi acceleration process. Nevertheless, specific distinctions exist from the standard Fermi mechanism, namely the small number of scatterings and the continuous evolution of the shock structure where the particle accelerates. When gaining energy, the particle usually performs a rotation to the transverse direction, and thus remains in the shock region, being susceptible to further acceleration.

An important consequence of the overall acceleration efficiency (∼10% of energy carried by the most energetic particles) is the increased relevance of nonlinear effects for the shock structure and for the acceleration process. In fact, when the accelerated particles yield ≳10% of the plasma energy, their dynamics become relevant and influences the evolution of the overall system (Jones & Ellison 1991). One of these nonlinear effects is the dynamic pressure of the accelerated particles that slows down the unshocked plasma before it reaches the sharp shock transition. This is indicated in the simulation results as the shock widens and the magnetic and electric field layers extend to the upstream (Keshet et al. 2008; Medvedev & Zakutnyaya 2008). Furthermore, the inclusion of the trapped particles population can have implications on the jump conditions, and current models can be extended to incorporate these effects into the equation of state in the shock region. This generalization has been made for the nonrelativistic and electrostatic shocks case (Forslund & Freidberg 1970; Sorasio et al. 2006).

The average energy growth of the particles can be reproduced by a multiscattering acceleration mechanism with ΔE = αE, assuming an increase in the time between collisions. The long-term acceleration implications of these effects cannot be inferred from our simulation spectra because of the exponential cutoff due to the finite simulated size that limits the injection of particles. Larger-scale simulations with longer propagation distances and larger acceleration times will elucidate the properties of the particle spectrum at higher energies, thus allowing for a detailed identification of the mechanisms responsible for α ≲ 1, and for an increase in time between collisions.

In summary, our results confirm the possibility of particle acceleration through a Fermi-like mechanism with a reduced number of energy gains, and generalized to reproduce the statistical data obtained with ab initio full PIC simulations that self-consistently resolve the turbulent and nonlinear evolution of the shock.

This work was partially supported by Fundação Calouste Gulbenkian and by Fundação para a Ciência e Tecnologia under grants SFRH/BD/35749/2007 and PTDC/FIS/66823/2006 (Portugal). The simulations presented were produced using the IST Cluster (IST/Portugal) and the Dawson cluster (UCLA). The authors would like to thank Profs. A. Spitkovsky, S. Schwartz, and M. Medvedev for useful discussions.