A data-driven MHD simulation of the 2011-02-15 coronal mass ejection from Active Region NOAA 11158

The data-driven simulation is carried out using the “Magnetic Flux Eruption” (MFE) code which solves the set of semi-relativistic MHD equations as described in Fan (2017, hereafter F17). The readers are referred to that paper for a description of the equations solved and the numerical methods. Here we only describe the changes made specifically for the setup of the current simulation.

As described in F17, the momentum equation includes the Boris correction with a reduced speed of light to limit the Lorentz force, and hence relax the stringent Courant condition on the numerical time-stepping due to the extremely high Alfvén speed present in the simulation domain which contains a strong active region. For this simulation, we have used a reduced speed of light of about $8000$ km/s, which remains significantly higher than the peak plasma flow speed and the the sound speed reached in the simulation. For the thermodynamics, we assume an ideal gas of fully ionized hydrogen, with $\gamma=5/3$, and solve the internal energy equation taking into account the following non-adiabatic effects that include the field aligned thermal conduction, optically thin radiative cooling, and coronal heating. However, we no longer include an empirical coronal heating (eq. (14) in Fan, 2017) for the heating term in the internal energy equation (term H in eq. (5) of Fan, 2017), but only include the resistive and viscous dissipation due to the numerical diffusion (as described on p.3 of Fan, 2017). For this simulation, we add to the lower boundary a random electric field (as described below in §2.2) representing the effect of turbulent convection that drives field line braiding, and the resultant (numerical) resistive and viscous heating provides the heating of the corona. This heating varies spatially and temporally self-consistently with the formation of the strong current layers in the 3D magnetic field.

The simulation is carried out in a spherical wedge domain whose lower boundary is centered on and tracked the active region. Figure 1 shows the simulation domain against the solar disk as viewed from the Earth perspective on 2011-02-15 at 02:00:00 UT.

It has a latitudinal width of 17.2 degrees, a longitudinal width of 18.8 degrees, and a radial range of 1 to 1.43 solar radius. The simulation domain is resolved with a grid of $360(r)\times 288(\theta)\times 315(\phi)$. The grid is stretched in the radial dimension (with the radial grid size $\Delta r$ increasing outward as a geometric series) and is uniform in the two horizontal dimensions. The peak radial resolution is 300 km at the bottom, and the horizontal resolution is 724 km at the bottom.

For the lower boundary magnetic flux transport driving conditions, we impose a time dependent horizontal electric field ${\bf E}_{h}$ that consists of three components:

The first component ${\bf E}_{h}^{P}$ corresponds to the horizontal component of the so-called “PTD” (Poloidal-Toroidal-Decomposition) electric field derived in Fisher et al. (2020, see the horizontal component of eq. (8) in that paper):

where $\dot{P}$ is computed from the observed normal magnetic field ($B_{r}$) evolution on the photosphere by solving the 2D Poisson equation (eq. (9) in Fisher et al., 2020):

due to Faraday’s law. Imposing ${\bf E}_{h}^{P}$ reproduces the observed photospheric normal magnetic field evolution on the lower boundary.

The second term ${\bf E}_{h}^{\rm Tw}$ on the right hand side of equation (1) is the “twisting” electric field. It is given as the horizontal gradient of a potential:

such that it does not alter the observed $B_{r}$ evolution (already reproduced by imposing the ${\bf E}_{h}^{P}$ component), and we assume that it corresponds to a vertical transport of a horizontal magnetic field into the domain, i.e.

where $v_{0}$ is the vertical transport speed, which we assume to be a constant, and ${\bf{\cal B}}_{h}$ is a horizontal magnetic field, which we note is not the same as the observed horizontal field and it can be shown from equation (5) that $\nabla_{h}\cdot{\bf{\cal B}}_{h}=0$ (which condition is generally not satisfied by the observed horizontal magnetic field). Further taking the horizontal divergence of equation (5) yields:

We let $(\nabla_{h}\times{\bf{\cal B}}_{h})_{r}$ be equal to $J_{r}=(\nabla_{h}\times{\bf B}^{\rm obs}_{h})_{r}$, which is the vertical electric current derived from the observed horizontal magnetic field ${\bf B}^{\rm obs}_{h}$ in the photospheric vector magnetograms. Thus we have

