EVOLUTION OF RELATIVISTIC PLASMOID CHAINS IN A POYNTING-DOMINATED PLASMA

Magnetic reconnection is a process that converts magnetic field energy into thermal and kinetic energy very efficiently (Biskamp 2000; Priest & Forbes 2000). Because of this fact, it is believed that magnetic reconnection plays an important role in various phenomena from laboratory plasmas to astrophysical plasmas. Recently, interest in the properties of relativistic magnetic reconnection has been growing, especially in Poynting-dominated plasmas, which are believed to be present in various high-energy astrophysical phenomena such as ultra-relativistic jets (Lovelace & Romanova 2003; Barkov & Baushev 2011), gamma-ray bursts (GRBs; Lyutikov & Blandford 2003; Zhang & Yan 2011), and pulsar winds (Kennel & Coroniti 1984a, 1984b; Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003). In those models, the Poynting energy of the plasma is assumed to be dissipated almost completely into thermal and kinetic energy at some distance from the central object. However, the cause of such an efficient dissipation process is still unknown. In the last decade, several studies have been performed with the goal of finding efficient dissipation processes (Lazarian & Vishniac 1999; Komissarov et al. 2007b; Takamoto et al. 2012; Inoue 2012; Amano & Kirk 2013; Mochol & Kirk 2013a, 2013b). Magnetic reconnection is one of the most promising candidates among these processes, and has been studied actively from analytical (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Lyubarsky 2005) and numerical points of view (Komissarov et al. 2007a; Zenitani et al. 2009b, 2009a; Dumbser & Zanotti 2009; Zenitani et al. 2010; Takahashi et al. 2011; Bessho & Bhattacharjee 2012).

