Unveiling the Initiation Route of Coronal Mass Ejections through Their Slow Rise Phase

In this work, we study the CME initiation by performing a 3D observationally inspired thermal MHD simulation labeled as "Simulation Driven-eruption" (abbreviated as "Simulation De") with the code MPI-AMRVAC (Xia et al. 2018), which solves the following equations in Cartesian coordinates:

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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \,(\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot \left[\rho {\boldsymbol{v}}{\boldsymbol{v}}+\left(\,p+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{2{\mu }_{0}}\right)\,{\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{B}}{\boldsymbol{B}}}{{\mu }_{0}}\right]\\ \,=\,\rho {\boldsymbol{g}}+2\mu {\rm{\nabla }}\cdot \left[{\boldsymbol{S}}-\displaystyle \frac{1}{3}\,({\rm{\nabla }}\cdot {\boldsymbol{v}})\,{\boldsymbol{I}}\right]\end{array}\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {e}_{\mathrm{int}}}{\partial t}+{\rm{\nabla }}\cdot ({e}_{\mathrm{int}}\,{\boldsymbol{v}})=-p{\rm{\nabla }}\cdot {\boldsymbol{v}}\\ +\,2\mu \,\left[{\boldsymbol{S}}:{\boldsymbol{S}}-\displaystyle \frac{1}{3}{\left(\,,{\rm{\nabla }}\cdot {\boldsymbol{v}}\right)}^{2}\right]+\eta {J}^{2}+{\rm{\nabla }}\cdot [{\kappa }_{| | }\,({\boldsymbol{b}}\cdot {\rm{\nabla }}T)\,{\boldsymbol{b}}]\end{array}\end{eqnarray} \tag{ 4 }$$

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

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

where ρ, p, T, e_int, v , B , J represent the mass density, thermal pressure, temperature, internal energy, velocity, magnetic field and current density, respectively. g = −g e _z is the gravity acceleration. μ is the dynamic viscosity coefficient and is set to 10⁻⁴ (in dimensionless unit) throughout the simulation, S is the strain tensor where S_ij = (1/2)(∂_i v_j + ∂_j v_i ), and I is the unit tensor. η is the resistivity coefficient whose value is set differently in different phases. κ_∣∣ = 8 × 10⁻⁷ T^5/2 erg cm⁻¹ s⁻¹ K⁻¹ is the parallel conductivity coefficient, and b = B /B is the normalized magnetic field. We introduce a generalized Lagrange multiplier ψ to maintain the ∇ · B = 0 condition (i.e., generalized Lagrange multiplier (GLM) method; Dedner et al. 2002). The evolution of ψ follows Equation (6), and c_h and c_p are constants. As a thermal MHD simulation, we consider compression heating, viscous dissipation, ohmic dissipation, and thermal conduction in the energy equation. In particular, we solve the internal energy equation instead of the total energy equation to avoid the heating from the numerical resistivity and numerical viscosity.

The equations are solved with the HLL scheme (Harten et al. 1983), the third-step Runge–Kutta time discretization method, and the fifth-order mp5-limited reconstruction (Suresh & Huynh 1997). We also use the magnetic field splitting method (Tanaka 1994; Xia et al. 2018) in this simulation: the magnetic field B is split into two parts, the invariable background field B ₀ and the deviation B ₁ . In this simulation, B ₀ is set as the initial magnetic field, and thus the initial B ₁ is 0. This method is conducive to eliminating numerical errors related to the divergence of the magnetic field and achieving more precise results, especially in regions where the magnetic field evolves with minor change from the initial value, although it is relatively less effective in regions where the change of magnetic field is considerable (e.g., the centers of polarities during the converging phase).

The equations are solved in the dimensionless form, with the magnetic permeability equal to 1. The atmosphere is set to a fully ionized ideal gas with a hydrogen-helium abundance ratio of 10:1 (Xia et al. 2012). The dimensionless unit of the length, time, mass density, thermal pressure, temperature, velocity, and magnetic field strength is 10 Mm, 67.89 s, 2.34 × 10⁻¹⁵ g cm⁻³, 0.51 erg cm⁻³, 1.6 MK, 147.30 km s⁻¹, and 2.53 G, respectively.

