© B. Cerutti and G. Giacinti 2021

1. Introduction

High-resolution X-ray 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 jet-like structure aligned along the pulsar rotation axis and of a bright torus-like 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, arc-like 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 X-rays. 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 X-ray 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 jet-torus 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 electron-positron pairs it contains; this gives rise to bright synchrotron radiation (Rees & Gunn 1974; Kennel & Coroniti 1984a,b). Two-dimensional (2D) axisymmetric magnetohydrodynamic (MHD) simulations have successfully reproduced the observed jet-torus 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 Doppler-boosted emission of the termination shock pointing toward the observer (Komissarov & Lyubarsky 2004; Komissarov & Lyutikov 2011). Three-dimensional (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 particle-in-cell (PIC) simulations, we have shown in a previous study that the anisotropy of the wind naturally results in the formation of a large-scale, 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 over-densities 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 semi-global 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 electron-positron pairs propagating outward with a bulk velocity of V₀/c = 0.99, or a bulk Lorentz factor of Γ₀ ≈ 7. This choice satisfies both the need to model an ultra-relativistic flow, as found in pulsar winds (i.e., Γ₀ ≫ 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 counter-propagating 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¹ (Cerutti & Giacinti 2020)

$$\begin{array}{c}\hfill B_{\varphi }=B_0\left(\frac{r_{\min }}{r}\right)tanh\left(\frac{\theta -\pi /2}{\chi }\right)sin\theta ,\end{array}$$(1)

B_r = 0 and B_θ = 0. The sin θ/r term is borrowed from the split-monopole model (Michel 1973; Bogovalov 1999), a good description of the asymptotic dissipation-less 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 non-vanishing component, E_θ = V₀B_ϕ/c.

In the steady state, the wind must carry the current density J = c/(4π)∇ × B. Using Eq. (1) yields

$$\begin{array}{c}\hfill J_r=\frac{c{B}_0}{4\pi }\frac{r_{\min }}{r^2}\{\frac{sin\theta }{\chi {cosh}^2\left(\frac{\theta -\pi /2}{\chi }\right)}+2tanh\left(\frac{\theta -\pi /2}{\chi }\right)cos\theta \},\end{array}$$(2)

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,

$$\begin{array}{c}\hfill \rho =\frac{1}{4\pi }\mathbf{\nabla }·\mathit{E}=\frac{V_0}{c^2}{J}_r={\rho }_{+}+{\rho }_{-},\end{array}$$(3)

where

$$\begin{array}{cc}& {\rho }_{+}=\frac{V_0{B}_0}{4\pi c}\frac{r_{\min }}{r^2}\frac{sin\theta }{\chi {cosh}^2\left(\frac{\theta -\pi /2}{\chi }\right)}\ge 0,\hfill \end{array}$$(4)

$$\begin{array}{cc}& {\rho }_{-}=\frac{V_0{B}_0}{2\pi c}\frac{r_{\min }}{r^2}tanh\left(\frac{\theta -\pi /2}{\chi }\right)cos\theta \le 0.\hfill \end{array}$$(5)

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₀. The pulsar wind must also fulfill quasi-neutrality (i.e., |ρ|≪ne). This condition is easily achieved by adding a uniform plasma density of n_w = n₀(r_min/r)² to both the electrons and positrons, that is drifting outward at V₀ 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 Courant-Friedrich-Lewy stability condition.

The fiducial plasma density, n₀, and magnetic field, B₀, set the upstream plasma magnetization parameter, the only free parameter in this study, which is defined as

$$\begin{array}{c}\hfill {\sigma }_0=\frac{B_0^2}{4\pi {\mathrm{\Gamma }}_0{n}_0{m}_{\mathrm{e}}{c}^2},\end{array}$$(6)

where m_e is the electron mass. In this work, we performed runs with σ₀ = 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σ₀. 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, $d_{\text{e}}^0={({\mathrm{\Gamma }}_0{m}_{\text{e}}{c}^2/8\pi n_0{e}^2)}^{1/2}$. The spatial resolution is about five cells per electron skin depth in all directions, therefore fitting about $L_{\varphi }=\pi r_{\text{min}}/2=320d_{\text{e}}^0$ 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 large-scale ring-like 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, skin-depth-scale 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 large-scale circulation pattern that takes place in the downstream medium and brings high-latitude 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 Kelvin-Helmholtz 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 Kelvin-Helmholtz instability may ultimately dominate and lead to the disruption and transformation of this large-scale 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. Three-dimensional rendering of the final state of the plasma density (r/rmin)2n/n0 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 space-time 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 Kelvin-Helmholtz instability is most severe for σ = 4.5 and 15.

Fig. 3. Two-dimensional slices through the 3D plasma density map shown in Fig. 2 (σ = 4.5, t = 8.44 rmin/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/rmin = 3.65, 6.13, and 9.17, zoomed-in 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 skin-depth scale, d_e, as it naturally occurs in a collisionless system, the fastest growing tearing mode has a wavelength of ${\lambda }_{\mathrm{T}}=2\pi \sqrt{3}{d}_{\mathrm{e}}$ (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, Nis, measured at a fixed radius, r = 8.35 rmin, 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 mono-energetic wind of pairs with a Lorentz factor of Γ ∼ 10⁶, the number density of particles at the shock is of order $n_{\text{Crab}}\backslash \kappa L_0/4\pi R_{\text{TS}}^2\mathrm{\Gamma }{m}_{\text{e}}{c}^3\approx 5\times {10}^{-8}$ cm⁻³, where L₀ ≈ 5 × 10³⁸ erg s⁻¹ is the Crab pulsar spin-down 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 skin-depth scale of order d_Crab ∼ 2 × 10¹² cm. Hence,

$$\begin{array}{c}\hfill N_{\mathrm{is},0}^{\mathrm{Crab}}\sim \frac{2\pi R_{\mathrm{TS}}}{{\lambda }_{\mathrm{T}}}=\frac{R_{\mathrm{TS}}}{\sqrt{3}{d}_{\mathrm{Crab}}}\sim {10}^5\mathrm{islands}\end{array}$$(7)

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₀, then their number at step n is

$$\begin{array}{c}\hfill N_n=2^{-n}{N}_0.\end{array}$$(8)

The duration of each merger is controlled by the reconnection timescale given by

$$\begin{array}{c}\hfill {\tau }_n=d_n/{\beta }_{\mathrm{rec}}c\end{array}$$(9)

at step n, where

$$\begin{array}{c}\hfill d_n=2^n{d}_0\end{array}$$(10)

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₀ = L_ϕ/2N₀, so the fiducial reconnection timescale is

$$\begin{array}{c}\hfill {\tau }_0=\frac{d_0}{{\beta }_{\mathrm{rec}}c}·\end{array}$$(11)

The total time to reach the nth generation is then

$$\begin{array}{c}\hfill t_n={\displaystyle \sum _{k=0}^{n-1}{\tau }_k={\tau }_0\sum _{k=0}^{n-1}{2}^k={\tau }_0(2^n-1).}\end{array}$$(12)

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

$$\begin{array}{c}\hfill N_{\mathrm{is}}(t)=\frac{N_0}{1+t/{\tau }_0}·\end{array}$$(13)

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 best-fit 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 = β_recc, 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_intanh((θ−π/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_uptanh((θ−π/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, ${\beta }_{\mathrm{rec}}^{\infty }\gtrsim 0.1$ (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, $V_{\mathrm{A}}=\sqrt{\sigma /(1+\sigma )}$, 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 tesc = RTS/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 large-scale dynamics of the post-shock 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 Kelvin-Helmholtz instability that is due to the strong velocity shear forming across the layer, which could disrupt the large-scale 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 drift-kink 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 light-crossing 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 light-crossing 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

$$\begin{array}{c}\hfill N_{\mathrm{is}}^{\mathrm{f}}=\frac{N_0}{1+t_{\mathrm{esc}}/{\tau }_0}=\frac{N_0}{1+{\beta }_{\mathrm{rec}}{R}_{\mathrm{TS}}/d_0}·\end{array}$$(14)

Given that β_recR_TS ≫ d₀, and noticing that N₀d₀ = L_ϕ/2 = πR_TS, this yields

$$\begin{array}{c}\hfill N_{\mathrm{is}}^{\mathrm{f}}=\pi {\beta }_{\mathrm{rec}}^{-1}\approx 30{\beta }_{0.1,\mathrm{rec}}^{-1}\mathrm{islands},\end{array}$$(15)

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 β_recR_TS ≫ d₀). 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 full-scale separation, this important conclusion implies that they are still relevant to capturing the high-end 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 larger-scale separation that is achievable in 2D, the initial number of islands is N₀ ≈ 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 light-crossing time of the layer, L_x/c, is $N_{\mathrm{is}}^{\mathrm{f}}\approx 20$, 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 Lx (solid blue line), and comparison with the hierarchical model with a best-fit 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 co-moving 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 co-moving 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 β_recc/Γ_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,

$$\begin{array}{c}\hfill \frac{2\pi }{N_{\mathrm{is}}}{v}_{\exp }\lesssim \frac{{\beta }_{\mathrm{rec}}c}{{\mathrm{\Gamma }}_{\exp }}·\end{array}$$(16)

Assuming a mildly relativistic expanding flow, Γ_exp ≈ 1, which is applicable to the nebula, the number of islands should not be smaller than

$$\begin{array}{c}\hfill N_{\mathrm{is}}^{\mathrm{f}}\gtrsim \frac{2\pi v_{\exp }}{{\beta }_{\mathrm{rec}}c}·\end{array}$$(17)

Taking v_exp/c = 0.5 for the Crab gives $N_{\mathrm{is}}^{\mathrm{f}}\gtrsim \pi {\beta }_{\mathrm{rec}}^{-1}$ 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 light-crossing 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 γ ∼ Γ₀σ, 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 loop-like 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 well-known 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 X-rays (amongst the hardest values found in the Crab Nebula; Mori et al. 2004). With a reported X-ray 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. Zoomed-in view on the current sheet sliced at r/rmin = 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, B2. Bottom panel: proxy for the total synchrotron power, nB2⟨γ⟩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 large-scale ring-like current sheet forms and reconnects in the equatorial plane downstream of pulsar wind termination shocks. Reconnection operates in a nearly ideal textbook-like 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 high-latitude magnetized downstream flow. Reconnection proceeds in the relativistic plasmoid-dominated 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 light-crossing time of the shock radius – after which Kelvin-Helmholtz 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, $N_{\mathrm{is}}^{\mathrm{f}}=\pi {\beta }_{\mathrm{rec}}^{-1}$. 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 well-defined 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 gamma-ray 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 under-luminous 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), high-resolution 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 large-scale 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.

Movie

Movie 1 associated with Fig. 8 (density_pwn3d_knots) Access here

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

All Figures

	Fig. 1. Diagram schematically representing the morphological features observed in the Crab Nebula in the optical and X-ray 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. Three-dimensional rendering of the final state of the plasma density (r/rmin)2n/n0 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. Two-dimensional slices through the 3D plasma density map shown in Fig. 2 (σ = 4.5, t = 8.44 rmin/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/rmin = 3.65, 6.13, and 9.17, zoomed-in 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, Nis, measured at a fixed radius, r = 8.35 rmin, 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 tesc = RTS/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 Lx (solid blue line), and comparison with the hierarchical model with a best-fit 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. Zoomed-in view on the current sheet sliced at r/rmin = 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, B2. Bottom panel: proxy for the total synchrotron power, nB2⟨γ⟩2. White contours show the magnetic field lines projected on this plane.
In the text