Particle Acceleration in Pulsar Wind Nebulae: PIC modelling

TL;DR

This paper reviews how particle-in-cell simulations help understand particle origins and acceleration mechanisms in pulsar wind nebulae, emphasizing magnetic reconnection's role alongside Fermi processes in explaining observed emissions.

Contribution

It highlights the importance of PIC simulations in revealing magnetic reconnection as a key particle acceleration process in PWNe.

Findings

Magnetic reconnection significantly contributes to particle acceleration in PWNe.

PIC simulations support a combined acceleration scenario involving Fermi processes and reconnection.

Reconnection processes explain multi-wavelength emission features of the Crab Nebula.

Abstract

We discuss the role of particle-in-cell (PIC) simulations in unveiling the origin of the emitting particles in PWNe. After describing the basics of the PIC technique, we summarize its implications for the quiescent and the flaring emission of the Crab Nebula, as a prototype of PWNe. A consensus seems to be emerging that, in addition to the standard scenario of particle acceleration via the Fermi process at the termination shock of the pulsar wind, magnetic reconnection in the wind, at the termination shock and in the Nebula plays a major role in powering the multi-wavelength signatures of PWNe.

Equations

$$\frac{\partial f}{\partial t}+\frac{\mathbf{p}}{\gamma m}\cdot \frac{\partial f}{\partial \mathbf{r}}+q\left(\mathbf{E}+\frac{\mathbf{v}\times \mathbf{B}}{c}\right)\cdot \frac{\partial f}{\partial \mathbf{p}}=0,$$

$$f\left(\mathbf{r},\mathbf{p},t\right)\approx \sum _{k=1}^N{w}_k\delta \left(\mathbf{r}-{\mathbf{r}}_k(t)\right)\delta \left(\mathbf{p}-{\mathbf{p}}_k(t)\right),$$

$$\frac{d\mathbf{u}}{dt}$$

$$\frac{d\mathbf{r}}{dt}$$

$$\frac{{\mathbf{u}}^{n+1/2}-{\mathbf{u}}^{n-1/2}}{\Delta t}=\frac{q{\mathbf{E}}^n}{mc}+\frac{q}{mc}\left(\frac{{\mathbf{u}}^n\times {\mathbf{B}}^n}{{\gamma }^n}\right).$$

$${\mathbf{u}}^{-}$$

$${\mathbf{u}}^{+}$$

$${\mathbf{u}}^{+}={\mathbf{u}}^{-}+{\mathbf{u}}^{-}\times \mathbf{s}+\left({\mathbf{u}}^{-}\times \mathbf{w}\right)\times \mathbf{s},$$

$$\mathbf{w}=\frac{q\Delta t{\mathbf{B}}^n}{2mc{\gamma }^n},\mathbf{s}=\frac{2\mathbf{w}}{1+{\mathbf{w}}^2},{\gamma }^n=\sqrt{1+{\left({\mathbf{u}}^{-}\right)}^2}.$$

$${\mathbf{r}}^{n+1}={\mathbf{r}}^n+c\Delta t\frac{{\mathbf{u}}^{n+1/2}}{{\gamma }^{n+1/2}},$$

$$\rho \left(\mathbf{r}\right)=\sum _{k=1}^N{q}_k{w}_k\delta \left(\mathbf{r}-{\mathbf{r}}_k\right),\mathbf{J}\left(\mathbf{r}\right)=\sum _{k=1}^N{q}_k{w}_k{\mathbf{v}}_k\delta \left(\mathbf{r}-{\mathbf{r}}_k\right),$$

$$\rho \left({\mathbf{r}}_i\right)=\sum _{k=1}^N{q}_k{w}_{k}S\left({\mathbf{r}}_i-{\mathbf{r}}_k\right),\mathbf{J}\left({\mathbf{r}}_i\right)=\sum _{k=1}^N{q}_k{w}_k{\mathbf{v}}_{k}S\left({\mathbf{r}}_i-{\mathbf{r}}_k\right),$$

$${\mathbf{J}}_{i,j}$$

$${\mathbf{J}}_{i+1,j}$$

$${\mathbf{J}}_{i,j+1}$$

$${\mathbf{J}}_{i+1,j+1}$$

$$a_k=\frac{x_k-x_i}{dx},b_k=\frac{y_k-y_j}{dy}$$

$$\frac{\partial \mathbf{E}}{\partial t}$$

$$\frac{\partial \mathbf{B}}{\partial t}$$

$$\frac{{\left(B_{\mathrm{z}}\right)}_{i+1/2,j+1/2,k}^{n+1/2}-{\left(B_{\mathrm{z}}\right)}_{i+1/2,j+1/2,k}^{n-1/2}}{\Delta t}=-c\frac{{\left(E_{\mathrm{y}}\right)}_{i+1,j+1/2,k}^n-{\left(E_{\mathrm{y}}\right)}_{i,j+1/2,k}^n}{\Delta x}$$

$$+c\frac{{\left(E_{\mathrm{x}}\right)}_{i+1/2,j+1,k}^n-{\left(E_{\mathrm{x}}\right)}_{i+1/2,j,k}^n}{\Delta y},$$

$${\left(\frac{c\Delta t}{\Delta x}\right)}^2(1D)$$

$${\left(c\Delta t\right)}^2\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)$$

$${\left(c\Delta t\right)}^2\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}+\frac{1}{\Delta z^2}\right)$$

$${\gamma }_{sync,e}\simeq 3.5\times 10^8{L}_{sd,38.5}^{1/6}{\dot{N}}_{40}^{-1/3}{ϵ}_{B,-2.5}^{-1/3}{R}_{\mathrm{TS},17.5}^{1/3}.$$

$${\gamma }_{conf,e}\simeq 1.9\times 10^7{L}_{sd,38.5}^{3/4}{\dot{N}}_{40}^{-1/2},$$

$$h{\nu }_{conf,e}\simeq 0.17L_{sd,38.5}^2{\dot{N}}_{40}^{-1}{ϵ}_{B,-2.5}^{1/2}{R}_{\mathrm{TS},17.5}^{-1}\mathrm{keV}$$

$$\frac{\lambda }{r_{L,hot}}\simeq 4\pi \kappa \frac{R_{LC}}{R_{TS}},$$

$${\mathbf{B}}_L(r\le r_{\mathrm{j}})\propto J_1(r{\alpha }_0){\mathbf{e}}_{\varphi }+J_0(r{\alpha }_0){\mathbf{e}}_z,$$

Full text

0.1 Introduction

