Data-constrained Magnetohydrodynamic Simulation of a Long-duration Eruptive Flare

Magnetic flux rope eruptions in the solar atmosphere are closely related to a series of solar phenomena, such as solar flares, coronal mass ejections (Chen 2011), and filament/prominence eruptions. The important physical processes include magnetic reconnection, particle acceleration, shocks, and so on. These eruptive phenomena could cause severe space weather when the eruptions are directed at the Earth. To prevent and reduce space weather disasters, we need to study the eruption mechanisms and forecast the eruptions. An effective way is to develop numerical models for solar eruptions and uncover these physical processes.

The mechanisms of magnetic flux rope eruptions are usually ascribed to magnetohydrodynamic (MHD) instabilities (Lin & Forbes 2000; Kliem & Török 2006; Démoulin & Aulanier 2010; Keppens et al. 2019) and magnetic reconnection (Antiochos et al. 1999; Moore et al. 2001). These two kinds of mechanisms can both be described by the MHD equations, specifically those for momentum conservation and magnetic induction. The difference between the two mechanisms relates to the different roles of the Lorentz force in the momentum equation and the dissipation term in the magnetic induction equation in initiating the eruptions. In the scenario of MHD instabilities, it is the Lorentz force that drives the system to further deviate from equilibrium. In the scenario of magnetic reconnection, however, magnetic reconnection initially changes the magnetic field connectivity, generates new Lorentz forces, lifts magnetic structures to a new equilibrium height, and then facilitates MHD instabilities. The connection between the two scenarios is that reconnection may also arise from resistive MHD instabilities, which may play a role in both linear and nonlinear MHD evolution.

Apart from some specific processes such as particle acceleration that need a treatment of particle dynamics, the model described by the MHD equations is good enough for studying magnetic flux rope eruptions. The key is to set up an initial condition as close as possible to the real case, which is usually a nonlinear force-free field (NLFFF) model, and to provide an appropriate boundary condition. Such an MHD simulation is called data-driven (Cheung & DeRosa 2012; Jiang et al. 2016; Guo et al. 2019a; Hayashi et al. 2019; Pomoell et al. 2019) or data-constrained (Kliem et al. 2013; Amari et al. 2018; Inoue et al. 2018). In the former, the boundary condition is updated in real time, while in the latter, it is fixed to the initial value or assigned numerically. Previous studies showed that data-constrained MHD simulations could reproduce magnetic flux rope eruptions if the initial condition is constructed using observational data not much earlier than the eruption time. Therefore, if we lack sufficient observational data covering the whole eruption process to perform data-driven MHD simulations, we can instead perform data-constrained ones that could reproduce solar eruptions.

Coronal nonpotential magnetic field is usually constructed by NLFFF models (Sakurai 1981; Wiegelmann & Sakurai 2012; Gilchrist & Wheatland 2014; Guo et al. 2017). However, most numerical NLFFF models cannot construct magnetic flux rope structures when the bottom magnetic field is weak or the flux ropes lie high in the atmosphere. A good way to overcome this drawback is to manually insert a flux rope a priori with the help of an analytical model. The Titov–Démoulin model (Titov & Démoulin 1999), which assumes that an active region comprises a major magnetic flux rope with an electric current channel and a background potential field, has been used in many pioneering efforts (e.g., Török et al. 2004; Török & Kliem 2005; Mei et al. 2017). In the Titov–Démoulin model, the path of the electric current channel is assumed to be a toroidal arc. Such a restriction has been relaxed in a newly developed method, the regularized Biot–Savart laws (RBSL; Titov et al. 2018), which allows one to prescribe an arbitrary shape of the electric current channel. Guo et al. (2019b) implemented the RBSL model in the Message Passing Interface Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al. 2003, 2012, 2021; Porth et al. 2014; Xia et al. 2018) to study the magnetic flux rope structure of an intermediate prominence. MPI-AMRVAC is an open-source and fully documented ⁴ numerical simulation code including multiple physics models and multiple numerical algorithms, powered by massive parallelization and adaptive mesh refinement capability.

There are several advantages in using RBSL to construct the initial condition for an MHD simulation. RBSL could reconstruct a flux rope in weak magnetic field regions, such as quiescent and intermediate prominence regions. RBSL could also reconstruct flux ropes high up in the corona, which are common structures usually manifested or associated with prominences, filaments, cavities, hot channels, and sigmoids (Cheng et al. 2017). Moreover, flux ropes reconstructed by RBSL and the magnetofrictional (MF) method are approximately force-free, and can be used as suitable initial conditions close to an equilibrium state before the eruption. We also note that such reconstructed flux ropes are in a nonpotential state and possess large free energy—a situation that is close to the real cases where eruptions are prone to occur. The above advantages make the method of RBSL plus MF a good way to initiate a data-constrained MHD simulation.