For magnetic reconnection to occur, the plasma should contain current sheets. Such structures evolve into the Sweet–Parker configuration when the Lundquist number S_L ≡ c_AL/η is small (Loureiro et al. 2005), where c_A is the Alfvén velocity, L is the sheet length, and η is the resistivity. It is well known that the reconnection rate of the Sweet–Parker sheet is very slow (Sweet 1958; Parker 1957, 1963), so a considerable number of studies have been devoted to finding a mechanism that enhances the magnetic reconnection. These mechanisms may include the anomalous resistivity in collisionless plasma (Ugai & Zheng 2005; Fujimoto 2011) and the effect of turbulence (Lazarian & Vishniac 1999; Kowal et al. 2009). Recently, it was found that spontaneous current sheet fragmentation in a non-relativistic plasma occurs via secondary tearing instabilities when the Lundquist number exceeds a critical value, leading to the so-called plasmoid chains. The critical value is thought to be about 10⁴ in non-relativistic plasmas (Shibata & Tanuma 2001; Loureiro et al. 2007; Samtaney et al. 2009; Uzdensky et al. 2010; Bárta et al. 2011b; Loureiro et al. 2012; Huang & Bhattacharjee 2013; Loureiro et al. 2013). In those works, it was shown that (1) the reconnection rate is enhanced by plasmoid chains and typically reaches v_R ∼ 10⁻²c_A and (2) the distribution of the plasmoid size is either a power law w^−p or an exponential function exp [ − w/α'] where p is the power law index, w is the plasmoid width, and α' is a constant. Since the plasma temperature in the plasmoid region is higher than that of the background plasmas, the plasmoid chains are expected to generate pulsed emission. Hence, the plasmoid chains are of interest from observational and theoretical points of view, especially in connection with solar flares (Bárta et al. 2011a).

The first study of relativistic plasmoid chains was done by Zanotti & Dumbser (2011). These authors performed two- and three-dimensional numerical simulations of relativistic magnetic reconnection using the relativistic resistive magnetohydrodynamic (RRMHD) approximation (Dumbser & Zanotti 2009). In those calculations, Zanotti & Dumbser (2011) assumed a background plasma with a high-Lundquist number S_L ∼ 10⁵–10⁸, a relativistic temperature k_BT ∼ mc², and a high magnetization parameter with respect to the mass density: $\sigma _m \equiv B_0^2 / 4 \pi \rho _0 \gamma _0^2 \sim 20,$ where B₀, ρ₀, and γ₀ are the background magnetic field in the laboratory frame, the rest mass density and the Lorentz factor, respectively. They also assumed the existence of a local anomalously large resistivity, so that their current sheet became very similar to the Petschek-type one. They found that relativistic magnetic reconnection is similar to Petschek-type reconnection with a critical Lundquist number ∼10⁸, which is much larger than the non-relativistic cases.

In this paper, we investigate the evolution of plasmoid chains in a cold Poynting-dominated background plasma with a large Lundquist number: S_L ∼ 10³–10⁵. In particular, we mainly investigate the statistical properties of plasmoid chains such as the distribution function of the plasmoid width and the dynamics of the X- and O-points along the current sheet. To study the evolution of the secondary tearing instability, we use a uniform, constant resistivity, and initialize the magnetic field with perturbations localized at the origin. This procedure allows us to understand the evolution of current sheets in which a tearing instability is triggered at a point.

In this section, we give a brief review of the non-relativistic plasmoid-chain theory. It is widely known that current sheets are unstable to the tearing instability. The maximum growth rate of this instability can be expressed as $\omega _{{\rm max}} = 1 / \sqrt{\tau _R \tau _A}$, where τ_R ≡ δ²/η is the resistive diffusion timescale and τ_A ≡ δ/c_A is the Alfvén crossing time across a current sheet, δ is the sheet thickness, η is the resistivity, and c_A is the Alfvén velocity in the background plasma (Furth et al. 1963; Low 1973; Komissarov et al. 2007a). This expression can be rewritten as follows:

$$\begin{equation} \omega _{{\rm max}} = 1 / \sqrt{\tau _R \tau _A} = \left(\frac{\delta }{c_A} \frac{\delta ^2}{\eta } \right)^{-1/2} = \frac{\tau _A^{-1}}{\sqrt{S_{\delta }}}, \end{equation} \tag{ 1 }$$

where S_δ = c_Aδ/η is the Lundquist number related to the sheet thickness δ. This equation shows that the tearing instability grows faster as the sheet thickness δ shrinks. The current sheet thickness behind plasmoids shrinks when the plasmoid grows along the current sheet and this process triggers the growth of other small plasmoids, which are called secondary plasmoids. Hence, we can expect that a current sheet will evolve into stochastic plasmoid chains in a few growth times of the largest plasmoid. Uzdensky et al. (2010) assumed the existence of a critical Lundquist number S_c at which current sheets become unstable to the plasmoid instability and discussed the physical nature of plasmoid chains. This critical value introduces the smallest elementary structure in the chain, called the "critical layer"; the related key parameters are the length scale L_c = S_cη/c_A, the thickness $\delta _c = L_c / \sqrt{S_c}$, and the reconnection rate $v_R = c_A / \sqrt{S_c}$. These authors also showed that the global reconnection rate is independent of the Lundquist number S_L when S_L > S_c and the plasmoid chains reach a statistical steady state. They obtained a global reconnection rate value of: v_R ∼ 10⁻²c_A by assuming that the value of the critical Lundquist number is S_c ∼ 10⁴, in accordance with the results of their numerical simulations.

As explained by Loureiro et al. (2007) and Bhattacharjee et al. (2009), the above expression for the growth rate of the tearing instability ω_max can be reinterpreted as follows. When we consider the Sweet–Parker current sheet, we can obtain a relation between sheet thickness δ and sheet length L: $\delta \sim L / \sqrt{S_L}$, where S_L ≡ Lc_A/η. Using this relation, the above equation can be rewritten as follows:

$$\begin{equation} \omega _{{\rm max}} \sim \frac{c_A}{L} \frac{L}{\delta } \sqrt{\frac{L}{\delta S_L}} \sim \frac{S_L^{1/4}}{\tau _{A,L}}, \end{equation} \tag{ 2 }$$

where τ_{A, L} = L/c_A. This equation means that the growth of the tearing instability becomes very fast when S_L reaches about 10⁴, which Loureiro et al. (2007) considered as the critical value of the Lundquist number. S_L depends on the current sheet length L and this fact means we need a very large numerical domain to study the effect of plasmoid chains.

A more complete derivation is presented in Bhattacharjee et al. (2009), Uzdensky et al. (2010), and Loureiro et al. (2013).

We model a very long current sheet using the RRMHD approximation. We use the RRMHD scheme developed by Takamoto & Inoue (2011) extended to the multi-dimensional case using the unsplit method (Gardiner & Stone 2005, 2008). To preserve the divergence-free constraint on the magnetic field, we use the constrained transport algorithm (Evans & Hawley 1988). We calculate the RRMHD equations in a conservative fashion; the mass density, momentum, and energy are also conserved within machine round-off error. For the equation of state, we assume a relativistic ideal gas with h = 1 + (Γ/(Γ − 1))(p/ρ), where Γ = 4/3, h is the specific relativistic enthalpy, ρ is the rest mass density, and p is the gas pressure.

For our numerical calculations, we prepare a square domain, [0, L_x] × [0, L_z] = [0, 20δ] × [0, 320δ], where δ is the current sheet thickness. We divide this domain into homogeneous numerical meshes with size Δ = 5δ/128 ∼ 0.04δ, which is equivalent to the mesh number N_x × N_z = 512 × 8192. Note that to reduce computational costs we solve only a quarter region of the current sheet and impose the point symmetric boundary condition about (x, z) = (0, 0) following Zenitani et al. (2009a). Hence, the above set up is equivalent to a square domain, [ − L_x, L_x] × [ − L_z, L_z] = [ − 20δ, 20δ] × [ − 320δ, 320δ] divided by the mesh number N_x × N_z = 1024 × 16384.¹ Along the boundaries x = L_x and z = L_z, inflows and outflows are allowed. For the initial condition, we consider the static relativistic Harris current sheet (Hoh 1966; Kirk & Skjæraasen 2003):

$$\begin{eqnarray} B_z(x) &=& B_0 \tanh (x / \delta), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} p(x) &=& p_{{\rm in}} + p_s / \cosh ^2 (x / \delta), \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \rho (x) &=& \rho _{{\rm in}} + \rho _s / \cosh ^2 (x / \delta), \end{eqnarray} \tag{ 5 }$$

where p and ρ are the gas pressure and the rest mass density, respectively, and other variables are set to zero except for a small perturbation of the magnetic field described later. For the upstream region of the current sheet, we consider a cold plasma ρ_in = 10p_in; for the inside of the sheet, we consider a relativistically hot plasma ρ_s = p_s where $p_s = B_0^2 / 8 \pi$. Note that the temperature of the sheet decreases with decreasing magnetic field strength. In this calculation, we use a constant resistivity η in order to concentrate on investigating the effects of plasmoid chains on the reconnection rate. Since it is easy to extend the law of resistivity, we will consider various kinds of resistivity in our future work.² To trigger the initial tearing instability at the origin (x, z) = (0, 0), we add the following small perturbation to the magnetic field:

