Numerical Simulation of a Fundamental Mechanism of Solar Eruption with Different Magnetic Flux Distributions

Following Amari et al. (2003a) and Jiang et al. (2021), we model the photospheric magnetogram of a bipolar field by the composition of two Gaussian functions,

where $\sigma_{x}$ and $\sigma_{y}$ control the extents of the magnetic flux distribution in x and y directions, respectively, and $y_{c}$ control the distance between the two magnetic polarities in the y direction. By adjusting these controlling parameters, we obtain different magnetograms with different flux distributions, some of which have long and strong-gradient PIL while other have short and weak ones. For example, by increasing $\sigma_{y}$ but fixing $\sigma_{x}$ and $y_{c}$, the shape of the polarities gradually changes from flat to round and the strong-gradient PIL becomes shorter, as can be seen in Fig. 1 with four values of $\sigma_{y}$. In another way, by increasing $y_{c}$ but fixing $\sigma_{x}$ and $\sigma_{y}$, the two opposite polarities become further away from each other and thus the field gradient across the PIL decreases, as shown in Fig. 2 for four values of $y_{c}$. For a reasonable comparison, we use the same value for the total unsigned magnetic flux of all the different magnetograms, as such the different magnetograms only differ in the distribution of the same amount of magnetic flux. Therefore, the larger the area of magnetic polarity is, the weaker the average magnetic field strength is, and vice versa. For the specific values of the parameters, we choose $[-10,10]\times[-10,10]$ with length unit of $14.4$ Mm as our computational domain. The values for the parameters as shown in Figs. 1 and 2 are given in their captions.

The BASIC mechanism is closely linked to the fully open field, because by ideally shearing a magnetic arcade, it asymptotically approaches the open field state. In Fig. 3, we show an example of the open field (and compared with the corresponding potential field) for a bipolar magnetogram. The procedure of the open field computation is as follows: first, we calculate the potential field $\textbf{B}_{\rm pot,u}(x,y,z)$ based on a unipolar magnetogram defined by $B_{z,\rm u}(x,y,0)={\rm sign}(y)B_{z}(x,y,0)$. The potential field is solved using the Green’s function method. Thus, all the field lines are open, going outwards from the bottom boundary to infinity. Then, we reverse the direction of the field lines that root in the original negative polarity, and in our case this can be simply done by defining $\textbf{B}_{\rm open}(x,y,z)={\rm sign}(y)\textbf{B}_{\rm pot,u}(x,y,z)$ according to the symmetry of the field. In this way, all the field is still open and current free, except at the interface (i.e., the $y=0$ plane) between the inverse-directed field lines, which forms the CS. Therefore this is the open field corresponding to the bipolar flux distribution.

The open field has important implications on the content of energy that a sheared bipolar field of the same magnetogram can store, as well as the intensity of CS that the sheared bipolar field can form. According to the Aly-Sturrock conjecture, the energy of such open field is the upper limit of the energy of all possible force-free fields with a given magnetic flux distribution on the bottom and a simply connected topology (Aly 1991; Sturrock 1991). Therefore, the upper limit of free magnetic energy that can be reached in a sheared arcade is the open field energy subtracted by the corresponding potential energy of the same magnetogram,

where $E_{\rm uf}$ denotes the upper limit of the free energy. If using the ratio of free magnetic energy to the potential energy as a measure of the non-potentiality of the field, namely $N=E_{\rm free}/E_{\rm pot}$, the upper limit of the non-potentiality of a sheared arcade can achieve is $N_{\rm max}=E_{\rm uf}/E_{\rm pot}$. The degree of non-potentiality $N$ of a field has been suggested as a critical factor in producing major eruptions (Moore et al. 2012; Sun et al. 2015), and thus $N_{\rm max}$ is an important indicator for the eruption capability of a magnetogram. For example, Moore et al. (2012) studied a large number of active regions and found that “there is a sharp upper limit to the free energy the field can hold that increases with the amount of magnetic field in the active region, the active region’s magnetic flux content, and that most active regions are near this limit when their field explodes in a coronal mass ejection/flare eruption.” In particular, using a proxy of magnetic shear for the free energy, they concluded that the non-potentiality $N$ is on the order of one for the active region’s core field, i.e., field rooted around the flare PIL, when the field is close to eruption. Therefore, for a sheared bipolar field that can produce major eruption, it should have its upper limit of non-potentiality somewhat close or above one, i.e., $N_{\rm max}\gtrsim 1$, and by calculating this parameter for a given magnetogram, one can charge whether it is capable of generating large eruption or not.

