EDP Sciences
Highlight
Free Access
Issue
A&A
Volume 607, November 2017
Article Number A134
Number of page(s) 13
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201731680
Published online 28 November 2017

© ESO, 2017

1. Introduction

It is commonly accepted that pulsars launch an ultrarelativistic, ultramagnetized wind of electron-positron pairs beyond the light-cylinder radius where co-rotation with the star is impossible (Kirk et al. 2009; Arons 2012; Pétri 2016; Cerutti & Beloborodov 2017). Observations of pulsar wind nebulae indicate that the plasma is purely non-thermal and weakly magnetized downstream of the wind termination shock (Rees & Gunn 1974; Kennel & Coroniti 1984; Begelman & Li 1992). Therefore, there must be an efficient transfer of the magnetic energy into particle kinetic energy somewhere between the pulsar magnetosphere and the nebula. This puzzle is well-known in the pulsar community as the “sigma problem”, but this is in fact a more generic problem of all relativistic magnetized outflows such as relativistic jets found in active galactic nuclei (AGN) and gamma-ray bursts. It is still unclear how and where magnetic dissipation happens in pulsars, but these objects are ideal testbeds for understanding this physical process.

The wind current sheet that forms between the two predominantly toroidal magnetic polarities is often suspected to be the main place for magnetic dissipation. If the magnetic axis is aligned with the pulsar rotation axis, the sheet is flat and fills the equatorial plane. If the magnetic axis is inclined at an angle χ ≠ 0, the sheet is oscillating and fills an equatorial wedge contained within π/ 2 ± χ and has the shape of a ballerina’s skirt of wavelength set by the light-cylinder radius (2πRLC) that we refer to herein as the “striped wind” (Coroniti 1990; Bogovalov 1999). Using a laminar model of magnetic dissipation, Coroniti (1990) and Michel (1994) found complete annihilation of the striped wind far before it reaches the termination shock. Lyubarsky & Kirk (2001) revisited this scenario and found that the wind accelerates as dissipation proceeds so that relativistic time dilation effects would not allow the wind to dissipate significantly, unless the particle density is higher than usually expected (Kirk & Skjæraasen 2003). This conclusion lead Lyubarsky (2003) to propose that the stripes survive up to the termination shock where they annihilate, leading to efficient plasma heating and non-thermal particle acceleration (Pétri & Lyubarsky 2007; Lyubarsky & Liverts 2008; Sironi & Spitkovsky 2011).

Regardless of where dissipation happens, the main physical mechanism at the origin of magnetic dissipation is reconnection. In the context of pulsars where the magnetic energy exceeds the rest mass energy of the plasma (high magnetization, σ ≫ 1), reconnection proceeds in the relativistic regime (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Lyubarsky 2005), that is, the plasma entering the reconnecting regions becomes necessarily relativistic. Recent ab-initio particle-in-cell (PIC) simulations have shown that relativistic reconnection is fast and efficient at dissipating the magnetic energy for both pair plasmas and electron-ion plasmas, and at producing broad hard particle spectra (see Kagan et al. 2015, for a review and references therein). These studies also revealed that fast reconnection is driven by the relativistic tearing instability (Zelenyi & Krasnoselskikh 1979) which fragments the current sheet into a hierarchical chain of magnetic islands separated by short current sheets (Uzdensky et al. 2010). In the close environment of pulsars, radiative losses are extreme (mostly curvature and synchrotron radiation) and effectively decrease the thickness of the layer, thus leading to an enhancement of the reconnection rate (Uzdensky & McKinney 2011; Uzdensky & Spitkovsky 2014). Pair creation at the base of the striped wind via photon-photon annihilation is another peculiarity of reconnection in pulsars (Lyubarskii 1996). The goal of this paper is to revisit the issue of magnetic dissipation and particle acceleration in the striped wind in light of the most recent advances in our understanding of relativistic reconnection.

Global PIC simulations of pulsar magnetospheres have already identified the important role of relativistic reconnection in the energy transfer between the fields and the particles (Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Philippov et al. 2015b; Cerutti et al. 2015; Belyaev 2015a; Philippov et al. 2015a), eventually leading to energetic pulsed emission (Cerutti et al. 2016; Philippov & Spitkovsky 2017). Simulations show efficient dissipation close to the light cylinder within 1−2 RLC, from 15−20% of the pulsar spindown power for an aligned rotator, down to a few percent for an orthogonal rotator (Philippov et al. 2015b). But these studies focus on the magnetosphere and hence are limited to a few RLC in radial extent, so the question of whether dissipation occurs at larger radii further in the wind but before the termination shock remains unanswered.

Here, we present the results of large box-size PIC simulations aimed at capturing the large-scale dynamics of the striped wind self-consistently. We pay little attention to the magnetosphere in this work, even though it is present in the simulation (see previous works). To reach large radial extensions, we perform two-dimensional (2D) simulations of the striped wind in the equatorial plane, starting from a split-monopole configuration. We begin with a detailed description of this numerical setup in Sect. 2 and present our results in Sect. 3, with an emphasis on the dynamics of reconnection along the striped wind. We also investigate how particle acceleration proceeds within the current sheet and we model the high-energy radiative signatures. Finally, we discuss the implications of these results in Sect. 4.

2. Numerical approach and setup

thumbnail Fig. 1

Numerical setup of a two-dimensional (2D) PIC simulation of the striped wind in the equator of an inclined split-monopole (Bogovalov 1999). Initially, magnetic field lines are purely radial and their polarity change at φ = Ωt ± π/ 2 (dotted line). The box extends from the stellar surface rmin = r (black disk) where e± pairs are created and where the fields are in co-rotation with the star, up to an absorbing layer between rabs and rmax where particles and fields vanish. The light-cylinder radius RLC = c/ Ω is shown by the dashed circle.

Open with DEXTER

