Numerical Simulations of Flare-productive Active Regions: δ-sunspots, Sheared Polarity Inversion Lines, Energy Storage, and Predictions

In this paper, we investigate the emergence of buoyant flux tubes that are initially set in the convection zone. We considered a rectangular computational domain with three-dimensional (3D) Cartesian coordinates $(x,y,z)$, where the z-coordinate increases upward. We solved the standard set of resistive magnetohydrodynamic (MHD) equations:

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(\rho {\boldsymbol{V}})-{\boldsymbol{\nabla }}\cdot \left[\rho {\boldsymbol{V}}{\boldsymbol{V}}+\left(p+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi }\right){\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{B}}{\boldsymbol{B}}}{4\pi }\right]+\rho {\boldsymbol{g}}=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=c{\boldsymbol{\nabla }}\times {\boldsymbol{E}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left(\rho U+\displaystyle \frac{1}{2}\rho {{\boldsymbol{V}}}^{2}+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi }\right)\\ &&\qquad +{\boldsymbol{\nabla }}\cdot \left[\left(\rho U+p+\displaystyle \frac{1}{2}\rho {{\boldsymbol{V}}}^{2}\right){\boldsymbol{V}}+\displaystyle \frac{c}{4\pi }{\boldsymbol{E}}\times {\boldsymbol{B}}\right]\\ &&\qquad -\,\rho {\boldsymbol{g}}\cdot {\boldsymbol{V}}=0,\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&U=\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{p}{\rho },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\boldsymbol{E}}=-\displaystyle \frac{1}{c}{\boldsymbol{V}}\times {\boldsymbol{B}}+\eta {\boldsymbol{J}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\boldsymbol{J}}=\displaystyle \frac{c}{4\pi }{\boldsymbol{\nabla }}\times {\boldsymbol{B}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&p=\displaystyle \frac{{k}_{{\rm{B}}}}{m}\rho T,\end{eqnarray} \tag{ 8 }$$

where ρ denotes the gas density, ${\boldsymbol{V}}$ the velocity vector, p the pressure, ${\boldsymbol{B}}$ the magnetic field, c the speed of light, ${\boldsymbol{E}}$ the electric field, and T the temperature, while U is the internal energy per unit mass, ${\boldsymbol{I}}$ the unit tensor, ${k}_{{\rm{B}}}$ the Boltzmann constant, $m(=\mathrm{const}.)$ the mean molecular mass, and ${\boldsymbol{g}}=(0,0,-{g}_{0})=(0,0,-1/\gamma )$ is the uniform gravitational acceleration. We assumed the medium to be an inviscid perfect gas with a specific heat ratio $\gamma =5/3$.

To make the above equations dimensionless, we introduced the following normalizing units: the pressure scale height ${H}_{0}=170\,\mathrm{km}$ for the length, the sound speed ${C}_{{\rm{s}}0}=6.8\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for the velocity, ${\tau }_{0}\equiv {H}_{0}/{C}_{{\rm{s}}0}=25\ {\rm{s}}$ for the time, and ${\rho }_{0}=1.4\times {10}^{-7}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ for the density, all of which are typical values in the photosphere. The gas pressure, temperature, and magnetic field strength were normalized by combinations of the units above, i.e., ${p}_{0}=\rho {C}_{{\rm{s}}0}^{2}\,=6.3\times {10}^{4}\,\mathrm{dyn}\,{\mathrm{cm}}^{-2}$, ${T}_{0}={{mC}}_{{\rm{s}}0}^{2}/(\gamma {k}_{{\rm{B}}})=\mathrm{5,600}\,{\rm{K}}$, and ${B}_{0}\,={({\rho }_{0}{C}_{{\rm{s}}0}^{2})}^{1/2}=250\,{\rm{G}}$, respectively.

We assumed an anomalous resistivity model with the form

$$\begin{eqnarray}\eta =\left\{\begin{array}{ll}0 & (J\lt {J}_{{\rm{C}}}\ \mathrm{or}\ \rho \gt {\rho }_{{\rm{C}}})\\ {\eta }_{0}(J/{J}_{{\rm{C}}}-1) & (J\geqslant {J}_{{\rm{C}}}\ \mathrm{and}\ \rho \lt {\rho }_{{\rm{C}}})\end{array}\right.,\end{eqnarray} \tag{ 9 }$$

where ${\eta }_{0}=0.1$, ${J}_{{\rm{C}}}=0.1$, and ${\rho }_{{\rm{C}}}=0.1$. The above treatment is intended to trigger magnetic reconnection in a low-density current sheet.

The initial background atmosphere consisted of three regions: an adiabatically stratified convection zone, a cool isothermal photosphere/chromosphere, and a hot isothermal corona (see Figure 2). We assumed $z/{H}_{0}=0$ to be the base height of the photosphere, and the initial temperature distribution of the convection zone ($z/{H}_{0}\leqslant 0$) was assumed to be

$$\begin{eqnarray}&&T(z)={T}_{\mathrm{ph}}-z{\left|\displaystyle \frac{{dT}}{{dz}}\right|}_{\mathrm{ad}},\end{eqnarray} \tag{ 10 }$$

where ${T}_{\mathrm{ph}}={T}_{0}$ is the photospheric temperature and

$$\begin{eqnarray}&&{\left|\displaystyle \frac{{dT}}{{dz}}\right|}_{\mathrm{ad}}=\displaystyle \frac{\gamma -1}{\gamma }\displaystyle \frac{{{mg}}_{0}}{{k}_{{\rm{B}}}}\end{eqnarray} \tag{ 11 }$$

is the adiabatic temperature gradient; i.e., the initial temperature profile of the convection zone is adiabatic. The temperature distribution of the atmosphere ($z/{H}_{0}\geqslant 0$) was expressed as

$$\begin{eqnarray}&&T(z)={T}_{\mathrm{ph}}+\displaystyle \frac{1}{2}({T}_{\mathrm{cor}}-{T}_{\mathrm{ph}})\left[\tanh \left(\displaystyle \frac{z-{z}_{\mathrm{cor}}}{{w}_{\mathrm{tr}}}\right)+1\right],\end{eqnarray} \tag{ 12 }$$

where ${T}_{\mathrm{cor}}=150{T}_{0}$ is the coronal temperature, ${z}_{\mathrm{cor}}=18{H}_{0}$ is the base of the corona, and ${w}_{\mathrm{tr}}=2{H}_{0}$ is the temperature scale height of the transition region. With the temperature profile above, the initial pressure and density profiles (Figure 2) were defined by the equation of static pressure balance:

$$\begin{eqnarray}&&\displaystyle \frac{{dp}(z)}{{dz}}+\rho (z){g}_{0}=0.\end{eqnarray} \tag{ 13 }$$

Zoom In Zoom Out Reset image size

A magnetic flux tube was embedded in the convection zone, and its longitudinal and azimuthal components of the flux tube are given by

$$\begin{eqnarray}&&{B}_{x}(r)={B}_{\mathrm{tube}}\exp \left(-\displaystyle \frac{{r}^{2}}{{R}_{\mathrm{tube}}^{2}}\right)\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{B}_{\phi }={{qrB}}_{x}(r),\end{eqnarray} \tag{ 15 }$$

where $r={[{(y-{y}_{\mathrm{tube}})}^{2}+(z-{z}_{\mathrm{tube}}^{2})]}^{1/2}$ is the radial distance from the tube axis, $({y}_{\mathrm{tube}},{z}_{\mathrm{tube}})$ the location of the tube axis, ${R}_{\mathrm{tube}}$ the radius, ${B}_{\mathrm{tube}}$ the magnetic field strength at the axis, and q the twist intensity. As the reference (typical) case, we considered $({y}_{\mathrm{tube}}/{H}_{0},{z}_{\mathrm{tube}}/{H}_{0})=(0,-30)$, ${R}_{\mathrm{tube}}/{H}_{0}=3$, ${B}_{\mathrm{tube}}/{B}_{0}=30$, and ${{qH}}_{0}=-0.2$. The total axial magnetic flux amounts to ${{\rm{\Phi }}}_{\mathrm{tube}}/({B}_{0}{H}_{0}^{2})=845$. These parameters indicate that the initial flux tube is located at a depth of 5.1 Mm and has a radius of 510 km; a central field strength of 7.5 kG (or the plasma $\beta \equiv 8\pi p/{B}^{2}\sim 10$), which yields an axial flux of $6\times {10}^{19}\ \mathrm{Mx};$ and a left-handed twist. The magnetic pressure, ${p}_{\mathrm{mag}}={B}^{2}/(8\pi )$, along the vertical axis is shown in Figure 2. The gas pressure inside the tube was defined as ${p}_{{\rm{i}}}=p(z)+\delta {p}_{\mathrm{exc}}$, with the pressure excess being

$$\begin{eqnarray}&&\delta {p}_{\mathrm{exc}}=\displaystyle \frac{{B}_{x}^{2}(r)}{8\pi }\left[{q}^{2}\left(\displaystyle \frac{{R}_{\mathrm{tube}}^{2}}{2}-{r}^{2}\right)-1\right](\lt 0).\end{eqnarray} \tag{ 16 }$$

To trigger the buoyant emergence, we reduced the density inside the flux tube, ${\rho }_{{\rm{i}}}=\rho (z)+\delta {\rho }_{\mathrm{exc}}$, where

$$\begin{eqnarray}&&\delta {\rho }_{\mathrm{exc}}=\rho (z)\displaystyle \frac{\delta {p}_{\mathrm{exc}}}{p(z)}\left[(1+\epsilon )\exp \left(-\displaystyle \frac{{(x-{x}_{\mathrm{tube}})}^{2}}{{\lambda }^{2}}\right)-\epsilon \right],\end{eqnarray} \tag{ 17 }$$

and ${x}_{\mathrm{tube}}$ and λ are the center and the length of the buoyant section, respectively, and is a factor that suppresses the emergence of both ends of the tube. The typical values are ${x}_{\mathrm{tube}}/{H}_{0}=0$, $\lambda /{H}_{0}=8$, and $\epsilon =0.2$. Depending on the simulation case, some parameters were modified as described in the next subsection.

The simulation domain was (−150, −150, −40) ≤ (x/H₀, y/H₀, z/H₀) ≤ (150, 150, 400), resolved by a 512 × 512 × 512 grid. The grid spacings for the x-, y-, and z-directions were Δx/H₀ = Δy/H₀ = 0.25 for (−20,−20) ≤ (x/H₀, y/H₀) ≤ (20, 20) and ${\rm{\Delta }}z/{H}_{0}=0.2$ for $-40\leqslant z/{H}_{0}\leqslant 15$. Outside this range, the spacings were smoothly increased up to ${\rm{\Delta }}x/{H}_{0}\,={\rm{\Delta }}y/{H}_{0}=0.8$ and ${\rm{\Delta }}z/{H}_{0}=1.8$. We assumed a periodic boundary condition for the x-direction and symmetric boundaries for both the y- and z-directions.

The simulation code we used is the same as that used by Takasao et al. (2015), which is based on the numerical scheme of Vögler et al. (2005): fourth-order central differences for calculating the spatial derivatives, and the four-step Runge–Kutta scheme for calculating the temporal derivatives. Artificial diffusivity, proposed by Rempel et al. (2009), was introduced to stabilize the calculation, while the ${\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}$ error was reduced by the iterative hyperbolic divergence cleaning technique based on the method described in Dedner et al. (2002).

Figure 3 shows the evolution of the magnetogram at $z/{H}_{0}=0$ and magnetic field lines for the reference case (a movie is attached to provide detailed evolution). The whole evolution is consistent with the previous 3D simulations by, e.g., Fan (2001), Archontis et al. (2004), and Toriumi & Yokoyama (2012): the horizontal flux tube makes an Ω-shaped arcade, which rises through the convection zone and eventually penetrates the photosphere, creating a magnetic dome in the corona with bipolar spots in the photosphere.

In order to model the four types of flare-productive ARs introduced in Section 1, we tested four simulation cases with initial conditions different from those of the reference case, which are summarized in the bottom row of Figure 1.

For the spot–spot case, the initial twist strength was intensified to ${{qH}}_{0}=-0.8$, which is higher than the critical value for the kink instability ($| q| {H}_{0}=0.33$: Linton et al. 1996). Owing to the stronger initial twist, the density deficit is higher for this case (Equations (16) and (17)), and thus the flux tube starts with a faster rising speed (see, e.g., Murray et al. 2006). At the same time, the kinking itself accelerates the flux tube: when the tube kinks, its axis is stretched, which enhances the buoyancy and increases the rise speed (Fan et al. 1999).

The second case, spot–satellite, was modeled by introducing a parasitic flux tube set in a direction perpendicular to the main flux tube. Perhaps this type can also be produced from a single flux tube that bifurcates. However, for simplicity, we here tested the two-tube scenario (main and parasitic tubes). The parasitic tube has parameters of ${R}_{\mathrm{tube}}/{H}_{0}=2$, ${B}_{\mathrm{tube}}/{B}_{0}=15$, and ${{qH}}_{0}=-0.2$ (directed to $+y$ with a left-handed twist). The tube center is located at $(x/{H}_{0},y/{H}_{0},z/{H}_{0})=(15,0,-14)$ and is kept in mechanical balance. In this model, a periodic boundary condition was applied in the y-direction.

The quadrupole flux tube has two buoyant sections along the axis and thus starts emergence at two locations. We changed the density perturbation of Equation (17) to

$$\begin{eqnarray}\delta {\rho }_{\mathrm{exc}} & = & \rho (z)\displaystyle \frac{\delta {p}_{\mathrm{exc}}}{p(z)}\left[(1+\epsilon )\exp \left(-\displaystyle \frac{{(x-{x}_{\mathrm{tube}1})}^{2}}{{\lambda }^{2}}\right)\right.\\ & & +\left.(1+\epsilon )\exp \left(-\displaystyle \frac{{(x-{x}_{\mathrm{tube}2})}^{2}}{{\lambda }^{2}}\right)-\epsilon \right],\end{eqnarray} \tag{ 18 }$$

where ${x}_{\mathrm{tube}1}/{H}_{0}=-3\lambda /{H}_{0}=-24$ and ${x}_{\mathrm{tube}2}/{H}_{0}=3\lambda \,/{H}_{0}=24$.

Finally, for the inter-AR case, we set two flux tubes in parallel in the convection zone. The two tubes have parameters of $({x}_{\mathrm{tube}1}/{H}_{0},{y}_{\mathrm{tube}1}/{H}_{0})=(3\lambda /{H}_{0},-3\lambda /{H}_{0})=(24,-24)$ and $({x}_{\mathrm{tube}2}/{H}_{0}\,,{y}_{\mathrm{tube}2}/{H}_{0})=(-3\lambda /{H}_{0},3\lambda /{H}_{0})=(-24,24)$.

In this work, for the purposes of comparing the simulations to observations, we refer to the emerged region as δ-spots if it is complex (qualitatively), compact (separation of the two polarities smaller than, say, 20H₀), and highly sheared (shear angle of the PIL ∼ 90°). Note that we take these values for just a threshold in the simulations and they are not actually measured from observations. In addition, although previous observations found that the δ-spots often show rotational motions and violate Hale's polarity rule (e.g., Kurokawa 1987; López Fuentes et al. 2000, 2003), these properties are not used as the definition here.

Figure 9 summarizes the detailed photospheric evolution of the spot–spot case. In this case, as the twisted flux tube emerges, a complex magnetic pattern is formed in the surface layer, and an elongated PIL, highlighted by the Y-axis is eventually created at the center between the sunspot pair P1 and N1. One prominent feature here is the counter-streaming shear flow along the PIL (see ${{\boldsymbol{V}}}_{{\rm{h}}}$ vector), with its orientation following the expansion of the magnetic arcades in the atmosphere. As a result of the shear flow, the horizontal magnetic field becomes highly inclined to the PIL direction (see ${{\boldsymbol{B}}}_{{\rm{h}}}$ vector) and the shear angle becomes almost 90° (see panel (k)). Here, the shear angle is measured from the direction of the potential field (the direction perpendicular to the PIL) and thus 90° is parallel to the PIL. The length of the highly sheared (∼90°) part along the PIL is ${L}_{\mathrm{PIL}}/{H}_{0}\sim 60$.

Zoom In Zoom Out Reset image size

For easy comparison of the PIL with other simulation cases and actual observations, we introduce the diameter of a reference sunspot. In the reference case in Section 2.2, the final photospheric magnetogram at $t/{\tau }_{0}=300$ shows a simple bipolar pair (Figure 3). We measure the area of this reference spot (region with $| {B}_{z}| /{B}_{0}\geqslant 0.5$) and define ${D}_{\mathrm{spot}}$ as the diameter of the circle with an area equivalent to the area of this spot. Then, we obtain ${D}_{\mathrm{spot}}=24.6$, which is used as a normalizing factor for the length scale in Figures 9(k) and (l). With this value, one can find that the length of the highly sheared PIL of the spot–spot case is as large as ${L}_{\mathrm{PIL}}/{D}_{\mathrm{spot}}=2.5$.

One may also note that the developed PIL has a strong horizontal field (see panel (i)). Moreover, this PIL reveals an alternating pattern of positive and negative polarities (see panels (h) and (l)). These features are highly reminiscent of the "magnetic channel" structure, which was introduced by Zirin & Wang (1993) as one of the key characteristics of the flare-producing PIL (Kubo et al. 2007; Wang et al. 2008). From the comparison with numerical simulations, Kusano et al. (2012) and Bamba et al. (2013) suggested the possibility that small-scale flux emergence at the magnetic channel in NOAA AR 10930 triggers the series of large flare eruptions.

In the spot–satellite case, the newly emerging field in close proximity to the main sunspot of the opposite polarity creates a compact sheared PIL within a δ-spot. Figure 10 shows that the small bipole P2-N2 appears immediately right of ($+x$-side of) the main spot N1. The horizontal flow field and the relative motion (middle column and panel (j)) indicate that as N1 proceeds to the right, P2 drifts along the lower edge of ($-y$-side of) N1 and produces a sheared PIL. Reflecting the scale of the parasitic tube and thus of the satellite spot, the length of the highly sheared part of the PIL at $t/{\tau }_{0}=250$ is only ${L}_{\mathrm{PIL}}/{H}_{0}\sim 5$ or about 20% of the typical spot diameter, ${D}_{\mathrm{spot}}$. Furthermore, in this case, the horizontal field is stronger at the PIL (see panel (f)).

Zoom In Zoom Out Reset image size

The newly emerging fields at the edge of preexisting sunspots have been reported to occasionally drive major flares. For example, Louis et al. (2014) found that emerging satellite spots ahead of the leading sunspot of NOAA AR 11515 produce a filament at the PIL, which eventually erupts at the onset of the M5.6 class flare that develops into a CME: compare especially their Figure 7 and Figure 10 of this paper. Similar behaviors have been reported by, e.g., Wang et al. (1991), Ishii et al. (1998), Schmieder et al. (1994), and Takasaki et al. (2004), while simulations have shown that CME eruptions can be triggered by newly emerged flux at the edges of ARs (e.g., Chen & Shibata 2000).

In this case, as the satellite polarities (P2 and N2) move away from the main spot (N1), the δ-configuration and sheared PIL eventually disappear (see $t/{\tau }_{0}=280$ in Figure 10). The δ-spots are only seen in the earliest phase of the satellite emergence.

Advected by horizontal flows (middle column of Figure 11), the two inner sunspots of the quadrupole case, N1 and P2, collide with each other at the center of the simulation domain. The distance between the two spots shows a monotonic decrease (panel (j)), and a strongly packed δ-spot is eventually created. The highly sheared PIL has a length of ${L}_{\mathrm{PIL}}/{H}_{0}\sim 15$, or ${L}_{\mathrm{PIL}}/{D}_{\mathrm{spot}}\sim 0.6$, with the strongest B_z gradient (see panel (l)). The horizontal field is best enhanced at the central PIL (panel (i)).

Zoom In Zoom Out Reset image size

Sun et al. (2012) showed that the quadrupole AR NOAA 11158, producing the first X-class event in Cycle 24, hosts a highly sheared PIL between the two colliding sunspots (Schrijver et al. 2011; Liu et al. 2012; Toriumi et al. 2013). They pointed out that the PIL has a strong horizontal field, which is in good agreement with the PIL simulated here: compare Figure 5 of Sun et al. (2012) and Figure 11 of this paper.

The final case, inter-AR, does not show the clear formation of a sheared PIL or a δ-spot (Figure 12). The two inner sunspots, N1 and P2, remain separated from each other and simply show a fly-by motion (see panel (j)). In the central region in the photosphere, the horizontal field exhibits a slight indication of a magnetic shear (panel (g): Y-axis). However, the shear angle and the B_z gradient are not significant (panels (k) and (l)).

Zoom In Zoom Out Reset image size

As mentioned in Section 1, the X1.2-class flare from NOAA ARs 11944 and 11943 produced a very fast CME. Möstl et al. (2015) pointed out the importance of AR magnetic structures in controlling the eruption of the CME. Although this flare event was not from the sheared PIL, it produced a fast CME that channeled through the open magnetic flux created between the two closed field systems, ARs 11944 and 11943.

In the four simulation cases of flare-productive ARs, we found that the quadrupole case produces the most highly sheared PIL with the largest B_z gradient in a well-developed δ-spot. In order to investigate the evolution of the sheared PIL, we take the quadrupole case as an example and plot the terms of the induction equation in Figure 13. Here we show the shear component of the photospheric horizontal field, ${B}_{Y}/{B}_{0}$, i.e., the magnetic field along the Y-axis in Figure 11(g), and each term of the induction equation,

$$\begin{eqnarray}\displaystyle \frac{\partial {B}_{Y}}{\partial t} & = & \mathop{\underbrace{({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){V}_{Y}}}\limits_{\mathrm{Stretching}}\mathop{\underbrace{-({\boldsymbol{V}}\cdot {\boldsymbol{\nabla }}){B}_{Y}}}\limits_{\mathrm{Advection}}\\ & & \mathop{\underbrace{-\left(\displaystyle \frac{\partial {V}_{X}}{\partial X}+\displaystyle \frac{\partial {V}_{Y}}{\partial Y}\right){B}_{Y}}}\limits_{\mathrm{Compression}\ (\mathrm{horizontal})}\mathop{\underbrace{-\displaystyle \frac{\partial {V}_{z}}{\partial z}{B}_{Y}}}\limits_{\mathrm{Compression}\ (\mathrm{vertical})}.\end{eqnarray} \tag{ 19 }$$

In Equation (19), we neglect the magnetic diffusion and divide the compression term, $-({\boldsymbol{\nabla }}\cdot {\boldsymbol{V}}){B}_{Y}$, in the horizontal and vertical components.

Zoom In Zoom Out Reset image size

Figure 13 shows that the shear field ${B}_{Y}/{B}_{0}$ appears at $t/{\tau }_{0}=230$, peaks around $t/{\tau }_{0}=275$, and then gradually decays. During this period, the advection term is initially dominant ($t/{\tau }_{0}\lesssim 260$), and as the advection weakens, the stretching term increases and becomes comparable to the advection term ($260\lesssim t/{\tau }_{0}\lesssim 280$). For most of the time, the two compression terms remain negative and do not contribute to the growth of the shear. In the final phase after B_Y attains its peak ($t/{\tau }_{0}\gtrsim 280$), the total value becomes negative and thus B_Y decreases. However, the horizontal compression turns positive and becomes the only term that works to sustain B_Y.

The above behavior can be explained in the following manner. In the quadrupole case, the sunspots of positive and negative polarities (N1 and P2 in Figure 11) approach the region center and produce a δ-like configuration. In the early phase, as the two spots come closer, they transport the horizontal field from both sides (see the horizontal field vector shown with red arrows in Figure 11). This effect enhances the advection term in the early phase (Figure 13). Then, after the two spots merge, they show a drifting motion (N1 to the right and P2 to the left: yellow arrows in Figure 11), which stretches the horizontal field along the PIL, leading to the enhancement of the stretching term in the later phase (Figure 13). The compression by the two approaching spots also becomes stronger. However, this is only true for the horizontal component (Figure 13). Since the emergence is a process of a nonlinear instability, the rising field drastically expands vertically and $\partial {V}_{z}/\partial z$ is positive (e.g., Shibata et al. 1989). The negative contribution of the vertical compression term in Figure 13 reflects this process.

Figure 14 summarizes the 3D magnetic field structures for the four simulation cases. For the spot–spot case, the green magnetic field lines, each connecting the main spot and the extended tail, approach the center of the domain from both sides and form an electric current sheet between them (indicated by an isosurface in the middle column: see the Appendix for the plotted values). As a result, the green field lines reconnect with each other and create the purple and yellow field lines. The newly created purple flux system is highly sheared and aligned almost parallel to the photospheric PIL. However, this purple flux is trapped by the overlying yellow flux that connects the two main sunspots.

Zoom In Zoom Out Reset image size

In the spot–satellite case, as the main flux tube (yellow) pushes the parasitic tube up (green), magnetic reconnection occurs between the two flux tubes, and the purple field lines are formed. One may find that in the subsurface layer, the purple field lines extend to both the main and parasitic flux tubes. The current sheet is developed between the two flux systems immediately above the sheared PIL in the photosphere. Reflecting the smaller scale of the PIL (Figure 10), the purple flux is compact and very low-lying compared to the other cases. In contrast to the spot–spot case, this purple flux is located at the edge of the AR and is not trapped by the overlying fields, i.e., it is exposed to the outer space.

The quadrupole and inter-AR cases somewhat resemble each other. As explained in Section 3, the two emerging flux systems of the quadrupole case (yellow and green) originate from the common single flux tube and thus are connected beneath the surface. Consequently, the two photospheric polarities show convergence motion and become tightly packed. Driven by this photospheric motion, magnetic reconnection between the two coronal loops occurs in the current layer with a sheet-like shape extending in parallel to the photospheric PIL. Eventually, the purple flux system is newly created, which short-circuits the two inner polarities.

Although the two coronal loops of the inter-AR case are not originally connected beneath the surface and consequently, the contact of the two flux systems is less vigorous, they in fact undergo magnetic reconnection in the atmosphere since the two bipoles expand above the surface. A vertically extending current sheet is seen in between. The two bipoles eventually form purple field lines that connect the two independent ARs, which may be related to the flux rope that erupted as a fast CME between NOAA ARs 11944 and 11943 (Möstl et al. 2015).

In order to examine the accumulation of magnetic energy in the atmosphere, we calculate the potential magnetic fields from the B_z map at ${z}_{{\rm{p}}}/{H}_{0}=2$ and measure the total magnetic energy

$$\begin{eqnarray}&&{E}_{\mathrm{mag}}={\int }_{z\geqslant {z}_{{\rm{p}}}}\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi }\,{dV},\end{eqnarray} \tag{ 20 }$$

potential energy

$$\begin{eqnarray}&&{E}_{\mathrm{pot}}={\int }_{z\geqslant {z}_{{\rm{p}}}}\displaystyle \frac{{{\boldsymbol{B}}}_{\mathrm{pot}}^{2}}{8\pi }\,{dV},\end{eqnarray} \tag{ 21 }$$

and free energy

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{mag}}\equiv {E}_{\mathrm{mag}}-{E}_{\mathrm{pot}}.\end{eqnarray} \tag{ 22 }$$

Figure 15 compares the time evolutions of ${E}_{\mathrm{mag}}$, ${E}_{\mathrm{pot}}$, and ${\rm{\Delta }}{E}_{\mathrm{mag}}$ for the four cases. The time ${\rm{\Delta }}t/{\tau }_{0}$ is measured since the flux appears at the photosphere. Here, the simulation case with the highest energy is the spot–spot case. Reflecting the large photospheric flux (Figure 8), it has a total energy and free energy that are about one order of magnitude greater than those of the other three cases.

Zoom In Zoom Out Reset image size

The free energy of the spot–satellite case is higher than that of the reference case because free energy is stored in the current layer between the main and parasitic tubes in the spot–satellite case (see Figure 14, spot–satellite). Since the current sheet lies lower in the atmosphere, where the density is higher and the reconnection is less effective (Section 2.1), the free energy is not significantly consumed and gradually increases over time.

The free energies of the remaining two cases exhibit an interesting oscillatory behavior. For example, the quadrupole case shows a large bump around ${\rm{\Delta }}t/{\tau }_{0}=50$, which corresponds to $t/{\tau }_{0}=270$. This is because while the potential energy (${E}_{\mathrm{pot}}$) follows the monotonous growth of the photospheric flux (Figure 8), coronal reconnection between the two magnetic loops occurs when the two inner polarities approach from $t/{\tau }_{0}=240$ (Figure 6) and the actual magnetic energy (${E}_{\mathrm{mag}}$) starts to reduce, leading to the drastic loss of free energy. The free energy decrease of the inter-AR case after ${\rm{\Delta }}t/{\tau }_{0}=40$ (corresponding to $t/{\tau }_{0}=260$) is also due to the coronal reconnection of the two magnetic systems, which occurs later than in the quadrupole case because the two systems are significantly separated from each other.

The bottom panel of Figure 15 is intended to compare the four cases under the condition that they have similar AR scales. For each simulation case, we normalize the free energy by the three-halves power of its total unsigned flux at the final stage ($t/{\tau }_{0}=300$), ${{\rm{\Phi }}}_{\mathrm{final}}^{3/2}$. Note that the AR area is approximately proportional to the total unsigned flux, Φ, because the photospheric field is almost uniquely determined by the pressure balance between the magnetic field and the external gas. Considering that the AR volume is roughly proportional to the area to the three-halves, here we normalize the free magnetic energies of the four different cases with ${{\rm{\Phi }}}^{3/2}$.

One can find that since the difference among the four cases in the bottom panel is less prominent than that in the middle panel, the free energy depends highly on the photospheric total flux of each AR (or equivalently area or volume). However, as seen from the bottom panel, even for ARs of the similar size scales, the stored free energy may differ by much, up to a factor of five, depending on the twist and geometrical configuration of subsurface emerging fields.

The prediction of flares and CMEs is currently one of the most important topics of solar-terrestrial physics. Since the measurement of the photospheric parameters from vector magnetic data is much easier than the reconstruction of full 3D magnetic fields, most of the current flare prediction schemes are based on such photospheric parameters. After Leka & Barnes (2003) and Barnes et al. (2007) made use of the vector magnetogram for flare prediction, Bobra & Couvidat (2015) extracted various parameters (including those suggested by Leka & Barnes 2003; Schrijver 2007, and Fisher et al. 2012) from the SDO/HMI vector magnetogram for each AR (SHARP data: Bobra et al. 2014) and obtained a good predictive performance for flares of ≥M1.0-class using a machine-learning algorithm. By adding flare history and ultraviolet observables to the SHARP parameters, Nishizuka et al. (2017) further developed flare prediction models with even higher performance (see also Muranushi et al. 2015; Liu et al. 2017). The SHARP parameters used by these authors are summarized in Table 1. The F-score (Fisher score) in the table indicates the scoring of the parameter given by Bobra & Couvidat (2015).

Table 1.  Properties of Flare Events

Keyword	Description	Formula	F-score	CC	Scaling	CME Rank
totusjh	Total unsigned current helicity	${H}_{{c}_{\mathrm{total}}}\propto \sum | {B}_{z}\cdot {J}_{z}|$	3560	0.922	E	16
totbsq	Total magnitude of the Lorentz force	$F\propto \sum {B}^{2}$	3051	0.925	E	⋯
totpot	Total photospheric magnetic free energy density	${\rho }_{\mathrm{tot}}\propto \sum {({{\boldsymbol{B}}}^{\mathrm{Obs}}-{{\boldsymbol{B}}}^{\mathrm{Pot}})}^{2}{dA}$	2996	0.952	E	8
totusjz	Total unsigned vertical current	${J}_{{z}_{\mathrm{total}}}=\sum | {J}_{z}| {dA}$	2733	0.933	E	12
absnjzh	Absolute value of the net current helicity	${H}_{{c}_{\mathrm{abs}}}\propto | \sum {B}_{z}\cdot {J}_{z}|$	2618	0.833	E	13
savncpp	Sum of the modulus of the net current per polarity	${J}_{{z}_{\mathrm{sum}}}\propto | {\sum }^{{B}_{z}^{+}}{J}_{z}{dA}| +| {\sum }^{{B}_{z}^{-}}{J}_{z}{dA}|$	2448	0.781	E	18
usflux	Total unsigned flux	${\rm{\Phi }}=\sum | {B}_{z}| {dA}$	2437	0.894	E	10
area_acr	Area of strong field pixels in the active region	$\mathrm{Area}=\sum \mathrm{Pixels}$	2047	0.865	E	14
totfz	Sum of the z-component of the Lorentz force	${F}_{z}\propto \sum ({B}_{x}^{2}+{B}_{y}^{2}-{B}_{z}^{2}){dA}$	1371	0.745	E	⋯
meanpot	Mean photospheric magnetic free energy	$\overline{\rho }\propto \tfrac{1}{N}\sum {({{\boldsymbol{B}}}^{\mathrm{Obs}}-{{\boldsymbol{B}}}^{\mathrm{Pot}})}^{2}$	1064	−0.406	I	5
r_value	Sum of the flux near the polarity inversion line	${\rm{\Phi }}=\sum | {B}_{{LoS}}| {dA}$ within R mask	1057	0.855	E	15
epsz	Sum of the z-component of the normalized Lorentz force	$\delta {F}_{z}\propto \tfrac{\sum ({B}_{x}^{2}+{B}_{y}^{2}-{B}_{z}^{2})}{\sum {B}^{2}}$	864.1	−0.774	I	⋯
shrgt45	Fraction of area with shear $\gt 45^\circ$	Area with shear $\gt 45^\circ$/total area	740.8	0.725	I	7
meanshr	Mean shear angle	$\overline{{\rm{\Gamma }}}=\tfrac{1}{N}\sum \arccos \left(\tfrac{{{\boldsymbol{B}}}^{\mathrm{Obs}}\cdot {{\boldsymbol{B}}}^{\mathrm{Pot}}}{| {B}^{\mathrm{Obs}}| | {B}^{\mathrm{Pot}}| }\right)$	727.9	0.813	I	6
meangam	Mean angle of field from radial	$\overline{\gamma }=\tfrac{1}{N}\sum \arctan \left(\tfrac{{B}_{h}}{{B}_{z}}\right)$	573.3	−0.535	I	11
meangbt	Mean gradient of the total field	$\overline{| {\rm{\nabla }}{B}_{\mathrm{tot}}| }=\tfrac{1}{N}\sum \sqrt{{\left(\tfrac{\partial B}{\partial x}\right)}^{2}+{\left(\tfrac{\partial B}{\partial y}\right)}^{2}}$	192.3	−0.530	I	4
meangbz	Mean gradient of the vertical field	$\overline{| {\rm{\nabla }}{B}_{z}| }=\tfrac{1}{N}\sum \sqrt{{\left(\tfrac{\partial {B}_{z}}{\partial x}\right)}^{2}+{\left(\tfrac{\partial {B}_{z}}{\partial y}\right)}^{2}}$	88.40	−0.197	I	19
meangbh	Mean gradient of the horizontal field	$\overline{| {\rm{\nabla }}{B}_{h}| }=\tfrac{1}{N}\sum \sqrt{{\left(\tfrac{\partial {B}_{h}}{\partial x}\right)}^{2}+{\left(\tfrac{\partial {B}_{h}}{\partial y}\right)}^{2}}$	79.40	−0.474	I	1
meanjzh	Mean current helicity (Bz contribution)	$\overline{{H}_{c}}\propto \tfrac{1}{N}\sum {B}_{z}\cdot {J}_{z}$	46.73	−0.094	I	2
totfy	Sum of the y-component of the Lorentz force	${F}_{y}\propto \sum {B}_{y}{B}_{z}{dA}$	28.92	0.456	E	⋯
meanjzd	Mean vertical current density	$\overline{{J}_{z}}\propto \tfrac{1}{N}\sum \left(\tfrac{\partial {B}_{y}}{\partial x}-\tfrac{\partial {B}_{x}}{\partial y}\right)$	17.44	−0.221	I	9
meanalp	Mean characteristic twist parameter, α	${\alpha }_{\mathrm{total}}\propto \tfrac{\sum {J}_{z}\cdot {B}_{z}}{\sum {B}_{z}^{2}}$	10.41	−0.161	I	3
totfx	Sum of the x-component of the Lorentz force	${F}_{x}\propto -\sum {B}_{x}{B}_{z}{dA}$	6.147	0.352	E	⋯
epsy	Sum the of y-component of the normalized Lorentz force	$\delta {F}_{y}\propto \tfrac{-\sum {B}_{y}{B}_{z}}{\sum {B}^{2}}$	0.647	−0.018	I	⋯
epsx	Sum of the x-component of the normalized Lorentz force	$\delta {F}_{x}\propto \tfrac{\sum {B}_{x}{B}_{z}}{\sum {B}^{2}}$	0.366	−0.108	I	⋯

Note. For descriptions and formulae of the SHARP parameters, we follow the original notations of Bobra & Couvidat (2015). In the analysis of the simulation results of this paper, we measured each parameter at the ${z}_{{\rm{p}}}/{H}_{0}=2$ plane every ${\rm{\Delta }}t/{\tau }_{0}=2$ after magnetic flux appears at ${z}_{{\rm{p}}}/{H}_{0}=2$. The threshold for the absolute field strength above which the total and mean values are calculated is selected to be $B/{B}_{0}=0.04$ (equivalently 10 G), while the threshold for the "strong field" (area_acr) and for measuring Schrijver (2007)'s R (r_value) is $B/{B}_{0}=0.2$ (500 G). The F-score here is the scoring of the parameters for predicting solar flares provided by Bobra & Couvidat (2015), while CC is the correlation coefficient obtained from the four simulation cases by comparing the free magnetic energy and SHARP parameter at each moment (see top panels of Figure 16). Following Welsch et al. (2009), we classified the parameters into extensive (E), where a given parameter increases with AR size, and intensive (I), where the parameter is independent of AR size. The rightmost column shows the ranking of features for predicting CME eruptions, as reported by Bobra & Ilonidis (2016).

Download table as:  ASCIITypeset image

Then, the question is as follows. Why do some of these parameters predict the flare eruptions well (indicated by higher F-scores), while the others do not (lower F)? The series of numerical simulations of the present work, which successfully reproduced a variety of complex, non-potential configurations of flare-productive ARs, is one of the best ways to resolve this mystery. It is worth noting that Guennou et al. (2017) recently examined various photospheric parameters using simulations proposed by Leake et al. (2013, 2014) and found that parameters related to the PIL are the best to describe the flare occurrence. However, their analysis was restricted by the limited size of ARs and complexity of the simulations, which could be compensated for by our simulations.

Considering that the flare occurrence is a releasing process of free magnetic energy stored in the atmosphere, we compare the SHARP parameters in Table 1 and the stored free energy, i.e., the maximum flare energy that could potentially be released, for the four simulation cases. The top six panels of Figure 16 show samples of the comparisons. In each diagram, the horizontal axis represents a SHARP parameter in question that is measured at the surface (${z}_{{\rm{p}}}/{H}_{0}=2$) every ${\rm{\Delta }}t/{\tau }_{0}=2$ after the flux appears at the surface, whereas the vertical axis represents the stored free magnetic energy at each moment (${\rm{\Delta }}{E}_{\mathrm{mag}}/{E}_{0}$, measured directly from the 3D computational domain: see Section 5.2). One can find that some SHARP parameters have strong proportionalities with the free energy (e.g., totusjh and totpod), while others do not (e.g., epsy and epsx). It is reasonable that parameters such as totusjh (total unsigned current helicity) show high correlations because the free energy is stored in the form of electric current.

Zoom In Zoom Out Reset image size

For each diagram, we compute the correlation coefficient, CC, in a log–log plot, which indicates how accurately a given SHARP parameter reflects the stored free energy, and show it in the bottom right of the diagram. For each plot, we assume each data point to be independent and simply derive CC from all data points, regardless of the simulation cases. Therefore, there is only one CC for each plot. The CC values for all 25 parameters are summarized in Table 1.

The bottom panel of Figure 16 is a scatter plot of the absolute correlation, $| {CC}|$, versus F-score for all SHARP parameters. It is clearly seen that parameters with larger F yield larger $| {CC}|$. In other words, the SHARP parameters that are excellent in predicting the flare events can predict the free energy in the atmosphere very accurately. On the other hand, the smaller-F parameters have a weaker to almost no correlation with the free energy, indicating that they are incapable of predicting the free energy. It should be noted that in each scatter plot, the correlation is strong (weak) because all four simulations show consistently strong (weak) correlations. For example, $| {CC}|$ of totpot is 0.95, which is due to high correlations of the four cases: 0.89 (spot–spot), 0.85 (spot–satellite), 0.82 (quadrupole), and 0.98 (inter-AR).

The above relationship between $| {CC}|$ and F confirms the suggestion by Welsch et al. (2009) that parameters strongly associated with the flare activity are extensive (scaling with AR size: indicated by "E" in Table 1) because the free energy is likely stored on large scales and non-local.

The even higher prediction rates obtained by Nishizuka et al. (2017) have probably been obtained because they not only added the flare history, but also included ultraviolet observables that are sensitive to chromospheric dynamics such as triggering processes (preflare brightenings) before the flares occur (see, e.g., Bamba et al. 2013 for a detailed observational analysis of the flare triggers). This suggests that although the SHARP parameters properly indicate the accumulation of free energy, this is not sufficient to accurately predict the exact occurrence of flares, and we need additional observables that represent the triggering of the flares.

The rightmost column of Table 1 shows the ranking of SHARP parameters for predicting CME eruptions, as reported by Bobra & Ilonidis (2016). In contrast to the flare prediction case, the parameters that successfully predict whether a given flare produces a CME are mostly intensive (independent of AR size: indicated by "I" in Table 1), not extensive. One can also see from this table that they show moderate to weak negative correlations (${CC}\lt 0$) between the SHARP parameters and the stored free energy. This is because these parameters are, in most cases, normalized by the AR size or some relevant factors. That is, the CME-predictive parameters are not perfectly independent of the AR size; rather, they have a negative dependence on the AR scale.

For example, according to Bobra & Ilonidis (2016), meangbt is one of the highest-ranking CME parameters, and it is the sum of the horizontal magnetic gradient normalized by the SHARP patch pixels N (representing AR area):

$$\begin{eqnarray}&&\overline{| {\rm{\nabla }}{B}_{\mathrm{tot}}| }=\displaystyle \frac{\sum \sqrt{{\left(\tfrac{\partial B}{\partial x}\right)}^{2}+{\left(\tfrac{\partial B}{\partial y}\right)}^{2}}}{N}.\end{eqnarray} \tag{ 23 }$$

Since the flare-causing PILs tend to have a high magnetic gradient (Section 1), the numerator of Equation (23) becomes larger for flaring ARs. However, the normalization by area cancels this trend, leading to a negative correlation with the free energy (${CC}=-0.530$: Figure 16) and thus a worse flare prediction rate (F = 192.3). Still, this parameter yields a high CME prediction performance. This may also be due to the normalization effect. That is, while the local flaring zone is characterized by the field gradient (numerator), the global scale of the AR is represented by the AR area (denominator), and the relationship between the two factors (their ratio) determines the CME productivity.

The above discussion is further supported by observational studies. Sun et al. (2015) investigated the flare-rich but CME-poor AR NOAA 12192 and found that this AR has a relatively weak non-potentiality with a relatively strong overlying field, and thus, the flux ropes fail to erupt as CMEs. They concluded that CME eruption is described by the relative measure of non-potentiality over the restriction of the background field. A statistical analysis by Toriumi et al. (2017) revealed that CME-less ARs show, on average, smaller ${S}_{\mathrm{ribbon}}/{S}_{\mathrm{spot}}$ (flare ribbon area normalized by total sunspot area) and $| {\rm{\Phi }}{| }_{\mathrm{ribbon}}/| {\rm{\Phi }}{| }_{\mathrm{AR}}$ (total unsigned magnetic flux in the ribbon normalized by the total unsigned flux of the entire AR). They argued that the magnetic relation between the large-scale structure of an AR and the localized flaring domain within it is a key factor determining the CME eruption. We shall leave the detailed numerical investigation of this relation for future work.