The physical domain of the simulation is in the range of −7 ≤ x ≤ 7, −7 ≤ y ≤ 7, and 0 ≤ z ≤ 14 (in dimensionless unit). We set a stretched grid (nx × ny × nz = 144 × 144 × 96), which is symmetrically stretched in x- and y-directions with a stretched ratio of 1.029, and unidirectionally stretched in z-direction with a stretched ratio of 1.028. The spatial resolution in three directions is in the range of 0.0297 ≤ dx ≤ 0.2262, 0.0297 ≤ dy ≤ 0.2262, and 0.0298 ≤dz ≤ 0.4103 (in dimensionless unit). In dimensional units, the finest resolutions in three directions are about 300 km, comparable to those of Solar Dynamics Observatory/Atmospheric Imaging Assembly (AIA; Lemen et al. 2012). This spatial resolution is lower than those in other simulations of CME eruptions (Aulanier et al. 2010; Jiang et al. 2021), but it is a compromise for the sake of saving computation resources for solving the thermal MHD equations.

The initial magnetic field (similar to that in Observationally driven High-order Magnetohydrodynamics code; OHM; Aulanier et al. 2010) is composed of two potential bipolar fields:

$$\begin{eqnarray}\begin{array}{rcl}{B}_{x}\,(t=0) & = & {{\rm{\Sigma }}}_{m=1}^{4}\,{c}_{m}\,(x-{x}_{m})\,{r}_{m}^{-3}\\ {B}_{y}\,(t=0) & = & {{\rm{\Sigma }}}_{m=1}^{4}\,{c}_{m}\,(\,y-{y}_{m}\,)\,{r}_{m}^{-3}\\ {B}_{z}\,(t=0) & = & {{\rm{\Sigma }}}_{m=1}^{4}\,{c}_{m}\,(z-{z}_{m})\,{r}_{m}^{-3}\\ {r}_{m} & = & \sqrt{{\left(x-{x}_{m}\,\right)}^{2}+{\left(\,y-{y}_{m}\,\right)}^{2}+{\left(z-{z}_{m}\,\right)}^{2}},\end{array}\end{eqnarray} \tag{ 7 }$$

where (c₁ = 60, x₁ = 0.9, y₁ = 0.3, z₁ = −1.1), (c₂ = −60, x₂ = −0.9, y₂ = −0.3, z₂ = −1.1), (c₃ = 45, x₃ = 9, y₃ = 3, z₃ = −11), and (c₄ = −45, x₄ = −9, y₄ = −3, z₄ = −11) in dimensionless unit. The first (last) two monopoles are close to (far from) each other and located shallow (deep), producing a strong (weak) potential field that decays rapidly (slowly) with altitude within the physical domain. The first potential field mimics the background field in an active region, while the second one mimics the background field of the Sun on a global scale. The second potential field is set to ensure that the Alfvén speed at the upper part of the domain is large enough and thus to avoid shocks caused by the driving motion. The maximum dimensionless field strength in the plane z = 0 (which is named as box bottom surface in the following) is about 45. Such a strength is smaller than that in active regions in observations, while this is also a compromise for saving computation resources.

The initial atmosphere is set to a hydrostatic corona with a uniform dimensionless temperature of 1. The box bottom surface corresponds to 3 Mm above the solar surface, and the initial dimensionless mass density at the box bottom surface is set to 1. The initial distribution of thermal pressure and mass density with altitude follows:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\nabla }}p & = & \rho g\\ p & = & \rho T\\ g & = & {g}_{0}\,\displaystyle \frac{{R}_{}^{2}}{{\left({R}_{}+0.3+z\right)}^{2}},\end{array}\end{eqnarray} \tag{ 8 }$$

where R_Sun = 69.55 is the dimensionless solar radius and g₀ = −0.126 is the dimensionless gravity acceleration at the solar surface. With such an initial atmosphere, the plasma β along the vertical slit through the origin is less than 1 in the region z < 3.7 (in dimensionless unit) at t = 0. While, as the magnetic field evolves, the plasma β further decreases: the plasma β along the vertical slit through the origin is less than 1 in the region z < 4.8 at t = 18, in the region z < 6.0 at t = 58, and in the region z < 7.5 at t = 68, which ensures that the modeled CME (progenitor) is always located in a low plasma β region which mimics the corona.

The initial velocity is set to zero, and the parameter ψ is also set to zero in the whole physical domain.

In this work, we impose line-tied motions on the bottom boundary to drive the magnetic field and set special bottom boundary conditions to ensure the line-tied condition. Here, the line-tied condition means that the footpoint of the field line can only move horizontally according to the prescribed motion in the case of ideal MHD (Aulanier et al. 2005).

