ACCELERATION OF PARTICLES AT THE TERMINATION SHOCK OF A RELATIVISTIC STRIPED WIND

The broadband spectrum of Pulsar Wind Nebulae (PWNe), extending from the radio up to the gamma-ray band, is usually modeled as synchrotron and inverse Compton radiation from a nonthermal population of electron–positron pairs. It is assumed that the emitting particles are accelerated to a power-law distribution at the so-called termination shock, where the momentum flux of the relativistic pulsar wind is balanced by the confining pressure of the nebula. However, it is still unclear how the wind termination shock can accelerate electrons and positrons (hereafter, simply "electrons") to the relativistic energies required to power the nebular emission.

A flat radio spectrum is the common observational feature of PWNe, whose spectral flux in the radio band scales as $F_{\nu }\propto \nu ^{-\alpha _r}$, with α_r = 0.0–0.3. For the Crab Nebula, the best studied example of a PWN, the synchrotron cooling time of radio-emitting electrons largely exceeds the age of the nebula, so one cannot exclude that such particles may be "primordial," i.e., accelerated at a very early stage of the PWN evolution (Atoyan 1999). However, observations of wisps in the radio band (Bietenholz et al. 2001) seem to suggest that radio-emitting electrons are being accelerated now, in the same region as those responsible for the optical and X-ray emission. In this case, the observed radio spectral index (α_r ∼ 0.25 for the Crab; Bietenholz et al. 1997) would imply a distribution of emitting electrons of the form dN/dE∝E^−p, with slope p = 2α_r + 1 ∼ 1.5. To explain the radio-to-optical emission of the Crab, the power law of shock-accelerated electrons should span at least three decades in particle energy (e.g., Lyubarsky 2003). Flat electron spectra with slope p < 2 below GeV energies are also required to model the radio emission of hotspots in radio galaxies (Stawarz et al. 2007) and X-ray observations of luminous blazar sources (Sikora et al. 2009).

For a power-law particle spectrum with 1 < p < 2, most of the particles lie at the low-energy end of the distribution, but the plasma energy content is dominated by particles at the high-energy cutoff. This contrasts with the energy spectrum expected from Fermi acceleration in relativistic shocks, which normally yields slopes p > 2 (e.g., Achterberg et al. 2001; Keshet & Waxman 2005), in which case low-energy particles dominate both by number and by energy. To produce broad spectra with p < 2, a mechanism is required that raises the mean electron energy, but leaves most of the particles at relatively low energies. In the context of pulsar winds, Hoshino et al. (1992) and Amato & Arons (2006) have proposed that, in a wind loaded with ions, resonant absorption of ion cyclotron waves by electrons and positrons could generate flat distributions downstream from the termination shock. However, the required ion injection rate largely exceeds what is expected from standard models of pulsar magnetospheres (e.g., Arons 2007).

An alternative possibility was discussed by Lyubarsky (2003), under the assumption that the flow upstream of the termination shock consists of alternating stripes of opposite magnetic polarity, separated by current sheets of hot plasma (from now on, a "striped wind"). For obliquely rotating pulsars, this is the configuration expected around the equatorial plane of the wind, where the sign of the toroidal field alternates with the pulsar period. If the wind remains dominated by Poynting flux until the termination shock (which is still an open question; see Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003), shock-compression of the stripes may drive annihilation of the alternating fields, a process known as "driven magnetic reconnection." The resulting transfer of energy from the fields to the particles could potentially generate flat electron distributions, with slope 1 < p < 2.

The physics of magnetic reconnection can be captured self-consistently only by means of multi-dimensional particle-in-cell (PIC) simulations. Fully kinetic PIC simulations provide a powerful tool to explore the microphysics of collisionless plasmas, since they can capture from first principles the fundamental interplay between charged particles and electromagnetic fields. In the context of relativistic magnetic reconnection in pair plasmas, most studies have explored the so-called undriven reconnection process (Zenitani & Hoshino 2001, 2005a, 2007, 2008; Jaroschek et al. 2004), where field annihilation is initiated by a transient seed perturbation to an otherwise stable current sheet. As discussed above, this is not the setup expected at the wind termination shock. Here, it is the shock-compression of the flow that steadily drives regions of opposite magnetic field polarity toward each other, causing reconnection.

Shock-driven reconnection in magnetically dominated flows has been studied by Pétri & Lyubarsky (2007) with analytical methods and one-dimensional (1D) PIC simulations.¹ They find a condition on the wind magnetization and the stripe wavelength needed for full dissipation of the alternating fields at the termination shock (roughly speaking, the post-shock Larmor radius should be larger than the separation between current sheets). However, in 1D all particles gain energy at the same rate, which results in a Maxwellian distribution downstream from the shock. The generation of nonthermal tails is possible only if different particles experience different energy gains, which requires multi-dimensional simulations. As a preliminary step toward multi-dimensional studies of the termination shock in a striped wind, Lyubarsky & Liverts (2008) have performed 2D simulations of driven magnetic reconnection. Two stripes of opposite magnetic polarity were compressed by means of an external force, which would imitate the effect of a shock. They found that driven magnetic reconnection can produce flat nonthermal tails, with p ∼ 1. Yet, their study was limited to a single current sheet, and most importantly it did not self-consistently account for the presence of a shock.

In this work, we explore via multi-dimensional PIC simulations the acceleration of particles at the termination shock of a striped relativistic electron–positron wind. The physics of a shock interacting with multiple current sheets is captured self-consistently in 2D and 3D, and the transfer of energy from the field to the particles via shock-driven reconnection is investigated from first principles. We find that the alternating fields are completely dissipated upon compression by the shock, and their energy is transferred to the particles, regardless of the properties of the flow. This contrasts with the conclusions of Pétri & Lyubarsky (2007), who found efficient field dissipation only for high magnetizations and short stripe wavelengths. We show that the reduced dimensionality of their 1D model could not correctly capture the growth of the so-called tearing-mode instability, which is of paramount importance for field annihilation in 2D and 3D configurations. We find that broad particle spectra with flat slopes (1 < p < 2) are a common by-product of shock-driven reconnection, but the extent of the power-law tail depends on the wind magnetization and the stripe wavelength.

This work is organized as follows. In Section 2, we discuss the setup of our simulations and the magnetic field geometry. In Section 3, the properties of the termination shock are investigated for one representative choice of wind parameters. We discuss both the shock structure, with focus on the physics of shock-driven magnetic reconnection, and the process of particle acceleration to nonthermal energies. In Section 4, we study how the particle spectrum downstream from the shock depends on the physical conditions in the wind, like magnetization and stripe wavelength. Although we mostly employ 2D simulations, in Section 5 we compare 2D and 3D results, finding good agreement. We summarize our findings in Section 6 and comment on the application of our results to astrophysical scenarios.

We use the 3D electromagnetic PIC code TRISTAN-MP (Buneman 1993; Spitkovsky 2005) to study the termination shock of a relativistic striped wind. The shock is set up by reflecting a magnetized electron–positron flow from a conducting wall located at x = 0 (Figure 1). The interaction between the incoming beam (that propagates along $-\mbox{\boldmath {$\hat{x}$}}$) and the reflected beam triggers the formation of a shock, which moves away from the wall along $+\mbox{\boldmath {$\hat{x}$}}$ (Figure 1). The simulation is performed in the "wall" frame, where the downstream plasma is at rest.

Zoom In Zoom Out Reset image size

To follow the shock evolution for longer times with fixed computational resources, we mainly utilize 2D computational domains in the x–y plane, with periodic boundary conditions in the y direction. In Section 5, we compare 2D and 3D simulations and show that 2D runs can capture most of the relevant physics. For both 2D and 3D domains, all three components of particle velocities and electromagnetic fields are tracked.

The incoming electron–positron stream is injected along $-\mbox{\boldmath {$\hat{x}$}}$ with bulk Lorentz factor γ₀. The flow carries a strong magnetic field, oriented along the y direction. It follows that the shock will be "perpendicular," i.e., with magnetic field orthogonal to the shock normal. The choice to initialize the field in the simulation plane (i.e., oriented along y, as opposed to z) is motivated by the agreement between 2D simulations with in-plane fields and 3D experiments, as we show in Section 5. The spatial profile of the injected magnetic field, sketched in Figure 1, is

$$\begin{eqnarray} \!\!\!B_{y}(x,t)&=&B_0 \tanh \left\lbrace \frac{1}{\delta }\left[\alpha +\cos \left(\frac{2\pi (x+\beta _0 c t)}{\lambda }\right)\right]\right\rbrace\nonumber \\ &\equiv & B_0 \tanh \zeta \end{eqnarray} \tag{ 1 }$$

and its associated motional electric field is E_z = β₀B_y. Here, β₀ is the three-velocity of the injected plasma. This particular choice for the magnetic field is motivated by analytical studies of the structure of pulsar winds (e.g., Pétri & Kirk 2005). At each given time, Equation (1) shows that the field profile is periodic with wavelength λ (see Figure 1). Current sheets, where the magnetic field vanishes, are identified by the condition ζ ≃ 0, and their half-thickness is Δ = λδ/2π ≪ λ. Across each current sheet, the field reverses its sign from +B₀ to −B₀ (or vice versa), where B₀ is the field strength in the region outside the current sheets, which we shall call "cold wind." In the following, we parameterize the field strength via the so-called magnetization parameter σ ≡ B²₀/4πγ₀mn_c0c², taken to be much larger than unity. Here, m is the electron (or positron) mass and n_c0 is the density of particles in the cold wind. By defining the relativistic skin depth in the cold wind $c/\omega _{p}\equiv \sqrt{\gamma _0 m c^2/ 4\pi e^2 n_{c0}}$, and the relativistic Larmor radius r_L ≡ γ₀mc²/eB₀, we can rewrite the magnetization parameter as σ = (c/ω_p)²/r²_L. Finally, the parameter α (−1 < α < 1) is a measure of the net field (i.e., averaged over one wavelength λ), such that 〈B_y〉_λ/B₀ = α/(2 − |α|). Although the magnetic field intensity in the cold wind is always B₀, the wavelength-averaged field is not necessarily zero, because two neighboring stripes generally have different widths (see Figure 1). In the pulsar wind, one expects 〈B_y〉_λ = 0 (which corresponds to α = 0) only in the equatorial plane (see the upper panel in Figure 1). We remark that a net field 〈B_y〉_λ ≠ 0 does not play the role of a guide field, which in our geometry would correspond to a uniform component of magnetic field in the z direction. No significant guide field is expected at the termination shock of pulsar winds, and the physics of shock-driven reconnection in the presence of a guide field will be discussed elsewhere.