In recent years, multi-dimensional magnetohydrodynamic (MHD) models of Pulsar Wind Nebulae (PWNe) have been able to reproduce the nebular morphology down to intricate details (see contributions by A. Mignone, L. del Zanna and O. Porth in this volume). In order to compare the results of MHD simulations to the multi-wavelength observations of PWNe (most notably, of the prototypical Crab Nebula), it is usually assumed that the termination shock of the pulsar wind — where the ram pressure of the ultra-relativistic wind emanating from the pulsar balances the thermal pressure of the surrounding nebula — is an efficient site of particle acceleration. This assumption cannot be tested directly within the MHD framework (which bears no information on the properties of the accelerated particles), but it requires fully-kinetic particle-in-cell (PIC) simulations. By capturing the interplay of charged particles and electromagnetic fields from first principles, PIC simulations allow to identify potential locations of particle acceleration in PWNe. On the other hand, due to computational constraints, PIC simulations are typically confined to a local description of the system, on scales much smaller than the nebular size. It is only by integrating the PIC results with a global MHD model of the nebula that we can properly reproduce the multi-wavelength signatures of PWNe.

In this chapter, we summarize the role of PIC simulations in understanding the origin of high-energy particles in PWNe. In Sect. 0.2 we describe the basics of the PIC method. In Sect. 0.3 we summarize the implications of PIC results on the quiescent and flaring emission from PWNe. We conclude in Sect. 0.4.

0.2 The particle-in-cell technique

This section is intended to provide a brief overview of the most common methods and algorithms used in explicit PIC codes. A detailed presentation of this technique can be found in Hockney and Eastwood (1988); Birdsall and Langdon (1991).

0.2.1 Collisionless plasmas

A necessary condition for non-thermal particle acceleration is the absence of Coulomb collision in the plasma of interest. This is the case for most high-energy astrophysical systems, and in particular pulsar wind nebulae, where plasmas are very diluted. Roughly speaking, a plasma can be considered as “collisionless” if the frequency of Coulomb collision ($\nu$) is much smaller than the plasma frequency, $\omega_{\rm pe}\gg\nu$. This condition implies that the number of particles per Debye sphere must be large, i.e., $N_{\rm D}\gg 1$. The dynamics of individual particle is driven by collection plasma phenomena rather than binary collisions at the sub-Debye length and plasma frequency scales, or simply referred below as the “kinetic” scale. As we will see in Section 0.3, these microscopic scales are involved in the particle acceleration processes and, thus, they must be well-resolved by simulations in contrast to the magnetohydrodynamic approach. To obtain meaningful astrophysical results, particle-based simulations must also capture large scale features, i.e., system size and long integration time.

The evolution of a collisionless plasma is governed by the Vlasov equation

[TABLE]

where $f\equiv d{\rm N}/d\mathbf{r}d\mathbf{p}$ is the particle distribution function defined in 6D phase space $\left(\mathbf{r},\mathbf{p}\right)$ and 1D in time, with $\mathbf{r}$ is the position and $\mathbf{p}=\gamma m\mathbf{v}$ is the momentum, and $q$ is the electric charge. Along with Maxwell’s equations for the fields ($\mathbf{E}$ and ${\mathbf{B}}$), this is the full set of equations to model a collisionless plasma from first principles.

0.2.2 The particle approach

Analytical solutions to the Vlasov equation are known for a few idealized situations only. In most cases, it must be solved numerically. They are at least two ways to solve this equation. In the first approach, phase space is treated as a continuous fluid and Vlasov equation is solved directly using semi-Lagrangian or Eulerian methods (Cheng and Knorr, 1976; Elkina and Büchner, 2006). This method has the advantage to be insensitive to particle noise, and hence can capture well weak plasma phenomena and broad particle distribution functions. In theory this is the most appropriate approach to follow, but in practice the usage of Vlasov codes is currently limited due to prohibitive numerical costs for multidimensional problems (6D). The second approach is the PIC method which is the main focus of this chapter.

In PIC, Vlasov equation is solved indirectly by integrating discrete particle trajectories. This approach is equivalent to the direct method, and an easy way to see this is to rewrite Vlasov equation as a usual advection equation: $\partial f/\partial t+\boldsymbol{\nabla}_{\mathbf{r,p}}\cdot\left(f\mathbf{U}\right)=\mathbf{0}$, where $\boldsymbol{\nabla}_{\mathbf{r,p}}=\left(\partial/\partial\mathbf{r},\partial/\partial\mathbf{p}\right)$ and $\mathbf{U}=\left(\mathbf{p}/\gamma m,q\left(\mathbf{E}+\mathbf{v}\times\mathbf{B}/c\right)\right)$. Thus, using the methods of characteristics, this first-order partial differential equation can be rewritten as a sets of ordinary differential equations (Newton’s law) along characteristic curves which corresponds here to particle trajectories. For point-like particles, the particle distribution function is then approximated as

[TABLE]

where $\delta$ is the Dirac delta function and $w_{k}$ is the particle weight. The number of particles must be very high for a good sampling of phase space and to be close to the exact solution of Vlasov equation. In practice, however, this number will be limited by computing resources and is always much smaller than the number of particles contained in real plasmas. To overcome this difficulty, a PIC particle represents a large number (given by the weight $w_{k}$) of physical particles that would follow the same trajectory in phase space (with the same $q/m$ ratio). For this reason, the simulation particles are usually called “macroparticles”.

Even though the plasma is collisionless, particles feels each other via long-range interactions. Summing over all particle-particle binary interactions, i.e. $N\left(N-1\right)/2\approx N^{2}$, is numerically expansive and hard to implement. Instead, in PIC, particles do not feel each other directly but via the electromagnetic fields known on the grid which result from the evolution of the plasma. In this case, the number of operations scales as the number of particles $N$ instead of $N^{2}$.

The PIC method has become increasingly popular in high-energy astrophysics to model non-thermal particle acceleration phenomena. PIC codes are much cheaper in comparison to Vlasov codes, and they are also conceptually simple, robust and easy to implement and parallelize efficiently to a large number of cores. This simplicity comes at the cost of significant particle noise which can lead to poor sampling of the particle distribution (e.g., steep power-law tails), difficulty in capturing subtle or weak phenomena, artificial collisions, and load-balancing issues in parallel computing.

0.2.3 Main computing procedures in PIC

Figure 1 describes the three main operations performed per timestep $\Delta t$ of an explicit PIC code: (i) Solve Newton’s equation for each particle to evolve velocities and positions (ii) Collect charge and current densities from all particles and deposit them on the grid, and (iii) Solve Maxwell’s equations to update the fields on the grid. Below is a brief technical description of each step:

Step 1: Particle push

The set of equations to solve are

[TABLE]