Inspired by the characteristics of the photospheric magnetic field evolution in observations (Yang et al. 2004; Green & Kliem 2009; Schrijver et al. 2011), two driving motions are considered in this simulation: a shearing motion that drives the initial potential field to a highly sheared state, and a converging motion that facilitates the flux cancellation near the polarity inversion line (PIL). The "Simulation De" is composed of two phases: the shearing phase (0 ≤ t ≤ 18) with only shearing motion applied and the converging phase (18 < t ≤ 68) with only converging motion applied. Such a setup differs from that in the previous work (Zuccarello et al. 2015), as here the converging motion is never switched off. We note that imposing only one type of motion at each phase makes the simulation easier to control and analyze. The shearing phase can be understood as a preparatory stage of the simulation to create a sheared magnetic field for the converging phase, during the latter of which we analyze the formation and eruption of the flux rope. The converging phase is designed to fulfill the flux cancellation model.

The motions are imposed at the cell-center bottom surface, which is a horizontal surface at the altitude of the cell center of the first layer of the physical domain (cell layers in the physical domain indexed by k = 1, 2, 3, ..., 96 from the bottom to the top) at each time step to directly control the evolution of the magnetic field in the layer k = 1. For "Simulation De," we set:

the shearing motion (${v}_{x}^{s},{v}_{y}^{s},{v}_{z}^{s}$):

$$\begin{eqnarray}\begin{array}{rcl}{v}_{x}^{s}\,(k=1;t) & = & \gamma \,(t)\,{V}_{x}^{s}\,(t)\\ {v}_{y}^{s}\,(k=1;t) & = & \gamma \,(t)\,{V}_{y}^{s}\,(t)\\ {v}_{z}^{s}\,(k=1;t) & = & 0\\ {V}_{x}^{s}\,(t) & = & {v}_{0}^{\max \;}{{\rm{\Psi }}}_{0}\,(t)\,{\partial }_{y}\,{\rm{\Psi }}\,(t)\\ {V}_{y}^{s}\,(t) & = & -{v}_{0}^{\max \;}{{\rm{\Psi }}}_{0}\,(t)\,{\partial }_{x}\,{\rm{\Psi }}\,(t)\\ \gamma \,(t) & = & \left\{\begin{array}{ll}\displaystyle \frac{1}{2}\;\tanh \,[3.75\,(t-1)]+\displaystyle \frac{1}{2} & 0\leqslant t\lt 16\\ -\displaystyle \frac{1}{2}\;\tanh \,[3.75\,(t-17)]+\displaystyle \frac{1}{2} & 16\leqslant t\leqslant 18,\end{array}\right.\end{array}\end{eqnarray} \tag{ 9 }$$

and the converging motion (${v}_{x}^{c},{v}_{y}^{c},{v}_{z}^{c}$):

$$\begin{eqnarray}\begin{array}{rcl}{v}_{x}^{c}\,(k=1;t) & = & \gamma \,(t)\,{V}_{x}^{c}\,(t)\\ {v}_{y}^{c}\,(k=1;t) & = & \gamma \,(t)\,{V}_{y}^{c}\,(t)\\ {v}_{z}^{c}\,(k=1;t) & = & 0\\ {V}_{x}^{c}\,(t) & = & {v}_{0}^{\max \;}{{\rm{\Psi }}}_{0}\,(t)\,{\partial }_{x}\,{\rm{\Psi }}\,(t)\\ {V}_{y}^{c}\,(t) & = & {v}_{0}^{\max \;}{{\rm{\Psi }}}_{0}\,(t)\,{\partial }_{y}\,{\rm{\Psi }}\,(t)\\ \gamma \,(t) & = & \displaystyle \frac{1}{2}\;\tanh \,[3.75\,(t-19)]+\displaystyle \frac{1}{2}\,18\lt t\leqslant 68,\end{array}\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\rm{\Psi }}\,(t)=\exp \;\left[-{{\rm{\Psi }}}_{1}{\,\left(\displaystyle \frac{{B}_{z}\,(k=1;t)}{{B}_{z}^{\max }\,(k=1;t)}\right)}^{2}\;\right],\end{eqnarray} \tag{ 11 }$$