$$\begin{equation} \delta A_y = - 0.03 B_0 \delta \exp [-(x^2 + z^2) / 4 \delta ^2]. \end{equation} \tag{ 6 }$$

The typical parameters used in our calculations are listed in Table 1. To model magnetic reconnection in high-energy astrophysical phenomena such as relativistic jets, the Y-point of a pulsar magnetosphere, and a GRB, we consider a magnetically dominated plasma with a magnetization parameter σ_in > 1. In the following sections, we present numerical results and consider the effects of plasmoid chains.

In this section, we present numerical results of the tearing instability and evolution of the plasmoid-chain.

4.1. Temperature Profile

Figure 1 shows snapshots of temperature profiles of runs B1–B4 at the time when the largest plasmoid resulting from the initial perturbation reaches the edge of numerical domain. Since plasmoids move at approximately the Alfvén speed of the upstream flow unless the plasmoid inertia is comparable to the magnetic field energy, the escape time is of the order of t_A.

Zoom In Zoom Out Reset image size

First, we find that many plasmoids evolve along the current sheet. As we discussed in the previous section, the evolution of a plasmoid induces a thinner current sheet behind it, leading to a secondary tearing instability and the generation of plasmoid chains. We also find that the thickness of the current sheet between plasmoids decreases and the apparent number of plasmoids increases with increasing magnetization parameter σ. We will discuss these effects in Section 4.2. At the origin (x, z) = (0, 0), we note the existence of a large hot region. This region is an artifact of our assumption of point symmetry about the origin, which means that plasmoids entering the region from above have counterparts entering from below with the same magnitude and opposite speed. The merger of these plasmoids results in a plasmoid with zero momentum at the origin; this plasmoid gradually accumulates matter as the simulation proceeds.

4.2. Reconnection Rate

Figure 2 shows the time evolution of the reconnection rate in units of the Alfvén crossing time, L_z/c_A ≡ t_A. The reconnection rate is defined as:

$$\begin{equation} v_R / c_A \equiv - \frac{c}{B_0 c_A L_z} \int ^{L_z}_0 dz E_y(x = 0, z). \end{equation} \tag{ 7 }$$

The top panel is the result of run B3 plotted using a logarithmic scale. Here, we see that the evolution of the reconnection rate can be divided into three phases, separated in the figure by vertical lines at t = t_A and t = 2.2t_A. To the left of the blue line, the reconnection rate shows exponential growth due to the initial tearing instability at the origin. Between the blue and green lines, the reconnection rate oscillates around a power law growth rate with an index of approximately 2. This region is where the plasmoid chains develop: small plasmoids start to appear, changing the growth rate from an exponential to a power law.³ Finally, the reconnection rate saturates to the right of the green line, which marks the time when the largest plasmoid escapes.

Zoom In Zoom Out Reset image size

In the bottom panel, we compare the reconnection rates of runs B1–B4. We find that runs B2–B4 show very similar evolution after the plasmoid instability is triggered. This result indicates that the reconnection rate in units of the Alfvén velocity, v_R/c_A, becomes nearly independent of the magnetization parameter σ in the strongly magnetized plasma, σ_in > 1, once the plasmoid instability starts. The reconnection rate grows until the largest plasmoid, which is initially triggered at the origin, escapes from the numerical domain, at which point the reconnection rate has increased up to ∼0.05c_A. After this point, the plasmoid chains reach a statistical steady state and the averaged reconnection rate is about 0.03c_A, which is approximately twice that of the relativistic tearing instability without plasmoid chains (Takahashi et al. 2011). Note that the reconnection rate of run B1 is lower than that of other runs. This fact is because in this case the plasmoid chains do not grow sufficiently, as can be seen in the left panel of Figure 1 where only one secondary plasmoid is generated. This effect reduces this run's reconnection rate compared with other runs including the fully evolved plasmoid chains. We will discuss later the reason why the plasmoid instability does not grow in this case.

Figure 3 shows the time-averaged reconnection rate 〈v_R/c_A〉 as a function of the Lundquist number S_L.⁴ The top panel is the relativistically strong magnetic field case, σ_in = 14, and the bottom panel is the non-relativistic magnetic field case, σ_in = 0.14. We calculate the time average of the reconnection rate curves over the plateau region where the plasmoid chains reach a statistical equilibrium state. As in the non-relativistic case, we find that the reconnection rate becomes independent of the Lundquist number when it is larger than a critical value S_c. For small Lundquist numbers, we find the Sweet–Parker sheet dependence $S_L^{-1/2}$ of the reconnection rate predicted by Lyutikov & Uzdensky (2003) and Lyubarsky (2005; see also Equation (A11) of the present paper). In our calculations, the critical value of the weakly magnetized case is S_c ∼ 10⁴, which is similar to the value indicated in the non-relativistic case; on the other hand, in the strongly magnetized case, the critical value is S_c ∼ 2–3 × 10³, which is a little less than that of the weak magnetic field case. This result can be explained as follows. After generating plasmoids, the current sheet between the plasmoids will become a Sweet–Parker current sheet. In this case, the sheet thickness can be obtained by Equation (A4). If we assume the reconnection jet velocity is the Alfvén velocity, the sheet thickness can be written as $\delta = L / \sqrt{2 \sigma _{{\rm in}} S_L}$, where we used Equation (A20) to estimate v_in. This result means that the sheet thickness decreases with increasing magnetic field strength. On the other hand, the growth time of the tearing instability is ${\sim} \sqrt{\delta ^3 / \eta c_A}$. Using these two expressions, the growth time of the tearing instability of the secondary current sheet is

