THREE-DIMENSIONAL RELATIVISTIC PAIR PLASMA RECONNECTION WITH RADIATIVE FEEDBACK IN THE CRAB NEBULA

All the simulations presented in this work were performed with Zeltron,⁶ a parallel 3D electromagnetic PIC code (Cerutti et al. 2013). Zeltron solves self-consistently Maxwell's equations using the Yee finite-difference time-domain (FDTD) algorithm (Yee 1966), and Newton's equation following the Boris FDTD algorithm (Birdsall & Langdon 2005). Unlike most PIC codes, Zeltron includes the effect of the radiation reaction force in Newton's equation (or the so-called Lorentz–Abraham–Dirac equation) induced by the emission of radiation by the particles (see also, e.g., Jaroschek & Hoshino 2009; Tamburini et al. 2010; Capdessus et al. 2012). In the ultra-relativistic regime, the radiation reaction force, g, is akin to a continuous friction force, proportional to the radiative power and opposite to the particle's direction of motion (e.g., Landau & Lifshitz 1975; Tamburini et al. 2010; Cerutti et al. 2012a). The expression of the radiation reaction force used in Zeltron is given by

$$\begin{equation} \mathbf {g}=-\frac{2}{3}r_{\rm e}^2\gamma \left[\left(\mathbf {E}+\frac{\mathbf {u}\times \mathbf {B}}{\gamma }\right)^2-\left(\frac{\mathbf {u}\cdot \mathbf {E}}{\gamma }\right)^2\right]\mathbf {u}, \end{equation} \tag{ 1 }$$

where r_e ≈ 2.82 × 10⁻¹³ cm is the classical radius of the electron, γ is the Lorentz factor of the particle, E and B are the electric and magnetic fields, and u = γv/c is the four-velocity divided by the speed of light. This formulation is valid if γB/B_QED ≪ 1, where B_QED = 4.4 × 10¹³ G is the quantum critical magnetic field. Because of the relativistic effects, the typical frequency of the expected radiation is ∼γ³ ≫ 1 times the relativistic cyclotron frequency. Hence, the radiation is not resolved by the grid and time step of the simulation. It must be calculated separately. Zeltron computes the emitted optically thin radiation (spectrum, and angular distributions) assuming it is pure synchrotron radiation. This is valid if the change of the particle energy is small, Δγ/γ ≪ 1, during the formation length of a synchrotron photon, given by the relativistic Larmor radius divided by γ. We checked a posteriori that this assumption is indeed correct. The code also models the inverse Compton drag force on the particle motion in an imposed photon field, but this effect is negligible in the context of the Crab flares (Cerutti et al. 2012a), hence this capability will not be utilized in the following. To perform the large 3D simulations presented in this paper, Zeltron ran on 97,200 cores on the Kraken supercomputer⁷ with nearly perfect scaling.

The simulation setup chosen here is almost identical to our previous 2D pair plasma reconnection simulations with radiation reaction force in Cerutti et al. (2013). The computational domain is a rectangular box of dimensions L_x, L_y and L_z, respectively along the x-, y- and z-directions, with periodic boundary conditions in all directions. We set up the simulation with two flat anti-parallel relativistic Harris current layers (Kirk & Skjæraasen 2003) in the xz-plane located at y = L_y/4 and y = 3L_y/4 (Figure 1). Having two current sheets is only a convenient numerical artifact that allows us to use periodic boundary conditions along the y-direction, but it does not have a physical meaning in the context of our model of the Crab flares where only one reconnection layer is involved. The electric current, J_z, flows in the ±z-directions, and is supported by electrons counter-streaming with positrons at a mildly relativistic drift velocity (relative to the speed of light) β_drift = 0.6. The plasma (electrons and positrons) is spatially distributed throughout the domain, with the following density profile

$$\begin{equation} n = \! \left\lbrace \begin{array}{@{}lll}n_0\left[\cosh \left(\frac{y-L_{\rm y}/4}{\delta }\right)\right]^{-2}+0.1n_0 &\mbox{if} & y<L_{\rm y}/2\\ n_0\left[\cosh \left(\frac{y-3L_{\rm y}/4}{\delta }\right)\right]^{-2}+0.1n_0 & \mbox{if} & y>L_{\rm y}/2\end{array} \right.. \end{equation} \tag{ 2 }$$