and Ψ₁ = 5.5 for the shearing motion and Ψ₁ = 27.5 for the converging motion. The maximum speed of the driving motion v₀ ^max is set to 0.16 (in dimensionless unit) for both the shearing and the converging motions. The function Ψ₀(t) is used to keep the maximum of the term $\sqrt{{V}_{x}^{2}+{V}_{y}^{2}}$ always v₀ ^max. As shown in Figure 1, the shearing motion is along the tangent direction of the contour of B_z (k = 1), while the converging motion is perpendicular to that tangent direction. Near the PIL, the shearing motion is antiparallel on two sides of the PIL, while the converging motion is toward the PIL.

Zoom In Zoom Out Reset image size

It should be noted that the maximum speed of the driving motions is set to roughly 20 times that in observations to speed up the simulation. While, it is still acceptable as it is less than the Alfvén speed and the sound speed in the core area affected by the driving motion: at t = 0, in the domain that −4.20 ≤ x ≤ 4.20, −4.20 ≤ y ≤ 4.20, and 0 ≤ z ≤ 8.36 (in dimensionless unit), the ratio of v₀ ^max to the average Alfvén speed is 0.101 and the ratio of v₀ ^max to the maximum Alfvén speed is 0.004; the ratio of v₀ ^max to the sound speed is constant of 0.124 in this domain.

Several settings are adopted to achieve the line-tied condition at the cell-center bottom surface. First, the GLM parameter ψ is fixed to zero in the layer k = 1, and we set a symmetric boundary condition centered at the cell-center bottom surface for ψ (see Section 2.5 for more details), so all components of the term ∇ · (ψ I ) are zero in the layer k = 1. Second, for "Simulation De," in the layer k = 1, the dissipation term in Equation (3) is modified into:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{v}}{\boldsymbol{B}}-{\boldsymbol{B}}{\boldsymbol{v}}+\psi {\boldsymbol{I}})\\ =\,\left\{\begin{array}{rl}0 & 0\leqslant t\leqslant 18\\ \eta \,\left(\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {y}^{2}}\right)\,{\boldsymbol{B}} & 18\lt t\leqslant 68.\end{array}\right.\end{array}\end{eqnarray} \tag{ 12 }$$

During the shearing phase of "Simulation De" (0 ≤ t ≤ 18), the dissipation term in this layer equals zero so that it has no effect on the magnetic field. The resistivity η is fixed to zero in the layer k = 1 and set to be uniform (η = 10⁻⁴; in dimensionless unit) in the whole region (including both the physical domain and the ghost cell layers) except the layer k = 1 during the shearing phase. However, during the converging phase of "Simulation De" (18 < t ≤ 68), we adopt a 2D dissipation term (Aulanier et al. 2010) in the layer k = 1 to allow the flux cancellation close to the PIL. A larger, uniform resistivity (η = 4 × 10⁻⁴; in dimensionless unit), which facilitates the flux cancellation, is set in the whole region, including the layer k = 1.

In addition, the shearing motion in Equation (9) satisfies that the z-component of the term ∇ · ( v B − B v ) is zero, so analytically B_z in the cell-center bottom surface should remain invariable during the shearing phase. To achieve this feature more accurately, we enforce that B_z remains invariable in the layer k = 1 during this phase.

The mp5-limited reconstruction method uses three ghost cell layers at each boundary. The mesh in the ghost cell layer is also stretched in the same way as that in the physical domain. In all ghost cell layers, the background magnetic field B ₀ remains invariable; the two components of the magnetic field B ₁ which are parallel to the boundary surface are derived by one-sided third-order equal-gradient extrapolation, and the component of B ₁ which is normal to the boundary surface is derived under the constraint of the ∇ · B = 0 condition. Such a setting will keep the divergence of the magnetic field zero in the layer k = 1, so it is reasonable to set ψ to zero in this layer.

We assume that at each time step, the atmosphere in both the bottom ghost cell layer and the layer k = 1 is hydrostatic. Thus, the thermal pressure in the bottom ghost cell layer can be derived with a centered difference. We further assume that the dimensionless temperature in the ghost cell layer is constant of 1, the same as the initial temperature, and then the mass density is derived by the ideal gas law. The thermal pressure and the mass density in the top ghost cell layers are derived with the same method. In the four side ghost cell layers, the thermal pressure and the mass density are fixed to their initial values.

