NUMERICAL EXPERIMENTS ON FINE STRUCTURE WITHIN RECONNECTING CURRENT SHEETS IN SOLAR FLARES

The system commences to evolve due to the attraction of the magnetic fields of opposite polarity at two sides of the current sheet, leading them to approach one another and to squeeze the sheet. The process is governed by the resistive MHD equations as follows:

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} & + & \nabla \cdot (\rho \mathbf {v})=0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\rho \left [\frac{\partial }{\partial t} + (\mathbf {v} \cdot \nabla) \right ] \mathbf {v}=-\nabla p + \mathbf {J} \times \mathbf {B}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\frac{\partial \mathbf {B}}{\partial t} = \nabla \times (\mathbf {v}\times \mathbf {B})+\eta _{m}\nabla ^{2} \mathbf {B}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\frac{\rho ^{\gamma }}{\gamma -1}\left (\frac{\partial }{\partial t} + \mathbf {v}\cdot \nabla\right)\frac{p}{\rho ^{\gamma }}= \mu _{0} \eta _{m}J^{2}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&p = \frac{1}{\tilde{\mu }}\tilde{R}\rho T, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\mathbf {J} = \frac{1}{\mu _{0}}\nabla \times \mathbf {B}, \end{eqnarray} \tag{ 6 }$$

which include the mass continuity equation, the momentum equation, the magnetic induction equation, the energy equation, the universal gas law, and Ampere's law. The quantities ρ, $\mathbf {v}$, $\mathbf {B}$, p, $\mathbf {J}$, and T are mass density, flow velocity, magnetic field, gas pressure, electric current density, and temperature, respectively; and η_m, γ, μ₀, and $\tilde{R}$ are the magnetic diffusivity, the ratio of specific heats, the magnetic permeability of free space, and the universal gas constant, respectively. In addition to the above equations, magnetic field $\mathbf {B}$ is subject to the divergence-free condition as well

$$\begin{equation} \nabla \cdot \mathbf {B} =0. \end{equation} \tag{ 7 }$$

Normalization of the above important quantities brings $\mathbf {B}$, ρ, p, and $\mathbf {v}$ to the following dimensionless form:

$$\begin{eqnarray} &&\mathbf {B}^{\prime }=\frac{\mathbf {B}}{B_{0}}, \quad \rho ^{\prime }=\frac{\rho }{\rho _{0}}, \quad p^{\prime }=\frac{p}{B_0^2/\mu _0}, \nonumber\\ && \mathbf {v}^{\prime }=\frac{\mathbf {v}}{v_{0}}, \quad T^{\prime }=\frac{\beta _{0} T}{2T_{0}}, \quad J^{\prime }=\frac{J}{J_{0}} \end{eqnarray} \tag{ 8 }$$

with

$$\begin{equation} v_{0}=B_{0}/\sqrt{\mu _{0} \rho _{0}}, \; J_{0}=\frac{B_{0}}{\mu _{0}L}, \; \mbox{and} \; \beta = 2\mu _{0} p_{0}/B_{0}^{2}, \end{equation} \tag{ 9 }$$

where ρ, $\mathbf {v}$, $\mathbf {B}$, and p are the variables with dimensions, those with prime (') are the corresponding dimensionless variables, and ρ₀, $\mathbf {B}_{0}$, T₀, L, and v₀ are the initial mass density, magnetic field, temperature, width of the simulation domain, and Alfvén speed of the ambient background plasma, respectively. The dimensionless time t' and spatial coordinates x' and z' are related to their dimensional counterparts by

$$\begin{equation} x^{\prime }=\frac{x}{L}, \quad z^{\prime }=\frac{z}{L}, \quad t^{\prime }=\frac{t}{t_{0}}, \end{equation} \tag{ 10 }$$

where t₀ = L/v₀ is the Alfvén timescale with respect to L. In our calculations below, the above characteristic values are taken as follows:

$$\begin{eqnarray} && B_{0}=10 \mbox{ G}, \quad T_{0}=10^{6} \mbox{ K}, \quad \rho _{0}=2.9\times 10^{-15} \mbox{ g cm}^{-3}, \nonumber\\ && \mbox{and} \quad L=10^{5} \mbox{ km}, \end{eqnarray} \tag{ 11 }$$

which lead to β = 0.1, v₀ ≈ 525 km s⁻¹, and t₀ ≈ 190 s.

Rewriting Equations (1) through (6) in dimensionless form, we have

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} &+& \nabla \cdot \rho \mathbf {v} =0, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&\rho \left [\frac{\partial \mathbf {v}}{\partial t} + (\mathbf {v} \cdot \nabla) \mathbf {v}\right ]=-\nabla p + \mathbf {J} \times \mathbf {B}, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \frac{\partial \mathbf {B}}{\partial t} &=& \nabla \times (\mathbf {v} \times \mathbf {B}) + {\frac{1}{R_{m}}}{\nabla }^{2} \mathbf {B}, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\frac{{\rho }^{\gamma }}{\gamma -1}\left (\frac{\partial }{\partial t} + \mathbf {v}\cdot \nabla\right) \frac{p}{{\rho }^{\gamma }}=\mu _{0} \frac{1}{R_{m}} {J}^{2}, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} p &=& \rho T, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \nabla \cdot \mathbf {B} &=& 0. \end{eqnarray} \tag{ 17 }$$