We use the relativistic electromagnetic PIC code ZELTRON (Cerutti et al. 2013) with spherical coordinates (r,θ,φ) as described in Cerutti et al. (2015, 2016). In order to reduce numerical cost and hence capture the large-scale evolution of the striped wind, we restrict the spatial grid to the equatorial plane (θ = π/ 2). The numerical domain is therefore a spherical wedge in the -plane and should not be confused with polar coordinates1. Time-dependent Maxwell’s equations are solved with the usual second-order accurate Yee-algorithm (Yee 1966) directly on the spherical mesh using cell-integrated formulae while particle trajectories are solved in Cartesian coordinates using the Boris push (Birdsall & Langdon 1991). Particle positions and velocities are then remapped in spherical coordinates at every time step for current deposition. Particle positions are confined within the -plane but all three components of the 4-velocity vector are retained. In addition to the Lorentz force, particles feel the radiation-reaction force due to the emission of curvature and synchrotron radiation (see Cerutti et al. 2016, for its implementation), the two main radiative processes within the close environment of pulsars.

The computational domain is a full disk (φ ∈ [0,2π]) of inner radius rmin and outer radius rmax (see Fig. 1 for an overview of the setup). We place the inner boundary at the surface of the star rmin = r and the outer boundary far away from the light-cylinder radius, rmax = 170 RLC, where RLC = 3r in our fiducial simulations. We have also investigated slower rotators with RLC = 6r and 9r. The outer boundary of the box is coated with a thick layer that smoothly absorbs all particles and electromagnetic waves leaving the box. This open boundary can be achieved by adding finite electric and magnetic conductivities within the absorbing layer (Cerutti et al. 2015; Belyaev 2015b). It implicitly assumes that no information will be able to return inward. We see below that this is a very good assumption in this context because the wind crosses the fast magnetosonic point well before it reaches the outer boundary. The full domain is composed of (4096 × 4096) cells with a logarithmic spacing in radius and a constant spacing in φ. This configuration is well adapted for modeling the striped wind because the fields and the density decrease with radius (as 1 /r and 1 /r2, respectively, for rRLC) and hence the particle Larmor radius (ρL) and skin-depth (de) increase with radius. It also allows to keep the cell aspect ratio constant with radius (Δr/rΔθ = const.).

We assume that the neutron star is a perfect conductor rotating at the constant angular velocity Ω = (c/RLC)ez. Magnetic field lines are frozen into the surface so that their solid rotation with the star induce the co-rotation electric field (1)at the inner boundary. Initially, the box is in vacuum (no particles) with a static magnetic field configuration. We first tried to create a striped wind from an inclined rotating dipole as it is usually done in the literature, but we found that this choice is not appropriate in this numerical setup because the assumption of a zero-gradient along the θ-direction is not valid for a dipole (∂B/∂θ ≠ 0): we found unphysical cyclic changes of the magnetic structure, even in vacuum. Instead, we chose a split-monopole configuration (Michel 1973; Bogovalov 1999), that is, purely radial field lines whose polarity reverses across the line of equation φ = Ωt ± π/ 2 (Fig. 1). The split-monopole has the advantage of depending only on radius (∂B/∂θ = 0) meaning that this solution does not exhibit the pathologies that we observed with the dipole. This configuration is also well motivated physically because it is a good approximation of the asymptotic structure of the pulsar wind, which is the main focus of this study.

As soon as the simulation begins and at every timestep, electron-positron pairs are uniformly created at the surface of the star to model the polar-cap discharge (Timokhin & Arons 2013; Chen & Beloborodov 2014; Philippov et al. 2015a), which was shown by Philippov et al. (2015b) to be more efficient for inclined rotators and hence most relevant to this study. Pair creation via photon-photon annihilation away from the star is neglected here. The injected plasma is neutral (one pair of macroparticles is injected per cell and per timestep) and has a high multiplicity, κn/nGJ, where nGJΩ·B/ 2πec is the fiducial Goldreich-Julian density (Goldreich & Julian 1969). In our set of simulations, κ = 5, 10, and 20. The high plasma multiplicity allows to reach the quasi force-free regime in the magnetosphere and wind (i.e., ρE + J × B/c = 0 with ρ,J the charge and current densities), except in the current sheet where reconnection operates. It is worth noting that the plasma will remain neutral everywhere in the box (ρ ≈ 0) because both species play a symmetric role by construction of the numerical setup used here. This is most appropriate to describe an orthogonal rotator. In the more general case where Ω·B ≠ 0, the plasma in the magnetosphere and the current sheet is charged.

The magnetization of the plasma at the surface of the star is very high, (2)where me is the rest mass of the electron. In this work, σ varies from 250 up to 2500. At r = r, the plasma skin-depth is smallest and is resolved by 4.4 down to 1.4 cells depending on the value of σ. The highest plasma frequency is well resolved in all cases, with , and down to 4 at the highest σ. The shortest radiative cooling timescale must also be well resolved by the simulations, , where re = e2/mec2 is the classical radius of the electron and frad is a numerical factor that scales the magnitude of the cooling strength. Here, is a fixed number, meaning that in all the simulations regardless of the magnetic field strength. Simulations ran until the wind reached the outer boundary, r = 150 RLC, which corresponds to 24 spin periods. At this point, the box contained about 2.75 billions particles which represents about 160 particles per cell if spread uniformly over the box. Table 1 lists the set of simulations and the physical parameters that we explored in this study.

Table 1

List of numerical simulations performed in this work.

thumbnail Fig. 2

Test simulation of a pure magnetic monopole compared with the analytical solution of Michel (1973; dashed lines). From top to bottom: non-vanishing components of the magnetic field, components of the β = E × B/B2 drift velocity, and spindown power L normalized to as a function of radius.

Open with DEXTER

3. Simulations of the striped wind

3.1. Test case: the monopole