The velocity in the bottom ghost cell layer k = n (layer indexed by k = −1, −2, −3 from the top to the bottom of the bottom ghost cells) is derived from the velocity in the layer k = 1 and that in the layer k = 1 − n by the linear extrapolation. This boundary condition is the same as the antisymmetric boundary condition for v_z and the pseudo antisymmetric boundary condition for v_x and v_y at the bottom boundary in the previous work (Aulanier et al. 2005). The velocity in top ghost cell layers is derived by the one-sided second-order zero-gradient extrapolation, and the z-component is forced to be no less than zero to avoid possible downward flows. The velocity in the side ghost cell layers is set by the antisymmetric boundary condition centered on the side surface of the box.

To achieve the goal of removing the contribution of ∂_z ψ in the layer k = 1, a symmetric bottom boundary condition for ψ is set centered on the cell center of this layer. The aim is fulfilled more precisely by replicating the flux of ψ from the top surface of the layer k = 1 to its bottom surface (box bottom surface). In the top and the side ghost cell layers, ψ is fixed to zero.

The evolution of the modeled CME progenitor and CME in "Simulation De" is shown in Figure 2. During the shearing phase (0 ≤ t ≤ 18), the initial potential field is driven to a highly sheared state under the shearing motion (Figure 2(b)). After that, during the converging phase (18 < t ≤ 68), a pre-eruptive flux rope composed of twisted field lines is first formed by magnetic reconnection under the converging motion, and then it erupts as a CME (Figure 2(b)). The pre-eruptive flux rope is hot and visible in synthetic EUV images mimicking AIA 335 Å observations (Figures 2(a) and 10).

The initial magnetic topology of the reconnection region is a bald patch (BP; Titov et al. 1993; during 20 ≤ t ≤ 45), characterized by ( B · ∇)B_z > 0 where B_z = 0. As shown by the quasi-separatrix layers (QSLs; Priest & Démoulin 1995), such a BP topology later gradually bifurcates into a hyperbolic flux tube (HFT; Titov et al. 2002; since t = 46; see Figures 5 and 8(f)). An obvious difference between the BP reconnection and the HFT reconnection is that the former forms flux rope field lines tangent to the bottom surface at their dips (Figure 3(a)), while the latter forms flux rope field lines whose dipped parts are detached from the bottom surface and low-lying loops (Figure 3(b)).

Zoom In Zoom Out Reset image size

The transition from the quasi-static phase to the slow rise phase occurs (at t = 52) just after the first appearance of HFT (at t = 46). This indicates that the HFT reconnection is likely the mechanism triggering the slow rise phase.

This speculation is supported by detailed analyses of the relation between the reconnection and the kinematics. Figures 5(a), (f) reveal that the upward outflow of the reconnection is very slow at the stage of BP (e.g., t = 45) and the early stage of HFT (46 ≤ t < 52), but it is markedly accelerated as the HFT is rapidly built up (e.g., from t = 52 to t = 58). Quantitatively, the growth rate of the reconnection upward outflow velocity is temporally consistent with the acceleration of the CME progenitor, and both of them start to increase shortly after the first appearance of HFT (Figure 5(g)). We consider that the velocity of the reconnection upward outflow reflects the ability of reconnection to contribute to the flux rope. Therefore, the above results indicate that the slow rise is triggered when the contribution from reconnection starts to quickly increase after the HFT forms. In other words, the slow rise is triggered by the HFT reconnection. We eliminate the possibility that the torus instability triggers the slow rise, since the failed eruption in "Simulation Ne" implies that the torus instability, as the triggering mechanism of the eruption, does not occur before t = 58 in the "Simulation De" (see Appendix A).

Zoom In Zoom Out Reset image size

We further investigate the mechanisms driving the slow rise of the flux rope. Considering that the torus instability occurs during the slow rise phase, we divide the slow rise phase into two stages, i.e., the earlier stage (52 ≤ t ≤ 58) and the later stage (59 ≤ t ≤ 65), and analyze the driving mechanisms in these two stages separately.

For the earlier stage, the motion of the flux rope is controlled by a net upward force on its lower part and a net downward force on its upper part (see Figure 6(b), taking an example at t = 58). The Lorentz force and the thermal pressure gradient force are the major contributors to the net force, while the gravity and the viscous force are of little importance (Figure 6). Especially, the maxima of the Lorentz force, the net force, and the velocity in the flux rope region are all concentrated near the HFT (Figures 6(a), (b)) during this stage. This implies that it is the HFT reconnection that powers the flux rope rise in the earlier stage of the slow rise phase. In addition, the synchronized evolutions of the flux rope acceleration and the growth rate of the reconnection upward outflow velocity in Figure 5(g) also indicate that the slow rise is driven by the HFT reconnection during the earlier stage.

