© B. Cerutti and G. Giacinti 2020

1. Introduction

Pulsar wind nebulae are archetypal cosmic particle accelerators. The most studied amongst them, the Crab Nebula, presents one of the best known examples of a purely nonthermal emission spectrum extending over 20 orders of magnitude in frequency range, from 100 MHz radio waves to 100 TeV gamma rays (Meyer et al. 2010). The bulk of the emission is almost certainly of synchrotron origin and extends up to the synchrotron burnoff limit, namely ≲100 MeV (de Jager et al. 1996), and slightly beyond during gamma-ray flares (Abdo et al. 2011; Tavani et al. 2011). The electron spectrum is a broad power-law distribution spreading over at least eight orders of magnitude in particle Lorentz factor, 10 ≲ γ ≲ 10⁹, with a major spectral break about half way, γ  ∼  10⁵. Below this break, the power law is hard and is responsible for the radio to infrared emission. Above the break, the spectrum steepens significantly and forms the optical to 100 MeV emission. Within the classical models of Rees & Gunn (1974) and Kennel & Coroniti (1984), the break as well as the high-energy component are interpreted as the injection of electron–positron pairs by the pulsar wind which are then accelerated at the wind termination shock front. This scenario is all the more promising as the slope of the injected particles above the break coincides with the first-order Fermi acceleration prediction, that is, dN/dγ ∝ γ^−2.2 (Bednarz & Ostrowski 1998; Kirk et al. 2000; Achterberg et al. 2001; Pelletier et al. 2017).

Nevertheless, particle acceleration is dramatically suppressed in the presence of a mean magnetic field transverse to the shock normal (Langdon et al. 1988; Begelman & Kirk 1990; Gallant et al. 1992), as found in pulsar wind nebulae where the magnetic field structure is mostly toroidal. If the plasma magnetization parameter σ, defined as the magnetic to particle enthalpy density ratio, is σ ≳ 10⁻², particles are unable to return to the shock front. Therefore, plasma turbulence is too weak to scatter particles back and forth multiple times across the shock as needed for the first-order Fermi process to operate (Lemoine & Pelletier 2010; Sironi et al. 2013; Plotnikov et al. 2018). Global three-dimensional (3D) magnetohydrodynamic (MHD) simulations of the Crab Nebula favor a mean plasma magnetization of order unity which can locally reach up to σ ≈ 10 at high latitudes (Porth et al. 2014). Thus, Fermi acceleration should be quenched, while at the same time these simulations indicate that particle acceleration most likely occurs within the equatorial regions of the shock front (Porth et al. 2014; Olmi et al. 2015).

The conclusion that magnetized relativistic shocks do not accelerate particles implicitly relies on the assumption that plasma turbulence must be self-generated within the flow, like in unmagnetized shocks where the Weibel instability seeds plasma turbulence (Medvedev & Loeb 1999; Spitkovsky 2008). This may not be the case and this is why recent studies have looked for an externally driven source of plasma turbulence, such as for example corrugations in the shock front (Lemoine 2016; Demidem et al. 2018) or the large-scale dynamics of the nebula driven by magnetic pitching and current-driven instabilities (Begelman 1998; Porth et al. 2014). Another possible solution is to consider the equatorial belt where the magnetic field vanishes by symmetry, and therefore where the Fermi process could operate (Giacinti & Kirk 2018).

The idea of driven magnetic reconnection within the large-scale pulsar wind current sheet at the shock front has also been considered to circumvent the above difficulties (Lyubarsky 2003; Pétri & Lyubarsky 2007), but this scenario requires an unusually high pair plasma supply (Sironi & Spitkovsky 2011) and assumes that negligible dissipation took place in the current sheet before the shock, which may not happen (Coroniti 1990; Cerutti & Philippov 2017). Another particle acceleration mechanism involves electron acceleration by the absorption of ion cyclotron waves emitted at the shock front (Hoshino et al. 1992; Amato & Arons 2006), but this model requires a high injection rate of ions in the wind. All things considered, the origin of particle acceleration in pulsar wind nebulae remains elusive (see reviews by Kirk et al. 2009; Amato 2020).

In this work, we revisit the model of particle acceleration in relativistic magnetized shocks, taking into account a realistic latitudinal dependence of the transverse magnetic field at the shock front as expected from the theory of pulsar winds (Michel 1973; Bogovalov 1999), in contrast with previous models which assume a uniform field. In essence, we extend the model proposed by Giacinti & Kirk (2018) to a larger latitudinal extent and use two-dimensional (2D) ab-initio particle-in-cell (PIC) simulations. Here, we focus our attention on the X-ray- and < 100 MeV gamma-ray-emitting electrons only. The remainder of the paper is organized as follows. Section 2 describes the physical model, the numerical setup, and the list of runs performed in this study. Simulation results are presented in Sects. 3–5, which are further discussed in Sect. 6, with particular emphasis on the Crab Nebula.

2. Numerical setup

Our setup is inspired from previous PIC simulations of relativistic shocks (Spitkovsky 2008; Sironi et al. 2013; Plotnikov et al. 2018). It is a Cartesian box initially filled with an ultra-relativistic cold and magnetized beam of electron–positron pairs propagating along the +x-direction, which mimics the radial direction in this case. The right boundary reflects the particles and the fields with no loss of energy in order to form two counter-streaming beams, which eventually leads to the formation of the shock. The key difference with previous studies is the anisotropic transverse field profile along the shock front, here along the y-direction, which mimics the latitude. This new setup leads to additional numerical complications that we describe below.

2.1. Fields

In a split-monopole configuration (Bogovalov 1999), an oscillatory current sheet forms within the wind at the interface between the two magnetic polarities and fills a spherical wedge in the equatorial regions. The wind zone containing the sheet is called the striped wind (Coroniti 1990; Kirk et al. 2009). Assuming that the sheet has fully dissipated before the wind enters the shock, we are left with the DC and axisymmetric component of the toroidal magnetic field, B_ϕ, and therefore the problem which was initially three-dimensional becomes essentially two-dimensional. Translated into Cartesian coordinates and assuming that the neutron star angular velocity and magnetic field vectors fulfill Ω ⋅ B >  0, a good proxy for the out-of-plane magnetic field profile is¹

