A global model of particle acceleration at pulsar wind termination shocks

Pulsar wind nebulae are efficient particle accelerators, and yet the processes at work remain elusive. Self-generated, microturbulence is too weak in relativistic magnetized shocks to accelerate particles over a wide energy range, suggesting that the global dynamics of the nebula may be involved in the acceleration process instead. In this work, we study the role played by the large-scale anisotropy of the transverse magnetic field profile on the shock dynamics. We performed large two-dimensional particle-in-cell simulations for a wide range of upstream plasma magnetizations. A large-scale velocity shear and current sheets form in the equatorial regions and at the poles, where they drive strong plasma turbulence via Kelvin-Helmholtz vortices and kinks. The mixing of current sheets in the downstream flow leads to efficient nonthermal particle acceleration. The power-law spectrum hardens with increasing magnetization, akin to those found in relativistic reconnection and kinetic turbulence studies. The high end of the spectrum is composed of particles surfing on the wake produced by elongated spearhead-shaped cavities forming at the shock front and piercing through the upstream flow. These particles are efficiently accelerated via the shear-flow acceleration mechanism near the Bohm limit. Magnetized relativistic shocks are very efficient particle accelerators. Capturing the global dynamics of the downstream flow is crucial to understanding them, and therefore local plane parallel studies may not be appropriate for pulsar wind nebulae and possibly other astrophysical relativistic magnetized shocks. A natural outcome of such shocks is a variable and Doppler-boosted synchrotron emission at the high end of the spectrum originating from the shock-front cavities, reminiscent of the mysterious Crab Nebula gamma-ray flares.


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 9 , with a major spectral break about half way, γ ∼ 10 5 . 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 −2 , 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).
Article number, page 1 of 12 arXiv:2008.07253v1 [astro-ph.HE] 17 Aug 2020 A&A proofs: manuscript no. shock_pwn 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 largescale 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 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 Sections 3-5, which are further discussed in Section 6, with particular emphasis on the Crab Nebula.

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 ultrarelativistic 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 counterstreaming 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.

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 1 where the z-direction plays the role of the toroidal direction, B 0 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 o ), L s /L y = 0.5 (χ = 45 o ), and for an orthogonal rotator L s /L y = 1 (χ = 90 o ). The background electric field is the ideal advection field 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 0 , and so E y (y) = V 0 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.

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 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 Figure 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 so that J x ≈ ρc in the ultrarelativistic limit.

Plasma density
Assuming that both species move along the +x-direction at the same speed V 0 , 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 where n + 0 = V 0 B 0 /4πecL s , e is the electron charge, and the electronic density profile where n − 0 = V 0 B 0 /8ecL 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 0 , so that n e = n 0 + n − for the total electron density and n p = n 0 + n + for the total positron density. In the simulations, this fiducial density is set by the chosen upstream magnetization parameter, where m e is the electron rest mass and Γ 0 = (1 − V 2 0 /c 2 ) −1/2 is the wind bulk Lorentz factor. Thus, the density contrast is where R 0 = Γ 0 m e c 2 eB 0 (10) is the fiducial particle Larmor radius.

Scale separation
The Crab Nebula spectrum suggests that the wind injects ∼1 TeV pairs at the shock front immersed in a B 0 ∼ 200µG background field. The fiducial particle Larmor radius is then of order R 0 = 1 TeV/eB 0 ∼ 1.6×10 13 cm. If one compares this gyroradius with the shock radius R sh ∼ 0.1pc, the scale separation is R sh /R 0 ∼ 2 × 10 4 such that we can verify that n ± /n 0 1 as expected even for high magnetization σ 0 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 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 −2 . 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, 0.1, 1, 10, 30, 100, which translates into physical box sizes ranging from 2590R 0 × 324R 0 for σ 0 = 0.1 to 44869R 0 × 5609R 0 for σ 0 = 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σ 0 ≈ 0.15σ 0 for L s /L y = 0.5 (σ 0 ≈ 0.4σ 0 for L s /L y = 0.1 andσ 0 ≈ 0.065σ 0 for L s /L y = 1). The cyclotron frequency ω c ∆t = 8.75 × 10 −2 √ σ 0 is well resolved in all runs, even for σ 0 = 100. The largest simulation is integrated for about Article number, page 3 of 12 A&A proofs: manuscript no. shock_pwn Table 1. List of all PIC simulations reported in this work.

Run
Size 4096 × 1024 0.1 0.5 0.015 S1 4096 × 1024 1 0.5 0.15 S10 8192 × 1024 10 0.5 1.5 S10_LS01 4096 × 1024 10 0.1 4 S10_LS1 4096 × 1024 10 1 0.65 S30 8192 × 1024 30 0.5 4.5 S100 4096 × 1024 100 0.5 15 7875ω pe t, or 43133ω c t for σ 0 = 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, 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, Γ 0 ∼ 10 2 −10 6 . In this work, we scale down the wind velocity to V 0 = 0.99c or Γ 0 ≈ 7, which fulfills the need to have Γ 0 1 and low numerical heating before the end of the simulation. A small temperature (kT 0 /m e c 2 = 10 −2 ) 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.

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. Figure 2 shows the shock structure at time ω pe t = 3920 in order of increasing magnetization, starting with a perfectly unmagnetized shock σ 0 = 0 (top panels) down to a strongly magnetized shock with σ 0 = 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 0 ≈ 3 as expected from the MHD jump conditions of an unmagnetized perpendicular shock (Kennel & Coroniti 1984;Plotnikov et al. 2018).