Zoom In Zoom Out Reset image size

The HFT reconnection drives the slow rise in the earlier stage, on the one hand, by adding twisted concave-upward field lines to the lower part of the flux rope (Figure B1). These field lines provide magnetic tension to the flux rope (Figure 6(a); known as the slingshot effect), enhance the force imbalance in the flux rope, and finally accelerate the flux rope. On the other hand, the reconnection upward outflow in the HFT is convected upwardly into the flux rope (Figures 5(a), (e)), and thus it can directly promote the flux rope rise. Specifically, the plasma in the reconnection upward outflow is accelerated by a net upward force which is dominated by the Lorentz force and the thermal pressure gradient force (Figure 5). The upward Lorentz force consists of magnetic tension force mainly at the HFT top and magnetic pressure gradient force mainly at the HFT center (Figure 7), the former of which is provided by the newly formed concave-upward field lines and the latter of which is associated to the accumulation of the sheared/guide component of the magnetic field. The upward thermal pressure gradient force is due to an increase in the gas pressure at the flanks of the flux rope bottom, which is a result of the density pileup pinched by the reconnection inflow (Figure 7).

Zoom In Zoom Out Reset image size

Furthermore, since the HFT appears, the two footprints of the HFT continuously separate from each other in the direction perpendicular to the PIL, and meanwhile they extend in the direction parallel to the PIL to both sides until they are completely separated (Figures 8(a)–(c); shown by the QSL footprint, which manifests as the flare ribbon in observations during the eruption (Zhao et al. 2016), but probably invisible before the eruption due to the weak reconnection rate). The reconnection electric field within the HFT, which represents the reconnection rate (Forbes & Priest 1984), also shows a gradually increasing trend since t = 46 and until t = 59 (Figure 8(h); also see Appendix C for more details). These results strongly suggest that the HFT and its associated reconnection are gradually built up during the earlier stage of the slow rise phase. Both the slingshot effect and the outflow effect of the HFT reconnection thus become stronger, which leads to the continuous increase of acceleration of the CME progenitor in this stage.

Zoom In Zoom Out Reset image size

For the later stage of the slow rise phase, we argue that the HFT reconnection still has a considerable contribution to the rising motion, since the HFT reconnection becomes even stronger in this stage, shown by the quickly increasing reconnection electric field (see that during 59 ≤ t ≤ 65 in Figure 8(h)). Moreover, the torus instability also promotes the slow rise in the later stage since it starts. On the one hand, the hoop force could directly accelerate the flux rope (Kliem & Török 2006); on the other hand, the torus instability could indirectly contribute to the slow rise by enhancing the HFT reconnection, evidenced by that the HFT reconnection electric field immediately shows a quick increase once the torus instability sets in (Figure 8(h)). This enhancement may be due to the phenomenon that the current structure in the HFT could become longer as the flux rope rises higher and faster with the contribution from the torus instability, which well reflects the close coupling of torus instability and HFT reconnection in the later stage.

A fast magnetic reconnection likely occurs during the impulsive acceleration phase, as a thin and long current sheet, which is usually the necessary condition for the fast magnetic reconnection, is formed at the start of the impulsive acceleration phase (t ∼ 66) and soon built up since then (Figure 9). The other evidence for the fast magnetic reconnection is the impulsively increasing reconnection electric field in the reconnection region during the impulsive acceleration phase (Figure 8(g)). Therefore, we consider that the impulsive acceleration phase (65 < t ≤ 68) is likely triggered as the magnetic reconnection coupled to the torus instability transforms from the relatively weaker HFT reconnection to the relatively stronger fast magnetic reconnection, and that the phase is continuously driven by the close coupling of the fast magnetic reconnection and the torus instability which induces an impulsive acceleration.

Zoom In Zoom Out Reset image size

In this work, we perform three control simulations to help us determine the onset time of the torus instability in "Simulation De." In the following, we first introduce the numerical setups in these three control simulations.