The potential field $\textbf{B}_{\rm pot}$ and the open field $\textbf{B}_{\rm open}$ are uniquely defined by

and an asymptotic decay at infinity. These fields have energies given, respectively, by the standard relations (e.g., Amari et al. 2003a)

The upper limit of free magnetic energy can be obtained by Eq. (2), and then the non-potentiality $N_{\rm max}$ of each magnetogram can be calculated. By the way, the magnetic field (potential field and open field) in volume can be obtained from the bottom magnetogram, and then the magnetic energy can be obtained by integration $E=\frac{1}{8\pi}\int B^{2}{\rm d}V$. However computing the magnetic field in the full volume is very time consuming. For example, with the Green’s function method, the computing time scales with the grid number as $N^{5}$ (assuming the volume is a cube and the length of each side is $N$). That is why we choose to use Eq. (4) which is much faster in calculation, since the magnetic energy can be obtained directly from the magnetogram without knowing the magnetic field in the full volume. The magnetic energy obtained by these two methods is basically the same. For the magnetogram in Fig. 1A, the potential field energy and open field energy obtained by using the first method are $8.923\times 10^{29}$ erg and $2.059\times 10^{30}$ erg, respectively, and the ones obtained by Eq. (4) is $8.603\times 10^{29}$ erg and $2.090\times 10^{30}$ erg, respectively, with relative errors of $3.586\%$ and $1.483\%$, respectively.

In the open field, all the magnetic field lines have one end rooted at the bottom surface and the other end extending up to infinity, so the field lines on the two side of the PIL run antiparallel, and in between a CS forms (see Fig. 3B and C). The current density outside CS is zero, and all free magnetic energy is stored through the CS (but not in the CS). In the 3D prospective view of the open field as shown in Fig. 3, the structure of the CS is denoted by the red iso-surface of current density $J=0.34\times 10^{-3}$ A m^-2. The central cross section of the open field and the profile of the current density along the central vertical line are also shown in Fig. 3. Since the magnetic field is discretized with a finite resolution, the current density in the CS is not infinite but rather changes with height; it first increases from nearly zero, reaches a peak value at a certain height, and then decrease towards zero again. Therefore, the maximum current density and its height can be obtained by calculating the current density in the CS. We consider that the CS of open field can be used as a proxy of the CS formed in the core field in our simulations, and particularly, the location of the maximum current density in the open field CS indicates the position where reconnection most likely starts to trigger an eruption, and the maximum current density itself should be related to the reconnection rate. Therefore, the two parameters; the peak value of current density of the open field CS and its height, are essential indicator of the strength of the CS, which is correlated with the strength of the eruption once being triggered. In order to calculate the current density of the open field CS, we first calculate the open field by the Green’s function method and then use the second-order central difference to get the current density in the CS.

Figure 4 shows the result for different magnetograms specified by the parameter $\sigma_{y}$ changing from $0.5$ to $2.5$ with increment of $0.1$ (while $\sigma_{x}=2.0$ and $y_{c}=0.8$ are fixed), and Fig. 5 shows the results for magnetograms with parameter $y_{c}$ changing from $0.1$ to $2.1$ with increment of $0.1$ (with $\sigma_{x}=2.0$ and $\sigma_{y}=2.0$ fixed). Specifically, Fig. 4A shows the variation of different magnetic energies and the non-potentiality $N_{\rm max}$ with parameter $\sigma_{y}$. With the increase of $\sigma_{y}$, all energies decrease, and $N_{\rm max}$ also shows an overall decrease from about $1.5$ to below $1.0$, which clearly shows that, overall, the more elongated the magnetic polarity is, the more the upper limits of free energy and non-potentiality the magnetogram can reach (except that a too small $\sigma_{y}$ actually reduces $N_{\rm max}$, which is discussed below). Figure 5A shows the variation of magnetic field energies and the non-potentiality $N_{\rm max}$ with parameter $y_{c}$. A similar decrease pattern is also seen (except that the potential energy increases mildly) with the increase of $y_{c}$, which shows that the closer the two magnetic poles are, the more free magnetic energy (and non-potentiality) the magnetogram can reach at most.

