Automated Driving for Global Nonpotential Simulations of the Solar Corona

To illustrate the electric field computation we will use an active region from the automatically generated HMI/SHARP database (Bobra et al. 2014), shown in Figure 1. A single line-of-sight magnetogram was selected when the region was close to Central Meridian. Following the methodology of Yeates (2020a), this has been (i) downsampled to the 180 × 360 grid in (s, ϕ) used throughout this paper, and (ii) corrected for flux balance. As discussed in the introduction, we use only line-of-sight magnetograms and hence choose a magnetogram near Central Meridian so as to reduce projection effects. This does have the consequence that any flux emergence after the region has passed Central Meridian will be missed, but this is a limitation of the input data rather than the emergence algorithm described in this section.

Zoom In Zoom Out Reset image size

Since our data comprise a single magnetogram for each active region, we assume that every region emerges over a fixed time interval T_em, typically chosen as 24 hr, although the effect of this parameter will be considered in Section 5.1. During this interval we impose a steady electric field E _⊥(s, ϕ) satisfying

$$\begin{eqnarray}&&-\hat{{\boldsymbol{r}}}\cdot {\rm{\nabla }}\times {{\boldsymbol{E}}}_{\perp }(s,\phi )=\displaystyle \frac{{B}_{r}(s,\phi )}{{T}_{\mathrm{em}}},\end{eqnarray} \tag{ 1 }$$

so that the required B_r (s, ϕ) distribution is superimposed at the end of this emergence on any pre-existing B_r distribution (which is usually weak in our simulations). The assumption of a steady emergence rate is reasonable for global simulations that focus on the postemergence evolution over longer timescales. In fact, piecewise-steady electric fields were also used to interpolate between subsequent magnetograms in the active region simulations of Mackay et al. (2011) and Cheung & DeRosa (2012).

Although any E _⊥ satisfying (1) will generate the same B_r on the surface, different choices could lead to different coronal magnetic fields when used as boundary conditions for three-dimensional simulations. Our approach is to decompose E _⊥ into two contributions:

$$\begin{eqnarray}&&{{\boldsymbol{E}}}_{\perp }(s,\phi )={{\boldsymbol{E}}}_{\perp }^{0}(s,\phi )-{\rm{\nabla }}\displaystyle \frac{{\rm{\Psi }}(s,\phi )}{{T}_{\mathrm{em}}}.\end{eqnarray} \tag{ 2 }$$

The first term ${{\boldsymbol{E}}}_{\perp }^{0}$ is the inductive component of the emergence electric field, and as such will generate a coronal magnetic field of minimum possible complexity. Our method for computing this term is described in Section 2.3. The other term is the noninductive component that does not change the corresponding magnetogram B_r but allows us to vary the free energy of the emerging regions. This twisting potential Ψ is discussed in Section 2.4.

From the original local region D—which may consist of several disconnected pieces—we first identify a simply connected region, $\overline{D}$, containing the newly emerging flux. Our aim is to compute an electric field that vanishes outside this region. As an illustration, consider the region shown in Figure 1. The black shading in Figure 1(b) shows the original region, D (chosen as ∣B_r ∣ > 10⁻³ G). This is not simply connected, having several disconnected components containing holes. However, our algorithm for calculating ${{\boldsymbol{E}}}_{\perp }^{0}$ requires a single simply connected region. From D, we apply successive steps of Gaussian smoothing to ∣B_r ∣ until such a region $\overline{D}$ is reached, shown by the wider gray area in Figure 1(b).

The electric field is computed only on $\overline{D}$, with the boundary condition E _⊥ × n = 0 on the boundary $\partial \overline{D}$, consistent with magnetic flux balance in $\overline{D}$. A staggered grid is used where ${B}_{r}^{i,j}$ is defined at the cell centers, ${E}_{s}^{i,j+1/2}$ is defined on the s-edges and ${E}_{\phi }^{i+1/2,j}$ is defined on the ϕ-edges (e.g., Figure 2). Here i and j denote the latitudinal (s) and azimuthal (ϕ) indices of the cell. Thus with Stokes' Theorem applied to Equation (1), ${E}_{s}^{i,j+1/2}$ and ${E}_{\phi }^{i+1/2,j}$ must satisfy

$$\begin{eqnarray}&&\begin{array}{l}{\left({{\ell }}_{s}{E}_{s}\right)}^{i,j-1/2}-{\left({{\ell }}_{s}{E}_{s}\right)}^{i,j+1/2}+{\left({{\ell }}_{\phi }{E}_{\phi }\right)}^{i+1/2,j}\\ -{\left({{\ell }}_{\phi }{E}_{\phi }\right)}^{i-1/2,j}=\displaystyle \frac{{R}_{\odot }^{2}{\rm{\Delta }}s{\rm{\Delta }}\phi {B}_{r}^{i,j}}{{T}_{\mathrm{em}}},\end{array}\end{eqnarray} \tag{ 3 }$$