The first term is the density of the drifting pairs carrying the initial current, concentrated within the layer half-thickness δ = λ_D/β_drift, where λ_D is the relativistic Debye length (Kirk & Skjæraasen 2003). This population is modeled with a uniform and isotropic (in the co-moving frame) distribution of macro-particles with variable weights to account for the density profile and to decrease the numerical noise in low-density regions. The second term is a uniform and isotropic background pair plasma at rest in the laboratory frame with a density chosen to be 10 times lower than at the center of the layers (i.e., 0.1n₀). The drifting and the background particles are distributed in energy according to a relativistic Maxwellian with the same temperature θ₀ ≡ kT/m_ec² = 10⁸, where k is the Boltzmann constant and m_e is the rest mass of the electron. The temperature of the drifting particles is defined in the co-moving frame. This temperature models the ultra-relativistic plasma already present in the Crab Nebula, prior to reconnection, whose particles could have been accelerated at the wind termination shock or even by other reconnection events throughout the nebula. However, observations show that in reality the background plasma is distributed according to a broad and steep power law, extending roughly between γ_min = 10⁶ and γ_max = 10⁹ (responsible for the UV to 100 MeV synchrotron spectrum). This large dynamic range of particle energies translates directly into an equally large dynamic range of relativistic Larmor radii and hence of length scales that must be resolved in the simulation, which is beyond the reach of our numerical capabilities.

Zoom In Zoom Out Reset image size

The initial electromagnetic field configuration is

$$\begin{eqnarray} \mathbf {B}&=&\left\lbrace \begin{array}{@{}lll}-B_0\tanh \left(\frac{y-L_{\rm y}/4}{\delta }\right)\mathbf {e_{\rm x}}+\alpha B_0 \mathbf {e_{\rm z}} &\mbox{if} & y<L_{\rm y}/2 \\ B_0\tanh \left(\frac{y-3L_{\rm y}/4}{\delta }\right)\mathbf {e_{\rm x}}+\alpha B_0 \mathbf {e_{\rm z}} & \mbox{if} & y>L_{\rm y}/2\end{array} \right., \qquad \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\nonumber \\ \mathbf {E}&=&\mathbf {0}, \end{eqnarray} \tag{ 4 }$$

where e_x, e_z are unit vectors along the x- and z-directions. B₀ is the upstream reconnecting magnetic field and α is a dimensionless parameter of the simulation that quantifies the strength of the guide field component B_z in units of B₀ (Figure 1). Observations constrain the magnetic field in the emitting region to about 1mG, which is much higher than the expected average quiescent field of order 100–200 μG (Meyer et al. 2010). In this work, we choose B₀ = 5 mG to be consistent with our previous studies of the flares (Cerutti et al. 2012a, 2013). Hence, the energy scale at which the radiation reaction force equals the electric force, assuming that E = B₀ = 5 mG, is

$$\begin{equation} \gamma _{\rm rad}m_{\rm e}c^2=\sqrt{\frac{3e m^2_{\rm e} c^4}{2r_{\rm e}^2 B_0}}\approx 1.3\times 10^9 m_{\rm e}c^2, \end{equation} \tag{ 5 }$$

where e is the fundamental electric charge. Below, we express lengths in units of the typical initial Larmor radius of the particles in the simulations, i.e., ρ₀ = θ₀m_ec²/eB₀ ≈ 3.4 × 10¹³ cm. In all the simulations, the layer half-thickness is then δ/ρ₀ ≈ 2.7 and the relativistic collisionless electron skin-depth $d_{\rm e}\equiv \sqrt{\theta _0 m_{\rm e}c^2/4\pi n_0 e^2}\approx 1.8\rho _0$. Similarly, timescales are given in units of the gyration time of the bulk of the particles in the plasma, i.e., $\omega _0^{-1}\equiv \rho _0/c\approx 1140$ s. The initial distribution of fields and plasma results in a low plasma-β or high magnetization of the upstream plasma (i.e., outside the layers). Here, the magnetization parameter is $\sigma \equiv B_0^2/4\pi (0.1 n_0)\theta _0 m_{\rm e} c^2\approx 16$.