$$\begin{array}{c}\hfill B_z(y)=B_{0}tanh\left(\frac{y}{L_{\mathrm{s}}}\right)sin\theta ,\end{array}$$(1)

where the z-direction plays the role of the toroidal direction, B₀ is the fiducial magnetic field strength in the upstream medium, and L_s is the spatial extent of the striped wind region set by the inclination angle between the magnetic axis and the pulsar spin axis, such that χ = πL_s/2L_y. The angle θ = π(y + L_y)/2L_y is the polar angle such that y = 0 (θ = π/2) should be understood as the equatorial plane, while y = ±L_y should be interpreted as the poles (θ = 0,π). These three places all have in common that the field vanishes exactly; this plays an important role in the following. Nevertheless, one should keep in mind that we neglect the curvature of the shock front with this Cartesian setup. Figure 1 (top panel) shows the dependence of B_z for a nearly aligned pulsar with L_s/L_y = 0.1 (χ = 9°), L_s/L_y = 0.5 (χ = 45°), and for an orthogonal rotator L_s/L_y = 1 (χ = 90°).

Fig. 1. Top: transverse magnetic field profile in the pulsar wind along the y-direction, Bz(y), according to Eq. (1) for Ls/Ly = 0.1, 0.5, 1. Bottom: corresponding electric current profile, Jx(y).

The background electric field is the ideal advection field

$$\begin{array}{c}\hfill \mathit{E}=-\frac{\mathit{V}\times \mathit{B}}{c},\end{array}$$(2)

where V is the bulk plasma velocity, and c the speed of light. For simplicity, we assume there is no latitudinal dependence of the wind velocity, V(y) = V₀, and so E_y(y) = V₀B_z(y)/c. The dimensions of the box in each direction are x ∈ [0, L_x] and y ∈ [ − L_y, L_y]. For numerical convenience, we apply periodic boundary conditions for both the fields and the particles along the y-directions. This choice has also a physical motivation because the toroidal field direction changes sign across the rotation axis.

2.2. Current and charge densities

Electron–positron pairs are continuously injected throughout the duration of the simulation. To save on computing time, new particles are uniformly created by an injector receding at the speed of light away from the right boundary, which is initially located at x = 0.95L_x. According to Ampère’s law, the pulsar wind must carry the following electric current density, 4πJ_x = cdB_z/dy, and therefore

$$\begin{array}{c}\hfill J_x=\frac{c{B}_0}{4\pi }[\frac{sin\theta }{L_{\mathrm{s}}{cosh}^2(y/L_{\mathrm{s}})}+\frac{\pi }{2L_y}tanh\left(\frac{y}{L_{\mathrm{s}}}\right)cos\theta ]=J_{+}+J_{-}.\end{array}$$(3)

The first term, J₊, is strictly positive; it peaks at the equator (y = 0) and vanishes at the poles (y = ±L_y), while the second term, J₋, is strictly negative with a minimum at the poles and vanishes at the equator (see the total current profile in Fig. 1, bottom panel). The total electric current passing through the yz-plane is zero as expected for a steady state wind. According to Gauss’s law, the motion electric field in Eq. (2) leads to the following distribution of electric charges in the wind

$$\begin{array}{c}\hfill \rho =\frac{1}{4\pi }\frac{\mathrm{d}{E}_y}{\mathrm{d}y}=\frac{V_0}{c^2}{J}_x,\end{array}$$(4)

so that J_x ≈ ρc in the ultrarelativistic limit.

2.3. Plasma density

Assuming that both species move along the +x-direction at the same speed V₀, the minimum amount of pairs required to model both the current and the charge densities in the wind is to consider the positronic density profile

$$\begin{array}{c}\hfill n_{+}=\frac{V_0{J}_{+}}{e{c}^2}=n_0^{+}\frac{sin\theta }{{cosh}^2(y/L_{\mathrm{s}})},\end{array}$$(5)

where $n_0^{+}=V_0{B}_0/4\pi ec{L}_{\text{s}}$, e is the electron charge, and the electronic density profile

$$\begin{array}{c}\hfill n_{-}=-\frac{V_0{J}_{-}}{e{c}^2}=-n_0^{-}tanh\left(\frac{y}{L_{\mathrm{s}}}\right)cos\theta ,\end{array}$$(6)

where $n_0^{-}=V_0{B}_0/8ec{L}_y$. Although polarized, pulsar winds are most likely quasi-neutral, meaning that the plasma density greatly exceeds the charge density, n ≫ |ρ|/e. Thus, on top of these minimum densities, we add a uniform neutral density of pairs n₀, so that n_e = n₀ + n₋ for the total electron density and n_p = n₀ + n₊ for the total positron density. In the simulations, this fiducial density is set by the chosen upstream magnetization parameter,

$$\begin{array}{c}\hfill {\sigma }_0=\frac{B_0^2}{8\pi {\mathrm{\Gamma }}_0{n}_0{m}_{\mathrm{e}}{c}^2},\end{array}$$(7)

where m_e is the electron rest mass and ${\mathrm{\Gamma }}_0={(1-V_0^2/c^2)}^{-1/2}$ is the wind bulk Lorentz factor. Thus, the density contrast is

$$\begin{array}{cc}& \frac{n_0^{+}}{n_0}=\frac{2V_0}{c}\frac{{\sigma }_0{R}_0}{L_{\mathrm{s}}},\hfill \end{array}$$(8)

$$\begin{array}{cc}& \frac{n_0^{-}}{n_0}=\frac{\pi V_0}{c}\frac{{\sigma }_0{R}_0}{L_y},\hfill \end{array}$$(9)

where

$$\begin{array}{c}\hfill R_0=\frac{{\mathrm{\Gamma }}_0{m}_{\mathrm{e}}{c}^2}{e{B}_0}\end{array}$$(10)