where ${R}_{\odot }^{2}{\rm{\Delta }}s{\rm{\Delta }}\phi$ is the area of a grid cell—uniform in these coordinates. The cell edge lengths are given by ${{\ell }}_{s}^{i,j+1/2}={R}_{\odot }[\arcsin ({s}^{i+1/2})-\arcsin ({s}^{i-1/2})]$ and ${{\ell }}_{\phi }^{i+1/2,j}={R}_{\odot }{\rm{\Delta }}\phi \sqrt{1-{\left({s}^{i+1/2}\right)}^{2}}$. The boundary condition is enforced by setting tangential components of ${E}_{s}^{i,j+1/2}$ or ${E}_{\phi }^{i+1/2,j}$ to zero on the boundary edges of $\overline{D}$, shown by dashed lines in Figure 2.

Zoom In Zoom Out Reset image size

When applied to each grid cell in $\overline{D}$, Equation (3) gives a system of linear equations for ${E}_{s}^{i,j+1/2}$ and ${E}_{\phi }^{i+1/2,j}$. The system is underdetermined, reflecting the freedom to add an arbitrary gradient term to E _⊥ without affecting $\hat{{\boldsymbol{r}}}\cdot {\rm{\nabla }}\times {{\boldsymbol{E}}}_{\perp }$. Our implementation solves for ${\left({{\ell }}_{s}{E}_{s}\right)}^{i,j+1/2}$ and ${\left({{\ell }}_{\phi }{E}_{\phi }\right)}^{i+1/2,j}$, since these weighted values are most useful in the magnetofrictional code used for our tests (Section 5). We choose the least-squares solution, ${{\boldsymbol{E}}}_{\perp }^{0}$, to this system, i.e., the solution that minimizes the L₂-norm. This is implemented with a sparse matrix and a sparse least-squares solver (in Python). This is mathematically equivalent to finding an electric field of the form ${{\boldsymbol{E}}}_{\perp }^{0}=-{\rm{\nabla }}\times \left({\rm{\Phi }}\hat{{\boldsymbol{r}}}\right)$, which is an inductive electric field (Yeates 2017). However, unlike the inductive electric fields computed by Weinzierl et al. (2016) or Yeates (2017), who solved Equation (3) over the full solar surface, here we have allowed nonzero electric field only within the region $\overline{D}$. Hence, we call this a "local inductive" electric field solution. The resulting E_s ⁰ and ${E}_{\phi }^{0}$ are shown in Figures 3(a) and (e).

Zoom In Zoom Out Reset image size

Although this solution ${{\boldsymbol{E}}}_{\perp }^{0}$ satisfies E _⊥ × n = 0 on $\partial \overline{D}$, Figures 3(a) and (e) show that the normal component ${{\boldsymbol{E}}}_{\perp }^{0}\cdot {\boldsymbol{n}}$ can be discontinuous at the edge of $\overline{D}$. While this does not affect $\hat{{\boldsymbol{r}}}\cdot {\rm{\nabla }}\times {{\boldsymbol{E}}}_{\perp }$, it could lead to the generation of spurious currents in coronal simulations. We remove the sharp discontinuity by smoothing ${{\boldsymbol{E}}}_{\perp }^{0}$ with an evolution of the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{E}}}_{\perp }^{0}}{\partial t}={\rm{\nabla }}\left({\rm{\nabla }}\cdot {{\boldsymbol{E}}}_{\perp }^{0}\right).\end{eqnarray} \tag{ 4 }$$

This leaves $\hat{{\boldsymbol{r}}}\cdot {\rm{\nabla }}\times {{\boldsymbol{E}}}_{\perp }^{0}$ unaffected provided the gradient is discretized consistently, which we ensure by using the staggered form

