Magnetohydrodynamic Simulations for Studying Solar Flare Trigger Mechanism

In this work, we used AR NOAA 10930 as a case study for our simulation. This AR was bipolar, with the negative polarity spot larger than the positive. It was very active, having produced at least 113 X-ray flares of different energy from 2006 December 4 to 2006 December 18 (Gopasyuk 2015). Here, we focus on the X3.4 class solar flare that occurred at 02:14 UT on 2006 December 13. Many studies of this AR have been extensively conducted on various aspects, i.e., sheared field (Kubo et al. 2007; Su et al. 2007), helicity and twist (Magara & Tsuneta 2008; Su et al. 2009; Inoue et al. 2011), rotating sunspot (Min & Chae 2009; Gopasyuk 2015), NLFFF extrapolation (Schrijver et al. 2008; Inoue et al. 2012), and MHD simulation (Fan et al. 2011; Amari et al. 2014). Sigmoidal structure has been reported by observations as well as NLFFF extrapolation (Min & Chae 2009; Inoue et al. 2012; Amari et al. 2014).

We used vector magnetic field of AR 10930, derived from Spectro Polarimeter (SP) data from the Solar Optical Telescope (SOT) instrument (Tsuneta et al. 2008), onboard the Hinode satellite (Kosugi et al. 2007) for the NLFFF extrapolation. We used Ca ii H line (3968.5 $\mathring{\rm A}$) data from the Broadband Filter Instrument in the SOT to examine how well the NLFFF and MHD simulation results agree with the structure of the flare ribbons. The X-ray image of the AR 10930 was obtained from the X-ray Telescope (XRT) on-board Hinode (Golub et al. 2007). This image is important to compare the NLFFF results with the coronal magnetic field configuration inferred from the X-ray image.

We inserted vector magnetogram data obtained from the Hinode/SP as a bottom boundary condition from the original 1000 × 512 pixels, in order to fit to the $240\times 128\times 128$ uniform grid used in the simulation box. The magnetogram's field-of-view is 297 × 163 arcsec, corresponding to 214 × 118 Mm on the Sun. The simulation box represents the rectangular domain of (−0.5L, −0, 25L, 0) ≤ (x, y, z) ≤ (0.5L, 0.25L, 0.5L), where L is the normalization of spatial length, which has an actual value of about 214 Mm.

We follow the MHD relaxation method of Inoue et al. (2014b) to reconstruct the coronal field of the AR we are interested in. We use vector magnetic field data obtained from the Hinode/SP magnetogram at 20:30 UT on 2006 December 12, which was about six hours before the X3.4 flare onset. The potential field of the AR is calculated as an initial condition from the normal component B_z of the vector magnetic field on the photosphere, using the Fourier method (Alissandrakis 1981). The initial density is chosen to be uniform. After inserting the observed tangential components (B_x and B_y) into the bottom boundary, the magnetic field in the whole domain is then evolved toward the force-free state. This evolution process is governed by a set of equations for zero plasma beta,