where $\mathbf{u}=\gamma\mathbf{v}/c$ is the particle 4-velocity vector divided by the speed of light and $\gamma=1/\sqrt{1-\left(v/c\right)^{2}}$ is the Lorentz factor. One of the most successful and most common method used in PIC to solve Newton’s equation is the Boris push111Other efficient methods exists as for instance the particle pusher developed by Vay (2008).. It has all the desirable numerical features one might think of: it is fast, stable, second order accurate, and conserves well the particle energy. The algorithm is based on the usual leapfrog integration method, i.e., 4-velocities $\mathbf{u}$ and positions $\mathbf{r}$ are staggered in time by half a timestep (Figure 2). If particle positions and fields are known at time $t^{n}$ ($\mathbf{r}^{n}$, $\mathbf{E}^{n}$, $\mathbf{B}^{n}$) and velocities at time $t^{n-1/2}$ ($\mathbf{u}^{n+1/2}$), the finite-difference time-centered expression of Eq. (3) is

[TABLE]

Now, the trick is to rewrite $\mathbf{u}^{n}$ appearing on the right-hand side of the equation as $\mathbf{u}^{n}=\left(\mathbf{u}^{n-1/2}+\mathbf{u}^{n+1/2}\right)/2$. Assuming that $\mathbf{E}^{n}$ and $\mathbf{B}^{n}$ are known, $\mathbf{u}^{n+1/2}$ can be extracted. It is convenient to define the following intermediate variables

[TABLE]

Then, using Eqs. (6)-(7) and after a few algebraic manipulations one finds

[TABLE]

where

[TABLE]

It is important to notice that the fields appearing in these equations are those felt at the particle position, not at the grid point where the fields are known. The fields must be interpolated to the particle positions. A linear interpolation scheme is usually sufficient. The final step is to update the particle positions

[TABLE]

where $\gamma^{n+1/2}=\sqrt{1+\left(\mathbf{u}^{n+1/2}\right)^{2}}$.

Step 2: Charge and current deposition

To solve Maxwell’s equations, we need the source terms $\rho$ and $\mathbf{J}$ that are given by the particles. In a continuous space, these macroscopic quantities can be recovered by summing over the contribution from all particles

[TABLE]

where $q_{k}$, $\mathbf{v}_{k}$ are respectively the electric charge and the 3-velocity of the particle $k$. In PIC, charges and currents from the particles must be collected and dispatched among the nearest grid points. Charge and current densities at the grid point $\mathbf{r}_{i}$ can be written as

[TABLE]

where $S$ is a shape function which depends on the desired deposition scheme. Even though the particles are point-­like, they have an effective size that is felt through the deposition of currents on the grid.

Figure 3 shows the example of a first order linear deposition method in a 2D Cartesian grid cell (or area-weighting method). The contributions from all the particles contained in the cell $(x_{i},y_{j})$ to the current $\mathbf{J}$ are given by

[TABLE]

where

[TABLE]

are the usual bilinear interpolation coefficients. In this particular example, particles have a triangular shape.

Step 3: Fields evolution

The last step is to update the fields on the grid. In principle, one needs solely to solve the time-dependent equations

[TABLE]

because the current density is already given by the particles as we have seen in the previous paragraph. The other two should be automatically satisfied, but this is not necessarily true due to truncation errors in the discretization of space and time derivatives. The current deposition procedure does not always guarantee charge conservation222The total particle charge is conserved, but not necessarily the charge deposited on the grid. (i.e., $\boldsymbol{\nabla}\cdot\mathbf{J}\neq-\partial\rho/\partial t$), but some solutions exist to enforce it to machine roundoff precision (Villasenor and Buneman, 1992; Esirkepov, 2001). Alternatively, Poisson equation should be solved to correct the electric field to make sure $\boldsymbol{\nabla}\cdot\mathbf{E}=4\pi\rho$. Parabolic and hyberbolic divergence cleaning methods also exist in the literature (Marder, 1987; Munz et al., 2000).

The finite difference time domain (FDTD) method proposed by Yee (1966) for solving the time-dependent Maxwell equations enforces $\boldsymbol{\nabla}\times\mathbf{B}=0$ to machine roundoff precision. This is the most commonly used method in explicit PIC codes. Like the Boris push, the FDTD method combines stability, efficiency and second order accuracy (here in both space and time). To achieve this, fields must be staggered in time and in space. Figure 4 shows the order in time (top panel), as well as the spatial configuration of the fields within a Cartesian cell333For a spherical geometry, see Holland 1983; Cerutti et al. 2015, 2016; Belyaev 2015. in 2D (bottom-left panel) and in 3D (bottom-right panel). For illustrative purposes, within this framework the $z$-component of Eq. (19) is

[TABLE]

where $\Delta x,~{}\Delta y$ are the spatial step size along $x$ and $y$. The FDTD method is stable under the usual Courant-­Friedrichs-­Lewy condition, i.e.,

[TABLE]

This is a purely numerical requirement, but physics imposes other constraints on the size of the steps, namely that the Debye length and the plasma frequency are well resolved by the code ($\Delta x/\Lambda_{\rm D}\ll 1$ and $\omega_{\rm pe}\Delta t\ll 1$), the latter condition being more stringent.

0.2.4 Boundary conditions

Periodic boundary conditions are robust, easy to implement and physically useful in many case studies but there are not always appropriate. Below is a brief description of some other boundary conditions usually employed in PIC simulations.

Particles

• •

A perfectly reflective wall scatters the particle with no loss of momentum and energy. It can be useful in the context of a perfectly conducting wall for the fields (see below).

• •

In the context of an absorbing wall or of an open boundary, it is appropriate to absorb particles at the boundary. In such a case, particles are simply removed from the simulation.

• •

Conversely, new particles can be injected in simulations. This can be motivated by the physics involved in the problem, as for instance in pulsars where pair creation is important (Timokhin and Arons, 2013; Chen and Beloborodov, 2014; Philippov et al., 2015), or by the numerical setup if for example there is an inflow from one side of the box as in PIC simulations of relativistic shocks (Spitkovsky, 2008b; Sironi and Spitkovsky, 2011a). In this case, an expanding box with an injector receding at the speed of light can be desirable to reduce numerical cost (see 0.3 below).

Fields

• •

Perfectly conducting walls allow to reflect electromagnetic waves. They are easily implemented by applying the usual boundary conditions, namely that the tangential component of $\mathbf{E}$ and the perpendicular component of $\mathbf{B}$ vanishes at the interface. Semi-reflective medium can also be easily coded using surface current and charge densities.

• •

It is sometimes useful to absorb all electromagnetic waves leaving the box to simulate an open boundary, as for instance in pulsar winds (Cerutti et al., 2015; Belyaev, 2015). In this case, the open boundary is coated with an absorbing medium of several cell thick where resistive terms are added to Maxwell’s equations