We compute the potential ${\psi}^{\rm Tw}$ by solving the 2D Poisson equation (7) given the measured vertical electric current at the photosphere $J_{r}$ from the HMI vector magnetograms and specifying $v_{0}$, which is an ad hoc parameter we can adjust. For the simulation presented here we have used $v_{0}=2.5\,{\rm km/s}$, which is small enough to ensure a quasi-static evolution for the build-up phase but high enough to compete with numerical diffusion to produce the eruptive behavior close to that observed. Once ${\psi}^{\rm Tw}$ is determined, the twisting electric field ${\bf E}_{h}^{\rm Tw}$ is given by equation (4). Imposing ${\bf E}_{h}^{\rm Tw}$ at the lower boundary corresponds to transporting a (divergence-free) horizontal magnetic field (${\bf{\cal B}}_{h}$) with the observed vertical electric current into the domain, without changing the normal magnetic field evolution at the lower boundary. It effectively transports twist into the corona based on the observed vertical electric current.

Figure 2 shows example snapshots of the observed normal magnetic field $B_{r}$, the vertical electric current density $J_{r}$, and the horizontal twisting electric field components $E_{\theta}^{\rm Tw}$ and $E_{\phi}^{\rm Tw}$ derived based on $J_{r}$ (from eqs. [4] and [7]), at the time of about 1.9 hour before the onset of the observed X-class flare.

From panels (a) and (b) we can see that significant $J_{r}$ which tends to be of the same sign as $B_{r}$ is present in the polarity concentrations P2 and N1 near the central polarity inversion line (PIL). On the other hand, in the polarity concentration N2, we find $J_{r}$ of predominantly opposite sign of $B_{r}$. As a result, the derived twisting electric field (panels (c) and (d)) transports positive twist into the corona at the central P2, N1 flux concentrations and negative twist at the N2 flux concentration. It can be seen from panels (c) and (d) that $E_{\theta}^{\rm Tw}$, and $E_{\phi}^{\rm Tw}$ transport a concentrated horizontal, shear magnetic field component into the corona at the central PIL. They also effectively drive a clockwise rotation of the P2 and N1 flux concentrations and a counter-clockwise rotation of the N2 flux concentration.

We note that Cheung & DeRosa (2012) have also used a horizontal electric field to drive twist into the corona in their magneto-frictional modeling of the evolution of AR 11158. In that case, effectively a uniform rotation of all the polarity concentrations is applied to drive the same sign of twist into the corona. Here with the twisting electric field, the twist injection varies spatially based on the distribution of the observed vertical electric current.

The last term on the right hand side of equation (1) is a random, horizontal electric field $\bf{E}_{h}^{\rm{random}}$ given by (see e.g. Fan, 2022):

Again, this ${\bf{E}}_{h}^{\rm{random}}$ is given as the gradient of a scalar field and thus does not alter the observed $B_{r}$ on the lower boundary (reproduced by imposing ${\bf{E}}_{h}^{P}$). The formulation of this electric field is inspired by that of the “STatistical InjecTion of Condensed Helicity” (“STITCH”) electric field (Mackay et al., 2014; Dahlin et al., 2022) used for modeling helicity condensation and filament channel formation. But here, $\xi$ used in equation (8) is a time dependent field that is made up of a superposition of 15721 randomly placed cells of opposite sign values. Each cell has a 2D Gaussian profile with a size scale of 1.74 Mm and a peak amplitude of $4.07\times 10^{13}\rm{cm^{2}/s}$, and varies temporally as a sinusoidal function with a period and life time of 11.89 min. Figure 2(e) shows a snapshot of the $\xi$ field, which illustrates its random cellular pattern, and the resulting random electric field components $E_{\theta}^{\rm random}$ (Figure 2(f)), and $E_{\phi}^{\rm random}$ (Figure 2(g)) produced at the lower boundary when this $\xi$ field and the $B_{r}$ given in Figure 2(a) are applied in equation (8). As is described in Fan (2022), ${\bf{E}}_{h}^{\rm{random}}$ effectively drives random rotations of the foot points of the $B_{r}$ flux concentrations, with positive (negative) $\xi$ corresponding to clockwise (counter-clockwise) rotation. Here the positive and negative signed cells in the $\xi$ field are statistically balanced, therefore no net twist or helicity is driven into the corona by ${\bf{E}}_{h}^{\rm{random}}$. It represents the effect of turbulent convection that drives field-line braiding which produces resistive and viscous heating in the corona. Similar ways of driving coronal heating by imposing random foot-point motions that represent turbulent convection at the photospheric lower boundary of MHD simulations have been widely used (e.g. Gudiksen & Nordlund, 2005; Warnecke & Peter, 2019).