In the cold wind, each computational cell is initialized with two electrons and two positrons, with a small thermal spread Θ_c ≡ kT_c/mc² = 10⁻⁴. We have also performed limited experiments with a larger number of particles per cell (up to 32 per species), obtaining essentially the same results. Such a background of cold particles is initialized also within the current sheets, but with the addition of a hot component whose density profile is n_h = n_h0/cosh ²ζ. Here, n_h0 ≡ η n_c0 is the density of hot particles at the center of the current sheet. The temperature of this hot component is determined by pressure balance between the gaseous part (inside the current sheet) and the magnetic part (in the cold wind), which yields a thermal spread Θ_h ≡ kT_h/mc² = σ/2η for a 3D relativistic Maxwellian. Finally, the variation of magnetic field across a current sheet needs to be sustained by a current flowing along $\mbox{\boldmath {$\hat{z}$}}$. In the wind frame, electrons and positrons in the current sheet should be drifting in opposite directions with the same speed $\beta _{h}=\sqrt{\sigma }/(\eta \gamma _0)\,(c/\omega _{p})/\Delta$. In summary, to initialize the distribution of hot electrons (and positrons) in the frame of the simulations, we need to make two successive Lorentz boosts: first, from their proper frame to the wind frame, where electrons and positrons counter-stream along $\mbox{\boldmath {$\hat{z}$}}$ with velocity ±β_h; then, from the wind frame to the simulation frame, where the wind propagates along $-\mbox{\boldmath {$\hat{x}$}}$ with Lorentz factor γ₀.²

In our study, we adopt γ₀ = 15 as our fiducial bulk Lorentz factor. As discussed in Section 4, we actually explore a wide range of Lorentz factors, from γ₀ = 3 to γ₀ = 375, and find basically the same results, modulo an overall shift in the energy scale. We vary the stripe wavelength by almost two orders of magnitude, from λ = 20 c/ω_p to λ = 1280 c/ω_p, and we discuss its effects on the structure of the shock and the process of shock-driven reconnection in Section 4.1. The dependence on the wind magnetization, though always in the regime of highly magnetized flows (σ ≫ 1), is presented in Section 4.2, where we vary σ from 10 to 100. Different values of the stripe-averaged magnetic field are investigated in Section 4.3, where we vary the parameter α = 2〈B_y〉_λ/(B₀ + |〈B_y〉_λ|) between 0.0 and 0.95. For pulsar winds, large values of |α| would be expected at high latitudes above the midplane, whereas α = 0 along the equator (see the upper panel in Figure 1).

The structure of the current sheet is completely determined once we fix the value of the density excess η = n_h0/n_c0, and of the sheet half-thickness Δ = λδ/2π. We choose η = 3, motivated by analytical studies of pulsar winds (Lyubarsky & Kirk 2001). However, we find that our results do not appreciably depend on the value of η (we have also tried with η = 1.5 and η = 6), provided that the contribution of current sheets to the upstream particle flux stays negligible (i.e., ηδ ≪ 1). Regarding the sheet half-thickness, we typically employ Δ = c/ω_p (for each choice of λ, this determines the value of δ = 2πΔ/λ), but we have verified that our results are not very sensitive to this parameter (we have tried up to Δ = 8 c/ω_p, for λ = 320 c/ω_p). The choice for the sheet thickness is guided, on the one hand, by the requirement that the contribution of current sheets to the upstream plasma energy flux should be small (i.e., σδ ≪ 1). On the other hand, we enforce Δ ≳ 0.5 c/ω_p, to ensure that the incoming current sheets are stable to the process of undriven reconnection. With this choice, we are confident that the trigger for magnetic dissipation will come entirely from the shock-driven compression of the flow.

We choose the relativistic electron (or positron) skin depth in the cold wind c/ω_p such that the smallest scale in the system, which for σ > 1 is the relativistic Larmor radius $r_{L}=(c/\omega _{p})/\sqrt{\sigma }$, is resolved with at least a few computational cells.³ We resolve c/ω_p with 7.5 computational cells for σ = 10, 10 cells for σ = 20, 13 cells for σ = 30, 17 cells for σ = 50, and 24 cells for σ = 100. For each value of the magnetization, we have verified that our results do not substantially change when doubling the resolution (for σ = 10, we have tried with up to c/ω_p = 30 cells). Our computational domain is typically ∼400 c/ω_p wide (along the y direction). In particular, we choose a transverse size of 3072 cells for σ = 10, of 4096 cells for σ = 20, of 5120 cells for σ = 30, of 7168 cells for σ = 50, and of 9216 cells for σ = 100. A large box is of paramount importance for the consistency of our findings. As we show in Appendix A, the results of our simulations converge only when the computational box is larger than the half-wavelength of the stripes. It follows that 1D simulations, as the ones presented by Pétri & Lyubarsky (2007), cannot capture all the relevant physics, especially for long wavelengths.

The results presented in the following sections refer to shocks that have already reached a steady state. For the typical parameters we explore, this happens before ω_pt = 3000, but in some cases (especially for long stripe wavelengths), we need to evolve our simulations up to longer times. In Appendix B, we follow the entire time evolution of a representative shock, from the earliest stages until it approaches a steady state.

In this section, we present the properties of a shock that propagates in a relativistic electron–positron flow with alternating fields (a "striped wind"). The case of uniform pre-shock fields has been investigated by Sironi & Spitkovsky (2009) and Sironi & Spitkovsky (2011) via multi-dimensional PIC simulations (for electron–positron and electron–ion plasmas, respectively). They studied how the magnetic field strength and obliquity can affect the efficiency of the Fermi process, where particles gain energy by bouncing back and forth across the shock (e.g., Blandford & Eichler 1987). They found that in perpendicular shocks, i.e., with field orthogonal to the shock normal, the Fermi mechanism can operate only if the flow is weakly magnetized. For magnetizations σ ≳ 10⁻³, the field is so strong that charged particles are forced to slide along the field lines, whose orientation prohibits repeated crossings of the shock, thus inhibiting the Fermi process.

In a medium with alternating fields, we find that the efficiency of Fermi acceleration is controlled by the magnitude of the stripe-averaged field 〈B_y〉_λ. Striped flows with 〈B_y〉_λ/B₀ ≪ 1 are equivalent to weakly magnetized shocks in uniform fields, with respect to particle acceleration via the Fermi process. In this section, we choose a value for the stripe-averaged field 〈B_y〉_λ/B₀ ≃ 0.05 (corresponding to α = 0.1) that is large enough to inhibit the Fermi process. This allows us to focus solely on the process of shock-driven reconnection and to isolate its effects on the shock structure and the particle spectrum. In pulsar winds, the special case 〈B_y〉_λ = 0 should be realized only along the equatorial plane of the wind. By choosing α = 0.1, we address the generic properties of the flow above (or below) the midplane.

In this section, we adopt γ₀ = 15 and σ = 10 as our fiducial values for the upstream bulk Lorentz factor and magnetization. To clarify the physics of shock-driven reconnection, we employ relatively large values for the stripe wavelength (λ = 320 c/ω_p or λ = 640 c/ω_p), and we choose α = 0.1, corresponding to a stripe-averaged field 〈B_y〉_λ/B₀ ≃ 0.05. The dependence of the shock properties on γ₀, σ, λ, and α is investigated in Section 4. There, we also explore the special case α = 0, where shock-driven reconnection coexists with Fermi acceleration.

3.1. Shock Structure