Here, we make data-constrained MHD simulations for an eruption event on 2011 June 21 using initial conditions computed by the method of RBSL plus MF with varying mass density. Our primary goal is to find an optimized simulation that reproduces as many observational features as possible. Observations and modeling methods are described in Section 2. We display an optimized model and compare it with observations in Section 3. Finally, we summarize the results and discuss the simulations in Section 4.

The event we model is a GOES C7.7 class flare that occurred on 2011 June 21. The flare was initiated at 01:18 UT and peaked at 03:26 UT. An unusual characteristic of the flare is that the 1.0–8.0 Å soft X-ray emission lasted until 12:00 UT as studied by Zhou et al. (2017), suggesting that the flare is an extremely long-duration event. The flare also rose slowly, which provides a good example for us to study the initial driving mechanisms of the flux rope eruption.

Figure 1 shows the evolution of the flare in multiple EUV wavelengths observed by the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) aboard the Solar Dynamics Observatory (SDO; Pesnell et al. 2012). The emission at 94 Å of Fe xviii (∼6.3 MK) reveals an expanding magnetic flux rope in the onset and rising phases at 01:12 (Figure 1(a)) and 02:12 UT (Figure 1(d)), respectively. The images at this wavelength also show the flaring loops at 03:12 UT close to the peak phase (Figure 1(g)). At the He ii 304 Å wave band (∼0.05 MK), we can see a filament along the path of the flux rope, where the filament material flowed in the flux rope, as shown in Figures 1(b) and (e) and the animation attached to Figure 1(b). The emission at this lower temperature was dominated by flare ribbons after the flare peak time. We also plot composite images by combining 304 Å, 211 Å of Fe xv (∼2 MK), and 171 Å of Fe ix (∼0.6 MK) emissions that originate in a wide region spanning the upper chromosphere, the transition region, and the quiet and active-region corona, which show the evolution of the flare (Figures 1(c), (f), and 1(i)). These observations provide necessary constraints on the input parameters for the following simulations.

We select a zero-β MHD model to simulate the evolution of the magnetic flux rope. This model neglects the gas pressure and gravity in the momentum equation, and neglects the energy equation by assuming the temperature to be zero. The basic equations are as follows:

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}{\boldsymbol{v}}-{\boldsymbol{B}}{\boldsymbol{B}})+{\rm{\nabla }}\left(\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{2}\right)=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{v}}{\boldsymbol{B}}-{\boldsymbol{B}}{\boldsymbol{v}})=-{\rm{\nabla }}\times (\eta {\boldsymbol{J}}),\end{eqnarray} \tag{ 3 }$$

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

where ρ represents the density, v the velocity, B the magnetic field, η the resistivity, and J the electric current density. This set of equations is the same as that used in Toriumi et al. (2020), but different from the one used in Guo et al. (2019a). Here, we omit the artificial density diffusion term in Equation (1) and the viscosity term in Equation (2). The resistivity term in Equation (3) is kept since it ensures a correct flux rope topology as shown in the following sections. Equations (1)–(4) are normalized such that the vacuum permeability, μ₀, equals 1. The normalization factors for the physical variables are listed as follows: L₀ = 1.0 × 10⁹ cm, t₀ = 85.9 s, ρ₀ = 2.3 × 10⁻¹⁵ g cm⁻³, v₀ = 1.2 × 10⁷ cm s⁻¹, and B₀ = 2.0 G. We note that these numbers are rounded off for clarity here, while they are provided in double precision in the code, and the relationships between these quantities, such as L₀ = t₀ × v₀ and $\displaystyle {B}_{0}^{2}/8\pi =\tfrac{1}{2}{\rho }_{0}{v}_{0}^{2}$, are satisfied. The resulting time step is always larger than 3.0 × 10⁻⁴ in the normalized unit, which is a reasonable value in numerical simulations. We usually do not have the problem of negative density. But we set a lower density limit to keep the Courant–Friedrichs–Lewy (CFL) condition to a reasonable value. Otherwise, it would drop too low and make the simulation unaffordable.

