MAGNETOHYDRODYNAMIC SIMULATION OF THE X2.2 SOLAR FLARE ON 2011 FEBRUARY 15. I. COMPARISON WITH THE OBSERVATIONS

To extrapolate the NLFFF, we use the vector magnetic field shown in the left panel of Figure 1(a), observed at 00:00 UT on 2011 February 15, approximately 120 minutes before the X2.2-class flare detected by SDO/HMI. The vector magnetic field covers a 216 × 216 Mm² region, divided by a 600 × 600 grid that is provided by the description of the cylindrical equal area projected and remapped vector magnetic field.⁴ It was obtained using the very fast inversion of the Stokes vector algorithm (Borrero et al. 2011) based on the Milne–Eddington approximation. A minimum energy method (Metcalf 1994; Metcalf et al. 2006; Leka et al. 2009) was used to resolve the 180° ambiguity in the azimuth angle of the magnetic field. We also use extreme-ultraviolet (EUV) images observed with 94 Å and 171 Å taken by Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board SDO and Ca ii H 3968 Å images during the flare taken by Hinode, in order to compare with the NLFFF and MHD simulation. Ca ii images were obtained using the Broadband Filter Imager of the SOT on board Hinode. The field of view corresponds to 11158 × 11158 with a resolution of 1024 × 1024 at 01:50:18 UT on February 15.

Zoom In Zoom Out Reset image size

We carry out the NLFFF and MHD simulation according to the following equations:

$$\begin{equation} \rho = |{\boldsymbol {B}}| \end{equation} \tag{ 1 }$$

or

$$\begin{equation} \frac{\partial \rho }{\partial t} = -{\boldsymbol {\nabla }} \cdot (\rho {\boldsymbol {v}}) + \xi \nabla ^2 (\rho - \rho _0), \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial {\boldsymbol {v}}}{\partial t} = - ({\boldsymbol {v}}\cdot {\boldsymbol {\nabla }}){\boldsymbol {v}} + \frac{1}{\rho } {\boldsymbol {J}}\times {\boldsymbol {B}} + \nu {\boldsymbol {\nabla }}^{2}{\boldsymbol {v}}, \end{equation} \tag{ 3 }$$

$$\begin{equation} \frac{\partial {\boldsymbol {B}}}{\partial t} = {\boldsymbol {\nabla }}\times ({\boldsymbol {v}}\times {\boldsymbol {B}} - \eta _{i}{\boldsymbol {J}}) - {\boldsymbol {\nabla }}\phi, \end{equation} \tag{ 4 }$$

$$\begin{equation} {\boldsymbol {J}} = {\boldsymbol {\nabla }}\times {\boldsymbol {B}}, \end{equation} \tag{ 5 }$$

$$\begin{equation} \frac{\partial \phi }{\partial t} + c^2_{h}{\boldsymbol {\nabla }}\cdot {\boldsymbol {B}} = -\frac{c^2_{h}}{c^2_{p}}\phi, \end{equation} \tag{ 6 }$$

where the NLFFF calculation employs Equation (1) as the density evolution in order to ease the relaxation by equalizing the Alfvén speed in space as shown in Inoue et al. (2013). ρ is the mass density, ρ₀ is an initial value of the density where the density evolution described by Equation (2) follows from Aulanier et al. (2005, 2010), ${\boldsymbol {B}}$ is the magnetic flux density, ${\boldsymbol {v}}$ is the velocity, ${\boldsymbol {J}}$ is the electric current density, and ϕ is a scalar potential. The last equation (Equation (6)) is used to reduce deviation from $\nabla \cdot {\boldsymbol {B}}=0$ and was introduced by Dedner et al. (2002).

The length, magnetic field, density, velocity, time, and electric current density are normalized by L* = 216 Mm, B* = 2500 G, ρ* = |B*|, $V_{A}^{*}\equiv B^{*}/(\mu _{0}\rho ^{*})^{1/2}$, where μ₀ is the magnetic permeability, $\tau _{A}^{*}\equiv L^{*}/V_{A}^{*}$, and J* = B*/μ₀L*, respectively. The non-dimensional viscosity ν is set as a constant (1.0 × 10⁻³), and the coefficients $c_h^2$, $c_p^2$ in Equation (6) also fix the constant values, 0.04 and 0.1, respectively. The subscript i of η corresponds to NLFFF or MHD. These equation sets, normalized values, and parameters are mostly identical to those in previous studies (Inoue et al. 2013), except for the density evolution described in Equation (2).