Shock structure and dynamics
Adding a nonzero but subdominant field (σ 0 = 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 highdensity 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 .
At even higher magnetization (σ 0 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 Figure 3 for σ 0 = 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 (Figure 4) which drives the formation of Kelvin-Helmholtz vortices well visible in the density maps as well as in the streamline pattern.
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 σ 0 = 30 (see zoomed-in view in the bottom panel in Figure 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 (ω pe t 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 Article number, page 4 of 12 Benoît Cerutti and Gwenael Giacinti: A global model of particle acceleration at pulsar wind termination shocks closes through the shock front and the spearhead cavities. At later stages (ω pe t 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 Figure 5 shows the total current and its schematic path within the numerical box at a late evolutionary stage (ω pe t = 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 ( Figure 5, bottom panel).

Total spectra and maximum energy
This rather complex shock structure leads to efficient nonthermal particle acceleration. The upper panel in Figure 6 shows the total particle spectrum, udN/du, where u = γβ is the di-mensionless particle momentum at time ω pe t = 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 = 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 −4 and extending to the same maximum energy as in the unmagnetized shock. In contrast, strongly magnetized shocks (σ 0 > 1) present a pronounced high-energy power-law tail extending to energies far beyond unmagnetized shocks, with a maximum Lorentz factor γ max /Γ 0 ∼ 500 for σ 0 = 30 compared with γ max /Γ 0 ∼ 20 for σ 0 = 0. The power law hardens as well with increasing magnetization, approaching the canonical firstorder Fermi acceleration spectrum dN/du ∝ u −2.2 for σ 0 = 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 Figure 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 magne-Article number, page 5 of 12 A&A proofs: manuscript no. shock_pwn  tized on average; that is,σ 0 ≈ 4.5 for L s /L y = 0.1, as compared withσ 0 ≈ 0.65 for L s /L y = 1. Figure 7 shows the time evolution of the maximum particle Lorentz factor of the total spectrum, γ max (t). For the unmagnetized shock (σ 0 = 0), the maximum energy approximately grows as γ max /Γ 0 ≈ 0.5 √ ω 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 = 0.1, 1, 10) show a similar acceleration rate at the early stages, but this is followed by a saturation at ω pe t 500 which is particularly visible for the σ 0 = 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).
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 σ 0 = 1, γ max increases again at ω pe t 2 × 10 3 . For σ 0 = 10, it occurs earlier, at ω pe t 10 3 , followed by a quasi-linear evolution γ max ∝ t. More highly magnetized solutions (σ 0 = 30, 100) only show a linear evolution of the maximum particle energy with time from the beginning of the simulations. The longest runs (σ 0 = 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).

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 Figure 8 at ω pe t = 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, Article number, page 6 of 12 Benoît Cerutti and Gwenael Giacinti: A global model of particle acceleration at pulsar wind termination shocks 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 Figure 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 = Γ 0 to γ max ≈ 10 4 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. Figure 9 shows two typical high-energy particle trajectories out of a randomly selected sample of 2000 particles for σ 0 = 30, which are meant to illustrate the different acceleration processes at work. Particle 1 shown in the left panels represents the ac-celeration 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 (ω pe t 2500), the particle is trapped inside the cavity with little acceleration. At ω pe t ∼ 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 Article number, page 8 of 12

Particle trajectories
Benoît Cerutti and Gwenael Giacinti: A global model of particle acceleration at pulsar wind termination shocks 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. In the 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.
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 Γ s = 1/ 1 − ∆V 2 /c 2 , and ∆V = (V 1 − V 2 )/(1 − V 1 V 2 /c 2 ) 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.
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 ω pe t ∼ 5000, the particle is captured by a current layer where it experiences an abrupt acceleration, from γ ∼ Γ 0 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 γ ∼ Γ 0 σ 0 = 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 secondorder Fermi acceleration. 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 ν/ν 0 Total 500 < x/de < 3000 5000 < x/de < 7500 ν −0.6 Fig. 10. Total synchrotron spectrum for σ 0 = 30 at ω pe t = 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 a , b , and c refer to Figure 11. Figure 10 shows the instantaneous synchrotron spectrum emitted by the pairs in the σ 0 = 30 simulation at ω pe t = 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 A&A proofs: manuscript no. shock_pwn Fig. 11. Spatial distribution of the synchrotron radiation flux integrated in the frequency bands defined in Figure 10: a (ν/ν 0 < 10 2 ), b (10 2 < ν/ν 0 < 10 6 ), and c (ν/ν 0 > 10 6 ).

Synchrotron radiation
is then given by where r e = e 2 /m e c 2 is the classical radius of the electron, and is the critical synchrotron frequency.B ⊥ 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) Frequencies are normalized by the fiducial synchrotron critical frequency The total synchrotron spectrum presents three main features: (i) a low-energy bump centered around ν/ν 0 ∼ 1, (ii) a plateau at intermediate frequencies (dE/dtdν ∝ ν 0 , 10 1 ν/ν 0 10 4 ) followed by, (iii) a steep power-law decline cutting off at ν/ν 0 ∼ 10 8 . 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 Figure 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, ν/ν 0 10 6 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-, ν/ν 0 < 10 2 , (ii) intermediate-, 10 2 < ν/ν 0 < 10 6 , and (iii) high-frequency band ν/ν 0 > 10 6 labeled as a , b and c 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 b , 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 Figure 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 ( c ), 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 becauseB ⊥ ≈ 0.

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 form-ing 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 powerlaw 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(Zhdankin et al. , 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 σ 0 = 30, or a mean upstream magnetizationσ 0 ≈ 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(Porth et al. , 2014. The spectrum extends from Γ 0 to ∼ Γ 0 σ 0 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 highestenergy part of the particle spectrum. The maximum particle energy increases nearly linearly with time, γ max ∝ t, in contrast to Weibel-dominated shock acceleration where γ max ∝ √ 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-ofplane 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.