Even though the driving electric field is derived from the photosphere magnetograms, we impose a fixed chromosphere temperature of 20,000 K and density of $10^{12}\,{\rm cm}^{-3}$ at the lower boundary. The boundary conditions for the side and top boundaries in the present simulation are the same as those used in Fan (2017).

To set up the initial state of the data-driven simulation, we start with a potential magnetic field extrapolated from the observed photospheric normal magnetic field $B_{r}$ at 2011-02-14 23:47 UT (Fig. 2(a), about 1.9 hours before the onset of the observed X-class flare), and numerically evolve the MHD state by driving at the lower boundary with only the random electric field ${\bf{E}}_{h}^{\rm{random}}$ until it reaches a quasi-steady state with a hot corona. Figure 3 shows the 3D magnetic field lines (panel (a)) colored with the twist rate defined as $\alpha\equiv{\bf J}\cdot{\bf B}/B^{2}$ where ${\bf J}=\nabla\times{\bf B}$ is the current density, and the radial profiles of the horizontally averaged temperature (panel (b)), density (panel (c)), and the Alfvén and sound speeds (panel (d)) of the relaxed, quasi-steady state reached. It can be seen that the magnetic field (panel (a)) for the relaxed state is no longer the potential field but contains some random twist $\alpha$ due to the driving random electric field (which does not drive a net helicity into the domain), although most field lines remain near zero $\alpha$ (white color), and the magnetic energy remains very close to that of the potential field energy (being about 1.0016 of the potential field energy). A quasi-steady atmosphere with a reasonable temperature and density stratification (panels (b) and (c)) that contains a chromosphere transitioning into a hot corona is established. The Alfvén speed is significantly higher than the sound speed throughout the domain (panel (d)).

We then use this relaxed state as the initial state (defined hereafter as $t=0$) for the simulation, driven at the lower boundary with all three electric fields given in the right hand side of equation (1) (derived from a time sequence of the vector magnetograms starting from 2011-02-14 23:47 UT) for a perid of over 2.8 hours. The following section describes the resulting evolution obtained from the boundary data-driven simulation.

Figure 4 shows the resulting evolution of the 3D magnetic field (left column images), the radial velocity in the central meridional cross-section (middle column), and several integrated quantities (right column) that include the free magnetic energy $E_{\rm free}$ (panel (g)), which is the excess of the total magnetic energy over that of the corresponding potential magnetic field, the total kinetic energy $E_{k}$ (panel (h)), and the various contributions to the rate of change of the total magnetic energy (panel (i)) which are described in the following.

(1) The black curve in panel (i) shows the input of magnetic energy by the Poynting flux integrated over the lower and upper boundaries: $\int(c/4\pi)({\bf E}_{h}\times{\bf B}_{h})_{r}\,dS$, where $(c/4\pi)({\bf E}_{h}\times{\bf B}_{h})_{r}$ is the radial component of the Poynting flux density and the integration $dS$ is over the area of the lower and upper boundaries (note no Poynting flux through the side boundaries due to the conducting wall side boundary condition). (2) The solid red curve shows the release of magnetic energy resulting from the Lorentz force work: $-\int(1/4\pi)((\nabla\times{\bf B})\times{\bf B})\cdot{\bf v}\,dV$, where $(1/4\pi)((\nabla\times{\bf B})\times{\bf B})$ is the Lorentz force, ${\bf v}$ is the velocity, and the integration $dV$ is over the volume of the domain. (3) The red dash-dotted curve shows the dissipation of magnetic energy by the numerical magnetic diffusion: $-\int(c/4\pi){\bf E}_{\rm num}\cdot(\nabla\times{\bf B})\,dV$, where ${\bf E}_{\rm num}$ is the electric field resulting from the numerical diffusion evaluated in the numerical code. The blue curve shows the sum of the above three, which is the net rate of change of the total magnetic energy. We see an overall continuous build up of the free magnetic energy $E_{\rm free}$ (panel (g)) due to the continuous Poynting flux input (black curve in panel (i)) at the lower boundary produced by the driving electric field (${\bf E}_{h}$ given by eq. [1]). This Poynting flux input is in excess of the dissipations of the magnetic energy due to the resistive dissipation (red dash-dotted curve in panel (i)) and the Lorentz force work (red solid curve in panel (i) ), resulting in a continuous net gain of the free magnetic energy (blue curve in panel (i)) for most of the time. But this continuous build up is punctuated by sudden free magnetic energy $E_{\rm free}$ releases and kinetic energy $E_{k}$ increases (at about $t=1.51$ hour, $t=2.1$ hour, and $t=2.7$ hour), due to the sudden enhancements of the Lorentz force work and resistive dissipation resulting from the loss of equilibrium of the magnetic field (see the movie for the whole course of the evolution).