[TABLE]

where $\lambda$ and $\lambda^{\star}$ are artificial electric and magnetic “conductivities”. The transition between the working domain and the absorbing layer should be gradual to avoid undesired reflections at the boundary. Conductivities usually are increasing function from the inner edge to the outer edge of the damping layer to make sure waves are completely absorbed. Eqs. (23-23) are valid for 1D layer. A perfectly matched layer is a generalization of these formulae to a multidimensional damping layer (Berenger, 1994, 1996). In this framework, fields must be split into two subcomponents and hence the number of equations to solve is doubled (up to 12 in 3D, for an application to pulsars see e.g., Kalapotharakos and Contopoulos 2009 in the context of force-free MHD simulations).

0.2.5 Parallelization

PIC codes must be efficiently parallelized to model large system size and long integration time to have meaningful astrophysical applications. A common practice is to use the domain decomposition technique. It consists in dividing the computational box into smaller domains where one or more cores are assigned. Each CPU goes through the main steps described in Sect. 0.2.3 every timestep and exchanges information with the neighbouring processes to send particle and field data at the interface between subdomains. Communications between an arbitrary number of processes are done thanks to the Message Passing Interface (MPI) library. PIC codes scale well to a large number of CPUs, today at least up to $\sim 10^{6}$ processes (see Fig. 5). These scaling plots are usually done under ideal conditions, and do not necessarily reflect problem-dependent loss of performance. In PIC, a poor load balancing severely slows down a simulation. If, for some reason, there is a concentration of particle in a few subdomains, only few processors will have to push a lot of particles while the others remain idle444Typically, pushing particles and depositing currents take $90\%$ of the computing time (without communications), this is the reason why load-balancing is so critical in PIC.. The way how the domain is decomposed for a given setup can usual make a big difference. Hybrid codes combining MPI and OpenMP, variable particle weighting, or dynamical changes of the domain decomposition are other solutions to have better performances.

0.3 Application to pulsar wind nebulae

In this section, we describe the main results obtained with PIC simulations on the efficiency of particle acceleration in PWNe. This section is divided into two parts: at first, we investigate particle acceleration at the termination shock of the pulsar wind, which is usually invoked to power the nebular quiescent emission (Sect. 0.3.1); then, we will discuss the origin of the gamma-ray flares detected from the Crab Nebula in the last few years (Sect. 0.3.2), focusing on the mechanisms that might explain such extreme particle acceleration events.

0.3.1 The quiescent emission

The Crab Nebula, our best laboratory for high energy astrophysics, has been observed over the entire electromagnetic spectrum from $\lesssim 100\,{\rm MHz}$ to $\gtrsim 1.5\,{\rm TeV}$ (see several contributions in this volume). Efficient acceleration of particles at the termination shock is required to explain its broadband spectrum. However, the flatness of the radio spectrum ($F_{\nu_{r}}\propto\nu^{-0.3}$, Bietenholz et al., 1997) is hard to reconcile with the steeper optical and X-ray slope ($F_{\nu_{X}}\propto\nu^{-1.1}$, Mori et al., 2004), unless the electron distribution is more complicated than a single power law. In fact, the radio band would require a distribution of emitting particles with a power-law slope $p=-d\log N/d\log\gamma\simeq 1.6$, whereas $p\gtrsim 2$ would be needed for the optical and X-ray emission. Even more fundamentally, how the pulsar wind termination shock can accelerate particles to the required “non-thermal” energies (i.e., well beyond the “thermal” peak of a Maxwellian distribution) is still an unsolved problem.

Particle acceleration in shocks is usually attributed to the Fermi process, where particles are energized by bouncing back and forth across the shock. Despite its importance, the Fermi process is still not understood from first principles. The highly nonlinear coupling between accelerated particles and magnetic turbulence — which is generated by the particles, and at the same time governs their acceleration — is extremely hard to incorporate in analytic models, and can be captured only with ab initio PIC simulations (for a review of the Fermi process in relativistic shocks, see Sironi et al., 2015a).

As we describe below, the efficiency of the Fermi process depends critically on the shock properties, e.g., composition, magnetization (i.e., the ratio $\sigma$ between the Poynting flux and the kinetic energy flux of the pre-shock flow) and magnetic obliquity (i.e., the angle $\theta$ between the upstream magnetic field and the shock direction of propagation).555In the limit $\gamma_{0}\gg 1$ of ultra-relativistic shocks, as appropriate for PWNe, the efficiency of the Fermi process does not depend on the shock Lorentz factor $\gamma_{0}$ (Sironi and Spitkovsky, 2009, 2011b; Sironi et al., 2013). Pulsar winds are thought to be dominated by electron-positron pairs (Bucciantini et al., 2011). MHD models of PWNe require $\sigma\gtrsim 0.01-0.1$ in order to reproduce the morphology of the Crab jet/plume. Finally, polarization measurements indicate that the nebular magnetic field should be toroidal around the symmetry axis of the system, so that the termination shock is “perpendicular” (i.e., with the field orthogonal to the flow direction).

PIC simulations of perpendicular magnetized shocks show negligible particle acceleration (Gallant et al., 1992; Hoshino, 2008; Sironi and Spitkovsky, 2009, 2011b; Sironi et al., 2013). Here, due to the lack of significant self-generated turbulence, charged particles are forced to slide along the background field lines, whose orientation prohibits repeated crossings of the shock. This inhibits the Fermi process, and in fact the particle distribution behind perpendicular shocks is purely thermal.

In summary, PIC simulations have shown that the shock configurations which are apparently most relevant for PWNe (i.e., ultra-relativistic magnetized perpendicular shocks) do not naturally result in efficient particle acceleration. This is in sharp contrast with the pronounced non-thermal signatures of the quiescent emission of PWNe. However, one key ingredient of the PIC results summarized above is that the pre-shock magnetic field direction stays uniform throughout the timespan of the simulations. This is generally not the case in pulsar winds. If the rotational and magnetic axes of the central pulsar are misaligned, around the equatorial plane the wind consists of toroidal stripes of alternating magnetic polarity, separated by current sheets of hot plasma. It is still a subject of active research whether the alternating stripes will dissipate their energy into particle heat ahead of the termination shock, or whether the wind remains dominated by Poynting flux till the termination shock (Lyubarsky and Kirk, 2001; Kirk and Skjæraasen, 2003; Sironi and Spitkovsky, 2011a). If the stripes are dissipated far ahead of the termination shock, the upstream flow is weakly magnetized and the pulsar wind reaches a terminal Lorentz factor (in the frame of the nebula) $\gamma_{0}\sim L_{sd}/m_{e}c^{2}\dot{N}\simeq 3.7\times 10^{4}L_{sd,38.5}\dot{N}_{40}^{-1}~{},$ where $L_{sd}\equiv 3\times 10^{38}L_{sd,38.5}\nobreak{\mathrm{\;ergs\,s^{-1}}}$ is the spin-down luminosity of the Crab, and $\dot{N}=10^{40}\dot{N}_{40}\nobreak{\mathrm{\;s^{-1}}}$ is the particle flux entering the nebula, including the radio-emitting electrons (Bucciantini et al., 2011).