Before doing the full split-monopole simulations, we performed the test of a pure monopole for which we have an exact analytical solution (Michel 1973) that we can compare with the numerical result. This solution is characterized by the remarkable feature that Eθ = Bφ at all radii, where Bφ = Br(r/RLC) and (the other components are all zero2). We recover this solution numerically with good accuracy. Figure 2 makes a direct comparison between the analytical (in dashed lines) and the numerical solution (in solid lines). The top panel shows the radial profile of the non-zero components of the magnetic field (Br and Bφ). The middle panel illustrates even better the good agreement with the analytical solution. This plot compares the E × B drift velocity with the expected radial profiles, that is, The bottom panel shows the radial Poynting flux going through a spherical surface as a function of radius. The expected spindown power is (6)The numerical simulation demonstrates not only the good agreement with the analytical result, but also the low numerical dissipation of the code since no physical dissipation is expected in this particular setup. This is, of course, in contrast with the split-monopole case where reconnection kicks in and results in magnetic dissipation which is the main physical effect investigated here (see Sect. 3.3 below).

3.2. Overall structure

thumbnail Fig. 3

Snapshot of the striped wind simulation zoomed in on the innermost regions (r< 3.5 RLC) at time t = 16.6 of rotation periods in run R3_S250_K10. From top left to bottom right: the plasma multiplicity n/nGJ in log 10 scale, the toroidal magnetic field (r/r)Bφ/B, the poloidal current density (r/r)2Jθ normalized to the fiducial Goldreich-Julian current density JGJ ≈ ΩB/ 2π, and the dissipation rate of the electromagnetic energy (E·J) normalized by (r/r)3(BJGJ). Solid lines show the magnetic field lines. They are omitted in the bottom panels for clarity. Radii are in units of the light-cylinder radius (shown here by a dashed circle). The magenta arrow indicates the direction of the north magnetic pole (i.e., φ = Ωt).

Open with DEXTER

thumbnail Fig. 4

Entire numerical domain at the end of the simulation (t = 25.76P) in run R3_S250_K10. The grey (log-) scale shows the plasma density. Oblique dashed lines are Archimedean spirals of polar equation r(φ) /RLC = Ωtφ + π/ 2 (blue) and r(φ) /RLC = Ωtφ + 3π/ 2 (red). The vertical black dashed line is at r = RLC.

Open with DEXTER

We now describe the results of the striped wind simulations obtained from an initial split-monopole configuration. Figure 3 shows a snapshot of the simulation (run R3_S250_K10) after 16.6 spin periods and zoomed in on the innermost regions (r< 3.5 RLC). Shortly after the onset of the simulation, the purely radial magnetic structure collapses into the more familiar structure of the pulsar magnetosphere, that is, with closed field lines in the magnetic equatorial regions (φ = Ωt ± π/ 2) and open field lines at the magnetic poles (φ = Ωt and φ = Ωt + π). The closed field line region is confined within the light cylinder (shown by the dashed circle in Fig. 3) and is in solid rotation with the star. The striped wind begins beyond the light cylinder. It is composed of two equatorial current sheets that form at the boundary between opposite magnetic polarities (mostly in Bφ, see top-right panel) and a cold relativistic wind of particles. The current sheets have the shape of two nested Archimedean spirals of wavelength 2πRLC separated from each other by πRLC (most visible in Fig. 4). The electric current flowing in each sheet (along θ only) has opposite sign (bottom-left panel). This can easily be understood as the opposite signs of × B across each layer. The thickness of the current layers at the light cylinder, δLC, is approximatively the local plasma skin-depth, that is, δLC ~ de ≈ 10-2RLC in this run. The layer thickness is resolved by at least 6 cells at the light cylinder.

The density map in the top-left panel of Fig. 3 reveals a high-density contrast between the wind zone (the inter-stripe region, with local multiplicities κ ~ 110) and the current sheets (κ ~ 10100). There is also a strong density gradient between both sides of the reconnection layer, the leading part of the wind being less dense than the trailing part. The magnetic field strength follows the same trend. This resembles the asymmetric reconnection configuration as found, for example, in the Earth’s magnetopause. As soon as the current sheets form, they become quickly unstable to the relativistic tearing mode (Zelenyi & Krasnoselskikh 1979) which leads to the initial fragmentation of the sheet into overdense elliptical structures (plasmoids or magnetic islands below) separated by secondary current sheets (or simply referred to as X-points in the following). The fragmentation of secondary current sheets then proceeds on smaller scales and leads to the formation of a hierarchical chain of plasmoids and short current sheets via the plasmoid instability (Uzdensky et al. 2010). The initial tearing instability was already observed in three-dimensional (3D) pulsar magnetosphere simulations (Philippov et al. 2015b; Cerutti et al. 2016; Philippov & Spitkovsky 2017) but they did not capture the subsequent plasmoid instability because of the lack of sufficient scale separation in these studies.

These instabilities are also well-known in classical plane-parallel configurations to mediate efficient conversion of the Poynting flux into non-thermal particle acceleration within secondary current sheets (e.g., Cerutti et al. 2013; Sironi et al. 2016). This point is illustrated in the bottom-right panel of Fig. 3 which shows the dissipation rate of the electromagnetic energy E·J. There is a net positive dissipation rate localized between plasmoids where the reconnection (non-ideal) electric field is aligned with the current. Inside magnetic islands, the wind (ideal) electric field (E = −V × B/c) which changes sign on each side of the islands dominates, resulting in a dipolar-like structure in the E·J map. There is therefore no net dissipation of the Poynting flux within the plasmoids or anywhere else in the wind outside of the secondary current sheets. We discuss the process of particle acceleration in greater detail in Sect. 3.4.

Magnetic reconnection is most vigorous and time-dependent in the innermost parts of the striped wind (r ≲ 10 RLC). Plasmoids quickly grow from the inflow of particles coming from X-points on either side and by merging with other plasmoids. Mergers between small plasmoids are frequent close to the light cylinder. At intermediate radii (~2−5 RLC) merging episodes are still possible, although much less frequent, but they involve bigger structures. At larger distances, the substructures in the sheets are frozen by the adiabatic expansion. At the outer edge of the simulation box, about a dozen plasmoids remain (Fig. 4). The overall picture does not change with the values of initial plasma magnetization and multiplicity, and pulsar angular velocity explored here.