$$\begin{eqnarray}&&\rho ={\rho }_{0}\displaystyle \frac{| {\boldsymbol{B}}| }{{B}_{0}},\end{eqnarray} \tag{ 1 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-{\boldsymbol{\nabla }}\times (-{\boldsymbol{v}}\times {\boldsymbol{B}})+{\eta }_{{}_{\mathrm{NLFF}}}{{\boldsymbol{\nabla }}}^{2}{\boldsymbol{B}}-{\rm{\nabla }}\phi ,\end{eqnarray} \tag{ 3 }$$

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

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

where ρ is the plasma density,  v is the plasma velocity, J is the current density, and B is the magnetic flux density. In this method, Equation (1) defines a pseudo-density (ρ) that is proportional to $| {\boldsymbol{B}}|$, in order to ease the relaxation by maintaining the Alfvén speed in space (Inoue et al. 2013). Equation (2) is the equation of motion for the zero plasma beta condition neglecting gravity. The last term in the induction Equation (3) includes the ${\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}$ cleaning potential (ϕ). The cleaning potential Equation (4) was introduced by Dedner et al. (2002) to reduce deviation from the solenoidal condition ${\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}=0$, where c_h and c_p are the coefficients related to advection and diffusion of ${\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}$, respectively.

The magnetic field (B) in the calculation is normalized by B₀, which equals 4000 G. Velocity, time, and electric current density are normalized by ${V}_{A}\equiv {B}_{0}/{({\mu }_{0}{\rho }_{0})}^{(1/2)}$, ${\tau }_{A}\equiv L/{V}_{A}$, and ${J}_{0}={B}_{0}/{\mu }_{0}L$, respectively. In the typical AR, ${\rho }_{0}=1.67\,\times 10$⁻¹² kg m⁻³, so that ${V}_{A}\approx 275$ Mm s⁻¹, ${\tau }_{A}\approx 0.8\,{\rm{s}}$, and ${J}_{0}\approx 1.5$ mA m⁻². We set the coefficients following Inoue et al. (2014b), where c_p² and c_h² are 0.1 and 0.04, respectively. The non-dimensional viscosity (ν) in Equation (2) is set as a constant ($1.0\times {10}^{-3}$). Magnetic diffusivity (η) in Equation (3) is defined as

$$\begin{eqnarray}&&{\eta }_{{}_{\mathrm{NLFF}}}={\eta }_{{}_{0}}+{\eta }_{{}_{1}}\displaystyle \frac{| {\boldsymbol{J}}\times {\boldsymbol{B}}| | {\boldsymbol{v}}{| }^{2}}{| {\boldsymbol{B}}{| }^{2}},\end{eqnarray} \tag{ 6 }$$

where ${\eta }_{{}_{0}}=5.0\times {10}^{-4}$ and ${\eta }_{{}_{1}}=1.0\times {10}^{-3}$ are non-dimensional parameters in the units of ${\mu }_{0}{V}_{A}L$ and ${({\mu }_{0}L)}^{2}/{V}_{A}$, respectively.

At the bottom boundary, once we run the program, the tangential components from the potential field are incrementally changed into the observed tangential components, whereas all physical values in the other boundaries are fixed. After the bottom boundary values of magnetic vector field are completely changed into the observed values, we set all the physical values for all boundaries (including the bottom) to be fixed during the calculation. The method for the NLFFF extrapolation and parameter setting in this work are almost identical to the NLFFF method of Inoue et al. (2014a).

The MHD simulation is performed in the same grid as the NLFFF extrapolation. It uses the NLFFF model and corresponding density as its initial conditions. The non-ideal zero-beta MHD equations are solved in the MHD simulation. Hence, the induction equation now takes the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-{\boldsymbol{\nabla }}\times (-{\boldsymbol{v}}\times {\boldsymbol{B}}+{\eta }_{{}_{\mathrm{MHD}}}{\boldsymbol{J}}),\end{eqnarray} \tag{ 7 }$$

and the continuity equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}}),\end{eqnarray} \tag{ 8 }$$

replaces Equation (1). The magnetic diffusion (η) in Equation (7) is defined as an anomalous resistivity following Inoue et al. (2014a),