A formulation of η_NLFFF is given by

$$\begin{equation} \eta _{{\rm NLFFF}} = 5\eta _0 + \eta _1 \frac{|{\boldsymbol {J}}\times {\boldsymbol {B}}||{\boldsymbol {v}}|^2}{\boldsymbol {|B|^2}}, \end{equation} \tag{ 7 }$$

which is slightly different from Inoue et al. (2013), where η₀ = 1.0 × 10⁻⁵ and η₁ = 1.0 × 10⁻³. The vector field shown in Figure 1(a) is preprocessed in accordance with Wiegelmann et al. (2006), and a potential field applied to an initial condition of the NLFFF calculation is extrapolated using the Green function method (Sakurai 1982). The observed tangential components of the magnetic field reside in the region enclosed by a dashed square, as shown in the left panel of Figure 1(a); the outside region is fixed by the potential field because we focused on the dynamics of the magnetic field in an early phase near the central area and constructed a field close to an equilibrium state numerically. Note that, according to the mathematical formula, this state mixing of the force-free field and the potential field cannot keep an equilibrium. Nevertheless, we employed this method in this study in order to avoid an effect from inconsistent force-free α residing outside areas where the weak magnetic fields are dominant.

Next, we perform MHD simulation under the zero β approximation (Mikic et al. 1989 or Amari et al. 1996) due to the low-β (0.01−0.1) coronal plasma where the time evolution of density is according to Equation (1), which is similar to Amari et al. (1996), i.e., the equations employed by the MHD simulation are identical to that of the NLFFF extrapolation method in Inoue et al. (2014) or Inoue et al. (2013), except the velocity limit is not imposed and the resistivity formulation is different.

The initial condition on the density is given by $\rho _0={\boldsymbol {|B|}}_{t=0}$ in all MHD (also NLFFF) calculations. In this case, the density is modeled by Equation (1); therefore, we cannot compare it with observational data exactly or discuss the timescale of the dynamics of the magnetic field, as the Alfvén timescale depends on this model. On the other hand, this density profile might be more reasonable than the constant density profile, e.g., ρ =1.0, because the real coronal density is stratified due to gravity, and the Alfvén wave propagates effectively in the upper corona where the magnetic field is driven. In addition to these, Inoue & Kusano (2006) showed the dynamics of a flux tube using two types of density. One is modeled by ρ = 1.0 at all times, i.e., the density variation is not considered. Another is obtained from a continuous Equation (2). Consequently, both dynamics of the magnetic field are nearly identical at the early eruptive phase. We will discuss this again later in this paper. A resistivity formula is given as an anomalous resistivity as follows:

$$\begin{equation} \eta _{{\rm MHD}} = \left\lbrace \begin{array}{@{}ll}\eta _{0} & \quad {J<j_{c}}, \\ \eta _{0}+\eta _{2} \displaystyle\left(\frac{J-j_{c}}{j_{c}}\right)^{2} & \quad {J>j_{c}}, \end{array}\right. \end{equation} \tag{ 8 }$$

where η₂ = 5.0 × 10⁻⁴ and j_c is the threshold current, set to 30 in this study.

In this simulation, a numerical box with dimensions of 216 × 216 × 216 Mm³ is given by 1 × 1 × 1 as a non-dimensional value. The grid number is assigned as 300 × 300 × 300, which is 2 × 2 binning from the original data. Regarding the boundary conditions of NLFFF (Run A) and MHD relaxation processes (Run C), the physical values of all boundaries are fixed (fixed boundary condition), except for the Neumann-type boundary condition, on which the convenient potential ϕ is imposed. On the other hand, in MHD simulations (Run B, and Runs D–F), tangential components of the magnetic field are released on all of the boundaries (released boundary conditions). The numerical scheme is quite identical to that in Inoue et al. (2014) or Inoue et al. (2013). The kind of simulation (NLFFF or MHD), equation of the density evolution, initial conditions, and boundary conditions used in each run are summarized in Table 1.