3.3. Magnetic reconnection and dissipation

3.3.1. Layer thickness and dissipation radius

Figure 4 clearly shows the broadening of the current sheets with radius. To quantify this, we measure the layer thickness δ at each X-point that we identify as a local density minimum along both spiral arms. We perform a Gaussian fit of the density profile along the direction perpendicular to the layer. We define the layer thickness as the full width at half maximum of the best-fit solution, δ = FWHM. Following Lyubarsky & Kirk (2001), we define the fraction of the striped wind wavelength filled by the current layer, Δ = δ/πRLC. We find that Δ increases linearly with radius (see top panel in Fig. 5). A linear fit gives (7)where ΔLC = δLC/πRLC ≈ 3.1 × 10-3 in this particular run. We find this correlation in all runs with very similar values of ΔLC (see Table 2). This result implies that by extrapolation the sheet would fill the entire wind region at (i.e., where Δ = 1), which would in turn lead to complete dissipation of the stripes. In this run, this would happen beyond 300 RLC which is unfortunately outside the numerical box.

Table 2

Current sheet filling factor at the light cylinder.

thumbnail Fig. 5

From top to bottom: current sheet filling factor (Δ = δ/πRLC), plasma drift velocity (βθ), dimensionless reconnection rate (), and electromagnetic energy flux toward the sheet (S′ = c|E′ × B′|/4π normalized by ) as function of radius in run R3_S250_K10. This analysis is done for each X-point (marked by a dot in these plots) along one spiral arm. The dashed line in the top panel is a linear fit to the numerical data of equation Δ(r) = ΔLC(r/RLC), where ΔLC ≈ 0.003. The horizontal dashed line gives the average dimensionless reconnection rate .

Open with DEXTER

The linear expansion of the sheet is consistent with the striped wind model of Coroniti (1990) and Michel (1994) which we briefly repeat here for completeness. Ampère’s law across the sheet for rRLC gives (8)where and are, respectively, the bulk drift velocity of the positrons and the electrons carrying the current3. Since Bφ ∝ 1 /r on the left-hand side and the plasma density n ∝ 1 /r2 on the right-hand side, the product βθδ should linearly increase with radius. Our analysis shown in Fig. 5 (middle panel) reveals that the plasma drift velocity is approximatively constant for rRLC, , which explains why Δ ∝ r. We notice that both species are perfectly counter-streaming in the sheet at all radii. The drift velocity is higher close to the light cylinder ( at rRLC) because Br cannot be neglected there and participates to the reconnecting field. This result also implies that the sheet always has enough charge to conduct the current required to maintain the discontinuity of the field (since βθ ≪ 1) as long as there is enough charge to begin with, as suggested by Arons (2012) but in contradiction with Usov (1975), Melatos & Melrose (1996). The same conclusion holds with the other values of σ, κ, and RLC explored here.

An important consequence of Eq. (7) is that the distance of complete dissipation is solely given by the thickness of the sheet at the light cylinder. Then, the relevant question is how does δLC depend on physical parameters. Using Eq. (8), and assuming that at the light cylinder (Fig. 5), gives (9)where ΓLC and are the bulk Lorentz factor and the plasma proper density at the light cylinder, respectively. The plasma multiplicity parameter is defined as . Then, Eq. (7) simply becomes (10)where the local plasma multiplicity is . This relation holds quite well at all radii in all the simulations within the range of the parameters explored here (see Fig. 6). The plasma multiplicity and bulk Lorentz factor at the light cylinder are about the same in all runs, κLC ~ 100 and ΓLC ≈ 1 (see Sect. 3.3.2), which explains why the layer thickness is very similar in all cases (even if κ and σ change). We reach the important conclusion that the striped wind would come to a complete dissipation at (where Δ = 1) (11)which is identical to Eq. (11) in Lyubarsky & Kirk (2001).

thumbnail Fig. 6

Distribution of the product ΔπΓLCκ computed at each X-point along the striped wind (run R3_S250_K10). The distribution peaks at ΔπΓLCκ ≈ 0.7 with a mean value of order unity (vertical red dashed line) as expected from Eq. (10).

Open with DEXTER

3.3.2. Reconnection rate and bulk velocities

Another key quantity to measure along the sheet is the dimensionless reconnection rate. It is usually defined as the ratio of the reconnection electric field induced at X-points and the upstream reconnecting magnetic field (i.e., just outside the layer). Equivalently, it is given by the speed at which fresh plasma and field enter and then leave an X-point, that is, as the ratio of the inflow to outflow plasma bulk velocities Vin/Vout. It is straightforward to measure this quantity in a 2D parallel configuration where the sheet has no bulk motion. In this context, the analysis is more complex because the relativistic bulk motion of the wind must be subtracted to extract the mildly relativistic motion induced by reconnection in the vicinity of each X-point. It is thus appropriate to measure the reconnection rate in the wind comoving frame (Lyubarskii 1996; Uzdensky & Spitkovsky 2014). Another difficulty is that the wind accelerates with radius and is not uniform along φ. Figure 7 shows the radial profile of the plasma bulk velocity averaged in φ. These profiles are consistent with the E × B drift velocity of a pure monopolar relativistic wind (see Fig. 2), but with some deviations close to the light-cylinder due to reconnection. The bulk Lorentz factor grows with radius, approximatively as (12)close to the light cylinder (Michel 1973; Buckley 1977; Contopoulos & Kazanas 2002), but quickly saturates to Γ ≈ 5−6 at large radii (r ≳ 20). The acceleration of the wind slows down after crossing the fast magnetosonic point located at r = rFMS ≈ 8.7 RLC and ΓFMS ≈ 3.4 for σ = 250 (indicated by the yellow star in Fig. 7), that is, where the wind velocity exceeds the Alfvén speed. This is given by the condition , where (13)(Michel 1969; Kirk et al. 2009). The plasma inertia becomes important beyond this point and the wind Lorentz factor increases only logarithmically with radius, approximatively as (14)in agreement with previous studies (Tomimatsu 1994; Beskin et al. 1998; Tchekhovskoy et al. 2009; Lyubarsky 2009). We see the same behavior in the other runs but, as expected, the distance of the fast point to the star increases with increasing magnetization: rFMS ≈ 33 RLC for σ = 1000 and rFMS ≈ 77 RLC for σ = 2500.