The two subsections below cover the two potential fates of the pulsar striped wind: at first, we investigate the physics of particle acceleration in a weakly magnetized shock (i.e., assuming that the alternating stripes have dissipated their magnetic energy far ahead of the termination shock); then, we assume that the stripes persist until the termination shock.

The termination shock of a weakly magnetized wind

Weakly magnetized ultra-relativistic shocks are mediated by electromagnetic plasma instabilities (the so-called Weibel instability, Weibel, 1959; Medvedev and Loeb, 1999; Gruzinov and Waxman, 1999). These instabilities build up a magnetic barrier, up to a level666The parameter $\epsilon_{B}$ denotes the magnetization of the turbulence, $\epsilon_{B}\,=\,\delta B^{2}/8\pi\gamma_{0}\rho_{0}c^{2}$, where $\delta B$ is the fluctuating magnetic field and $\rho_{0}$ is the mass density of the pre-shock flow. This should not be confused with the magnetization $\sigma=B_{0}^{2}/4\pi\gamma_{0}\rho_{0}c^{2}$, which quantifies the strength of the pre-existing ordered upstream field $B_{0}$. $\epsilon_{B}\,\sim\,10^{-2}-10^{-1}$, sufficient to deflect strongly the incoming particles and thus mediate the shock transition. The instability — triggered by a stream of shock-reflected particles propagating ahead of the shock — generates filamentary magnetic structures in the upstream region (Fig. 6), which in turn scatter the particles back and forth across the shock, mediating Fermi acceleration.

Such shocks do self-consistently accelerate particles up to nonthermal energies, via the Fermi process (Spitkovsky, 2008a, b; Martins et al., 2009; Haugbølle, 2011; Sironi et al., 2013). The accelerated particles populate in the downstream region a power-law tail $dN/d\gamma\propto\gamma^{-p}$ with a slope $p\sim 2.5$, that contains $\sim 3\%$ of the particles and $\sim 10\%$ of the flow energy.

The particle energy spectrum extends over time to higher and higher energies, as shown in Fig. 7. For electron-positron flows, as appropriate for pulsar winds, the maximum post-shock particle Lorentz factor increases with time as $\gamma_{max}\sim 0.5\,\gamma_{0}\,(\omega_{\rm p}t)^{1/2}$ (Sironi et al., 2013).777This scaling is shallower than the so-called (and commonly assumed) Bohm limit $\gamma_{max}\propto t$, and it naturally results from the small-scale nature of the Weibel turbulence generated in the shock layer (see Fig. 6). The plasma frequency $\omega_{\rm p}$ can be computed from the number density ahead of the termination shock, which is $n_{{\rm TS}}=\dot{N}/(4\pi R_{\rm TS}^{2}c)$, assuming an isotropic particle flux. Here, $R_{\rm TS}\equiv 3\times 10^{17}R_{\rm TS,17.5}\nobreak{\mathrm{\;cm}}$ is the termination shock radius. Balancing the acceleration rate with the synchrotron cooling rate in the self-generated Weibel fields, the maximum electron Lorentz factor is

[TABLE]

A stronger constraint comes from the requirement that the diffusion length of the highest energy electrons be smaller than the termination shock radius (i.e. a confinement constraint). Alternatively, the acceleration time should be shorter than $R_{\rm TS}/c$, which yields the critical limit

[TABLE]

which is generally more constraining than the cooling-limited Lorentz factor $\gamma_{sync,e}$. The corresponding synchrotron photons will have energies

[TABLE]

which are apparently too small to explain the X-ray spectrum of the Crab, extending to energies beyond a few tens of MeV.

At face value, Fermi acceleration at the termination shock of PWNe is not a likely candidate for producing X-ray photons via the synchrotron process. Yet, the steady-state hard X-ray and gamma-ray spectra of PWNe do look like the consequences of Fermi acceleration — particle distributions with $p\simeq 2.4$ are a natural prediction of the Fermi process in ultra-relativistic shocks (Kirk et al., 2000; Achterberg et al., 2001; Keshet and Waxman, 2005). In this regard, we argue that the wind termination shock might form in a macroscopically turbulent medium, with the outer scale of the turbulence driven by the large-scale shear flows in the nebula (Komissarov and Lyubarsky, 2004; Del Zanna et al., 2004; Camus et al., 2009). If the large-scale motions drive a turbulent cascade to shorter wavelengths, back-scattering of the particles in this downstream turbulence, along with upstream reflection by the transverse magnetic field of the wind, might sustain Fermi acceleration to higher energies.

An alternative mechanism leading to particle acceleration to higher energies may be connected to the accelerator behind the recently discovered gamma-ray flares in the Crab Nebula (see Sect. 0.3.2). Runaway acceleration of electrons and positrons at reconnection X-lines, a linear accelerator, may inject energetic beams into the shock, with the mean energy per particle approaching the whole open field line voltage, $\gtrsim 10^{16}\nobreak{\mathrm{\;V}}$ in the Crab (Arons, 2012), as required to explain the Crab GeV flares. This high-energy population can drive cyclotron turbulence when gyrating in the shock-compressed fields, and resonant absorption of the cyclotron harmonics can accelerate the electron-positron pairs in a broad spectrum, with maximum energy again comparable to the whole open field line voltage (Hoshino et al., 1992; Amato and Arons, 2006).

The termination shock of a strongly magnetized striped wind

Assuming that the stripes survive until the termination shock, we now describe the physics of particle acceleration if the pre-shock flow carries a strong magnetic field of intensity $B_{0}$, oriented perpendicular to the shock direction of propagation and alternating with wavelength $\lambda$.888The wavelength $\lambda$ of the striped wind equals $c\,P$, where $P$ is the pulsar period. Although the magnetic field strength in the wind is always $B_{0}$, the wavelength-averaged field $\langle B_{\phi}\rangle_{\lambda}$ can vary from zero up to $B_{0}$, depending on the relative widths of the regions of positive and negative field (see the sketch in Fig. 8). In pulsar winds, one expects $\langle B_{\phi}\rangle_{\lambda}=0$ only in the equatorial plane (where the stripes are symmetric), whereas $|\langle B_{\phi}\rangle_{\lambda}|/B_{0}\rightarrow 1$ at high latitudes. As a proxy for latitude, we choose $\alpha=2\langle B_{\phi}\rangle_{\lambda}/(B_{0}+|\langle B_{\phi}\rangle_{\lambda}|)$, which varies between zero and unity.