$$\begin{eqnarray}{\eta }_{{}_{\mathrm{MHD}}}(t)=\left\{\begin{array}{ll}{\eta }_{{}_{2}}, & J\leqslant {j}_{c}\\ {\eta }_{{}_{2}}+{\eta }_{{}_{3}}\left(\tfrac{J-{j}_{c}}{{j}_{c}}\right), & J\gt {j}_{c},\end{array}\right.\end{eqnarray} \tag{ 9 }$$

where ${\eta }_{{}_{2}}=1.0\times {10}^{-5}$, ${\eta }_{{}_{3}}=5.0\times {10}^{-3}$, and the threshold current density, ${j}_{c}=300$. This anomalous resistivity can be expected to enhance the reconnection of the field lines in the regions of strong current (Inoue et al. 2014a).

Kusano et al. (2012) suggested that the trigger structure is located near the photospheric PIL. Accordingly, we expect the area near the PIL of the core field to be particularly effective for triggering a flare. According to Bamba et al. (2013), the trigger structure of the X3.4 flare studied here was situated in the area marked by the yellow circle in Figure 1(a). They showed that a highly sheared structure existed along the PIL, and a small positive polarity magnetic island grew near the PIL, as is shown in Figure 1(b). This location was obtained from their study of the topological features of the flare ribbons and their associated highly sheared structure. They found that the emerging flux of the magnetic island, which was located between the flare ribbons, triggered the X3.4 flare six hours after this magnetogram was taken. The orientation of the bipole flux was opposite to that of the large-scale magnetic field of the AR. Thus, it could lead to the eruption of the sheared or twisted magnetic field lines by introducing a reconnection that formed and destabilized the flux rope. Therefore, we chose this location as the place where the small bipole structure is injected as flux that emerges into the initial field in our simulation, to trigger the eruption.

Zoom In Zoom Out Reset image size

The emerging flux model follows the method of Kusano et al. (2012), where the small bipole is made from a magnetic torus that ascends from below the simulation box. The bipole structure is a sphere with radius r_e filled with a purely toroidal field of uniform strength B_e. An electric field ${{\boldsymbol{B}}}_{e}\times {{\boldsymbol{v}}}_{e}$ is imposed in the cross-section of the bottom plane to let the torus ascend with velocity v_e, chosen to be constant during the period $0\leqslant t\leqslant {\tau }_{e}(={r}_{e}/{v}_{e})$. The injected bipole structure has the azimuthal orientation angle ${\phi }_{e}$, defined as shown in Figure 1(c). The bipole is injected at the coordinate $P(x=294,y=-98)$ arcsec; it has magnetic intensity ${B}_{e}\,=15$ and radius ${r}_{e}=0.01$. It starts to ascend with the constant velocity ${v}_{e}=0.02$ at t = 0, and is stopped at $t={\tau }_{e}=0.5$ when the center of the sphere reaches the bottom plane. This velocity is higher than the typically observed photospheric velocities, but still slower than the coronal Alfvén velocity. Therefore, it is appropriate for the problem studied here. We perform simulations with various angles ${\phi }_{e}$. Eight cases are run as summarized in Table 1. Case C (${\phi }_{e}=110^\circ$) is displayed in Figure 1(d).

Table 1.  Azimuth Angle (${\phi }_{e}$) and Total Reconnected Flux (${{\rm{\Phi }}}_{\mathrm{rec}}$) in the Area of the Red Square in Figure 1 for Different Cases Performed in the Simulations, Estimated from the Flux Covered by Field Lines with a Large Displacement ${\rm{\Delta }}x({x}_{0})$

Run Case	Orientation (${\phi }_{e}$)	${{\rm{\Phi }}}_{\mathrm{rec}}$ a
A	100	2.44 × 10−4
B	500	4.68 × 10−4
C	1100	6.2 × 10−4
D	1350	6.27 × 10−4
E	1800	8.12 × 10−4
F	2250	4.81 × 10−4
G	2700	2.14 × 10−4
H	3150	1.29 × 10−4

Note.

^aThe values are normalized by ${{\rm{\Phi }}}_{0}\equiv {B}_{0}{L}^{2}=1.83\times {10}^{24}$ Mx.

Download table as:  ASCIITypeset image

Based on its orientation with respect to the pre-existing field, the bipole configurations imposed in our simulations can be classified as right polarity (${\phi }_{e}\approx 0^\circ$), RS (${\phi }_{e}\approx 90^\circ$), OP (${\phi }_{e}\approx 180^\circ$), and normal shear (${\phi }_{e}\approx 270^\circ$) types, using the terms introduced by Kusano et al. (2012). RS-type configuration is defined as the bipole flux whose orientation is almost oppositely directed to the shear (non-potential) field component. Here, we define the RS type to be the bipole with ${\phi }_{e}\approx 90^\circ$, because the shear field in the area around the PIL has a left-handed twist such that the magnetic helicity is negative. The left-handed shear and twist can be seen from the reverse S shape of the sigmoid and the angle between the threads of the sigmoid and the PIL. This was confirmed by a computation of twist map by Inoue et al. (2012). OP type, on the other hand, is defined as the bipole structure with the orientation almost opposite to the averaged potential field.

The constraint for the tangential components of magnetic field on the top and bottom boundaries is set to be released during the simulations, whereas the normal components are fixed except for the area where the bipole flux is emerging. At the side boundaries, all physical values are fixed during the simulations. Due to the relatively small size of the numerical box and fixed side boundary conditions, we cannot expect that the simulation will produce a large expansion of the field such as coronal mass ejection. We only focus on the dynamics of the beginning phase of the flare process.

The top view of field lines in the NLFFF model, plotted over normal component of the magnetogram data, is shown in Figure 2(a). It shows that open magnetic field dominates the AR in the area within and surrounding the negative polarity. This is due to the imbalance of the flux between the negative and positive polarities in the AR. The coronal magnetic field is closed in the area surrounding the PIL. We call this the core field of the AR. As shown in Figure 2(a), the core field shows a strong shear, as can also be seen in the photosphere (Figure 1(b) and (d)).

Zoom In Zoom Out Reset image size

The extrapolated NLFFF is strongly sheared, particularly on the lower part of the corona. This core field may contain a large amount of free energy, because it differs strongly from the potential field. A comparison with X-ray images taken by the XRT instrument onboard Hinode (Figure 2(b)) shows that the NLFFF model agrees well with the observation, in terms of the presence of high shear at the PIL. Moreover, the NLFFF model implies that the sigmoidal structure at the PIL consists of short arcade-type field lines (Figure 2(a)). This sigmoidal structure is important because it shows that the core field of the AR is highly sheared (Su et al. 2007; Min & Chae 2009). The reverse-S shape is well reproduced in the NLFFF extrapolation used here, as well as in several previous works by Inoue et al. (2012) and Amari et al. (2014).

As a reference case, we first carry out a simulation without imposing any external perturbation. This simulation is performed to show the nature of the system if there is no emerging flux imposed. The residual Lorentz force in the NLFFF model is verified to be too weak to trigger an eruption (Figures 3(a) and (b)). The simulation shows that the shear of the magnetic field slightly weakens under the condition of the released tangential field components in the bottom boundary, as shown in Figure 3(b). It is easy to understand that the magnetic field will naturally relax to a lower-energy state, which is closer to the potential field configuration. This verifies that any eruption of the NLFFF must be driven by an external disturbance. The relaxation to a stable equilibrium state can also be seen from the evolution of the kinetic energy in the box in Figure 3(c), where a brief initial rise (due to the residual Lorentz force in the initial condition) is followed by a monotonic decrease (similar to Run B in Inoue et al. 2014a).

Zoom In Zoom Out Reset image size

We find that several configurations of the trigger structure can lead to an eruption; these are the structures in cases C, D, E, and F. However, each type of triggering structure creates a different dynamic and topology in the erupting flux rope. Here, we carefully analyze the eruption of the flux ropes in our simulations to clarify the typical dynamics involved in the erupting process. Based on the relation between flare reconnection and flux rope formation, all eruptive cases in our simulation can be categorized into two distinct groups: eruption-induced reconnections and reconnection-induced eruptions. The former is the case when the flux rope is formed before the flare reconnection occurs. In this case, the emerging bipole flux triggers the creation of an unstable flux rope through pre-flare reconnection. Subsequently, the flare reconnection is generated below the flux rope during its eruption. As for the latter, the role of the emerging bipole is to trigger the reconnection between pre-existing magnetic field by reducing the shear of the overlying field, which then creates an unstable flux rope. These two types of dynamical process were also observed in the simulations conducted by Kusano et al. (2012).

Eruption-induced reconnection features are clearly observed in the simulation results for Case E (${\phi }_{e}=180^\circ$) and Case F (${\phi }_{e}=225^\circ$). However, although they share some common features of eruption-induced reconnection, each case has its own topological structure characteristics, due to the difference in the azimuthal angle of the bipole. Figure 4(a)–(d) shows the evolution of the eruption-induced reconnection from the initial condition until the expansion of sigmoidal flux ropes. At the beginning (Figure 4(a)), the magnetic field lines do not form any large twisted flux ropes. After the small bipole structure emerges in the photosphere, it starts to reconnect with the pre-existing field and forms a flux rope with sigmoidal shape (Figure 4(b)). The flux rope then acquires a higher twist (Figure 4(c)) and expands outward (Figure 4(d)). Figure 5(a)–(d) shows the detailed dynamics of the eruption-induced reconnection, from a different point of view. When the bipole flux emerges, a flux rope occurs immediately after the reconnection of the magnetic field near the PIL via the OP-type structure, and high electric current regions (shown by red shade) are formed (Figure 5(b)). This reconnection tends to create a large flux rope with high twist (yellow lines in Figure 5(c)) that quickly erupts, as can be seen in Figure 5(d), where the flux rope lifts the overlying field and finally induces flare reconnection below the flux rope.

Zoom In Zoom Out Reset image size

In Cases C and D, the process of eruption follows reconnection-induced eruption scenarios found in Kusano et al. (2012). Figure 6(a)–(d) shows the evolution of the reconnection-induced eruption for Case C. At the first step, the core magnetic field near the PIL reconnects with the RS structure bipole (Figure 6(b)). Subsequently, the overlying field at the center of the RS structure collapses (Figure 6(c)) and starts to create a large flux rope (Figure 6(d)). Figure 7(a)–(d) shows the detailed dynamics of flux rope formation in the reconnection-induced eruption, from a different point of view. Just after the bipole emerges, some of the pre-existing field lines near the PIL (blue lines in Figure 7(b)) are in contact with the bipole field and create a current sheet. This reconnection reduces the sheared field (blue lines), which causes the overlying field to collapse toward the center. Part of the collapsed field finally reconnects with the bipole structure (green lines in Figure 7(c) and (d)). However, the higher overlying field tends to form a flux rope that then erupts upward, as shown by the yellow lines in Figure 7(d). The topological structure of these steps in the present simulation is consistent with the previous simulation by Kusano et al. (2012) although it is more difficult to observe than in the previous simulation. This is due to the real coronal magnetic field having a more complex structure than the symmetric one seen in the initial boundary condition of LFFF in the Kusano et al. (2012) simulation. All of these steps are also observed in Case D.

Zoom In Zoom Out Reset image size

From the simulation results, we find that the azimuthal angle of the bipole structure plays a very important role in determining the overall dynamics of the magnetic field. Some imposed bipole structures clearly are not effective in triggering a flare, whereas some others are very effective in triggering a flare, with some aspects of the evolution also depending on the azimuth angle. A summary of the simulation results is shown in Figure 8. The events that do not show any eruptive characteristics are the cases where the imposed emerging flux is oriented relatively parallel to the potential or shear (non-potential) components of the pre-existing magnetic core field at PIL. In our simulations, these are the cases A, B, G, and H, with ${\phi }_{e}=10^\circ$, 50°, 270°, and 315°, respectively. Otherwise, an eruption is triggered.

Zoom In Zoom Out Reset image size

Flare ribbons can represent the topology of the reconnecting magnetic field so well that they can be used to check the results of the simulations. The ribbons mark the footpoints of magnetic field lines that reconnect during the flare. A proxy for flare ribbons in the simulations is made by following the method introduced by Toriumi et al. (2013), and applied also in Inoue et al. (2014a). We calculate the total displacement of the footpoint for each field line for a given time and consider field lines with a large footpoint displacement to be reconnected ones. We trace each magnetic field line from each point $({{\boldsymbol{x}}}_{0})$ in the bottom plane and identify the end point of the field line ${{\boldsymbol{x}}}_{1}({{\boldsymbol{x}}}_{0},{t}_{0})$. Here, the end point position ${{\boldsymbol{x}}}_{1}$ as a function of start point ${{\boldsymbol{x}}}_{0}$ and time t_n is denoted as ${{\boldsymbol{x}}}_{1}({{\boldsymbol{x}}}_{0},{t}_{n}$). After some time, we trace again the end point of the field line $({{\boldsymbol{x}}}_{1}({{\boldsymbol{x}}}_{0},{t}_{n+1}))$ from the start point $({{\boldsymbol{x}}}_{0})$. Therefore, the displacement of the end point position for one start point ${{\boldsymbol{x}}}_{0}$ is given by

$$\begin{eqnarray}&&\delta ({{\boldsymbol{x}}}_{0},{t}_{n})=| {{\boldsymbol{x}}}_{1}({{\boldsymbol{x}}}_{0},{t}_{n+1})-{{\boldsymbol{x}}}_{1}({{\boldsymbol{x}}}_{0},{t}_{n})| .\end{eqnarray} \tag{ 10 }$$

By integrating δ for a given time $t={t}_{N}$, we can obtain the total displacement of the end point from the initial state to an arbitrary time step,

$$\begin{eqnarray}&&{\rm{\Delta }}x({x}_{0},t)=\displaystyle \sum _{n=0}^{N}\delta ({x}_{0},{t}_{n}).\end{eqnarray} \tag{ 11 }$$

We assume that the high value of total displacement (${\rm{\Delta }}x$) is due to reconnection. Synthetic flare ribbons are constructed by plotting the value of the total displacement of all field lines as a function of their footpoint position on the bottom plane. The results are shown in Figure 9 for eight different cases. In order to emphasize the distribution of the flare ribbons, we only plot footpoints with total displacement exceeding 0.04. We also calculate the total reconnected flux of B_z from the areas within the red square in Figure 1(a) that show a large displacement of field lines, and denote this as ${{\rm{\Phi }}}_{\mathrm{rec}}$ in Table 1. These total fluxes show how much flux reconnects when the emerging flux meets the condition of a flare trigger.

Zoom In Zoom Out Reset image size

It is obvious from Figure 9 that the flare ribbons constructed from the simulations are different in each case. This suggests that the topology of the reconnected magnetic field is unique for each case. Therefore, the dynamics of the magnetic field due to the interaction between the pre-existing magnetic field and the emerged bipole structure depend considerably on the orientation of the bipole. The results also show that the topology of the magnetic field involved in the eruption process triggered by the RS- and OP-type structures differs, even though both cases result in eruption. The RS-type structure (Case C and Case D) tends to produce local flare ribbons in the area of the core-field of the AR. On the other hand, the OP-type structure (Case E) can generate more extended flare ribbons of more complex structure.