The system is initially set at an equilibrium, i.e., there is a force balance across the reconnection layers between the magnetic pressure and the drifting particle pressure. This equilibrium is unstable to two competing instabilities, namely the relativistic tearing and kink instabilities, as well as oblique modes that combine tearing and kink modes (Zenitani & Hoshino 2005, 2008; Daughton et al. 2011; Kagan et al. 2013). In contrast to Cerutti et al. (2013), we choose here not to apply any initial perturbation in order to avoid any artificial enhancement of one type of instability over the other. Instabilities are seeded with the numerical noise only. This choice has a direct computational cost because the lack of perturbation significantly delays the onset of reconnection (see Sections 3.1 and 4.1), but it enables a fair comparison between the growth rates of both instabilities (Sections 3.2 and 4.2). Another important consequence of this choice specific to this study is the significant radiative cooling of the particles before reconnection can accelerate them (Section 3.3).

In this work, we performed a series of 14 simulations. This set includes 12 2D simulations, and 2 3D simulations. The 2D simulations are designed to study the effect of the guide field strength, α = 0, 0.25, 0.5, 0.75 and 1 on the developments of instabilities (kink and tearing) and on particle acceleration/emission. To analyze the developments both instabilities separately, we follow the same approach as Zenitani & Hoshino (2007, 2008), i.e., we consider the dynamics of the Harris current layers in the xy-plane where the tearing modes alone develop similarly to our previous study in Cerutti et al. (2013), and in the yz-plane where the kink modes alone develop (the kink and the tearing modes are perpendicular to each other). The box is square of size L_x × L_y = (200ρ₀)² and L_y × L_z = (200ρ₀)² with 1440² cells and 16 particles per cell (all species together). The spatial resolution is ρ₀/Δx ≈ 7.2, where Δx is the grid spacing in the x-direction (Δx = Δy = Δz), which ensures the conservation of the total energy to within ≲ 1% error throughout the simulation. From this 2D scan, we identify the best/worst conditions for efficient particle acceleration and emission above the radiation reaction limit in 3D. The 3D box is cubical of size L_x × L_y × L_z = (200ρ₀)³ with 1440³ grid cells and 16 particles per cell (all species together). The simulation time step is set at 0.3 times the critical Courant–Friedrichs–Lewy time step, $\Delta t=0.3\Delta t_{\rm CFL}\approx 0.029\omega _0^{-1}$ in 2D and ${\approx }0.024\omega _0^{-1}$ in 3D, in order to maintain satisfactory total energy conservation in the presence of strong radiative damping. Table 1 enumerates all the simulations presented here.

Figure 2 shows the time evolution of the total plasma density and field lines at four characteristic stages of 2D magnetic reconnection in the xy-plane with α = 0 (run 2DXY0). Because there is no initial perturbation, the layers remain static until tω₀ ≈ 120 when the layer tears apart into about 7 plasmoids per layer separated by X-points where field lines reconnect. The noise of the macro-particles in the PIC code is sufficient to seed the tearing instability. The reconnection electric field E_z is maximum at X-points and is responsible for most of particle acceleration. The high magnetic tension of freshly reconnected field lines pushes the plasma toward the ±x-directions and drives the large scale reconnection outflow that forces magnetic islands to merge with each other. Reconnection proceeds until there is only one big island per layer remaining in the box. At the end of the simulation (tω₀ = 353), about 70% of the initial magnetic energy is dissipated in the form of particle kinetic energy. All the energy gained by the particles is then lost via the emission of synchrotron radiation. Adding a guide field does not suppress the tearing instability, but it creates a charge separation across the layer that induces a strong E_y electric field (see also Zenitani & Hoshino 2008; Cerutti et al. 2013).

Zoom In Zoom Out Reset image size

Figure 3 presents the time evolution of the 2D simulation in the yz-plane with no guide field (run 2DYZ0). The initial setup of fields and particles is identical to run 2DXY0, except that the reconnecting field (B_x) is now perpendicular to the simulation plane. Hence, reconnection and tearing modes cannot be captured by this simulation. Instead, we observe the development of the kink instability as early as tω₀ ≈ 100 in the form of a small sinusoidal deformation of the current sheets with respect to the initial layer mid-plane. The sinusoidal deformation proceeds along the z-direction, with the deformation amplitude in the y-direction increasing rapidly up to about a quarter of the simulation box size (about 50ρ₀). At this stage, the folded current layers are disrupted, leading to fast and efficient magnetic dissipation. About 55% of the total magnetic energy is dissipated by the end of the simulation. The guide field has a dramatic influence on the stability of the layers. The amplitude of the deformation as well as the magnetic energy dissipated decreases with increasing guide field. For α ≳ 0.75, the layers remain flat during the entire duration of the simulation and no magnetic energy is dissipated. In this case, the only noticeable time evolution is a slight decrease of the layer thickness due to synchrotron cooling. To maintain pressure balance across the layers with the unchanged upstream magnetic field, the layer must compress to compensate for the radiative energy losses (Uzdensky & McKinney 2011). We observed also a compression of the reconnection layer in the xy-reconnection simulations.