thumbnail Fig. 7

Top: radial profile of the plasma bulk velocity components V/c of the striped wind averaged in φ (solid lines) and the analytical E × B drift velocity of a pure monopole (dashed lines) at t = 25.76P in run R3_S250_K10. Bottom: corresponding φ-averaged bulk Lorentz factor of the striped wind (blue line). The wind first accelerates approximatively as (dashed line) and then grows as (dot-dashed line) passed the fast point (yellow star) where the (red line) and the Γ curves intersect.

Open with DEXTER

We use the φ-averaged bulk velocity profile shown in Fig. 7 to perform a Lorentz boost on each point of the striped wind, and neglect the effects of acceleration of the wind. Fields are changed according to the usual Lorentz transformation formulae, that is, where primed quantities are defined in the comoving frame. Far from the sheets, the electric field vanishes E′ ≈ 0 and the magnetic field shrinks as B′ ≈ B/ Γ, as expected if E·B = 0. We compute the dimensionless reconnection rate, , as the ratio of the electric field induced inside and around each X-point by reconnection, , with the upstream magnetic field strength B. We obtain a dimensionless reconnection rate around 0.4 throughout the wind (see bottom panel in Fig. 5). These values are high compared with previous local plane parallel simulations of reconnection (with a rate of order 0.10.2). A possible explanation is that reconnection is driven by large-scale plasma motions in the wind rather than being spontaneous as is almost always the case in local simulations. Intense synchrotron losses may also be at the origin of this enhancement as anticipated by Uzdensky & McKinney (2011) and as observed in Cerutti et al. (2013). Figure 5 (bottom panel) shows the module of the Poynting vector defined in the comoving frame, S′ = c|E′ × B′|/4π, normalized by . This quantity measures the amount of electromagnetic energy flux flowing towards X-points along the striped wind. Even though the reconnection rate remains constant, the available energy flux sharply decreases with radius because if Γ is constant.

thumbnail Fig. 8

Top: radial profile of the toroidal magnetic field strength r|Bφ| (blue line) and its striped-average (red line) for φ = 0 at t = 25.76P in run R3_S250_K10. Bottom: radial profile of the Poynting flux (L/L, blue line) and the plasma magnetization parameter (σ = B2/ 4πΓnmec2, red line).

Open with DEXTER

Relativistic magnetic reconnection leads to efficient magnetic dissipation. Figure 8 shows the decay with radius of the stripe-averaged toroidal magnetic field strength consumed by reconnection in the sheet and converted into heat and particle acceleration in run R3_S250_K10. This translates into a monotonic decrease of the Poynting flux with radius, from L ≈ 0.4 L at the light cylinder down to L ≈ 0.05 L at r = 150 RLC (i.e., a drop of almost 90%), without any sign of saturation. The plasma magnetisation drops accordingly with increasing radius down to σ of order unity at r = 150 RLC, due to both increasing wind Lorentz factor and decreasing magnetic field strength. This is also consistent with the increase of the sheet filling factor with radius. We notice a more efficient dissipation close to the light cylinder where the sheet is thinner and where the tearing instability is most active. At large distances, reconnection proceeds in a smoother, more laminar way which results in an approximate constant dissipation rate.

3.4. Particle acceleration and radiative signatures

thumbnail Fig. 9

Spatial distribution of the particle mean Lorentz factor (top) and synchrotron flux (bottom, in arbitrary units) zoomed in on the innermost regions of the striped wind in run R3_S250_K10. The plasma density isocontours (black lines) show the location of plasmoids.

Open with DEXTER

thumbnail Fig. 10

Trajectories of 30 simulation particles whose Lorentz factor (color-coded) exceeded γ = 30 at least once over their history in run R3_S250_K10. Plasma density at t = 5.5P is shown in gray.

Open with DEXTER

The magnetic energy dissipated via reconnection is channeled into particle kinetic energy. As was already noticed in the map of E·J in Fig. 3, dissipation of the Poynting flux happens primarily between plasmoids, that is, within secondary current layers where the reconnection electric field accelerates particles. This picture is consistent with local plane parallel PIC simulations (Zenitani & Hoshino 2001; Cerutti et al. 2012b; Sironi & Spitkovsky 2014; Nalewajko et al. 2015) which identified X-points as the main particle accelerating regions.

Figure 9 (top) shows the spatial distribution of the mean particle Lorentz factor zoomed in on the innermost regions of the striped wind. This figure. clearly demonstrates that close to the light cylinder high-energy particles are concentrated within secondary layers only and not in plasmoids, in contrast with previous reconnection studies. This important difference is attributed to the strong radiative losses present in this study. Particles are first accelerated at X-points where the effective perpendicular magnetic field, given by (17)is small and, hence, radiative losses are small (β is the particle 3-velocity divided by c). But as soon as the particles reach a plasmoid, they feel an abrupt increase of which results in catastrophic radiative energy losses (where γ = (1−β2)− 1/2 is the particle Lorentz factor) and an intense emission of energetic synchrotron radiation (Cerutti et al. 2013, 2014). In other words, plasmoids act as beam dumps and are the main emitting regions of the striped wind. This is well visible in Fig. 9 which shows an almost perfect spatial anti-correlation between the synchrotron flux and the mean Lorentz factor maps.