is the fiducial particle Larmor radius.

2.4. Scale separation

The Crab Nebula spectrum suggests that the wind injects ∼1 TeV pairs at the shock front immersed in a B₀ ∼ 200 μG background field. The fiducial particle Larmor radius is then of order R₀ = 1 TeV/eB₀ ∼ 1.6 × 10¹³ cm. If one compares this gyroradius with the shock radius R_sh ∼ 0.1 pc, the scale separation is R_sh/R₀ ∼ 2 × 10⁴ such that we can verify that n_±/n₀ ≪ 1 as expected even for high magnetization σ₀ ≲ 100 (except perhaps in the equatorial plane for a nearly aligned pulsar).

The main numerical challenge is to reach a sufficiently large separation of scales between the microscopic Larmor radius and plasma skin-depth scales, where particle acceleration processes take place, and the global shock size scale. The most stringent constraint in PIC simulations is to resolve the plasma skin depth d_e and plasma frequency ω_pe scales and therefore these quantities determine the minimum spatial and time resolution of the simulations. In all runs, the fiducial skin-depth scale defined in the upstream flow is resolved by eight cells in all directions, where

$$\begin{array}{c}\hfill d_{\mathrm{e}}=\sqrt{\frac{{\mathrm{\Gamma }}_0{m}_{\mathrm{e}}{c}^2}{8\pi n_0{e}^2}}=\sqrt{{\sigma }_0}{R}_0.\end{array}$$(11)

The plasma density and the mean particle Lorentz factor in the downstream medium can differ significantly from the upstream parameters due to the compression of the flow and particle acceleration. A posteriori, we found that the plasma skin depth in the downstream flow is resolved by at least five cells in all the simulations.

The simulation time-step Δt is determined by the usual Courant-Friedrich-Lewy condition, such that ω_peΔt ≈ 8.75 × 10⁻². The largest simulation contains 65536 × 8192 cells along the x- and y-directions, which corresponds to a 8192d_e × 1024d_e box size. In this work, we do not perform a systematic study of the effect of the transverse size of the shock in order to focus our attention on the largest possible sizes. We inject 16 particles per cell per time-step. We ran simulations for σ₀ = 0, 0.1, 1, 10, 30, 100, which translates into physical box sizes ranging from 2590R₀ × 324R₀ for σ₀ = 0.1 to 44869R₀ × 5609R₀ for σ₀ = 30, which is close to the scale separation we are seeking for the Crab Nebula. Due to the strong anisotropy of the transverse magnetic field profile, the average magnetization in the wind is ${\overline{\sigma }}_0\approx 0.15{\sigma }_0$ for L_s/L_y = 0.5 (${\overline{\sigma }}_0\approx 0.4{\sigma }_0$ for L_s/L_y = 0.1 and ${\overline{\sigma }}_0\approx 0.065{\sigma }_0$ for L_s/L_y = 1). The cyclotron frequency ${\omega }_{\mathrm{c}}\mathrm{\Delta }t=8.75\times {10}^{-2}\sqrt{{\sigma }_0}$ is well resolved in all runs, even for σ₀ = 100. The largest simulation is integrated for about 7875ω_pet, or 43133ω_ct for σ₀ = 30. Scaled to the Crab Nebula, this represents a total simulation time of about 260 days which is of the order of a few times the dynamical timescale of the nebula, t_neb ∼ R_sh/c ∼ 100 days.

Aside from the separation of scale, we must also mitigate the effect of numerical Cherenkov radiation that tends to slow down and heat up ultra-relativistic beams (Greenwood et al. 2004). This instability grows with the amplitude of the Lorentz factor of the beam. The wind Lorentz factor in pulsar winds is uncertain but is most likely very high, Γ₀  ∼  10² − 10⁶. In this work, we scale down the wind velocity to V₀ = 0.99c or Γ₀  ≈  7, which fulfills the need to have Γ₀  ≫  1 and low numerical heating before the end of the simulation. A small temperature (kT₀/m_ec² = 10⁻²) is added and spatial filtering is applied to the current density to even further delay the onset of the instability. Radiative cooling (primarily synchrotron and inverse Compton) is neglected in this work.

2.5. Summary of all runs

All runs in this study were performed with the Cartesian version of the Zeltron PIC code (Cerutti et al. 2013; Cerutti & Werner 2019). Table 1 gives the list of all runs reported in this work.

3. Shock structure and dynamics

Figure 2 shows the shock structure at time ω_pet = 3920 in order of increasing magnetization, starting with a perfectly unmagnetized shock σ₀ = 0 (top panels) down to a strongly magnetized shock with σ₀ = 100 (bottom panel) and the transition in between these two extreme regimes. The unmagnetized case serves as a control simulation we can compare with previous studies (e.g., Spitkovsky 2008; Keshet et al. 2009; Sironi et al. 2013; Plotnikov et al. 2018). As expected for an unmagnetized shock, the interaction between the incoming and reflected beams leads to the formation of Weibel filaments which then mediate magnetic turbulence and ultimately the formation of the shock. These filaments are highly visible as self-generated magnetic structures of alternative polarities, whose strength peaks at the shock front and slowly decays downstream. The filamentation proceeds in the wind ahead of the shock front due to the reflected beam of particles propagating back upstream. This region, usually called the precursor, effectively decelerates the incoming flow and sustains magnetic turbulence (see e.g., Lemoine et al. 2019), which is the key ingredient to bring particles back to the shock and initiate the Fermi process. The downstream plasma density increases by n/n₀ ≈ 3 as expected from the MHD jump conditions of an unmagnetized perpendicular shock (Kennel & Coroniti 1984; Plotnikov et al. 2018).

Fig. 2. Total plasma density (left panels) and transverse magnetic field strength (right panels) for σ0 = 0, 0.1, 1, 10, 30, 100 (from top to bottom) and Ls/Ly = 0.5, at ωpet = 3920.