At the termination shock, the compression of the flow forces the annihilation of nearby field lines, a process known as driven magnetic reconnection (Lyubarsky, 2003; Sironi and Spitkovsky, 2011a, 2012). As shown in Fig. 9, magnetic reconnection erases the striped structure of the flow (panel (a)), and transfers most of the energy stored in the magnetic fields (panel (d)) to the particles, whose distribution becomes much hotter behind the shock (see panel (f), for $x\lesssim 1000\,c/\omega_{\rm p}$). As a result of field dissipation, the average particle energy increases by a factor of $\sigma$ across the shock, regardless of the stripe width $\lambda$ or the wind magnetization $\sigma$ (as long as $\alpha\lesssim 0.1$). The reconnection process manifests itself as characteristic islands in density (panel (c)) and magnetic energy (panel (e)), separated by X-points where the magnetic field lines tear and reconnect.

The incoming particles are accelerated by the reconnection electric field at the X-points and, in the post-shock spectrum, they populate a broad distribution (red line in panel (g)), extending to much higher energies than expected in thermal equilibrium (dotted line). For the parameters studied in Fig. 9, the slope of the non-thermal tail is $p\simeq 1.5$ (dashed line in panel (g)), harder than what the Fermi process normally gives in relativistic shocks.999Hard particle spectra are found to be a generic by-product of magnetic reconnection in the relativistic regime appropriate for pulsar winds (e.g., Sironi and Spitkovsky, 2014; Guo et al., 2014; Melzani et al., 2014; Werner et al., 2016; Sironi et al., 2015b, 2016). While efficient field dissipation (and so, efficient transfer of field energy to the particles) occurs irrespective of the wind properties (if $\alpha\lesssim 0.1$), the width of the downstream particle spectrum is sensitive to the stripe wavelength and the wind magnetization through the combination $\lambda/r_{L,hot}$, namely the stripe wavelength measured in units of the post-shock particle Larmor radius (i.e., after dissipation has taken place, and the mean particle energy has increased by a factor of $\sigma$). A Maxwellian-like spectrum is obtained for $\lambda/r_{L,hot}\lesssim$ a few tens, whereas in the limit $\lambda/r_{L,hot}\gg 1$ the spectrum approaches a broad power-law tail of index $1<p<2$, extending from $\gamma_{\rm min}\simeq\gamma_{0}$ up to $\gamma_{\rm max}\simeq\gamma_{0}\sigma^{1/(2-p)}$.

The particles are accelerated primarily by the reconnection electric field at the X-points, rather than by bouncing back and forth across the shock, as in the standard Fermi mechanism. Quite surprisingly, the Fermi process can still operate along the equatorial plane of the wind, where the stripes are quasi-symmetric ($\alpha\lesssim 0.01$). Here, the highest energy particles accelerated by the reconnection electric field can escape ahead of the shock, and be injected into a Fermi-like acceleration cycle. In the post-shock spectrum, they populate a power-law tail with slope $p\simeq 2.5$, that extends beyond the hard component produced by reconnection.

At higher latitudes, the presence of a non-negligible stripe-averaged field $\langle B_{\phi}\rangle_{\lambda}$ inhibits the Fermi process, in analogy to the case of perpendicular magnetized shocks (with uniform fields) discussed at the beginning of this section. The efficiency of particle acceleration via shock-driven reconnection is also affected, as we show in Fig. 10. For $\alpha\lesssim 0.1$ (i.e., relatively close to the equatorial plane), the field is efficiently dissipated, and the shape of the spectrum is nearly independent of latitude. For $\alpha\gtrsim 0.3$, the post-shock particle spectrum consists of two components. The low-energy peak comes from cold plasma with mean Lorentz factor $\sim\gamma_{0}$, whereas the high-energy part is populated by hot particles that gained energy from field dissipation, so that their mean Lorentz factor is now $\sim\gamma_{0}\sigma$. As $\alpha$ increases, the fraction of upstream Poynting flux available for dissipation decreases, which explains why the high-energy component in the spectra of Fig. 10 gets de-populated, at the expense of the low-energy part. The limit $|\alpha|\rightarrow 1$ (yellow curve for $\alpha=0.95$) approaches the result expected for an unstriped wind (purple line).

Based on these findings, one could interpret the optical and X-ray signatures of the Crab, which require a particle spectrum with $p\simeq 2.5$, as synchrotron emission from the particles that are Fermi-accelerated close to the equatorial plane of the wind. In addition, the spectral index required for the radio spectrum of the Crab ($p\simeq 1.5$) could naturally result from the broad hard component of particles accelerated by shock-driven reconnection. However, the particle spectrum in the simulations approaches the hard tail required by the observations only when the combination $\lambda/r_{L,hot}$ exceeds a few tens (for smaller values, the spectrum is a narrow thermal-like distribution). At the termination shock of the pulsar wind ($R=R_{TS}$) we have

[TABLE]

where $R_{LC}=c/\Omega$ is the light cylinder radius ($\Omega=2\pi/P$ is the pulsar rotational frequency), and $\kappa$ is the so-called multiplicity in the wind (i.e., the ratio of the actual density to the Goldreich-Julian density, Goldreich and Julian (1969)). For the Crab, $R_{TS}\simeq 5\times 10^{8}R_{LC}$ (Hester et al., 2002) and most available models estimate $\kappa\simeq 10^{4}-10^{6}$ (Bucciantini et al., 2011). Based on our findings, the resulting value of $\lambda/(r_{L}\sigma)\lesssim 0.01$ would yield a Maxwellian-like spectrum, at odds with the wide flat spectrum required by observations. If radio-emitting electrons are accelerated at the termination shock of pulsar winds via magnetic reconnection, a revision of the existing theories of pulsar magnetospheres is required.

0.3.2 The flaring emission

In recent years, the Fermi and AGILE satellites have detected a number of hours- to week-long flares at GeV energies, which surprisingly falsify the widely-believed “standard candle” nature of the high-energy Crab emission. During these events the Crab nebula gamma-ray flux above 100 MeV exceeded its average value by a factor of several or higher (Tavani et al., 2011; Abdo et al., 2011; Buehler et al., 2012), while at other wavelengths nothing unusual was observed (Weisskopf et al., 2013). The observed gamma-ray flares happen with a cadence of $\sim 1\,{\rm year}$ (e.g., Buehler and Blandford, 2014) and there are no associated pulsar timing glitches. Variability on timescales as short as a few hours has been reported. The peak isotropic luminosity is roughly $10^{36}\,{\rm ergs/s}$ and the energy radiated is $\sim 10^{41}\,{\rm ergs}$. The flares and secular observations (Wilson-Hodge et al., 2011) demonstrate that the energy conversion is intermittent and that the mechanism can be locally cataclysmic.

