Dissipation of the striped pulsar wind and non-thermal particle acceleration: 3D PIC simulations

The formation of a large-scale current sheet is a generic feature of pulsar magnetospheres. If the magnetic axis is misaligned with the star rotation axis, the current sheet is an oscillatory structure filling an equatorial wedge determined by the inclination angle, known as the striped wind. Relativistic reconnection could lead to significant dissipation of magnetic energy and particle acceleration although the efficiency of this process is debated in this context. In this study, we aim at reconciling global models of pulsar wind dynamics and reconnection in the stripes within the same numerical framework, in order to shed new light on dissipation and particle acceleration in pulsar winds. To this end, we perform large three-dimensional particle-in-cell simulations of a split-monopole magnetosphere, from the stellar surface up to fifty light-cylinder radii away from the pulsar. Plasmoid-dominated reconnection efficiently fragments the current sheet into a dynamical network of interacting flux ropes separated by secondary current sheets which consume the field efficiently at all radii, even past the fast magnetosonic point. Our results suggest there is a universal dissipation radius solely determined by the reconnection rate in the sheet, lying well upstream the termination shock radius in isolated pair producing pulsars. The wind bulk Lorentz factor is much less relativistic than previously thought. In the comoving frame, the wind is composed of hot pairs trapped within flux ropes with a hard broad power-law spectrum, whose maximum energy is limited by the magnetization of the wind at launch. We conclude that the striped wind is most likely fully dissipated when it enters the pulsar wind nebula. The predicted wind particle spectrum after dissipation is reminiscent of the Crab Nebula radio-emitting electrons.


Introduction
Large-scale current sheets are generic features of planetary and stellar magnetospheres. Their formation can be externallydriven as in the Earth magnetotail shaped by the Solar wind, or internally-driven by the intrinsic magnetic activity of the star or if the magnetosphere is rapidly-rotating like in Jupiter. Perhaps the most extreme example of rotationally-driven current sheets in an astrophysical environment is found in the vicinity of pulsars. The short rotation period of the star (P ∼ 1 − 10 3 ms) combined with strong surface magnetic fields (B ∼ 10 9 − 10 12 G) lead to significant field line winding and opening beyond the light cylinder, a virtual cylindrical surface of radius R LC = cP/2π ∼ 50 − 50, 000 km, beyond which the corotation velocity becomes superluminal. A large-scale current sheet forms outside the light cylinder where both magnetic polarities of the star meet (Michel 1971;Coroniti 1990). The magnetic and current structures are supported by a plasma of relativistic electron-positron pairs self-generated near the stellar surface via pair production. This plasma flows along open field lines in the form of a radially expanding, relativistic magnetized wind, simply referred as the pulsar wind in the following (Rees & Gunn 1974;Kennel & Coroniti 1984).
If the magnetic axis is aligned with the rotation axis, the magnetosphere is axisymmetric and the current sheet is a flat disk lying in the equatorial plane. If the magnetic axis is inclined, the current sheet has the shape of an oscillatory structure of wavelength 2πR LC confined within an equatorial wedge of half-opening angle set by the magnetic inclination angle (Bogovalov 1999). A cut through this surface at a constant latitude gives rise to a succession of narrow stripes of currents separated by smooth wind regions ( Figure 1). For this reason, this structure is usually referred to as the "striped wind" (Coroniti 1990;Kirk et al. 2009). Away from this region, the wind is smooth and is well described by a rotating monopole-like configuration (Michel 1973).
One fundamental question refers to the fate of the stripes as the wind propagates outward, and so far this issue has led to contradictory conclusions (Coroniti 1990;Lyubarsky & Kirk 2001;Lyubarsky 2003;Kirk & Skjaeraasen 2003;Cerutti & Philippov 2017). It is generally accepted that relativistic reconnection occurs within the current layer (e.g., Kagan et al. 2015) leading to a transfer of magnetic energy into particle kinetic energy. The main uncertainty lies in the rate of dissipation and its feedback on the global dynamics of the wind. In a collisionless plasma, the current layer thickness is determined by the plasma kinetic scales, i.e., of order the plasma skin-depth and particle Larmor radius scales. In pulsars, this scale is microscopic such that the aspect ratio of the current sheet is very Article number, page 1 of 11 arXiv:2008.11462v1 [astro-ph.HE] 26 Aug 2020 A&A proofs: manuscript no. striped3d large meaning that the layer will most likely reconnect into the plasmoid-dominated regime, a regime where reconnection is fast (Uzdensky et al. 2010). Recent particle-in-cell (PIC) simulations of plane-parallel reconnection have confirmed the efficiency at dissipating the field and at accelerating particles of relativistic reconnection mediated by the plasmoid instability (e.g, Zenitani & Hoshino 2001;Cerutti et al. 2012;Sironi & Spitkovsky 2014;Werner et al. 2018). In a parallel effort, global PIC simulations of pulsar magnetospheres have shown the major role of reconnection at dissipating a sizeable fraction of the Poynting flux into high-energy particles and pulsed gamma-ray emission Philippov & Spitkovsky 2018;Kalapotharakos et al. 2018). These simulations were focused on the inner magnetospheric regions and restricted to a few light-cylinder radii only so that the large-scale evolution of dissipation was not probed.
In this study, we aim at reconciling global models of pulsar wind dynamics, reconnection and particle acceleration in the stripes within the same numerical framework, using large threedimensional (3D) PIC simulations. The latter are supplemented by a series of two-dimensional (2D) simulations restricted to the equatorial plane to explore the parameter space and the effect of numerical resolution. This work is the logical continuation of our previous effort in this direction (Cerutti & Philippov 2017). We begin first by introducing the numerical setup in Sect. 2. Simulation results are presented in Sect. 3 and discussed in Sect. 4 with an emphasis on dissipation and particle acceleration. Radiative signatures will not be discussed in this paper and will be left to a future study.