3.1.1. 3D Magnetic Structure of the NLFFF

3.1.2. Physical Properties of the NLFFF

Next, we show the physical properties of the NLFFF, focusing in particular on how much the reconstructed field is in a force-free state, its stability, and how much magnetic twist is accumulated in it. Figure 2(a) shows a distribution of the force-free α where the vertical and horizontal axes represent the values of the force-free α estimated on opposite footpoints of each field line that are calculated within the dashed square in Figure 1(a). Like Inoue et al. (2013), these force-free α values are estimated at 720 km approximately, i.e., 1 grid size in this case above the photosphere, and the field lines are traced from the region where B_z < −250 G. This result seems to show slightly larger scatter than that in our previous result (Inoue et al. 2013, ; see Figure 2(a)) even though many points are aligned close to the y = x line plotted in green. This might be due to the fine force-free α compared to the previous one residing in the vector field reduced by an 8×8 binning process.

Zoom In Zoom Out Reset image size

We perform MHD simulations using a reconstructed NLFFF as an initial condition to find a state in the details of the NLFFF where the tangential components of the magnetic field are released on the bottom boundary, employing constant resistivity. These calculations are designated as Run B. Figure 2(b) shows the temporal evolution of the kinetic energy ($\int(\rho |{\boldsymbol {v}}|^2/2)dV$, where dV is a volume element) for Run B. We first observe an increase of the kinetic energy because the NLFFF cannot meet an equilibrium state exactly in addition to the configuration further deviating from the NLFFF during the simulation due to the release boundary condition at the bottom surface. On the other hand, after that, this profile decreases as time goes on, which suggests that the reconstructed field moves back to the potential field even on the disturbances derived from the velocity. Two insets in Figure 2(b) represent the 3D view of the field lines at t = 0 and t = 15, marked by the circles. We clearly see that the twisted field lines at t = 15 become looser than the initial state t = 0, i.e., the NLFFF is in a stable state and never shows the drastic dynamics as seen in the observations.

We also investigate the magnetic twist that is defined here by

$$\begin{equation} T_n = \frac{1}{4\pi }\int \alpha dl, \end{equation} \tag{ 9 }$$

where α and dl are the force-free α(= ${\boldsymbol {J}} \cdot {\boldsymbol {B}}/|{\boldsymbol {B}}|^2$) and a line element of a field line, respectively. We trace the field lines starting from the photospheric region with |B_z| > 0.2 (500 G) within the dashed square in the left panel of Figure 1(a), in order to focus on the twisted lines in the strong magnetic field. Figure 2(c) shows the 3D field lines over the B_z distribution. Orange lines represent the strongly twisted lines over the half-turn twist that are traced inside of the red contours (T_n = 0.5) shown in Figure 2(d). Figure 2(e) represents the twist distribution on the positive polarity. This result indicates that most of the twists have a value of less than one turn, suggesting that they are stable against the kink mode instability. These results are quite similar to our previous result in Inoue et al. (2013) even though some parts of the NLFFF model are different.

From the above section showing the physical properties of the NLFFF, we found that the NLFFF never shows dramatic dynamics or eruption away from the solar surface, which is contrary to the observation (Schrijver et al. 2011). Therefore, there might be a mechanism that excites the stable NLFFF into a dynamically unstable phase, causing the observed eruptive phenomena. Thus, we perform the MHD relaxation process using the NLFFF as an initial condition, which is similar to the NLFFF calculation in Run A, except a velocity limit is released and an anomalous resistivity described in Equation (8) is set in this calculation, in order to form further strongly twisted lines with more than a one-turn twist. Because the anomalous resistivity plays a role in enhancing the reconnection in the strong current region, we can expect it to generate through reconnection more strongly twisted lines than that formed in the NLFFF. Because the magnetic field is fixed at the bottom boundary during this calculation, more strongly twisted lines are formed without changing the original horizontal components on the vector field, i.e., the map of the observed force-free α on the vector field is not changed. We call this calculation Run C.