$$\begin{equation} \tau _{{\rm tearing,2nd}} \sim \frac{\tau _{A,L}}{(2 \sigma _{{\rm in}})^{3/4} S_L^{1/4}} \propto \sigma _{{\rm in}}^{- 3/4} c_A^{-5/4}. \end{equation} \tag{ 8 }$$

This result means that as the magnetic field strength becomes stronger, the secondary tearing instability grows faster and the plasmoid instability occurs much more readily, especially along the reconnection jet that resulted from the initially triggered plasmoid. Similarly, using the characteristic wavelength of the tearing instability, λ_tearing ∼ δ[δc_A/η]^1/4, the characteristic wavelength of the secondary tearing instability can be obtained as:

$$\begin{equation} \lambda _{{\rm tearing,2nd}} \sim L / \left[(2 \sigma _{{\rm in}})^{5/8} S_L^{3/8} \right] \propto \sigma _{{\rm in}}^{-5/8} c_A^{-3/8}. \end{equation} \tag{ 9 }$$

This result also indicates that the plasmoid instability evolves more readily as the background magnetization parameter becomes larger. Note that Equation (9) means that a background plasma with a larger magnetization parameter necessitates a smaller Lundquist number with respect to the sheet length for the plasmoid instability due to the smaller characteristic wavelength of the secondary tearing instability. This result also supports the findings shown in Figure 3, which indicate that the critical Lundquist number becomes smaller as the magnetization parameter of the background plasma becomes larger.

Zoom In Zoom Out Reset image size

As pointed out by Uzdensky et al. (2010), the reconnection rate of plasmoid chains can be written as $v_R / c_A \sim 1 / \sqrt{S_c}$ using the relation of the Sweet–Parker sheet. If we use the above critical value, S_c = 3 × 10³, in the strongly magnetized case, the reconnection rate is ∼0.02c_A, which agrees with the values indicated in the top panel of Figure 2.

4.3. The Evolution of the Plasmoid Structure

Figure 1 shows that the aspect ratios of plasmoids take different values, depending on the magnetization parameter σ_in; the aspect ratio seems to decrease as the magnetization parameter σ_in increases. This result can be explained as follows. The left panel of Figure 4 is the density profile of a plasmoid at t = t_A. This figure shows that its aspect ratio is about 14:1 and that the inner structure of the plasmoid is very similar to that of the Petschek reconnection case that was investigated by Zenitani & Miyoshi (2011). The right panel of Figure 4 is the density profile of the same plasmoid at t = 1.2t_A. We find that the plasmoid size in the z-direction shrinks because of the appearance of slow shocks. These shocks are generated by the steepening of slow waves that are induced by collisions with other plasmoids. In the example shown in Figure 4, slow waves are generated by collisions with the plasmoid at z ∼ 48δ in the left panel. As these slow shocks propagate across the plasmoid, the upstream plasma in the plasmoid is compressed and the plasmoid size shrinks in the z-direction. Figure 5 shows the density configuration of the plasmoid triggered by the initial perturbation of runs B1 and B2 at a time just before the plasmoid escapes from the numerical domain. In run B1, σ_in = 0.14 and we find that the aspect ratio of the plasmoid maintains its initial value of approximately 14:1. This result occurs because the plasmoid instability in run B1 does not grow sufficiently, as explained in the previous sections; the largest plasmoid does not experience a collision with a smaller plasmoid. On the other hand, in run B2, many collisions with smaller plasmoids reduce the aspect ratio of the plasmoids to about 6:1. In our calculations, the aspect ratio does not show any rapid time evolution after it reaches the above ratio of 6:1. Although we cannot be certain that this ratio is the final state, it seems that the aspect ratio depends very weakly on time. Finally, we cannot find any strong dependence of the above aspect ratio on the magnetization parameter σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To understand the physical nature of plasmoid chains, it is helpful to trace the trajectories of the X- and O-points that are the magnetic null points: B_x = B_z = 0. At the X-points, the magnetic configuration around them is of the X-type and these are points where magnetic reconnection occurs; at O-points, the magnetic configuration around them is of the O-type and these points are usually the sites of plasmoids. Figure 6 shows a plot of the trajectories of the X- and O-points of run B3. This figure shows that in the initial phase there is only one O-point, which is generated by the initial perturbation at the origin. Around t = t_A, small plasmoids start to develop behind the initial O-point and the number of O-points and X-points gradually increases with time. Around t = 2.2t_A, the initial plasmoid reaches the boundary of the numerical domain and escapes from the domain. After that, the current sheet is filled with X- and O-points and the plasmoid chains are fully evolved. This picture is consistent with the temporal evolution of the reconnection rate shown in Figure 2. Figure 6 shows that most points, particularly those close to the initial plasmoid, move steadily toward larger z. Their velocity is approximately 0.8c. Note that the X- and O-points that are not close to the initial plasmoid gradually start to move in both directions along the z-direction. This region is confined around the origin initially and expands with time as the plasmoid chains evolve. Finally, this region covers the entire simulation domain and the plasmoid chains reach a statistical equilibrium state around t = 3t_A.

Zoom In Zoom Out Reset image size