The initial condition for the magnetic field is derived by combining the RBSL method (Titov et al. 2018) and the MF method (Guo et al. 2016a, 2016b). Here, we use the same event and the same procedure as used in Guo et al. (2019b). There are four important parameters in the RBSL method, namely, the axis path, the minor radius, the magnetic flux, and the electric current of the flux rope. The axis path is reconstructed by triangulation from observations by SDO/AIA and the Extreme Ultraviolet Imager (EUVI; Wuelser et al. 2004) aboard the twin Solar Terrestrial Relations Observatory (STEREO). The minor radius, a, of the flux rope is constrained by the width of the filament, and we adopt a = 0.02 R_⊙ = 14 Mm in this case. The typical magnetic flux is assumed to be the average of the unsigned values on the two footprints of the flux rope within the minor radius a, namely, F₀ = (∣F₊∣ + ∣F₋∣)/2 = 3.7 × 10²⁰ Mx in this event. Nevertheless, we also test different magnetic fluxes to check how they affect the eruption behavior of the flux rope in the following simulations. Finally, the typical electric current is calculated to be ${I}_{0}=(5\sqrt{2}{F}_{0})/(3{\mu }_{0}a)\,=5.0\times {10}^{11}$ A. Of course, the electric current also changes if we adopt different magnetic fluxes. Figure 2(a) shows a magnetic flux rope with F = 3F₀ and I = 3I₀, which represent optimized parameters to yield reasonable simulation results comparable to observations. The MHD simulation is performed in a local Cartesian coordinate system, where x-, y-, and z-axes point in the local west, north, and radial directions, respectively. The simulation domain is $[{x}_{\min },{x}_{\max }]\times [{y}_{\min }$, ${y}_{\max }]\times [{z}_{\min },{z}_{\max }]=[-26.5,26.5]\times [-26.5,26.5]\times [0.1,53.2]$, which is resolved by a stretched grid with 360 × 360 × 180 cells. The grid size is 0.15 per cell in the x and y directions, and ranges from 0.03 to 1.07 geometrically from bottom to top in the z direction. We note that the field of view in the current MHD simulations is larger than that for RBSL models in Guo et al. (2019b). The height extends to h = 53.2 (i.e., 5.32 × 10⁵ km), which is comparable to the solar radius of 6.95 × 10⁵ km. This advantage allows us to follow a flux rope to a sufficiently high altitude.

The initial condition for the density is provided by a plane-parallel atmosphere model, which is solved by dp/dz = −ρg and p = ρ T in the normalized form. We prescribe a temperature profile with the same form and parameters as used in Guo et al. (2019a). The temperature profile comprises two constant values of T_pho = 0.006 below h = 0.35 in the photosphere and T_cor = 1.0 above h = 1.0 in the corona. They are connected by a linear function between h = 0.35 and h = 1.0. The normalization unit of the temperature is T₀ = 1.0 × 10⁶ K. It is noted that the temperature profile is only used in the derivation of the initial density profile of the atmosphere model, while it does not explicitly appear in the MHD equations of the zero-β model (since gravity is ignored and temperature and pressure are zero by definition). The idea is to initialize the density in a realistically stratified way. The normalized density drops from 1.0 × 10⁸ at the bottom to 0.8 at the top of the simulation box. Since the zero-β model allows the use of any density profile as the background atmosphere, we can simulate the effect of mass inertia on the evolution of the flux rope. After a large number of tests, we find that the observed height–time profile of the filament, in particular its two phases of evolution, can only be reproduced by adopting two different density profiles. We thus split the MHD simulation into two stages, one with an initial density higher than usual in the corona and the other with an initial density given by the atmospheric model. The former case corresponds to the first evolution stage from 01:12 UT, where the initial density is 1.0 × 10⁴ ρ₀, about 100 times the background density at the height of the prominence (the average value is 1.0 × 10² ρ₀). Such a case represents the added mass density within the filament. The latter case corresponds to the second evolution stage from 02:10 UT, when the filament starts to accelerate. We describe both cases in detail in Section 3.2. Note that the initial velocity is zero everywhere in the computation box.

We adopt the same boundary conditions as those used in the data-constrained simulation in Guo et al. (2019a). Namely, we assign the initial density values to two ghost layers on all the six boundaries. The velocity on these boundary layers is set to be zero during the whole simulation. The observed vector magnetic field at 01:12 UT on 2011 June 21 is assigned to the inner ghost layer of the bottom boundary. For other ghost layers, the magnetic field is determined by a zero-gradient extrapolation, where the normal component of the vector magnetic field is reset to satisfy ∇ · B = 0. We note that all the physical quantities, including density, velocity, and magnetic field, are assigned at the centers of the ghost cells.

Equations (1)–(4) are solved by MPI-AMRVAC. Specifically, we use a three-step method for time progression, the HLL Riemann solver to derive the fluxes of the physical variables, and the Koren limiter to reconstruct variables from cell center to cell face. The resistivity in Equation (3) is chosen to be η = 0 when J < J_c,r, and $\eta ={\eta }_{0}{[(J-{J}_{c,{\rm{r}}})/{J}_{c,{\rm{r}}}]}^{2}$ when J ≥ J_c,r, where J_c,r is the critical current density to invoke the resistivity, and η₀ is the coefficient of the anomalous resistivity. In the following, we first present a case in Section 3 with J_c,r = 5.0 in dimensionless units (equivalent to 8.0 × 10⁻⁹ A cm⁻² in cgs units) and η₀ = 0.1 (equivalent to 1.2 × 10¹⁶ cm² s⁻¹ in cgs units). The simulation time ranges from 01:12 to 03:12 UT. Control parameters for this optimized Model 1 are summarized in Table 1. We also discuss some test cases by varying the control parameters in Section 3.2 and an ideal MHD case with η = 0 everywhere in Section 4.