Figure 3 shows the temporal evolution of the kinetic energy for Run C. We clearly see that the kinetic energy starts to accelerate after t = 0.5 and then increases as time goes on. The result in Run C indicates that the solution is deviating from the NLFFF, suggesting that it converts into the dynamic phase through the reconnection due to the anomalous resistivity.

Zoom In Zoom Out Reset image size

In order to confirm this suggestion that we note in the previous section, we observe the 3D magnetic structure and estimate the twist of the magnetic field lines at t = 1, marked by the vertical dashed line in Figure 3, which is just beginning to convert into the dynamic phase. Figures 4(a) and (b) show the 3D view of the magnetic field at t = 0, i.e., NLFFF and t = 1. The twisted lines in orange at t = 1 look slightly different from the initial NLFFF state, even though the overlying field lines surrounding these twisted lines almost keep their original configurations during this period in Run C. We further see the distribution of the magnetic twist and 3D magnetic structure in Figures 4(c) and (d) to understand the more detailed magnetic structure close to the PIL in the central area. Figure 4(c) shows that most of the twist values in NLFFF are less than one turn as discussed in the above section; on the other hand, the strongly twisted lines with more than one turn are generated at t = 1 as shown in Figure 4(d), indicating the locations of their footpoints marked by dashed circles in the upper panel, and showing the field lines traced from there in the lower panel. Interestingly, the locations of these footpoints are close to those of the sunquake detected by helioseismology (Zharkov et al. 2011), which was triggered at the footpoint of the erupting flux rope at an early flare phase. In addition, the change of the field line topology from t = 0 to t = 1 is reminiscent of the tether-cutting reconnection (Liu et al. 2012; Wang et al. 2012; Liu et al. 2013).

Zoom In Zoom Out Reset image size

Therefore, this result might suggest that the tether-cutting reconnection plays an important role in breaking the equilibrium state before the flare.

3.3.1. 3D Dynamics of Eruptive Magnetic Fields

We suggested in the above section that the magnetic field with a strongly twisted line over one turn can break the equilibrium condition and convert the NLFFF into the dynamic phase. However, its time evolution is not consistent in that simulation because the three components of the magnetic field are fixed on the boundaries, that is clearly inconsistent with the boundary condition to Equation (4), meaning that the dynamics is not shown properly. In this section, we therefore perform the MHD simulation using the magnetic field shown in Figure 4(b), which is breaking the equilibrium state, as an initial condition where the tangential components are released on all the boundaries to discuss its time evolution. This calculation is called Run D.

We first show an overview of the 3D dynamics of the magnetic field lines for Run D with vertical velocity (v_z) in Figure 5. This result shows that the twisted field lines with more than a half-turn twist, being formed at t = 1 in Figure 4(b) (in orange), erupt away from the solar surface even though the tangential components of the magnetic field are released on all the boundaries. Because of the release boundary condition, the tangential components are deviating from that of the original vector field over time, which is relaxing toward them close to the potential field. The NLFFF obtained from Run A never shows these dynamics; this result thus indicates that the strongly twisted lines over the one-turn twist formed in the NLFFF are required in this AR to break the quasi-equilibrium state.

Zoom In Zoom Out Reset image size

3.3.2. Comparison with Observational Data

We compare these simulation results with observational data to confirm the reliability of our simulation. Figure 6(a) shows the EUV images in 171 Å taken by SDO/AIA. Figure 6(b) plots the contour of |B_z| = 0.25 over the AIA images. EUV is strongly enhanced at the PIL between the sunspots located in the central area; on the other hand, we can also observe another enhanced location, marked by a red circle, away from the central area. This location corresponds to the edge of the negative polarity on the east side in Figure 6(b). From Figure 5, we can see that one footpoint of the twisted field lines also jumps from the central area to the edge of the negative polarity after t = 2.5, as suggested from Figure 6. Therefore, these observational results might support our simulation.

Zoom In Zoom Out Reset image size

Next we compare our simulation results with the Ca ii image taken by Hinode. We focus on investigating the reconnection of only the twisted lines with more than one-third turn (T_n = 0.3) because Inoue et al. (2011, 2012, 2013) suggested that these twisted lines close to the PIL are related to the two-ribbon flares where the Ca ii image is strongly enhanced. Furthermore, in this study we focus on the early eruptive phase as shown in Figure 7(a) because in the late phase the tangential components of the magnetic field on the bottom boundary are deviating gradually from the observational one due to the released boundary condition.