$$\begin{eqnarray}\begin{array}{rcl}{\left({{\ell }}_{s}{E}_{s}\right)}_{k+1}^{i,j+1/2} & = & {\left({{\ell }}_{s}{E}_{s}\right)}_{k}^{i,j+1/2}\\ & & +{\rm{\Delta }}t\left({D}_{k}^{i+1/2,j+1/2}-{D}_{k}^{i-1/2,j+1/2}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\left({{\ell }}_{\phi }{E}_{\phi }\right)}_{k+1}^{i+1/2,j} & = & {\left({{\ell }}_{\phi }{E}_{\phi }\right)}_{k}^{i+1/2,j}\\ & & +{\rm{\Delta }}t\left({D}_{k}^{i+1/2,j+1/2}-{D}_{k}^{i+1/2,j-1/2}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

where ${D}_{k}^{i+1/2,j+1/2}$ is a finite-difference approximation to ${\rm{\nabla }}\cdot {{\boldsymbol{E}}}_{\perp }^{0}$ at time step k. The step size Δt is chosen for numerical stability and we integrate for a total time ≈7 × 10⁻⁵ T, where T is the time that would be taken to diffuse a distance R_⊙. We choose this short integration time to minimize the amount of smoothing, which is a balance between reducing currents at the edge of $\overline{D}$ and maintaining localization of ${{\boldsymbol{E}}}_{\perp }^{0}$. If the amount of smoothing were increased, ${{\boldsymbol{E}}}_{\perp }^{0}$ would eventually tend toward the global inductive solution. The result with our chosen amount of smoothing is illustrated in Figures 3(b) and (f). We denote this smoothed local inductive electric field by $\overline{{{\boldsymbol{E}}}_{\perp }^{0}}$, and use it in practice for the ${{\boldsymbol{E}}}_{\perp }^{0}$ component in Equation (2).

It is worth noting that $\overline{{{\boldsymbol{E}}}_{\perp }^{0}}$ differs significantly from electric fields that have previously been used to drive magnetofrictional simulations. For example, Figures 3(c) and (g) show the components of the global inductive solution, which satisfies ${{\boldsymbol{E}}}_{\perp }=-{\rm{\nabla }}\times \left({\rm{\Phi }}\hat{{\boldsymbol{r}}}\right)$ on the whole spherical surface. By definition, this minimizes the L₂-norm of E _⊥ (70 V m⁻¹ compared to 111 V m⁻¹ for the local inductive solution), but at the expense of spreading nonzero E _⊥ far outside of $\overline{D}$, decaying only as the inverse of the distance from $\overline{D}$ (Yeates 2017). Although no additional B_r (R_⊙, s, ϕ) is generated by this nonzero E _⊥ outside $\overline{D}$, additional electric fields are generated and these can change the evolution of the coronal magnetic field outside of the emergence region (illustrated in Section 5.3 below). Thus our local inductive $\overline{{{\boldsymbol{E}}}_{\perp }^{0}}$ is more appropriate for incorporating newly emerging regions in global simulations.

As a further comparison, Figures 3(d) and (h) show a sparse E _⊥ computed over the whole solar surface for this region, using the L₁-minimization method of Yeates (2017). Because it is essentially a compromise between sparsity (localization) and smoothness, the E _⊥ generated by the sparse algorithm is less smooth than our L₂-based electric fields. For a fairer comparison, the solution shown in Figures 3(d) and (h) has had the same smoothing applied via Equation (4). It is clear that the global sparse E _⊥ shares some of the broad features of the local inductive solution, with ∣E_s ∣ > ∣E_ϕ ∣ overall and localization around $\overline{D}$. However, there are two disadvantages of the sparse solution. First, although it is localized, it is not necessarily localized specifically to the region $\overline{D}$, as in this example. Second, the sparse solution is more expensive to compute because it requires multiple least-squares solutions over the whole solar surface, in an iterative loop, rather than a single solution over only the local region $\overline{D}$. Thus we find that the local inductive $\overline{{{\boldsymbol{E}}}_{\perp }^{0}}$ is an improvement. Note that, although it does not concern us here, the sparse solution is also found to lead to significant spurious coronal currents when applied with higher-cadence sequences of observational magnetograms (Mackay & Yeates 2021). Since minimizing the L₁ norm is equivalent to finding the sparsest solution, we can compare the relative sparseness of the different E _⊥ by their L₁ norms, which for this example are 2256 V m⁻¹ for the local inductive solution, 6115 V m⁻¹ for the global inductive solution, and 1978 V m⁻¹ for the global sparse solution. Thus—as is seen in Figure 3—the local inductive solution does not have a significantly larger footprint than the sparse solution itself, even after the smoothing step.

The noninductive electric field, resulting from −∇(Ψ/T_em) in Equation (2), typically increases both the Poynting flux of magnetic energy and the flux of magnetic helicity into the corona (Kazachenko et al. 2014; Pomoell et al. 2019). Indeed, the existence of nonzero free magnetic energy is an important feature of many real active regions. While we do not follow the detailed injection of this free energy during the emergence process, it is nevertheless important to account for the net emergence of any nonpotential structure. Thus we include a nonzero "twisting" component −∇(Ψ/T_em). In line with observations and the theory of helicity condensation, we choose to concentrate this twisting near polarity inversion lines, similar to Mikić et al. (2018; see also Mackay et al. 2018; Yeates et al. 2018). The choice of twisting potential Ψ(s, ϕ) will not affect $\hat{{\boldsymbol{r}}}\cdot {\rm{\nabla }}\times {{\boldsymbol{E}}}_{\perp }$, but it will lead to the generation of additional B_s and B_ϕ components in the corona since it is imposed only at r = R_⊙. Note that Mackay et al. (2018) instead modified E by adding an E_r component (in their case, through modifying A ), but the effect on B is the same, the difference being only in choice of noninductive component.

The potential Ψ(s, ϕ) for a given region is constructed as follows, illustrated in Figure 4:

• 1.  
  Polarity inversion lines are identified using the function

  $$\begin{eqnarray}&&{f}_{\mathrm{pil}}(s,\phi )={\left(\left|{\rm{\nabla }}{\overline{B}}_{r}\right|-\left|{\rm{\nabla }}\left|{\overline{B}}_{r}\right|\right|\right)}^{2},\end{eqnarray} \tag{ 7 }$$

  where ${\overline{B}}_{r}$ is a smoothed B_r distribution, computed by mild Gaussian smoothing with a standard deviation of less than one grid cell (Figure 4(a)). The function f_pil is nonzero near polarity inversion lines and zero away from them. It is strongest where the gradient of B_r is strongest (Figure 4(b)).
• 2.  
  A somewhat broader twisting region is defined by Gaussian smoothing of f_pil with a standard deviation of two grid cells to give 〈f_pil〉. This is chosen approximately so as to localize the twisting within the interior of $\overline{D}$ (Figure 4(c)).
• 3.  
  The potential is then

  $$\begin{eqnarray}&&{\rm{\Psi }}(s,\phi )=\tau {b}_{0}\langle {f}_{\mathrm{pil}}\rangle {\overline{B}}_{r},\end{eqnarray} \tag{ 8 }$$

  as shown in Figure 4(d). Here b₀ is a constant normalization factor, chosen so that $\max | {\rm{\nabla }}{\rm{\Psi }}| =\tau L\max | \overline{{B}_{r}}|$. The dimensionless "twist" parameter τ is signed such that τ > 0 leads to positive helicity (sinistral chirality) and τ < 0 to negative helicity (dextral chirality). We set L = 0.015R_⊙, which is the height of the first coronal grid cell in our magnetofrictional calculations. With this normalization, the horizontal field strength generated by the twisting should be roughly of the order ∣τ B_r ∣.

Zoom In Zoom Out Reset image size

The coupled surface flux transport and magnetofrictional approach were first introduced by van Ballegooijen et al. (2000), and later Yeates et al. (2008) extended it to cover the global corona. The basic idea is that the coronal magnetic field B evolves continuously through a sequence of nonpotential states in response to the changing surface boundary. Large-scale plasma motions and magnetic diffusion on the surface drive the coronal field evolution in a quasi-static manner. In the corona, we solve for the magnetic vector potential B = ∇ × A ,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{A}}}{\partial t}=-{\boldsymbol{E}},\end{eqnarray} \tag{ 9 }$$