Concerning the X-points, we find that they are located near the midpoint between two O-points, as is expected since they are generated by the tearing instability that ejects two plasmoids away from the X-point. Figure 6 shows that many X-points move along the current sheet and that most of them disappear after a short time due to the merger of two neighboring plasmoids or the collapse of the X-points (Loureiro et al. 2005). In addition, we find that the X-points move in a way similar to that of the nearest plasmoid, as reported by Bárta et al. (2011b). Since the X-points are thought to play an important role in the particle acceleration, their dynamical time along the current sheet will impose an upper limit on the acceleration time. For example, if we consider a current sheet with plasmoid chains in statistical equilibrium, with a critical Lundquist number S_c, the sheet length between them can be estimated as L_c ∼ S_cη/c_A and the dynamical time can be estimated as

$$\begin{equation} t_{{\rm acc}} \sim L_c / c_A \sim S_c \eta / c_A^2. \end{equation} \tag{ 10 }$$

Direct acceleration by the electric field at the X-points will be limited by this time scale. Using our parameters, the value of the acceleration time is t_acc ∼ 1.5 × 10⁻²t_A. Note that Figure 6 includes X-points whose lifetimes are much longer than the above value. Their typical lifetime is about 1.5 × 10⁻¹t_A and some of them survive for a much longer time. Figure 6 indicates that the X-points accompany large plasmoids that have a somewhat large amount of space around them.

Note that sometimes large holes appear in the current sheet in Figure 6, such as at z = 0 or z = 100δ. This effect is due to the "monster plasmoids" that result from the merger of many smaller plasmoids. In particular, the monster plasmoid at z = 100δ around t = 3t_A shows an interesting behavior. In the initial phase, it behaves in nearly the same way as the other plasmoids. In the later phase, its inertia becomes much larger than that of the surrounding plasmoids and its dynamics start to resemble Brownian motion since it moves stochastically around an average trajectory that has a low velocity.

As discussed in the previous sections, the evolution of a plasmoid induces a secondary tearing instability and generates small plasmoids behind it; the small plasmoids in turn induce more tearing instabilities. As a result, the current sheet evolves into a plasmoid chain. Since the distribution of plasmoid sizes is potentially important for high-energy astrophysical phenomena, we investigate this subject using our numerical results.

The statistical behavior of plasmoid chains was investigated in Uzdensky et al. (2010), Fermo et al. (2010, 2011), Loureiro et al. (2012), and Huang & Bhattacharjee (2012, 2013). In those papers, the authors discussed the time evolution of the distribution function of plasmoid chains using the following model kinetic equation:

$$\begin{equation} \frac{\partial f}{\partial t} + \alpha \frac{\partial f}{\partial \Psi } = \zeta \delta (\Psi) - \frac{f N}{\tau _A} - \frac{f}{\tau _A}, \end{equation} \tag{ 11 }$$

where f(Ψ) is the distribution function, Ψ is the magnetic flux of a plasmoid, $N(\Psi) \equiv \int ^{\infty }_{\Psi } f(\Psi ^{\prime }) d \Psi ^{\prime }$ is the cumulative distribution function, $\alpha \sim B_0 c_A / \sqrt{S_c}$ is the plasmoid growing rate of a plasmoid, τ_A ∼ L/c_A is the Alfvén crossing time of a plasmoid across the plasmoid chains with scale L, and ζ is the magnitude of the source of plasmoids. Thus, the second term on the left-hand side describes the growth of plasmoids; the first term on the right-hand side is the source of plasmoids. The second term on the right-hand side is the loss of plasmoids due to mergers with larger plasmoids; the third term is the advection loss. Some analytical steady-state solutions of Equation (11) in large Ψ regions can be obtained as follows. When the loss of plasmoids is mainly by advection, N ≪ 1, we obtain f∝exp [ − Ψ/ατ_A]; when the loss of plasmoids is mainly due to plasmoid mergers, N ≫ 1, we obtain f ∼ 2ατ_AΨ⁻². In this derivation, we assumed that the speed of plasmoids is of the order c_A, corresponding to the assumption that the plasmoid crossing time is τ_A ∼ L/c_A. Recently, Huang & Bhattacharjee (2013) showed that dropping this assumption allows a solution f∝Ψ⁻¹. Since the magnetic flux can be expressed as Ψ ∼ B₀w, where w is the plasmoid size perpendicular to their current sheet (Uzdensky et al. 2010), the above distribution function of the magnetic flux can be used to find the plasmoid size distribution.

The top panel of Figure 7 shows the time-averaged distribution of the plasmoid size of run B3 with error bars. The time average is taken between t = t_A and t = 2.2t_A, which are equivalent to the starting time of the plasmoid instability and the escaping time of the initially triggered plasmoid, respectively, as indicated in Figure 2. Figure 7 shows that the plasmoid size distribution is consistent with a size distribution with a power law index of −2 for small plasmoids in the range [0.1δ, δ] of the plasmoid width, as predicted by previous works for the non-relativistic case. From the above discussion, this result means that the plasmoid loss is mainly due to plasmoid mergers. This finding is a natural consequence because we consider the distribution at the escape time of the initially triggered plasmoid; no plasmoids can escape from the plasmoid chains at that time due to the presence of the initially triggered plasmoid. In addition, larger plasmoids, around δ < w < 5δ, deviate from the power law index of −2 and tend rather to an index of −1. This result indicates that the velocity of large plasmoids deviates from the Alfvén velocity. This fact is because the large plasmoids have large inertias, which reduce their velocity, as in the case of monster plasmoids. Note that the distribution function of the largest plasmoid size region, around w > 5δ, drops rapidly and clearly deviates from a power law. This fact is because the number of plasmoids is too small to show statistically significant results, as can also be seen from the large error bars of this region.

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 7, we plot the plasmoid size distribution of runs B1–B4. We find that the distribution of the strong magnetic field case, run B4, shows a very similar behavior to run B3; it shows a power law with an index of −2 in the range [0.1δ, δ] of the plasmoid width and changes to an index of −1 for larger plasmoids. In the weak magnetic field cases, runs B1 and B2, the distributions also have an index of −1 in the range of larger plasmoids. However, the distribution of the smaller plasmoid size region seems to have an index of −1, too. We propose that this result occurs because small plasmoids are not sufficiently evolved in these runs to show a clear size dependence.