Zoom In Zoom Out Reset image size

Following Toriumi et al. (2013), the reconnected field lines are estimated using the spatial variance of the field line connectivity, which is allowed only by the magnetic reconnection. The formulation is defined by the following equation:

$$\begin{equation*} \delta ({\boldsymbol {x}}_0,t_n) = |{\boldsymbol {x}}_{1}({\boldsymbol {x}}_0, t_{n+1}) - {\boldsymbol {x}}_{1}({\boldsymbol {x}}_0, t_{n})| \quad ({|T_n| \geqq 0.3}), \end{equation*}$$

where ${\boldsymbol {x}}_{1}({\boldsymbol {x}}_0, t_{n})$ is a location of one footpoint of each field line at time t_n, which is traced from another footpoint at ${\boldsymbol {x}}_0$. Eventually, we calculate

$$\begin{equation} \Delta ({\boldsymbol {x}}_0,t) = \int _{0}^{t} \delta ({\boldsymbol {x}}_0,t_n)dt_n, \end{equation} \tag{ 10 }$$

where the enhanced region in $\Delta ({\boldsymbol {x}}_0, t)$ indicates memories in which the reconnection took place dramatically in the twisted lines. The left panel in Figure 7(b) exhibits the Ca ii image taken by Hinode/SOT at the early flare phase, and the middle and right panels show the distribution of $\Delta ({\boldsymbol {x}}_0)$ at t = 2.0 and t = 4.0, respectively. Interestingly, the distribution of $\Delta ({\boldsymbol {x}}_0)$ is similar to that of the observed two-ribbon flares seen in the left panel of Figure 7(b). This result clearly shows that the two-ribbon flare is enhanced due to the reconnection of twisted lines, which is consistent with classical flare models (see Shibata & Magara 2011).

Finally, we show the field line structure after the flare. Figure 8(a) shows the field lines over the flare ribbons obtained from the MHD simulation from which all of the field lines are traced. Our simulation shows the post-flare loops on the flare ribbons as often seen in observed EUV or Ca ii images. We compare the simulation results with the EUV image with 94 Å taken by SDO/AIA shown in Figure 8(b). Before showing the comparison with the result of the MHD simulation, Figure 8(c) shows the NLFFF after the flare, which is reconstructed from the vector field at 03:00 UT on February 15, plotted over the EUV image. The NLFFF almost captures the post-flare loops in the EUV image shown in Figure 8(b). On the other hand, the field line structure at t = 10 obtained from the MHD simulation slightly deviates from the observed image and looks like a relaxation compared to those of the NLFFF because the tangential components of the magnetic field are released on the bottom boundary, relaxing toward the state close to the potential field. Nevertheless, they roughly capture the EUV image, and this result would also support the reliability of our simulation.

Zoom In Zoom Out Reset image size

In the previous section, the MHD simulation showed that the twisted lines initially embedded in the lower corona successfully erupt from the solar surface. However, it is still unclear whether these twisted lines grow to become coronal mass ejection (CME) or not. In order to clarify this, we discuss the decay index profile related to the torus instability (Kliem & Török 2006), which is an important factor in determining whether the flux tube can escape from the lower corona or not. The decay index refers to an estimation of a criterion for torus instability and is given by the following equation:

$$\begin{equation} n(z) = - \frac{z}{|{\boldsymbol {B}}|}\frac{\partial |{\boldsymbol {B}}|}{\partial z}, \end{equation} \tag{ 11 }$$

where the threshold is at approximately n_c = 1.5 (Török & Kliem 2007). Because it is difficult to separate an external field from the NLFFF, following Fan & Gibson (2007) or Aulanier et al. (2010), the decay index on the potential field is plotted in Figure 9(a), marked by a cross in the inset, where the strong twisted lines of the NLFFF are embedded. From this result, the threshold is located at approximately h = 0.2, corresponding to 43.2 Mm. Figure 9(b) shows the time evolution of selected twisted field lines with more than a half turn where at t = 0.0 the twisted line exists in the lower area corresponding to less than h = 0.2 at the top of the box. This indicates that the twisted line is stable against torus instability. On the other hand, as time goes on, these cross over the critical height as a result of the initial eruption. This result might give an important scenario related to the dynamics of the flare CME.