The flare properties suggest that an extreme accelerator is at work. The typical decay time of the flaring episodes, which is attributed to synchrotron cooling, together with the $\sim$ GeV peak frequency, allows to solve simultaneously for the magnetic field strength $\sim 5$ mG (as compared to the nebula-averaged $\sim 200\,{\rm\mu G}$) and for the extreme energy of the emitting particles $\sim$ PeV. From the $\sim 10$-hour rise time of the flares, one can estimate the size $\sim 10^{15}$ cm of the emission region. In order to accelerate up to PeV energies within this length, the accelerating electric field needs to be comparable to the inferred magnetic field (i.e., $E\sim B$).

Fermi acceleration at the termination shock of the Crab nebula fails to explain the observed GeV flares (Sironi et al., 2013). In contrast, the requirement that $E\sim B$ is naturally satisfied in reconnection layers, in the relativistic regime where the magnetic energy per particle exceeds its rest mass, or equivalently where the magnetization $\sigma=B_{0}^{2}/4\pi\rho_{0}c^{2}\gg 1$. The reconnection scenario would work best in the most magnetized regions of the nebula, i.e., near the poles and possibly in the jets (Cerutti et al., 2012; Lyubarsky, 2012; Komissarov, 2013; Mignone et al., 2013). Unfortunately, current gamma-ray telescopes do not have the angular resolution to pin down the precise location of the flares within the Nebula. Below, we discuss how PIC simulations help unveiling the role of magnetic reconnection in the Crab Nebula, as the underlying particle accelerator that powers the GeV flares.

Plane-parallel reconnection

In the simplest geometry of magnetic reconnection, the field lines are parallel to a pre-existing current sheet, with opposite polarity on the two sides of the current sheet. We shall call this setup as “plane-parallel reconnection” (see Kagan et al., 2015, for a review). Using 2D and 3D PIC simulations, it has been recently shown that most of the features of the Crab flares can be explained with relativistic plane-parallel reconnection (timescale, energetics, particle and photon spectra). The key arguments in favor of reconnection for the Crab flares are:

• –

The flare spectrum requires an electron power-law population with hard slope $p=-d\log N/d\log\gamma\lesssim 2$, which is not attainable in shocks, but it naturally results from relativistic reconnection (Sironi and Spitkovsky, 2014; Guo et al., 2014; Melzani et al., 2014; Werner et al., 2016; Sironi et al., 2015b, 2016). As shown in Fig. 11, the power-law slope depends on the flow magnetization, being harder for higher $\sigma$ ($p\sim 1.5$ for $\sigma=50$, compare solid and dotted green lines). The slope is steeper for lower magnetizations ($p\sim 4$ for $\sigma=1$, solid and dotted black lines), approaching the result of non-relativistic reconnection, yielding poor acceleration efficiencies (Drake et al., 2010).

• –

The $\sim\,$GeV peak energy of the flares is well above the classical synchrotron “burnoff” limit of $\sim 236\,\eta$ MeV (as measured in the fluid rest frame), which is obtained by balancing acceleration due to an electric field $E=\eta\,B$ with synchrotron cooling losses; unlike in shocks, where $\eta<1$, in the reconnection layer one finds $\eta>1$, thus boosting the synchrotron limit to the observed $\sim$ GeV peak. In a reconnection scenario, this requires the accelerating particles to stay confined within the reconnection layer, where $\eta>1$. Uzdensky et al. (2011) showed analytically that as the particle energy increases, the trajectory gets more and more focused along the electric field, with vanishing cooling losses. This has now been confirmed with PIC simulations (Cerutti et al., 2013, 2014b). In particular, these studies demonstrated that reconnection can accelerate particles above the synchrotron radiation burn-off limit (Guilbert et al., 1983; de Jager et al., 1996) deep inside the reconnection layer where the electric field overcome the magnetic field (see Fig. 12). This result is crucial because it can explain the emission of $>100~{}$MeV synchrotron radiation emitted during every Crab flare, which would be impossible to achieve in ideal MHD.

• –

The short rise time ($\sim 10$ hours) of the flaring episodes can naturally result from the inhomogeneity of the reconnection layer, which is fragmented into a chain of magnetic islands, or plasmoids, as shown in Fig. 13 (in 3D, these plasmoids appear as elongated magnetic flux ropes). The plasmoids are overdense (Fig. 13a), filled with energetic particles and confined by strong fields. The plasma flows into the reconnection layer at $v_{\rm rec}\simeq 0.15\,c$ for $\sigma=10$ (Fig. 13b). The inflow speed is nearly independent of $\sigma$ for larger magnetizations (Sironi et al., 2016), in agreement with analytical models (Lyubarsky, 2005). After entering the sheet, the flow is advected out by the tension force of the reconnected field. The motion in the reconnection exhausts is ultra-relativistic (Fig. 13c), approaching a bulk four-velocity $\Gamma v_{\rm out}\sim\sqrt{\sigma}\,c$, in agreement with the theory (Lyubarsky, 2005). The relativistic bulk motion of the plasmoids in the reconnection layer plays a critical role in enhancing — via Doppler boosting — their emission signatures.

Aside from bulk Doppler beaming, an energy-dependent “kinetic” beaming has also been proposed to explain the extreme time variability of the Crab flares (Cerutti et al., 2013, 2014b). In particular, while low-energy particles are nearly isotropic, at high energies ($\gamma\gtrsim\sigma$) the particles exhibit clear sign of anisotropy with two beams pointing roughly towards the $\pm x$-directions, i.e., along the reconnection exhausts. Hence, the beams are not necessarily pointing along the direction $z$ of the reconnection electric field because the tension of the reconnected field lines pushes the particles away from the X-points in the form of a reconnection outflow towards the magnetic islands. Nonetheless, the direction of the beam of energetic particles is not static: it wiggles rapidly within the $xz$-plane, which results in rapid flares of energetic radiation when the beam crosses the line of sight of a distant observer. Since in the Crab the particles emitting $>100$ MeV synchrotron radiation should be accelerated and radiating over a sub-Larmor timescale, one expects that the highest energy radiation should keep the imprint of the particle anisotropy (regardless of the acceleration process), while the low-energy radiation should be more isotropic (Cerutti et al., 2013, 2014b).