3.1. Simulation Results and Comparison with Observations

Figure 2 shows the simulation results with the parameters described in Section 2. Three snapshots in each column correspond to the flare onset, rising, and peak times, which are exactly the same times for observed images as shown in Figure 1. Animations attached to Figure 2 display a more detailed evolution process with a time cadence of 2 minutes from 01:12 to 03:12 UT. The flux rope evolves slowly before 02:10 UT as shown in Figures 2(a)–(d), and the attached movies. It rises and expands drastically after 02:10 UT as seen in Figures 2(e) and (f). The electric current density shows a fluctuating behavior during the whole period: it is initially high in the flux rope (Figure 2(b)) and then dissipates to a low level in the slowly evolving process (Figure 2(d)); however, it increases to a high level again when the flux rope starts to erupt and finally decreases again to an even lower level (watch the attached animation after 02:14 UT).

To compare the simulation results with the observations, we back-project the magnetic field vector and its location to the viewing angles of STEREO_A/B and SDO, where the new X-, Y-, and Z-axes point to the west, north, and the observer, respectively. The back-projection is realized by three elementary rotations ${{ \mathcal R }}_{x}(-{B}_{0}){{ \mathcal R }}_{y}(-L){{ \mathcal R }}_{x}({B}_{1})$, where the rotation matrices can be found in Guo et al. (2017). The latitude of the disk center, B₀, the longitude of the active region reference point measured from the local central meridian, L, and the latitude of the reference point, B₁, for STEREO_A/B and SDO are (69, −90.9, 142), (− 72, 981, 142), and (17, 56, 142), respectively.

We overlay the back-projected flux ropes on the 304 Å images in Figure 3. The morphology of the flux rope resembles that of the filament very well (Figures 3(a), (b), and (c)). The eruption path deviates from the radial direction, leaning toward the south (Figures 3(d), (e), and (f)), which is consistent with the observations obtained by EUVI aboard STEREO_A/B. The apex of the flux rope reaches a large altitude at 03:12 UT (Figures 3(g) and (i)). And the flux rope keeps rising after 03:12 UT, since the velocity is still large as shown in Section 3.2. This result is consistent with the observations of the Large Angle and Spectrometric Coronagraph (LASCO; Brueckner et al. 1995) C2 and C3 on the Solar and Heliospheric Observatory (SOHO), which display a halo coronal mass ejection associated with the flux rope eruption. Here, we only simulate the flux rope eruption until 03:12 UT. After that, with the flux rope rising high and the filament material draining down, the flux rope becomes very weak in emission. Thus, the observations obtained by STEREO_A/B and SDO provide few constraints on the simulation parameters. In the future, we expect to study a flux rope eruption that appears in the field of view of SOHO/LASCO C2 and C3, and envisage a corresponding simulation with a much larger scale. In such a case, we can study the propagation behavior of the flux rope in more detail. Here, we only focus on the onset and initial rising phases.

Zoom In Zoom Out Reset image size

We then compute the quasi-separatrix layers (QSLs) based on the magnetic field in the simulation results, and compare them with the observations. QSLs are 3D thin layers where the field linkage has a huge gradient (Démoulin et al. 1996). They delineate the boundaries of special magnetic domains, such as flux ropes, spine-fan structures, hyperbolic flux tubes, and reconnecting current sheets. QSLs are quantified by the squashing factor, Q, and are usually defined as regions where Q ≫ 2 (Titov et al. 2002). We compute the value of Q on two slices at each snapshot with the open-source code written by Kai E. Yang ⁵ . This code implements the method proposed by Scott et al. (2017). Note that another similar method has been proposed by Tassev & Savcheva (2017).