Methodology and setup
We use the relativistic electromagnetic PIC code ZELTRON (Cerutti et al. 2013;Cerutti & Werner 2019) in its full 3D spherical coordinates (r, θ, φ) version first introduced in Cerutti et al. (2016). The numerical grid is logarithmically spaced along the r-axis. This choice is well-suited for this problem where the plasma density and the field strength present a sharp gradient in the vicinity of the star, and at the same time allows us to probe large physical distances, the key objective in this work. The grid along the θ-direction follows a cos θ-spacing. This is a natural choice for a 3D spherical grid as it keeps the volume of the cell the same at a given radius, i.e., the grid is refined at the equator but it is coarser at the poles. This choice is also motivated by the fact that the pulsar wind power is concentrated within the equatorial regions (for a monopole it scales as ∝ sin 2 θ). The grid is uniformly spaced along the φ-direction.
The full numerical grid is composed of (2016 × 1024 × 512) cells along the r-, θand φ-directions respectively. The box extends from the stellar surface r min = r up to r max = 100 R LC where we fixed R LC /r = 3, θ = [0.03π, 0.97π] and φ = [0, 2π]. A damping layer absorbs all outgoing electromagnetic waves and particles to mimic an open boundary . We apply reflective boundary conditions along the θ-boundaries for the particles and axial symmetry to the fields (Holland 1983). As for the φ-direction, we apply standard periodic boundary conditions to the fields and the particles. The rotation axis of the star is aligned with the axis of the spherical domain, i.e. θ = 0. The magnetic axis is inclined at an angle χ with respect to the rotation axis and is rotating at the angular velocity of the star Ω. Following Cerutti & Philippov (2017), we chose a split-monopole magnetic configuration (Michel 1973;Bogovalov 1999) which is a good proxy for the asymptotic structure of the pulsar wind which is the main region of interest in this work. At t = 0, the magnetic field configuration is purely radial, where ±B is for the northern/southern hemisphere. At any instant t, the plane which separates both magnetic polarities at the stellar surface is given by the following condition sin θ sin χ cos (Ωt − φ) + cos χ cos θ = 0.
The solid rotation of field lines is enforced by applying the corotation electric field at the stellar surface at every time step. Assuming a perfectly conducting neutron star yields Fresh plasma exclusively composed of electron-positron pairs is uniformly injected at the stellar surface. The plasma is in corotation with the star and has a net radial velocity determined by the E × B drift velocity of the monopole solution. Normalized by the speed of light, the initial particle velocity components are (Cerutti & Beloborodov 2017) where R = r sin θ is the cylindrical radius. The injected plasma is neutral and has a multiplicity κ = n /n GJ = 10, where n GJ = ΩB /2πec is the fiducial plasma density (Goldreich & Julian 1969) and e the elementary electric charge. This prescription is a simple numerical recipe to fill efficiently the magnetosphere with plasma and therefore reach a quasi force-free configuration which is most appropriate to model active pulsars . In this work, we are not aiming at modeling pair production which is most likely happening in the inner magnetosphere but focus our numerical resources on the wind region instead, assuming that plenty of pairs are produced along all field lines. On average, one new pair is injected per cell at the neutron star surface every 8 time step. At the end of the simulations when most of the numerical box has been filled with plasma, there are about ∼ 10 10 particles which represents about 10 particles per cell on average. In addition to the Lorentz force, particles feel the radiationreaction force due to both curvature and synchrotron radiation (see Cerutti et al. 2016 for its implementation). The simulations time step is fixed at half the Courant-Friedrichs-Lewy time step, ∆t = 0.5∆t CFL . One pulsar period is 3.6 × 10 4 ∆t. All simulations were evolved for about 10 rotation periods, i.e., for about 3.6 × 10 5 ∆t. The fiducial plasma magnetization at the star surface is set at where m e is the electron rest mass. The shortest plasma time scale, ω −1 pe , as well as the shortest radiative cooling time scale, ω −1 c , are well resolved in all simulations, with ω pe ∆t ≈ 20, and Notes. 2D runs are in the rφ-plane for θ = 90 • . 'log' and 'uni' stands respectively for logarithmic and uniform grid spacing along the rdirection. r max is expressed in units of R LC and the obliquity angle χ is in degrees.
ω c ∆t ≈ 14. The smallest plasma scale, the local skin depth d e = c/ω pe , is also very well resolved everywhere in the simulation box with d e /∆r 10 − 30, where ∆r is the radial size of a cell. At the light cylinder, the local plasma skin depth is d e /R LC ≈ 3 × 10 −2 . We have performed 3D simulations with three different magnetic obliquity angles: χ = 30 • , 60 • and 85 • , all other parameters remaining identical. We ran another 3D run for χ = 60 • and r max = 50 R LC with a constant spacing along the radial direction to evaluate the role of the grid on magnetic dissipation. We have also performed a series of 2D runs limited to the equatorial plane as in Cerutti & Philippov (2017) to explore the parameter space and to investigate the role of numerical resolution on the dissipation of the striped wind. To this end, we performed a first series of 6 simulations with κ = 10 and σ = 250 using a logarithmic or a constant radial grid spacing with different numerical resolution both in r and φ. In a second set, we aim at assessing the role of the wind magnetization σ = 125, 250, 750 while keeping κ = 10. In the last set, we explore the role of the plasma multiplicity while keeping the stellar magnetic field the same, with κ = 0.2, 0.6, 2, 6, 20, 60. Table 1 draws the full listing of all runs used in this work.