Explosive reconnection

Despite its successes, plane-parallel reconnection does not seem to be fast enough to explain the short rise time of the Crab flares. In particular, recent PIC studies of relativistic reconnection have demonstrated that the reconnection rate (inflow velocity) in 3D simulations of plane-parallel reconnection is significantly lower than in 2D. For a reference magnetization $\sigma=10$ the reconnection rate in 2D is $v_{rec}/c\sim 0.1$, whereas in 3D it is only $v_{rec}/c\sim 0.02$ (Sironi and Spitkovsky, 2014). The slower reconnection rate leads to a weaker accelerating electric field. Moreover, for a given flare duration it translates into a smaller utilised magnetic energy.

To overcome this difficulty, it has recently been proposed that the Crab flares might result from explosive reconnection episodes (a process that has been called “magnetoluminescence” by Blandford et al. (2017), for the rapid conversion of magnetic energy into high-energy particles and then into radiation), accompanying the relaxation of force-free equilibria on dynamical timescales (i.e., corresponding to an effective reconnection rate of $v_{rec}/c\sim 1$). In particular, Nalewajko et al. (2016); Lyutikov et al. (2016); Yuan et al. (2016) have carried out PIC simulations of the relaxation of force-free equilibria in application to the Crab flares.

As a representative case, we consider the configuration of two Lundquist’s force-free cylinders surrounded by uniform magnetic field (Lyutikov et al., 2016),

[TABLE]

Here, $J_{0},\,J_{1}$ are Bessel functions of zeroth and first order and the constant $\alpha_{0}\simeq 3.8317$ is the first root of $J_{0}$. This solution is terminated at the first zero of $J_{1}$, which we denote as $r_{j}$ and hence continued with $B_{z}=B_{z}(r_{j})$ and $B_{\phi}=0$ for $r>r_{j}$. Since the total current of the flux tube is zero, the azimuthal field vanishes at the boundary of the rope, and so the evolution is initially very slow (i.e., the initial configuration is dynamically stable). To speed things up, the ropes are pushed towards each other.

In Fig. 14, we present the 2D pattern of the out-of-plane field $B_{z}$ (left column) and of the in-plane magnetic energy fraction $\epsilon_{B,\rm in}=(B_{x}^{2}+B_{y}^{2})/8\pi nmc^{2}$ (right column; with superimposed magnetic field lines), from a PIC simulation with $\sigma_{\rm in}=43$ (only defined with the in-plane fields) and $\,r_{\rm j}=61\,r_{\rm L,hot}$ (where $r_{\rm L,hot}$ is the Larmor radius of particles heated by reconnection). As the two magnetic ropes slowly approach, driven by the initial velocity push, reconnection is triggered in the plane $x=0$, as indicated by the formation and subsequent ejection of small-scale plasmoids. As a result of reconnection, an increasing number of field lines, that initially closed around one of the ropes, are now engulfing both magnetic islands. Their tension force causes the two ropes to approach and merge on a quick (dynamical) timescale, starting at $ct/\,r_{\rm j}\sim 4.5$ and ending at $ct/\,r_{\rm j}\sim 7.5$ (see that the distance of the rightmost island from the center rapidly decreases, as indicated by the black line in the middle panel of Fig. 15). The tension force drives the particles in the flux ropes toward the center, with a fast reconnection speed peaking at $v_{\rm rec}/c\sim 0.3$ (red line in the middle panel of Fig. 15).101010The reconnection rate is measured to be in the range $v_{\rm rec}/c\sim 0.2-0.5$, which increases with the magnetization and saturates at around 0.5 at high magnetization limit (Lyutikov et al., 2016). The reconnection layer at $x=0$ stretches up to a length of $\sim 2\,r_{\rm j}$, and secondary plasmoids are formed. In this phase of evolution, the fraction of initial energy released to the particles is small ($\epsilon_{\rm kin}/\epsilon_{\rm tot}(0)\sim 0.1$, top panel in Fig. 15), but the particles advected into the central X-point experience a dramatic episode of acceleration. As shown in the bottom panel of Fig. 15, the cutoff Lorentz factor $\gamma_{\rm max}$ of the particle spectrum presents a dramatic evolution, increasing up to $\gamma_{\rm max}/\gamma_{\rm th}\sim 10^{3}$ within a couple of dynamical times (here, $\gamma_{\rm th}$ is the initial “thermal” Lorentz factor).

This phase of extremely fast particle acceleration on a dynamical timescale is analogous to the relaxation of unstable “ABC” force-free structures discussed by Nalewajko et al. (2016); Lyutikov et al. (2016), and it constitutes the most promising scenario to explain the Crab flares. The particle acceleration efficiency and the hardness of the power-law slope depend on the mean magnetization of the configuration, in a similar fashion as in plane-parallel reconnection scenarios. The particle spectrum gets harder as the mean magnetization increases (Fig. 16); both the non-thermal particle fraction and the maximum particle energy increase with the magnetization (Lyutikov et al., 2016; Nalewajko et al., 2016).

Since the highest energy particles are first accelerated in the current layers by the parallel electric field, they do not radiate much when they are inside the sheet, because the curvature of their trajectory is small. Most of the radiation is produced when particles are ejected from the current layers — their trajectories start to bend significantly in the ambient magnetic field (Yuan et al., 2016). Such a separation of acceleration site and radiative loss site facilitates acceleration beyond the synchrotron radiation reaction limit, as required by the Crab flares. Fast variability of the observed photon flux can be produced when compact plasmoids that contain high-energy particles are ejected from the ends of the current layers and get destroyed. These give beamed radiation. An observer sees high intensity radiation when the beam happens to be aligned with the line of sight. Such peaks in emission are accompanied by an increase in the polarization degree and rapid change of polarization angle in the high-energy band (Yuan et al., 2016). The variability timescale is determined by the spatial extent of the emitting structure, e.g. the plasmoids, thus can be much shorter than the light crossing time of the region that collapses.

0.4 Conclusions

In this chapter, we have discussed the role of PIC simulations in unveiling the origin of the emitting particles in PWNe. After describing the basics of the PIC technique, we have summarized its implications for the quiescent and flaring emission of the Crab Nebula, as a prototype of PWNe. A consensus seems to be emerging that, in addition to the standard scenario of particle acceleration via the Fermi process at the termination shock of the pulsar wind, magnetic reconnection in the wind, at the termination shock and in the nebula plays a major role in powering the multi-wavelength emission signatures of PWNe.

Acknowledgements.

LS acknowledges support from DoE DE-SC0016542 and NASA Fermi NNX16AR75G. BC acknowledges support from CNES and Labex OSUG@2020 (ANR10 LABX56).