We note that, in Fig. 4A, $N_{\rm max}$ has a maximum at $\sigma_{y}\sim 0.8$ and a smaller $\sigma_{y}$ actually reduces $N_{\rm max}$. This is because as the polarity distance is fixed at $y_{c}=0.8$, a too small value of $\sigma_{y}\leq 0.8$, i.e., a too much concentration of the polarity in the $y$ direction, will actually make strong magnetic flux distributed farther away from the PIL, and thus the magnetic field gradient near the PIL is reduced. It hints that the gradient on the PIL is a more important indicator for $N_{\rm max}$. To confirm this, we calculate the line integral of gradient of the vertical field $B_{z}$ across the PIL for each magnetogram, named as $L$, which is defined by

Figure 6A shows the diagram of $L$ versus $N_{\rm max}$ for all the different sets of parameters shown in Figs. 4 and 5. Strikingly, $N_{\rm max}$ increases monotonically with the increase of $L$, which conforms our argument. Hereafter we will refer to the parameter $L$ as simply the strength of the PIL. Note that $L$ for $\sigma_{y}=0.5$ is less than that of $\sigma_{y}=0.8$, and thus $N_{\rm max}$ of the former is less than that of the latter. We also note that with decreasing of $L$, $N_{\rm max}$ decreases more and more fast.

For the characteristic of the CS, Fig. 4B shows that with the increase of $\sigma_{y}$, the peak value of current density decreases and its height increases. This indicates that the more elongated the magnetic polarity is, the stronger the CS can form. In Fig. 5B, with the increase of $y_{c}$, the maximum current value decreases and its height increases. This shows the closer the magnetic polarities are, the stronger the CS can form. Again, the PIL strength $L$ is crucial. As can be seen in Fig. 6, with the increase of $L$, the peak value of current density increases monotonically overall (except some points due the discretization errors), and the height decreases systematically, which means that the PIL strength is positively correlated with the strength of the CS.

For each magnetogram, we carried out an MHD simulation which begins with the potential field and is driven continually by a rotational flows at each magnetic polarity at the lower boundary, which creates magnetic shear along the PIL. The rotational flow is defined as

with $\psi$ given by

where $B_{z,{\rm max}}$ is the largest value of the photosphere $B_{z}$, and $v_{0}$ is a constant for scaling such that the maximum of the surface velocity is $4.4$ km s^-1, close to the typical flow speed in the photosphere ($\sim$$1$ km s^-1). The flow speed is smaller than the sound speed by two orders of magnitude and the local Alfv$\acute{\text{e}}$n speed by three orders, respectively, thus representing a quasi-static stress of the coronal magnetic field. The flow pattern has been shown in Figs. 1 and 2. This is an incompressible, anti-clockwise rotational flow that does not change with time, and it will not modify the flux distribution at the bottom surface.

We numerically solve the full MHD equations with both coronal plasma pressure and solar gravity included, in a 3D Cartesian geometry by an advanced conservation element and solution element (CESE) method implemented on an adaptive mesh refinement (AMR) grid (Jiang et al. 2010; Feng et al. 2010; Jiang et al. 2016, 2021). Since the controlling equations, the numerical code and essentially all the setting of initial and boundary conditions are the same as used in Jiang et al. (2021), the readers are referred to that paper (in particular, the Method part) for more details. We note that no explicit resistivity is used in the magnetic induction equation, but magnetic reconnection can still occur due to numerical diffusion when a current layer is sufficiently narrow such that its width is close to the grid resolution.

The computational volume spans a Cartesian box of approximately $(-270,-270,0)$ Mm $\leqslant(x,y,z)\leqslant(270,270,540)$ Mm (where $z=0$ represents the solar surface). The volume is large enough such that the simulation runs can be stopped before the disturbance from the simulated eruption reaches any of these boundaries. The full volume is resolved by a block-structured grid with AMR in which the base resolution is $\Delta x=\Delta y=\Delta z=\Delta=2.88$ Mm, and the highest resolution of $\Delta=360$ km is used to capture the formation of CS and the subsequent reconnection.

We have run in total nine experiments for a selected set in the parameter space as listed in Table 1. Specifically, CASE I to CASE IV represent $\sigma_{y}=0.5$, $1.0$, $1.5$ and $2.0$, respectively (with fixed $\sigma_{x}=2.0$ and $y_{c}=0.8$). CASE V to CASE VIII represent $y_{c}=0.1$, $0.6$, $1.1$ and $1.6$, respectively (with fixed $\sigma_{x}=2.0$ and $\sigma_{y}=1.0$). In order to verify the analysis of free magnetic energy maximum in Fig. 4, we set CASE IX, in which $\sigma_{x}=2.0$, $\sigma_{y}=0.8$ and $y_{c}=0.8$.