where E = − v × B + N . Here, E and N represent the electric field, and the nonideal part of Ohm's law, respectively. Although the corona is highly conducting, the nonideal term reflects the fact that we are modeling the large-scale, mean magnetic field; thus N describes the effect of unresolved smaller-scale turbulent motions (van Ballegooijen et al. 2000). The velocity v in the corona—discussed in the next paragraph—is a combination of frictional relaxation and a solar wind outflow. On the lower boundary, we apply a different electric field which is a superposition of (i) the emergence electric fields E _⊥ for individual active regions, as per Section 2, and (ii) a global electric field ${{\boldsymbol{E}}}_{\perp }=-{{\boldsymbol{v}}}_{s}\times ({B}_{r}\hat{{\boldsymbol{r}}})+{\eta }_{0}{\rm{\nabla }}\times ({B}_{r}\hat{{\boldsymbol{r}}})$. Here the large-scale velocity v _s (θ) on the lower surface mimics two large-scale plasma flows: meridional circulation and differential rotation (for more details, see Section 2.2 of Yeates 2014). The surface diffusivity η₀ = 455 km² s⁻¹ represents the net effect of small-scale (supergranular) flows. The computational domain extends radially within R_⊙ ≤ r ≤ 2.5 R_⊙ and includes the full extent of colatitude 0°–180° and longitude 0°–360°. We solve Equation (9) using a finite-difference method on an equally spaced grid of 60 × 180 × 360 cells in $\mathrm{log}(r/{R}_{\odot })$, sine latitude, and longitude.

This nonpotential coronal model is a simplified version of full-scale MHD models. We do not solve the full momentum equation to simulate the velocity evolution, but rather employ a "frictional" velocity v proportional to the Lorentz force. Such an artificial velocity field drives the coronal magnetic field toward a force-free equilibrium, J × B = 0 where J = ∇ × B . The velocity field within the corona is modeled accordingly as

$$\begin{eqnarray}&&{\boldsymbol{v}}=\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{\nu }+{v}_{\mathrm{out}}(r)\hat{{\boldsymbol{r}}}\end{eqnarray} \tag{ 10 }$$

where $\nu ={\nu }_{0}\,| {\boldsymbol{B}}{| }^{2}/({r}^{2}{\sin }^{2}\theta )$, is the friction coefficient, with ν₀ = 2.8 × 10⁵ s. On the photosphere, the frictional velocity is set to zero. The second term in Equation (10), with ${v}_{\mathrm{out}}(r)={v}_{0}{\left(r/{R}_{\odot }\right)}^{11.5}$ and v₀ = 100 km s⁻¹, models (crudely) the effect of the solar wind in the upper corona (see Rice & Yeates 2021). It ensures that magnetic field lines become essentially radial by 2.5 R_⊙, except temporarily during eruptions. All components of B are periodic in ϕ. At the inner (1.0 R_⊙) and outer (2.5 R_⊙) boundaries the transverse current density, J × $\hat{{\boldsymbol{r}}}$, is set to zero, such that the Lorentz force is always tangential to the boundaries. The nonideal term N is modeled using fourth-order hyperdiffusion, which preserves magnetic helicity density A · B in the volume (van Ballegooijen & Cranmer 2008). We use, N = − ( B /∣ B ∣²) ∇ (η_h∣ B ∣²∇α), where α = J · B /∣ B ∣² and η_h = 10³¹ cm⁴ s⁻¹.