Figure 4 displays the evolution of QSLs at three selected snapshots. The evolution of lgQ distributions and 304 Å images in Figures 4(a), (d), and (g) shows that the QSLs on the bottom boundary (photosphere) reproduce the flare ribbons very well, both in their shapes and evolution behavior (also see the animation attached to Figure 4). The most prominent feature shown in Figures 4(a) and (b) is the QSLs surrounding the flux rope at its footprints. Figure 4(c) shows the circular QSLs surrounding and in the body of the flux rope. The time for Figures 4(d)–(f) is different from that for the middle panels of Figures 1–3, since the QSLs coinciding with the separating flare ribbons have not fully developed at 02:12 UT. Instead, we show the QSLs at 02:20 UT in Figure 4(d) to highlight the separating QSLs within the orange polygon. A hyperbolic flux tube, where the QSLs intersect with themselves, develops under the flux rope as shown in Figure 4(f). The QSLs coinciding with the flare ribbons further separate from each other (Figures 4(g) and (h)), and the QSLs at the boundary and in the body of the flux rope expand drastically at 03:12 UT (Figure 4(i)). We also note that QSLs in Figures 4(f) and (i) display some tree-ring patterns, which represent the complexity of the flux rope and the precision of the MHD simulations in this study. Both the ability to reproduce the evolution of the flare ribbon and the ability to resolve the complex structure of the QSLs indicate the advantage of the present MHD simulations.

Download figure:

Video Standard image High-resolution image

3.2. Flux Rope Kinematics

The kinematics of the flux rope obtained from the observations provides a quantitative constraint on the simulations. Therefore, we first measure the evolution of the flux rope, which is manifested by the hot channel shown at SDO/AIA 94 Å and the filament at STEREO_B/EUVI 304 Å. Figures 5(a) and (b) display such observations at two selected times, when the hot channel and filament are clearly observed. We select two slices along the direction of the flux rope eruption as shown in Figures 5(a) and (b). We then align the time series of the slice to get the time–distance diagrams at 94 Å and 304 Å as shown in Figures 5(e) and (f), respectively. From them we visually determine the apexes of the hot channel and filament. The measurements are repeated ten times to reduce the errors. We then get the average positions as shown by the red and green diamonds and the standard deviations as shown by the error bars in Figures 5(e) and (f). Most importantly, the time–distance profiles reveal two phases of evolution, namely, a slow evolution phase and a fast eruption phase. Finally, we perform linear fittings to the fast eruption phase of the time–distance profiles. The velocities are 76.6 and 78.8 km s⁻¹ for the observations at 94 Å and 304 Å, respectively. Comparatively, the velocities in the slow evolution phase are small (below a few km s⁻¹).

Zoom In Zoom Out Reset image size

Flux rope eruptions comprising two phases of evolution are now thought to be a common phenomenon (Cheng et al. 2020). Therefore, any MHD model has to reproduce the two phases of kinematic evolution as shown in the time–distance profiles. As mentioned in Section 2, our simulations consider two different density distributions for the two phases: a high density in the background corona (also in the flux rope) for the slow evolution phase, and a low density for the fast eruption phase. As already mentioned in Section 2, we first determine the density distribution of the corona by solving the hydrostatic stratification. The average density in the corona at the height of the flux rope apex is about 1.0 × 10² ρ₀. Since the flux rope is filled with filament material, its density is much higher than the average density in the corona. We have to artificially set a higher initial density for the simulation of the slow evolution phase. Here, we change the density to 1.0 × 10⁴ ρ₀ at the places where the density is initially lower than this value. It is worth noting that a factor 100 in density is the typical value quoted for any prominence material. At 02:10 UT, however, the density is reset to the initial value calculated from the atmosphere model.

This setup is motivated by drainage of the filament material, which is universally observed before an eruption (Jenkins et al. 2019; Chen et al. 2020). Such drainage is present in this event during the filament evolution (see the animation attached to Figure 1). In the slow evolution phase, the magnetic flux rope lies lower in the atmosphere with higher density. When it starts to rise, filament material drains down and the mass density in the flux rope decreases. However, the zero-β model cannot deal with the density evolution accurately. That is why we adopt two distinct initial density distributions for the two evolution phases. Figure 6(a) shows a typical snapshot for Model 1 at 02:32 UT. It is used for comparison with other test cases discussed below.

Zoom In Zoom Out Reset image size

The eruption velocity of the flux rope is sensitively controlled by the magnetic field strength and the density within the flux rope. We have tested different cases with different magnetic fluxes of the flux rope, namely, F = 1F₀, 2F₀, 3F₀, and 5F₀, but the same initial density from the atmosphere model. The results show that the time–distance profiles are all linear in the four cases. The flux rope starts to rise at the beginning of the simulation, and the two phases of the eruption, as shown in the observations, cannot be reproduced. Quantitatively, the velocities in the cases of F = 1F₀ and 2F₀ are too small to be compared to the observations in the fast eruption phase, but the velocity in the case of F = 5F₀ is too large. Figure 6(b) shows a snapshot for the test case of F = 5F₀ at 01:52 UT, whose control parameters are listed in Table 1 as Model 2. It is found that the height of the flux rope of Model 2 at 01:52 UT is close to that of Model 1 at 02:32 UT, which implies that the flux rope of Model 2 erupts faster than that of Model 1.