Figure 10 presents the trajectories of 30 simulation particles which were accelerated at least once over their history. A closer look at individual trajectories reveal multiple episodes of acceleration and cooling which depend on the evolution of the particles within the layer, and on the dynamics of plasmoids and mergers. Kinks in the trajectories are almost always associated with the capture of the particles by a plasmoid and by rapid radiative losses. We note that this connection cannot be directly observed in Fig. 10 because the particle trajectories are over-plotted onto the density map shown at a given time (5.5 rotation periods).

thumbnail Fig. 11

Top: mean particle Lorentz factor (γ, blue line) and mean synchrotron photon energy (ν ⟩ /ν0, where ν0 = 3eB/ 4πmec, red line) as functions of radius in run R3_S250_K10. Bottom: radial profile of the synchrotron flux normalized by its maximum value.

Open with DEXTER

thumbnail Fig. 12

Particle energy spectra (γdN/ dγ, top) and spectral energy distributions (νFν, bottom) averaged in φ as a function of radius (color-coded) for σ = 2500 (run R3_S2500_K10). The vertical dashed line in the top panel shows γ = σLC ≈ 500.

Open with DEXTER

thumbnail Fig. 13

Particle spectra at the outer edge of the simulation box r = 150 RLC for σ = 250 (green), 1000 (red), and 2500 (blue). The inset plot shows the mean particle Lorentz factor γ as function of σLC. The oblique dashed line is γ ⟩ = σLC.

Open with DEXTER

Further away (rRLC), particles continue to accelerate because magnetic reconnection is still active and synchrotron radiative losses decrease with radius (Fig. 11, top). Figure 12 (top panel) shows the φ-averaged particle spectra γdN/ dγ as a function of radius for σ = 2500 (run R3_S2500_K10). At small radii, the particle spectra are dominated by a narrow component at low energy (γ a few) which corresponds to the mild energization that experienced the particles at the surface of the star after their injection. Then, the spectrum develops a steep power-law tail that hardens with increasing radius. At large radii (r ≳ 100 RLC), the spectrum converges towards a narrow distribution centred around γ ≈ 500, which coincides with the value of the magnetization parameter at the light-cylinder radius σLCσ(RLC) (Philippov & Spitkovsky 2014; Cerutti et al. 2015, 2016; Philippov & Spitkovsky 2017), regardless of the strength of radiative cooling (Kirk 2004; Cerutti et al. 2012a; Kagan et al. 2016). This correlation is confirmed by the other simulations where γ ⟩ ≈ σLC (see Fig. 13). This relation can be recast as (Cerutti et al. 2015) (18)where φpc = eΦpc/mec2 and Φpc = μΩ2/c2 is the vacuum potential drop across the polar cap. In comparison, the radiation-reaction-limited particle Lorentz factor at the light cylinder is always much smaller than σLC in our simulations. It can be estimated by the following expression (e.g., Cerutti et al. 2013): (19)Here, varies from 12 for σLC = 50 to for σLC ≈ 500. We note that in Crab-like pulsars this quantity can be as high as .

It is worth noticing here that in contrast to plane parallel reconnection studies at very high magnetization (Sironi & Spitkovsky 2014; Werner et al. 2016; Guo et al. 2014), the particle spectra do not converge to a broad, hard (p ≈ −1) power law. We see hints of such a hard power-law tail at intermediate radii only. Instead, all particles end up gaining a very similar amount of energy regardless of their histories. The characteristic energy is naturally set by σLC because it represents the total available magnetic energy per particle. The fact that the sheets become so close to each other (and eventually merging with each other) at large radii could be a possible explanation for this important difference with plane parallel studies where the sheet is isolated or as far away as possible from its neighbor thus preventing any interaction between them. This result is also consistent with the linear expansion of the sheet with radius found in the previous Section (Sect. 3.3.1, Eqs. (7)(10)). Assuming that the layer thickness scales with the mean particle Larmor radius in the sheet δ ~ ⟨ γmec2/eB and with γ ⟩ ≈ σLC ~ φpc/ ΓLCκLC (we dropped the factor 4) yields , which coincides with the expression found in Eq. (10). One can also recover this relation by setting the layer thickness to the local plasma skin-depth, δ ~ de = ( ⟨ γmec2/ 4πne2)1/2.

The bottom panel in Fig. 12 shows the resulting spectral energy distribution (νFν) averaged in azimuth as function of radius. Synchrotron radiation is the dominant process at the origin of the high-energy emission. In contrast with the particles, the mean synchrotron photon frequency decreases with radius (Fig. 11). The highest-energy photons are found close to but outside of the light cylinder where there are both energetic particles and a strong magnetic field. The frequency-integrated flux is also concentrated at the light cylinder within a few RLC (Fig. 11, bottom panel). It is therefore sufficient to consider the innermost parts of the striped wind to model pulse profiles. Following Cerutti et al. (2016), we compute the synchrotron flux as received by an observer at infinity looking in the plane of the simulation. Figure 14 presents a series of individual light curves computed after each pulsar period and folded over the pulsar phase whose origin is defined by the direction of the north magnetic pole. The obtained light curves present two pulses at phase 0.25 and 0.75, in agreement with previous 3D simulations for an observer looking along the equator (Cerutti et al. 2016; Philippov & Spitkovsky 2017).

thumbnail Fig. 14

Top: series of individual pulse profiles computed after each pulsar period. Bottom: pulse profile averaged over all periods. The origin of phases is defined by the direction of the north magnetic pole.

Open with DEXTER

Thanks to the long integration time of the simulations, we can study the pulse-to-pulse variability. We find bright substructures randomly distributed within each pulse, reminiscent of micropulses observed in some pulsars in radio (Lyne & Graham-Smith 2006). These substructures are emitted by the particles confined within the plasmoids and shining towards the observer. Since most of the emission comes from the base of the striped wind, the synchrotron flux probes the most active reconnecting regions, that is, where there is the largest number of plasmoids and plasmoid mergers. However, due to the low photon flux at high energies, instruments are only sensitive to the average pulse profile. Averaging over all pulses washes out the intra-pulse variability and leaves the two main peaks which are very stable from one pulse to another (bottom panel in Fig. 11).