In Figure 8, we plot the time-averaged distribution functions after the initially triggered plasmoids escaped: t > 2.2t_A. The top panel of Figure 8 shows the time-averaged distribution of the plasmoid size of run B3 with error bars. In the small plasmoid region, the distribution function has an index of −2, similar to the previous case. However, the distribution function of the larger plasmoid region, w > δ, drops rapidly and clearly cannot be approximated by a power law. We consider this fact to be due to the effects of plasmoid loss by advection. Since the initially triggered plasmoid already escaped from the simulation domain in this case, the plasmoids can freely escape from the domain and this fact results in the exponential decay of the distribution function, as indicated by the above discussion using the kinetic equation.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 8 shows a plot of the plasmoid size distribution of runs B1–B4. The behavior of the distribution functions in the small plasmoid region, w < δ, is very similar to that shown in Figure 7. However, in the large plasmoid region, rapid decay is also observed, similar to the strongly magnetized case, σ = 14. Note that the distribution function of the weakly magnetized case, σ = 0.14, seems to be a power law in the large plasmoid region, w > δ. Unfortunately, our data do not encompass a sufficiently large number of plasmoids in this region, so we cannot conclude that this result is statistically correct.

In this paper, we investigated the evolution of plasmoid chains in a high-σ plasma. We modeled the relativistic current sheet with a cold background plasma using the relativistic resistive magnetohydrodynamic approximation and solved for its temporal evolution numerically. We performed various calculations using different magnetization parameters of the background plasma from σ_in = 0.14 to σ_in = 29 and different Lundquist numbers with respect to the sheet length from S_L ∼ 10³ to S_L ∼ 10⁵. The numerical results show that the initially induced plasmoid triggers a secondary tearing instability and that the current sheet is gradually filled with many plasmoids. That is, the current sheet evolves into a plasmoid chain, as predicted by non-relativistic work. As expected, this plasmoid instability enhances the reconnection rate, which grows up to ∼0.05c_A until the initially triggered plasmoid escapes from the simulation domain. Subsequently, the plasmoid chain reaches a statistical equilibrium state and the temporally averaged reconnection rate in a steady state becomes ∼0.03c_A. Since the maximum value of the Alfvén velocity is the speed of light c, our numerical calculation indicates that the maximum reconnection rate of the plasmoid chain is 0.03c. In our calculations, the evolution of the reconnection rate shows a similar behavior in strongly magnetized cases (σ_in > 1). Although the weakly magnetized case, σ_in = 0.14, shows a different behavior, we consider this result to be due to the larger wavelength of the secondary plasmoid instability indicated by Equation (9). Note that the above critical value is much smaller than that obtained recently by Zanotti & Dumbser (2011), who found S_c ∼ 10⁸. We believe this difference comes from their assumption of a relativistically hot background plasma. A high temperature reduces the magnetization parameter $\sigma _{{\rm in}} = B_0^2 / 4 \pi \rho _0 h_0 \gamma _0^2$ and the critical Lundquist number becomes large when the σ parameter is small, as shown in Section 4.2.

We also investigated the behavior of the O-points and the X-points. In our simulations, the initial perturbation is confined to the origin. The triggered plasmoid shrinks the current sheet behind it, inducing secondary tearing instabilities. Those O- and X-points that are close to the triggered plasmoid move in the same direction, but the other points start to move in both directions along the sheet, reflecting the final state of statistical equilibrium of the plasmoid chain. Most of the X- and O-points disappear by merging, which limits their lifetime and therefore limits the time over which particles can be accelerated by the electric field at such points. We estimate this lifetime using the parameters of the plasmoid chain. As predicted for the non-relativistic case, we noted the appearance of "monster plasmoids." Interestingly, our calculations show that monster plasmoids slow down as they evolve because of their increasing inertia. Ultimately, they display Brownian-like motion around fixed points.

Finally, we investigated the plasmoid size distribution. Our numerical results show that in strongly magnetized cases the distribution becomes a power law with an index of −2 in the small plasmoid region and a power law with an index of−1 in the large plasmoid region before the initially triggered plasmoid escapes. These results indicate that the plasmoid loss is mainly due to mergers; the plasmoid velocity is of order the Alfvén velocity in the small plasmoid region, but is lower in the large plasmoid region. This fact is because the plasmoid inertia increases with increasing size, preventing large plasmoids from moving at the Alfvén speed. After the escape of the initial plasmoid, the distribution function in the large plasmoid region shows an exponential decay because of the free advective escape of plasmoids from the domain.