We also fix the magnetic flux as F = 3F₀ and test different density models. One group of tests is to restrict the high-density region to a small area along the flux rope, roughly assumed to be the location of the filament, but to reset the density to a normal value at 02:10 UT. In practice, we use the electric current density as a good physical parameter to restrict the high-density region. Five cases are tested, with the critical electric current density J_c,d = 0.0, 0.01, 0.1, 5.0, and 20.0. The density is initially set as 1.0 × 10⁴ ρ₀ when the current density is larger than J_c,d and the background density smaller than 1.0 × 10⁴ ρ₀. When J_c,d = 0.0, the density in the whole corona is set to a high value; however, when J_c,d = 20.0, only the density in the flux rope is set to the high value. The simulation results show that the cases with J_c,d = 0.01, 0.1, 5.0, and 20.0 cannot suppress the velocity to a rather small value as shown in observations before 02:10 UT. Figures 6(c) and (d) show Model 3 with J_c,d = 5.0 (see Table 1 for other control parameters) at two snapshots, 01:12 UT and 01:38 UT. The density is higher than the background corona in the vicinity of the flux rope at 01:12 UT as shown in Figure 6(c), but the flux rope has risen higher than its initial height as shown in Figure 6(d). Note that it is supposed to be steady before the density profile is reset at 02:10 UT. With all the test parameters, we find that only the case with J_c,d = 0.0 could reproduce the observed time–distance profile. We also make a test by changing the time of resetting the density to normal from 02:10 UT to 01:50 UT. The results show that the fast eruption phase indeed occurs earlier correspondingly.

With the aforementioned tests, the observed time–distance profile can be reproduced in the simulation with a magnetic flux of F = 3F₀ and a higher initial density distribution before 02:10 UT but a normal one after it. We measure the apex of the flux rope axis as shown in Figures 5(c) and (d). The measurements are repeated ten times for each view angle, one from SDO and the other from STEREO_B. Then, the time–distance profiles are overlaid on that derived from observations (Figures 5(e) and (f)). The errors are estimated as the standard deviation of the measurements. Therefore, by adopting appropriate parameters, the simulation can indeed reproduce both the slow evolution and fast eruption phases of the flux rope kinematics. Note that the flux rope in the simulation rises to a higher altitude than the observations could trace. The simulation thus makes up the missing story of the flux rope when it moved outside the field of view in observations.

3.3. Eruption Mechanism

Magnetic flux rope eruptions are often thought to be driven by ideal MHD instabilities, or equivalently, the loss of an equilibrium mechanism (Démoulin & Priest 1988; Forbes & Isenberg 1991; Lin & Forbes 2000; Démoulin & Aulanier 2010; Kliem et al. 2014). MHD instabilities include helical kink instability (Török et al. 2004) and torus instability (Kliem & Török 2006) in one single current channel, tilt-kink instability between two parallel current channels (Keppens et al. 2014), and coalescence-kink instability between two antiparallel current channels (Makwana et al. 2018). Our simulation is most relevant to the kink and torus instabilities since there is only one major flux rope in this case. The kink instability is invoked when the twist angle is larger than a critical value, which depends on the specific equilibrium parameters, such as the magnetic field distribution, plasma β, and the geometric parameters of the flux rope. Some typical critical twist angles for the kink instability have been proposed to be 2.5π (Hood & Priest 1981), 3.3π (Hood & Priest 1979), and 3.5π (Török et al.2004). However, it can be as large as ∼5π (Mikic et al. 1990; Baty & Heyvaerts 1996) or 6π (Hood & Priest 1979) for different magnetic configurations.

The torus instability occurs when the decay index, $n\,=-d\mathrm{ln}{B}_{\mathrm{ex},\mathrm{pol}}/d\mathrm{ln}R$, exceeds a critical value n_c (e.g., Liu 2008), where B_ex,pol is the poloidal component of the external magnetic field and R is the major radius of the torus-shaped flux rope. This model relies on a force balance between the outward hoop force, generated by the toroidal current and the internal poloidal magnetic field, and the inward strapping force, generated by the toroidal current and the external poloidal magnetic field. Therefore, when the external poloidal magnetic field decreases fast enough, the hoop force dominates over the strapping force, and the torus instability occurs. A classical instability analysis shows that the critical decay index is n_c = 1.5 (Bateman 1978; Kliem & Török 2006) when the current channel is semicircular and without any internal structures. It has also been suggested that n_c = 1 if the current channel is thin and straight (van Tend & Kuperus 1978). Démoulin & Aulanier (2010) found that n_c ∈ [1.1, 1.3] when the current channel has a deformable shape and a finite thickness. Olmedo & Zhang (2010) concluded that n_c depends on the ratio of the height of the flux rope apex, Z, to the half-distance of the flux rope footprints, S₀. When Z/S₀ ≥ 1 (Z/S₀ = 1 is equivalent to a half-torus), they found that n_c ∈ [0, 2] depending on the distribution of the external magnetic field. Fan (2010) studied the conditions for the dynamic eruption of a flux rope by means of an isothermal MHD simulation. The torus instability is also regarded as the driving mechanism of the eruption and the critical decay index is found to be n_c = 1.74 in their case.