Adding a nonzero but subdominant field (σ₀ = 0.1) quenches the formation of Weibel filaments in most of the flow, except within the equatorial and polar regions where the effective magnetization is small (only a few filaments are visible in each of these regions). The upstream magnetic field is compressed downstream by a factor of approximately three as expected in the low-σ limit. Similarly, the plasma density is compressed by approximately the same amount, with a noticeable depletion in the equatorial and polar regions. We observe a well-defined high-density rim at the shock front at intermediate heights. This is a known feature of magnetized shocks which results from the magnetic reflection of the incoming flow of particles (see, e.g., Plotnikov et al. 2018). When the particles cross the shock front, they lose the support of the transverse electric field which allowed them to go in straight lines upstream. As a result, the incoming particles gyrate coherently at the shock front resulting in this characteristic plasma density bump and the emission of an electromagnetic precursor upstream (Gallant et al. 1992).

At even higher magnetization (σ₀ ≳ 1), the shock structure changes dramatically and we are now leaving known territory for a new phenomenology. As the magnetization increases, what appears to be a shock front travels faster upstream. In magnetized regions located at intermediate heights, the jumps in the plasma density and the magnetic field strength across the shock front also decrease with increasing magnetization. While these features are consistent with a weak shock, the equatorial and polar regions behave very differently. The flow is strongly compressed into a highly turbulent state driven by kinks (current-driven) and Kelvin-Helmholtz vortices (shear-driven). The compression of the flow into low-field regions is the result of the magnetic pressure force. The global bulk flow then quickly converges towards the pattern shown in the upper panel in Fig. 3 for σ₀ = 30. This figure shows the dimensionless bulk momentum of the flow along the x-direction, U_x = ΓV_x/c, and plasma velocity streamlines. In magnetized regions (y/d_e = ±256), the flow decelerates down to about ∼ + 0.3c. In the low-field regions, the flow velocity returns back to the shock with a net bulk velocity ≈ − 0.5c in the downstream medium. This is a major difference with uniform shocks where the plasma is at rest in the downstream region. In the transition region (y/d_e ∼ ±100), there is a strong velocity shear ΔV ∼ 0.8c (Fig. 4) which drives the formation of Kelvin-Helmholtz vortices well visible in the density maps as well as in the streamline pattern.

Fig. 3. Top: plasma bulk momentum along the x-direction, Ux = ΓVx/c, for σ0 = 30 over the full simulation box at time ωpet = 3920. Bottom: zoomed-in view of the equatorial region of the shock front shown by the rectangular box in the upper panel. Solid black lines with arrows are the plasma velocity streamlines.

Fig. 4. Plasma bulk velocity profile in the downstream flow along the x- and y-directions at x/de = 1500 for σ0 = 30 at time ωpet = 3920.

An intriguing and robust feature of the shock is the spearhead-shaped structures developing at the shock front in the equatorial and polar regions, which are elongated along the direction of the flow. These structures are low-field, low-density regions characterized by a mildly relativistic backflow motion up to U_x ≈ −4 at x/d_e ≈ 700 for σ₀ = 30 (see zoomed-in view in the bottom panel in Fig. 3). They are also characterized by a large and abrupt velocity shear at their boundaries. The sheath-like structures around them gradually deflect the incoming flow sideways, such that there is no clear sign of a standard shock pattern here. Away from this triangular-shaped precursor drilling through the upstream medium, the incoming flow is perfectly laminar with no sign of plasma turbulence. The size of these cavities continuously grows with time without any sign of saturation. The kink-like motion of the plasma concentrated in the midplane seems to depart from the base of these structures. The plasma carries away the electric current within high-density filaments, simply referred to as “current layers” in the following. The current flows along the equator to sustain the jump in the magnetic field polarity. In the early phases of the simulation (ω_pet ≲ 4000), the current then flows along the ±y-directions at the x = L_x boundary and reaches the polar regions where it flows in the opposite direction. This electric circuit gradually closes through the shock front and the spearhead cavities. At later stages (ω_pet ≳ 4000), Kelvin-Helmholtz vortices combined with the kink lead to an efficient mixing of the downstream flow into a highly turbulent state. The top panel in Fig. 5 shows the total current and its schematic path within the numerical box at a late evolutionary stage (ω_pet = 7840), when a turbulent mixing state has been reached far downstream. It is important to notice that the downstream flow, and in particular the current layers and their associated cavities, are electrically charged with a net negative charge in the equator and a net positive charge at the poles (and vice versa if Ω ⋅ B <  0). Before they mix and reconnect far downstream, each layer is surrounded by a low-density background plasma with the opposite sign of charge, but of the same sign as the upstream flow (Fig. 5, bottom panel).

Fig. 5. Top: magnitude of the total electric current density, |J|, normalized by the fiducial upstream current J0 = cB0/4πLy for σ0 = 30 at ωpet = 7840. The white arrows schematically show how the current flows and closes within the simulation box. Bottom: plasma charge density, ρ, normalized by the fiducial upstream charge density J0/c.

4. Particle acceleration

4.1. Total spectra and maximum energy

This rather complex shock structure leads to efficient nonthermal particle acceleration. The upper panel in Fig. 6 shows the total particle spectrum, udN/du, where u = γβ is the dimensionless particle momentum at time ω_pet = 3920 for all the magnetizations simulated here. As expected, the unmagnetized shock produces a high-energy power-law tail extending beyond the thermal bath. For σ₀ = 0.1, particle acceleration is quenched in most of the shock front (thermal spectrum), except in the low-field regions which produces a weak excess starting below udN/du ≲ 10⁻⁴ and extending to the same maximum energy as in the unmagnetized shock. In contrast, strongly magnetized shocks (σ₀ >  1) present a pronounced high-energy power-law tail extending to energies far beyond unmagnetized shocks, with a maximum Lorentz factor γ_max/Γ₀ ∼ 500 for σ₀ = 30 compared with γ_max/Γ₀ ∼ 20 for σ₀ = 0. The power law hardens as well with increasing magnetization, approaching the canonical first-order Fermi acceleration spectrum dN/du ∝ u^−2.2 for σ₀ = 30. At this magnetization, the spectrum breaks and steepens at high energy (γ ≳ 300) and a new component emerges whose nature will become clear later. The bottom panel in Fig. 6 shows the dependence of the particle spectrum with the striped wind filling factor, L_s/L_y. Particle acceleration is more pronounced at low L_s, which is consistent with the trend reported earlier that particle acceleration is more effective when the wind is more magnetized on average; that is, ${\overline{\sigma }}_0\approx 4.5$ for L_s/L_y = 0.1, as compared with ${\overline{\sigma }}_0\approx 0.65$ for L_s/L_y = 1.