Plasma structures
Shortly after the onset of the simulations, the plasma blown from the star establishes a force-free split-monopole-like con- figuration as it propagates out. It is characterized by a toroidaldominated magnetic structure whose polarity reverses across the current sheet. The shape of this sheet is consistent with the prediction of Bogovalov (1999): it is a 3D spherical Archimedean spiral of wavelength 2πR LC defined within the spherical wedge π/2 − χ ≤ θ ≤ π/2 + χ (Figures 1, 2). This region, which includes both the sheet and the plasma in between, is the striped wind (Michel 1971;Coroniti 1990;Kirk et al. 2009). Away from the stripes, the wind nearly resembles a single rotating magnetic monopole (Michel 1973).
The current sheet fragments soon after its formation near the light cylinder. It is unstable to the relativistic tearing instability (Zelenyi & Krasnoselskikh 1979;Pétri & Kirk 2007) which mediates fast magnetic reconnection. It results in the formation of a dynamical chain of plasma overdensities trapped in magnetic loops, or plasmoids, separated by secondary current sheets where the field reconnects. These features are clearly visible on the zoomed-in view of the plasma density in both the toroidal and poloidal planes in Figure 2. In full 3D, these structures form a network of interconnected flux ropes reminiscent of plane parallel reconnection simulations (Kagan et al. 2013;Cerutti et al. 2014;Werner & Uzdensky 2017). Secondary flux ropes are being produced within secondary current sheets which then merge with others to form bigger structures. The dynamical nature of reconnection is most pronounced within 10 R LC and is quenched by the expansion of the wind at larger radii. Although severely fragmented, the striped wind structure retains its global coherent structure. Hence, the 2D picture drawn in Cerutti & Philippov (2017) qualitatively holds in full 3D.
We observe that the sheet is also mildly kink unstable (Zenitani & Hoshino 2005;Cerutti et al. 2014). It appears as wiggles most visible in the zoomed-in view of the poloidal plane in Figure 2, but these distortions do not lead to the complete disruption of the sheet. The kink is more effective at low magnetic obliquities Cerutti et al. 2016) and can lead Fig. 2. 2D slices of the plasma density r 2 (n/n GJ ) for χ = 60 • . Top: rθ-plane containing the magnetic axis at the phase φ = Ωt. Bottom: rφ-plane at the equator (θ = 90 • ). The radius is expressed in units of R LC . The right panels show zoomed-in views of the regions delimited by the white boxes drawn on the left panels. The striped wind region is contained within π/2 − χ ≤ θ ≤ π/2 + χ (white dashed lines).
to a significant latitudinal spreading of the striped wind region, well outside of its natural boundaries. At χ = 30 • , the kink instability may be responsible for a widening of the stripes of about 10 • , whereas there are no noticeable deviations for χ = 60 • and above ( Figure 2, top panel). The plasma density is distributed in a highly inhomogeneous and anisotropic manner. The unstriped wind is composed of a low-density uniform plasma, of multiplicity κ = n/n GJ of order unity, except close to the axis where numerical plasma fluctuations are artificially enhanced by the spherical grid. In the striped zone, there is a strong plasma density contrast between, in order of increasing density, the wind (κ ∼ 2-3), secondary current layers (κ ∼ 4-5) and the far more denser flux ropes (κ ∼ 10-10 3 ). The wind zone itself between two stripes is inhomogeneous, in contrast with the unstriped region. There is a clear plasma depletion on the leading edge of the spiral which was already reported in previous studies Cerutti & Philippov 2017;Philippov & Spitkovsky 2018). In this sense, reconnection proceeds in a highly asymmetric way.