4. Conclusion

We have investigated the formation and the dissipation of the striped pulsar wind structure using large 2D PIC simulations of a rotating split-monopole. We observe that shortly after its formation, the wind current sheet breaks up into a highly dynamical chain of magnetic islands separated by secondary current sheets, a configuration reminiscent of previous plane parallel reconnection studies. Close to the light cylinder, the dynamics of plasmoids is very time-dependent: they form, expand and merge with each other multiple times. This is the most dynamical, turbulent-like phase of reconnection which is accompanied by efficient dissipation of the Poynting flux (Zrake 2016). At larger radii, adiabatic expansion takes over and no more plasmoids form or merge. Reconnection proceeds in a more laminar way and continues to dissipate the Poynting flux without any sign of saturation up to at least 150 RLC where the numerical box ends. The current sheet thickness expands linearly with radius in accordance with the increase of the plasma skin-depth in the wind.

Complete dissipation does not occur within the numerical domain (up to 90%) but simulations suggest that the dissipation radius is set by the plasma multiplicity at the light cylinder, rdiss/RLC = πΓLCκLC. This corresponds to the location where the layer thickness is equal to the wavelength of the stripes. This remarkably simple result has important consequences. With multiplicities of order κ ~ 103−104 (Hibschman & Arons 2001; Timokhin & Arons 2013) and ΓLC ≲ 100, our results imply that the stripes should dissipate far before the wind reaches the termination shock in isolated systems, even in the Crab pulsar where Rshock/RLC ~ 109 ≫ ΓLCκ ~ 105−106 (Coroniti 1990; Michel 1994). We note that this conclusion is not consistent with Lyubarsky & Kirk (2001) and Kirk & Skjæraasen (2003). A possible explanation for this discrepancy could be that reconnection is very effective at dissipating magnetic energy in the inner zone of the striped region (about half is dissipated), where the wind is not super-Alfvenic. This transition region is not considered by the analytical models but it is well captured by the simulations. Another reason could be that we find no significant acceleration of the flow in the super-fast wind region as the stripes dissipate and hence no time-dilation effect to slow down the process in contrast to what Lyubarsky & Kirk (2001) found. This could be an artefact of the relatively low magnetization explored in our simulations.

What happens beyond the dissipation radius is unknown. What is clear, however, is that the current sheet will not lack electric charge, as long as there is enough charge at the light cylinder, because the particle drift velocities needed to carry the current are subluminal and decrease with radius (Arons 2012). Magnetic dissipation results in efficient particle acceleration and emission of high-energy pulsed synchrotron emission. Even though magnetic dissipation is spread over hundreds of light-cylinder radii, the synchrotron flux piles up at the light cylinder. This implies that the pulsed high-energy emission probes the innermost parts of the pulsar wind where reconnection is most time-dependent, resulting in significant pulse-to-pulse variability (Cerutti et al. 2016). The particle spectra asymptotically converge to a narrow energy distribution with no sign of a high-energy power-law tail in contrast with plane parallel reconnection studies (Sironi & Spitkovsky 2014; Guo et al. 2014; Werner et al. 2016). The characteristic particle Lorentz factor in the wind after dissipation is set by the magnetization parameter at the light cylinder γσLC which also depends on the plasma multiplicity parameter showing once again how critical this parameter is in this problem. It may be determined by both the plasma supply from the star via the polar-cap discharge and photon-photon annihilation within the current sheet (Lyubarskii 1996; Philippov & Spitkovsky 2017). While polar-cap cascade was extensively studied in the past, pair creation within the current sheet is still poorly understood. Dedicated studies of this phenomenon would be very valuable.

To conclude, the striped component of the pulsar wind is most likely fully annihilated far before it reaches the nebula. The wind is left with its aligned DC component (not present in this work) which may remain highly magnetized up to the termination shock. Three-dimensional simulations would be required to investigate the full latitude-dependent evolution of the striped wind. Note that even though the striped component disappears, a large-scale current sheet survives in the equatorial regions which might be a promising location for extra particle acceleration and magnetic dissipation in pulsar wind nebulae. Pulsars in binary systems such as gamma-ray binaries (Dubus 2013) or transitional millisecond pulsars (Archibald et al. 2009; Stappers et al. 2014) where the pulsar wind zone is truncated at much smaller radii than in isolated systems (by the companion star wind or by an accretion disk) may provide ideal environments to probe the innermost regions of the wind where we expect dissipation to be most active.


1

It is worth noting that a relativistic pulsar wind cannot be achieved in a cylindrical geometry, more specifically we find that EzBφ in this case.

2

We note that we use the same notation for the spherical and cylindrical radii here (both written r) since it makes no difference in the equator.

3

We note that this expression assumes a neutral plasma, ne+ne, which is the case in this setup (see Sect. 2).

Acknowledgments

We thank A. Beloborodov, G. Dubus, A. Levinson, Y. Lyubarsky, E. Nakar, and A. Spitkovsky for discussions. B.C. acknowledges support from CNES and Labex OSUG@2020 (ANR10 LABX56). This research was also supported by the NASA Earth and Space Science Fellowship Program (grant NNX15AT50H to AP), Porter Ogden Jacobus Fellowship awarded by Princeton University to AP, NASA through Einstein Postdoctoral Fellowship grant number PF7-180165 awarded by the Chandra X-ray Center to AP, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. The simulations presented in this article used computational resources of TGCC under the allocation t2016047669 made by GENCI, of the NASA/Ames HEC Program (SMD-16-7816), as well as NSF through an XSEDE computational time allocation TG-AST100035 on TACC Stampede.