The typical structure of a relativistic shock propagating in a striped wind is shown in Figures 2 and 3, which refer to ω_pt = 3750, after the flow structure has reached a steady state. The longitudinal profile of the transition region, along the direction x of shock propagation, is presented in Figure 2, whereas Figure 3 shows the flow properties in the x–y plane of the simulation, zooming in on a region around the shock (as delimited by the vertical dashed red lines in Figure 2(a)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The incoming flow (x ≳ 2600 c/ω_p) carries an alternating magnetic field of wavelength λ = 640 c/ω_p (Figure 2(c); black line for B_y, red for E_z). The field reverses its polarity across current sheets filled with hot overdense plasma (see the spikes in density and average particle energy of Figures 2(a) and (d), respectively), where the magnetic energy is much lower than in the cold wind (see the troughs in the magnetic energy fraction _B = B²/8πγ₀mn_c0c², black line of Figure 2(b)). For x ≳ 2600 c/ω_p, the only perturbation to the incoming flow comes from the so-called electromagnetic precursor, a coherent train of transverse electromagnetic waves propagating into the upstream at the speed of light (see the short-scale ripples at x ≃ 3850 c/ω_p in Figures 2(b) and (c)). The electromagnetic precursor wave is generated in the initial stages of shock evolution, before shock-driven reconnection appreciably changes the structure of the shock. At such early times, the flow resembles a magnetized shock propagating in uniform fields (see Appendix B), which is well known to radiate electromagnetic waves into the upstream (Gallant et al. 1992; Hoshino 2008; Sironi & Spitkovsky 2009, 2011).

The steady-state structure of the flow in Figure 2 shows the presence of two shocks. The main shock (which we shall call the "hydrodynamic shock" from now on, for reasons that will become clear below) corresponds to the jump in density occurring at x ≃ 1000 c/ω_p. At x ≃ 2600 c/ω_p, well ahead of the hydrodynamic shock, the incoming flow crosses a fast magnetohydrodynamic (MHD) shock, whose location is highlighted in Figure 2 by a vertical dotted blue line. At the fast shock, the wind is decelerated and compressed, as shown by the modest increase in density (Figure 2(a)), magnetic energy (black line in Figure 2(b)), and transverse magnetic field B_y (black line in Figure 2(c)). The deceleration of the flow appears most clearly in Figure 2(c), as a mismatch between the motional electric field E_z and the magnetic field B_y behind the fast shock (red and black line, respectively). Irrespective of the initial bulk Lorentz factor of the flow, the plasma drift velocity ($=-E_z/B_y\, \mbox{\boldmath {$\hat{x}$}}$) downstream from the fast shock is only moderately relativistic, since E_z/B_y ≃ 0.8 there (Figure 2(c)). The flow deceleration at the fast shock is accompanied by a partial isotropization of the particle distribution, mostly in the x–z plane orthogonal to the field. Overall, the flow kinetic energy decreases across the fast shock from γ₀ = 15 to 〈γ〉 ≃ 10 (Figure 2(d)), to compensate for the growth in magnetic energy (Figure 2(b)).

The fast MHD shock has a dramatic impact on the structure of the striped wind. Across the fast shock, the stripe wavelength shrinks from λ = 640 c/ω_p to λ ≃ 550 c/ω_p (see Figure 2(c)). Most importantly, the compression induced by the passage of the fast shock through a given current sheet triggers the onset of magnetic reconnection. Small-scale islands develop inside the current sheet as a result of the tearing-mode instability. The growth and evolution of reconnection islands is clearly shown in the 2D plots of density and magnetic energy of Figures 3(a) and (b), respectively. The current sheet at x ≃ 2500 c/ω_p, just downstream from the fast shock (located at x ≃ 2600 c/ω_p), appears partitioned into a number of overdense islands with strong magnetic fields (red point-like regions in Figure 3(b), at x ≃ 2500 c/ω_p). A magnetic X-point exists in between each pair of neighboring islands, where field lines of opposite polarity break and reconnect. As the flow recedes from the fast shock, the small-scale islands created inside each current sheet coalesce to form bigger islands. This temporal evolution corresponds, in the snapshot taken at ω_pt = 3750 in Figure 3, to a pattern of bigger (and so, fewer) islands for current sheets that are farther behind the fast shock.

As reconnection proceeds with distance downstream from the fast shock, more and more magnetic energy is released and transferred to particles. In fact, the spikes in average particle energy of Figure 2(d) become stronger for current sheets farther from the MHD shock (see also the longitudinal phase space of positrons in Figure 2(e), for 1000 c/ω_p ≲ x ≲ 2600 c/ω_p). The transfer of magnetic energy to the particles continues up to the point where reconnection islands fill the entire region between neighboring current sheets. Now, the striped structure of the flow is erased, and a hydrodynamic shock forms (located at x ≃ 1000 c/ω_p in Figures 2 and 3).

Behind the hydrodynamic shock, the flow comes to rest (in the wall frame of our simulations), and the particle distribution is isotropic in three dimensions (compare longitudinal and transverse phase spaces in Figures 2(e) and (f), respectively). The post-shock number density approaches n_d ≃ 4 n_u (Figure 2(a), where n_u is the density ahead of the fast shock), and the shock velocity is β_sh ≃ 1/3 (see Appendix B), in agreement with the jump conditions of a relativistic unmagnetized 3D plasma. In fact, most of the energy per particle downstream from the hydrodynamic shock is in kinetic form (as opposed to the dominant electromagnetic component of the incoming striped flow). This becomes clear when comparing the profile of the mean kinetic energy per particle (Figure 2(d)) with the 1D plot of the mean magnetic energy per particle (red line in Figure 2(b)).⁴ Behind the hydrodynamic shock, the average kinetic energy per particle increases from 〈γ〉/γ₀ ≃ 1 up to 〈γ〉/γ₀ ≃ σ + 1 ≃ 11, as expected in the case of full dissipation of magnetic fields. Correspondingly, the mean magnetic energy per particle drops to zero (red line in Figure 2(b)). As best seen in experiments with a smaller stripe wavelength, in the case α ≠ 0 a net magnetic field survives downstream from the shock, as predicted by Pétri & Lyubarsky (2007), and shown by the 1D profile of B_y (black line in Figure 2(c)). This residual magnetic field results from shock-compression of the stripe-averaged upstream field (recall that α = 0.1 here, corresponding to 〈B_y〉_λ/B₀ ≃ 0.05).

We point out that the picture described above accurately represents the long-term behavior of the shock. In particular, the distance between the fast and the hydrodynamic shock (≃ 1600 c/ω_p in Figures 2 and 3) stays approximately constant in time (see Appendix B). The hydrodynamic shock moves with β_sh ≃ 1/3, and the fast MHD shock remains far enough ahead of the main shock such that to guarantee that reconnection islands will fill the space between neighboring current sheets by the time the flow enters the hydrodynamic shock. For a fixed reconnection rate, this implies that the distance between the fast and the hydrodynamic shock will scale linearly with the stripe wavelength (as will the size of reconnection islands just ahead of the main shock), a prediction that we have directly verified in our simulations.

Finally, we refer to Appendix A for a detailed discussion of the differences between our 2D results and the 1D simulations presented by Pétri & Lyubarsky (2007). As apparent in Figure 3, the process of island coalescence, which is essential for the formation of the hydrodynamic shock, can be captured correctly only with multi-dimensional simulations. Indeed, for the parameters employed in this section, the 1D model by Pétri & Lyubarsky (2007) would predict negligible field dissipation, in sharp contrast with our findings. This clearly emphasizes the importance of multi-dimensional physics for our understanding of shock-driven reconnection.

3.2. Particle Spectrum and Acceleration

The particle energy spectrum at different locations through the flow, as marked by arrows at the bottom of Figure 3(d), is shown in Figures 3(e) and (f). Figure 3(e) shows the particle distribution in the cold wind, whereas Figure 3(f) focuses on the hot plasma within current sheets. Downstream from the hydrodynamic shock, the striped structure is completely erased, and no difference persists between current sheets and cold wind. Here, the particle spectrum (black line in both panels) is in the form of a flat power-law tail (with spectral index p ≃ 1.4, dashed line) extending from γ_min ≃ 30 to γ_max ≃ 500, where it cuts off exponentially. Since 1 < p < 2, most of the particles are found at γ ≃ γ_min, but most of the energy is contributed by particles at γ ≃ γ_max. For comparison, a 3D Maxwellian with the same average Lorentz factor 〈γ〉 = γ₀(σ + 1) ≃ 165 is plotted in Figure 3(e) as a dotted line, to show that the actual particle spectrum is much broader than a thermal distribution.⁵ By following the flow from the fast to the hydrodynamic shock, we now clarify how such a particle spectrum is generated.

Most of the incoming plasma is contained in the cold striped wind exterior to current sheets, and its spectrum is initially a cold Maxwellian drifting with bulk Lorentz factor γ₀ = 15 (light blue line in Figure 3(e)). After crossing the fast shock, the flow decelerates and it partly isotropizes, which explains why the particle spectrum is now broader and it peaks at lower energies (dark blue line in Figure 3(e), but see also the 2D plot of average particle energy in Figure 3(d), and its 1D projection in Figure 2(d)).

The evolution of the particle population within current sheets is much more dramatic, as shown in Figure 3(f).⁶ Ahead of the fast shock, the spectrum at x ≃ 2850 c/ω_p (yellow line in Figure 3(f)) shows, in addition to the cold plasma in the striped wind, the hot particles in the current sheet. Their distribution is initialized as a drifting 3D Maxwellian with thermal spread Θ_h = σ/2η, so that in the simulation frame the spectrum should peak at 〈γ〉 = 3γ₀Θ_h ≃ 75, as observed (we remind that σ = 10 and η = 3). Behind the fast shock (red line for x ≃ 2200 c/ω_p, brown line for x ≃ 1350 c/ω_p), the particle spectrum consists of two components: a low-energy peak, reminiscent of the particle distribution far from current sheets (compare with the dark blue line in Figure 3(e)); and a high-energy component, which grows in number and extends to increasingly higher energies, as the flow propagates toward the hydrodynamic shock. This component is populated by particles that were initially outside the current sheet, and have entered the sheet in the course of the reconnection process. They mostly end up in the reconnection islands of Figure 3, which act as reservoirs of particles (Figure 3(a)) and particle energy (the high-energy regions of Figure 3(d), colored in red, match very well the overdense islands of Figure 3(a)). In fact, the main contributions to the high-energy end of the spectrum at x ≃ 1350 c/ω_p (brown line in Figure 3(f)) come from the two major islands of that current sheet.

As magnetic reconnection proceeds with distance behind the fast shock, particles are accelerated to higher and higher energies, as we discuss below in more detail. This continues until reconnection islands grow as big as the distance between neighboring sheets, at which point the current sheets occupy the entire flow (and the hydrodynamic shock forms). This explains the smooth transition from the spectrum within current sheets to the distribution downstream of the hydrodynamic shock (yellow to black lines in Figure 3(f)). Also, since most of the particles that end up with high energies were initially part of the cold wind, the resulting spectrum behind the hydrodynamic shock is almost independent of the initial setup of the current sheet (most notably, sheet thickness Δ and overdensity η).

We now investigate in more detail the physics of particle acceleration within current sheets. In Figure 4, we follow the trajectories of a representative sample of positrons, extracted from the simulation of a striped wind with λ = 320 c/ω_p (but our conclusions hold for all wavelengths). All the selected positrons are initially in the cold wind (i.e., far from current sheets), at roughly the same x-location (within a range of 5 c/ω_p). At first, they propagate toward the fast shock, which they cross at ω_pt ≃ 700, as shown by the spreading of their energies in Figure 4(a). They encounter the closest current sheet (in this case, the one that was to the right of their initial position) at ω_pt = 1008, as shown in Figure 4(b) together with the corresponding 2D plot of particle number density. Starting from this time, their histories can diverge significantly, depending on their y-location at the moment of interaction with the current sheet. The particles that will end up with relatively low energies (γ < 3γ₀, black lines in Figure 4(a)) are found at ω_pt = 1008 exclusively around magnetic islands (their locations are indicated by filled circles in Figure 4(b)). At later times, they will get trapped on the outskirts of the growing islands, without appreciable changes in energy.

Zoom In Zoom Out Reset image size

In contrast, the positrons that will eventually reach high energies (γ > 6γ₀, red and blue lines in Figure 4(a)) are concentrated at ω_pt = 1008 in the vicinity of magnetic X-points (the particle locations are shown by open circles in Figure 4(b)). It is at X-points that magnetic field lines of opposite polarity are converging and tearing (white contours in Figure 4(b) show the magnetic field lines), and magnetic energy is transferred to particles. Starting from ω_pt = 1008 (vertical dashed line in Figure 4(a)), the energies of these positrons grow simultaneously and explosively (red and blue lines in Figure 4(a)), due to the reconnection electric field parallel to the X-line (Zenitani & Hoshino 2001; Jaroschek et al. 2004; Lyubarsky & Liverts 2008). We find that X-point acceleration represents the first stage of energization for most of the particles that will end up with high energies. Particles pre-accelerated at X-points continue to gain energy by compression (Lyubarsky & Liverts 2008; Giannios 2010), as they are advected from the X-point into the closest island. No major systematic changes of energy are observed after the particles get trapped within magnetic islands. Different X-points may provide different energy gains, as shown by the red (respectively, blue) curves in Figure 4(a), which refer to high-energy positrons that get accelerated at the X-point above (respectively, below) the small magnetic island in the middle of the current sheet. In summary, the energy gain of a given particle is determined on the one hand by the strength of the reconnection electric field at the X-point, and, on the other hand, by the time available for acceleration, as the particle flows from the X-point into the closest island.

In retrospect, the importance of X-points for the generation of nonthermal particles is not surprising. In magnetically dominated flows, particles can untie from magnetic field lines and be injected into the acceleration process only in regions where the electric field is comparable or exceeds the magnetic field, i.e., green or blue regions in the current sheets of Figure 3(c) (where we plot the difference _B − _E between magnetic and electric energy). Based on this criterion, X-points are natural sites of particle energization, since the magnetic field vanishes there. As discussed above, the main agent of particle acceleration at X-points is the reconnection electric field, which for the magnetic geometry employed in this work is oriented along z (more precisely, along $+\mbox{\boldmath {$\hat{z}$}}$ if the magnetic field B_y is negative to the left of the current sheet and positive to the right; along $-\mbox{\boldmath {$\hat{z}$}}$ in the opposite case). The sign of the reconnection electric field leaves an imprint on the particle z-momentum. As shown by the positron transverse phase space in Figure 2(f), current sheets where the reconnection field is expected to be along $+\mbox{\boldmath {$\hat{z}$}}$ are mostly populated by high-energy positrons with γβ_z > 0 (and electrons with γβ_z < 0), and the opposite holds where the reconnection field is along $-\mbox{\boldmath {$\hat{z}$}}$.

Accelerated by the reconnection electric field, high-energy particles stream at first along the z direction, and then, by effect of the reconnected magnetic field B_x, they are advected along y into the closest island. In our 2D simulations, where the z dimension is not resolved, the relevant timescale for acceleration is the time to drift from an X-point into the closest magnetic island. As the flow propagates downstream from the fast shock, magnetic islands grow and merge (see Figure 3(a)), so that the typical distance between an X-point and the closest island increases, allowing particles to be accelerated to higher energies. This explains why the upper cutoff in the particle spectra of Figure 3(f) shifts to higher energies for current sheets farther from the fast shock.⁷

Since our 2D simulations cannot resolve potential perturbations to the particle motion along z, the validity of our conclusions may be questioned, if the structure of X-points in a realistic 3D scenario is much different from the 2D picture analyzed here. In particular, in 3D we expect the current sheet to be folded, by effect of the so-called drift-kink instability (Zenitani & Hoshino 2005a, 2007). The characteristic wavelength of the kink mode in the x–z plane will introduce a different length scale that could potentially compete with the distance between X-points and islands (in the x–y plane) in constraining the maximum energy of accelerated particles. In Section 5, we address this important point and show that the latter constraint is more restrictive than the former, at least for the parameters explored in this work. It follows that our 2D simulations with in-plane fields can capture very well the main 3D properties of the shock, and in particular the physics of particle acceleration.

Although we find that the reconnection electric field at X-points is the main agent of particle acceleration, we have verified, by tracing the orbits of a large number of particles extracted directly from the simulation, that a variety of other mechanisms can play a role in particle energization. During the coalescence of magnetic islands, particles may be accelerated by the anti-reconnection electric field at the X-point located in between the two merging islands (Jaroschek et al. 2004; Oka et al. 2010). Alternatively, particle energization may be due to a kind of first-order Fermi process within magnetic islands, when particles reflect from the ends of contracting islands following coalescence (Drake et al. 2006). Also, the second-order Fermi process can operate between randomly moving islands just downstream from the hydrodynamic shock. The two latter mechanisms primarily change the components of particle momentum along the plane where the magnetic field lies (the x–y plane, in our simulations). In contrast, for reconnection and anti-reconnection at X-points, it is the momentum component along the reconnection electric field (along z, in our case) that receives the greatest boost.

The energization mechanisms discussed so far are peculiar to the process of magnetic reconnection, being seeded by the field dissipation itself, or by the resulting fluid motions. In shock-dominated systems, two more acceleration processes can be relevant, whose free energy comes from the converging motion of the pre-shock and post-shock flows. In first-order Fermi acceleration (or diffusive shock acceleration, DSA), particles diffuse back and forth across the shock front and gain energy by scattering from magnetic turbulence embedded in the converging flows (e.g., Blandford & Ostriker 1978; Bell 1978; Blandford & Eichler 1987). In magnetized flows, particles may also gain energy directly from the large-scale motional electric field while they gyrate around the shock front, a process known as shock-drift acceleration (SDA; e.g., Chen & Armstrong 1975; Webb et al. 1983; Begelman & Kirk 1990). As we describe in the following section, the SDA process is most efficient in shocks where 〈B_y〉_λ/B₀ ≳ 0.015, whereas the DSA mechanism can only operate if the stripe-averaged field is small, 〈B_y〉_λ/B₀ ≲ 0.015.

In this section, we investigate how the physical conditions in the wind can affect the structure of the shock and the physics of magnetic reconnection, with particular focus on the particle spectrum downstream from the hydrodynamic shock. We will explore the dependence of our results on the stripe wavelength (in Section 4.1), the wind magnetization (in Section 4.2), the strength of the stripe-averaged field (in Section 4.3), and the wind bulk Lorentz factor (in Section 4.4).

The 1D model presented by Pétri & Lyubarsky (2007) predicts that magnetic dissipation should not appreciably affect the dynamics of the flow, provided that $\lambda _{\, c/\omega _{p}}/\sigma \gtrsim 8\,\xi _1/5$. Here, $\lambda _{\, c/\omega _{p}}=\lambda /(\, c/\omega _{p})$ is the stripe wavelength in units of the relativistic skin depth in the cold wind, and ξ₁ = 6–10 is the ratio of downstream sheet thickness to downstream Larmor radius (its value is calibrated with 1D PIC simulations). In this case, the striped structure of the wind will be preserved downstream, just compressed, and the flow will satisfy the MHD jump conditions for a strongly magnetized fluid. Alternatively, full dissipation of the alternating fields should occur if $\lambda _{\, c/\omega _{p}}/\sqrt{\sigma }\lesssim 4\,\xi _1/5$, when the thickness of downstream current sheets becomes comparable to the distance between neighboring sheets. In this case, Pétri & Lyubarsky (2007) find that the post-shock flow is purely hydrodynamical, since all the magnetic field energy has been transferred to the particles.

As we discuss below, the results of our 2D and 3D simulations suggest a different picture. We find that complete dissipation of the alternating fields by magnetic reconnection is a common by-product of the shock evolution, even in the cases when the 1D model of Pétri & Lyubarsky (2007) predicts that the stripes should be preserved downstream from the shock. As suggested by Figure 3, multi-dimensional simulations are of key importance to correctly capture the evolution of the tearing-mode instability, which is ultimately responsible for field dissipation and particle acceleration.

Even though complete field dissipation (and transfer of field energy to the particles) occurs irrespective of the wind properties, the shape of the downstream particle spectrum is sensitive to the stripe wavelength and the wind magnetization. In most cases, the particle spectrum can be described as a power law in energy, with a flat slope of index 1 < p < 2. The slope is moderately sensitive to the pre-shock parameters, but it is the width of the spectrum that depends most dramatically on the properties of the wind. If the spectrum can be modeled as a power law of index p between γ_min and γ_max, we can derive a relation between the lower and upper cutoffs by requiring that the mean particle Lorentz factor be 〈γ〉 = γ₀(σ + 1), as expected for full dissipation of magnetic energy. We obtain

$$\begin{eqnarray} &&\left(\frac{\gamma _{{\rm min}}}{\gamma _0}\right)^{p-1} \left(\frac{\gamma _{{\rm max}}}{\gamma _0}\right)^{2-p}=\frac{2-p}{p-1}\,(\sigma +1). \end{eqnarray} \tag{ 2 }$$

For 1 < p < 2 and fixed σ, we see that this expression includes both the case of narrow distributions with γ_min ≃ γ_max ≃ γ₀σ and the case of wide spectra with γ_min ≃ γ₀ and γ_max ≃ γ₀σ^{1/(2 − p)}. As we show in Section 4.1, the stripe wavelength is the main factor, at fixed wind magnetization, that controls the transition from narrow Maxwellian-like spectra to broad power-law tails.

In Figures 5–7, 9 and 12, we show the dependence of the downstream particle spectrum on the conditions in the wind. In most cases, the spectrum is measured between 700 c/ω_p and 300 c/ω_p downstream from the hydrodynamic shock.⁸ In the subpanel of each figure, we estimate the mean particle kinetic energy 〈γ〉/γ₀ (black line) and the width of the particle spectrum, measured as the ratio γ_max/γ_min of upper to lower energy cutoff (red line). More precisely, γ_min will be defined as the Lorentz factor where most of the particles lie (i.e., where γ dN/dγ peaks), whereas γ_max will correspond to the location where most of the energy resides (i.e., where γ²dN/dγ is at maximum).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Dependence on the Stripe Wavelength

4.2. Dependence on the Wind Magnetization

4.3. Dependence on the Stripe-averaged Field

In this section, we explore the dependence of the shock properties on the strength of the stripe-averaged field 〈B_y〉_λ. The magnitude B₀ of the field in the cold wind is fixed by the value of the magnetization σ. A non-zero stripe-averaged field is established if the portion of each stripe with B_y = +B₀ has a different width than the region with B_y = −B₀. If λ₊ is the length of the former and λ₋ of the latter (with the constraint λ₊ + λ₋ = λ), then we have 〈B_y〉_λ/B₀ = (λ₊ − λ₋)/λ. Alternatively, by defining α ≡ (λ₊ − λ₋)/λ₊, we have 〈B_y〉_λ/B₀ = α/(2 − |α|). In pulsar winds, |α| → 1 at latitudes approaching the inclination angle between magnetic and rotational axes of the pulsar (see the upper panel of Figure 1). The special case α = 0, appropriate for the wind midplane, will be discussed in more detail in Section 4.3.1.

The results presented in the previous sections referred to shocks with a moderate value of α = 0.1. As a representative case of high-α flows, Figure 8 shows the shock structure for α = 0.75, corresponding to 〈B_y〉_λ/B₀ ≃ 0.6. The stripe wavelength λ = 320 c/ω_p and the wind magnetization σ = 10 are chosen to facilitate comparison with the results discussed in Section 3.1. As in the case α = 0.1, the pre-shock flow interacts at first with a fast MHD shock (at x ≃ 2170 c/ω_p in Figure 8, indicated by a vertical dotted blue line), that triggers magnetic reconnection. With increasing α (and so, with increasing 〈B_y〉_λ), a smaller fraction of the incoming Poynting flux is in the form of alternating fields, and as such available for dissipation. For α > 0 (which implies λ₊ > λ₋), reconnection behind the fast shock will eventually annihilate the field where B_y = −B₀, but a finite region where the field stays B_y = +B₀ will survive. Dissipation of the alternating component is completed when the size of reconnection islands approaches λ₋ (or, in general, the smallest between λ₊ and λ₋), as shown by the 2D plots of density and magnetic energy in Figures 8(a) and (b), respectively. At this point, a second shock forms (at x ≃ 2050 c/ω_p in Figure 8), which would correspond in the regime of small α (i.e., α ≲ 0.1) to the hydrodynamic shock discussed in Section 3.1. For larger α, the speed of this second shock increases, from β_sh ≃ 1/3, as appropriate for an unmagnetized post-shock plasma, up to β_sh → 1, the shock velocity in a σ ≫ 1 flow. In the limit |α| → 1, the shock that in the regime α ≲ 0.1 was identified as "hydrodynamic" becomes eventually degenerate with the fast MHD shock. In the case α = 0.75 shown in Figure 8, the fast shock stays at a distance of ≃ 100 c/ω_p ahead of the main shock, and they both move at β_sh ≃ 0.6.

Zoom In Zoom Out Reset image size

Downstream from the main shock (x ≲ 2050 c/ω_p), the residual magnetic flux that has not been annihilated on the way from the fast shock survives within stripes of high magnetic field, separated by regions of relatively low magnetic content (see the 2D and 1D plots of magnetic energy in Figures 8(b) and (d), respectively; and the 1D profile of B_y, black line in Figure 8(e)). The plasma inside the high-field stripes is relatively cold, not having gained much energy from field dissipation. In contrast, the low-field regions are populated with hot particles energized by the reconnection field, as shown in the 1D profile of mean kinetic energy per particle of Figure 8(f) and in the longitudinal positron phase space of Figure 8(g). We find that the post-shock flow retains a striped structure, with magnetically dominated regions separated by relatively wide sheets of hot plasma, for values of |α| ≳ 0.3 (or equivalently, for |〈B_y〉_λ|/B₀ ≳ 0.15). As |α| increases, the longitudinal extent of the high-field stripes grows, at the expense of the regions of hot dense plasma, and magnetic pressure dominates over particle pressure in the energetics of the downstream fluid.

The striped structure of the downstream flow leaves a clear imprint on the post-shock particle spectrum, as Figure 9 shows. There, we present how the particle distribution changes when the parameter α is varied between α = 0.0 and α = 0.95, for fixed wavelength λ = 320 c/ω_p, magnetization σ = 10, and bulk Lorentz factor γ₀ = 15.

Zoom In Zoom Out Reset image size

For α ≲ 0.3, the average particle Lorentz factor downstream from the hydrodynamic shock approaches 〈γ〉/γ₀ ≃ σ + 1 ≃ 11, suggesting full dissipation of magnetic fields (black line in the subpanel of Figure 9). The shape of the spectrum is independent of α, in the regime α ≲ 0.1. The spectra for α = 0.0 (red curve) and α = 0.1 (black curve) differ only at the upper end, where the case α = 0.1 shows the high-energy bump of shock-drift-accelerated particles. This component disappears for α = 0.0, when the stripe-averaged field 〈E_z〉_λ = β₀〈B_y〉_λ powering the SDA process vanishes.

For α ≳ 0.3, the post-shock particle spectrum consists of two components. The low-energy peak comes from cold plasma residing in the high-field stripes described above, whereas the high-energy part is populated by hot particles that gained energy from field dissipation. As α increases, the fraction of upstream Poynting flux available for dissipation decreases, which explains why the high-energy component in the spectra of Figure 9 gets de-populated, at the expense of the low-energy part. As a consequence, the average particle energy downstream from the shock monotonically drops from 〈γ〉 ≃ γ₀(σ + 1), corresponding to the case of full dissipation, down to 〈γ〉 ≃ γ₀, as expected for negligible field annihilation (black line in the subpanel of Figure 9). In the limit |α| → 1 (yellow curve for α = 0.95), most of the downstream plasma is contained within the high-field regions, and the low-energy component of the spectrum approaches the result expected for an unstriped wind (purple line).¹²

With increasing α, the high-energy component of the downstream spectrum recedes to lower energies (see the trend in Figure 9, from the blue line for α = 0.45, to the yellow line for α = 0.95). As described in Section 4.1, the maximum energy of particles accelerated at X-points depends primarily on the amount of time available for acceleration, which is set by the characteristic distance between a given X-point and the closest island (in the same current sheet). In turn, this scales with the maximum size of reconnection islands, which is determined by the smallest between λ₊ and λ₋ (here λ₋ < λ₊, since α > 0). For larger α at fixed λ, the extent of λ₋ decreases, so the maximum energy to which particles can be accelerated by the reconnection electric field decreases. This explains the trend in the upper cutoff of the spectra in Figure 9. Given that the location of the high-energy cutoff is set by the value of λ₋ (or equivalently, by the maximum size of reconnection islands), one would expect that simulations with the same λ₋, but different choices of λ and α > 0, should yield similar spectra, at the high-energy end. We have tested this, and we find that flows with the same λ₋ show magnetic islands of similar sizes, and their downstream spectra almost overlap at high energies. Of course, this scaling does not hold for the low-energy spectral component, since it is populated by particles that were initially in the region with B_y = +B₀, and are therefore insensitive to the value of λ₋.

4.3.1. The Case of Zero Stripe-averaged Field

The special case 〈B_y〉_λ = 0 (or equivalently, α = 0), which is relevant for the equatorial plane of pulsar winds, is worth a deeper investigation. As anticipated in Section 4.1, hot downstream particles can propagate back from the shock into the upstream, and cross at least one wavelength of the striped wind, provided that $\lambda _{\, c/\omega _{p}}/\sqrt{\sigma }\lesssim 2\, \xi _2$, where ξ₂ = 1–10 is a numerical factor that depends weakly on λ and σ. For wind parameters that satisfy such criterion, the motion of hot particles moving into the upstream (the "returning" particles, from now on) will be regulated, to zeroth order, by the stripe-averaged fields 〈B_y〉_λ and 〈E_z〉_λ = β₀〈B_y〉_λ. For α ≳ 0.03 (or equivalently, 〈B_y〉_λ/B₀ ≳ 0.015), the returning particles will be quickly advected downstream by the stripe-averaged fields, after being accelerated via the SDA mechanism discussed in Section 4.1. On the other hand, for α ≲ 0.03 they can propagate as if the upstream flow were unmagnetized. The special case α = 0.0 is taken here as a representative example of the shock structure for α ≲ 0.03.

Figure 10 compares the structure of the flow between the cases α = 0.1 (left panels) and α = 0.0 (right panels), for the same wavelength λ = 80 c/ω_p and magnetization σ = 10. The positron phase space of the case α = 0.0 (Figure 10(c), right column) shows a diffuse population of hot particles moving back from the shock into the upstream (x ≳ 1050 c/ω_p), which were absent in the case α = 0.1 (Figure 10(c), left column). For both α = 0.0 and α = 0.1, as the flow propagates from the fast shock (at x ≃ 1050 c/ω_p) to the hydrodynamic shock (at x ≃ 1000 c/ω_p), magnetic reconnection transfers energy from the field to the particles, which in the downstream region end up with an average Lorentz factor 〈γ〉 ≃ γ₀(σ + 1) (see Section 4.1). In other words, the physics of magnetic reconnection is not sensitive to the value of α, in the regime α ≲ 0.1. However, it is only in the case α = 0.0 (or α ≲ 0.03, generally speaking), when the upstream stripe-averaged fields are small, that the hot particles resulting from shock-driven reconnection can stream for large distances ahead of the shock (see the phase space in the right panel of Figure 10(c)).

Zoom In Zoom Out Reset image size

When a sufficient number of returning particles build up in the pre-shock region, their counter-streaming with respect to the incoming flow can trigger the filamentation (or Weibel) instability (Weibel 1959; Medvedev & Loeb 1999; Gruzinov & Waxman 1999). In fact, the filamentary structures seen ahead of the shock in the 2D plot of B_z (panel (d), right column) are generated by the Weibel instability. A similar pattern appears in the 2D plots of electric field E_y (=−β₀B_z), particle density, and magnetic energy (not shown). In contrast, no signatures of filamentation are seen for α = 0.1 (Figure 10(d), left column), since in this case no particles can propagate far enough ahead of the shock to seed the Weibel instability.

For both α = 0.0 and α = 0.1, the hydrodynamic shock satisfies the jump conditions appropriate for an unmagnetized 3D relativistic fluid. However, the filamentation instability appreciably changes the structure of the flow for α = 0.0, especially at late times. The density profile in the transition region (Figure 10(a)) is much smoother in the case α = 0.0 (right column) than it is for α = 0.1 (left column). More importantly, the magnetic fields generated ahead of the shock by the Weibel instability survive in the downstream region of the α = 0.0 shock, as shown by the 1D profile of magnetic energy in Figure 10(b) (right column, at x ≲ 1000 c/ω_p). Instead, for α = 0.1, the only magnetic field left behind the shock comes from compression of the pre-shock stripe-averaged field, which results in a much lower level of post-shock magnetic energy (see the left panel of Figure 10(b)).

The shape of the downstream energy spectrum is also significantly different between the cases α = 0.0 and α = 0.1. Figure 11 follows the time evolution of the particle distribution behind a shock with α = 0.0, compared with the case α = 0.1. As further discussed in Appendix B, the downstream particle spectrum in a flow with α ≳ 0.03 does not appreciably change with time, so the dotted blue line in Figure 11, which refers to ω_pt = 3000 for the case α = 0.1, is a good representation of the downstream energy spectrum at all times (for α = 0.1). In contrast, as shown by the solid lines, the evolution of the particle spectrum in the case α = 0.0 is much more dramatic. At early times (black line for ω_pt = 1500), the spectrum for α = 0.0 is remarkably similar to the (time-independent) particle distribution of the case α = 0.1 (dotted blue line). This confirms once more that the physics of shock-driven reconnection, which is ultimately responsible for shaping the post-shock spectrum at early times, does not appreciably depend on α, in the regime α ≲ 0.1. As mentioned above, the minor high-energy bump seen in the spectrum of α = 0.1 (dotted blue line), but not for α = 0.0 (black line), results from particles accelerated at the hydrodynamic shock via the SDA mechanism.

Zoom In Zoom Out Reset image size

Starting from ω_pt = 3000 (solid blue line), the spectrum for α = 0.0 changes significantly. The high-energy part of the broad peak that was present at earlier times (at γ ≃ 300) gets de-populated, at the expense of a prominent tail of nonthermal particles (at γ ≳ 10³). The tail contains approximately a few percent of particles and a few percent of the total particle energy. With time, the slope of the nonthermal tail becomes flatter (for comparison, the dashed line in Figure 11 has p = 2.5), and its high-energy cutoff extends to higher Lorentz factors, as tracked by the solid line in the subpanel of Figure 11. In contrast, no appreciable evolution in the value of the maximum particle Lorentz factor is observed for the case α = 0.1 (dotted line in the subpanel). The particles in the nonthermal tail seen for α = 0.0 have been accelerated via a kind of first-order Fermi mechanism (or diffusive shock acceleration, DSA), by scattering off the turbulence generated by the Weibel instability (see Figure 10(d), right column). Their energization proceeds in the same way as for a strictly unmagnetized shock, which was discussed by Spitkovsky (2008). By comparing the black line (for ω_pt = 1500) and the blue line (for ω_pt = 3000) at γ ≳ 10², we infer that injection into the Fermi process is primarily due to the so-called thermal leakage, i.e., the hottest particles from the downstream region can escape ahead of the shock, where they take part in the Fermi acceleration cycles. Of course, once the fields produced by the Weibel instability grow stronger than the background field B₀ (see the scale in Figure 10(d), right column), the Weibel-generated turbulence in the shock layer may mediate further injection of particles into the Fermi process.

Finally, we remark that, based on the condition $\lambda _{\, c/\omega _{p}}/\sqrt{\sigma }\lesssim 2\, \xi _2$ discussed above, injection into the Fermi process is expected to be more favorable for short stripe wavelengths and highly magnetized winds. In fact, we observe that the efficiency of Fermi acceleration, meaning the fraction of particles participating in the Fermi cycles, is higher for smaller λ (at fixed σ) and for larger σ (at fixed λ). This explains why, for λ = 320 c/ω_p (so, a wavelength longer than the value λ = 80 c/ω_p used here), Figure 9 shows that the spectra for α = 0.0 and α = 0.1 are almost identical, in contrast with the conclusions of Figure 11. For λ = 320 c/ω_p, fewer particles can be injected into the Fermi process, and the structure of the flow at ω_pt = 3000 is practically the same for all values of α ≲ 0.1.

4.4. Dependence on the Bulk Lorentz Factor

Figure 12 shows the particle spectrum downstream from the hydrodynamic shock, for different bulk Lorentz factors of the striped wind (from γ₀ = 3, in red, up to γ₀ = 375, in blue). The spectra are shifted along the x-axis by 15/γ₀, to facilitate comparison with the reference case γ₀ = 15 (black curve) discussed in the previous sections. Once normalized in this way, the spectra overlap almost perfectly, suggesting that the physics of shock-driven reconnection is not sensitive to the upstream bulk Lorentz factor. The value of γ₀ enters the definition of the relativistic plasma skin depth c/ω_p ($\propto \sqrt{\gamma _0}$) and of the relativistic Larmor radius ($r_L=\, c/\omega _{p}/\sqrt{\sigma }\propto \sqrt{\gamma _0}$), but no residual dependence on γ₀ remains if length and timescales are normalized with respect to c/ω_p and ω⁻¹_p, respectively. More specifically, the conditions that control the physics of shock-driven reconnection (see Section 4.1 and Pétri & Lyubarsky 2007) depend only on the stripe wavelength λ and the magnetization σ, but they do not explicitly depend on the wind Lorentz factor. This explains why our results are nearly the same (modulo an overall shift in energy scale) across a wide range of γ₀.

Zoom In Zoom Out Reset image size

The results presented so far are based on 2D experiments with magnetic field initialized in the plane of the simulations. This is the best configuration to study the tearing-mode instability, whose wavevector lies in the plane of the field, but it artificially inhibits the growth of cross-plane instabilities, which may also be important for the structure of the shock. In this section, we study the 3D physics of shock-driven reconnection, with particular emphasis on the process of particle energization. The simulation setup parallels very closely what we described in Section 2, with the magnetic field oriented initially along the y direction. We adopt our fiducial values for σ = 10 and α = 0.1, and we investigate two choices for the stripe wavelength, λ = 80 c/ω_p and λ = 160 c/ω_p. For each value of λ, the extent of the simulation domain along y and z is chosen such that our results are not artificially affected by the periodicity of our boundary conditions (see Appendix A). We adopt L_y = L_z = 50 c/ω_p for λ = 80 c/ω_p, and L_y = L_z = 100 c/ω_p for λ = 160 c/ω_p.

Figure 13 shows the structure of the shock at ω_pt = 990, for a striped wind of wavelength λ = 160 c/ω_p. In agreement with the results presented in Section 3.1 for 2D simulations, the structure of the flow in Figure 13 shows the presence of two shocks, a fast MHD shock located at x ≃ 630 c/ω_p and a hydrodynamic shock at x ≃ 250 c/ω_p. At the fast shock, the flow is decelerated and compressed, as shown by the increase in magnetic energy in Figure 13(b). The passage of the fast shock through the incoming current sheets initiates magnetic field annihilation. The dissipation of alternating fields proceeds as the flow propagates from the fast to the hydrodynamic shock, and little magnetic energy remains downstream from the hydrodynamic shock (see Figure 13(b) at x ≲ 250 c/ω_p), in agreement with our 2D results. As observed in 2D, the particle density downstream from the hydrodynamic shock (Figure 13(a)) and the shock velocity are in agreement with MHD jump conditions for a relativistic unmagnetized fluid. In fact, most of the post-shock energy is in kinetic form, rather than electromagnetic form.¹³

Zoom In Zoom Out Reset image size

As the flow propagates from the fast to the hydrodynamic shock, the structure of the current sheets is affected by the growth of two competing modes. In the x–y plane of the magnetic field, the tearing-mode instability breaks the current sheet into a sequence of high-density islands (see Figure 13(a)), separated by X-points where the field lines tear and reconnect. In the x–z plane orthogonal to the field, the folding of the current sheet seen in Figure 13 is governed by the drift-kink instability, driven by the current of fast-drifting plasmas in a thin sheet (Daughton 1998; Zenitani & Hoshino 2005a, 2007). The 2D geometry employed in the previous sections was chosen to select the tearing mode and suppress the drift-kink mode. In the regime of undriven reconnection, the drift-kink instability is generally thought to grow faster than the tearing mode (Zenitani & Hoshino 2005b, 2007), but for very thin sheets the opposite hierarchy is observed, as shown via 3D PIC simulations by Liu et al. (2011). In our experiments of shock-driven reconnection, the current sheet just behind the fast shock (x ≃ 580 c/ω_p) shows the pattern of magnetic islands typical of the tearing mode, whereas the drift-kink mode appears only farther downstream (x ≲ 500 c/ω_p).

The nature of the instability that dominates the process of field annihilation leaves an imprint on the resulting particle distribution. Zenitani & Hoshino (2007) have shown with PIC simulations that field dissipation due to the drift-kink instability does not result in nonthermal particle acceleration, the plasma just being heated. This is because the field lines remain straight, so that all particles gain energy at the same rate. Nonthermal particles are produced only by the tearing mode (e.g., Zenitani & Hoshino 2001). In our 2D simulations with in-plane fields, i.e., the geometry required to capture the tearing mode, we have shown that particles are accelerated to nonthermal energies by the reconnection electric field. They initially move along the z direction parallel to the reconnection field, and then they drift along y from a given X-point into the closest magnetic island. In 3D, one could argue that the folding of the current sheet introduced by the drift-kink instability in the x–z plane (not resolved in 2D) may deflect the accelerating particles out of the current sheet, thus suppressing their energization.

In Figure 14, we address this important issue, by comparing for λ = 80 c/ω_p the post-shock particle spectrum of a 3D simulation (red line) to the results of 2D experiments with in-plane or out-of-plane magnetic fields (blue and green lines, respectively). The excellent agreement between 3D and 2D in-plane results (red and blue lines, respectively) suggests that the physics of particle acceleration by shock-driven reconnection is captured extremely well by our 2D simulations with in-plane fields.¹⁴ In turn, this implies that, at least for λ = 80 c/ω_p, the maximum energy of accelerated particles is constrained by the distance between X-points and islands (in the x–y plane), rather than by the current sheet folding in the x–z plane. Just upstream of the hydrodynamic shock, the typical spacing between X-points and islands scales linearly with the stripe wavelength. A similar trend may be expected in the characteristic scale of the drift-kink mode, for the current sheet just ahead of the hydrodynamic shock (see Figure 13(b) at x ≃ 350 c/ω_p). If this is the case, our conclusions will also hold for longer stripe wavelengths, as we have directly verified in the case λ = 160 c/ω_p. A detailed 3D study of particle acceleration by shock-driven reconnection will be presented elsewhere.

Zoom In Zoom Out Reset image size

Finally, we point out that the tension force of a guide field (along z, in our geometry) can easily stabilize the drift-kink mode, with little or no effect on the tearing mode (Zenitani & Hoshino 2008). In the presence of a guide field, the 3D physics of shock-driven reconnection should be described very accurately by 2D simulations with alternating fields lying in the simulation plane. A guide field component, though, is not usually expected in the context of pulsar winds (but see Pétri & Kirk 2005).

We have explored by means of 2D and 3D PIC simulations the internal structure and acceleration properties of relativistic shocks that propagate in an electron–positron striped wind, i.e., a flow consisting of stripes of alternating field polarity separated by current sheets of hot plasma. This work extends the study of Sironi & Spitkovsky (2009), which investigated relativistic shocks in a pair plasma carrying uniform magnetic fields.

We find that a fast MHD shock propagates into the pre-shock striped flow, compressing the incoming current sheets and initiating the process of driven magnetic reconnection, via the tearing-mode instability. Reconnection islands seeded by the passage of the fast shock grow and coalesce, while magnetic energy is dissipated at X-points located in between each pair of islands. When reconnection islands grow so big to occupy the entire region between neighboring current sheets, the striped structure of the flow is erased, and a hydrodynamic shock forms. Downstream from the shock, the average particle energy is larger than in the pre-shock flow by a factor of ≃ σ, where σ is the wind magnetization. In other words, the energy stored in the alternating fields has been entirely transferred to the particles via shock-driven magnetic reconnection, and the downstream fluid behaves like an unmagnetized plasma. The only field component surviving in the downstream region comes from shock-compression of the stripe-averaged pre-shock field (i.e., the non-alternating component), if non-zero. In our 2D and 3D simulations, we find that complete annihilation of the alternating fields occurs irrespective of the stripe wavelength λ or the magnetization σ, in contrast with the results of the 1D analysis by Pétri & Lyubarsky (2007). In Appendix A we show that multi-dimensional studies are of paramount importance to correctly capture the physics of shock-driven reconnection.

The main agent of particle energization is the reconnection electric field, as particles drift from a given X-point into the closest island. Whether all particles have comparable energy gains, or only a few particles are accelerated to high energies, and the majority stay cold, depends sensitively on the properties of the wind. We have explored the dependence of the downstream particle spectrum on the stripe wavelength λ, the wind magnetization σ (or equivalently, the field strength B₀), and the stripe-averaged magnetic field 〈B_y〉_λ (or equivalently, the parameter α, with α = 0 for 〈B_y〉_λ = 0, and |α| → 1 in the limit |〈B_y〉_λ| → B₀). For fixed α = 0.1, we find that the shape of the post-shock spectrum depends primarily on the combination λ/(r_Lσ), where r_L is the relativistic Larmor radius in the striped wind. This is just the stripe wavelength measured in units of the post-shock plasma skin depth. For small values of λ/(r_Lσ) (≲ a few tens), the spectrum resembles a Maxwellian distribution with mean energy 〈γ〉 ≃ γ₀σ, where γ₀ is the bulk Lorentz factor of the pre-shock flow. In the limit of very large values of λ/(r_Lσ) (≳ a few hundreds), the spectrum approaches a broad power-law tail, with a flat spectral slope p ≃ 1.5. In this regime, the tail extends from γ_min ≃ γ₀ up to γ_max ≃ γ₀σ^{1/(2 − p)}. In terms of the flow structure, the former case (i.e., a Maxwellian-like spectrum) is realized when all particles pass in the vicinity of an X-point, thus gaining energy from field dissipation. In contrast, the latter case occurs when most of the particles remain far from X-points, thus staying cold, and the energy of annihilating fields is transferred to only a few particles, which end up dominating the energy content of the flow.

As |α| increases from 0.1 up to 1, the fraction of upstream Poynting flux in the form of alternating fields becomes smaller. As the amount of field energy available for dissipation decreases, the average particle Lorentz factor downstream from the shock recedes from 〈γ〉 ≃ γ₀σ, as expected for complete field annihilation, down to 〈γ〉 ≃ γ₀, the result of an unstriped wind. In the opposite limit |α| ≲ 0.1, i.e., in the case of nearly symmetric stripes, the particles accelerated by the reconnection electric field to the highest energies can escape ahead of the shock, unaffected by the stripe-averaged field that would tend to advect them back downstream. Their counter-streaming with respect to the incoming flow can seed the filamentation (or Weibel) instability, as it happens for shocks propagating in unmagnetized plasmas (Spitkovsky 2005, 2008). In turn, the turbulence generated by the Weibel instability can mediate particle acceleration to even higher energies, via a Fermi-like diffusive process. The nonthermal tail produced by the Fermi mechanism is steeper than that given by magnetic reconnection, with characteristic slopes p ≳ 2.5. Typically, a few percent of particles are injected into the Fermi process, and the injection efficiency scales in inverse proportion to the ratio λ/(r_Lσ). With time, the nonthermal tail of Fermi-accelerated particles extends to higher energies. For |α| ≳ 0.1, the propagation of particles far ahead of the shock is inhibited by the stripe-averaged field, and Fermi acceleration is suppressed. In this case, a few percent of particles can get accelerated by the stripe-averaged electric field 〈E_z〉_λ while gyrating around the shock front, a process known as shock-drift acceleration (i.e., Begelman & Kirk 1990). The tail resulting from shock-drift acceleration has a limited energy extent, and it does not evolve in time.

We have confirmed our findings with an extensive set of convergence tests. Most importantly, as we have shown in Section 5, 2D simulations with in-plane magnetic fields—the setup that we have primarily adopted in this work—can reproduce very accurately the shock physics and particle spectrum observed in 3D simulations. In other words, the current sheet folding caused by the drift-kink instability in the plane orthogonal to the field (Zenitani & Hoshino 2007) does not represent a serious constraint for the maximum particle energy achieved in the course of the reconnection process. Rather, the maximum energy of accelerated particles is limited by the time to drift from a given X-point into the closest island.

These findings can place important constraints on current models of PWNe, bubbles of synchrotron-emitting plasma powered by the relativistic wind of young pulsars. If the magnetic and rotation axes of the pulsar are not aligned, the wind consists of stripes of opposite magnetic field polarity, alternating with the pulsar period. The wind wavelength will be λ = 2πR_LC, where R_LC = c/Ω is the so-called light-cylinder radius and Ω is the angular frequency of the pulsar. The stripe-averaged field vanishes along the equatorial plane (i.e., α = 0 there), and the striped structure persists up to a latitude equal to the inclination angle between the pulsar's magnetic and rotational axes. At higher latitudes, the wind carries a uniform (i.e., non-alternating) field, so the physics discussed in the present work does not apply (see, instead, Gallant et al. 1992; Hoshino et al. 1992; Amato & Arons 2006; Sironi & Spitkovsky 2009, 2011).

Our results suggest that efficient dissipation of the alternating fields should occur at the termination shock of pulsar winds, irrespective of the properties of the flow. It follows that current observational constraints on the upstream magnetization parameter (Rees & Gunn 1974; Kennel & Coroniti 1984a; Del Zanna et al. 2004) should be attributed not to the total Poynting flux, but to the Poynting flux associated with the stripe-averaged magnetic field, the only component surviving downstream from the shock. The upstream flow may be Poynting-dominated, provided that most of the Poynting flux is in the form of an alternating field, which is annihilated at the shock. Thus, our results provide a solution to the so-called sigma problem, i.e., that pulsar winds should be magnetically dominated at the light-cylinder radius, but kinetically dominated downstream from the termination shock. The idea of shock-driven reconnection as a solution to the sigma problem was originally proposed by Lyubarsky (2003). Our work demonstrates that complete annihilation of the alternating fields occurs regardless of the properties of the wind, in the regime of relativistic magnetically dominated flows.

In addition, the particle energy spectra extracted from our simulations can be used directly to interpret the radiative signature of PWNe. The radio spectrum of the Crab Nebula, the prototype of the class of PWNe, requires a population of nonthermal particles with a flat spectral slope (p ≃ 1.5), extending at least across three decades in energy, from 10² MeV up to 10⁵ MeV. The particle spectrum should be steeper at higher energies, with p ≃ 2.5, to explain the optical and X-ray flux.

Our results suggest that, for a slope p ≃ 1.5, the width of the electron power-law tail resulting from reconnection will be ≃ σ². To produce a particle distribution extending over three decades in energy, the wind magnetization at the termination shock needs to be σ ≳ 30. Most importantly, broad particle spectra are produced by shock-driven reconnection only if λ/(r_Lσ) ≳ a few tens. For pulsar winds, the particle number density scales with distance from the pulsar as n = n_LC(R_LC/R)², where the density at the light-cylinder radius is usually written as n_LC = κ n_GJ. Here, n_GJ = ΩB_LC/2πec is the Goldreich–Julian density (Goldreich & Julian 1969), and κ is the so-called multiplicity. Conservation of energy along the flow streamlines implies that at the termination shock γ₀(1 + σ)κ = ω_LC/2 Ω, where ω_LC = eB_LC/mc. It follows that at the termination shock (R = R_TS)

$$\begin{eqnarray} &&\frac{\lambda }{r_L \sigma }\simeq 4\pi \kappa \frac{R_{{\rm LC}}}{R_{{\rm TS}}}. \end{eqnarray} \tag{ 3 }$$

The constraint λ/(r_Lσ) ≳ a few tens required to produce broad particle spectra is extremely challenging for current models of PWNe and pulsar magnetospheres. The radius of the Crab termination shock in the equatorial plane is inferred from X-ray observations to be R_TS ≃ 0.1 pc ≃ 5 × 10⁸R_LC (Hester et al. 2002). Most available models estimate κ ≃ 10⁴–10⁶ (Bucciantini et al. 2011), depending on whether or not radio-emitting electrons are included in the analysis. Based on our findings, the resulting value of λ/(r_Lσ) ≲ 0.01 would yield a Maxwellian-like spectrum, at odds with the wide flat spectrum required by observations.

If radio-emitting electrons are produced in the equatorial plane by shock-driven reconnection, a revision of the existing estimates of κ is required. In this respect, we point out that the values of κ quoted above are averages over latitude, and one cannot exclude that particle injection into the pulsar wind is highly anisotropic, with multiplicity as large as 10⁸ along the equatorial plane. Of course, any given choice of κ will constrain γ₀ and σ via γ₀(1 + σ)κ = ω_LC/2 Ω, a relation that holds at each latitude. From the equatorial plane, radio-emitting electrons will quickly fill the whole nebula, transported by the strong fluid motions observed in MHD models of PWNe downstream from the termination shock (e.g., Komissarov & Lyubarsky 2004; Del Zanna et al. 2004; Camus et al. 2009). This would provide a natural explanation for the lack of gradients in the radio spectral slope of the Crab (Bietenholz & Kronberg 1992). Alternatively, the radio part of the spectrum may be produced at higher latitudes, where the termination shock is closer to the pulsar, and the ratio in Equation (3) becomes larger. However, in this case one should also take into account the fact that, at high latitudes, the fraction of upstream Poynting flux available for dissipation is lower (since |α| is larger), which results in a narrower spectrum, everything else being fixed.¹⁵

On the other hand, the small value of λ/(r_Lσ) expected in the equatorial plane on the basis of current models of the Crab Nebula is the most favorable setup for particle acceleration via the Fermi diffusive process, with efficiency of a few percent by number. In this respect, one could interpret the optical and X-ray signature of the Crab, which requires a particle spectrum with p ≃ 2.5, as synchrotron emission from such Fermi-accelerated particles. Based on our results, optical and X-ray emitting electrons should be produced close to the equatorial plane of the wind, where |α| ≲ 0.1. If the same equatorial wedge is responsible for accelerating the radio-emitting electrons, our findings predict that they should outnumber the optical and X-ray-emitting particles by a factor of few hundred. This is in agreement with current models of PWNe spectra, which require κ ≃ 10⁶ if the radio band is included in the analysis (Bucciantini et al. 2011), a value which is two orders of magnitude larger than if only higher frequency bands are considered (κ ≃ 10⁴; Kennel & Coroniti 1984b).

Although our work focuses primarily on the physics of PWNe, shocks propagating into striped flows may be present also in blazar jets and gamma-ray bursts (e.g., Thompson 2006), where they may provide an appealing explanation, based on strong physical grounds, for any flat particle distribution inferred from observations. The broadband emission of hotspots in radio galaxies, usually interpreted as the termination shocks of relativistic jets, indicates that the spectra of accelerated electrons need to be flat (1 < p < 2) below GeV energies (Stawarz et al. 2007). In the case of luminous blazar sources, very flat electron spectra below GeV energies are inferred directly from X-ray observations (Celotti & Ghisellini 2008; Sikora et al. 2009). The steeper particle spectrum required at ≳GeV energies, with slope p ≃ 2.5, could instead result from Fermi acceleration across the shock, if the stripe-averaged field satisfies |α| ≲ 0.1.

We thank J. Arons, R. Blandford, N. Bucciantini, A. Loeb, J. McKinney, R. Narayan, and E. Quataert for insightful comments. L.S. gratefully thanks D. Giannios for many inspiring discussions. This work was supported by NSF grant AST-0807381 and NASA grant NNX10AI19G. The simulations presented in this article were performed on computational resources supported by the PICSciE-OIT High Performance Computing Center and Visualization Laboratory, and at National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US Department of Energy under contract No. DE-AC02-05CH11231.

In Figure 15 we investigate the dependence of our findings on the transverse size L_y of the computational domain, for λ = 640 c/ω_p, σ = 10, and α = 0.1. The agreement between the black curve (L_y = 400 c/ω_p) and the red line (L_y = 800 c/ω_p) suggests that our results become insensitive to the transverse size of the box (i.e., they converge with respect to L_y), for boxes larger than 400 c/ω_p. Generally speaking, we find that a good criterion for convergence is L_y ≳ λ/2. When this condition is fulfilled, the growth of magnetic islands behind the fast shock is not artificially inhibited by the periodicity of our boundaries in the y direction. Magnetic reconnection then proceeds up to the point when islands from two neighboring current sheets will merge, which happens when their size is ≃ λ/2 (for α ≲ 0.1).

Zoom In Zoom Out Reset image size

On the other hand, for L_y ≲ 300 c/ω_p, the growth and coalescence of reconnection islands within a given current sheet artificially stops when only one island is left in the sheet, and its size approaches L_y. At this point, a second shock forms, located behind the fast shock. For L_y ≳ 400 c/ω_p, this would correspond to the hydrodynamic shock discussed in Section 3.1, where the striped structure of the flow is erased, and field energy is entirely transferred to the particles. Instead, in the case L_y ≲ 300 c/ω_p, complete dissipation of the alternating fields is prohibited by the fact that the growth of islands is artificially inhibited by the transverse extent of the box. Here, the post-shock flow retains a striped structure, with regions of hot plasma separated by highly magnetized walls of cold particles. With decreasing L_y, less magnetic energy is transferred to the particles behind the fast shock, and the post-shock flow gets more dominated by magnetic field energy, with respect to particle kinetic energy (see the decrease in average particle energy with decreasing L_y, in the subpanel of Figure 15). This also explains why in Figure 15 the normalization of the low-energy peak, populated by the cold particles in the high-field regions, grows with decreasing L_y, at the expense of the high-energy component of hot particles that gained energy from field dissipation.

As discussed in the main body of the paper (see Section 4.1), the maximum size of magnetic islands (which here is set by the smallest between λ/2 and L_y) correlates with the highest energy at which particles can be accelerated by the reconnection electric field. The trend in the upper spectral cutoff of Figure 15 can then be explained by the decrease with L_y in the maximum size of reconnection islands. In this respect, the spectra in Figure 15, interpreted as a monotonic trend in the size of reconnection islands, can be directly compared to the spatial sequence of spectra in Figure 3(f), which follow the growth of magnetic islands as the flow propagates from the fast to the hydrodynamic shock. A similar trend as in Figure 15 is also observed in the large-α limit of Figure 9 (yet for a different wavelength, λ = 320 c/ω_p as opposed to λ = 640 c/ω_p). However, in that case the maximum size of magnetic islands was set by a physical constraint, namely the smallest between λ₊ (where B_y = +B₀) and λ₋ (where B_y = −B₀), whereas here it is the transverse size of the box which artificially limits the growth of magnetic islands for L_y ≲ 300 c/ω_p.

As the downstream region gets more magnetized with decreasing L_y, the main shock, which would be moving at β_sh ≃ 1/3 for L_y ≳ 400 c/ω_p, propagates at a faster velocity, eventually catching up with the fast MHD shock, in the limit L_y ≪ λ. In this regime, we recover the 1D results by Pétri & Lyubarsky (2007). The wind parameters employed in this section, λ = 640 c/ω_p and σ = 10, satisfy the condition $\lambda _{\, c/\omega _{p}}/\sigma \gtrsim 8\,\xi _1/5$ (where ξ₁ = 6–10). In this case, the 1D model by Pétri & Lyubarsky (2007) would predict negligible field dissipation, in agreement with our results for a very narrow box (yellow line for L_y = 12.5 c/ω_p; see also the subpanel, for small values of L_y). However, as clarified by Figure 15, this is just an artificial consequence of the reduced dimensionality of their model, which cannot correctly capture the development of the tearing-mode instability, and the resulting growth and coalescence of magnetic islands. In a more realistic 2D scenario, complete field dissipation is achieved by the time the flow enters the hydrodynamic shock. This stresses that multi-dimensional simulations are essential for our understanding of shock-driven reconnection.

In Figure 16 we follow the time evolution of the shock, from the earliest stages until it reaches a steady state. We choose our fiducial values for λ = 640 c/ω_p, σ = 10, and 〈B_y〉_λ ≃ 0.05 (corresponding to α = 0.1), but our main conclusions hold basically for the whole parameter space we have investigated. At very early times (uppermost row, for ω_pt = 750), only the fast MHD shock is present (its location at different times is shown as a vertical dotted blue line in Figure 16). The incoming flow decelerates at the fast shock, forming a ring in momentum space lying in the x–z plane orthogonal to the magnetic field. Since the spread in y-momentum is much smaller than in the other two directions, the fluid behind the fast shock can be treated as a 2D relativistic plasma. In fact, the increase in density, magnetic energy, and transverse magnetic field across the fast shock, as well as the drop in average kinetic energy per particle, is in good agreement with the MHD jump conditions for a 2D relativistic fluid.

Zoom In Zoom Out Reset image size

Behind the fast shock, the energy content of the flow at ω_pt = 750 is dominated by magnetic fields. However, even at such early times, we see that current sheets are getting wider (left panel) and populated with hotter particles (see the mean kinetic energy per particle, black line in the right panel), as we proceed farther downstream from the fast shock. As described in Section 3, a larger and larger fraction of the cold wind is being channeled into reconnection islands, seeded by the passage of the fast shock through the upstream current sheets. At X-points located in between each pair of neighboring islands, magnetic energy is being transferred to the particles, which explains why the peaks in mean kinetic energy per particle are getting higher with distance behind the fast shock (black line in the right panel), at the expense of the mean magnetic energy per particle (red line). Yet, the striped structure, with magnetized regions of cold plasma separated by hot current sheets, still persists downstream from the fast shock.

The shock structure changes significantly at ω_pt = 2250 (second row from the top). Behind the fast shock, which is now much weaker (see the modest density jump at x ≃ 1850 c/ω_p), reconnection islands have enough time to grow and fill the entire region in between neighboring current sheets. At x ≃ 700 c/ω_p, the striped structure is erased, and a hydrodynamic shock forms (see the density profile in the left panel). The density plateau behind the hydrodynamic shock is in agreement with the expected jump conditions for a relativistic unmagnetized 3D plasma, namely n_d/n_u ≃ 4. In fact, the energy balance presented in the right panel shows that, by the time the flow crosses the hydrodynamic shock, most of the field energy (red line) has been converted to particle kinetic energy (black line).

A similar picture holds for later times (ω_pt = 3750 in the third row from the top, ω_pt = 5250 in the lowermost row). The density behind the hydrodynamic shock stabilizes at n_d ≃ 4 n_u, and the transfer of energy from magnetic to kinetic form becomes more and more complete. The change of adiabatic index due to the transition from a magnetically dominated fluid to a kinetically dominated plasma results in a lower pressure to drive the fast shock, which decelerates (compare the locations of the vertical dotted blue lines at different times). Between ω_pt = 3750 and ω_pt = 5250, the overall structure of the transition region reaches a steady state. The fast shock propagates ahead of the hydrodynamic shock, roughly at a fixed distance (i.e., the speed of the fast shock approaches the velocity β_sh ≃ 1/3 of the hydrodynamic shock). The separation between the two shocks obeys the requirement that reconnection islands seeded by the fast shock should fill the entire region in between current sheets, by the time the flow arrives at the hydrodynamic shock. When the flow structure has reached the steady state, we find that the particle spectrum downstream from the hydrodynamic shock does not significantly change with time.