Poynting flux and dissipation
We now turn our attention to the central question of magnetic dissipation. In a dissipationless steady state split-monopole magnetosphere, the outgoing Poynting flux integrated over a spherical radius is conserved in virtue of the Poynting flux theorem. The predicted value is (8) Figure 3 shows the radial and latitudinal dependence of the Poynting flux, L, for all 3D simulations. The first element to notice is that the numerical values are closer to L ≈ L 0 /5 at the star surface. This discrepancy is explained by the fact that the analytical split-monopole solution used in Eq. (8)   cylinder, and thus do not contribute to the outflowing Poynting flux. Only polar field lines remain open so that the magnetospheric structure simulated here qualitatively resembles a forcefree dipole inside the light cylinder. Nonetheless, the Poynting flux does not present a strong dependence with magnetic obliquity as expected from the split-monopole solution (Bogovalov 1999 so that Eq. (8) should be changed into the more general expression (e.g., Tchekhovskoy et al. 2016) Thus, physically L/L 0 ≈ 0.2 is a measure of the fraction of the solid angle squared filled by open field lines in the inclined splitmonopole simulations.
The presence of reconnection in the current sheet leads to significant dissipation which translates into a decay of the radial Poynting flux. Here, dissipation is by no means spurious numerical dissipation but it has a real physical origin: the reservoir of Poynting flux provided by the star is gradually consumed by reconnection which converts magnetic free energy into particle kinetic energy and radiation, such that the total energy in the system is conserved to a good accuracy as shown in Figure 4. Irregularities in the total power curve reflect the intermittent nature of reconnection. The radial profile of the Poynting flux first drops by about 20% within a few light-cylinder radii after launching, a consistent number with past studies focusing on the magnetosphere and inner wind zone (Parfrey et al. 2012;Philippov & Spitkovsky 2014;Cerutti et al. 2015;Belyaev 2015). Dissipation continues further at a much slower but steady rate beyond r 10 R LC until it reaches about 40-50% of the total initial flux at r = 50 R LC (Figure 3). This rate seems independent of the magnetic obliquity. It does not seem to be sensitive on the choice of the grid spacing either: after a brief overshoot at low radii, likely due to the low numerical resolution in comparison with the log-spacing run, the constant r-spacing simulation asymptotically converges towards the same amount and rate of dissipation (see the solid and the dashed red lines in Figure 3).
To strengthen this point further, 2D simulations restricted to the equatorial plane with different numerical resolutions, box sizes and grid spacing present the same evolution of Poynting flux and amount of dissipation (top panel in Figure 5). Figure 3 also shows the latitudinal dependence of the Poynting flux and of its dissipation. As expected, the Poynting flux is preferentially distributed within the equatorial regions and the θ-profiles closely resemble the split-monopole solution, i.e., dL/dΩ ∝ sin 2 θ. We note that the profiles are slightly sharper than predicted although not as much as reported in Tchekhovskoy et al. (2013) where dL/dΩ ∝ sin 4 θ for an initially dipolar magnetic field configuration. The transition between the striped and the unstriped wind regions is smooth except at χ = 30 • where small dips are visible. An interesting feature is that the same fraction of power is dissipated at all latitudes within the striped-wind region. This fraction reaches about 40-50% depending on the obliquity, and vanishes within the dissipationless unstriped-wind region (see lower-right panel in Figure 3) 1 . Incidentally, the kink-induced widening of the sheet at χ = 30 • described in Sect. 3.1 results in a broader angular dissipation rate profile, well outside the boundaries 60 • ≤ θ ≤ 120 • .
In Figure 5 (middle and bottom panels), we report on the dependence of dissipation with the plasma magnetization, σ , and multiplicity, κ , based on our large set of 2D simulations. While there is no noticeable dependence with magnetization (at least as long as σ 1), dissipation monotonically increases with increasing plasma multiplicity. Low-multiplicity solutions (κ < 1) present 10-20% dissipation rate at r/R LC = 50. In contrast, high-multiplicity solution (κ 1) show up to 70-80% dissipation without any sign of saturation with radius, indicating that the striped wind structure most likely disappear far before the pulsar wind terminates in isolated systems (Sect. 4). Putting ev-  Table 1 for grid parameters). The 3D run with χ = 85 • is reported here for comparison with 2D runs (black solid line). Middle: Dependence of dissipation with the plasma magnetization for κ = 20. Bottom: Dependence of dissipation with the plasma multiplicity for σ = 250. The dashed lines are bestfit models to the analytical solution proposed in Sect. 4,Eq. (20). In all panels, the Poynting flux is normalized to its light-cylinder value, L/L LC . erything together, simulations reported in this work suggest that reconnection proceeds efficiently and homogeneously within the striped wind, regardless of magnetic inclination, numerical resolution, grid spacing and plasma magnetization, but there is a clear distinction between charge-starved (slow dissipation) and high-multiplicity winds (fast dissipation).