The BASIC scenario as demonstrated in Jiang et al. (2021) is that, by continually shearing a bipolar coronal field, a CS forms slowly within the arcade and once reconnection sets in, the whole arcade explodes and forms a fast-ejecting magnetic flux rope (i.e., CME). Here we first show briefly an example of such process in CASE I, and then analyze the different runs based on the different magnetograms.

Figure 7 and Supplementary Movie 1 show evolution of magnetic field lines, current density, and velocity during the whole simulation process. The time unit is $\tau=105$ s (all the times mentioned in this paper are expressed with the same time unit). After a period of surface flow driving, the magnetic field structure has evolved from the initial potential field to a configuration with a strong shear immediately above the PIL, where the current density is enhanced, while the envelope field is still current-free ($t\leq 76$). It can be clearly seen that the entire magnetic field structure inflates during the energy injection stage, since the magnetic pressure of the core field increases gradually by the continuous shear. As a result, it stretches outward the envelope field, making the bipolar magnetic arcade tend to approach an open field configuration (Fig. 7B). During this quasi-static evolution process, the current is squeezed from the volumetric distribution into a vertical, narrow layer extending above the PIL, forming a vertical CS (Fig. 7C).

As a critical point, when the thickness of the CS decreases down to the grid resolution, magnetic reconnection sets in and triggers an eruption. This transition from the pre-eruption to eruption onset is clearly manifested in the evolution of energies as shown in Fig. 8 (see the curves colored in magenta), which have a sharp transition at $t=78$. The kinetic energy increases impulsively to nearly $7\%$ original magnetic potential energy in a time duration of $\Delta t=5$. Meanwhile, the magnetic energy releases quickly during the eruption. The onset of the eruption can be more clearly shown by the time profiles of the magnetic energy release rate and the kinetic energy increase rate, and both of them have a sharp increase at the beginning of the eruption (Fig. 8B).

With the onset of reconnection, a plasmoid (i.e., MFR in 3D) originates from the tip of the CS and rises quickly, leaving behind a cusp structure separating the reconnected, post-flare loops from the un-reconnected field (Fig. 7C and D). The plasmoid expands quickly and meanwhile, an arc-shaped fast magnetosonic shock is formed in front of the plasmoid. All of these evolving structures are proof of the typical coronal magnetic eruption leading to CME. The shock marks the front edge of the CME, and its average speed is about $603.4$ km s^-1 (Figs. 7D and 10).

Figures 8 and 9 show the temporal evolution of magnetic and kinetic energies (and their changing rate) of all the cases from CASE I to IX. Note that in each CASE, the energies are normalized by the corresponding potential field energy, i.e., the magnetic energy at $t=0$. Overall, all the different runs (except CASE VIII) show a similar evolution pattern: magnetic energy first increases monotonically for a long time, approaching the open field energy, while the kinetic energy remains to be a very low level; at a critical point, the magnetic energy begins to decrease rapidly along with an impulsive rise of the kinetic energy and the evolutions of the two energies are closely correlated in time, which indicates that the free magnetic energy is released to accelerate the plasma. In Supplementary Movies 2 and 3, we show the evolution of current density on the central cross section for all the cases. As can be seen, they all follow the same BASIC scenario; that is, first a CS forms during the magnetic energy increasing phase and then reconnection sets in and triggers eruption. Therefore, these different runs demonstrate the robustness of the BASIC mechanism.

Nevertheless, the magnitudes, or intensities, of the eruptions in the different cases are different, as can be seen by comparing the energy conversion rates during the impulsive phase. For example, in the five experiments with increasing $\sigma_{y}$ as shown in Fig. 8, the maximum release rate of magnetic energy increases first, reaching the largest at $\sigma_{y}=0.8$ and then decreases with higher $\sigma_{y}$, which is exactly consistent with the dependence of $N_{\rm max}$ on $\sigma_{y}$ as shown in Fig. 4. Again, Fig. 9 shows that with the increase of $y_{c}$, the eruption intensity decreases, consistent with the dependence of $N_{\rm max}$ on $y_{c}$ as shown in Fig. 5. In Fig. 11A and B, we further show how the strength of the PIL, i.e., $L$, is related with the intensity of the eruption, as quantified by the peak values of the kinetic energy increasing rate and the magnetic energy releasing rate as well as the speed of the leading edge of the CME (i.e., the shock). It clearly shows that the eruption intensity is correlated positively with the PIL strength. And in Fig. 11C, we show the non-potentiality $N$ at the onset time of the eruption, as compared with the corresponding $N_{\rm max}$. The non-potentiality increases overall with increase of the PIL strength, consistent with (but a bit slower than) that of the $N_{\rm max}$. Interestingly, we note that the non-potentiality at the eruption onset is mostly close to one, which strikingly agrees with the statistical study of observed active regions (Moore et al. 2012), and hints that the BASIC mechanism is responsible for those eruptions.

