TY - JOUR
AU - Baty, H.
AB - Abstract We study the magnetohydrodynamic tearing instability occurring in a double current sheet configuration when a guide field is present. This is investigated by means of resistive relativistic magnetohydrodynamic simulations. Following the dynamics of the double tearing mode (DTM), we are able to compute synthetic synchrotron spectra in the explosive reconnection phase. The pulsar-striped wind model represents a site where such current sheets are formed, including a guide field. The variability of the Crab nebula/pulsar system, seen as flares, can be therefore naturally explained by the DTM explosive phase in the striped wind. Our results indicate that the Crab GeV flare can be explained by the DTM in the striped wind region if the magnetization parameter σ is around 105. magnetic reconnection, MHD, plasmas, methods: numerical 1 INTRODUCTION In a recent study (Baty, Petri & Zenitani 2013, hereafter Paper I), a model has been proposed to address the origin of the strong flares observed in X-rays and gamma-rays from pulsars and magnetars environments (Palmer et al. 2005; Striani et al. 2011; Takamoto et al. 2014). These flares are short and powerful. Indeed, the observations of Crab pulsar nebula show flare duration between a few hours and up to several days with a rising/falling time-scale of a few hours/days (Buehler et al. 2012; Striani et al. 2013). This duration appears to be too large to be associated with the rotation of the neutron star (period of 33 ms), and on the other hand too short for typical nebula dynamical time-scale (≈1 yr). Relativistic magnetic reconnection is usually believed to be an efficient mechanism to explain these flares (Clausen-Brown & Lyutikov 2012), possibly in a way similar to solar flares where magnetic energy is suddenly released and converted into other forms of energy. Indeed, the highly magnetized wind in the nebula entails a magnetic energy reservoir that is large enough to explain the typical observed flare energy (≈1034 J). However, achieving a fast enough time-scale for releasing the magnetic energy is more tricky. For example, the standard view of tearing instabilities developing in a single current sheet leads to a magnetic reconnection (with growth of magnetic islands) on a relatively slow time-scale that is also very dependent on the (unknown) resistive dissipation (Lyutikov & Uzdensky 2003; Komissarov, Barkov & Lyutikov 2007). Recently, Uzdensky, Cerutti & Begelman (2011) and Cerutti, Uzdensky & Begelman (2012) discussed that photons emitted by an extremely accelerated particle in an X-line may be able to explain the Crab flares, and Cerutti et al. (2014) calculated the evolution of a relativistic current sheet in the pulsar wind nebula using three-dimensional particle-in-cell (PIC) simulations. However, they did not succeed in reproducing the observed synchrotron spectra, in particular the 2011 April flare event. In Paper I, an explosive mechanism is described, that is based on a magnetic reconnection process associated with the double tearing mode (DTM). The favoured site for the emission corresponds to a region situated in the stripped wind close to the light cylinder radius (r ≈ 50rL). Indeed, the DTM is a magnetohydrodynamic (MHD) instability that is well known to develop in multiple current sheets, as expected from the magnetic structure of pulsars magnetosphere with the presence of a current sheet wobbling around the equatorial plane (Coroniti 1990). The explosive character of the mechanism is due to the sudden and fast development of a secondary instability that is structurally driven by the interaction of the initial magnetic islands situated on two successive current layers. It also manifests by a merging of the magnetic islands and is characterized by a bulk plasma flow reaching a magnitude close to the Alfvén speed. More recently, the robustness of the mechanism has been shown by performing new relativistic MHD simulations of the DTM of magnetically dominated plasmas, shifting to more realistic high magnetization parameter (σ) and magnetic Lundquist number (S which is a measure of the dissipation) than previously explored in Paper I (Pétri et al. 2015, hereafter Paper II). The results indicate a weak dependence of the time-scale of the explosive phase on the magnetic Lundquist number and the magnetization parameter, S0.2−0.3 and σ0.3, respectively. In this paper, with the aim to improve our model, we extend the previous results obtained in the two-dimensional (2D) approximation, by adding the effect of a magnetic guide field component perpendicular to the main plane. Such structure is indeed expected from a pulsar current sheet structure when the rotation and magnetic axis do not coincide. Indeed, the magnetic field structure in the striped wind can be described by an analytical expression given by Pétri (2013). The main magnetic component is toroidal, along the azimuthal direction Bφ, whereas the poloidal part is directed along the radius Br. Because of this particular magnetic configuration, we get Bφ ∝ rL/r and Br ∝ (rL/r)2 both having the same magnitude at the light cylinder. Therefore, the effect of a guide field in the pulsar striped wind should become negligible for distances much larger than rL. We re-examine the scaling laws dependence with σ and S, with a special emphasis on the effect of plasmoid chain, that can develop in a transient way in the high S regime. We also investigate the energetics of the reconnection process on the basis of our approximate magnetized fluid model. More precisely, the change of the temperature distribution during the explosive phase is determined in order to deduce an effective thermal particles spectrum, and the ensuing synchrotron emission. The results are then compared to the observations. The paper is organized as follows. The numerical setup is described in Section 2. The results of the DTM simulations are presented in Section 3. Section 4 is devoted to the computation of the synchrotron spectra deduced from the temperature profile during the explosive reconnection event. Before ending with the conclusions in Section 6, we discuss the implications of our work for the relevant parameters in the striped wind (sheet thickness and other parameters) in Section 5. 2 NUMERICAL SETUP The evolution of the 2D double current sheets is modelled by using the relativistic resistive magnetohydrodynamic (RRMHD) approximation. The RRMHD equations are solved by using a numerical method developed by Takamoto & Inoue (2011) which solves the equations in a conservative form using a finite volume method. Magnetic field is updated with the constrained transport algorithm (Evans & Hawley 1988) which allows us to treat the divergence-free magnetic field. In this paper, the numerical setup follows Paper I, with an extra magnetic field component BG added to model the guide field. The initial current sheets are described by assuming the static relativistic Harris current sheets. The initial magnetic field profile is \begin{equation} {\boldsymbol B} = B_0 \left[1 + \tanh [(y - y_0)/ l] - \tanh [(y + y_0)/ l] \right]{\boldsymbol x} + B_{\rm G} {\boldsymbol z} ,\!\! \! \end{equation} (1) and the initial density profile is \begin{equation} \rho = \rho _{\rm b} + \rho _0 \left[1 / \cosh ^2 [(y - y_0) / l] + 1 / \cosh ^2 [(y - y_0) / l] \right] ,\!\! \! \end{equation} (2) where y0 = 3l is the half separation of the current sheets, and l is the half thickness of the sheets. The gas pressure in the sheets is determined by the pressure balance. Note also that, for the sake of simplification, the guide magnetic field component is assumed to be uniform. The initial temperature T is assumed to be uniform and constant, kBT/mc2 = 1, where kB, m, and c are the Boltzmann constant, particle mass, and light velocity, respectively. The corresponding sound velocity is cs ≃ 0.516c. The upstream background density ρb is determined by the magnetization parameter σ ≡ B2/4πρhc2, where h is the specific enthalpy. The relativistic ideal equation of state is assumed, in normalized units h = 1 + Γ/(Γ − 1)p/ρc2, where Γ = 4/3. The resistivity coefficient η is assumed to be uniform and constant, and is related to the Lundquist number S through S = 4πcl/η (4π is necessary due to the Gauss unit). Note that in this paper, we use the light speed as the characteristic velocity to define our Lundquist number. The simulation domain is bounded by a rectangular domain of dimensions [0, Lx] × [−Ly, Ly], with a fixed size value Lx = 20l. As the boundary condition of x-direction is taken to be periodic, the Lx value is precisely chosen in order to get a linearly unstable mode with a wavelength close to the fastest one according to MHD stability analysis (see Papers I and II for more detail). The Ly value is taken to be large enough to prevent the boundary effects on the sheets evolution, as free boundary conditions are used in this lateral y-direction. We divide this simulation domain into homogeneous numerical square meshes with size Ly/Ny. A small divergence-free magnetic field perturbation is also added (at t = 0) as a vector potential δAz at the cell-edges in order to trigger the system evolution \begin{equation} \delta A_z = - B_1 \frac{L_x}{2 \pi } \cos \left[\frac{2 \pi }{L_x} (x - L_x/2) \right] \exp [- (|y| - y_0)^2] , \end{equation} (3) where B1 = 3 × 10−4B0 is the perturbation amplitude. The parameters σ, S, Ly, Ny, BG used in the simulations are given in Table 1.1 The actual numerical value of the parameters, such as ρ0, ρb, B0, follow Paper I. Table 1. List of the parameters. Name . σ . S . Ly . Ny . BG/B0 . cA, x/c . SA1 0.2 200 10l 1024 0 0.408 SA2 0.4 0.535 SA3 0.6 0.612 SA4 1.2 0.739 SA5 12 0.961 SA6 120 0.996 SB1 0.2 3200 10l 2048 0 0.408 SB2 0.4 10l 2048 0.535 SB3 0.6 20l 4096 0.612 SB4 1.2 20l 4096 0.739 SB5 12 20l 4096 0.961 SB6 120 20l 4096 0.996 RA1 0.4 400 10l 1024 0 0.535 RA2 800 1024 RA3 1600 1024 RA4 3200 2048 RC1 12 400 10l 1024 0 0.961 RC2 800 1024 RC3 1600 1024 RC4 3200 2048 GA1 12 1600 20l 2048 0.01 0.961 GA2 0.1 0.956 GA3 0.5 0.859 GA4 1 0.679 GA5 1.5 0.533 GB1 0.2 1600 20l 2048 0.5 0.365 GB2 1.2 0.661 GB3 12 0.859 GB4 120 0.891 GC1 12 200 20l 2048 0.859 GC2 800 Name . σ . S . Ly . Ny . BG/B0 . cA, x/c . SA1 0.2 200 10l 1024 0 0.408 SA2 0.4 0.535 SA3 0.6 0.612 SA4 1.2 0.739 SA5 12 0.961 SA6 120 0.996 SB1 0.2 3200 10l 2048 0 0.408 SB2 0.4 10l 2048 0.535 SB3 0.6 20l 4096 0.612 SB4 1.2 20l 4096 0.739 SB5 12 20l 4096 0.961 SB6 120 20l 4096 0.996 RA1 0.4 400 10l 1024 0 0.535 RA2 800 1024 RA3 1600 1024 RA4 3200 2048 RC1 12 400 10l 1024 0 0.961 RC2 800 1024 RC3 1600 1024 RC4 3200 2048 GA1 12 1600 20l 2048 0.01 0.961 GA2 0.1 0.956 GA3 0.5 0.859 GA4 1 0.679 GA5 1.5 0.533 GB1 0.2 1600 20l 2048 0.5 0.365 GB2 1.2 0.661 GB3 12 0.859 GB4 120 0.891 GC1 12 200 20l 2048 0.859 GC2 800 Open in new tab Table 1. List of the parameters. Name . σ . S . Ly . Ny . BG/B0 . cA, x/c . SA1 0.2 200 10l 1024 0 0.408 SA2 0.4 0.535 SA3 0.6 0.612 SA4 1.2 0.739 SA5 12 0.961 SA6 120 0.996 SB1 0.2 3200 10l 2048 0 0.408 SB2 0.4 10l 2048 0.535 SB3 0.6 20l 4096 0.612 SB4 1.2 20l 4096 0.739 SB5 12 20l 4096 0.961 SB6 120 20l 4096 0.996 RA1 0.4 400 10l 1024 0 0.535 RA2 800 1024 RA3 1600 1024 RA4 3200 2048 RC1 12 400 10l 1024 0 0.961 RC2 800 1024 RC3 1600 1024 RC4 3200 2048 GA1 12 1600 20l 2048 0.01 0.961 GA2 0.1 0.956 GA3 0.5 0.859 GA4 1 0.679 GA5 1.5 0.533 GB1 0.2 1600 20l 2048 0.5 0.365 GB2 1.2 0.661 GB3 12 0.859 GB4 120 0.891 GC1 12 200 20l 2048 0.859 GC2 800 Name . σ . S . Ly . Ny . BG/B0 . cA, x/c . SA1 0.2 200 10l 1024 0 0.408 SA2 0.4 0.535 SA3 0.6 0.612 SA4 1.2 0.739 SA5 12 0.961 SA6 120 0.996 SB1 0.2 3200 10l 2048 0 0.408 SB2 0.4 10l 2048 0.535 SB3 0.6 20l 4096 0.612 SB4 1.2 20l 4096 0.739 SB5 12 20l 4096 0.961 SB6 120 20l 4096 0.996 RA1 0.4 400 10l 1024 0 0.535 RA2 800 1024 RA3 1600 1024 RA4 3200 2048 RC1 12 400 10l 1024 0 0.961 RC2 800 1024 RC3 1600 1024 RC4 3200 2048 GA1 12 1600 20l 2048 0.01 0.961 GA2 0.1 0.956 GA3 0.5 0.859 GA4 1 0.679 GA5 1.5 0.533 GB1 0.2 1600 20l 2048 0.5 0.365 GB2 1.2 0.661 GB3 12 0.859 GB4 120 0.891 GC1 12 200 20l 2048 0.859 GC2 800 Open in new tab 3 RESULTS As previously reported in Paper I, the evolution of the system exhibits four distinct stages: (1) an initial linear growth of two magnetic islands arranged asymmetrically in the two current layers and developing on a time-scale S1/2; (2) an ensuing saturation on the slow diffusion time-scale S1; (3) an explosive growth of the islands interacting with a triangular deformation and driving a fast magnetic reconnection event characterized by a merging between the two islands; (4) a relaxation towards a new stable state. This overall behaviour does not change in highly magnetized plasma regime with large Lundquist numbers, as reported in Paper II. An important result for the robustness of our model is that, a relatively weak dependence of the time-scale of the explosive phase on the magnetic Lundquist number (S0.2−0.3) and on the magnetization parameter (σ0.3) is obtained. Interestingly, the explosive stage is however characterized by a transient growth of smaller magnetic islands identified as plasmoid-chain instability. The plasmoids appear when a critical local Lundquist number is reached, as it is the case for runs using S = 3.2 × 103 (see fig. 1 in Paper II). Thus, in order to quantify the effect of the plasmoid chain, we re-examine in somewhat more details the scaling law dependence of the explosive stage. 3.1 Scaling law dependence and plasmoid-chain effect Following Papers I and II, we use the maximum four-velocity in x-direction, ux = γVx, as a characteristic variable of the DTM evolution. First, we examine the scaling law of ux on the magnetization parameter for different Lundquist number values. The results are shown in panels of Fig. 1. The top and bottom panels are for the low and large Lundquist number cases, respectively. In the low-Lundquist number case, runs SA1-6, we found that the maximum four-velocity scales as ∼σ0.5 and ∼σ0.07 in the small σ and large σ regimes, respectively.2 This indicates that the maximum four-velocity in the low-σ region is proportional to the Alfvén four-velocity: |$\gamma _{\rm A} c_{\rm A} = \sqrt{\sigma }$|⁠. However, we note that the obtained maximum velocity is not the Alfvén velocity but approximately just half of it. This is consistent with the outflow velocity of the single relativistic tearing instability with small Lundquist number, indicating the deceleration by the large enthalpy in the outflow region (Takahashi et al. 2011). The bottom panel of Fig. 1 shows the maximum four-velocity γVx in the case of S = 3200, runs SB1-6. Differently from the low-Lundquist number case, the maximum four-velocity can be well reproduced by the Alfvén four-velocity, |$\sqrt{\sigma }$|⁠. This is because, in the large Lundquist number case, the sheets become plasmoid chain as reported in Paper II, which includes the Petschek-type structure (Baty 2012; Takamoto 2013). Figure 1. Open in new tabDownload slide Maximum four-velocity γVx of DTM with respect to the magnetization parameter σ. Top: Low-Lundquist number cases, S = 200: runs SA1–6. Bottom: Large Lundquist number cases, S = 3200: runs SB1–6. Figure 1. Open in new tabDownload slide Maximum four-velocity γVx of DTM with respect to the magnetization parameter σ. Top: Low-Lundquist number cases, S = 200: runs SA1–6. Bottom: Large Lundquist number cases, S = 3200: runs SB1–6. Next, we examine the scaling law of γVx on the Lundquist number. The numerical results are reported in panels of Fig. 2. The top and bottom panels are in the small and high σ cases, respectively. In the small σ case, the maximum velocity scales approximately as S0.3 as reported in Paper I. In the large σ case, the maximum four-velocity scales as S0.3−0.4 as in the small Lundquist number regime. However, the scaling becomes ∼S0.8 and gradually saturates in the large Lundquist number regime. This is because the maximum velocity becomes very close to the Alfvén four-velocity in the upstream region, |$\sqrt{\sigma } \simeq 3.5$|⁠, when σ = 12. Figure 2. Open in new tabDownload slide Maximum four-velocity γ Vx of DTM with respect to the Lundquist number. Top: small magnetization cases, σ = 0.4: runs RA1–4, SA2. Bottom: high magnetization cases, σ = 12: runs RC1–4, SB5. Figure 2. Open in new tabDownload slide Maximum four-velocity γ Vx of DTM with respect to the Lundquist number. Top: small magnetization cases, σ = 0.4: runs RA1–4, SA2. Bottom: high magnetization cases, σ = 12: runs RC1–4, SB5. 3.2 Guide magnetic field effects In order to investigate the effect of a magnetic guide field component, the direction of the total magnetic field is varied with a fixed magnitude, that is, |${\boldsymbol B} = B \cos \Phi {\boldsymbol x} + B \sin \Phi {\boldsymbol z} \equiv B_0 {\boldsymbol x} + B_{\rm G} {\boldsymbol z}$|⁠, where B is the constant total magnetic field strength, and Φ is the angle between x-axis and the magnetic field direction. The corresponding runs are GA1–GA5 in Table 1. Fig. 3 is a snapshot of temperature profile of run GA4. The basic spatial structure is roughly similar to the case without guide field. However, the resulting temperature from the burst phase is much smaller than the case without guide field (see Fig. 10). This is because the compression of the guide field reduces the increase of gas pressure. Note also that the sheet region does not have plasmoid chain, differently from the case without guide field as reported in Paper II. This is due to the stabilization of tearing instability by guide field (Somov & Verneta 1993), which prohibits the secondary-tearing instability leading to the plasmoid chain. Fig. 4 is the temporal evolution of the maximum four-velocity in x-direction. It shows that the initial single tearing instability is slowed down by guide field as is predicted by Somov & Verneta (1993) in non-relativistic work. It also shows the evolution of secondary phase is slowed down and the maximum velocity becomes slower as increasing the guide field component. Figure 3. Open in new tabDownload slide A snapshot of temperature profile at the burst phase of run GA4. Figure 3. Open in new tabDownload slide A snapshot of temperature profile at the burst phase of run GA4. Figure 4. Open in new tabDownload slide Temporal evolution of the maximum four-velocity with different guide field. Figure 4. Open in new tabDownload slide Temporal evolution of the maximum four-velocity with different guide field. Fig. 5 plots the guide field dependence of the maximum four-velocity. In the figure, the maximum four-velocity exhibits a decrease with increasing guide field as reported by Zenitani, Hesse & Klimas (2009) and Zanotti & Dumbser (2011). This can be interpreted by a decrease of the reconnection magnetic field component B0, which results in the decrease of the Alfvén velocity, |$c_{\rm A} = \sqrt{\sigma _x/(1 + \sigma )}$|⁠, along the sheets.3 This is due to the decrease of the reconnection magnetic field strength as: B0 = B cos Φ. Fig. 6 shows the Rutherford time (see Papers I and II), defined as the time at which the explosive phase is triggered. The results show that the Rutherford time increases linearly. Note that the slope is very steep, and the Rutherford time becomes even twice when BG = 0.5B0. Note that these properties are very similar to the behaviour of the relativistic Petschek slow shocks with guide field (Lyubarsky 2005). Figure 5. Open in new tabDownload slide The maximum four-velocity γ Vx: runs GA1–GA5. Figure 5. Open in new tabDownload slide The maximum four-velocity γ Vx: runs GA1–GA5. Figure 6. Open in new tabDownload slide The Rutherford time (see text for the definition). The horizontal axes are the ratio of the guide field to the reconnecting magnetic field component. Figure 6. Open in new tabDownload slide The Rutherford time (see text for the definition). The horizontal axes are the ratio of the guide field to the reconnecting magnetic field component. Fig. 7 is the maximum four-velocity dependence including guide field which is fixed as BG/B0 = 0.5. The top panel shows the magnetization parameter dependence. It shows the dependence is completely different from the no guide field case in Fig. 1, and the guide field drastically reduces the maximum four-velocity, although the maximum four-velocity slightly increases with σ. This is because the guide field hinders the evolution of plasmoid chain as indicated in Fig. 7, and the sheets become usual relativistic Sweet–Parker sheet whose outflow saturates an upper limit before reaching the Alfvén velocity (Takahashi et al. 2011). The bottom panel shows the Lundquist number dependence. It shows the maximum four-velocity increases with Lundquist number similar to no guide field case shown in Fig. 2. However, the power-law index of the Lundquist number is around 0.6 which is an intermediate value indicated in the bottom panel of Fig. 2, that is, between 0.4 and 0.8. These panels clearly show that the guide field basically weakens the energy conversion resulted from DTM, and the burst phase becomes less explosive even in highly magnetized plasma with high-Lundquist number. Figure 7. Open in new tabDownload slide Top: maximum four-velocity dependence on σ, run GB1–GB4. Bottom: maximum four-velocity dependence on S, GC1, GC2, GB3. In these calculations, the ratio BG/B0 is fixed. Figure 7. Open in new tabDownload slide Top: maximum four-velocity dependence on σ, run GB1–GB4. Bottom: maximum four-velocity dependence on S, GC1, GC2, GB3. In these calculations, the ratio BG/B0 is fixed. Fig. 8 is the energy distribution to each component in the reconnected region of runs GA2 and GA4. It shows most magnetic field energy is converted into the thermal energy in the both cases. However, the guide field energy occupies around 30 per cent of the total energy in GA4. This indicates the guide field magnetic field pressure compensates the compression by DTM, and reduces the conversion into thermal energy. This is basically consistent with the two-fluid result by Zenitani et al. (2009). However, the thermal energy is still dominant in the case of GA4, and this reflects the fact that DTM is not related to slow shocks. Figure 8. Open in new tabDownload slide Energy density in the thermalized region after DTM burst with guide field component, corresponding to runs: GA2 and GA4. Figure 8. Open in new tabDownload slide Energy density in the thermalized region after DTM burst with guide field component, corresponding to runs: GA2 and GA4. 4 DTM THERMAL ENERGY SPECTRUM AND CRAB FLARES In our previous works (Paper I; Paper II), we did not discuss in detail the radiation properties associated with the DTM evolution. In order to obtain more quantitative results about possible radiative signatures of the reconnection phase, in this section, we compute the thermal synchrotron spectrum expected to emanate from the relativistic and highly magnetized plasma during its explosive phase. The procedure is explained in the following lines, our aim being to compare our synthetic spectra with the variability of the Crab gamma-ray flares. 4.1 Temperature profile As already discussed in Papers I and II, the DTM shows an explosive reconnection phase during which temperature drastically increases in the reconnected region even in the case of moderate value of the magnetization parameter σ. In this paper, we explore ultrarelativistic Poynting-dominated flows assuming a high-Lundquist number for the plasmas, and show the resulting explosive reconnection becomes more violent, that is, the flow speed becomes faster and the temperature becomes higher, as the magnetization σ-parameter and the Lundquist number are increased. Since the Crab pulsar wind is considered to be a high-σ plasma with very large Lundquist number, we expect the DTM to be sufficiently energetic to explain the flares in the Crab pulsar/nebula system. Fig. 9 shows the profile of temperature just before and after the explosive phase of DTM in a highly magnetized plasma, σ = 120, with high-Lundquist number, S = 3200. As is indicated in Fig. 9, the high-Lundquist number plasma gives very narrow reconnection sheets in which several plasmoids can be observed (see also our Paper II). The plasmoid region becomes very dense and also very hot with temperature as high as kBT ≳ σmc2 because of the compression by nearly Alfvénic reconnection flows along the sheet (Takamoto 2013). After the explosive phase, the magnetic field between the two original sheets is forced to reconnect, and dissipate its energy into kinetic bulk flow and thermal energy. Since the motion inside of the reconnected region becomes highly stochastic, the resulting kinetic bulk flow energy rapidly dissipated into the thermal energy. In particular, the random motion induces strong compression, and this also increases the temperature in the plasmoids. Figure 9. Open in new tabDownload slide The temperature profiles log10[kBT/mc2] during the double tearing burst phase of the run SB6. The top and bottom panels show just before and after the burst phase, respectively. Figure 9. Open in new tabDownload slide The temperature profiles log10[kBT/mc2] during the double tearing burst phase of the run SB6. The top and bottom panels show just before and after the burst phase, respectively. 4.2 Synchrotron energy spectrum Next, we investigate the synchrotron energy spectrum using our numerical results. The typical photon energy for the synchrotron emission can be written as \begin{equation} \epsilon _{\rm sync} = \frac{3}{2} \gamma ^2 \frac{B}{B_{\rm q}} m c^2 , \end{equation} (4) where B is the magnetic field strength as measured in the frame where ϵsync is detected, Bq = m2c3/eℏ is the critical magnetic field, m is the electron mass, e is it electric charge, γ is the typical Lorentz factor of the electrons, and ℏ is the reduced Planck constant. Assuming the synchrotron radiation is mainly emitted by the pair plasma which is in local thermal equilibrium due to the MHD approximation, the Lorentz factor in equation (4) can be expressed for an ultrarelativistic plasma equation of state as: γ  ∼  3kBT/mc2 ≡ 3Θ. For later convenience, we introduced the normalized temperature Θ. Note that the temperature T depends on the location in the simulation box. In the explosive phase, a strong gradient is formed leading to drastic variation in the temperature profile from point to point as seen in Fig. 9. In the pulsar-striped wind region, the dominant background magnetic field component is toroidal and propagates as an entropy wave thus decreasing with radius r according to, B0 = BL rL/r, where rL  ∼  1.5 × 106 m is the radius of the light cylinder and BL  ∼  100 T is the magnetic field strength at rL for the Crab pulsar. The synchrotron photon energy in the fluid comoving frame becomes \begin{eqnarray} \bar{\epsilon }_{\rm sync} &\sim & 1.73 \times 10^{-4\,} {\rm eV} \left(\frac{3 k_{\rm B} T}{m c^2}\right)^2 \left( \frac{\bar{B}}{1\,{\rm T}} \right) \nonumber \\ &&= 0.156\, {\rm eV}\, \Theta ^2 \Gamma _{\rm W}^{-1} \left(\frac{r}{r_{\rm L}}\right)^{-1} \left( \frac{\bar{B}}{\bar{B}_0} \right) \left( \frac{B_L}{100\,{\rm T}} \right) , \end{eqnarray} (5) where |$\bar{B}_0 = B_0 / \Gamma _{\rm W}$| is the background magnetic field measured in the fluid comoving frame, or the simulation frame and ΓW is the Lorentz factor of the pulsar wind. Finally, the Lorentz transformation into the pulsar rest frame (observer frame) gives us \begin{equation} \epsilon _{\rm sync,L} = \frac{\bar{\epsilon }_{\rm sync}}{\Gamma _{\rm W}\,(1-\beta _{\rm W}\,\cos \vartheta )} \equiv \delta _{\rm W} \bar{\epsilon }_{\rm sync}, \end{equation} (6) where ϑ is the angle between the velocity of the reconnecting blob and the observer line of sight, βW is the velocity of the wind in the unit of light velocity, and δW is the ‘Doppler factor’. If we assume that this blob is pointing towards Earth in the most favourable case, δW ≃ 2ΓW, the peak energy of the photons will be \begin{equation} \epsilon _{\rm sync,L} \sim 0.311\, {\rm eV} \Theta ^2 \left( \frac{r}{r_{\rm L}} \right)^{-1} \left( \frac{\bar{B}}{\bar{B}_0} \right) \left( \frac{B_L}{100\,{\rm T}} \right). \end{equation} (7) The energy spectrum can be obtained by integrating equation (7) over the numerical spatial domain, ∫dVnϵsync, where n is the number density of emitting particles in the laboratory frame. The top panel of Fig. 10 shows the energy spectra with several value for the magnetization parameter σ corresponding to the runs SB1–SB6 at the burst phase when reaching maximum four-velocity. Typical values for the Crab pulsar are: r = 50rL4 and ΓW = 300.5 This panel clearly shows the energy in their body part increases with σ-parameter, approximately 100 times increase as σ-parameter increases by one order of magnitude. This can be understood from the bottom panel of Fig. 10. The panel is the maximum temperature during the explosive phase in terms of the σ-parameter. It shows the maximum temperature increases linearly with the σ-parameter. Since the photon energy depends on the temperature as T2, as indicated in equation (7), the resulting energy spectrum can be estimated as T2 ∝ σ2. The observed Crab flares have their mean energy around 102–103 MeV energy region (Buehler et al. 2012; Bühler & Blandford 2014), so that the DTM can explain the Crab flares if the Crab pulsar wind has σ ≳ 105 which is in agreement with current estimates from theoretical expectations (Kirk, Lyubarsky & Petri 2009) about a trans-Alfvénic flow. We also note that the energy spectrum in Fig. 10 increases with σ-parameter. This is because, in high σ cases, the number of plasmoids along the sheets increases as reported in Takamoto (2013). In plasmoid region, the plasma is characterized by a high temperature and a high density, which increases the high-energy photons and make the energy spectrum harder.6 Note that the plasma is far below the radiation reaction limit because in the wind frame the energy of the synchrotron photons is much less that the limit of 240 MeV. Therefore, the dynamics of the plasma presented in this paper is not significantly perturbed by the radiative losses. Our RRMHD code does not include any such losses so far. We computed the spectra by a post-processing algorithm. Figure 10. Open in new tabDownload slide Top: the energy spectrum with various magnetization parameter σ. Bottom: maximum temperature with respect to background magnetization parameter. In both cases, S = 3200 is assumed, and the pulsar parameters are r = 50rL and ΓW = 300. The runs are corresponded to runs SB1–SB6, and the data at the burst phase reaching maximum four-velocity are used. Figure 10. Open in new tabDownload slide Top: the energy spectrum with various magnetization parameter σ. Bottom: maximum temperature with respect to background magnetization parameter. In both cases, S = 3200 is assumed, and the pulsar parameters are r = 50rL and ΓW = 300. The runs are corresponded to runs SB1–SB6, and the data at the burst phase reaching maximum four-velocity are used. Finally, the synchrotron cooling time in the wind comoving frame can be written as \begin{eqnarray} \bar{\tau }_{\rm sync} &=& \frac{3}{4} \frac{m c}{\gamma \sigma _{\rm T} U_{\rm B}} \nonumber \\ &\simeq & 2.3 \times 10^7 {\rm s} \nonumber \\ &&\times \left(\frac{\Gamma _{\rm W}}{300}\right)^2 \Theta ^{-1} \left( \frac{\bar{B}}{0.05 \bar{B}_0} \right)^{-2} \left( \frac{r}{50 r_{\rm L}} \right)^2 \left( \frac{\bar{B}_{\rm L}}{100\, {\rm T}} \right)^{-2} ,\nonumber\\ \end{eqnarray} (8) where σT is the Thomson scattering cross-section and |$\bar{U}_{\rm B}$| is the magnetic field energy density in the fluid comoving frame. Note that this time-scale is estimated in the merged sheet where the magnetic field release occurs. In our simulations, the dynamical time of the explosive phase in the simulation frame, |$\bar{\tau }_{\rm dyn}$|⁠, is about 50l/c. Using l = 2παΓWrL,7 it reduces to \begin{equation} \bar{\tau }_{\rm dyn} = 25\, {\rm s} \left(\frac{\alpha }{0.05} \right) \left(\frac{\Gamma _{\rm W}}{300}\right) \left( \frac{P}{33\,{\rm ms}} \right) , \end{equation} (9) where P = 2πrL/c is the rotation period of the central neutron star.8 Comparing these two time-scales, we obtain \begin{eqnarray} \frac{\bar{\tau }_{\rm sync}}{\bar{\tau }_{\rm dyn}} &=& 9.4 \times 10^5 \left(\frac{\Gamma _{\rm W}}{300}\right) \Theta ^{-1} \left( \frac{\bar{B}}{0.05 \bar{B}_0} \right)^{-2} \left( \frac{r}{50 r_{\rm L}} \right)^2 \nonumber \\ &&\times \left( \frac{\bar{B}_{\rm L}}{100\, {\rm T}} \right)^{-2} \left(\frac{\alpha }{0.05} \right)^{-1} \left( \frac{P}{33\,{\rm ms}} \right)^{-1} . \end{eqnarray} (10) This indicates the maximum energy by DTM in the Crab pulsar wind would be obtained when Θ  ∼  106. Note that we obtain a cut-off energy scale around 300 MeV by substituting this temperature into equation (7) assuming |$\bar{B}\, {\sim }\, 0.05 \bar{B}_0$| and r = 50rL, which roughly reproduces the observed cut-off energy, 350 MeV. Fig. 10 also indicates this temperature can be obtained when σ  ∼  105. 5 APPLICABILITY TO CRAB PULSAR WIND In this paper, we discuss the DTM to explain the Crab GeV flares. For the simulation setup, we assume a separation of two consecutive current sheets to be around three times of their thickness. Although the sheet separation being fairly confidently constrained to be of the order of the light cylinder radius using the geometry of the striped wind model (Michel 1971), there is still no accepted theoretical estimate giving the thickness of one sheet. However, from an observational point of view, the Fermi Large Area Telescope (Fermi/LAT) detected more than hundred gamma-ray pulsars showing small scale details in their light curves in the GeV band (Abdo et al. 2010). If we assume that the striped wind current sheets are responsible for the gamma-ray pulses, the width of each pulse should reflect the thickness of the sheet itself. Indeed, for sufficiently relativistic outflows, the spiral structure combined to relativistic beaming effects gives rise to pulsed emission. If the current sheet has a half thickness of l, then the half width of the corresponding pulse should be l/π β rL, where β ≈ 1. For the Crab pulsar, the pulse width is about 10 per cent of the pulsar period thus the underlying current sheet thickness should be around 10 per cent of the striped wind wavelength (Pétri & Kirk 2005; Pétri 2012). We consider this supports our MHD treatment of the double current sheets since the indicated sheet width is much larger than the kinetic scales, such as the electron's Larmor radius. Next, we discuss the conditions required to trigger the DTM in the Crab pulsar wind. The first step is the triggering of the simple tearing mode, which needs kl ≲ 0.6 so that l ≲ λ/10, where λ = 2π/k is the wavelength of the perturbation.9 The second step is that the DTM requires ky0 ≲ 1 such that 2πy0 ≲ λ. In our case y0 = 3l, thus 6πl ≲ λ which is approximately comparable to the first condition. According to the pulse profile of the Crab, we have l ≈  2πrL/10; therefore, simple tearing is triggered for wavelengths λ ≳ 2πrL which will be easily satisfied since the sheet length in the striped wind is much longer than the sheet separation 2πrL. The DTM also demands the distance to a boundary from sheets should be larger than the sheet separation. Although the sheet separation in the striped wind is considered to be equal around the equator, in a high-latitude region, it gradually becomes a combination of long separation λl and short separation λs which satisfies λl + λs = 2πrL. This means the DTM responsible for the GeV flare should occur around high-latitude region in the wind as indicated in Fig. 11. Since the inclination angle of Crab pulsar is considered to be larger than 45° (Rankin 1993; Moffett & Hankins 1999), we expect it is always possible to find a proper latitude in order for the growth of the DTM. Figure 11. Open in new tabDownload slide A meridional profile of current sheets in the striped wind. The inclination angle is assumed π/4, and the light cylinder radius rL is used for the unit of each axis. A solid line is z = x tan θ where θ = π/6 along which the distance between the sheets becomes imbalanced. Figure 11. Open in new tabDownload slide A meridional profile of current sheets in the striped wind. The inclination angle is assumed π/4, and the light cylinder radius rL is used for the unit of each axis. A solid line is z = x tan θ where θ = π/6 along which the distance between the sheets becomes imbalanced. Finally, as is indicated in Fig. 5, the DTM becomes less explosive as the guide field increases. We consider the guide field in the current sheets will be small in the striped wind, and our scenario can still be applicable in the actual Crab pulsar wind. This is because, first, the radial magnetic field component decreases with increasing the radial coordinate as: Br ∝ (r/rL)−2, so that it is negligible at r ≃ 50rL. Secondly, the magnetic field in the striped wind is basically generated by the oscillation motion of the pulsar magnetosphere around the light cylinder where the magnetic field is theoretically considered to be antiparallel across the current sheet in the equator. Hence, it is likely that the magnetic field in the wind region is also nearly antiparallel across the sheets of striped wind. One possibility resulting in a strong guide field component in the current sheets is the existence of a monopole-like radial magnetic field in the sheet at the light cylinder. As indicated in Fig. 12, such a magnetic field will immediately change into a poloidal component in the sheet due to the oscillation motion, and survive globally in the striped wind region. However, the theoretical studies of pulsar magnetospheres do not support this monopole-like magnetic field structure, and this will be unlikely to occur. For these reasons, we expect the guide field to be very small, and does not reduce the energy release. Figure 12. Open in new tabDownload slide A schematic picture of generating guide field in the striped wind current sheets via a radial magnetic field component in the pulsar magnetosphere. The green line is the radial magnetic field responsible for the guide field in the wind region. The blue lines are the poloidal magnetic field responsible for the reconnection magnetic field in the wind region. Figure 12. Open in new tabDownload slide A schematic picture of generating guide field in the striped wind current sheets via a radial magnetic field component in the pulsar magnetosphere. The green line is the radial magnetic field responsible for the guide field in the wind region. The blue lines are the poloidal magnetic field responsible for the reconnection magnetic field in the wind region. 6 CONCLUSIONS The striped wind structure represents a natural magnetic field configuration arising from the rotation of a magnetized neutron star surrounded by a relativistic pair plasma. It generates a current sheet wobbling around the equatorial plane. From simple analytical models, it is known that the magnetic field lines become mainly toroidal far from the neutron star with a poloidal component decreasing much faster than the toroidal component. Therefore, locally, the striped wind can be depicted as a current sheet with a guide field for which the strength depends on the distance from the light cylinder. Our new RRMHD simulations of the DTM including guide field effect confirm our previous works. DTM is a good candidate to explain the variability of the gamma-ray emission of the Crab flares. Indeed, the explosive phase exhibits a magnetic reconnection event releasing magnetic energy into heat and therefore also into synchrotron radiation. The comparison between the synchrotron cooling time and the dynamical time-scale shows that the plasma regime is far below the radiation reaction limit. The spectrum peaks at an energy εp simply related to the magnetization parameter σ by εp ∝ σ2. This allows us to constrain the magnetization of the Crab to lie around σ ≈ 105. This is in addition to the previous constraints on the location of the origin of the flares estimated to be around r ≈ 50 rL where rL = c/Ω is the light cylinder radius, c is the speed of light, and Ω is the rotation speed of the pulsar and on the maximum Lorentz factor of the striped wind, Γ ≲ 300 (Paper II). We would like to thank Yasunobu Uchiyama and Dmitry Khangulyan for many fruitful comments and discussions. Numerical computations were carried out on the Cray XC30 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. Calculations were also carried out on SR16000 at YITP in Kyoto University. This work is supported in part by the Postdoctoral Fellowships for Research Abroad programme by the Japan Society for the Promotion of Science No. 20130253 and also by the Research Fellowship for Young Scientists (PD) by the Japan Society for the Promotion of Science No. 20156571 (MT). JP and HB acknowledge financial support from the French National Research Agency (ANR) through the grant no. ANR-13-JS05-0003-01 (project EMPERE). 1 Note that in this paper, we treated high-Lundquist number plasmas. In general, it is very difficult to find the numerical convergence in the presence of plasmoid chain due to its turbulent nature. In this paper, we have checked the convergence of the maximum velocity and temperature at the explosive phase of DTM, which is given in Appendix A. 2 Note that the obtained scaling law on σ-parameter is a little different from our previous work which indicated ∼σ0.3. This is mainly due to including σ = 0.2 which is actually included in a new scaling regime we found for the first time in this paper. 3 Here, |$\sigma _x \equiv B_x^2/4 \pi \rho h c^2 \gamma ^2$| is the magnetization parameter in terms of the reconnection magnetic field component. 4 Here, the emission region is assumed to start approximately at a fixed radius outside the light cylinder following the prescription given by Pétri & Kirk (2005). This is a necessary requirement for both fitting the synchrotron spectrum to the observed value (as the synchrotron emissivity decreases with distance) and adjusting the time-scale of the Crab flares taking into account time dilation. Although this makes it difficult to discuss time-scale of DTM, we use the DTM time-scale as the observed flare time-scale in equation (9) since this is one of the characteristic time-scales of DTM model in the striped wind. The validity of the above assumption will be investigated in a future work. 5 This is twice larger than the value inferred in Paper II, ΓW ≲ 150. This is due to changing of the values of some parameters in the dynamical time of the explosive phase, equation (9), in order to explain the cut-off energy of the observed flare. 6 Our numerical results in the case of σ = 120 gave us an effective index in the high-energy region γF, DT ∼ 1.3, which is not so far away from the observed energy spectral index of the flaring component γF, obs = 1.27 ± 0.12. Comparing with the non-thermal PIC simulation results obtained by Cerutti et al. (2014), the DTM spectrum is harder in high-energy region due to the plasmoids. 7 α is a ratio between the sheet width and the sheet separation (see also Fig. 11). 8 Note that if we consider the Lorentz-time dilation, this gives a flare time-scale: τdyn, lab  ∼  2 h in the pulsar rest frame, which is within the observed flare time-scale, 8 h. 9 Although the initial tearing instability needs somewhat long time to evolve in high-Lundquist number plasma, recent PIC simulations (Cerutti et al. 2015; Philippov, Spitkovsky & Cerutti 2015) show that the current sheets in pulsar winds suffer from strong perturbation around light cylinder radius, and we expect such a strong perturbation helps the initial tearing instability to grow before reaching r = 50rL. REFERENCES Abdo A. A. et al. ApJ 2010 708 1254 Crossref Search ADS Baty H. Phys. Plasmas 2012 19 092110 Crossref Search ADS Baty H. Petri J. Zenitani S. MNRAS 2013 436 L20 (Paper I) Crossref Search ADS Buehler R. et al. ApJ 2012 749 26 Crossref Search ADS Bühler R. Blandford R. Rep. Prog. Phys. 2014 77 066901 Crossref Search ADS PubMed Cerutti B. Uzdensky D. A. Begelman M. C. ApJ 2012 746 148 Crossref Search ADS Cerutti B. Werner G. R. Uzdensky D. A. Begelman M. C. ApJ 2014 782 104 Crossref Search ADS Cerutti B. Philippov A. Parfrey K. Spitkovsky A. MNRAS 2015 448 606 Crossref Search ADS Clausen-Brown E. Lyutikov M. MNRAS 2012 426 1374 Crossref Search ADS Coroniti F. V. ApJ 1990 349 538 Crossref Search ADS Evans C. R. Hawley J. F. ApJ 1988 332 659 Crossref Search ADS Kirk J. G. Lyubarsky Y. Petri J. Becker W. Astrophysics and Space Science Library, Vol. 357, The Theory of Pulsar Winds and Nebulae 2009 Berlin Springer-Verlag 421 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Komissarov S. S. Barkov M. Lyutikov M. MNRAS 2007 374 415 Crossref Search ADS Lyubarsky Y. E. MNRAS 2005 358 113 Crossref Search ADS Lyutikov M. Uzdensky D. ApJ 2003 589 893 Crossref Search ADS Michel F. C. Comm. Astrophys. Space Phys. 1971 3 80 Moffett D. A. Hankins T. H. ApJ 1999 522 1046 Crossref Search ADS Palmer D. M. et al. Nature 2005 434 1107 Crossref Search ADS PubMed Pétri J. MNRAS 2012 424 2023 Crossref Search ADS Pétri J. MNRAS 2013 434 2636 Crossref Search ADS Pétri J. Kirk J. G. ApJ 2005 627 L37 Crossref Search ADS Pétri J. Takamoto M. Baty H. Zenitani S. Plasma Phys. Control. Fusion 2015 57 014034 (Paper II) Crossref Search ADS Philippov A. A. Spitkovsky A. Cerutti B. ApJ 2015 801 L19 Crossref Search ADS Rankin J. M. ApJS 1993 85 145 Crossref Search ADS Somov B. V. Verneta A. I. Space Sci. Rev. 1993 65 253 Crossref Search ADS Striani E. et al. ApJ 2011 741 L5 Crossref Search ADS Striani E. et al. ApJ 2013 765 52 Crossref Search ADS Takahashi H. R. Kudoh T. Masada Y. Matsumoto J. ApJ 2011 739 L53 Crossref Search ADS Takamoto M. ApJ 2013 775 50 Crossref Search ADS Takamoto M. Inoue T. ApJ 2011 735 113 Crossref Search ADS Takamoto M. Kisaka S. Suzuki T. K. Terasawa T. ApJ 2014 787 84 Crossref Search ADS Uzdensky D. A. Cerutti B. Begelman M. C. ApJ 2011 737 L40 Crossref Search ADS Zanotti O. Dumbser M. MNRAS 2011 418 1004 Crossref Search ADS Zenitani S. Hesse M. Klimas A. ApJ 2009 705 907 Crossref Search ADS APPENDIX A: CONVERGENCE OF THE CALCULATIONS In this paper, we dealt with high-Lundquist number plasmas. In general, high resolution is necessary to simulate high-Lundquist number plasma to reduce numerical dissipation. Fig. A1 shows temporal evolutions of the maximum velocity of run RC2 with different resolutions. It shows the resolution approximately Nx ≳ S/2 is sufficient for reproducing the maximum four-velocity. Note that our code does not include viscosity and thermal conductivity but only resistivity, so that it is in general impossible to obtain perfect numerical convergence. Figure A1. Open in new tabDownload slide Temporal evolutions of the maximum four-velocity of run RC2 using different numerical resolution. Figure A1. Open in new tabDownload slide Temporal evolutions of the maximum four-velocity of run RC2 using different numerical resolution. © 2015 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society © 2015 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - Thermal synchrotron radiation from RRMHD simulations of the double tearing mode reconnection – application to the Crab flares
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stv2163
DA - 2015-12-11
UR - https://www.deepdyve.com/lp/oxford-university-press/thermal-synchrotron-radiation-from-rrmhd-simulations-of-the-double-zMuGNRBVtV
SP - 2972
EP - 2980
VL - 454
IS - 3
DP - DeepDyve
ER -