Zoom In Zoom Out Reset image size

To compare the relative strength of the tearing instability versus the kink instability, we perform a spectral analysis of the fastest growing modes that develop in the simulations. To study the kink instability, we do a fast Fourier transform (FFT) along the z-direction of the small variations of the reconnecting magnetic field in the bottom layer mid-plane, δB_x(z, t) = B_x(y = L_y/4, z, t) − B_x(y = L_y/4, z, 0), during the early phase of the 2D simulations in the yz-plane. For the tearing modes, we follow the same procedure for the fluctuations in the reconnected field along the x-direction, δB_y(x, t) = B_y(x, y = L_y/4, t) − B_x(y = L_y/4, z, 0), in the 2D simulations in the xy-plane.

We present in Figure 4 the time evolution of the fastest growing modes as well as the dispersion relations for the tearing and kink modes, with no guide field. In the linear regime (tω₀ ≲ 125), we infer the growth rates by fitting the amplitude of each mode with |FFT(δB_{x, y}/B₀)|∝exp (γ_gr(k)t), where γ_gr is the growth rate of the mode of wave-number k. We find that the fastest growing tearing mode is at k_xδ ≈ 0.58, which coincides with the analytical expectation of $k_{\rm x}\delta =1/\sqrt{3}$ (Zelenyi & Krasnoselskikh 1979) as found by Zenitani & Hoshino (2007). The wavelength of this mode is $L_{\rm x}/\lambda _{\rm x}=L_{\rm x}/2\pi \sqrt{(}3)\delta \approx 7$; this explains the number of plasmoids formed in the early stages of reconnection (see top right panel in Figure 2). The fastest growing kink mode has a wavelength L_z/λ_z ≈ 8 (or k_zδ ≈ 0.67) which is consistent with the deformation of the current layers observed in Figure 3, top right panel. The corresponding growth rate is $\gamma _{\rm KI}\omega _0^{-1}\approx 0.055$, which is comparable with the fastest tearing growth rate, $\gamma _{\rm TI}\omega _0^{-1}\approx 0.045$ (Figure 4, bottom panel). This is expected for an ultra-relativistic plasma (kT ≫ m_ec²) with a drift velocity β_drift = 0.6 (Zenitani & Hoshino 2007). It is worth noting that the dispersion relations for both the kink and the tearing instabilities are not sharply peaked around the fastest growing modes; a broad range of low-frequency modes is almost equally unstable (i.e., for 0 < k_{x, z}δ ≲ 1).

Zoom In Zoom Out Reset image size

In agreement with Zenitani & Hoshino (2008) and as pointed out in Section 3.1, we find that the kink instability depends sensitively on the guide field strength. Figure 5 shows that the fastest growth rate decreases rapidly between α = 0.25 and α = 0.75 from $\gamma _{\rm KI}\omega _0^{-1}\approx 0.055$ to undetectable levels. Thus, the guide field stabilizes the layer along the z-direction. In contrast, as mentioned earlier in Section 3.1, the tearing growth rate depends only mildly on α; we note a decrease from $\gamma _{\rm TI}\omega _0^{-1}\approx 0.045$ for α = 0 to $\gamma _{\rm TI}\omega _0^{-1}\approx 0.025$ for α = 1. For α ≳ 0.5, the tearing instability dominates over the kink.

Zoom In Zoom Out Reset image size