Here, R_m = v₀L/η_m is the magnetic Reynolds number based on the characteristic length L and the Alfvén speed v_A. We note here that the prime (') of each dimensionless variable has been ignored in the above equations for simplicity, and any variable appearing hereafter is dimensionless. Equations (12) through (15) will be solved by the SHASTA code in which magnetic vector potential A is involved such that

$$\begin{equation} \frac{\partial {\bf A}}{\partial t}+\nabla \cdot ({\bf A v})={\bf A} \nabla \cdot {\bf v} + R_{m}^{-1}\nabla ^{2}{\bf A} \end{equation} \tag{ 18 }$$

is used to replace Equation (14) to ensure that condition (17) always holds in our calculations. Here A is related to B in the way of B = ∇ × A. In the two-dimensional case, this becomes ${\bf B} = \nabla \times A \hat{y}$ with A being the y-component of A and $\hat{y}$ being the unit vector in the y-direction. Specifically, we have replaced the two separate diffusive equations for B_x and B_z in Weber's SHASTA code by a single diffusive equation for A, and there are two subroutines which compute A from $\mathbf {B}$ and $\mathbf {B}$ from A (see also the discussions of Forbes & Malherbe 1991).

Under such an arrangement for the magnetic field and vector potential, we are able to write the initial and boundary conditions for solving Equations (12) through (17). The initial configuration shown in Figure 1 includes a vertical current sheet that separates two magnetic fields of opposite polarity, and is line-tied to the bottom boundary surface. It is described mathematically by

$$\begin{eqnarray} B_z(x, z, t=0) &=& x/|x|, \quad \mbox{for} \quad |x|> w, \quad \mbox{and} \nonumber \\ B_z(x, z, t=0) &=& \sin(\pi x/2w), \quad \mbox{for} \quad |x|\le w. \nonumber \\ B_x(x, z, t=0) &=& 0, \nonumber \\ \mathbf {v}(x, z) &=& 0,\nonumber \\ \rho (x, z, t=0) &=& (\beta + 1- B^2)/\beta \nonumber \\ p(x, z, t=0) &=& T \rho (x, z, t=0)=\beta \rho (x, z, t=0)/ 2,\qquad \end{eqnarray} \tag{ 19 }$$

where w is the half-width of the current sheet, and all parameters in Equation (19) are dimensionless.

The boundary conditions are arranged in this way: the right, left, and top sides of the simulation domain are the open or the free boundary on which the plasma and the magnetic flux are allowed to enter or exit freely and the boundary at the bottom ensures that the magnetic field is line-tied,

$$\begin{equation} \frac{\partial A(x, 0, t)}{\partial t} =0, \end{equation} \tag{ 20 }$$

and the plasma does not slip and is fixed to the wall,

$$\begin{equation} {\bf v}(x, 0, t) = {\bf 0}, \end{equation} \tag{ 21 }$$

respectively. Combining boundary conditions (20) and (21) with the initial conditions in Equation (19), we found

$$\begin{eqnarray*} \frac{\partial p(x, z, t)}{\partial z} \bigg | _{z=0} &=& 0 \quad \mbox{ and} \\ \frac{\partial \rho (x, z, t)}{\partial z} \bigg | _{z=0} &=& 0. \end{eqnarray*}$$

Since ∂A/∂t = −v_xB_z − η_mJ, Equations (20) and (21) require that J vanish on the bottom boundary as well, namely,

$$\begin{eqnarray*} &&J(x, 0, t) = 0 \quad \mbox{for} \quad t > 0. \end{eqnarray*}$$

The evolution in the system of interest starting with the configuration shown in Figure 1 can be studied via solving Equations (12) through (17) with conditions (19) through (21). The modified SHASTA code is used to solve these equations after several necessary tests. We describe the code and the test in the next section.

The SHASTA code is well suited for studying shock waves since it can sustain a sharp shock transition over only two or three mesh points, and shock transitions in the code are free from spurious oscillations caused by Gibb's phenomenon. The original SHASTA code solves the resistive MHD equations in a conservative form (see below) so that quantities such as mass, momentum, magnetic flux, and total energy are normally conserved to high accuracy (<1% error). However, large errors are sometimes generated in the pressure and temperature if the plasma β is greatly different from unity. These errors can be reduced if the conservative form of the energy equation is replaced by an energy equation for the pressure. Therefore, the conservation of the total energy was not maintained as well as other parameters.

A revised version made by Weber et al. (1979) computes diffusion implicitly using a standard tridiagonal matrix procedure for parabolic equations. When the diffusion is computed the advective terms are excluded, and similarly when the advection is computed the diffusive terms are excluded. As long as the time steps are small compared to the physical timescales, the results are the same as if both sets of terms had been solved simultaneously.

Using the SHASTA code revised by Weber et al. (1979), Forbes & Malherbe (1991) performed numerical experiments of magnetic reconnection and radiative cooling that may take place in the two-ribbon flare process. They noticed that the tearing mode did occur in the sheet and played an apparent role in diffusing the magnetic field. Many plasma blobs resulting from the tearing mode were recognized, and the development of the instability in the linear stage duplicated the results of Furth et al. (1963) in theory, but when the evolution began to develop non-linearly, deviations were seen. However, the low resolution of the grid prevented them from further investigating the case of high R_m and low β.

With the improvement in computational hardware and the available resources, Shen & Lin (2009, 2010) further modified and improved SHASTA by implementing AMR so that a case of a large magnetic Reynolds number (up to 5 × 10⁴) and small plasma β (∼0.1) could be studied smoothly. In the present work, we perform more tests and choose an appropriate approach for investigating the reconnection process in the two-ribbon flare.

Several cases are considered with different combinations of some parameters in addition to those arranged above. These parameters include the initial grid number of the experiment domain in each case, N_g, magnetic Reynolds number, R_m, and the number of refinements of the initial grids via AMR, N_AMR. They are listed in Table 1 for six cases, and the initial values for the plasma β and the half-width of the current sheet, w, are taken as 0.1 for all cases. In these cases, refinement by AMR was conducted for cases 1, 3, 4, and 6, and the refining times N_AMR = 2 for all these four cases (see values of N_AMR listed in Column 4 of Table 1). For cases 2 and 5, on the other hand, no AMR refinement was conducted. In calculations below R_m ranges from 5 × 10² to 5 × 10⁵, and calculations are performed only in the region of x > 0 because of the symmetry in the configuration. Figure 2 displays some results of our tests for different cases.

Zoom In Zoom Out Reset image size

For the low R_m case, for example, case 1 (R_m = 500), reconnection quickly sets in through almost the whole current sheet even when the sheet is still thick, and formation of an X-point in the sheet becomes apparent after a group of magnetic arches (loops) has been created (see also Shen & Lin 2009). Figure 2(a) shows the current density distribution and the magnetic configuration at time t = 24.0 for this case. The AMR technique (N_AMR = 2) is used and calculations in the experiment run smoothly. As a result of low R_m, the magnetic field can be dissipated at a fast rate in the whole current sheet, so no magnetic island develops in the sheet. As the AMR technique was applied, some fine features, including a termination shock (see also discussions of Forbes & Malherbe 1991 for more details), on the top of the loop system can be well recognized.

As R_m increases, reconnection apparently becomes difficult in a thick current sheet. In the case of higher R_m and lower grid resolution without refinement, for example, case 2 (N_g = 151 × 301 and N_AMR = 0 as shown in Table 1), the calculation quickly breaks down as the sheet gets thinner, and the experiment stops. With application of the AMR technique (N_AMR = 2), on the other hand, the calculation was immediately stabilized (case 3), and fast reconnection commences when the sheet gets thin enough and magnetic islands form (see Figure 2(b)). This suggests that the AMR technique helps stabilize the calculation to a certain extent, and formation of the magnetic islands inside the sheet indicates the beginning of the tearing mode instability (e.g., see Furth et al. 1963). The tearing mode produces many magnetic islands, or plasmoids or plasma blobs, inside the current sheet, and the reconnection outflow moves these small structures either upward or downward depending on when and where individual blobs form. We shall see later that magnetic islands are more likely to appear in the reconnection process occurring in high R_m environments.

With much bigger R_m, for instance, R_m = 5 × 10⁴ in cases 4 and 5, a much smaller grid size and shorter time step are required to maintain the stability of calculations. We perform our tests in different approaches for these two cases. In case 4, we start calculations with a lower grid resolution (N_g = 151 × 301 and N_AMR = 2) and apply the AMR technique to calculations and in case 5, we start with a higher but fixed grid resolution (N_g = 301 × 601 and N_AMR = 0). As expected, the experiments in both cases run smoothly and give similar evolutionary pictures of the magnetic reconnection process: the current sheet continues to become thinner as evolution in the system commences. Not until the thickness w decreases to about 1/20 of its initial value do the characteristics of the tearing mode instability begin to appear, namely, the appearance of multiple magnetic islands in the sheet. Figure 2(c) displays the magnetic configuration and the current contours at time t = 34.6 for case 4, which clearly manifests two magnetic islands in the sheet.

Furthermore, we increase R_m to 5 × 10⁵ with higher grid resolution, N_g = 301 × 601, and AMR, N_AMR = 2, in case 6. Fast reconnection becomes more difficult to start, and the sheet gets much thinner than previous cases before magnetic islands form. Figure 2(d) shows the magnetic configuration and distribution of J for case 6 at t = 39.6. A lot of magnetic islands are seen to form along the current sheet, but numerical errors also create obvious disturbances to both the field and J. This suggests that the power of AMR may have reached its maximum at R_m = 5 × 10⁵ for the SHASTA code. However, we also tried increasing the grid resolution to N_g = 601 × 901 and found that the experiment breaks down midway through calculations. We thus realize that the SHASTA code may not be able to deal well with cases of R_m > 10⁵, regardless of whether AMR or very high grid resolution is used. Hence, we stay with the parameter combination for case 5 in our calculations below in order to avoid extra errors and to be sure that the experiment results are reliable.

Before ending this section, we should note two things. First of all, numerical diffusion is inevitable in any numerical experiment, whether or not AMR is used. Even different levels of AMR may effectively cause different levels of numerical diffusion, yielding numerical noise in our results. To identify the numerical noise, we perform an estimate as below.

The magnetic induction equation (Equation (14)) can be written as

$$\begin{equation} \frac{\partial A}{\partial t} = ({\bf v}\times {\bf B})_{y} - \frac{J}{R_{m}} \end{equation} \tag{ 22 }$$

in two dimensions. In the absence of numerical diffusion, both sides of Equation (22) should ideally balance with each other. In the realistic numerical experiment, on the other hand, an error exists such that the two sides do not balance each other as a result of the numerical diffusion. The behavior of this error can be studied by comparing the difference of the values of two sides near the current sheet for specific cases. We conduct this comparison for the parameters adjacent to the current sheet because the numerical diffusion in this region reaches maximum.

Suppose we define a ≡ ∂A/∂t and b ≡ (v × B)_y − J/R_m, and we calculate a and b near the principal X-point (PX-point) for cases 4 and 5 over the time interval between t = 10 and t = 25, then estimate three kinds of relative difference between a and b: |a − b|/|a|, |a − b|/|b|, and 2|a − b|/|a + b|, respectively. The results are plotted in Figure 3, in which dotted curves are for case 4, solid curves for case 5, and various sizes of each symbol are for various relative differences between a and b. As expected, none of these results is zero, which apparently suggests existence of the error due to numerical diffusion. However, the behavior of these errors varies from case to case. In case 4, the grid resolution is low initially, so the errors are fairly large at the beginning. As AMR functions properly, we see that all three errors decrease to about 20% continuously until the pre-set upper limit of the AMR level is reached. In case 5, the AMR function is turned off, but the initial resolution is high. Therefore, the error due to numerical diffusion is controlled at a low level of around 20% and is roughly kept at this level over the whole period of the test.

Zoom In Zoom Out Reset image size

Overall, Figure 3 suggests that the error caused by numerical diffusion is not large when the calculation is in steady progress in both cases 4 and 5. For the fixed grid size of high resolution, this error is controlled at a low level from the very beginning. If the initial resolution of the grid is low but AMR functions properly, however, the error is rather large at the beginning but is eventually suppressed to a low level. In our cases, this error is about 20% as the experiment progresses steadily, namely, the numerical noise included in our results is at the level of about 20% of the physical signal.

Another aspect we need to note is that the information revealed by the above tests indicates that the performance of a code for the given computational resources, with or without AMR, is sensitive to the magnetic Reynolds number R_m. In the case of small R_m in an apparent way, for example, case 1, the code always runs smoothly with low grid resolution (small N_g) and a low level of the AMR technique (small N_AMR or even N_AMR = 0), and no significant features of turbulence occur inside the current sheet (see also Figure 4(a) as well as the results of Forbes & Malherbe 1991). As R_m increases, however, running the code properly becomes tricky and difficult, and high initial grid resolution or a high AMR level is required. However, comparing the results of cases 2 through 5 with one another suggests the trade-off between N_g and N_AMR such that large N_g is required if the level of AMR is low (namely, small N_AMR) and selecting N_g and N_AMR or their combination is determined by R_m for the specific case.

Zoom In Zoom Out Reset image size

The current sheet shown in Figure 1 is about 10 times denser than its surroundings. At an early stage of evolution, the equilibrium in the system breaks down as a result of the initial perturbation imposed on the system. Plasma and magnetic fields on the two sides of the sheet commence to move to one another gradually, causing the sheet to become thinner, but the current sheet does not fragment until evolution reaches the stage at t = 26.8 when the sheet is already very thin because of the high value of R_m (5 × 10⁴). As discussed in Section 2, the thinning process is sensitive to the value of R_m. In the low R_m condition, the current sheet diffuses fast enough and no plasmoid forms inside the sheet. Previous numerical experiments indicated that the threshold ratio of the length-to-thickness of the sheet, 2π, for the tearing mode given by Furth et al. (1963) should be the lower limit. In these experiments, the value of this ratio may exceed 50 before the plasmoids begin to form (e.g., see Ni et al. 2010 and references therein). On the other hand, our experiments show that the evolution of the current sheet represented by the decrease in the sheet thickness 2w is not uniform along the current sheet. Instead squeezing of the sheet due to the plasma inflow mainly occurs at a specific location at around z = 0.48 from which the first X-point forms at t = 23.0. We shall notice later that this X-point eventually becomes the principal one, namely, the PX-point.

We focus on the evolution of the current sheet at the PX-point. We plot values of w near the PX-point versus t within the time interval between t = 15.0 and t = 40.0 in Figure 5 (see the solid curve). In the early phase, the half-width of the current sheet w, which has the initial value 0.1, continues to decrease slowly. At t = 26.8 when the first magnetic island starts to form, w has decreased to 0.0075. After time t = 26.8, w gradually decreases to the minimum value 0.0025 and it subsequently fluctuates around this minimum value. We shall see later that this oscillating behavior is related to the formation of multiple X- and O-points in the sheet.

Zoom In Zoom Out Reset image size

At the same time, the electric current density J inside the sheet evolves as well. Initially, its maximum value at the center of the sheet is around 26. It increases gradually as w decreases. As we did for w, we plot J near the PX-point versus time in Figure 5 (see the dotted curve). The curve shows that J starts increasing significantly after t = 26.8 when the first O-point commences to form. Therefore, time t = 26.8 separates slow and fast evolution stages. By time t = 29.4, J at the PX-point increases rapidly to around 265, which is an order of magnitude higher than its initial value. After that, J decreases slowly as a result of the diffusion of the magnetic field through reconnection.

Following the formation of the first X-point (or the PX-point), magnetic reconnection successively creates closed field lines (flare loops) below the PX-point. At t = 26.8, the first O-type neutral point (or magnetic island) starts to form above the PX-point at z = 0.61, representing the ignition of the tearing mode instability in the current sheet. Figure 4 shows the distribution of several variables in the x–z plane at time t = 26.8 (first row) and t = 34.0 (second row), respectively. These variables include the magnetic field together with ρ (left panel), J (left middle panel), v_x (right middle panel), and v_z (right panel). Carefully tracking the evolution of each of these variables, we notice that J continues to increase in the progress of the PX-point formation, and the maximum value of J occurs at the PX-point. This indicates that magnetic reconnection has begun and the fastest reconnection takes place at the PX-point.

However, what happens to the mass is different. The mass density inside the sheet remains almost unchanged in this stage. With the decrease in the sheet thickness, the mass was continuously expelled out of the sheet in two directions (see the contours of v_z in the right panel of Figure 4). The process continues until the first O-point forms above the PX-point at z = 0.61 after t = 26.8. We identify time t = 26.8 with the moment when the tearing mode begins. Following its formation, the first O-point moves upward and leaves the PX-point quickly due to the reconnection outflow. The mass density inside the sheet reaches a dynamical balance with the outside, and the conjugate reconnected magnetic field line at the opposite side of the X-point forms a new flare loop and leaves the X-point quickly as well because of the loop shrinkage that has long been reported in both observations (e.g., Švestka et al. 1987; Forbes & Acton 1996; Vršnak et al. 2004, and references therein) and theories (Lin et al. 1995, 2004, and references therein).

In the same stage, as shown by velocity contours in the right middle panel of the first row in Figure 4 as well as in the associated movies, plasma inflows outside the sheet are faster near the PX-point than at the other locations, and they are always symmetric with respect to the PX-point. The region where the inflow speed value (|v_x|) reaches its maximum rapidly expands in both x- and z-directions as the first O-point forms, and it bifurcates into two as the O-point leaves the PX-point. One of the two regions moves with the O-point, and another one stays on the top of flare loops. In the process, the maximum value of the inflow velocity continues to increase from 0.01 to about 0.1. Simultaneously, the outflow speed v_z increases faster than the inflow speed. Unlike the inflow speed, however, the speed of the outflow along the sheet is not symmetric with respect to the PX-point such that the upward outflow is roughly two times faster than the downward outflow until the first O-point forms. For example, the speeds of both outflows increase from zero as the evolution begins. That of the upward flow reaches 0.24 and that of the downward flow reaches 0.12 as the O-point appears. Then the speed of the downward flow quickly increases, and it eventually catches up with that of the upward flow as the first O-point fully forms.

The second O-point starts to form below the PX-point at around t = 29.8, and the third one appears above the PX-point after t = 30.2. By the time t = 34.0, a series of O-points (plasmoids, plasma blobs, or magnetic islands) has formed in the current sheet. With these new plasmoids forming, multiple X-points are created between every pair of plasmoids. The lower row of Figure 4 shows that some plasma blobs with denser material are located at heights z = 0.57, 0.70, 1.03, 1.11, and 1.53. Meanwhile, the PX-point is also moving upward at a speed apparently lower than that of the other X- and O-points.

As we can see from velocity panels in the lower row of Figure 4, motions of plasma blobs and the reconnection outflow are completely mixed with one another, and they manifest rich kinematic features. For the convenience of showing their motions, we identify several blobs in the sheet and plot their locations along the sheet (or on the z-axis due to the symmetry) versus time in Figure 6. The background of this plot is black and each bright dot represents a blob. Blobs A and B are two specific examples of these blobs such that Blob A is moving upward and Blob B is moving downward.

Zoom In Zoom Out Reset image size

In addition, we also plot locations of the PX- and S-points at various times in Figure 6, in which the solid curve is for the PX-point and the dashed one is for the S-point. We noticed that these two special points are not co-located in space although they are very close to one another, as suggested by Seaton (2008), Oka et al. (2008), T. G. Forbes (2009, private communications), and Murphy (2010). Figure 6 shows clearly that motions of blobs bifurcate at the region where these two points are located. Those moving upward eventually leave the domain of the experiment through the open top boundary. For example, Blob A leaves the domain at around t = 35. Because the bottom boundary condition is line-tied and the end of each field line there is fixed, reconnected magnetic flux is accumulated above the boundary creating flare loops, and plasma blobs moving downward will thus eventually stop their motions at the top of flare loops. For example, Blob B collides with the closed flare loop beneath and stops its motion at t = 33.

Motions of blobs also manifest accelerations, and that of the upward motion is higher than that of the downward motion. For example, Blob B forms near z = 0.28, starts moving at a small velocity at around t = 30.5, and continues to speed up until it hits the closed field lines at t = 32.8 when its velocity reaches v = 0.4. At the same time, the height of the loop top is about z = 0.18. From its formation at t = 30.5 to its collision with the loop top at t = 32.8, the lifetime of Blob B is short and is only around Δt = 2.3.

Blobs moving upward are accelerated more rapidly. More information is revealed by them since they have significantly more room to move. Both Figure 6 and the associated movie clearly indicate that merging occurs frequently among these blobs. For example, in the motion of Blob A, at least three blobs created later successively catch up and merge with it at t = 31.8, t = 33.5, and t = 35.4, respectively. The velocity of the resultant blob following the merging remains almost unchanged although the blobs catching up with the initial Blob A moved significantly faster than Blob A before merging. From the plot in Figure 6, we can deduce the average speed of Blob A before and after merging, which is around 0.41. This time interval in which merging occurs roughly covers the process of formation and disappearance of Blob B as shown by Figure 6.

Besides the motion of blobs in the z-direction, we also notice the motion of the reconnected plasma outside blobs. It is very interesting such that, along the motion direction, the plasma behind a blob may flow faster than this blob, and it quickly slows down when it gets close to the blob. In Figure 7, we plot the plasma velocity distribution along the z-axis at time t = 26.8 (solid curve) and t = 34.0 (dashed curve), respectively, corresponding to the plots in Figure 4. The two curves in Figure 7 show that the plasma speeds vary smoothly along the z-axis before the formation of blobs. Following the formation of blobs, however, flow patterns fluctuate and change, especially around individual blobs, and plasma outside blobs may flow faster than blobs at some specific locations.

Zoom In Zoom Out Reset image size

To describe this phenomenon more clearly, we also specify the locations of some blobs at t = 34.0 (see the arrows) on the dashed curve and their velocities (see the positive and negative numbers). For the upward flow, for example, we identified three blobs that are located near z = 1.0 and z = 1.5. The two blobs near z = 1.0 have velocities v = 0.46 and v = 0.28, but the plasma behind (or below) them flows significantly fast. For example, the blob near z = 1.5 moves at speed v = 0.30, and the reconnected plasma flow behind could move at speeds up to v = 0.8. The same situation occurs in the downward flow. Two blobs moving downward have been identified near z = 0.5 and z = 0.7, which move at speeds of v = −0.17 and v = −0.30, respectively, and the plasma flow behind (or above) them moves at speeds up to v = −0.65 and v = −0.50, respectively.

According to the characteristic values given in Equation (11), the above values of the speed of the plasma blobs inside the reconnection outflows range from 147 to 242 km s⁻¹ for those in the upward flow, and from 89 to 159 km s⁻¹ for those in the downward flow. Obviously, the upward motion is faster than the downward motion, which is consistent with that deduced from observations. Table 2 lists all the results about the speed of the plasma blob observed in various situations that we are able to collect. These results indicate that, in general, the anti-sunward blobs move faster than the sunward ones, and apparently the former is much faster than the latter. This is probably due to the fact that the closed flare loop below the current sheet constitutes an obstacle to the motion of the plasma blob. Comparing our results with observational ones indicates that the present experiment displays an intermediate case as a result of a relatively weak background magnetic field (10 G) and high plasma β (∼0.1).

As we mentioned earlier, Figure 6 also plots positions of the S- and PX-points on the z-axis as functions of time. The dashed curve is for the S-point, near which the reconnection outflow is roughly static, and the solid curve is for the PX-point, near which the magnetic field almost vanishes and the reconnection rate reaches maximum. These two curves display the overall motion pattern of the S- and PX-points, which are upward, and a couple of smooth fluctuations are also seen. The two points move apparently slower than blobs, but the formation of blobs causes their motions to speed up. This scenario is consistent with the fact that the formation of blobs inside the current sheet accelerates the reconnection process and yields quicker production, accumulation, and growth of flare loops below the two points. It is the growth of the flare loop that pushes the two points upward.

As seen in Figure 6, the S- and PX-points are in general not co-located, and their relative positions switch back and forth as reconnection progresses. The initial motion of the two points is smooth, with the S-point above the PX-point. Following the appearance of blobs at t = 26.8, the relative positions of the two points begin to fluctuate. This pattern of motion is seen more clearly in Figure 8, which shows the distance between the PX-point and the S-point versus time. From Figures 6 and 8, it is apparent that the direction in which a plasma blob moves is related to the relative positions of the S- and PX-points. When the S-point is above the PX-point, newly appearing blobs move upward. When the S-point is below the PX-point, newly appearing blobs move downward. Because new blobs are not observed to appear in the short distance between the S-point and PX-point, this implies that plasma blobs never go through the PX-point.

Zoom In Zoom Out Reset image size

Another important parameter for the reconnection process is the rate of magnetic reconnection, M_A, which is the Alfvén Mach number of the reconnection inflow and is defined as the inflow velocity, v_i, near the current sheet compared to the local Alfvén speed, v_A:

$$\begin{equation} M_{A}=\frac{v_{i}}{v_{A}}. \end{equation} \tag{ 23 }$$

After multiple X-points appear in the nonlinear phase of the evolution, the magnetic reconnection process displays violent and nonlinear characteristics, and M_A is not necessarily a constant along the current sheet. Figure 4 suggests that M_A reaches its maximum near the PX-point. This is consistent with observations of the reconnection processes taking place in the magnetosphere and other numerical experiments in various environments (e.g., see also discussions of Lin et al. 2008 and references therein for more details). Therefore we focus on the behavior of M_A near the PX-point. To obtain M_A, we locate a small square with a side length of 0.1 right to the PX-point and evaluate v_i and v_A on average over the square. As suggested in Figure 4, v_i is actually identified with v_x near the PX-point.

Figure 9 plots M_A at the PX-point as a function of time. An overall feature of the curve is that M_A increases slowly in the early phases and manifests significant increases after plasma blobs appear in the current sheet. As Figure 9 shows, M_A remains less than 0.01 until t = 28.0 when the first blob forms. Then it jumps to 0.05 within a very short period from t = 28.0 to t = 29.6. Consequently, it evolves in a violent fluctuating fashion such that its maximum value reaches 0.17 and its minimum value drops to 0.023. This is probably due to the impact of the development of the blob. No observation or in situ measurement of the reconnection process in space manifested this behavior of M_A. That may be due to poor spatial resolution in the relevant observations and measurements, but recent numerical experiments display such a feature quite often (e.g., see Skender & Lapenta 2010; Shepherd & Cassak 2010, and references therein). In this period, its average value gradually increases from 0.05 at t = 29.6 to 0.085 at t = 37.0 (see the dotted curve in the plot). After t = 37.0, reconnection becomes slow again because the magnetic fields weaken after most of the magnetic energy has been released.

Zoom In Zoom Out Reset image size

To further improve our understanding and knowledge about the impact of the tearing mode and the associated turbulence on various properties of the reconnecting current sheet, we check the scaling law for the current sheet thickness deduced by Lin et al. (2007, 2009) on the basis of the results obtained in the present work. According to the classical theory of the tearing mode instability (Furth et al. 1963), the tearing mode develops within the timescale τ which is shorter than the classical diffusion timescale τ_d and is longer than the Alfvén timescale τ_A (see also Priest & Forbes 2000, pp. 177–189). This directly yields the following constraints on the wave number k and the sheet half-thickness w (see also Lin et al. 2007):

$$\begin{equation} S^{-1/4} < \bar{k} < 1, \end{equation} \tag{ 24 }$$

where $\bar{k} = kw$, and S = τ_d/τ_A is the Lundquist number (or the magnetic Reynolds number; in fact R_m = S) near the current sheet. Lin et al. (2007) further showed that M_A = S⁻¹. To help determine the lower limit to the observed w in order to remove projection effects to the greatest extent from the observational results for w, Lin et al. (2007, 2009) related w_min to M_A and the wavelength of the tearing mode λ, which is identified with the distance between two adjacent plasma blobs flowing along the current sheet in the same direction. The relationship can be summarized in the form below

$$\begin{equation} w_{{\rm min}} = M_{A}^{\alpha }\frac{\lambda }{2\pi } \end{equation} \tag{ 25 }$$

for different types of tearing modes and different boundary conditions. Here α is a dimensionless parameter between 0 and 1.

Loureiro et al. (2007) and Ni et al. (2010) found that a Sweet–Parker current sheet with a large aspect ratio is unstable to the long wave perturbation in a high Lundquist number environment, and the maximum growth rate of the instability in the linear stage is related to the Lundquist number in terms of S^1/4v_A/L_cs, where L_cs is the length of the current sheet used in these works. The related numerical simulations suggest that the critical value of the aspect ratio for the tearing mode is about 10² in the case of S = 10⁴ (see Loureiro et al. 2005), and the number of plasmoids appearing inside the current sheet in a given time interval at the linear stage is found to scale as S^3/8 (Loureiro et al. 2007; Ni et al. 2010). In the nonlinear regime, the reconnection rate is nearly independent of S, and the number of plasmoids is roughly proportional to S (Bhattacharjee et al. 2009; Huang & Bhattacharjee 2010).

Huang & Bhattacharjee (2010) also found that the thickness of the sheet is related to the scale of the system, L, as LR^1/2_mc/R_m, where R_mc is the critical value of the magnetic Reynolds number at which the plasmoids start to appear in the reconnecting current sheet. In their experiments, R_m is between 5 × 10⁵ and 2 × 10⁶, and R_mc = 4 × 10⁴. This brings the scale of the sheet thickness into the range 10⁻⁴L to 4 × 10⁻⁴L. Comparing our results shown in Figure 5 and those deduced from Equation (25) indicates that the sheet thickness of Huang & Bhattacharjee (2010) is at least an order of magnitude smaller than what we have obtained in the present work. Therefore, the difference between the results of the two works is apparent.

However, we need to note that the analyses and results of other works mentioned above were conducted and obtained in symmetric configurations. Obviously, the reconnection occurring in either the solar eruption or in the geomagnetic tail is not in a symmetric configuration. Instead the evolution in the magnetic field is subject to the condition of line-tying to the bottom surface and no restriction is imposed on the field in other directions. Therefore, both the tearing mode and the plasmoid instabilities studied here should behave differently from those taking place in the cases investigated by Loureiro et al. (2007), Bhattacharjee et al. (2009), Huang & Bhattacharjee (2010), and Ni et al. (2010), and we do not expect the same scaling law found in this work as those for the cases of the symmetric configuration to exist.

Following the practice of Lin et al. (2007, 2009), in this work we take λ as the distance between two adjacent plasma blobs after they are fully formed and recognizable, M_A is measured near the PX-point so that it can be estimated from Figure 9, and w_min is measured at the middle point between the two corresponding blobs. From Figures 4 and 6, we track the blobs moving upward because Lin et al. (2007, 2009) tracked those leaving the Sun. We identified three pairs of such blobs, and plot their locations, z_b, on the z-axis versus time in Figure 10 (see symbol "+" in the figure).

Zoom In Zoom Out Reset image size

Figure 10(a) displays the locations of a pair of blobs in the time interval between t = 31.2 and t = 33.3. In this period, their distance λ gradually increases from 0.29 to 0.61, the corresponding values of M_A vary from 0.03 to 0.088, and w_min varies from 0.0048 to 0.007. Substituting these values into Equation (25) brings α into the range 0.70–1.08 (see the squares in Figure 10(a); the dashed line is the fit to the results). Figures 10(b) and (c) are for the cases of two blobs approaching one another as they move in the sheet. Figure 10(b) covers the period from t = 31.8 to t = 33.2, M_A ranges from 0.034 to 0.086, λ is from 0.16 to 0.3, and w_min ranges from 0.004 to 0.007. We find that α is between 0.50 and 0.86, and it tends to decrease as two blobs merge with one another. Figure 10(c) shows a situation similar to that in Figure 10(b) and covers the period from t = 33.6 to t = 35.0. In this period, λ is between 0.29 and 0.45, M_A changes in the range from 0.054 to 0.105, and w_min varies from 0.005 to 0.01. Therefore, we deduce α according to Equation (25) such that 0.65 < α < 1.12.

In the above three cases, Figure 10(a) is for the case in which two blobs are moving apart, and the largest distance between them is λ = 0.61, which is the largest one of two adjacent blobs that we can track within the calculation domain. Figures 10(b) and (c) are the example of two blobs merging with one another. In these two cases, we note that we chose smaller values of λ, 0.29 and 0.16, respectively, in an arbitrary fashion, instead of λ = 0 when two of them completely merge together. This is because the turbulence developed by the tearing mode is still in the linear stage prior to the merging of the two blobs and Equation (25) holds. Whether Equation (25) still remains valid as two blobs totally merge or even as they become very close to one another is an open question, and will be investigated in future work.

What we have found so far from three examples indicates that 0.50 < α < 1.12. Our results also show that α tends to decrease with λ. This probably has something to do with the fact that the value of w between two adjacent blobs in the present experiment remains almost unchanged. Consulting the work by Lin et al. (2009) which gave a value of α between 1/7 and 1, we conclude that the range of α values deduced in this work is consistent with theirs, which can be considered a confirmation of their theoretical results by the numerical experiment performed here.