Wind kinematics and particle spectra
In this section, we exploit the particle data to reconstruct the bulk motion of the wind and the particle spectrum in the context of magnetic dissipation and particle acceleration. Figure 6 focuses on the radial and latitudinal dependence of the wind bulk Lorentz factor Γ. Here, all quantities are averaged over azimuth so that the wind can be referred as a single homogeneous entity. Within r/R LC 4, the wind accelerates quasi-linearly with cylindrical radius as expected from the monopole prediction, i.e., Γ = (1 + r 2 sin 2 θ/R 2 LC ) 1/2 if the plasma follows the E × B drift velocity (see Eqs. 4-6). Then, the wind quickly reaches the fast is Michel's magnetization parameter (Kirk et al. 2009). With σ ≈ 15 at the fast point, the wind remains highly magnetized. Past this point, the wind acceleration significantly slows down and saturates at Γ ∞ ≈ 11 for r/R LC 30 within the equatorial regions. This result is not sensitive to the magnetic obliquity angle and confirms theoretical expectations (Tomimatsu 1994; Beskin et al. 1998), and our previous findings in 2D although we obtain here significantly larger asymptotic values in full 3D (Γ 2D ∞ ≈ 6 reported in Cerutti & Philippov 2017). The latitudinal dependence of the bulk Lorentz factor is consistent with the monopole solution in the unstriped region with Γ(θ) ∝ sin θ. At large radii in the striped region, the profile is much flatter and significantly departs from the ideal solution (bottom panel in Figure 6). These deviations are most likely due to dissipation and particle acceleration within the current sheets whose effects are not negligible at large radii.
The wind acceleration is an ideal process solely governed by the E × B drift velocity. It should not be confused with nonthermal particle acceleration due to reconnection or other nonideal dissipative processes which comes on top of this ideal process. Figure 7 shows the radial and latitudinal evolution of the individual particle Lorentz factor spectrum, (γ/N)dN/dγ. The r-dependence highlights the formation of a broad non-thermal particle spectrum as magnetic energy is consumed and transferred to the particles via reconnection. The saturated-looking state of the particle spectrum at large radii hides a strong latitude  Fig. 7. Top panel: Normalized particle energy spectrum, (γ/N)dN/dγ, within the spherical shells of radius r and r + R LC as function of radius (color-coded) for χ = 60 • . Middle panel: Latitudinal dependence (colorcoded) of the particle energy spectrum averaged between 40 ≤ r/R LC ≤ 50 for χ = 60 • . Bottom panel: Total particle spectra beyond r = 40 R LC and 80 • ≤ θ ≤ 100 • for χ = 30 • , 60 • and 85 • . dependence shown in the middle panel of Figure 7. The spectrum is narrow and peaks at low energies near the rotation axis ( γ ∼ 3). The mean energy shifts to higher energy without any significant spectral broadening in the unstriped zone. This evolution follows the ideal bulk acceleration of the wind which is not accompanied by non-thermal acceleration. In contrast, the particle spectrum dramatically widens inside the striped region where non-thermal particle acceleration pushes particles up to γ Γ ∞ . This evolution is not very sensitive to the magnetic obliquity.
The asymptotic spectrum in the striped wind can be decomposed into two parts: a thermal bath peaking at γ ≈ Γ ∞ and a hard power-law spectrum dN/dγ ∝ γ −1 cutting at µ LC To uncover the origin of these two components, it is useful to look at the φresolved particle spectrum (see Figure 8). This analysis clearly shows that the low-energy particles are associated with the wind region located far upstream the current layers where the ideal monopole wind acceleration mechanism applies. A closer look reveals that the wind itself is composed of two distinct parts, of slightly different energies, a low-energy component at the leading edge of the spiral sheet corresponding to the low-density region with γ ≈ 8 and a high-energy component at the trailing edge of the spiral with γ ≈ 15. There is a sharp spatial segregation about half way in between two current sheets (around φ = 0 • and φ = 180 • in Figure 8) which marks the zone of influence Article number, page 7 of 11 A&A proofs: manuscript no. striped3d The locations of the current sheets and wind regions are labelled "Layer 1", "Layer 2" and "wind" respectively. Middle: Corresponding φ-resolved particle energy spectrum. Bottom: Decomposition of the total particle spectrum in the striped region (black dashed line) into a "wind" component (solid red line) and a "Layers" component (solid blue line). The vertical dotted lines show γ = Γ ∞ ≈ 11 and γ = σ LC ≈ 50. of each reconnecting layers on the plasma in the wind. As for the hard power-law component, it is unambiguously associated with non-thermal particle acceleration within the reconnection layers. The spectral index of ≈ −1 is consistent with high-σ reconnection simulation studies (Zenitani & Hoshino 2001;Sironi & Spitkovsky 2014;Guo et al. 2015;Werner et al. 2016). Therefore, to a first order, the asymptotic particle spectrum within the striped wind can be approximatively described as The asymptotic spectrum does not seem to collapse into a narrower distribution at large radii, in contrast to our previous 2D runs suggesting that 3D effects may be important in producing a power-law spectrum.