Fig. 6. Top: total particle spectrum as a function of the dimensionless total particle momentum, u = γβ for all the magnetizations explored in this study and at the same time ωpet = 3920. Bottom: effect of the striped wind filling factor on the final particle spectrum, Ls/Ly for σ0 = 10 at the same time as in the top panel. The black dotted line shows a power law dN/du ∝ u−p of index 2.2. The black dashed line is a 2D relativistic Maxwellian distribution for comparison.

Figure 7 shows the time evolution of the maximum particle Lorentz factor of the total spectrum, γ_max(t). For the unmagnetized shock (σ₀ = 0), the maximum energy approximately grows as ${\gamma }_{\max }/{\mathrm{\Gamma }}_0\approx 0.5\sqrt{{\omega }_{\mathrm{pe}}t}$ without any sign of saturation, in agreement with Sironi et al. (2013). The square-root dependence reflects the microscopic nature of the Weibel-driven turbulence. Finite but mildly magnetized solutions (σ₀ = 0.1, 1, 10) show a similar acceleration rate at the early stages, but this is followed by a saturation at ω_pet ≳ 500 which is particularly visible for the σ₀ = 0.1 solution. This saturation is related to the finite size of the turbulent region in the upstream flow, which approximately scales as the particle Larmor radius in the background field (Lemoine & Pelletier 2010; Sironi et al. 2013; Plotnikov et al. 2018).

Fig. 7. Time evolution of the maximum particle Lorentz factor of the total spectrum, γmax, in all the runs. The dashed black lines show ${\gamma }_{\max }=0.5{\mathrm{\Gamma }}_0\sqrt{{\omega }_{\mathrm{pe}}t}$ and γmax = 0.5Γ0ωpet.

For uniform shocks, there would be no more evolution of the particle spectrum. In contrast, we observe here that the maximum energy then increases again and at a much faster rate in the late evolution of the simulation. For σ₀ = 1, γ_max increases again at ω_pet ≳ 2 × 10³. For σ₀ = 10, it occurs earlier, at ω_pet ≳ 10³, followed by a quasi-linear evolution γ_max ∝ t. More highly magnetized solutions (σ₀ = 30, 100) only show a linear evolution of the maximum particle energy with time from the beginning of the simulations. The longest runs (σ₀ = 10 and 30) show no sign of saturation. The linear increase of the particle energy with time is evidence for efficient particle acceleration, compatible with the Bohm regime but in sharp contrast with shock acceleration mediated by self-generated microturbulence (Sironi et al. 2013).

4.2. Phase space and local spectra

To gain further physical insight into the origin and location of particle acceleration, we compute the mean particle Lorentz factor, ⟨γ⟩, in each cell of the simulation as reported in the top panel in Fig. 8 at ω_pet = 7840. High-energy particles are located in high-density regions within the current layers and inside the spearhead-shaped cavities at the shock front. A look into the x–p_x phase-space distribution within the midplane reveals that the latter contain the highest energy particles in the simulations, γ ∼ 10⁴, in the form of a strongly beamed distribution traveling against the incoming flow. Electrons (positrons) are preferentially accelerated in the equatorial (polar) cavities and vice versa if Ω ⋅ B <  0. The asymmetry between both species as well as the anisotropy of their momentum distribution decrease downstream and even fully disappear where the flow becomes turbulent (x/d_e ≳ 4000). The bottom panels in Fig. 8 show the particle spectra measured in four areas defined in the top panel. Areas 1 and 2 are restricted to the shock-front cavities. The asymmetric acceleration between electrons and positrons is clearly visible here. In contrast to the total spectrum, the spectrum measured in these cavities is consistent with a single power law extending from γ_min = Γ₀ to γ_max ≈ 10⁴ with an index close to but slightly harder than −2.2. Area 3 is limited to the shock front at intermediate latitudes where particle acceleration is quenched. Area 4 focuses on the far-downstream region where the current layers reconnect and merge in a turbulent manner. The spectrum is hard at low energies, but cuts off noticeably below the maximum energy measured in the cavities. This difference explains the high-energy break at γ ≳ 300 in the total spectrum, beyond which the spectral component from the shock-front cavities takes over.

Fig. 8. Top panel: mean particle Lorentz factor at ωpet = 7840 for σ0 = 30. Middle panels: electron and positron x–px phase-space distribution dN/dxdydpx along the midplane integrated within y/de ± 75. Bottom panels: local particle spectra measured at four locations shown in the top panel and labeled from ① to ④. The dashed line is the total spectrum, and the dotted line is dN/du ∝ u−2.2 for reference.

4.3. Particle trajectories