We perform three-dimensional magnetofriction simulations for 213 days, covering 2011 June 1 to 31 2011 December 31, around the Solar Cycle 24 Maximum. A three-dimensional magnetic field distribution is required to initialize the simulation, which we obtain from a potential field source surface extrapolation (Yeates 2018) from the HMI radial-component pole-filled Carrington map 2210 (Sun 2018).

The emerging regions during this period are extracted automatically from the HMI/SHARP database, using the methodology of Yeates (2020a), which is fully automated and repeatable (the code is freely available at Yeates 2022). In summary, a single line-of-sight magnetogram is selected at the time when the region was closest to Central Meridian. The extraction code discards those SHARPs that are either too small to resolve at the chosen 180 × 360 resolution, have highly imbalanced flux (near unipolar), or represent a repeat observation of an earlier region surviving on the solar surface for more than one Carrington rotation. This leaves 102 active regions with realistic shapes to emerge in the magnetofrictional simulations, using our new emergence technique described in Section 2.

Since our simulations model solar maximum with a hundred active regions emerging, the coronal magnetic field evolution is quite dynamic, with numerous nonpotential structures forming and erupting. Thus, unlike the solar minimum study by Bhowmik & Yeates (2021), where individual nonpotential structures were tracked and analyzed over a few days, here we primarily focus on global measures of nonpotentiality in the three-dimensional coronal magnetic field. The temporal evolution of these measures is depicted in Figure 5 for a particular magnetofrictional simulation where active regions have emergence duration T_em = 24 hr and all have the same twist τ = ±0.1 (positive in the southern hemisphere, negative in the northern hemisphere).

Zoom In Zoom Out Reset image size

Figure 5(a) shows the evolution of the photospheric magnetic flux calculated by integrating ∣B_r ∣ over the surface r = R_⊙. The sudden rises in photospheric flux denote the epochs of active region emergence. Besides the newly emerged active regions, any significant changes in the coronal magnetic field distribution, especially at the outer boundary (2.5 R_⊙), are caused by eruptions of coronal structures. These need not be directly caused by active region emergence, as seen in Figure 5(b), which shows the open magnetic flux (integral of ∣B_r ∣ over the surface r = 2.5 R_⊙). Figure 5(c) shows the total magnetic energy in the corona, calculated by integrating ∣ B ∣²/(8π) over the whole simulation volume. Another important measure is the mean current density ∣ J ∣ per unit volume within the corona, shown in Figure 5(d). The current density is directly linked with the nonpotential part of the magnetic field.

We also study the evolution of the average horizontal magnetic field, ${B}_{\perp }={\left({B}_{\theta }^{2}+{B}_{\phi }^{2}\right)}^{1/2}$, at r = 2.5 R_⊙. The coronal magnetic field at the outer boundary is primarily radial due to the solar wind. However, B_⊥ becomes significantly enhanced when any nonpotential coronal magnetic structure erupts and passes through the outer boundary (Yeates & Hornig 2016; Lowder & Yeates 2017; Bhowmik & Yeates 2021). The temporal evolution of the average 〈B_⊥〉 over the surface r = 2.5 R_⊙ is depicted in Figure 5(e), where we separate the averages in the northern and southern hemispheres.

The bottom row, Figure 5(f), shows the evolution of (relative) helicity flux through the outer boundary. Helicity is injected into the coronal magnetic field through large-scale shearing by differential rotation as well as the emergence of twisted active regions. With hyperdiffusion, the volume dissipation of helicity is negligible, with only an extremely small amount from numerical diffusion (e.g., Figure 1 of Bhowmik & Yeates 2021). Thus the ejection of unstable nonpotential structures with high helicity content through the outer boundary is the prime means by which the coronal model sheds helicity. The total (signed) flux of relative helicity through the outer boundary is gauge independent and computed by integrating $\tilde{{\boldsymbol{A}}}\times \left(2{\boldsymbol{E}}+\partial \tilde{{\boldsymbol{A}}}/\partial t\right)$ over the surface r = 2.5 R_⊙, where $\tilde{{\boldsymbol{A}}}$ is an appropriately chosen vector potential that must be computed from B (for details, see Yeates & Hornig 2016).

As discussed already, B_⊥ is a useful indicator for eruptions of nonpotential coronal structures, as well as the amount of magnetic flux involved. We calculate 〈B_⊥〉 at 2.5 R_⊙ at a cadence of 20 s; thus, individual peaks with significant amplitudes can be considered to be associated with separate eruptions. As an indicator of the number of eruptions, we count the number of peaks in 〈B_⊥〉 (for each hemisphere separately). The total number of peaks will depend on the particular threshold of B_⊥ imposed for selecting the peaks. Figures 6(a)–(c) show the selected peaks with asterisks for a reasonable choice of the threshold. Figure 6(d) depicts how the number of peaks (in the northern hemisphere) varies with an increasing threshold value, and for simulations with different active region twists τ (to be discussed in the following section). The simulation that is shown in Figures 5 and 6(a)–(c) is the dark blue curve in Figure 6(d) (labeled τ = 0.1).