The critical quantities of interest here are the particle energy distributions, γ²dN/dγ, and the instantaneous optically thin synchrotron radiation spectral energy distribution (SED) emitted by the particles, $\nu F_{\nu }\equiv E^2 d{\rm N_{\rm ph}}/dt {dE}$, where E is the photon energy. Figure 6 shows the time evolution of the total particle spectra at different times with no guide field in the xy-plane (top panel) and in the yz-plane (bottom panel). In the early stage (tω₀ ≲ 132), both simulations are subject to pure synchrotron cooling (i.e., with no acceleration or heating) of the plasma that results in a decrease of the typical Lorentz factor of the particles from γ/γ_rad ≈ 0.3 at tω₀ = 0 to γ/γ_rad ≈ 0.08 at tω₀ = 132. The decrease of the mean particle energy within the layer explains the shrinking of the layer thickness described in Section 3.1.

Zoom In Zoom Out Reset image size

At tω₀ ≳ 132, the instabilities trigger magnetic dissipation and particles are energized, but the particle spectra differ significantly in both cases. In run 2DXY0, where the tearing instability drives reconnection, the particle spectrum extends to higher and higher energy with time until the end of the simulation, where the maximum energy reaches γ_max/γ_rad ≈ 2.5, i.e., well above the nominal radiation reaction limit. The spectrum above γ/γ_rad = 0.1 cannot be simply modeled with a single power law, but it is well contained between two steep power laws of index −2 and −3. We know from our previous study that the high-energy particles are accelerated via the reconnection electric field at X-points and follow relativistic Speiser orbits (Cerutti et al. 2013). The maximum energy is then given by the electric potential drop along the z-direction (neglecting radiative losses), i.e.,

$$\begin{equation} \gamma _{\rm max}\sim \frac{e E_{\rm z} L_{\rm x}}{m_{\rm e}c^2}=\frac{e \beta _{\rm rec}B_0 L_{\rm x}}{m_{\rm e}c^2}\approx 3\gamma _{\rm rad}, \end{equation} \tag{ 6 }$$

for a dimensionless reconnection rate β_rec ≈ 0.2. Particles above the radiation reaction limit (γ > γ_rad) account for about 5% of the total energy of the plasma at tω₀ = 318 (Figure 7, top panel), and are responsible for the emission of synchrotron radiation above 160 MeV. Figure 7 (bottom panel) shows the resulting isotropic synchrotron radiation SED at tω₀ = 318, where about 11% of the radiative power is >160 MeV. The SED peaks at E = 10 MeV and extends with a power law of index −0.42 up to about 300–400 MeV before cutting off exponentially.

Zoom In Zoom Out Reset image size

In contrast, in run 2DYZ0, where the kink instability drives the annihilation of the magnetic field, the particles are heated up to a typical energy γ/γ_rad ≈ 0.3. The particle spectrum is composed of a Maxwellian-like distribution on top of a cooled distribution of particles formed at tω₀ ≲ 132 (Figures 6 and 8). The mean energy of the hot particles corresponds to a nearly uniform redistribution of the total dissipated magnetic energy to kinetic energy of background particles, i.e.,

$$\begin{equation} \langle \gamma \rangle \sim 0.55\times \frac{B_0^2/8\pi }{(0.1n_0)m_{\rm e}c^2}=0.55\times \frac{\sigma \theta _0}{2}\approx 0.34\gamma _{\rm rad}, \end{equation} \tag{ 7 }$$

where the numerical factor 0.55 accounts for the fraction of the total magnetic energy dissipated at the end of the simulation. Hence, the development of the kink prevents the acceleration of particles above γ_rad and the emission of synchrotron photons above 160 MeV. Figure 8 (bottom panel) shows that the total synchrotron radiation SED peaks and cuts off at E = 10 MeV, far below the desired energies >160 MeV.

Zoom In Zoom Out Reset image size

Because a moderate guide field suppresses the effect of the kink instability, hence magnetic dissipation, the particles are not heated for α ≳ 0.5, and the initial spectrum continues cooling until the end of the yz-plane simulation where the particles radiate low-energy (∼1 MeV) synchrotron radiation (Figure 8). In the xy-plane reconnection simulations, the guide field tends to decrease the maximum energy of the particles and of the emitted radiation (Figure 7, see also Cerutti et al. 2013). The guide field deflects the particles outside the layer, reducing the time spent by the particle within the accelerating region.