Magnetic reconnection is one of the most efficient mechanisms of magnetic field dissipation and is expected to play an important role in many astrophysical phenomena. As shown in this paper, once the tearing instability evolves and generates plasmoids, the plasmoid chain always evolves in the current sheet between the initially generated plasmoids if the Lundquist number of the current sheets is beyond the critical value, especially in a Poynting-dominated plasma. Since plasmoids are associated with high temperature plasma and accelerated particles, they can be used to explain intermittent observational signals from high-energy astrophysical objects such as pulsed emission from the Y-point of the Crab pulsar magnetosphere (Uzdensky & Spitkovsky 2012) and the multi-timescale TeV flares in blazars (Giannios 2013). In this paper, we assumed a constant resistivity and used an approximate equation of state corresponding to a relativistic, adiabatic gas. Nevertheless, we believe that our results have revealed the general properties of plasmoid chains in a Poynting-dominated background plasma.

We thank John Kirk, Iwona Mochol, Simone Giacche, Seiji Zenitani, Keizo Fujimoto, Takaaki Yokoyama, and Tsuyoshi Inoue for many fruitful comments and discussions. We also thank our referees for many useful comments on our paper. Numerical computations were carried out on SR16000 at YITP at Kyoto University. Calculations were also carried out on the Cray XT4 at the Center for Computational Astrophysics, CfCA, of the National Astronomical Observatory of Japan. This work is supported by Max-Planck-Institut für Kernphysik and the Postdoctoral Fellowships for Research Abroad program by the Japan Society for the Promotion of Science (JSPS) No. 20130253 (M.T.).

In this Appendix, we derive the reconnection rate of the relativistic Sweet–Parker current sheet. The basic relations have already been presented by several authors (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Lyubarsky 2005). Here, we clarify the dependence on the external pressure, following the non-relativistic approach of Priest & Forbes (2000). A schematic picture of the Sweet–Parker current sheet is shown in Figure 9.

Zoom In Zoom Out Reset image size

We assume a steady state plasma that can be well described by a relativistic magnetohydrodynamic approximation at all points other than the X-point. We also assume that the plasma is homogeneous in the y-direction. The background magnetic field of the inflow region is B_in = B_ine_x if z > δ/2 and B_in = −B_ine_x if z < −δ/2; the field in the sheet region is B = B_in ± B_se_z, where |B_in| ≫ |B_s| and is a very small constant. We assume that B = 0 at the X-point. In this case, the electric field E_y is constant and we can obtain the following relation:

$$\begin{equation} B_{{\rm in}} v_{{\rm in}} = B_s v_s, \end{equation} \tag{ A1 }$$

where v is the fluid 3-velocity. In addition, we can obtain the following relation at the X-point where the magnetic field and the flow velocity are zero:

$$\begin{equation} B_{{\rm in}} v_{{\rm in}} = \eta j, \end{equation} \tag{ A2 }$$

where η is the resistivity and j is the current vector described by Ampére's law:

$$\begin{equation} j \sim B_{{\rm in}} / \delta. \end{equation} \tag{ A3 }$$

From Equations (A2) and (A3), the sheet thickness δ can be expressed as:

$$\begin{equation} \delta \sim \eta / v_{{\rm in}}. \end{equation} \tag{ A4 }$$

The mass and the energy conservation equation are given by

$$\begin{eqnarray} \rho _{{\rm in}} \gamma _{{\rm in}} v_{{\rm in}} L &=& \rho _s \gamma _s v_s \delta, \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} (\rho _{{\rm in}} h_{{\rm in}} \gamma _{{\rm in}}^2 + B_{{\rm in}}^2) v_{{\rm in}} L &=& (\rho _s h_s \gamma _s^2 + B_s^2) v_s \delta, \end{eqnarray} \tag{ A6 }$$

where ρ is the rest mass density, γ is the Lorentz factor, L is the curvature scale of the background magnetic field, and δ is the current sheet thickness. h = 1 + Γ/(Γ − 1)p/ρ is the specific entharpy of the ideal gas where Γ = 4/3 is the relativistic heat ratio and p is the gas pressure. Here, we also assume the cold upstream plasma, p_in = 0; in the sheet region we assume a hot plasma ρ_s ≪ p_s whose pressure can be determined through pressure equilibrium, $p_{s} = B_{{\rm in}}^2 / 2 \gamma _{{\rm in}}^2$. Then, the energy equation can be rewritten as

$$\begin{equation} \rho _{{\rm in}} \gamma _{{\rm in}}^2 (1 + \sigma _{{\rm in}}) v_{{\rm in}} L = \left[ \frac{2 B_{{\rm in}}^2 \gamma _s^2}{\gamma _{{\rm in}}^2} + B_s^2 \right] v_s \delta, \end{equation} \tag{ A7 }$$

where σ ≡ B²/ρhγ² is the magnetization parameter. Using Equation (A1), the above equation reduces to

$$\begin{equation} (1 + \sigma _{{\rm in}}) \gamma _{{\rm in}}^2 v_{{\rm in}} = \left[ 2 \sigma _{{\rm in}} \gamma _s^2 + \frac{\sigma _{{\rm in}}}{v_s^2} \gamma _{{\rm in}}^2 v_{{\rm in}}^2 \right] v_s \frac{\delta }{L}. \end{equation} \tag{ A8 }$$

Using Equation (A4),

$$\begin{equation} (1 + \sigma _{{\rm in}}) \gamma _{{\rm in}}^2 v_{{\rm in}}^2 \sim \left[ 2 \gamma _s^2 v_s^2 + \gamma _{{\rm in}}^2 v_{{\rm in}}^2 \right] \frac{\sigma _{{\rm in}}}{S_l v_s}, \end{equation} \tag{ A9 }$$