Figure 9 shows two typical high-energy particle trajectories out of a randomly selected sample of 2000 particles for σ₀ = 30, which are meant to illustrate the different acceleration processes at work. Particle 1 shown in the left panels represents the acceleration history of the highest energy particles found in the simulations. As already pointed out in Sect. 4.2, these particles inflate the elongated cavities near the shock front and in the upstream medium. In the early phases (ω_pet ≲ 2500), the particle is trapped inside the cavity with little acceleration. At ω_pet ∼ 2500, the shock front catches up with the particle, which leads to a rapid and uninterrupted acceleration up to at least γ ≳ 2000. During this phase, the particle is trapped in a region that looks like a wake produced behind the cavities where the current layer forms and departs from. The particle trajectory moves back and forth across the equatorial plane where the magnetic field reverses, such that it is well described by the relativistic analog of Speiser orbits (Speiser 1965; Cerutti et al. 2012). As the particle surfs on the wake, its Larmor radius becomes sufficiently large to experience the strong macroscopic bulk-velocity shear, and therefore we associate the acceleration of these particles with the tangential shear-flow acceleration mechanism, which in essence is another form of the Fermi process. In this regime, the energy gain is due to Lorentz-frame transformation as the particle is scattered back and forth across the velocity-shear layer, and is of order Δγ/γ ∼ Γ_s − 1 after each crossing, where ${\mathrm{\Gamma }}_{\mathrm{s}}=1/\sqrt{1-\mathrm{\Delta }{V}^2/c^2}$, and ΔV = (V₁ − V₂)/(1 − V₁V₂/c²) is the velocity shear sampled by the particle between frames 1 and 2 (Ostrowski 1990; Rieger & Duffy 2004). In this simulation, the velocity shear is mildly relativistic, ΔV/c ∼ 0.8, leading to Δγ/γ ∼ 0.7. The acceleration proceeds until the particle is kicked out and advected downstream, which plays the role of an escape mechanism.

Fig. 9. Typical high-energy particle trajectories accelerated at the shock front (left panels) and in the turbulent downstream medium (right panels) for σ0 = 30. Top panels: the Lorentz factor is color-coded along the particle trajectories, itself plotted on top of the plasma density map (gray scale) at the final time of the particle tracking. Bottom panels: time evolution of the particle Lorentz factor.

Particle 2 is more representative of particles accelerated in the turbulent flow further downstream, and therefore of the bulk of the energetic particles in the simulations. It is injected at intermediate latitudes and flows in the laminar magnetized medium between the equatorial and polar current layers without significant energy gain. At ω_pet ∼ 5000, the particle is captured by a current layer where it experiences an abrupt acceleration, from γ ∼ Γ₀ to γ ∼ 200. We associate this impulsive phenomenon as direct acceleration via relativistic reconnection occurring within the current sheet, which naturally boosts the particle energy to γ ∼ Γ₀σ₀ = 200 (e.g., Werner et al. 2016). This event is followed by a much slower stochastic acceleration. At this stage, the particle has reached the turbulent flow where current layers are mixed together and collide at nearly random velocities. This environment favors multiple particle scattering which leads to a stochastic increase or decrease of its energy, but with a net positive gain. In this sense, this process is reminiscent of a second-order Fermi acceleration.

5. Synchrotron radiation

Figure 10 shows the instantaneous synchrotron spectrum emitted by the pairs in the σ₀ = 30 simulation at ω_pet = 7840, assuming the plasma is optically thin everywhere. To reconstruct the spectrum, we use a delta-function approximation. Each particle emits a single photon radiating away the total power lost by the parent particle, P_sync, with a frequency set by the synchrotron critical frequency. We recall here that radiative losses are neglected in the simulations for simplicity, and therefore the synchrotron spectrum computed here is not meant to be compared with, for instance, the observed Crab Nebula spectrum, which most likely results from a cooled particle distribution. Instead, our goal here is to characterize the main emission pattern at the shock in different frequency bands. The power spectrum of a single photon is then given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}t\mathrm{d}\nu }\approx P_{\mathrm{sync}}\delta (\nu -{\nu }_{\mathrm{c}}),\end{array}$$(12)

Fig. 10. Total synchrotron spectrum for σ0 = 30 at ωpet = 7840 (blue solid line). The dashed lines show the spectra emitted near the shock front (red) and in the far-downstream region (black). The dotted line is a pure power law which would be emitted by a p = 2.2 power-law electron spectrum, such that dE/dtdν ∝ ν( − p + 1)/2 = ν−0.6. The frequency bands labeled ⓐ, ⓑ, and ⓒ refer to Fig. 11.

where

$$\begin{array}{c}\hfill P_{\mathrm{sync}}=\frac{2}{3}{r}_{\mathrm{e}}^{2}c{\gamma }^2{\stackrel{\sim }{B}}_{\perp }^2,\end{array}$$(13)

r_e = e²/m_ec² is the classical radius of the electron, and

$$\begin{array}{c}\hfill {\nu }_{\mathrm{c}}=\frac{3e}{4\pi m_{\mathrm{e}}c}{\gamma }^2{\stackrel{\sim }{B}}_{\perp },\end{array}$$(14)

is the critical synchrotron frequency. ${\stackrel{\sim }{B}}_{\perp }$ is the effective perpendicular (to the particle velocity vector) magnetic field in the presence of a strong electric field given by (e.g., Cerutti et al. 2016)

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathit{B}}}_{\mathbf{\perp }}=\mathit{E}+\mathit{\beta }\times \mathit{B}-(\mathit{\beta }·\mathit{E})\mathit{\beta }.\end{array}$$(15)

Frequencies are normalized by the fiducial synchrotron critical frequency

$$\begin{array}{c}\hfill {\nu }_0=\frac{3e}{4\pi m_{\mathrm{e}}c}{\mathrm{\Gamma }}_0^2{B}_0.\end{array}$$(16)

The total synchrotron spectrum presents three main features: (i) a low-energy bump centered around ν/ν₀ ∼ 1, (ii) a plateau at intermediate frequencies (dE/dtdν ∝ ν⁰, 10¹ ≲ ν/ν₀ ≲ 10⁴) followed by, (iii) a steep power-law decline cutting off at ν/ν₀ ∼ 10⁸. The latter is slightly steeper than the canonical dE/dtdν ∝ ν^{( − p + 1)/2} = ν^−0.6 synchrotron spectrum emitted by a p = 2.2 power-law electron spectrum immersed into a uniform magnetic field.