Magnetic flux rope eruptions can be either successful or confined to the lower corona. A well-known mechanism for the confined eruption is the failed kink scenario (Török & Kliem 2005), when the kink instability occurs while the torus instability does not. In this case, a flux rope undergoes a writhing motion by converting part of its twist to writhe and the eruption stops at the lower corona. Myers et al. (2015) found another mechanism to confine an eruption, namely, the failed torus scenario, which occurs when the magnetic tension force generated by the toroidal magnetic field and poloidal electric current is large enough to counteract the outward hoop force. Large magnetic tension force requires a low twist number of a flux rope. Therefore, confinement of an eruption requires a condition that either the twist or the decay index is low, while a successful eruption is generated when both the twist and decay index are large enough. This conclusion is clearly presented in the phase diagram of twist (reciprocal of the safety factor) versus decay index as proposed by Myers et al. (2015).

We provide the phase diagram of twist versus decay index in the evolution process of the magnetic flux rope in Figure 7. The twist is computed by using the Berger–Prior formula (Berger & Prior 2006). Two key parameters in the calculation of the twist are the axis curve and sample field lines of the flux rope. The axis at each snapshot of the MHD models is determined as the field line most perpendicular to the vertical plane as shown in Figure 4(b). We select 100 sample field lines randomly distributed within the border of the closed QSLs. The decay index is computed along an eruption path as shown in Figures 5(c) and (d), where the two slices determine the 3D eruption path with a unit vector of e _pro. Then, we assume the flux rope axis points to e _axi, which is assumed to be a unit vector perpendicular to e _pro and parallel to the horizontal plane. There are two solutions for e _axi. We choose the one yielding a positive poloidal field component. Finally, we determine the poloidal direction with a unit vector, e _pol, by requiring e _pol × ( e _pro × e _axi) = 0. The poloidal component of the external magnetic field, B_ex,pol, is determined as the scalar product of the potential magnetic field, B _pot, and the poloidal unit vector, e _pol, namely, B_ex,pol(R) = B _pot(R) · e _pol. The distance R is measured along the eruption path starting from the bottom boundary, which is also assumed to be the major radius of the erupting flux rope. The decay index is then computed as $n=-d\mathrm{ln}{B}_{\mathrm{ex},\mathrm{pol}}/d\mathrm{ln}R$. Many previous studies only used a stationary potential field model to study the decay index (e.g., Guo et al. 2010; Xu et al. 2012; Zuccarello et al. 2014; Li et al. 2019). Our study has both similarities and differences compared to these previous studies. We all adopt the potential field as the external magnetic field, even in the MHD simulation, since the normal magnetic field on the bottom boundary does not change too much in the eruptive process. However, in the present study, we use the evolution of the flux rope in the MHD simulation to determine the position and poloidal direction, which are important parameters to compute the decay index.

Zoom In Zoom Out Reset image size

Figure 7 shows that the twist varies in a small range of 2.69 to 2.75 and the decay index changes from 0.84 to 0.64 from 01:12 to 02:10 UT. Then, the twist remains almost constant from 02:10 to 02:28 UT, while the decay index increases drastically from 0.64 to 2.37 during this period. After 02:28 UT, both the twist and decay index show a decreasing trend, except for a slight increase in the decay index in the very late phase after 02:50 UT. Figure 2 and the attached movies show that the flux rope does not rotate and writhe obviously during the eruption process. This indicates that the kink instability might not play a significant role in this case. Instead, the torus instability is likely the major mechanism to drive the eruption. Since the twist is high, the magnetic tension force would not be large enough to confine the eruption, which is consistent with the successful eruption of this event. Besides, as implied by our simulations in Section 3.2, material draining due to the rising motion of the flux rope plays an important role in promoting the eruption, which reduces the inertia of the flux rope and allows it to more easily enter the torus instability regime.

We performed a data-constrained MHD simulation for the long-duration flare of C7.7 class on 2011 June 21. The zero-β MHD model is adopted in this study. The initial condition is provided by a magnetic flux rope structure constructed with the RBSL and MF methods. The bottom boundary is prescribed with the vector magnetic field observed at the initial time of the eruption and with the velocity fixed to zero. The simulation results reproduce the observations in several aspects. The simulated magnetic flux rope propagates in a direction that deviates from the radial one and leans toward the south, which is consistent with what is revealed by the SDO/AIA 94 and 304 Å observations. The simulated flux rope continues to rise when it moves close to the computation boundary, which indicates a successful eruption as seen from the SOHO/LASCO observations. In particular, we reproduce both the slow evolution and fast eruption phases of the flux rope kinematics as observed at 94 Å and 304 Å by assuming two different density distributions in two stages: a higher initial density in the first stage and a lower (normal) initial density in the second. Such a change in the initial density roughly imitates the draining process of filament material, which cannot be fully incorporated in zero-β simulations.