During the first 1.5 hour period, $E_{\rm free}$ increases steadily, while $E_{k}$ and the Lorentz force work remain small with the magnetic field being in quasi-equilibrium. We can see from the movie of Figure 4 that during this time the magnetic field is being sheared and twisted as indicated by the color of the field lines, forming positive-twisted (red field lines) sigmoidal shaped loops above the central (PIL) and negative-twisted (blue field lines) inverse-S shaped loops connecting to the following negative polarity sunspot (panel (a)). At about $t=1.506$ hour, the free magnetic energy reaches a peak value of about $2.25\times 10^{32}{\rm erg}$ (panel (g)), and then we see the onset of the first eruption with a sudden outward acceleration of the central sigmoid field (panels (b) and (e)), forming a positive-twisted erupting flux rope (red twisted field lines) which also pushes out and accelerates an outer negative-twisted field (blue field lines) (panels (c) and (f) and the movie).

Figure 5 shows the pre-eruption magnetic field just before the onset of the eruption at $t=1.506$ hour. The top two panels show a set of 3D magnetic field lines colored with the current density $J$ (panel (a)) and the twist rate $\alpha$ (panel (b)).

Comparing the pre-eruption field in Figure 5(b) to the nearly potential field of the initial state in Figure 3(a), we find that the twisting electric field has built up highly sheared, forward S-shaped (sigmoid) loops with positive (right-handed) twist rates $\alpha$ above the central PIL connecting polarities P2, N1, and also positive-twisted loops connecting polarities P1, N1. Furthermore, we find negative-twisted (left-hand-twisted) loops (blue with negative $\alpha$) connecting P2, N2 polarities, and some higher loops with negative twist connecting P1, N2 polarities. As pointed out earlier in Figure 2, the twisting electric field imposed in N2 (with opposite signs of $J_{r}$ and $B_{r}$) transports negative twist into the corona, opposite to the twisting electric field dominating in the other polarity concentrations that drives positive twist into the corona.

It can be seen in Figure 5(a) that the field lines with the strongest current density $J$ show morphology in qualitative agreement with the brightening loops in the observed hot channel images of SDO/AIA (Figs. 5(e) and 5(f)). This indicates that the model captures the structure of the non-potential (energized) magnetic field. Figures 5(c) and 5(d) show respectively the synthetic AIA 131 Å and 94 Å images at the same time of the pre-eruption magnetic field shown in Figures 5(a) and 5(b). The synthetic images are computed by integrating along the line of sight (from the Earth perspective) through the simulation domain:

where $l$ denotes the length along the line of sight through the simulation domain, $I$ denotes the intensity at each pixel in units of DN/s/pixel (shown in LOG scale in the images), $n_{e}$ is the electron number density, and $f(T)$ is the temperature response function of the corresponding AIA filter (obtained by running the SolarSoft routine get_aia_response.pro). The synthetic AIA images (Figs. 5(c) and 5(d)) also show some qualitatively similar emission features as those in the observed ones (Figs. 5(e) and 5(f)), i.e. the central forward S-shaped sigmoid emission along the central PIL, and the surrounding large loops connecting the major polarity concentrations. We find that the central sigmoid emission is due to the heating produced by the formation of a strong current layer in the strongly sheared field above the central PIL.