A toy model for dissipation
We propose a simple analytical model for the evolution of magnetic dissipation inspired by simulations in an attempt to extrapolate our results to realistic pulsar winds. In steady state, the radial variations of the Poynting flux is solely governed by Joule dissipation between two spherical shells that we will choose here to be of radius r = R LC and an arbitrary radius r > R LC , such that To a very good accuracy, we have To make further progress, we need to have a closer look at the distribution of currents and fields in the vicinity of current sheets where Joule dissipation is localized. Figure 9 shows the radial profile of B φ , J θ and the product of the two for κ = 0.2, 2, and 20 at a pulsar phase φ chosen such that it crosses only currents sheets and avoids magnetic islands. The spatial distribution of J · E within islands is dipolar such that there is no net contribution integrated over their volume, in contrast with secondary current sheets. The magnitude of the current carried by the particles depends on the plasma multiplicity and is localized in the form of thin sheets as expected. At low multiplicities (κ < 1), the conduction current is small, too small in fact to explain the magnetic reversal and it is not localized at the magnetic nulls. In this regime, the current is mostly carried by the displacement current (∂E/∂t, only possible for an oblique rotator), which is a logical consequence of charge starvation in the wind. With increasing multiplicities, the magnetic profile changes from a nearly sinusoidal shape at low-κ , to an asymmetric square shape at high-κ with sharp gradients where the field reverses and where a strong conduction current is localized. Thus, the electric field in Joule's term may be understood as the reconnection electric field, E rec , which represents a fraction of the upstream magnetic field, B up φ , such that where β rec is the dimensionless reconnection rate (Lyubarskii 1996;Uzdensky & Spitkovsky 2014). Another important observation we can make from Figure 9 is that the magnetic field strength which contributes to dissipation is always on the trailing edge of the spiral because of the asymmetric nature of reconnection in the striped wind, the most pronounced effect being visible at high multiplicities. This conclusion is also compatible with the asymmetry in the phase-resolved particle spectrum reported in Sect. 3.3. We find that the upstream magnetic field strength at the trailing edge of the sheet, B up φ , is of order the ideal split monopole field at all radii, even though the striped-averaged field strength decreases with radius due to dissipation, i.e., such that the sheet is always fed with fresh, unreconnected magnetic field. The electric current is given by Ampère's law. Assuming perfect symmetry on both sides of the current layer for the sake of simplicity, we have (e.g., Cerutti & Philippov 2017) 4π where δ is the layer thickness. Putting everything together, Joule's term within a single current sheet can be approximately estimated as The integral over φ at constant radius and latitude vanishes everywhere, except at two phases where the layers are located. The angular width of each layer can be estimated as ∆φ ∼ δ/r, such that Eq. (13) becomes i.e., dissipation is independent of the width of the current layer. The integrals over θ and r lead to the final result,  Fig. 10. Evolution of the dissipation rate β rec (blue line, dots and axis) and the extrapolated full dissipation radius of the stripes R diss /R LC = exp(β −1 rec ) (red line, dots and axis) with the typical plasma multiplicity measured at secondary current sheets, κ X . The estimated radius of the Crab pulsar wind termination shock is shown by the horizontal dotted line (R TS /R LC ∼ 10 9 ) for comparison.
where we absorbed the constant 2/π of order unity into the reconnection rate, i.e., β rec ← 2β rec /π. Applying this model to simulation data provides a good description of the Poynting flux decay and a direct measure of the rate of dissipation. Figure 10 shows the evolution of the dissipation rate with the plasma multiplicity. It is determined from the best-fit model to the dissipation curves in Figure 5 (bottom panel), assuming the evolution given in Eq. (20). As expected, at low multiplicities dissipation is low because the current is mostly carried away by the displacement current. The contribution from the conduction current increases with plasma supply until it reaches an approximate saturation at high-multiplicities where β rec ∼ 0.1 -0.2. One should keep in mind that the rate measured from the simulations includes additionnal physical effects neglected in the above toy model which are of order unity, such as asymmetries, deviations from the split monopole solutions, or the filling factor of plasmoids which do not contribute to dissipation. These effects may account for the slow increase of the dissipation rate at high multiplicities.
Another effect not captured by this simple model is dissipation due to the formation of vacuum gaps in the low-multiplicity solutions, as reported in Cerutti et al. (2015) for an aligned rotator. In the 2D equatorial plane simulations used here, Ω · B = 0 and therefore no gap can form in this special configuration. We would expect an additional source of dissipation at low, but finite multiplicities in a full 3D setup, which we did not explore in this work. In this sense, the model gives a lower limit of dissipation in this regime. This being said, the real solution must smoothly connect to the fully dissipationless vacuum solution (Deutsch 1955), meaning that dissipation should necessarily cease as we asymptotically approach the vacuum regime, as reported here.