Zoom In Zoom Out Reset image size

In the first set of simulations, we fix ∣τ∣ to the same magnitude for all emerging regions. As an illustration of the effect of this parameter, Figure 7 shows the magnetic field structure after the emergence of an individual region in the magnetofrictional code, for four different values τ = 0, −0.1, −0.4, 0.1. Each magnetic field line L is colored by its helicity contribution ${ \mathcal A }| {B}_{r}|$, where

$$\begin{eqnarray}&&{ \mathcal A }={\int }_{L}\displaystyle \frac{\tilde{{\boldsymbol{A}}}\cdot {\boldsymbol{B}}}{| {\boldsymbol{B}}| }\,{\rm{d}}l\end{eqnarray} \tag{ 11 }$$

is the field line helicity of the field line (Yeates & Hornig 2016), and ∣B_r ∣ is averaged between the two endpoints of L. We use the vector potential $\tilde{{\boldsymbol{A}}}$ (as in Section 3.3), in order that integrating ${ \mathcal A }| {B}_{r}|$ over both boundaries r = R_⊙ and r = 2.5R_⊙ would give (twice) the relative helicity (Yeates 2020b). The quantity ${ \mathcal A }| {B}_{r}|$ is a robust measure of the distribution of magnetic helicity within the corona, as both ${ \mathcal A }$ and ∣B_r ∣ would be invariant during an ideal evolution in which the field line endpoints remained fixed. From Figure 7, we see that the amount and sign of helicity generated within the core of the active region are clearly controlled by the τ parameter. Moreover, Figure 7(c) suggests that τ = −0.4 generates unphysically high twisting in the corona. We therefore consider four full test simulations with ∣τ∣ = 0.0, 0.05, 0.1, 0.2.

Zoom In Zoom Out Reset image size

Depending on the sign of the associated τ, an emerging region will have either positive or negative helicity. Observational works (e.g., Pevtsov et al. 1995; Wang 2013; Park et al. 2020) have found that helicity of active regions, on average, follows the hemispheric rule, such that active regions in the northern and southern hemispheres are more likely to have negative and positive helicity, respectively. This is known as the hemispheric sign rule and is observed in other coronal features, too (Pevtsov et al. 2014). Thus, for this simplest set of simulations, we impose that all emerging regions strictly follow the hemispheric rule. Note that this does not preclude the formation of structures of the opposite sign of helicity as active regions spread out and interact. For example, in simulations with simplified bipolar active regions, Yeates & Mackay (2009) looked at the chirality of filament channels as active regions decay and interact. This chirality is a proxy for the sign of magnetic helicity and is observed to follow an even stronger hemispheric rule than active region helicity (Martin 1998; Mackay et al. 2010). It was found that reversing the sign of helicity in the emerging regions caused only 32% of the simulated filament channels to reverse their chirality, highlighting how active region emergence is only one source of magnetic helicity in the corona.

As mentioned above, observations show clearly that there are deviations from the hemispheric helicity sign rule in active regions (Wang 2013; Park et al. 2020). To explore the effect of such variations, we try a second set of simulations where ∣τ∣ is still the same for all regions, but the sign of τ is chosen for each individual region according to vector magnetogram observations. These vector data are not used in determining ${{\boldsymbol{E}}}_{\perp }^{0}$ in Equation (2), but they can, in principle, give us physical information to constrain Ψ. In particular, the HMI/SHARP metadata provide the parameter ${\alpha }_{\mathrm{ob}}=\sum {J}_{z}\,{B}_{z}/(\sum {B}_{z}^{2}$) (in Mm⁻¹) for each region (Table 3 of Bobra et al. 2014). This does not correspond directly to the magnetic helicity of the region (which cannot be computed without further data) but is a proxy for the current helicity (integral of J · B ), which represents the net local twisting of the magnetic field. The current helicity often, though not always, has the same sign as the magnetic helicity itself (Russell et al. 2019). The distributions of α_ob in each hemisphere are depicted in Figure 8 for our test simulation period. We can see that among the 102 active regions during the simulation period, the majority follow the hemispheric sign rule (negative/positive in the North/South). However, we notice a significant number that deviate from the sign rule. Thus, in our second set of magnetofrictional simulations, we consider constant values of ∣τ∣, but with the sign chosen according to the observed sign of α_ob for that region.

Zoom In Zoom Out Reset image size

In a more realistic simulation, emerging regions would not all have the same amplitude of twist. From the distribution of observed α_ob (Figure 8), it is clear that active regions have diverse values of α_ob, suggesting that diverse values of τ are likely appropriate. This is also supported by computations of photospheric helicity flux in the literature (Georgoulis et al. 2009). As a first attempt to take this variation into account, we perform a third set of simulations where the τ for each region is proportional to its observed value of α_ob from the HMI/SHARPs metadata. We set τ = c α_ob, where the proportionality constant c has units of length.