References

All Tables

Table 1

List of numerical simulations performed in this work.

Table 2

Current sheet filling factor at the light cylinder.

All Figures

thumbnail Fig. 1

Numerical setup of a two-dimensional (2D) PIC simulation of the striped wind in the equator of an inclined split-monopole (Bogovalov 1999). Initially, magnetic field lines are purely radial and their polarity change at φ = Ωt ± π/ 2 (dotted line). The box extends from the stellar surface rmin = r (black disk) where e± pairs are created and where the fields are in co-rotation with the star, up to an absorbing layer between rabs and rmax where particles and fields vanish. The light-cylinder radius RLC = c/ Ω is shown by the dashed circle.

Open with DEXTER
In the text
thumbnail Fig. 2

Test simulation of a pure magnetic monopole compared with the analytical solution of Michel (1973; dashed lines). From top to bottom: non-vanishing components of the magnetic field, components of the β = E × B/B2 drift velocity, and spindown power L normalized to as a function of radius.

Open with DEXTER
In the text
thumbnail Fig. 3

Snapshot of the striped wind simulation zoomed in on the innermost regions (r< 3.5 RLC) at time t = 16.6 of rotation periods in run R3_S250_K10. From top left to bottom right: the plasma multiplicity n/nGJ in log 10 scale, the toroidal magnetic field (r/r)Bφ/B, the poloidal current density (r/r)2Jθ normalized to the fiducial Goldreich-Julian current density JGJ ≈ ΩB/ 2π, and the dissipation rate of the electromagnetic energy (E·J) normalized by (r/r)3(BJGJ). Solid lines show the magnetic field lines. They are omitted in the bottom panels for clarity. Radii are in units of the light-cylinder radius (shown here by a dashed circle). The magenta arrow indicates the direction of the north magnetic pole (i.e., φ = Ωt).

Open with DEXTER
In the text
thumbnail Fig. 4

Entire numerical domain at the end of the simulation (t = 25.76P) in run R3_S250_K10. The grey (log-) scale shows the plasma density. Oblique dashed lines are Archimedean spirals of polar equation r(φ) /RLC = Ωtφ + π/ 2 (blue) and r(φ) /RLC = Ωtφ + 3π/ 2 (red). The vertical black dashed line is at r = RLC.

Open with DEXTER
In the text
thumbnail Fig. 5

From top to bottom: current sheet filling factor (Δ = δ/πRLC), plasma drift velocity (βθ), dimensionless reconnection rate (), and electromagnetic energy flux toward the sheet (S′ = c|E′ × B′|/4π normalized by ) as function of radius in run R3_S250_K10. This analysis is done for each X-point (marked by a dot in these plots) along one spiral arm. The dashed line in the top panel is a linear fit to the numerical data of equation Δ(r) = ΔLC(r/RLC), where ΔLC ≈ 0.003. The horizontal dashed line gives the average dimensionless reconnection rate .

Open with DEXTER
In the text
thumbnail Fig. 6

Distribution of the product ΔπΓLCκ computed at each X-point along the striped wind (run R3_S250_K10). The distribution peaks at ΔπΓLCκ ≈ 0.7 with a mean value of order unity (vertical red dashed line) as expected from Eq. (10).

Open with DEXTER
In the text
thumbnail Fig. 7

Top: radial profile of the plasma bulk velocity components V/c of the striped wind averaged in φ (solid lines) and the analytical E × B drift velocity of a pure monopole (dashed lines) at t = 25.76P in run R3_S250_K10. Bottom: corresponding φ-averaged bulk Lorentz factor of the striped wind (blue line). The wind first accelerates approximatively as (dashed line) and then grows as (dot-dashed line) passed the fast point (yellow star) where the (red line) and the Γ curves intersect.

Open with DEXTER
In the text
thumbnail Fig. 8

Top: radial profile of the toroidal magnetic field strength r|Bφ| (blue line) and its striped-average (red line) for φ = 0 at t = 25.76P in run R3_S250_K10. Bottom: radial profile of the Poynting flux (L/L, blue line) and the plasma magnetization parameter (σ = B2/ 4πΓnmec2, red line).

Open with DEXTER
In the text
thumbnail Fig. 9

Spatial distribution of the particle mean Lorentz factor (top) and synchrotron flux (bottom, in arbitrary units) zoomed in on the innermost regions of the striped wind in run R3_S250_K10. The plasma density isocontours (black lines) show the location of plasmoids.

Open with DEXTER
In the text
thumbnail Fig. 10

Trajectories of 30 simulation particles whose Lorentz factor (color-coded) exceeded γ = 30 at least once over their history in run R3_S250_K10. Plasma density at t = 5.5P is shown in gray.

Open with DEXTER
In the text
thumbnail Fig. 11

Top: mean particle Lorentz factor (γ, blue line) and mean synchrotron photon energy (ν ⟩ /ν0, where ν0 = 3eB/ 4πmec, red line) as functions of radius in run R3_S250_K10. Bottom: radial profile of the synchrotron flux normalized by its maximum value.

Open with DEXTER
In the text
thumbnail Fig. 12

Particle energy spectra (γdN/ dγ, top) and spectral energy distributions (νFν, bottom) averaged in φ as a function of radius (color-coded) for σ = 2500 (run R3_S2500_K10). The vertical dashed line in the top panel shows γ = σLC ≈ 500.

Open with DEXTER
In the text
thumbnail Fig. 13

Particle spectra at the outer edge of the simulation box r = 150 RLC for σ = 250 (green), 1000 (red), and 2500 (blue). The inset plot shows the mean particle Lorentz factor γ as function of σLC. The oblique dashed line is γ ⟩ = σLC.

Open with DEXTER
In the text
thumbnail Fig. 14

Top: series of individual pulse profiles computed after each pulsar period. Bottom: pulse profile averaged over all periods. The origin of phases is defined by the direction of the north magnetic pole.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.