The first control simulation, referred to as "Simulation Half-driven-eruption" (abbreviated as "Simulation He"), is composed of two phases: the shearing phase (0 ≤ t ≤ 18) and the converging phase (18 < t ≤ 70). The setups in "Simulation He" are the same as those in "Simulation De," except that the amplitude of the converging flow is gradually switched off to half during 59 < t ≤ 60 and then kept at half after that (60 < t ≤ 70) by modifying the function of γ(t) in Equation (10) in these periods:

$$\begin{eqnarray}\gamma \,(t)=\left\{\begin{array}{rl}-\displaystyle \frac{1}{4}\;\tanh \,[6.0\,(t-59.5)]+\displaystyle \frac{3}{4} & 59\lt t\leqslant 60\\ \displaystyle \frac{1}{2} & 60\lt t\leqslant 70,\end{array}\right.\end{eqnarray} \tag{ A1 }$$

The second control simulation, referred to as "Simulation Undriven-eruption" (abbreviated as "Simulation Ue"), is composed of three phases: the shearing phase (0 ≤ t ≤ 18), the converging phase (18 < t ≤ 60), and the relaxation phase (60 < t ≤ 71). Before and at t = 59, the setups of "Simulation Ue" are the same as those of "Simulation De." After t = 59, the converging flow is gradually switched off to zero during 59 < t ≤ 60, and then the velocity on the cell-center bottom surface is fixed to zero during the relaxation phase (60 < t ≤ 71), also by modifying the function of γ(t) in Equation (10) in these periods:

$$\begin{eqnarray}\gamma \,(t)=\left\{\begin{array}{rl}-\displaystyle \frac{1}{2}\;\tanh \,[6.0\,(t-59.5)]+\displaystyle \frac{1}{2} & 59\lt t\leqslant 60\\ 0 & 60\lt t\leqslant 71.\end{array}\right.\end{eqnarray} \tag{ A2 }$$

Correspondingly, during 59 < t ≤ 71, the dissipation term in the Equation (3) in the layer k = 1 is modified into:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{v}}{\boldsymbol{B}}-{\boldsymbol{B}}{\boldsymbol{v}}+\psi {\boldsymbol{I}})\\ =\,\left\{\begin{array}{rl}\eta \,\left(\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {y}^{2}}\right)\,{\boldsymbol{B}} & 59\lt t\leqslant 60\\ 0 & 60\lt t\leqslant 71.\end{array}\right.\end{array}\end{eqnarray} \tag{ A3 }$$

During 59 < t ≤ 60, the setup of η is the same as that during 18 < t ≤ 59. During the relaxation phase (60 < t ≤ 71), η is set to zero in the layer k = 1 and set to be uniform (η = 4 × 10⁻⁴; in dimensionless unit) in the whole region except the layer k = 1. All other setups during 59 < t ≤ 71 are the same as those before t = 59.

The third control simulation, referred to as "Simulation No-eruption" (abbreviated as "Simulation Ne"), is also composed of three phases: the shearing phase (0 ≤ t ≤ 18), the converging phase (18 < t ≤ 59), and the relaxation phase (59 < t ≤ 72). The setups in the three phases of "Simulation Ne" are the same as those in the three phases of "Simulation Ue," respectively, except that the converging motion is switched off a little earlier, during 58 < t ≤ 59.

The kinematics of CME progenitors and CMEs in these three control simulations is estimated by a method same as that applied to "Simulation De." For these three simulations, the point, from which the overlying field line is traced, is the same as that used in "Simulation De." The velocity at this fixed point is extremely small, very close to or equal to zero, in all phases of all three simulations.

The kinematics of CME progenitors and CMEs in control simulations and "Simulation De" is shown in Figures A1(a)–(c). We note that, for three control simulations, the decreasing acceleration or even negative acceleration in a period shortly after switching off the driving motion (e.g., 59 ≤ t < 61 in "Simulation He," 59 ≤ t < 61 in "Simulation Ue," 58 ≤ t < 62 in "Simulation Ne") could be due to the viscous and resistive effects, as well as the relative relaxations of magnetic tension and pressure around the tipping point of the equilibrium curve (Démoulin & Aulanier 2010).

Zoom In Zoom Out Reset image size