To determine the best choice of c, we compare α_ob for each of the 102 regions with the equivalent quantity α_sm computed in the simulated region after emergence. Specifically, we define ${\alpha }_{\mathrm{sm}}=\sum {J}_{r}\,{B}_{r}/(\sum {B}_{r}^{2})$ (in Mm⁻¹) where J_r and B_r are calculated at r = 1.0154 R_⊙. To minimize the effect on α_sm of any nonpotentiality in the surrounding (pre-existing) magnetic field, this computation is done not during the full simulation but independently for each region, emerging it into a potential field initialized on the day before emergence. The average quantity α_sm is then computed at the end time of the emergence.

To compare α_ob and α_sm for each region, we check two factors: first, whether they have the same sign, and second, whether the ratio ${r}_{\mathrm{ob}}^{\mathrm{sm}}={\alpha }_{\mathrm{sm}}/{\alpha }_{\mathrm{ob}}$ is close to unity. Table 1 shows the results of evaluating the ratio when we use three alternative values of c: 2, 3, and 4 Mm. The first row shows the number of regions where α_sm and α_ob have opposite signs. Note that this number decreases from 19 to 15 with increasing c. The remaining rows depict the statistics where α_sm and α_ob have the same sign—true for more than 80% of regions. Those with ratios $0.5\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 1.5$ best match observation, and we call them "good-fit" regions. Note that the maximum number of good-fit regions are obtained for c = 3 Mm. Also, the number of regions with a very large ratio (${r}_{\mathrm{ob}}^{\mathrm{sm}}\geqslant 6.0$) is minimized when c = 3.0 Mm. This suggests c = 3.0 Mm is the optimum choice for the proportionality constant.

Table 1. Ratio between α_sm and α_ob