Zoom In Zoom Out Reset image size

We present a more detailed temporal evolution of the same selected twisted lines in Figure 9(c), projecting them at each time on a 2D plane with a red line tracing the temporal evolution of the loop top. The dashed line indicates the critical height of the torus instability. Figure 9(d) shows the time variation of the loop top, i.e., its ascending velocity upon appearance, with the vertical dashed line where the twisted line is passing marking the critical height. This result does not exhibit the sudden acceleration of the twisted line in this study, as shown in Aulanier et al. (2010) and Fan (2010), even though the top of the twisted line exceeds the critical height and its velocity increases slowly after t = 5.0.

We point out two possible causes of this problem. One is that the whole box is not large enough to trace the twisted lines for a long time. In addition, we set the rigid wall, which is not able to expel the overlying field lines to the outside of the numerical domain, at all boundaries; consequently, the magnetic pressure accumulates close to the boundary because the twisted lines keep pushing the overlying field lines. Another is due to the bottom boundary condition, i.e., due to the release of the tangential components. Torus instability is driven by the poloidal field of the flux tube generated by the toroidal current in it. However, the release condition for the tangential components relaxes the twisted lines, which means that the hoop force driving the twisted lines upward is loosening gradually. These two factors preventing the twisted lines from accelerating might work effectively on their dynamics. In other words, these twisted lines might have the potential for a solar eruption if settling these problems on the size of the numerical box and the boundary condition. Therefore, for future works, it is important to understand how to bring the three observed components of the magnetic field as a boundary condition into an induction equation, i.e., we would require a data-driven simulation as presented in Cheung & DeRosa (2012).

We check the dynamics of the density formula dependence by replacing Equation (1) by Equation (2) where ξ is given by 1.0 × 10⁻⁴. This calculation is called Run E. Figure 10(a) shows the temporal evolution of the kinetic energy for Run D and Run E, respectively. Both of them have a lot of similar profiles, except we can see a little bit of deviation at the late phase. On the other hand, Figure 10(b) shows the temporal evolution of the iteration number, which shows that Run E requires much more calculation time than Run D. The 3D field line structures marked by each circle in Figure 10(a) are shown in Figures 10(c) and (d). As a consequence, we conclude that the dynamics is almost the same between them, but the cost performance of Run D is much better than that of Run E, as shown in Inoue & Kusano (2006). This conclusion is also in agreement with the previous work by Amari et al. (1996) and their series of papers.

Zoom In Zoom Out Reset image size

We further perform an MHD calculation (Run F) replacing the anomalous resistivity with a constant one to also show the dynamics of the dependence on the resistivity formula. Figure 11(a) shows the temporal evolution of the maximum value of the current density (=$|{\boldsymbol {J}}|_{{\rm max}}$) measured above 3600 km for Run D and Run F, respectively, which can exclude the strong current density existing close to the photosphere. Although the value of Run D initially undergoes a sudden decrease and remains less than the one in Run F during the calculation, this sudden decrease is probably due to the diffusive effect derived from the anomalous resistivity. On the other hand, because the critical current set in the anomalous resistivity corresponds to a value of 30, it works by t = 3 approximately. Nevertheless, temporal evolutions of the kinetic energy for Run D and Run F shown in Figure 10(b) have quite similar profiles, although we can see a little bit of deviation between them. Furthermore, 3D field line structures and color maps of $|{\boldsymbol {J}}|$ for Run D and Run E shown in Figures 10(c) and (d) are almost the same, too, even at t = 5.0. From these results, in this case the anomalous resistivity does not work effectively on the dynamics while the twisted lines are ascending. Note that the reconnection is very important in the early phase to form the strongly twisted lines in the lower corona that launch from the solar surface. However, after that launching, the reconnection would take place only accompanied by the ascension of the twisted lines but not play an important role in the acceleration of the twisted line, as pointed out by Inoue & Kusano (2006).

Zoom In Zoom Out Reset image size