To examine how close to force-free the pre-eruption magnetic field is, we have evaluated $\sigma_{J}$, the current-weighted mean of the sine of the angle $\theta$ between the current density ${\bf J}$ and the magnetic field ${\bf B}$ (Wheatland et al., 2000):

where

the subscript “i” denotes each grid point and the sum is over all the grid points in the simulation domain. We found $\sigma_{J}\approx 0.06\ll 1$, corresponding to a (current-weighted) mean angle of about $3.4^{\circ}$ between the current density and the magnetic field vectors, i.e. close to being force-free. Furthermore we found that the vertical current density obtained at the simulation lower boundary is close to the observed $J_{r}$ from the HMI vector magnetogram at the corresponding time, for the large-scale main current patches, as can be seen in Figure 6 (compare panels (b) and (c)).

This shows that the modeled pre-eruption magnetic field is close to the solution of a NLFFF extrapolation from the observed vector magnetogram. Indeed, the peak free energy of the pre-eruption magnetic field reached just before the first eruption ($2.25\times 10^{32}{\rm erg}$) is close to the result obtained in Sun et al. (2012), who carried out the NLFFF extrapolation of the coronal magnetic field of AR 11158 and found that the free magnetic energy reaches a maximum of $\sim 2.6\times 10^{32}{\rm erg}$ just befored the onset of the observed X2.2 flare. As described in §2.2, the twisting electric field we impose does not enforce the observed horizontal magnetic field at the lower boundary. Nevertheless, we find that the resulting horizontal magnetic field at the lower boundary for the pre-eruption magnetic field that is built up is very similar to that of the observed.

Figure 7 shows a comparison of the horizontal components of the magnetic field at the lower boundary at the onset of the simulated eruption with those from the photospheric vector magnetogram at the same time. We see very similar patterns for both the $B_{\theta}$ and $B_{\phi}$ components comared to the observed ones.

Figure 8 shows the 3D evolution of a set of erupting magnetic field lines as viewed from 2 different perspectives (top 2 rows) during the first eruption. The field lines are traced from a set of Lagrangian tracer points tracked in the velocity field. The eruption begins with the acceleration of the highly sheared sigmoid field with positive twist above the central PIL (see the red field lines in panel (a) and the sudden onset of the rise velocity from panel (c) to panel (f) of Fig. 8). To examine what might have triggered the sudden acceleration of the sigmoid field above the central PIL and whether it contains a magnetic flux rope, we plotted in Figure 9 the distribution of twist rate $\alpha$ and the radial curvature ${\bf B}\cdot\nabla(B_{r})/B^{2}$ in the vertical cross-section, whose location is indicated by the black line across the central PIL in Figure 6(a), at the time of the onset of the eruption ($t=1.506$ hour, top row) and at a time during the eruption($t=1.545$ hour, bottom row).

It can be seen in the top left panel that most of the current of the positive-twisted (red) sigmoid field region is still below the contour of the critical decay index of $d({\rm ln}B_{p{\rm h}})/d({\rm ln}h)=-1.5$ for the onset of the torus instability (e.g. Kliem & Török, 2006). Here $B_{p{\rm h}}$ is the horizontal field strength of the corresponding potential field and $h$ denotes the height above the surface. In the mean time, it can also be seen that a thin current layer (thin layer of strong positive $\alpha$) has developed just above the PIL (located at about $Y=-0.01$ in Figure 9) in the sigmoid field region, which indicates the onset of the tether-cutting reconnection at the current layer may be the trigger of the onset of the eruption instead of the torus instability. Furthermore, we find no dipped field with positive radial curvature above the central PIL (at $Y=-0.01$) in the top right panel of Figure 9, indicating that at the onset of the eruption, the sigmoid field above the PIL is a sheared arcade instead of a flux rope. The bottom panels of Figure 9 corresponding to shortly after eruption onset show that a flux rope with dipped field lines with positive radial curvature forms during the eruption (bottom right panel) as a result of the tether-cutting reconnection in the vertical current sheet above the central PIL (bottom left panel).