The bulk of the synchrotron spectrum originates from the far-downstream region where the flow isotropizes and turbulent reconnection accelerates the particles (see black dashed line in Fig. 10). The spectrum emitted in the vicinity of the shock front is composed of the low-energy bump and a hard power-law tail (dE/dtdν ∝ ν^−0.3) spanning over nearly seven orders of magnitude in frequency range. The contribution from the shock front is subdominant at almost all frequencies except at the high-end of the spectrum, ν/ν₀ ≳ 10⁶ where it rises above the steeper spectrum of the far-downstream region. Figure 11 shows the spatial distribution of the emitted synchrotron flux integrated in three frequency bands: (i) low-, ν/ν₀ <  10², (ii) intermediate-, 10² <  ν/ν₀ <  10⁶, and (iii) high-frequency band ν/ν₀ >  10⁶ labeled as ⓐ, ⓑ and ⓒ respectively. The low-energy flux is rather uniformly distributed in the downstream flow, that is, from the shock front to the back end of the numerical box, from the most magnetized regions to the current layers, but with the notable exception of the shock front cavities. At intermediate frequencies in band ⓑ, the emission in concentrated within the outer edges of the current layers up to the shock front cavities which are lighting up in this band. Although high-energy particles are also located inside the layers (see Fig. 8), they do not emit significant radiation because the fields almost vanish in there. Away from the current layers, in smooth magnetized regions, the field is high but energetic particles are absent leading to no flux from these regions. In the high-energy band (ⓒ), the emission is dominated by the edges of the shock-front cavities piercing through the upstream, as well as by the base of the current sheets in the near downstream medium. The rest of the upstream flow remains dark in all bands as expected because ${\stackrel{\sim }{B}}_{\perp }\approx 0$.

Fig. 11. Spatial distribution of the synchrotron radiation flux integrated in the frequency bands defined in Fig. 10: ⓐ (ν/ν0 <  102), ⓑ (102 <  ν/ν0 <  106), and ⓒ (ν/ν0 >  106).

6. Discussion and conclusion

Here, we show that the anisotropic nature of the pulsar wind has a dramatic effect on the structure and evolution of the shock front. The usual local plane-parallel approximation does not apply because of the critical role of the global dynamics of the downstream flow, and therefore a latitudinally broad simulation box is essential. A salient feature is the formation of a sharp velocity shear between strongly- and weakly magnetized regions, which combined with current-driven instabilities leads to strong plasma turbulence in the downstream flow. Current sheets forming in the equatorial plane and at the poles mix and reconnect leading to efficient nonthermal particle acceleration. The efficiency of particle acceleration increases with σ, which is similar to relativistic reconnection but opposite to Fermi acceleration in uniform shocks. Turbulent reconnection leads to a power-law electron spectrum with a slope that hardens with σ, as expected from recent studies of relativistic reconnection (Sironi & Spitkovsky 2014; Werner et al. 2016) and kinetic turbulence (Zhdankin et al. 2017, 2018). Scaled up to the Crab Nebula, with dN/dγ ∝ γ^−2.2, the injected (uncooled) X-ray electron spectrum is consistent with the high-σ shock solution σ₀ = 30, or a mean upstream magnetization ${\overline{\sigma }}_0\approx 5$. This result is compatible with global 3D MHD models which advocate high-σ solutions to explain the morphology of the Crab Nebula (Porth et al. 2013, 2014). The spectrum extends from Γ₀ to ∼Γ₀σ₀ beyond which it steepens significantly and where another component takes over and extends the total spectrum to even higher energies. It is tempting to connect this result with the mysterious high-energy break in the Crab Nebula spectrum, where the electron spectral index decreases by ∼0.5 (Meyer et al. 2010).

The high-energy component originates from another robust and extraordinary feature of anisotropic magnetized shocks, which is the formation of elongated cavities at the base of the polar and equatorial regions drilling through the upstream medium. These structures are low-field, low-density regions moving with mildly relativistic speeds against the incoming flow, and their sizes continuously grow with time. They are inflated by the highest energy particles in the box, which follow relativistic Speiser orbits. These special trajectories, typically found in reconnection layers, are captured by the midplane where the magnetic field polarity reverses. This trapping mechanism provides stability to the cavities themselves, whose sizes constantly adjust to the particle Larmor radius. These particles are energized via shear-flow acceleration at the interface between the cavities and the incoming flow. This component alone explains the highest-energy part of the particle spectrum. The maximum particle energy increases nearly linearly with time, γ_max ∝ t, in contrast to Weibel-dominated shock acceleration where ${\gamma }_{\max }\propto \sqrt{t}$ (Sironi et al. 2013; Plotnikov et al. 2018), meaning that the acceleration process is very efficient, close to the Bohm limit. Another important difference with Weibel-dominated shock acceleration is that the particle spectrum does not show signs of saturation, the maximum energy grows steadily until the end of the simulations. Presumably, the particle energy will be limited by the transverse size of the shock or by radiative losses such as those in the Crab Nebula where the electron maximum energy is limited by the synchrotron burn-off limit (Guilbert et al. 1983; de Jager et al. 1996). These cavities also have the peculiarity to preferentially accelerate one sign of charge, electrons in the equatorial region, positrons (ions) at the poles, and vice versa if Ω ⋅ B <  0, as also reported by Giacinti & Kirk (2018).

The modeling of synchrotron radiation indicates that while the bulk of the emission is produced quasi-isotropically in the downstream region, the high-end of the synchrotron spectrum is concentrated within the edges and the wakes of the shock front cavities where both the field and the particle energy are the highest. Translated in terms of the Crab Nebula features, we predict that the 100 MeV emission is preferably localized at the inner ring, and specifically on the side receding away from us as the emission in the cavities is Doppler boosted in the direction opposite to the incoming pulsar wind. There might also be a weaker contribution from the base of the counter jet because of the smaller volume involved in comparison with the equatorial region. The wake behind the cavities is highly dynamical. Due to the kink instability, the high-energy beam sweeps a wide angular range. Therefore, gamma-ray flares at the high end of the synchrotron spectrum come out as a natural consequence of particle acceleration at the pulsar wind termination shock, and more generally in any relativistic, magnetized, and anisotropic shocks with a possible application to the hotspots of relativistic jets and gamma-ray bursts. The mildly relativistic bulk motion of the backflow, with Γ ∼ 3 − 4, would naturally push the radiation energy above the rest-frame 160 MeV synchrotron burn-off limit, a persistent feature of the Crab Nebula gamma-ray flares (Tavani et al. 2011; Abdo et al. 2011). Although these results suggest a connection between particle acceleration in the cavities and gamma-ray flares, more work is needed for a solid conclusion.