Figure 9 (left panels) shows the time evolution of the plasma density⁸ in the zero-guide field simulation at tω₀ = 0, 173, 211 and 269. The initial stage where the layer remains apparently static lasts for about tω₀ = 144, i.e., half of the whole simulation time. At tω₀ ≳ 144, overdensities appear in the layers in the form of 7–8 tubes (flux ropes) elongated along the z-direction. These structures are generated by the tearing instability and are the 3D generalization of the magnetic islands observed in 2D reconnection. As the simulation proceeds into the nonlinear regime, the flux ropes merge with each other creating bigger ones, as magnetic islands do in 2D reconnection. However, in 3D this process does not happen at the same time everywhere along the z-direction, which results in the formation of a network of interconnected flux ropes at intermediate times (173 ≲ tω₀ ≲ 211).

Zoom In Zoom Out Reset image size

In parallel to this process, the kink instability deforms the two layers along the z-direction in the form of sine-like translation of the layers' mid-planes in the ±y-directions. During the most active period of reconnection (tω₀ ≳ 173), the kink instability takes over and eventually destroys the flux ropes formed by the tearing modes (see left bottom panel in Figure 9). Only a few coherent structures survive at the end of the simulation (tω₀ = 269). In particular the reconnection electric field, which is strongest along the X-lines between two flux ropes, loses its initial coherence. This results in efficient particle heating but poor particle acceleration (see below, Section 4.3). At the end of this run about 52% of the total magnetic energy is dissipated, although the simulation does not reach the fully saturated state.

The right panels in Figure 9 shows the time evolution of the plasma density for α = 0.5 guide field. One sees immediately that the guide field effectively suppresses the kink deformations of the layers in the ±y-directions, as expected from the 2D simulations in the yz-plane (See Section 3) and from Zenitani & Hoshino (2008). In contrast, the tearing instability seems undisturbed and breaks the layer into a network of eight flux tubes. Toward the end of the simulation, there are about three well-defined flux ropes containing almost all the plasma that went through reconnection. At this point in time, 20% of the total magnetic energy (i.e., including the reconnecting and the guide field energy) has dissipated, in agreement with the 2D run 2DXY050.