For all four simulations, in the period 57 ≤ t ≤ 64, the flux rope is basically along the y-direction (Figure A2; taking examples at t = 57 and t = 64 in "Simulation De"). Therefore, the magnetic field at the intersection of the flux rope axis and the plane y = 0 is considered in the y-direction in this period, and the intersection is determined with the contour of B_x = 0 and the contour of B_z = 0 in the plane y = 0. The axis field line is then traced from the intersection. The decay index, $n=-d\,(\mathrm{ln}\;{B}_{t})/d\,(\mathrm{ln}\;h)$, describes the decay rate of the component (B_t ) of the potential field, which is transverse and perpendicular to the flux rope, with the height (h). At each analyzed moment, the potential field is extrapolated by the Green's function method (performed by Solar Software (SSW; Freeland & Handy 2012) package code "optimization_fff.pro") with B_z at the cell-center bottom surface as the input. During 57 ≤ t ≤ 64, the component B_t (h) is defined as the absolute value of B_x of the potential field at (0, 0, h), considering the flux rope axis is along the y-direction, and the horizontal coordinate of its apex is around (0, 0).

Zoom In Zoom Out Reset image size

We first analyze the kinematic evolution in "Simulation Ue" (see red curves in Figures A1(a)–(c)). We consider that the onset time of the CME eruption in the physical sense (i.e., the time when the stable equilibrium of the flux rope is broken) is between t = 59 and t = 64 for the following two reasons: (a) The CME progenitor fails to erupt in "Simulation Ne," in which all setups are the same as those in "Simulation Ue" but the converging motion is switched off during 58 < t ≤ 59. This means that the CME in "Simulation Ue" erupts after t = 58. The earliest limit of its eruption onset time is therefore set at t = 59. (b) The acceleration in "Simulation Ue" changes from negative to positive at t ∼ 64 and increases continuously after that (Figure A1(c)). This indicates that the flux rope is out of equilibrium after that, and thus the latest limit of the eruption onset time is set at t = 64. For "Simulation Ue," the decay index of the potential field at the apex of the flux rope axis is 1.5227 at t = 59 and 1.6239 at t = 64, very close to the threshold of 1.5 for the occurrence of torus instability (Kliem & Török 2006), indicating that the CME eruption is highly probably triggered by the torus instability. The HFT reconnection that starts earlier (t = 46) is considered unable to lead to the eruption as the eruption fails in the presence of HFT in "Simulation Ne" (Aulanier et al. 2010).

Next, we analyze the kinematic evolution in "Simulation De." We consider that the (pre-)eruptive flux rope in "Simulation De" is quite similar to that in "Simulation Ue" in the period 59 ≤ t ≤ 64, as the height of the flux rope axis apex in "Simulation Ue" is only 7.6% lower than that in "Simulation De" at t = 64. Therefore, if the torus instability also sets in "Simulation De," the critical decay index for the torus instability onset should be in a range quite similar to 1.5227 ≤ n ≤ 1.6239. In the following, we take the range 1.5227 ≤ n ≤ 1.6239 for the critical decay index in "Simulation De." In Figure A1(d), it is clear that the height of the flux rope axis apex in "Simulation De" is already higher than the height of the upper limit of the critical decay index at t = 63, indicating that the torus instability could set in "Simulation De" and that its onset time should be before t = 63. The earliest onset time of the torus instability in "Simulation De" is still t = 59, as the apex of the flux rope axis reaches the height of the lower limit of the critical decay index at this moment. Therefore, we summarize that the torus instability sets in "Simulation De" during 59 ≤ t < 63.

We note that it is not an accident that the time range of the torus instability onset in "Simulation De" is earlier than that in "Simulation Ue." In the period 60 < t ≤ 64, the driving motion becomes stronger and stronger from "Simulation Ue" to "Simulation He" to "Simulation De," with the maximum speed of the driving motion varying from 0 to 0.08 to 0.16 (in dimensionless unit). In Figures A1(e), (f), we show that, on the one hand, for these three simulations, the flux rope in the simulation with a stronger driving motion could rise to a higher altitude at the same moment (taking examples at t = 63 (see panel (e)) and at t = 64 (see panel (f))). On the other hand, in comparison among the three simulations, the one with a stronger driving motion has a higher height corresponding to the same critical decay index at the same moment (see panels (e) and (f)). However, it is clear that such a difference in the height of the flux rope axis is much larger than the difference in the height of a certain decay index, when comparing each two of these three simulations. This suggests that in comparison among different simulations, the stronger the imposed driving motion the earlier the flux rope axis reaches the height of the critical decay index, and thus the earlier the torus instability starts, regardless of the change in the critical decay index itself.