${r}_{\mathrm{ob}}^{\mathrm{sm}}={\alpha }_{\mathrm{sm}}/{\alpha }_{\mathrm{ob}}$	c = 2.0	c = 3.0	c = 4.0
${r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 0.0$	19	15	15
$0.0\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 0.5$	18	15	12
$0.5\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 1.5$	31	33	31
$1.5\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 2.5$	12	12	12
$2.5\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 3.5$	7	11	12
$3.5\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 4.5$	3	3	5
$4.5\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 6.0$	1	3	3
$6.0\leqslant {r}_{\mathrm{ob}}^{\mathrm{sm}}$	11	10	12

Download table as:  ASCIITypeset image

Even for the optimum choice of c, there are a considerable number of SHARP regions with a mismatch between α_sm and α_ob. We have investigated all of the regions where ${r}_{\mathrm{ob}}^{\mathrm{sm}}$ deviated significantly from unity (either wrong sign or ${r}_{\mathrm{ob}}^{\mathrm{sm}}\geqslant 1.5$), and find that they fall into several classes. For example, 16 were tiny active regions where the value and sign of the simulated and observed α will be heavily biased by a few grid points. Seven regions emerged very near to pre-existing active regions. A total of 11 regions had complex, mixed-polarity B_r distributions rather than simple bipoles, and six had two separate active regions within the same SHARP. The remaining 29 active regions were influenced by the current in the surrounding field. In particular, 15 among those 29 regions had the correct sign within the core of the region, but the overall measured α_sm had either negative ${r}_{\mathrm{ob}}^{\mathrm{sm}}$ or ${r}_{\mathrm{ob}}^{\mathrm{sm}}\lt 0.5$. For these regions, increasing the value of c reduced the discrepancy. The effect of "wrongly" imposed twist in the complex regions and SHARPs with two spots (total 17 regions) on the global scenario will be discussed later (in Section 5.2).

As prescribed in Section 2.2, each active region emerges over a time interval of duration T_em, during which a steady electric field E _⊥(s, ϕ) is imposed. Note that, in any given simulation, our code uses the same value of T_em for all 102 active regions. While keeping the twist parameter fixed at τ = ±0.1 (negative for all regions in the northern hemisphere, positive for all regions in the southern hemisphere), we considered three different values of T_em: 6, 12, 24 hr. The effect is illustrated for an individual emerging region in Figure 9. We see that there is no qualitative difference in the magnetic field structure, the field line helicities (the colors of the field lines in the plot), or the current density distribution. However, we do see an increase in the peak current density when the emergence is carried out more rapidly, due to the reduced time for relaxation available. Another reason for the difference in current is the fact that supergranular diffusion continues to smooth out the photospheric field during the emergence process, at a rate independent of T_em.

Zoom In Zoom Out Reset image size

The effect of choosing different T_em on the global diagnostics is shown in Figure 10, where we focus on a subinterval (September 2011) for clarity. In Figures 10(a) and (b), we find that decreasing the emergence duration, in general, corresponds to a little more magnetic energy and current density in the corona at the epochs of emergence (marked by the asterisks). However, these short-lived increments decay quickly, and the overall variation in energy and current density show significant similarities even with different T_em. Evolution of the other quantities, like open magnetic flux, helicity flux, and B_⊥ in the northern and southern hemispheres (Figures 10(c)–(f)) do exhibit relatively higher differences. In particular, predictions of specific eruption timings from the model must be viewed as uncertain and likely require a more detailed assimilation of input data. Still, we notice no systematic overall effect of T_em on the global evolution when we consider the total time-integrated amplitude of open flux, helicity flux, the total B_⊥ (in each hemisphere), and the total number of peaks at the end of each simulation. The differences remain <2.5% for all of the time-integrated quantities when comparing the three simulations. Thus in the remainder of the paper, we fix the standard value T_em = 24 hr, which is consistent with typical observations of sunspot emergence (van Driel-Gesztelyi & Green 2015).

Zoom In Zoom Out Reset image size

We have performed 13 global magnetofrictional simulations with the same 102 emerging regions but different τ, as described in Section 4. For each of these simulations, T_em was set at 24 hr. These comprise three sets of simulations with constant ∣τ∣: 0.0, 0.05, 0.1, 0.2 (sign according to the hemispheric rule), constant ∣τ∣: 0.05 s, 0.1 s, 0.2 s (sign according to observed sign of α_ob), and τ = c α_ob: 2α, 2.5α, 3α, 3.5α, 4α, 10α. (The coefficients denote the proportionality constant c in Mm.) As stated previously, we primarily focus on the evolution of the global coronal magnetic field. Thus we compute the total time-integrated amplitude of open flux, total energy, average current density, helicity flux, average B_⊥ (in each hemisphere), and the number of peaks in B_⊥ (in each hemisphere). A detailed comparison of those quantities for individual simulations is depicted in Figure 11 where the horizontal axis represents different choices of τ.

Zoom In Zoom Out Reset image size

We can see that all of the measures have minimum amplitude when τ = 0.0, which corresponds to ${{\boldsymbol{E}}}_{\perp }={{\boldsymbol{E}}}_{\perp }^{0}$. Comparing the pair of simulations [0.05 & 0.05s], [0.1 & 0.1s], and [0.2 & 0.2s], we find that the simulations with signs of τ according to observed α_ob result in lower amplitudes of the total open flux, magnetic energy, current density, number of peaks, etc. (for example, see the pairs of data points highlighted by gray boxes Figure 11). The interpretation is that emerging regions with both signs of helicity tend to reduce the net coronal helicity in each hemisphere, allowing the coronal field to relax to a state with less free energy and current, leading to fewer and/or weaker eruptions.

The simulations where τ = c α_ob have increasing current, magnetic energy, and helicity flux with increasing c, and generally increasing numbers of eruptions. This increasing trend is visible in the data points encircled in Figure 11. Comparing the cases with constant ∣τ∣ and τ = c α_obs, we can say that the case with c = 3.0 Mm is globally similar to ∣τ∣ = 0.05 (with observed sign). On the other hand, the simulation with c = 10.0 Mm is comparable with the ∣τ∣ = 0.1 simulation, although in this case the match between α_sm and α_ob is poorer (Section 4.3).

Not shown in Figure 11, we have performed additional global simulations with a higher constant twist amplitude, ∣τ∣ = 0.4 (with both uniform and observed signs). With such strong twists, most of the active regions erupt immediately after emergence, and we feel that the observed twisting of magnetic field lines is unrealistically high compared to typical observations of coronal loops. This excessive twist is illustrated in Figure 7(c). As such, we feel that the simulations with ∣τ∣ = 0.1 or τ = 10α_ob provide a rough upper limit to the reasonable amount of helicity that should be injected into the emerging regions.

As we noted in Section 4.3, the ratios ${r}_{\mathrm{ob}}^{\mathrm{sm}}$ for a number of regions with complex B_r distributions differed substantially from the ideal value of unity. Of these regions, 17 were either complex multipolar regions or SHARPs containing two separate active regions. Hence to evaluate the effect of the "wrong" twist of these regions on the global scenario, we ran a further global simulation with the twist of these regions set to zero and the remaining 85 regions left as τ = 3.0 α_ob. This modification causes slight decreases in the total open flux (2%), magnetic energy (1%), and current density (2.7%). Interestingly, quantities closely associated with eruptions through the upper boundary showed somewhat more significant effects: total helicity flux was reduced by 4.2%, and the number of peaks in B_⊥ decreased by 5.2% and 12.3% in the northern and southern hemispheres, respectively (see the data points within individual rectangles and corresponding to τ: 3α m in Figure 11). This gives an indication of the importance of choosing τ carefully.

Lastly, we performed two more simulations with constant ∣τ∣ = 0.0 and 0.1 (with signs according to the hemispheric rule) but used the global inductive solution to compute E _⊥ for each of the active regions (see Figures 3(c) and (g)). Even for the run with τ = 0.0, the global inductive simulation generated higher amounts of total open flux (5.3%), magnetic energy (1.6%), current density (5.4%), and helicity flux (4.9%), as well as more peaks in B_⊥ (North: 6.9% and South: 5.5%) as compared to the local solution (see the data points corresponding to τ = 0.0 G, which are encircled separately in Figure 11). Running another such simulation with τ = 0.1 (hemispheric signs) only for 2011 June revealed that the differences increased further when the chosen twist amplitude was higher.