An obvious limitation of the proposed model is its Cartesian plane-parallel geometry. Any curvature of the shock front and the radial expansion of the downstream flow are therefore neglected in this work. A more realistic configuration with a spherical geometry and a finite curvature of the shock front would be the logical next step to confirm our findings. This would also break the symmetry between the poles and the equator, which play nearly identical roles in this study. In particular, the cavity in the equator may then have a dominant contribution in the acceleration (charge asymmetry) and radiation pattern (gamma-ray flares) because of its larger volume in comparison with the poles. An extension to 3D simulations would also be desirable to fully capture magnetic reconnection that occurs here in the out-of-plane direction and therefore some of the important features of reconnection are missing here (e.g., tearing mode and plasmoid formation). Three-dimensional simulations can also capture possible departures from axisymmetry, which may explain for instance the knotty nature of the inner ring in the Crab Nebula. Synchrotron cooling could play an important role in the dynamics of the shock front cavities where it is most severe. Although they are probably highly subdominant in number compared with pairs, ions are most likely present in the wind and in the nebula. If the efficient particle acceleration mechanism revealed in this work also applies to them, pulsar wind nebulae could be an important source of the Galactic cosmic-ray population. Dissipation of the oscillating current sheet in the pulsar wind may lead to efficient particle acceleration ahead of the termination shock (Cerutti & Philippov 2017). An excess of energetic particles in the equatorial regions may affect the shock dynamics in return. The exploration of all of the above effects provides a wide array of possible future investigations.

Movie

Simulation movies

all trajectories sigma30 Access here

Bz sigma30 Access here

densities sigma30 Access here

Acknowledgments

We are grateful to the anonymous referee for his/her very supportive report. This work has been supported by the Programme National des Hautes Énergies of CNRS/INSU and CNES. We acknowledge PRACE and GENCI (allocations A0050407669 and A0070407669) for awarding us access to Joliot-Curie at GENCI@CEA, France.

References

All Tables

All Figures

	Fig. 1. Top: transverse magnetic field profile in the pulsar wind along the y-direction, Bz(y), according to Eq. (1) for Ls/Ly = 0.1, 0.5, 1. Bottom: corresponding electric current profile, Jx(y).
In the text

	Fig. 2. Total plasma density (left panels) and transverse magnetic field strength (right panels) for σ0 = 0, 0.1, 1, 10, 30, 100 (from top to bottom) and Ls/Ly = 0.5, at ωpet = 3920.
In the text

	Fig. 3. Top: plasma bulk momentum along the x-direction, Ux = ΓVx/c, for σ0 = 30 over the full simulation box at time ωpet = 3920. Bottom: zoomed-in view of the equatorial region of the shock front shown by the rectangular box in the upper panel. Solid black lines with arrows are the plasma velocity streamlines.
In the text

	Fig. 4. Plasma bulk velocity profile in the downstream flow along the x- and y-directions at x/de = 1500 for σ0 = 30 at time ωpet = 3920.
In the text

	Fig. 5. Top: magnitude of the total electric current density, |J|, normalized by the fiducial upstream current J0 = cB0/4πLy for σ0 = 30 at ωpet = 7840. The white arrows schematically show how the current flows and closes within the simulation box. Bottom: plasma charge density, ρ, normalized by the fiducial upstream charge density J0/c.
In the text

Fig. 6. Top: total particle spectrum as a function of the dimensionless total particle momentum, u = γβ for all the magnetizations explored in this study and at the same time ωpet = 3920. Bottom: effect of the striped wind filling factor on the final particle spectrum, Ls/Ly for σ0 = 10 at the same time as in the top panel. The black dotted line shows a power law dN/du ∝ u−p of index 2.2. The black dashed line is a 2D relativistic Maxwellian distribution for comparison.

In the text

	Fig. 7. Time evolution of the maximum particle Lorentz factor of the total spectrum, γmax, in all the runs. The dashed black lines show ${\gamma }_{\max }=0.5{\mathrm{\Gamma }}_0\sqrt{{\omega }_{\mathrm{pe}}t}$ and γmax = 0.5Γ0ωpet.
In the text

Fig. 8. Top panel: mean particle Lorentz factor at ωpet = 7840 for σ0 = 30. Middle panels: electron and positron x–px phase-space distribution dN/dxdydpx along the midplane integrated within y/de ± 75. Bottom panels: local particle spectra measured at four locations shown in the top panel and labeled from ① to ④. The dashed line is the total spectrum, and the dotted line is dN/du ∝ u−2.2 for reference.

In the text

Fig. 9. Typical high-energy particle trajectories accelerated at the shock front (left panels) and in the turbulent downstream medium (right panels) for σ0 = 30. Top panels: the Lorentz factor is color-coded along the particle trajectories, itself plotted on top of the plasma density map (gray scale) at the final time of the particle tracking. Bottom panels: time evolution of the particle Lorentz factor.

In the text

	Fig. 10. Total synchrotron spectrum for σ0 = 30 at ωpet = 7840 (blue solid line). The dashed lines show the spectra emitted near the shock front (red) and in the far-downstream region (black). The dotted line is a pure power law which would be emitted by a p = 2.2 power-law electron spectrum, such that dE/dtdν ∝ ν( − p + 1)/2 = ν−0.6. The frequency bands labeled ⓐ, ⓑ, and ⓒ refer to Fig. 11.
In the text

	Fig. 11. Spatial distribution of the synchrotron radiation flux integrated in the frequency bands defined in Fig. 10: ⓐ (ν/ν0 <  102), ⓑ (102 <  ν/ν0 <  106), and ⓒ (ν/ν0 >  106).
In the text