Implications
If the model presented above is a fair description of real astrophysical pulsar winds, it has important astrophysical implications. The first observation to make is that dissipation depends only on the reconnection rate. Studies of local plane-parallel reconnection consistently find a high reconnection rate, β rec ∼ 0.1 -0.2 (see, e.g., Kagan et al. 2015 and references therein), similar to what is reported here. This rate does not seem to depend on anything as long as the layer is thin, meaning that it is of the order the plasma kinetic scales, δ ∼ d e . In particular, it should not depend on the plasma multiplicity as long as κ 1 as reported in Figure 10, and thus, the radius at which the striped wind has fully dissipated should be similar for all pulsars producing pairs, i.e., like gamma-ray pulsars. In this sense, our results suggest that there is a universal dissipation radius for all pulsar winds loaded with pairs. Using Eq. (20) and assuming a naive extrapolation of the model to any radius, we predict a complete decay of the Poynting flux at a radius for β rec = 0.1 -0.2. Although this radius is quite sensitive to the exact value of the reconnection rate, simulations show that β rec is high meaning that the striped wind will most certainly decay entirely far before reaching its termination (Coroniti 1990;Cerutti & Philippov 2017), unless it is truncated at short distances by an accretion disk or a companion star wind. At this stage, it is important to emphasize that the expansion of the current layer plays no role at dissipating the field, as long as the layer thickness remains small compared with the stripe half-wavelength, δ/πR LC 1 (Lyubarsky & Kirk 2001;Kirk & Skjaeraasen 2003;Zrake & Arons 2017), which is the case in all of the simulations reported here. In contrast to Eq. (21), this condition explicitly depends on the plasma multiplicity: a denser sheet is also a thinner one, thus the condition for which two consecutive sheets would overlap is pushed further away as the multiplicity increases (Cerutti & Philippov 2017). Therefore, R diss should be seen as an upper limit for the radius of full dissipation. It is also worth noting that we do not find any evidence for a significant bulk acceleration of the wind due to magnetic dissipation as anticipated by Lyubarsky & Kirk (2001) and Kirk & Skjaeraasen (2003). Here, reconnection proceeds at a similar rate even past the fast point where the bulk Lorentz factor remains constant ( Figure 6). As discussed further below, dissipation does not perform work on the wind as a bulk but rather benefits to disorganized energetic pairs trapped within plasmoids.
The wind kinematics and the shape of the particle spectrum reported here also have important astrophysical consequences. The commonly accepted picture is that the wind is composed of cold, nearly monoenergetic pairs travelling at ultra-relativistic velocities such that Γ ∼ 10 4 − 10 6 (Rees & Gunn 1974;Wilson & Rees 1978;Kennel & Coroniti 1984). While most of the predicted wind properties are recovered here, we find that the bulk Lorentz factor may not be as relativistic as previously thought. A good proxy is given by the wind Lorentz factor at the fast magnetosonic point, Γ ∞ ∼ µ 1/3 M . The magnetization at the light cylinder is poorly constrained mostly because of the uncertainty on the plasma multiplicity, but it could of order where P 100 = P/100 ms is the pulsar period, B 5 = B LC /10 5 G and κ 4 = κ/10 4 , for a typical young gamma-ray pulsar. Hence, we expect Γ ∞ ∼ µ 1/3 M ∼ 100 at most, i.e., more similar to what is usually inferred in gamma-ray burst jets 2 . This gives a good estimate of the bulk particle Lorentz factor throughout the pulsar wind. It is also a fair estimate of the individual particle Lorentz factor located within ideal regions, i.e., the unstriped polar regions and the inter-stripe medium where the particle spectrum remains narrow and cold. Within the current layers, relativistic reconnection leads to the formation of a broad power law of index p ∼ −1, with a low-energy cut off at γ ≈ Γ ∞ ∼ 10 2 , and a high-energy cut off at γ ≈ µ LC M ∼ 10 6 . The fundamental difference with the standard picture is that the pairs accelerated in the sheet are not cold with γ = Γ, but instead the particles remain hot and trapped within the flux ropes in the wind frame. Thus, we expect the pulsar wind to inject ultrarelativistic pairs with a hard spectrum into the nebula and not just at a single energy scale. Interestingly, applying this model to the Crab Nebula with the typical magnetization scale quoted above µ M ∼ 10 6 could provide a natural explanation for the mysterious radio-electron component responsible for the hard, low-energy emission from the nebula. In the spectral modeling by Meyer et al. (2010), this population is well fitted by a single power law of index −1.6 cutting off at γ min ∼ 20 and γ max ≈ 2 × 10 5 which fits well within our results. Although our spectrum is slightly too hard, a larger and therefore more realistic separation of scales between the layer thickness and the light-cylinder scale could lead to a significant spectral steepening (Petropoulou & Sironi 2018). Pairs injected into the shock could be further accelerated into the nebula by some other mechanisms to form the softer X-ray to gamma-ray electron component which would naturally result in a smooth transition between both components.