where S_l ≡ Lc/η is the Lundquist number using the speed of light as the characteristic velocity. From this equation, we can obtain the following relation between v_in and v_s:

$$\begin{equation} \gamma _{{\rm in}} v_{{\rm in}} \sim \sqrt{\frac{2}{S_l v_s - c_A^2}} \gamma _s v_s c_A, \end{equation} \tag{ A10 }$$

where $c_A \equiv \sqrt{\sigma / (1 + \sigma)}$ is the Alfvén velocity. If we consider a plasma with a high Lundquist number S_l ≫ 1, the above equation reduces to

$$\begin{equation} \gamma _{{\rm in}} v_{{\rm in}} \sim \sqrt{\frac{2}{S_l}} \gamma _s \sqrt{v_s} c_A. \end{equation} \tag{ A11 }$$

This equation shows that the inflow velocity, which is equivalent to the reconnection rate, is inversely proportional to $\sqrt{S_l}$, which is the same conclusion drawn from the non-relativistic Sweet–Parker current sheet model. To obtain an explicit solution of the upstream velocity, we have to add another equation to the above equations. Here, we consider the equation of motion along the x-direction. The relativistic hydrodynamical equation of motion in the current sheet is given by

$$\begin{equation} \frac{\rho _s h_s \gamma _s^2 v_s^2}{L} \sim j B_s - \frac{p_o - p_N}{L} \sim \frac{B_{{\rm in}}}{\delta } B_s - \frac{p_o - p_N}{L}, \end{equation} \tag{ A12 }$$

where we used Equation (A3), p_o is pressure at the edge of the current sheet, and p_N ∼ p_s is the pressure at the X-point. Note that the scale L is a characteristic scale length in the above equation and it should be the curvature scale of the background magnetic field. Using the pressure equilibrium, the above equation reduces to

$$\begin{equation} 2 \frac{B_{{\rm in}}^2}{\gamma _{{\rm in}}^2} \frac{\gamma _s^2 v_s^2}{L} \sim \frac{B_{{\rm in}}}{\delta } B_s - \frac{p_o - p_N}{L}. \end{equation} \tag{ A13 }$$

From the mass conservation equation (Equations (A5) and (A1)), we can obtain the following relation

$$\begin{equation} \frac{B_s}{\delta } = \frac{B_{{\rm in}}}{L} \frac{\rho _s \gamma _s}{\rho _{{\rm in}} \gamma _{{\rm in}}}. \end{equation} \tag{ A14 }$$

Substituting this relation into Equation (A13), we obtain

$$\begin{equation} 2 \frac{B_{{\rm in}}^2}{\gamma _{{\rm in}}^2} \frac{\gamma _s^2 v_s^2}{L} \sim \frac{B_{{\rm in}}^2}{L} \frac{\rho _s \gamma _s}{\rho _{{\rm in}} \gamma _{{\rm in}}} - \frac{p_o - p_N}{L}. \end{equation} \tag{ A15 }$$

This equation reduces to

$$\begin{equation} \gamma _s^2 v_s^2 \sim \frac{\rho _s \gamma _s \gamma _{{\rm in}}}{2 \rho _{{\rm in}}} - \frac{\gamma _{{\rm in}}^2}{B_{{\rm in}}^2} (p_o - p_N). \end{equation} \tag{ A16 }$$

From the mass conservation equation (Equations (A5) and (A4)), we can obtain the following relation:

$$\begin{equation} \frac{\rho _s \gamma _s}{\rho _{{\rm in}}} \sim \frac{\gamma _{{\rm in}} v_{{\rm in}}}{v_s} \frac{L}{\delta } \sim \frac{\gamma _{{\rm in}} v_{{\rm in}}^2}{v_s} S_l. \end{equation} \tag{ A17 }$$

Using this equation, Equation (A16) can be rewritten as follows:

$$\begin{equation} \gamma _s \sim \gamma _A \sqrt{1 + \frac{1}{2} \left[ 1 - \frac{p_o}{p_N} \right]} \equiv \gamma _A \alpha, \end{equation} \tag{ A18 }$$

where γ_A is the Lorentz factor of the Alfvén velocity in the upstream region and α is the effect of the pressure gradient. Note that when p₀ > 3p_N, α becomes an imaginary number. This result is because, in this case, the reconnection outflow is prevented by the pressure gradient force and this fact results in the break down of the steady state assumption. Using this equation and Equation (A11), we obtain the following form of the reconnection rate:

$$\begin{equation} \gamma _{{\rm in}} v_{{\rm in}} / c_A \sim \sqrt{\frac{2}{S_L}} \sqrt{\gamma _A c_A \alpha \sqrt{\gamma _A^2 \alpha ^2 - 1}}, \end{equation} \tag{ A19 }$$

where S_L ≡ Lc_A/η is the Lundquist number relating to the Alfvén velocity. When p_o = p_N, the above equation reduces to

$$\begin{equation} \gamma _{{\rm in}} v_{{\rm in}} / c_A \sim \sqrt{\frac{2}{S_L}} \sqrt{\sigma _{{\rm in}}}, \end{equation} \tag{ A20 }$$

which is equivalent to the relation obtained in Lyutikov & Uzdensky (2003). When p_o = 0 or the reconnection outflow is ejected into a very cold region, Equation (A19) means that the reconnection rate is enhanced due to the pressure gradient force. Finally, when p_o > p_N, which is, for example, the case where plasmoids exist at the edge of the current sheet, Equation (A19) means that the reconnection rate is reduced by the pressure gradient force.