Issue 
A&A
Volume 656, December 2021



Article Number  A91  
Number of page(s)  8  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202142178  
Published online  06 December 2021 
Formation of giant plasmoids at the pulsar wind termination shock: A possible origin of the innerring knots in the Crab Nebula^{⋆}
^{1}
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: benoit.cerutti@univgrenoblealpes.fr
^{2}
MaxPlanckInstitut für Kernphysik, Postfach 103980, 69029 Heidelberg, Germany
email: gwenael.giacinti@mpihd.mpg.de
Received:
8
September
2021
Accepted:
4
November
2021
Context. Nearby pulsar wind nebulae exhibit complex morphological features: jets, torus, arcs, and knots. These structures are well captured and understood in the scope of global magnetohydrodynamic models. However, the origin of knots in the inner radius of the Crab Nebula remains elusive.
Aims. In this work, we investigate the dynamics of the shock front and downstream flow with a special emphasis on the reconnecting equatorial current sheet. We examine whether giant plasmoids produced in the reconnection process could be good candidates for the knots.
Methods. To this end, we perform large semiglobal threedimensional particleincell simulations in a spherical geometry. The hierarchical merging plasmoid model is used to extrapolate numerical results to pulsar wind nebula scales.
Results. The shocked material collapses into the midplane, forming and feeding a largescale, but thin, ringlike current layer. The sheet breaks up into a dynamical chain of merging plasmoids, reminiscent of threedimensional reconnection. Plasmoids grow to a macroscopic size. The final number of plasmoids predicted is solely governed by the inverse of the dimensionless reconnection rate.
Conclusions. The formation of giant plasmoids is a robust feature of pulsar wind termination shocks. They provide a natural explanation for the innerring knots in the Crab Nebula, provided that the nebula is highly magnetized.
Key words: acceleration of particles / magnetic reconnection / radiation mechanisms: nonthermal / methods: numerical / pulsars: general / stars: winds, outflows
Movie is available at https://www.aanda.org
© B. Cerutti and G. Giacinti 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Highresolution Xray images of nearby pulsar wind nebulae, as best illustrated by the iconic Chandra images of the Crab Nebula (Weisskopf et al. 2000; Hester et al. 2002; Mori et al. 2004) and of the Vela pulsar wind nebula (Helfand et al. 2001; Pavlov et al. 2001a; Durant et al. 2013), have uncovered complex morphological features. This played a crucial role in the recent development of pulsar wind nebula theory (see, e.g., Kargaltsev et al. 2015 for a review and references therein) and, more generally, the study of relativistic magnetized outflows.
At the largest scales, a pulsar wind nebula can be well described as being composed of a jetlike structure aligned along the pulsar rotation axis and of a bright toruslike structure lying in the equatorial plane. The pulsar is seemingly located at the center of the torus and near the base of the jet. At smaller scales, and currently only visible in the Crab and Vela pulsar wind nebulae, one can distinguish bright thin, arclike structures within the torus that are moving away from the pulsar, also known as the “wisps” (Scargle 1969), and an “inner ring” located at the inner edge of the torus (Weisskopf et al. 2000). A closer look at the Crab inner ring reveals that it is in fact composed of a series of about two dozen bright compact knots (Weisskopf et al. 2000; Hester et al. 2002), and at least two knots are also reported in Vela (Pavlov et al. 2001b). These knots are highly variable, although they seem stationary and confined within the inner ring. They appear as unpolarized structures in the optical band (Hester 2008), and they should not be confused with the “inner knot”. The latter is another salient feature of the Crab Nebula discovered with the Hubble Space Telescope in the optical (Hester et al. 1995); it is highly variable and polarized (Moran et al. 2013), but it is too close to the pulsar projected on the sky to be singled out as a separate component in Xrays. Figure 1 schematically represents all of the aforementioned morphological features.
Fig. 1.
Diagram schematically representing the morphological features observed in the Crab Nebula in the optical and Xray bands. The MHD model captures all these features well, except the knots along the inner ring (in red), whose origin is unknown. Solving this mystery is the main focus of this work. 
It is now well understood that the jettorus structure stems from the interaction of the anisotropic relativistic magnetized wind blown by the pulsar with the surrounding medium (Bogovalov & Khangoulian 2002; Lyubarsky 2002). A shock forms from this interaction, leading to compression and heating of the pulsar wind and (re)acceleration of the relativistic electronpositron pairs it contains; this gives rise to bright synchrotron radiation (Rees & Gunn 1974; Kennel & Coroniti 1984a,b). Twodimensional (2D) axisymmetric magnetohydrodynamic (MHD) simulations have successfully reproduced the observed jettorus structure, which was identified as the plasma downstream of the pulsar wind termination shock (Komissarov & Lyubarsky 2003, 2004; Del Zanna et al. 2004). In this framework, the moving wisps result from the feedback of the highly dynamical downstream flow on the shock front location (Camus et al. 2009), while the inner knot originates from the Dopplerboosted emission of the termination shock pointing toward the observer (Komissarov & Lyubarsky 2004; Komissarov & Lyutikov 2011). Threedimensional (3D) MHD simulations confirmed and accurately reproduced the morphology and time evolution of the Crab Nebula down to the smallest details (Porth et al. 2013, 2014), giving us high confidence in this physical model.
The morphological features that are not captured in the scope of the 3D MHD model are the brightness of the inner ring and its knots (Porth et al. 2014). This region is often regarded as the location of the termination shock itself in the equatorial plane (Weisskopf et al. 2000). Yet, it stands out from the rest of the nebula because this is where both polarities of the (mostly toroidal) magnetic field meet and reconnect (Giacinti & Kirk 2018). Using 2D particleincell (PIC) simulations, we have shown in a previous study that the anisotropy of the wind naturally results in the formation of a largescale, but thin, reconnecting current sheet in the midplane, which plays a key role in particle acceleration in the nebula (as well as shear flow acceleration; see Cerutti & Giacinti 2020). Therefore, the equatorial current sheet is a good candidate for explaining the brightness of the inner ring, but its knotty nature could not be addressed with these 2D simulations.
It is now well established that a long, thin current sheet is prone to the tearing and plasmoid instabilities (Zelenyi & Krasnoselskikh 1979; Loureiro et al. 2007; Uzdensky et al. 2010) that lead to the fragmentation of the sheet into multiple secondary current layers separated by plasma overdensities trapped in magnetic islands, or flux ropes in 3D. As reconnection proceeds, magnetic islands merge to form bigger structures, eventually giving rise to giant, macroscopic plasmoids (Loureiro et al. 2012) filled with energetic particles. In this work, we investigate the dynamics of the shock front and downstream flow in light of more realistic semiglobal 3D PIC simulations in a spherical geometry. We draw special attention to the development of flux ropes and examine whether these structures could be good candidates for explaining the observed knots along the inner ring of the Crab Nebula. In Sect. 2 we describe our numerical setup, which was inspired from Cerutti & Giacinti (2020) and adapted to a spherical geometry. In Sect. 3 we present the dynamics of the shock, with an emphasis on the formation of the equatorial current sheet and the flux ropes it contains. We provide a detailed analysis of their evolution and propose a toy model that can be used to apply and extrapolate our results to any pulsar wind nebula.
2. Numerical setup
We worked with the Zeltron PIC code using 3D spherical coordinates (r, θ, ϕ; Cerutti et al. 2013, 2016). The numerical domain is a spherical wedge of the inner radius, r_{min}, and outer radius, r_{max} = 10 r_{min}, limited by the polar angle ranging from θ_{min} = π/4 to θ_{max} = 3π/4, and azimuth ranging from ϕ_{min} = 0 to ϕ_{max} = π/2. The polar axis (θ = 0) is aligned along the pulsar rotation axis. We removed the polar regions to concentrate our computational resources on the equatorial regions. We verified that this choice does not significantly affect the dynamics of the flow in the equatorial belt and the results presented in this work. All runs were composed of (1600 × 1024 × 1024) cells along the r, θ, and ϕ directions, respectively. The grid cells were equally spaced in log r, cos θ, and ϕ to optimize the numerical resolution. The outer radius was a perfectly reflecting boundary with no loss of energy for both the electromagnetic fields and the particles. The θ boundaries were periodic for the particles but reflecting for the fields, and the ϕ boundaries were periodic for both the fields and the particles.
The numerical box was initialized with a perfectly radial wind composed exclusively of electronpositron pairs propagating outward with a bulk velocity of V_{0}/c = 0.99, or a bulk Lorentz factor of Γ_{0} ≈ 7. This choice satisfies both the need to model an ultrarelativistic flow, as found in pulsar winds (i.e., Γ_{0} ≫ 1), and the need to limit spurious heating of the beam due to numerical Cherenkov radiation, which limits the maximum Lorentz factor in the simulation to Γ ≲ 10. The shock forms when the initial beam of plasma hits and reflects on the radial outer boundary, giving rise to two identical but counterpropagating beams (e.g., Sironi et al. 2013). This setup is numerically convenient and easy to implement but has the disadvantage of producing a reverse shock that propagates inward; this is not the case in the Crab Nebula, where it is approximately stationary in the frame of the observer. The outer boundary can be seen as the contact discontinuity separating the pulsar wind material from the ambient medium, which is not simulated here. Therefore, the numerical experiment presented here is valid in the frame of the contact discontinuity of the flow, which is presumably still moving outward in young pulsar wind nebulae such as the Crab.
This far from the pulsar magnetospheric regions where the wind is launched, the magnetic field in the upstream medium entering the shock can be considered as purely toroidal (Michel 1973; Bogovalov 1999). We further assumed that the oscillating component of the field shaped by the pulsar rotation when the rotation and magnetic axes of the star are misaligned at an angle χ ≠ 0 (i.e., the striped wind; Coroniti 1990; Michel 1994) has fully dissipated when it reaches the termination shock (Coroniti 1990; Cerutti & Philippov 2017; Cerutti et al. 2020). A good proxy for the field is given by^{1} (Cerutti & Giacinti 2020)
B_{r} = 0 and B_{θ} = 0. The sin θ/r term is borrowed from the splitmonopole model (Michel 1973; Bogovalov 1999), a good description of the asymptotic dissipationless pulsar wind structure, while the tanh((θ − π/2)/χ) term translates the full dissipation of the striped wind component. The electric field in the wind is set by the ideal Ohm’s law, E = −V × B/c, giving a single nonvanishing component, E_{θ} = V_{0}B_{ϕ}/c.
In the steady state, the wind must carry the current density J = c/(4π)∇ × B. Using Eq. (1) yields
J_{θ} = 0 and J_{ϕ} = 0. Ohm’s law combined with Gauss’ law indicates that the wind has a net charge, given by the following charge density,
where
Following Cerutti & Giacinti (2020), we assumed that the positive charge density, ρ_{+}, is carried away by positrons only of number density n_{+} = ρ_{+}/e (where e is the elementary electric charge) and ρ_{−} by electrons of density n_{−} = −ρ_{−}/e, both species propagating outward in the wind at the speed V_{0}. The pulsar wind must also fulfill quasineutrality (i.e., ρ≪ne). This condition is easily achieved by adding a uniform plasma density of n_{w} = n_{0}(r_{min}/r)^{2} to both the electrons and positrons, that is drifting outward at V_{0} such that the total electron density is n_{e} = n_{w} + n_{−} and the total positron density is n_{p} = n_{w} + n_{+}. The injected plasma density profile in the wind is modeled by four particles per cell. Particles are continuously added to the wind zone from an injector that is initially located at r = 0.95 r_{max} and then recedes at the speed of light toward the inner boundary as the simulation continues. The run stops once the injector has reached the inner boundary, that is, t_{max} = 8.5r_{min}/c; this represents about 12 000 time steps, which are fixed at half the CourantFriedrichLewy stability condition.
The fiducial plasma density, n_{0}, and magnetic field, B_{0}, set the upstream plasma magnetization parameter, the only free parameter in this study, which is defined as
where m_{e} is the electron mass. In this work, we performed runs with σ_{0} = 3, 10, 30, and 100. Taking the polar profile of B_{ϕ} into account and fixing χ = π/4 gives an average magnetization, σ, of the order σ ≈ 0.15σ_{0}. Therefore, we effectively explored σ = 0.45, 1.5, 4.5, and 15, meaning from mildly to highly magnetized nebulae. The 3D MHD model favors magnetization of order unity and perhaps higher (Porth et al. 2013, 2014). The physical size of the box is identical in all runs in terms of electron plasma skin depth, . The spatial resolution is about five cells per electron skin depth in all directions, therefore fitting about in the simulation box along the ϕ direction.
3. Results
3.1. Dynamics and evolution of the shock
The shock forms immediately after the reflection of the beam on the outer boundary. The shock front recedes faster inward as one moves away from the equatorial plane, where the flow is more magnetized and therefore where the shock is weaker, leading to a significant departure from the spherical symmetry of the shock front. The downstream flow is decelerated to mildly relativistic speeds. The most striking feature is the development of a largescale ringlike current sheet in the equatorial plane due to the sharp magnetic discontinuity forming behind the shock. The downstream flow collapses into the midplane and conducts the electric current via a thin, skindepthscale layer. This layer then fragments into multiple interconnected flux ropes along the azimuthal direction, as shown in the 3D rendering of the final plasma density state for σ = 4.5 (Fig. 2). The current layer is continuously fed with fresh plasma thanks to a largescale circulation pattern that takes place in the downstream medium and brings highlatitude magnetized plasma down to the midplane (see Fig. 3, panel a). This mechanism even leads to a reversal of the radial flow velocity in the vicinity of the layer, which in turn excites KelvinHelmholtz modes, as reported in Cerutti & Giacinti (2020) and as conjectured by Begelman (1999). Due to the limited dynamical range probed in these 3D runs, only the early stages of this instability can be captured in the outermost regions of the box (r/r_{min} ≳ 8; see panel a and the bottom of panel c in Fig. 3). At later stages, the KelvinHelmholtz instability may ultimately dominate and lead to the disruption and transformation of this largescale current sheet into a turbulent flow (Cerutti & Giacinti 2020), hence limiting the radial extent of the layer and flux ropes. The shock front cavity reported in Cerutti & Giacinti (2020) is not observed in these simulations, most likely due to the small system size and integration time.
Fig. 2.
Threedimensional rendering of the final state of the plasma density (r/r_{min})^{2}n/n_{0} for the σ = 4.5 run. For the sake of clarity, axes are not shown; see Fig. 3 for 2D slices with axes and length scales. 
This network of interconnected flux ropes, reminiscent of 3D plane parallel relativistic reconnection studies (e.g., Kagan et al. 2013; Cerutti et al. 2014; Werner & Uzdensky 2017, 2021), is clearly visible in the equatorial slice of the plasma density shown in Fig. 3, panel b. This figure is remarkable in the sense that it captures in a single snapshot the time evolution and the expected merging tree pattern of flux ropes forming in the layer (see for instance the merging tree in the spacetime diagram presented in Nalewajko et al. 2015). In the θϕ plane, the evolution of the layer closely resembles 2D reconnection (panel c in Fig. 3). Near the shock front, multiple identical islands form, and this phenomenon corresponds to the linear phase of the tearing instability (at r/r_{min} ∼ 3 in Fig. 3, top of panel c). Islands are pushed toward one another by the reconnection outflows and merge, leaving behind a chain of bigger plasmoids separated by secondary layers. These shorter layers become themselves unstable to the secondary tearing mode, and more, but smaller, plasmoids are created (r/r_{min} ≳ 6; Fig. 2 and panel b of Fig. 3), which eventually merge with the big islands. The overall dynamics of the shock is robust against magnetization. The main differences are the plasma density contrast between the islands, the secondary layers, and the ambient downstream flow that increases with increasing magnetization, meaning that islands are much more visible as structures from the background flow at high magnetization. The strength of the velocity shear across the layer also increases with increasing magnetization, such that the KelvinHelmholtz instability is most severe for σ = 4.5 and 15.
Fig. 3.
Twodimensional slices through the 3D plasma density map shown in Fig. 2 (σ = 4.5, t = 8.44 r_{min}/c). Panel a: poloidal cut at ϕ = 0°. White solid lines with arrows represent the plasma bulk velocity streamlines. Panel b: equatorial cut (θ = π/2). Panel c: cuts through the θϕ plane at r/r_{min} = 3.65, 6.13, and 9.17, zoomedin on the current layer. These radii are indicated by dashed black lines in panels a and b. Panel c: white contours with arrows show the magnetic field lines. 
3.2. Hierarchical merging and formation of giant plasmoids
The merging process inevitably leads to a decrease in the number of plasmoids and to an increase in their sizes (e.g., Sironi et al. 2016). To assess whether major plasmoids can grow to macroscopic sizes at the physical scales of pulsar wind nebulae, and therefore examine whether they are viable candidates for explaining the knots, we need to characterize this evolution accurately in the simulations. Figure 4 shows the number of plasmoids as a function of time and at a fixed radius (r/r_{min} = 8.35). The center of a plasmoid is effectively determined as a local maximum in the magnetic flux function defined in the θϕ plane (local minima being the centers of secondary current layers, i.e., X points). In all runs, the initial number of plasmoids is about 30. Assuming that the layer thickness, δ, is on the order of the electron skindepth scale, d_{e}, as it naturally occurs in a collisionless system, the fastest growing tearing mode has a wavelength of (Zelenyi & Krasnoselskikh 1979; Zenitani & Hoshino 2007; Cerutti et al. 2014). The ratio of the box size to the tearing wavelength predicts that L_{ϕ}/λ_{T} ≈ 30 identical islands should form in the linear phase, in agreement with our results. This number should also be independent of σ, as reported here, since the box size is a fixed number of electron skin depths. The number of islands then decreases with time and reaches about a dozen toward the end of the simulation. The final number decreases with increasing magnetization from about 16 islands for σ = 0.45 down to about 8 islands for σ = 15.
Fig. 4.
Time evolution of the number of plasmoids, N_{is}, measured at a fixed radius, r = 8.35 r_{min}, for each value of σ (histograms). The dashed lines represent the best analytical fits to the hierarchical merging model (Zhou et al. 2019). 
In the Crab Nebula, the scale separation between the skin depth and the shock size is much larger than what is achievable with current numerical resources. Assuming a monoenergetic wind of pairs with a Lorentz factor of Γ ∼ 10^{6}, the number density of particles at the shock is of order cm^{−3}, where L_{0} ≈ 5 × 10^{38} erg s^{−1} is the Crab pulsar spindown power, R_{TS} ∼ 0.1 parsec is the radius of the termination shock (i.e., the inner ring), and κ = 3 is the shock compression ratio assuming a strong shock. This gives a skindepth scale of order d_{Crab} ∼ 2 × 10^{12} cm. Hence,
should be created initially in the Crab Nebula, that is to say, it is about a thousand times larger than what is achieved by 3D PIC simulations (∼4 × 30 = 120 islands, if considering a full 2π domain size). We now address the crucial question of how many islands will remain at the end of the merging process.
To answer this question and fill the gap between our simulation results and the Crab Nebula scale, we made use of the hierarchical model proposed by Zhou et al. (2019, 2020, 2021), which we reproduce here for completeness. This model assumes that the merging mechanism happens in discrete steps, during which identical, equally spaced islands merge by pairs. It also neglects the generation of secondary plasmoids as reconnection proceeds. If the initial number of plasmoids is N_{0}, then their number at step n is
The duration of each merger is controlled by the reconnection timescale given by
at step n, where
is the half separation between the plasmoids and β_{rec} is the dimensionless reconnection rate (i.e., the velocity at which reconnection operates, normalized by the speed of light). The initial half distance between islands is d_{0} = L_{ϕ}/2N_{0}, so the fiducial reconnection timescale is
The total time to reach the nth generation is then
Putting everything together and dropping the subscript n in the continuous limit of multiple successive merging events, the hierarchical model predicts that the number of islands should decay as
Figure 4 shows that the hierarchical model describes the overall evolution of islands in the simulation well. Leaving the reconnection rate and the initial number of islands as free parameters, we obtain the bestfit solutions, shown in Fig. 4 by dashed lines. The inferred reconnection rate measured at this particular radius ranges from β_{rec} ≈ 0.05 for σ = 0.45 up to β_{rec} ≈ 0.15 for σ = 15. These values fall well within the usual rates reported in the literature (e.g., Werner et al. 2018).
To verify that we are indeed measuring the reconnection rate with this model, we performed a more direct measurement using two diagnostics. The first method consists in measuring the asymptotic inflow velocity toward the major X point, V_{in} = β_{rec}c, at multiple radii. Once the major X point is identified (using the global minimum of the magnetic flux function), we fit the plasma bulk velocity profile across the layer as V_{θ} = V_{in}tanh((θ−π/2)/δ_{θ}), where δ_{θ} ≈ δ/r is the angular layer thickness. The second method consists in measuring the ratio of the reconnection electric field inside the layer, here the radial component E_{rec}, with the upstream reconnecting magnetic field, here the toroidal component assuming the same profile as for the velocity, B_{ϕ} = B_{up}tanh((θ−π/2)/δ_{θ}), so that β_{rec} = E_{rec}/B_{up}. Both methods give similar results, but the rate measured with the fields is cleaner (less sensitive to the particle noise) and is reported in Fig. 5. It ranges from about β_{rec} = 0.03 for σ = 0.45 up to β_{rec} = 0.12 for σ = 15, in good agreement with the rates inferred from the hierarchical model. The increase in the rate with increasing magnetization was also reported in Werner et al. (2018). Beyond σ ≳ 10, the rate should reach a saturation, (see Fig. 9 in Werner et al. 2018). We note that if the rate is defined as the reconnection speed normalized by the Alfvén velocity, , it would be nearly independent of σ.
Fig. 5.
Dimensionless reconnection rate, β_{rec}, as a function of σ (blue dots). This rate is measured at multiple locations in the simulation domain, and the error bars show the standard deviation from the mean value. The final number of giant plasmoids is inferred from the hierarchical merging model at the equatorial plane of pulsar wind nebulae, assuming an escape time of t_{esc} = R_{TS}/c (red dots). 
The reconnection rate having been characterized, the last parameter to fix before deducing the final number of giant plasmoids is the total duration of the reconnection process. If the lifetime of the sheet were unlimited, reconnection would saturate and all of the plasmoids would ultimately end up in a single giant island, as seen for instance in simulations with periodic boundary conditions. In pulsar wind nebulae, the lifetime of the layer will be limited by the largescale dynamics of the postshock flow. A naive estimate can be given by the advection time of the flow, moving at a speed of order 0.5c in the Crab Nebula (Hester et al. 2002). As mentioned earlier, another possible mechanism is the KelvinHelmholtz instability that is due to the strong velocity shear forming across the layer, which could disrupt the largescale coherent structure of the sheet, converting it into a turbulent flow (Cerutti & Giacinti 2020). The timescale is on the order of ∼0.1 R_{TS}/ΔV ≲ R_{TS}/c (Miura & Pritchett 1982). Another viable source of disruption is the driftkink instability (Begelman 1998; Mizuno et al. 2011; Porth et al. 2014). In this case, the growing time of this instability is also on the order of the lightcrossing time of the shock size. All in all, it is an interplay of at least these three effects that will limit the duration of the merging process to about one lightcrossing time of the shock size. Thus, defining the escape time as t_{esc} = R_{TS}/c, and using the hierarchical model (Eq. (13)), the final number of islands is
Given that β_{rec}R_{TS} ≫ d_{0}, and noticing that N_{0}d_{0} = L_{ϕ}/2 = πR_{TS}, this yields
where β_{0.1, rec} = β_{rec}/0.1, meaning that the final number of giant plasmoids is solely governed by the reconnection rate, and thus it is remarkably independent of the system size, but also of the initial number of plasmoids (as long as β_{rec}R_{TS} ≫ d_{0}). Figure 5 reports the corresponding number of islands as a function of σ derived from Eq. (15). The reconnection rate is reported in the same figure. Although the simulations presented here cannot capture the fullscale separation, this important conclusion implies that they are still relevant to capturing the highend part of the distribution and therefore predicting the correct number of giant plasmoids.
To further strengthen this point, we performed a standard 2D Cartesian PIC simulation of a Harris reconnection current layer of total length L_{x} for σ = 5 and with no guide field (i.e., the magnetic component perpendicular to the reconnection plane). Thanks to the largerscale separation that is achievable in 2D, the initial number of islands is N_{0} ≈ 175. This number decays with time in good agreement with the hierarchical model (Fig. 6). The inferred rate is β_{rec} ≈ 0.075, which is compatible with the rate measured in the σ = 4.5 3D run. The number of plasmoids after one lightcrossing time of the layer, L_{x}/c, is , which is comparable to the 3D shock simulation despite their initial number differing by a factor of ≈6.
Fig. 6.
Time evolution of the number of plasmoids measured in a 2D Cartesian PIC simulation initialized with a Harris layer of total length L_{x} (solid blue line), and comparison with the hierarchical model with a bestfit reconnection rate β_{rec} = 0.075 (dashed red line). Small secondary plasmoids forming in the late stages of reconnection between giant plasmoids are not considered in this measurement. 
3.3. The effect of expansion
In the above discussion, the effect of expansion has not been explicitly taken into account. Yet, measurements of proper motions in the Crab Nebula show that the downstream flow is expanding at mildly relativistic speeds, as also expected from the shock jump conditions (Kennel & Coroniti 1984b). Simulations presented here do not capture this effect well, despite the spherical geometry, because of the static reflecting outer boundary condition chosen for numerical convenience. Flux ropes are produced nearly at rest downstream of the shock as the latter recedes inward, and therefore they do not experience significant expansion over time.
To assess the effect of expansion, we considered a collection of N_{is} equally spaced islands comoving in the radial direction at a constant velocity, v_{exp}. Each closest pair of islands is separated by an angle α = 2π/N_{is}. In the comoving frame, islands only have a motion along the azimuthal direction (i.e., perpendicular to the flow); this is set by the reconnection speed, which in the lab frame is β_{rec}c/Γ_{exp}, where Γ_{exp} is the Lorentz factor of the expanding flow. Islands will be able to merge if the speed at which they are pushed away from each other by expansion (i.e., αv_{exp}) is smaller than the reconnection speed,
Assuming a mildly relativistic expanding flow, Γ_{exp} ≈ 1, which is applicable to the nebula, the number of islands should not be smaller than
Taking v_{exp}/c = 0.5 for the Crab gives islands, which is the same result we get from Eq. (15). In other words, even if the current layer had an infinite lifetime, it would never be able to reach a fully saturated state where only a single giant plasmoid would remain, in contrast to a static and finite system. Therefore, the spherical expansion of the downstream flow truncates the tail of the distribution predicted by the hierarchical model. This causality argument gives an additional and robust constraint in favor of the expected number of giant flux ropes formulated in Eq. (15), which was obtained by limiting the lifetime of the hierarchical merging process to the lightcrossing time of the shock.
3.4. Synchrotron hotspots
Giant magnetic islands are filled with all of the particles energized by reconnection up to a particle Lorentz factor γ ∼ Γ_{0}σ, and hence they are intense sources of nonthermal synchrotron radiation and appear as bright knots or, more precisely, as filaments downstream of the shock front along the equatorial plane (Fig. 7). Although synchrotron emission is highly polarized and the magnetic field in islands well ordered, the polarization signal averaged over the surface of the island for a looplike magnetic structure should be canceled out, leaving a mostly unpolarized signal. This would provide a natural explanation as to why knots are not visible in polarized optical light (Hester 2008), in contrast to the inner knot. The observed flux variability of the knots can also be accounted for by the wellknown dynamical nature of the flux rope network. With a few dozen knots in the inner ring in the Crab, this model thus favors a magnetized σ ≳ 5 pulsar wind nebula (Lyutikov 2010; Porth et al. 2014; Cerutti & Giacinti 2020). This conclusion is independently supported by the hard spectral index measured along the inner ring in Xrays (amongst the hardest values found in the Crab Nebula; Mori et al. 2004). With a reported Xray spectral index of α ∼ 0.8 in the ring (F ∝ ν^{−α}, where ν is the photon frequency and F the flux of photons per second), this yields an injected particle distribution of dN/dγ ∝ γ^{−1.6}. Such a hard particle distribution can be achieved by reconnection for σ ≳ 5 (Guo et al. 2014; Werner et al. 2016).
Fig. 7.
Zoomedin view on the current sheet sliced at r/r_{min} = 6.13 for the σ = 4.5 run. Top panel: square of the mean particle Lorentz factor, ⟨γ⟩^{2}. Middle panel: square of the total magnetic field strength, B^{2}. Bottom panel: proxy for the total synchrotron power, nB^{2}⟨γ⟩^{2}. White contours show the magnetic field lines projected on this plane. 
4. Conclusion
Pulsar wind nebulae decidedly offer a formidable and clean environment for studying relativistic magnetic reconnection in nature. In this work we have shown that a largescale ringlike current sheet forms and reconnects in the equatorial plane downstream of pulsar wind termination shocks. Reconnection operates in a nearly ideal textbooklike configuration: a thin, flat layer with periodic boundary conditions (along the azimuthal direction) and no guide field that is continuously reforming behind the shock. It is fed by the highlatitude magnetized downstream flow. Reconnection proceeds in the relativistic plasmoiddominated regime, leading to an efficient fragmentation of the current layer, dissipation of the toroidal magnetic field, and nonthermal particle acceleration. A large number of small, nearly identical plasmoids are first created in the immediate vicinity of the shock front. They merge in a hierarchical manner, leaving behind a few macroscopic structures that contain all the energetic particles processed by reconnection. Assuming that reconnection is unperturbed for at most a lightcrossing time of the shock radius – after which KelvinHelmholtz and kink instabilities will mix and disrupt the layer, leading to a turbulent state – the final number of giant plasmoids is solely determined by the inverse of the dimensionless reconnection rate, . A similar conclusion can be drawn from a causality argument: plasmoids can merge until the expansion of the downstream flow freezes their evolution. The minimum final number of large flux ropes depends on the ratio of the expansion to the reconnection speeds.
This remarkable result predicts that only magnetized nebulae such as the Crab (σ ≳ 5) should present discernible large islands (a few dozen at most); otherwise, they would appear as a continuous ring rather than a collection of welldefined discrete entities. We argue that the presence of these giant plasmoids (i.e., flux ropes) is a robust feature of pulsar wind termination shocks and provides a natural explanation as to the mysterious knots observed along the inner ring in the Crab Nebula (which, in fact, should appear more like filaments rather than circular structures). Hence, counting the knots gives a direct measure of the reconnection rate and average plasma magnetization injected by the wind into the nebula. In this interpretation, the presence of knots represents a real smoking gun of the magnetic reconnection process, which could well be involved in the Crab gammaray flare activity in the reconnection scenario (Uzdensky et al. 2011; Cerutti et al. 2013). Another important implication of this work is that the shock front should not coincide with the inner ring; instead, it should be located closer inward and should be underluminous in comparison with the giant plasmoids forming further downstream.
To make more progress and capture these knots in a fully global setup along with the other features of the nebula (such as the jet, the wisps, and the inner knot, which are missing here), highresolution 3D resistive MHD simulations with test particles could be used to provide an adequate answer before this becomes feasible with full PIC simulations. Last, we would like to point out that the same argument as the one developed in this work can be applied to other astrophysical environments where a largescale current sheet forms, such as the striped wind current sheet in pulsars (e.g., Cerutti et al. 2020) and the ergospheric current sheet in black hole magnetospheres (Parfrey et al. 2019; Crinquand et al. 2021), with perhaps significant observational consequences to come.
Assuming that χ ≤ π/2. If χ > π/2, then B_{ϕ} → −B_{ϕ} in Eq. (1).
Acknowledgments
The authors would like to thank Benjamin Crinquand, Lorenzo Sironi, and the referee Andrei Bykov for valuable comments regarding this work, and Gilles Henri for drawing our attention to the effect of expansion. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 863412). Computing resources were provided by TGCC and CINES under the allocation A0090407669 made by GENCI.
References
 Begelman, M. C. 1998, ApJ, 493, 291 [Google Scholar]
 Begelman, M. C. 1999, ApJ, 512, 755 [Google Scholar]
 Bogovalov, S. V. 1999, A&A, 349, 1017 [Google Scholar]
 Bogovalov, S. V., & Khangoulian, D. V. 2002, MNRAS, 336, L53 [Google Scholar]
 Camus, N. F., Komissarov, S. S., Bucciantini, N., & Hughes, P. A. 2009, MNRAS, 400, 1241 [Google Scholar]
 Cerutti, B., & Giacinti, G. 2020, A&A, 642, A123 [Google Scholar]
 Cerutti, B., & Philippov, A. A. 2017, A&A, 607, A134 [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147 [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2014, ApJ, 782, 104 [Google Scholar]
 Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401 [Google Scholar]
 Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204 [Google Scholar]
 Coroniti, F. V. 1990, ApJ, 349, 538 [Google Scholar]
 Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., & Philippov, A. 2021, A&A, 650, A163 [Google Scholar]
 Del Zanna, L., Amato, E., & Bucciantini, N. 2004, A&A, 421, 1063 [Google Scholar]
 Durant, M., Kargaltsev, O., Pavlov, G. G., Kropotina, J., & Levenfish, K. 2013, ApJ, 763, 72 [Google Scholar]
 Giacinti, G., & Kirk, J. G. 2018, ApJ, 863, 18 [Google Scholar]
 Guo, F., Li, H., Daughton, W., & Liu, Y.H. 2014, Phys. Rev. Lett., 113, 155005 [Google Scholar]
 Helfand, D. J., Gotthelf, E. V., & Halpern, J. P. 2001, ApJ, 556, 380 [Google Scholar]
 Hester, J. J. 2008, ARA&A, 46, 127 [Google Scholar]
 Hester, J. J., Scowen, P. A., Sankrit, R., et al. 1995, ApJ, 448, 240 [Google Scholar]
 Hester, J. J., Mori, K., Burrows, D., et al. 2002, ApJ, 577, L49 [Google Scholar]
 Kagan, D., Milosavljević, M., & Spitkovsky, A. 2013, ApJ, 774, 41 [Google Scholar]
 Kargaltsev, O., Cerutti, B., Lyubarsky, Y., & Striani, E. 2015, Space Sci. Rev., 191, 391 [Google Scholar]
 Kennel, C. F., & Coroniti, F. V. 1984a, ApJ, 283, 710 [Google Scholar]
 Kennel, C. F., & Coroniti, F. V. 1984b, ApJ, 283, 694 [Google Scholar]
 Komissarov, S. S., & Lyubarsky, Y. E. 2003, MNRAS, 344, L93 [Google Scholar]
 Komissarov, S., & Lyubarsky, Y. 2004, Ap&SS, 293, 107 [Google Scholar]
 Komissarov, S. S., & Lyutikov, M. 2011, MNRAS, 414, 2017 [Google Scholar]
 Loureiro, N. F., Schekochihin, A. A., & Cowley, S. C. 2007, Phys. Plasmas, 14, 100703 [Google Scholar]
 Loureiro, N. F., Samtaney, R., Schekochihin, A. A., & Uzdensky, D. A. 2012, Phys. Plasmas, 19, 042303 [Google Scholar]
 Lyubarsky, Y. E. 2002, MNRAS, 329, L34 [Google Scholar]
 Lyutikov, M. 2010, MNRAS, 405, 1809 [NASA ADS] [Google Scholar]
 Michel, F. C. 1973, ApJ, 180, L133 [Google Scholar]
 Michel, F. C. 1994, ApJ, 431, 397 [Google Scholar]
 Miura, A., & Pritchett, P. L. 1982, J. Geophys. Res., 87, 7431 [Google Scholar]
 Mizuno, Y., Lyubarsky, Y., Nishikawa, K.I., & Hardee, P. E. 2011, ApJ, 728, 90 [Google Scholar]
 Moran, P., Shearer, A., Mignani, R. P., et al. 2013, MNRAS, 433, 2564 [Google Scholar]
 Mori, K., Burrows, D. N., Hester, J. J., et al. 2004, ApJ, 609, 186 [NASA ADS] [CrossRef] [Google Scholar]
 Nalewajko, K., Uzdensky, D. A., Cerutti, B., Werner, G. R., & Begelman, M. C. 2015, ApJ, 815, 101 [Google Scholar]
 Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101 [Google Scholar]
 Pavlov, G. G., Zavlin, V. E., Sanwal, D., Burwitz, V., & Garmire, G. P. 2001a, ApJ, 552, L129 [Google Scholar]
 Pavlov, G. G., Kargaltsev, O. Y., Sanwal, D., & Garmire, G. P. 2001b, ApJ, 554, L189 [Google Scholar]
 Porth, O., Komissarov, S. S., & Keppens, R. 2013, MNRAS, 431, L48 [Google Scholar]
 Porth, O., Komissarov, S. S., & Keppens, R. 2014, MNRAS, 438, 278 [Google Scholar]
 Rees, M. J., & Gunn, J. E. 1974, MNRAS, 167, 1 [Google Scholar]
 Scargle, J. D. 1969, ApJ, 156, 401 [Google Scholar]
 Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54 [Google Scholar]
 Sironi, L., Giannios, D., & Petropoulou, M. 2016, MNRAS, 462, 48 [Google Scholar]
 Uzdensky, D. A., Loureiro, N. F., & Schekochihin, A. A. 2010, Phys. Rev. Lett., 105, 235002 [Google Scholar]
 Uzdensky, D. A., Cerutti, B., & Begelman, M. C. 2011, ApJ, 737, L40 [Google Scholar]
 Weisskopf, M. C., Hester, J. J., Tennant, A. F., et al. 2000, ApJ, 536, L81 [Google Scholar]
 Werner, G. R., & Uzdensky, D. A. 2017, ApJ, 843, L27 [Google Scholar]
 Werner, G. R., & Uzdensky, D. A. 2021, ArXiv eprints [arXiv:2106.02790] [Google Scholar]
 Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8 [Google Scholar]
 Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840 [Google Scholar]
 Zelenyi, L. M., & Krasnoselskikh, V. V. 1979, Sov. Astron., 23, 460 [Google Scholar]
 Zenitani, S., & Hoshino, M. 2007, ApJ, 670, 702 [Google Scholar]
 Zhou, M., Bhat, P., Loureiro, N. F., & Uzdensky, D. A. 2019, Phys. Rev. Res., 1, 012004 [Google Scholar]
 Zhou, M., Loureiro, N. F., & Uzdensky, D. A. 2020, J. Plasma Phys., 86, 535860401 [Google Scholar]
 Zhou, M., Wu, D. H., Loureiro, N. F., & Uzdensky, D. A. 2021, J. Plasma Phys., submitted [arXiv:2104.13757] [Google Scholar]
Movie
Movie 1 associated with Fig. 8 (density_pwn3d_knots) (Access here)
All Figures
Fig. 1.
Diagram schematically representing the morphological features observed in the Crab Nebula in the optical and Xray bands. The MHD model captures all these features well, except the knots along the inner ring (in red), whose origin is unknown. Solving this mystery is the main focus of this work. 

In the text 
Fig. 2.
Threedimensional rendering of the final state of the plasma density (r/r_{min})^{2}n/n_{0} for the σ = 4.5 run. For the sake of clarity, axes are not shown; see Fig. 3 for 2D slices with axes and length scales. 

In the text 
Fig. 3.
Twodimensional slices through the 3D plasma density map shown in Fig. 2 (σ = 4.5, t = 8.44 r_{min}/c). Panel a: poloidal cut at ϕ = 0°. White solid lines with arrows represent the plasma bulk velocity streamlines. Panel b: equatorial cut (θ = π/2). Panel c: cuts through the θϕ plane at r/r_{min} = 3.65, 6.13, and 9.17, zoomedin on the current layer. These radii are indicated by dashed black lines in panels a and b. Panel c: white contours with arrows show the magnetic field lines. 

In the text 
Fig. 4.
Time evolution of the number of plasmoids, N_{is}, measured at a fixed radius, r = 8.35 r_{min}, for each value of σ (histograms). The dashed lines represent the best analytical fits to the hierarchical merging model (Zhou et al. 2019). 

In the text 
Fig. 5.
Dimensionless reconnection rate, β_{rec}, as a function of σ (blue dots). This rate is measured at multiple locations in the simulation domain, and the error bars show the standard deviation from the mean value. The final number of giant plasmoids is inferred from the hierarchical merging model at the equatorial plane of pulsar wind nebulae, assuming an escape time of t_{esc} = R_{TS}/c (red dots). 

In the text 
Fig. 6.
Time evolution of the number of plasmoids measured in a 2D Cartesian PIC simulation initialized with a Harris layer of total length L_{x} (solid blue line), and comparison with the hierarchical model with a bestfit reconnection rate β_{rec} = 0.075 (dashed red line). Small secondary plasmoids forming in the late stages of reconnection between giant plasmoids are not considered in this measurement. 

In the text 
Fig. 7.
Zoomedin view on the current sheet sliced at r/r_{min} = 6.13 for the σ = 4.5 run. Top panel: square of the mean particle Lorentz factor, ⟨γ⟩^{2}. Middle panel: square of the total magnetic field strength, B^{2}. Bottom panel: proxy for the total synchrotron power, nB^{2}⟨γ⟩^{2}. White contours show the magnetic field lines projected on this plane. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.