The zero-β MHD model omits the gas pressure, gravity, and the energy equation. This is acceptable if the magnetic pressure is much larger than the gas pressure and gravity, which is plausible in the solar corona. There are also some difficulties in the full thermodynamic MHD simulations, which require for very high spatial resolutions and small time steps to resolve the pressure gradient and energy evolution, making them too expensive to be affordable in simulating flux rope dynamics from the photosphere to the corona. In this sense, the zero-β MHD model is economic and able to simulate the magnetic topology, flux rope velocity, and effects of resistivity reasonably in many circumstances. On the other hand, the thermodynamics of the plasma cannot be properly revealed in such a simplified model. Therefore, a zero-β MHD model is unable to simulate the filament formation and dynamics in a consistent way. In our case, the mass density within the filament is supposed to be very high at the initial time in order to stabilize the magnetic flux rope for a relatively long period, an hour or so before the eruption. After the rising motion of the flux rope is triggered, the density gradually decreases due to material draining. In order to take into account such a change in the density, we adopt two different initial density values in the two stages of our simulations. Although this is a simplified treatment, our results do indicate that the material draining can play an important role in promoting the flux rope eruptions. This further confirms the finding by Zhou et al. (2006) and Fan (2020). In the future, we need a full MHD model to more accurately treat the thermodynamics of the flux rope that is omitted by the zero-β model.

In our simulations, the resistive term in Equation (3) is determined such that η = 0 when J < J_c,r and $\eta \,={\eta }_{0}{[(J-{J}_{c,{\rm{r}}})/{J}_{c,{\rm{r}}}]}^{2}$ when J ≥ J_c,r, where J_c,r = 5.0 and η₀ = 0.1. In an ideal case, the resistive term is zero everywhere, as adopted in many simulations (e.g., Török & Kliem 2003; Fan & Gibson 2007). We also tested the ideal case by adopting η = 0 everywhere. Figures 6(a) and (e) compare the simulation results of the resistive and ideal cases (Models 1 and 4 in Table 1) at two typical times, when the eruption fully develops. It is found that the middle part of the flux rope kinks down to the bottom boundary in the ideal case (Figure 6(e)), while the flux rope erupts entirely in the resistive case (Figure 6(a)). Both cases produce current sheets under the erupting flux rope, where stretched envelope field lines of opposite directions meet and reconnect (also refer to Yokoyama & Shibata 1997). Note that in the ideal case, unavoidable numerical resistivity controls the reconnection processes, driven by round-off and truncation errors. We have also tested a few cases by varying the parameters J_c,r and η₀, namely, J_c,r = 0.1, 5.0, and 20.0, and η₀ = 0.001, 0.01, 0.1, and 1.0. The results show that a large value of J_c,r (≥20.0) or a small value of η₀ (≤0.01) cannot eliminate the kink in the middle part of the flux rope. In contrast, small J_c,r (≤0.1) or large η₀ (≥1.0) would dissipate the flux rope significantly. Figure 6(f) shows the result of Model 5 with η₀ = 1.0 at 01:14 UT. The evolution of the flux rope is fully controlled by the resistive time step, which is much smaller than the dynamic time step. And the flux rope is dissipated rapidly. We therefore conclude that this specific event requires an anomalous resistivity prescription as adopted, because it would otherwise show different topological evolutions. In most simulations, η is determined empirically and changes from case to case.

To better constrain future studies, we need more observations to construct the initial magnetic field model, including the axis path, minor radius, magnetic flux, and electric current of the flux rope. We also need a complete parameter survey in order to construct an equilibrium model including the thermal conditions. Finally, full MHD simulations are required to include more physics such as gravity, gas pressure, thermal conduction, radiation loss, and background heating, which are missing from the zero-β models.

The SDO data are provided by NASA/SDO and the AIA and HMI science teams. The SECCHI data are provided by STEREO and the SECCHI consortium. Y.G., Z.Z., M.D.D., and P.F.C. were supported by NSFC (11773016, 11733003, 11961131002, and 11533005) and 2020YFC2201201. C.X. was supported by NSFC (11803031). R.K. was supported by a joint FWO-NSFC grant G0E9619N and received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 833251 PROMINENT ERC-ADG 2018) and by Internal funds KU Leuven, project C14/19/089 TRACESpace. The numerical computation was conducted in the High Performance Computing Center (HPCC) in Nanjing University.