Figure 12 shows the central vertical cross section of current density $J$ (normalized by magnetic field strength $B$) at time immediately close to the eruption onset time for the different experiments. The location of the CS and the maximum current density in the CS are also shown in Fig. 11D. These results are consistent with the calculation of the open field in Section 2.2. Furthermore, we can find that the lower the position of the CS and the higher the current density are, the greater the explosive intensity is (Fig. 11).

The only case in our experiments that did not reach an eruption is CASE VIII. This is because during the quasi-static shearing process, the field expands fast in the later phase (e.g., $t>100$) and strongly presses the numerical boundaries before a CS is formed (or before a sufficient amount of free magnetic energy is accumulated to approach a open field), therefore we have stopped the simulation run to avoid a too much influence from the numerical boundaries on the results. Ideally, with a larger computational box, a free expansion of the field driven by the surface shearing motion will create a CS, but its height is too large (and the current density is too low) to trigger an efficient eruption.

It is interesting that the ratio $E_{\rm open}/E_{\rm pot}$ of around $1.7$ is apparently a threshold to start an eruption. Such a similar ratio has been found in Amari’s simulations on flux rope formation and instability (Amari et al. 2000, 2003a, 2003b). Actually, for all the different distributions of magnetic flux as we have considered in this paper, the lowest ratio of $E_{\rm open}/E_{\rm pot}$ is $1.7$ (see Fig. 6A, where $N_{\rm max}=E_{\rm open}/E_{\rm pot}-1$). The reason why $E_{\rm open}/E_{\rm pot}$ of $1.7$ appears to be a threshold to start an eruption is that in our scenario the CS can only form (and thus to trigger an eruption) when the magnetic field is sufficiently sheared such that its energy is close to the open field energy. Furthermore, as we have analyzed, the $E_{\rm open}/E_{\rm pot}$ (or $N_{\rm max}$) should be large enough to let the CS to form at a low height and with a large current density, such that the reconnection can be efficient to produce an eruption. Otherwise, if the $E_{\rm open}/E_{\rm pot}$ is too small, the free energy that can be attained is small, and the CS will form at a too large height (and with too low current density) to trigger an efficient eruption, as our experiment CASE VIII shows. The smaller the ratio $E_{\rm open}/E_{\rm pot}$ is, the higher the CS will form, and the smaller the free energy can be reached, and thus the less efficiently the eruption can be produced. For example, in the extreme case when $E_{\rm open}/E_{\rm pot}$ reaches it a lower limit of $E_{\rm open}/E_{\rm pot}=1$, which corresponds to the flux distribution with two opposite polarities being infinitely far away from each other, it cannot produce any eruption no matter how large the polarities are rotated.

Although using a different scenario (i.e., first a flux rope is created by surface converging and/or cancellation and then the flux rope runs into instability to initiate an eruption), Amari’s simulations also show the same behavior that their eruption occurs only when the magnetic energy is close to the open field energy. Thus, both of our and Amari’s simulations using a single bipolar flux distribution suggest that the ratio $E_{\rm open}/E_{\rm pot}$ should be larger than $1.7$ to start an eruption. Amari’s simulations consist of two important phases of energizing. Their first phase is the same as ours, i.e., by rotating the polarities to inject free energy into the field. The key difference is that in Amari et al. (2000, 2003a, 2003b) simulations the rotation is stopped before a CS is formed; then, in the second phase, they modified the flux content by opposite flux emerging (Amari et al. 2000) or surface diffusion (Amari et al. 2003a), and/or modified the flux distribution by surface converging flow (Amari et al. 2003b) at the bottom boundary. In this phase, the magnetic topology will be changed from a sheared arcade to a flux rope, through the slow reconnection near the bottom boundary. More importantly, the corresponding open field of the evolving magnetic flux distribution will change, being faster than that of the total magnetic energy, and can eventually lead to an eruption when the total magnetic energy is close to the open field energy.

Finally, we note that the eruption onset times of the different experiments are different. This is related to the magnetic energy injection rate, which depends on the surface flow distribution. As shown in Fig. 8, the magnetic energy injection rate decreases when $\sigma_{y}$ increases, and thus the eruption onset time is systematically postponed, because longer time is needed for free magnetic energy accumulation.