Summary
We have performed large 3D PIC simulations of pulsar winds with a focus on magnetic dissipation and particle acceleration within the striped region. The global structure of the striped wind is consistent with the split-monopole prediction (Michel 1973;Bogovalov 1999). The current sheet is prone to the plasmoid instability shortly after its launching at the light cylinder. This instability leads to an efficient fragmentation of the sheet into a network of interconnected flux ropes separated by secondary thin current layers where the field reconnects, a structure reminiscent of 3D plane-parallel reconnection studies. This chain is highly dynamical, flux ropes form and merge to form bigger structures, which may result in bright short bursts of radio emission (Philippov et al. 2019;Lyubarsky 2019) aligned with the incoherent gamma-ray pulsed emission from the sheet . The sheet is also kink unstable which leads to a significant widening of the striped wind region at low magnetic obliquity. Reconnection in this environment proceeds in a highly asymmetric way, in a qualitatively similar manner as in the dayside of the Earth magnetopause. More efforts are needed to better understand asymmetric reconnection in the relativistic regime using local studies.
Relativistic reconnection gradually consumes the oscillatory component of the toroidal magnetic field. This feature is robust against numerical resolution, grid spacing and plasma magnetization. The Poynting flux monotonically decreases with the distance to the light cylinder, reaching up to about 40% dissipation at the outer parts of the box, i.e., r = 50 R LC . The dissipation rate weakly depends on latitude within the striped region and on the magnetic obliquity angle. Based on a large set of 2D simulations restricted to the equatorial plane, we can establish that the fate of the striped wind is not sealed by the sheet width as previously thought, but it is rather controlled by the dimensionless reconnection rate β rec ≈ 0.2 which lies well within reported values in plane parallel reconnection studies where β rec ∼ 0.1 − 0.2. This rate is known to weakly depend on the system size and plasma magnetization in the ultra-relativistic regime (σ 1, Werner et al. 2018), meaning here that dissipation should not depend on other parameters such as the plasma multiplicity as long as there is a large supply of pairs (κ 1). Therefore, we propose there is a universal dissipation radius valid in all pair producing pulsars of order R diss /R LC = exp β −1 rec ∼ 10 2 − 10 4 , meaning that the stripes should disappear far before reaching the wind termination shock radius in isolated systems like the Crab pulsar where R TS /R LC ∼ 10 9 .
The wind bulk Lorentz is not as relativistic as previously imagined. After crossing the fast magnetosonic point, the wind speed quickly saturates to Γ ∞ ≈ µ 1/3 M and thus does probably not exceed Γ ∞ 100 which is closer to gamma-ray burst jet Lorentz factor than the standard Crab pulsar wind model where Γ ∼ 10 3 − 10 6 (Rees & Gunn 1974;Wilson & Rees 1978;Kennel & Coroniti 1984). On top of the wind acceleration driven by ideal processes, non-thermal particle acceleration proceeds in the striped region driven by reconnection. The particle spectrum is consistent with high-σ reconnection simulations, meaning a hard power-law distribution of index p ∼ −1 between γ min ∼ µ 1/3 M and γ max ∼ µ M . Scaled up to realistic Crab-like parameters yields the pulsar wind would be composed of ultra-relativistic pairs distributed as dN/dγ ∝ γ −1 between 10 2 γ 10 5 . Injected at the shock front, the wind particles could then naturally explain the mysterious hard radio component in the Crab Nebula (Meyer et al. 2010).