Following the analysis presented in Section 3.2, we perform a Fourier decomposition of the magnetic fluctuations in the bottom layer mid-plane, (x, y = L_y/4, z), to study the most unstable modes that develop in the 3D simulations. Figure 10 presents the growth rate of each modes in the (k_x × k_z)-plane estimated from the variations of B_x (left panel) and B_y (right panel), for α = 0. As pointed out in Section 3.2 and by Zenitani & Hoshino (2008) and Kagan et al. (2013), we find that the reconnecting field B_x effectively captures the kink-like modes along k_z whereas the reconnected field B_y is most sensitive to tearing-like modes along k_x. The dispersion relations show that pure kink (along k_z for k_x = 0) and pure tearing (along k_x for k_z = 0) modes grow at rates in very good agreement with the corresponding 2D simulations. With a growth rate γ_GR ≈ 0.06ω₀, the fastest growing mode in the simulation is a pure kink mode of wavenumber k_zδ ≈ 0.7, or L_z/λ_z ≈ 8 consistent with the deformation of the layer observed in the earlier stage of reconnection (Figure 9, left panels) and with the 2D run 2DYZ0. The fastest tearing mode has a growth rate γ_GR ≈ 0.045ω₀ at k_xδ ≈ 0.5 and generates the ≈7 initial flux ropes obtained in the simulation. The (k_x × k_z)-plane is also filled with oblique modes, i.e., waves with a non-zero k_x- and k_z-component, with growth rates comparable to the fastest tearing and kink modes. The existence of these modes is reflected by the flux ropes being slightly tilted in the xz-plane. Adding an α = 0.5 guide field decreases the amplitude of the low-frequency (k_zδ ≲ 1) growth rates of the kink modes (Figure 11). In particular, the growth rate of the fastest mode for α = 0, k_zδ = 0.7, decreases from 0.06ω₀ to 0.03ω₀. As a result, the fastest growing kink mode is now at k_zδ = 0.8 with a rate ≈0.04ω₀, while the fastest growing tearing modes is approximatively unchanged, in excellent agreement with the 2D runs (Figure 11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 presents the particle and photon energy distributions averaged over all directions at tω₀ = 265, with no guide field. The distributions are remarkably similar to the 2D run 2DYZ0 ones, and differ significantly from run 2DXY0. The high-energy part of the particle spectrum peaks at γ/γ_rad = 0.3 which is the signature of particle heating via magnetic dissipation rather than particle acceleration through tearing-dominated reconnection (Section 3.3, Equation (7)). We note that the spectrum extends to higher energy than the pure magnetic dissipation scenario, slightly above γ_rad, suggesting that there is a non-thermal component as well. On the contrary, in the α = 0.5 guide field case (run 3DG050, Figure 13), things are closer to the pure tearing reconnection case of run 2DXY050. The particle energy distribution is almost flat in the 0.04 ≲ γ/γ_rad ≲ 0.4 range, but barely reaches above γ_rad as in the zero-guide field case. Nevertheless, the synchrotron emission >160 MeV is more intense than in the zero-guide field case. Even though there is clear evidence for particle acceleration above the radiation reaction limit, the effect remains slightly weaker than in 2D with no guide field. A bigger box size would help to improve the significance of this result.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The angular distribution of the particles is of critical interest for determining the apparent isotropic radiation flux seen by a distant observer who probes one direction only. In 2D reconnection, we expect a pronounced beaming of the particles that increases rapidly with their energy (Cerutti et al. 2012b, 2013). We confirm here that this phenomenon exists also in 3D, even with a finite guide field. Figure 14 presents energy-resolved maps of the angular distribution of the positrons (left panels) and their optically thin synchrotron radiation (right panel) in run 3D050. The direction of motion of the particles is measured with two angles: the latitude, ϕ, varying between −90° and 90°, defined as

$$\begin{equation} \phi =\sin ^{-1}\left(\frac{u_{\rm y}}{\sqrt{u_{\rm x}^2+u_{\rm y}^2+u_{\rm z}^2}}\right), \end{equation} \tag{ 8 }$$

and the longitude, λ, defined between −180° and 180° given by

$$\begin{equation} \lambda = \left\lbrace \begin{array}{@{} lll}\cos ^{-1}\left(\frac{u_{\rm z}}{\sqrt{u_{\rm x}^2+u_{\rm z}^2}}\right) & \mbox{if} & \sin \lambda > 0 \\ {-}\cos ^{-1}\left(\frac{u_{\rm z}}{\sqrt{u_{\rm x}^2+u_{\rm z}^2}}\right) & \mbox{if} & \sin \lambda <0\end{array} \right., \end{equation} \tag{ 9 }$$

where u_x, u_y, and u_z are the components of the particle 4-velocity vector.

Zoom In Zoom Out Reset image size

We find that the low-energy particles (γ/γ_rad ≲ 0.1) nearly conserve the initially imposed isotropy, because they are still upstream and have not been energized by reconnection. In contrast, the high-energy particles (γ/γ_rad ≳ 0.1) are significantly beamed along the reconnection plane (at X-lines and with flux ropes) within ϕ = ±15° and λ = ±60°. The λ = ±60° angle is of special interest here because it coincides with the direction of the undisturbed magnetic field lines outside the reconnection layers for a α = 0.5 guide field (λ₀ = ±tan ⁻¹(1/α) ≈ ±63°). The particles are accelerated along the z-direction by the reconnection electric field, and move back and forth across the layer mid-plane following relativistic Speiser orbits (Cerutti et al. 2013). At the same time, the particles are deflected away by the reconnected field and the guide field creating a characteristic "S" shape in the angular maps. To a lesser extent, the zero-guide field case also presents some degree of anisotropy, but the deformation and then the disruption of the layer by the kink instability effectively broaden the beams.

The synchrotron angular distribution closely follows the particle one, essentially because relativistic particles radiate along their direction of motion within a cone of semi-aperture angle ∼1/γ ≪ 1. However, there is a noticeable offset between the distribution of the highest-energy particles with γ ≳ γ_rad and the radiation above 100 MeV. This discrepancy is due to the different zones where particles accelerate and where particles radiate. In the accelerating zone, the electric field is intense and leads to linear particle acceleration along the z-axis, whereas the perpendicular magnetic field, B_⊥, is weak deep inside the reconnection layers, yielding little synchrotron radiation. These high-energy particles then radiate ≳ 100 MeV emission abruptly, i.e., within a fraction of a Larmor gyration, only when they are deflected outside the layer where B_⊥ ∼ B₀. The beam dump is well localized at λ = ±60°, i.e., along the upstream magnetic field lines (see hot-spots in Figure 14, bottom-right panel).

As a consequence of this strong anisotropy, the observed spectra of particles and radiation depend sensitively on the viewing angle. Figure 15 compares the isotropic particle (top panel) and photon energy distributions (apparent intrinsic isotropic luminosities, νL_ν, bottom panel) with the distributions along one of the directions dominated by the highest-energy particles (ϕ = −92, λ = 345) and radiation (ϕ = 55, λ = −636) as it appears to a distant observer at tω₀ = 288, for α = 0.5 (see Section 4.4 and Figure 14). The bottom panel in Figure 15 also compares the Fermi-LAT measurements of the gamma-ray spectra of the 2009 February, 2010 September and 2011 April flares and the average quiescent spectrum of the Crab Nebula (Abdo et al. 2011; Buehler et al. 2012) with the simulated gamma-ray spectra. The particle spectrum in the (ϕ = −92, λ = 345) direction is very hard, and peaks at γ/γ_rad ≈ 0.6 with a total apparent isotropic energy ≈7 × 10⁴⁰ erg, which corresponds to about half the total magnetic energy dissipated by the end of the simulation, or about 10 times more energy than in the isotropic distribution. The photon spectrum in the (ϕ = 55, λ = −636) direction peaks at about 100 MeV but extends up to ∼1 GeV. The resulting gamma-ray luminosity above 100 MeV is L_γ ≈ 3 × 10³⁵ erg s⁻¹, which is about 40 times brighter than the isotropically averaged flux from the layer, and about 3 times brighter than the observed average (quiescent) luminosity of the Crab Nebula ($L_{\gamma }^{\rm Crab}\approx 10^{35}$ erg s⁻¹, assuming a distance of 2 kpc between the nebula and the observer). Thanks to beaming, the level of gamma rays expected from the model is consistent with the moderately bright flares, such as the 2009 February or 2010 September ones (see Figure 15, bottom panel), but the simulation cannot reproduce the flux of the brightest flares, such as the 2011 April event. Presumably, increasing the box size would be enough to increase the density of the emitting particles >100 MeV and account for the most intense flares.

Zoom In Zoom Out Reset image size

The beam of high-energy radiation is also time variable, both in direction and in intensity. Figure 16 presents the computed time evolution of the synchrotron photon flux integrated above 100 MeV in the directions defined by ϕ = 55, λ = −636 and ϕ = −55, λ = 491, as well as the computed time evolution of the isotropically averaged flux and the observed average flux in the Crab Nebula measured by the Fermi-LAT (Buehler et al. 2012) for comparison. This calculation assumes that all the photons emitted at a given instant throughout the box reach the observer at the same time, i.e., it ignores the time delay between photons emitted in different regions with respect to the observer. Along the directions probed here, the >100 MeV flux doubling timescale is of order $10\hbox{--}20\omega _0^{-1}$ or 3–6 hr, for both the rising and the decaying time, which is compatible with the observations of the Crab flares (Buehler et al. 2012; Mayer et al. 2013), as well as with our previous 2D simulations in Cerutti et al. (2013). Shorter variability timescales may still exist in the simulation but our measurement is limited by the data dumping period set at $T_{\rm dump}\approx 10\omega _0^{-1}$. These synthetic lightcurves also clearly illustrate the effect of the particle beaming on the observed flux discussed in Section 4.5. Along the direction of the beam, the >100 MeV flux can be ≳ 10 times more intense than the isotropically averaged one, and even exceeds the measured quiescent gamma-ray flux of the Crab Nebula by a factor ≈2–3 at the peak of the lightcurves 3.5–4 days after the beginning of the simulation. Each time the beam crosses our line of sight, we see a rapid bright flare of the most energetic radiation emitted in the simulation.

Zoom In Zoom Out Reset image size

The high-energy particles are strongly bunched within the magnetic flux ropes (within magnetic islands in 2D, see Cerutti et al. 2012b, 2013). As a result, the typical size of the emitting regions is comparable to the dimensions of the flux ropes, i.e., of order L_x/10 ≈ 20ρ₀ along the x- and y-directions (Figure 9), which corresponds to about six light-hours. We conclude that particle bunching is at the origin of the ultra-short time variability (3–6 hr) found in the reconstructed lightcurve (Figure 16). Particle bunching and anisotropy help to alleviate the severe energetic constraints imposed by the Crab flares.