The onset of the tether-cutting reconnection is facilitated by the numerical resistivity, triggered by the thinning of the current sheet in the strongly sheared arcade field above the PIL. Figure 10 shows a 3D view of example field lines before and after undergoing tether-cuttig reconnections. We see that a pair of highly sheared arcade field lines above the central PIL in the earlier time instance (left panel) have reconnected and transformed into a longer dipped field line that erupts upward and a lower short loop that shrinks back down in the later time instance (right panel). The formation of a flux rope with dipped field lines can also be seen in panels (d) and (g) in Figure 8. As the positive-twisted inner field erupts, it encounters and reconnects with the negative-twisted field above (see top row of Fig. 8). Eventually, the erupting field develops a complex structure that consists of an outer flux rope containing mainly negative twist (left-handed twist) and an inner flux rope containing mainly positive twist (right-handed twist) (see panels (m) and (n) in Fig 8), with the outer flux rope accelerating to a speed of about $800$ km/s as it approaches the top of the domain (panels (l) and (o) in Fig. 8). It can also be seen in panels (l) and (o) in Fig. 8 (and the associated movie) that the rise velocity of the inner flux rope appears to decelerate towards the end. This is due to its collision and reconnection with the overlying negative twisted flux, which forms the outer flux rope that continues to accelerate and exits the domain while the inner flux rope slows down as it loses twist via reconnection with part of the outer magnetic flux. The top boundary is an open boundary allowing plasma to flow through. However, the side boundaries are closed wall boundaries where both the magnetic and velocity fields become parallel to the walls. Thus the lateral expansion and deflection of the erupting flux ropes are artificially constrained by the side boundaries (see result in the next section).

To study the observational signatures of the simulated erupting field, we have computed the synthetic EUV images in the AIA 171 Å channel but viewed from the STEREO-B perspective. The right panel of Figure 11 shows the running difference of two synthetic images at times $t=1.555$ hour and $t=1.565$ hour during the eruption. It shows an outgoing double shell structure produced by the two outgoing flux ropes. A running difference image from the observed STEREO-B EUVI 171 passband images during the eruption (left panel of Fig. 11) shows a similar outgoing double-shell structure, although the observed out-going structure expands side-way significantly more than the simulated one because of the constraining side wall boundary of the simulation domain. Simulation with a widened domain is needed for a quantitative comparison of the kinematics of the erupting structure.

CME-associated coronal dimmings are thought to correspond to footpoints of the expanding CME magnetic structures and caused by a density decrease in these structures due to their rapid expansion (e.g. Thompson et al., 1998; Dissauer et al., 2018). Here we use our data-driven simulation to identify the footpoints of the erupting field lines and compare them with the observed EUV dimmings of the CME event. The bottom panels in Figure 12 show the lengths of the field lines traced from each locations at the base of the corona at the onset of the eruption (panel (d)), and at a time when the outer erupting flux rope is exiting the domain (panel (e)), and the difference between the two times (panel (f)) which shows the change of the field line lengths due to the eruption.

The regions of large field line length increases (as represented by the red patches in panel (f)) correspond to the foot points of the stretched-out erupting field lines. We find that these regions of the erupting field line foot points correspond well with both the two core dimming regions (indicated by the two short arrows) as well as the the more diffuse secondary dimming region (indicated by the long arrow) seen in the AIA 211 Å channel observations (panels (b) and (c)). This agreement adds validation to the modeled erupting field.

Observations have also shown that there is a rapid and irreversible change of the photospheric magnetic field associated with a solar flare (e.g Wang et al., 1994; Sun et al., 2012; Liu et al., 2012) as a result of the “implosion” of the coronal magnetic field just above the flare PIL due to the release of the magnetic energy (Hudson, 2000; Fisher et al., 2012). In Figure 13 we show the change of the horizontal magnetic field strength $\Delta B_{h}$ at the lower boundary 12 min after the onset of the eruption in the simulation compared to the observed change of the horizontal magnetic field strength in the photosphere vector magnetogram 12 min after the onset of the observed X-class flare.

We find that the change of $B_{h}$ in the simulation shows a similar pattern as the observed change in the central flaring region, where there is a significant enhancement of $B_{h}$ over the central PIL due to the implosion of the post reconnection loops, while there is a systematic decrease of $B_{h}$ in the periphery of the flare region away from the central PIL, where the field becomes more vertical as it erupts. This agreement is further evidence that our simulated eruption qualitatively reproduces the magnetic field reconfiguration during the flare.