STITCH: A Subgrid-scale Model for Energy Buildup in the Solar Corona

Solar eruptions are the most energetic explosive events in the solar system. The largest of these events, known as coronal mass ejections (CMEs), release vast quantities of solar plasma and magnetic fields into the solar wind, generating shocks that disturb the entire heliosphere, driving hazardous energetic particles and magnetospheric storms at Earth. A major goal of solar modeling efforts is operational prediction of the events that present the greatest threats to human assets and technology. A critical obstacle to achieving this goal, however, is the paucity of quantitative measurements of the coronal magnetic field. Further complicating matters are the remarkable diversity of eruptions in size, from CMEs that fill the heliosphere down to coronal jets that span only a few megameters (Raouafi et al. 2016); the morphology of the erupting material in coronagraphs, including both so-called stealth CMEs (Ma et al. 2010) and halo CMEs (Howard et al. 1982); and the accompanying signatures in the low atmosphere, in particular the chromospheric flare ribbons, which can be as simple as the classic, linear two-ribbon structure or as complex as circular (Masson et al. 2009) or multiribbon (Wang et al. 2014) configurations, and the overlying arcades of flare loops in the corona.

Despite such diversity in eruptive manifestations, the energy source is generally agreed (at least for the largest events) to be filament channels, highly stressed magnetic structures that form above reversals in the photospheric line-of-sight field known as polarity inversion lines (PILs; Gaizauskas 1998; Martin 1998). Filament channels are force-balanced structures consisting of a highly sheared magnetic field (that is, approximately parallel to the PIL) restrained by an overlying unsheared field that links the opposite-polarity fluxes to either side of the PIL. Many of these channels are observed to exhibit S-shaped sigmoids in their hot coronal loops and/or cool filament threads (if present), with associated strong electric currents signaling the storage of magnetic free energy (Komm et al. 2015). An eruption occurs when this force balance is disrupted by some trigger, either an ideal instability (Fan 2001; Linker et al. 2003; Kliem & Torok 2006) or magnetic reconnection (Antiochos 1998; Antiochos et al. 1999; Moore et al. 2001), which initiates ejection of the built-up shear field. Understanding how the eruption will occur and the hazard it presents reduces to understanding how the generic filament channel formation and eruption processes occur in the complex and diverse magnetic field configurations that house the PILs.

Several models have been proposed for filament channel formation. The simplest example, perhaps, is large-scale shear flows or sunspot/active region rotation (e.g., Cheung & DeRosa 2012; Sun et al. 2012). However, such large-scale coherent motions are not always observed and are unlikely to explain all examples of filament channel formation (Gaizauskas 1998; Martin 1998). For the quiet Sun especially, it has often been suggested that a combination of differential rotation and flux cancellation generates coronal flux ropes that constitute filaments and their channels (van Ballegooijen & Martens 1989; van Ballegooijen 2004; Mackay & Yeates 2012; Yeates & Mackay 2012). However, the chirality of these structures, if arising from these processes alone, has been shown (Mackay et al. 2014, 2018, and references therein) to be inconsistent with the global helicity patterns observed at all scales within eruptive coronal structures. For example, Ouyang et al. (2017) found that 92% of their sample of 571 erupting filaments followed the empirical rule that those in the northern/southern hemisphere exhibit negative/positive magnetic helicity.

The recently developed helicity condensation model (Antiochos 2013) posits injection of helicity/energy into the corona via the small-scale flows associated with photospheric convective motions. In a sufficiently "turbulent" system supporting widespread reconnection across ubiquitous small-scale current sheets, this structure and associated energy will be transferred to increasingly larger scales via the well-known inverse cascade of helicity (Berger 1984; Berger & Field 1984; Finn & Antonsen 1985; Biskamp 1993). In the context of the Sun, this cascade causes small-scale twist injected throughout an area enclosed by a PIL to eventually "condense" at that PIL in the form of a large-scale shear. This theory has been rigorously explored through MHD simulations that demonstrate the formation of filament channel structures consistent with those observed on the Sun (Knizhnik et al. 2015, 2017a, 2017b, 2018; Zhao et al. 2015) and evolution of the structure through the generation of a fast CME with an eruptive flare (Dahlin et al. 2019). Importantly, the filament channels generated by helicity condensation, when it is added to the global-scale processes of differential rotation, meridional flow, flux cancellation, and newly emerged active region flux, are consistent with the observed chirality patterns on the Sun (Mackay et al. 2014, 2018).

All of the previous full MHD studies of helicity condensation, cited above, have used close-packed (i.e., tessellated), boundary-driven vortical flows to represent helicity injection via convective granule/supergranule motions. Such calculations are computationally expensive, as they must resolve all of the cellular flows and the small-scale reconnection dynamics that drive the helicity cascade. Additionally, the flows cannot overlap with the PILs lest they trigger excessive small-scale diffusion across structures generated at the grid scale, a challenging constraint for driving the highly complex magnetic configurations typical of eruptive events. An important implication of the helicity condensation model, however, is its relative insensitivity to the particular mechanism or spatial scale of the helicity injection (Knizhnik et al. 2017a, 2019). So long as there is a net helicity injection with sufficiently complex turbulent dynamics, the shear will be transferred to the PILs. Hence, while resolving the small-scale reconnection is critically important for understanding nanoflare-type coronal heating, ultimately the large-scale energy buildup is insensitive to the underlying details of the helicity injection and cascade processes. This motivates the search for a method that captures the essential large-scale features of the energy buildup without explicitly calculating all of the small-scale structure and dynamics. Such a method, if sufficiently flexible in its implementation, should be agnostic with respect to the particular mechanism for helicity injection. Hence, it could be used more generally, even in cases where turbulent convective motions are not expected to be the dominant source of energy and helicity.

In Mackay et al. (2014, 2018), a statistically averaged approximation to the full helicity condensation (FHC) process was derived and employed to investigate the formation of filament channels across the full Sun and over multiple solar rotations. Those studies used a magnetofrictional model for the corona that accurately reproduces the large-scale and long-term formation of filament channels, but which cannot self-consistently capture explosive eruptive dynamics because the model output is a sequence of force-free coronal magnetic field equilibria. In this article, we extend and test this statistical approach for general use in full MHD simulations. The resulting model, which we call STatistical InjecTion of Condensed Helicity (STITCH), encapsulates the injection of helicity into the corona by vortical convection cells, the formation of current sheets at boundaries between adjacent cells, and the reconnection of the induced horizontal magnetic fields across the sheets. STITCH has three critical advantages over the MHD treatment of the FHC process with all its attendant small-scale dynamics and structures. First, it is numerically efficient in both memory and computations, requiring substantially lower spatial resolution and temporal cadence. Second, it is simple and ease to use. Third, it is highly flexible, and it is readily applicable to the complex magnetic field distributions observed on the Sun.

We discuss the STITCH model and its numerical implementation in Section 2. Our magnetohydrodynamics code, ARMS, and initial configuration are described in Section 3. We compare STITCH to the comprehensive helicity condensation model in Section 4. Thereafter, we show examples of the flexibility of STITCH in localizing the filament channel (Section 5) and in being used to energize more complex and realistic magnetic configurations (Section 6). STITCH is then placed in the larger context of methods used to model solar eruptive events in Section 7. To conclude, we summarize our results and discuss the implications of using STITCH in simulations of coronal energy buildup in Section 8.

In his exposition of the helicity condensation model, Antiochos (2013) described how an ensemble of vortical motions at the photosphere builds up magnetic shear in the corona. Current sheets form at the interfaces between pairs of neighboring vortical cells that have the same sense (clockwise or counterclockwise) of rotation. Magnetic reconnection across the sheets depletes the twist component of magnetic flux between the cells, leaving behind the twist flux that envelopes both cells jointly. This process repeats in an inverse cascade that transports the twist flux to ever-larger spatial scales. The cascade halts at the boundary of the flux system—defined by a reversal of either the sign of the vertical magnetic field or the sense of rotation of the vortical cells—because the twist flux has the same direction on both sides of such a boundary and, hence, cannot reconnect away. The twist flux encircling the flux system is said to "condense" at that boundary, and the linkage between the twist flux and the enclosed vertical flux imparts magnetic helicity to the structure. The term "helicity condensation" had been coined previously to describe this process when it occurs in laboratory experiments and numerical simulations (Biskamp 1993).

Several subsequent investigations of helicity condensation as applied to the corona have been carried out (Knizhnik et al. 2015, 2017a, 2017b, 2018; Zhao et al. 2015; Dahlin et al. 2019) in which numerous small-scale vortical cells—as many as 100+—were resolved on the computational grid. The needed helicity-preserving simulations were performed with the Adaptively Refined Magnetohydrodynamics Solver (ARMS; DeVore & Antiochos 2008). The processes of twist-flux accumulation, current-sheet formation and reconnection, and helicity condensation that were predicted theoretically were confirmed quantitatively and in detail. However, the required well-resolved, long-duration calculations of this type are very resource intensive.

In a separate line of investigation that applied the helicity condensation concept to the full-Sun magnetic field over many weeks of evolution, we developed and employed a subgrid-scale representation of the model (Mackay et al. 2014, 2018). The essential assumptions that were made are as follows: (1) the vortical cells are small, numerous, and nominally identical, although their properties may exhibit large-scale variations; and (2) the reconnection across current sheets between the cells is so efficient that the oppositely directed twist fields can be treated as simply canceling algebraically. As derived in the Appendix of Mackay et al. (2014, Equation (A1)), the resulting subgrid-scale model is expressed as an additional term in the induction equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}{\boldsymbol{\times }}\displaystyle \frac{\partial {{\boldsymbol{A}}}_{s}}{\partial t}=-{\boldsymbol{\nabla }}{\boldsymbol{\times }}{{\boldsymbol{\nabla }}}_{s}\left(\zeta {B}_{n}\right).\end{eqnarray} \tag{ 1 }$$

Here B is the magnetic field; A is the vector potential; subscripts n and s denote components normal and parallel to the surface, respectively; and ζ is a parameter characterizing the vortical flows. Specifically,

$$\begin{eqnarray}&&\zeta =\left\langle {{\ell }}^{2}{\omega }_{l}\right\rangle /2,\end{eqnarray} \tag{ 2 }$$

where ℓ is the radius and ω_l is the angular rotation rate of the cells, and the angle brackets denote a local average in space and time. The parameter ζ has dimensions of diffusivity, and, as will be seen below, the local rate of injection of helicity (per unit area and time) into the corona is determined directly by the product $\zeta {B}_{n}^{2}$ (Mackay et al. 2014). From the form of Equation (1), it is clear that the changes in B are strongest where the product ζ B_n varies most rapidly across the surface. In particular, this occurs at the flux system boundaries described previously. We refer to the subgrid-scale model represented by Equation (1) as STITCH.

The STITCH term has several notable properties. First and most fundamentally, by construction it preserves the divergence-free character of the magnetic field,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}=-{\boldsymbol{\nabla }}\cdot \left[{\boldsymbol{\nabla }}{\boldsymbol{\times }}{{\boldsymbol{\nabla }}}_{s}\left(\zeta {B}_{n}\right)\right]=0,\end{eqnarray} \tag{ 3 }$$

subject only to the requirement that ζ B_n be sufficiently smooth (i.e., differentiable) everywhere. Second, it prescribes an injection of twist flux into the corona that leaves the normal magnetic field component B_n unchanged,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{n}}{\partial t}=-\hat{{\boldsymbol{n}}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{\times }}{{\boldsymbol{\nabla }}}_{s}\left(\zeta {B}_{n}\right)=-\hat{{\boldsymbol{n}}}\cdot {{\boldsymbol{\nabla }}}_{s}{\boldsymbol{\times }}{{\boldsymbol{\nabla }}}_{s}\left(\zeta {B}_{n}\right)\equiv 0.\end{eqnarray} \tag{ 4 }$$

The change in the horizontal field component B _s is

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{B}}}_{s}}{\partial t}=-{{\boldsymbol{\nabla }}}_{n}{\boldsymbol{\times }}{{\boldsymbol{\nabla }}}_{s}\left(\zeta {B}_{n}\right)=+{{\boldsymbol{\nabla }}}_{s}{\boldsymbol{\times }}{{\boldsymbol{\nabla }}}_{n}\left(\zeta {B}_{n}\right);\end{eqnarray} \tag{ 5 }$$

we note that the exchange of curl and gradient operators in Equation (5) is valid in both spherical and Cartesian coordinates. Third, the total rate of injection of relative helicity, H_r , into the atmosphere is given by the surface integral of the local rate due to the individual vortical cells (for details, see Mackay et al. 2014),

$$\begin{eqnarray}&&\displaystyle \frac{{{dH}}_{r}}{{dt}}=-2\int \,\int \zeta {B}_{n}^{2}d\sigma .\end{eqnarray} \tag{ 6 }$$

This is required in order for the STITCH term to be consistent with helicity conservation throughout the processes of reconnection and transport of twist flux to the flux system boundaries.

STITCH is implemented numerically by computing the rate at which horizontal magnetic flux is injected into the domain according to Equation (5). This flux is calculated on the cell faces oriented perpendicular to the two horizontal coordinates s₁ and s₂, where $\hat{{{\boldsymbol{s}}}_{{\bf{1}}}}{\boldsymbol{\times }}\hat{{{\boldsymbol{s}}}_{{\bf{2}}}}=\hat{{\boldsymbol{n}}}$. Integrating Equation (5) over differential cell areas ${{dA}}_{{s}_{1}}$ and ${{dA}}_{{s}_{2}}$, respectively, we obtain the flux change rates

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{1}}}{\partial t} & = & \int \,\int {{dA}}_{{s}_{1}}\hat{{{\boldsymbol{s}}}_{{\bf{1}}}}\cdot \displaystyle \frac{\partial {{\boldsymbol{B}}}_{s}}{\partial t}\\ & = & +\int \,\int {{ds}}_{2}{dn}\displaystyle \frac{\partial }{\partial {s}_{2}}\displaystyle \frac{\partial }{\partial n}\left(\zeta {B}_{n}\right),\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{2}}}{\partial t} & = & \int \,\int {{dA}}_{{s}_{2}}\hat{{{\boldsymbol{s}}}_{{\bf{2}}}}\cdot \displaystyle \frac{\partial {{\boldsymbol{B}}}_{s}}{\partial t}\\ & = & -\int \,\int {{ds}}_{1}{dn}\displaystyle \frac{\partial }{\partial {s}_{1}}\displaystyle \frac{\partial }{\partial n}\left(\zeta {B}_{n}\right).\end{array}\end{eqnarray} \tag{ 8 }$$

These expressions are coordinate system independent. We localize the angular rotation of the convection cells near the coronal base, which is positioned at height n = n₀, by assuming that ζ falls to zero at the top of the first layer of grid cells at height n = n₁ = n₀ + h. Integrating Equations (7) and (8) over n from n₀ to n₁, they simplify to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{1}}}{\partial t}=-\int {{ds}}_{2}\displaystyle \frac{\partial }{\partial {s}_{2}}\left(\zeta {B}_{n}{| }_{{n}_{0}}\right),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{2}}}{\partial t}=+\int {{ds}}_{1}\displaystyle \frac{\partial }{\partial {s}_{1}}\left(\zeta {B}_{n}{| }_{{n}_{0}}\right),\end{eqnarray} \tag{ 10 }$$

where $\zeta {B}_{n}{| }_{{n}_{0}}$ denotes the value of ζ B_n evaluated at the surface, n = n₀. These expressions integrate immediately to yield the flux changes

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{1}}}{\partial t}=-{{\rm{\Delta }}}_{{s}_{2}}\left(\zeta {B}_{n}{| }_{{n}_{0}}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{2}}}{\partial t}=+{{\rm{\Delta }}}_{{s}_{1}}\left(\zeta {B}_{n}{| }_{{n}_{0}}\right),\end{eqnarray} \tag{ 12 }$$

where ${{\rm{\Delta }}}_{{s}_{1}}$ (${{\rm{\Delta }}}_{{s}_{2}}$) denotes a finite difference whose argument is evaluated at cell vertices along coordinate s₁ (s₂).

The flux changes on the four faces of any grid cell involve pairwise differences of the argument $\left(\zeta {B}_{n}{| }_{{n}_{0}}\right)$ evaluated at the four vertices of the (s₁, s₂) grid. Summing these changes algebraically, to calculate the change in the net flux leaving the cell through its four faces, yields zero; each vertex contribution appears twice in the sum with opposite signs (+, −). This is simply a discrete expression of Gauss's law, Equation (3), integrated over the cell volume.

The final step in the calculation is to update the horizontal magnetic field components using the flux changes in Equations (11) and (12),

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{{s}_{1}}}{\partial t}=\displaystyle \frac{1}{{A}_{{s}_{1}}}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{1}}}{\partial t}=-\displaystyle \frac{1}{{A}_{{s}_{1}}}{{\rm{\Delta }}}_{{s}_{2}}\left(\zeta {B}_{n}{| }_{{n}_{0}}\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{{s}_{2}}}{\partial t}=\displaystyle \frac{1}{{A}_{{s}_{2}}}\displaystyle \frac{\partial {{\rm{\Phi }}}_{{s}_{2}}}{\partial t}=+\displaystyle \frac{1}{{A}_{{s}_{2}}}{{\rm{\Delta }}}_{{s}_{1}}\left(\zeta {B}_{n}{| }_{{n}_{0}}\right),\end{eqnarray} \tag{ 14 }$$

where ${A}_{{s}_{1}}$ (${A}_{{s}_{2}}$) is the area of the cell face perpendicular to the coordinate s₁ (s₂). Evaluating the $\left(\zeta {B}_{n}{| }_{{n}_{0}}\right)$ vertex values and multiplying them by the finite-difference time increment, Δt, yields the required changes to the field components, ${\rm{\Delta }}{B}_{{s}_{1}}$ and ${\rm{\Delta }}{B}_{{s}_{2}}$.

Using STITCH in a numerical simulation requires specifying the magnitude of the parameter ζ to be used. By appeal to observations of supergranules (Duvall & Gizon 2000; Gizon & Duvall 2003; Komm et al. 2007), together with modeling results from flux transport simulations, Mackay et al. (2014) estimated the solar value to fall in the range (6 ± 4) × 10² km²s⁻¹. As expected, this is also the magnitude of the corresponding diffusion coefficient with which supergranules disperse the surface magnetic flux. The value is set by the supergranules' relatively large size (ℓ ≈ 2 × 10⁴ km) but very low vortical speed (ℓ ω_ℓ ≈ 6 × 10⁻² km s⁻¹). In MHD calculations that are performed over much shorter time intervals, of necessity, this value is impractically small to have any significant effect. Fortunately, adopting substantially larger values in MHD simulations is no barrier, so long as the resulting shear flux injection occurs more slowly than the governing fast MHD processes in the system, such as Alfvén-wave propagation. Accelerating the helicity condensation process by enhancing the STITCH coefficient is akin to accelerating the surface flows imposed in many other simulations to form and energize filament channels (e.g., Karpen et al. 2012). So long as the injection speed is sub-Alfvénic—preferably, if practical, also subsonic—the system evolves quasi-statically as the magnetic free energy builds up, until it transitions to fast, dynamical behavior that may include a violent eruption. The only caveat is that the duration of the energy buildup phase is artificially shortened compared to that required on the Sun; the critical point, however, is that both the real and model buildup phases are much more slowly evolving than the eruption phase.

We calculate a constraint on the allowed value of the STITCH parameter ζ by limiting the speed of propagation across the grid of any induced horizontal field component δ B _s to the local Alfvén speed, as follows. We assume that the grid spacings are Δs horizontally and h vertically (in practice, these three are similar if not identical), the magnetic field components are δ B_s horizontally and δ B_n vertically (these two may differ substantially in magnitude), and δ B_n varies at the smallest accessible scale (i.e., the horizontal grid spacing Δs). The STITCH induction Equation (5) then yields directly

$$\begin{eqnarray}&&\displaystyle \frac{\partial \delta {B}_{s}}{\partial t}\approx \displaystyle \frac{\zeta \delta {B}_{n}}{h{\rm{\Delta }}s}.\end{eqnarray} \tag{ 15 }$$

We can associate the rate of change of δ B_s with an effective shear speed V_s applied to the footpoints of the magnetic field at the line-tied surface, with the result

$$\begin{eqnarray}&&\displaystyle \frac{\partial \delta {B}_{s}}{\partial t}\approx \displaystyle \frac{{V}_{s}\delta {B}_{n}}{h}.\end{eqnarray} \tag{ 16 }$$

Equating the above two rates, we find the relationship between V_s and ζ,

$$\begin{eqnarray}&&{V}_{s}\approx \displaystyle \frac{\zeta }{{\rm{\Delta }}s}.\end{eqnarray} \tag{ 17 }$$

With this identification of the grid-dependent (due to Δs) effective shear speed V_s , we now define a flux injection Courant number _fi according to

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{fi}}\equiv \displaystyle \frac{{V}_{s}}{{V}_{{\rm{A}}}}\equiv \displaystyle \frac{\zeta }{{V}_{{\rm{A}}}{\rm{\Delta }}s},\end{eqnarray} \tag{ 18 }$$

where V_A is the characteristic Alfvén speed in the active region. As noted below, we found empirically that _fi ≈ 1/6 was a satisfactory upper limit on ζ for the STITCH simulations reported in this paper. This result is entirely in line with the numerical time increments Δt (relative to the minimum signal crossing time on the grid) that are adopted in MHD simulations, where fractional values on the order of a few tenths are acceptable and are used routinely to obtain numerically stable solutions.

The net result of the above considerations is that we can set an upper limit on the allowed STITCH parameter ζ simply by adopting a suitable value for the flux injection Courant number _fi,

$$\begin{eqnarray}&&{\zeta }_{\max }={\epsilon }_{\mathrm{fi}}{V}_{{\rm{A}}}{\rm{\Delta }}s.\end{eqnarray} \tag{ 19 }$$

This numerically constrained upper bound depends on the physical Alfvén speed, V_A, as should be expected, but also on the assumed grid spacing, Δs. The physical spatial scales in the simulation, such as the width w of the filament channel or the diameter L of the active region, are larger than the grid spacing. Consequently, the Courant numbers for processes at these scales are progressively smaller than _fi; for example, that for the shear injection of horizontal flux at the channel-width scale w is (Δs/w)_fi, and that for injecting an amount of horizontal flux equal to the vertical flux in the active region at its diameter L is (Δs/L)_fi. The physically intuitive requirements that these two subsidiary parameters be small are well satisfied by our constraint on ζ expressed in Equation (19).

Our numerical calculations were performed with ARMS (DeVore & Antiochos 2008), which has been used extensively to model both CMEs/eruptive flares (see also Lynch et al. 2008, 2009, 2016, 2021; Karpen et al. 2012; Masson et al. 2013) and the formation of filament channels (Knizhnik et al. 2015, 2017a, 2017b, 2018; Zhao et al. 2015). Our fiducial calculations adopt the magnetic field configuration shown in Figure 1. It consists of an elliptical, bipolar active region centered at 225N latitude with peak radial field ∣B_r ∣ ≈ 50 G, embedded in a background 10 G solar dipole field (for full description see Dahlin et al. 2019). The resulting topology is the well-known embedded bipole with a separatrix dome and a pair of spine lines emanating from a 3D null point (Antiochos 1990; Lau & Finn 1990; Priest & Titov 1996). A region of maximally refined grid enclosed the entire separatrix dome (gray/silver field lines arching over the arcade of red loops; Figure 1) and a substantial volume above it, in order to resolve the eventual small-scale reconnection dynamics. The initial atmosphere was a spherically symmetric hydrostatic equilibrium with an inverse-r temperature profile at a base temperature T_s = 2 × 10⁶ K and pressure P_s = 4 × 10⁻¹ dyn cm⁻². We solved the ideal MHD equations with an adiabatic temperature equation. The domain extents were r ∈ [1R_s , 30R_s ], θ ∈ [π/16, 15π/16], and ϕ ∈ [ − π, + π].

Zoom In Zoom Out Reset image size

For the FHC calculation, the system was driven by 107 tessellated, vortical, B_r -conserving cellular flows within the black shaded minority-polarity region. The flows were subsonic and sub-Alfvénic, attaining a maximum speed ∣V_⊥∣ ≈ 50 km s⁻¹ after an initial sinusoidal ramp-up interval 1000 s in duration. Elsewhere on the surface, the magnetic field was line-tied at rest ( V _⊥ = 0) for the entire duration of the simulation. The resultant rates of change of the horizontal magnetic field components are shown in Figures 2(a) and (d). An extended energy buildup phase occurred over about 90,000 s, which required approximately two continuous months of computational time on the NCCS discover supercomputer at NASA GSFC.

Zoom In Zoom Out Reset image size

For our first STITCH calculation described below, we assumed that ζ was uniform at ζ₀ = 1.4 × 10⁶ km² s⁻¹ within the minority-polarity region and zero everywhere outside of that region. This corresponds to a flux injection Courant number _fi ≈ 1/6 (Equation (18)), given the characteristic Alfvén speed V_A = 4 × 10³ km s⁻¹ and minimum grid spacing Δs = 2 × 10³ km.

The resulting profiles of injected horizontal magnetic field, whose maximum rate was about 1 G s⁻¹, are shown in Figures 2(b) and (e) (Figures 2(c) and (f) show a STITCH case with localized ζ that is discussed in detail in Section 5). In this and other STITCH calculations, the magnetic field was line-tied at rest ( V _⊥ = 0) over the whole surface for the entire duration of the simulation. To achieve a smooth start-up, we ramped up the tangential field injection sinusoidally over the first 1000 s of elapsed time in the simulation. The energy buildup phase was much shorter with STITCH, in this case being only about 6000 s in total duration and requiring approximately 4 days of computational time.

Although the maximum rates of horizontal field injection are similar in the FHC and STITCH cases, as shown by Figure 2, the timescales for accumulating the helicity and free energy necessary to drive eruption differ by a substantial factor of 15. One reason for this is that the width of the injection region is much larger, by a factor of about 4, for STITCH versus FHC; the rates of shear flux accumulation differ by the same factor. Whereas the STITCH profiles are smooth and space filling, the FHC profiles are patchy owing to the cellular structure; hence, FHC has a filling factor of about 0.5, compared to unity for STITCH. Lastly, the changing direction of the vortical flows in the FHC case lowers the average effective injection rate parallel to the PIL; this provides yet another reduction factor of about 2 in efficacy of FHC compared to STITCH. These cumulative effects account for the major difference in the required drive times to eruption in the two cases.

The stated rate of field injection and time interval in the STITCH case suggest that a tangential field strength on the order of 6 kG should build up; however, this does not occur. The tangential field is injected just above the surface, but then it is redistributed all along the length of the filament channel field lines by Alfvén waves. This reduces the accumulated maximum tangential field strength to approximately 100 G, roughly twice the peak radial field strength inside the minority polarity. Expressed another way, we find that the tangential flux injected by STITCH over the simulation's duration is roughly equal to the normal flux within the minority-polarity region.

Figures 3 and 4 show snapshots of the filament channel formation and eruption phases for the FHC and STITCH cases, respectively. Both generate low-lying, highly sheared orange field lines (the longest of which stretch more than halfway along the PIL) that are restrained by overlying, weakly sheared red field lines (Figures 3(a) and 4(a)). These snapshots are taken about two-thirds of the way through the filament channel formation phase for the two cases. Eventually, the entire filament channel erupts as a quasi-circular filament with maximum speed V > 1000 km s⁻¹, indicated by dark red shading in the plane of the sky (Figures 3(b) and 4(b)). The kinetic energy is near its maximum value at the times shown for both cases.

One significant point of contrast between the two cases is the shear imparted to the red loops overlying the PIL of the active region under FHC; under STITCH, the same loops are essentially shear-free. The helicity-injecting flows in the former fill the minority-polarity portion of the active region, inducing twist everywhere that subsequently condenses at the PIL, but with a finite time delay and with a spatial distribution that is smoothed across the region of flow. In contrast, the STITCH model with uniform ζ instantaneously injects tangential field adjacent to the PIL, across a zone whose characteristic width w = ∣B_r ∣/∣∇_s B_r ∣ is determined solely by the magnetic field profile. The constant stirring of the field within the filament channel by the helicity condensation flows means that the effective minimum width of the FHC channel will be the characteristic diameter of the vortical cells, d = 2ℓ (recall Equation (2)). Due to the limited resolution available for use in our simulations, these two measures of the width are comparable, d ≈ w, as can be seen directly in Figure 2. Therefore, the resulting filament channels exhibit similar widths of their strongly sheared core fields, as can be seen in Figures 3(a) and 4(a).

Figure 5 shows the evolution of the volume-integrated magnetic and kinetic energies, along with the maximum value of ∣B_ϕ ∣, a proxy for the strength of the shear magnetic field. Throughout the entire evolution, except for a brief period lasting about 100 s at the beginning of the driving phase, this maximum occurs at the center of the southern half of the PIL near the inner boundary.

Zoom In Zoom Out Reset image size

In order to compare the macroscopic evolution most conveniently, we have scaled the respective time axes to align the approximate eruption times in the two cases, t ≈ 90,000 s (FHC) and t ≈ 6000 s (STITCH). The reduction in the nominal driving time is matched by a corresponding decrease by a factor of 15 in the computational cost of the STITCH run. Qualitatively, all of the trends in the displayed global variables are shared between FHC and STITCH. Quantitatively, on the other hand, there are both close similarities and significant differences, as we now discuss in detail.

The accumulation of magnetic energy from its initial value, E_M ≥ E_M0, and its later reduction in conjunction with eruption, are shown in Figure 5(a). We find that the accumulated magnetic free energy is some 20% higher in the STITCH case than in the FHC case, E_M /E_M0 ≈ 1.24 and 1.20, respectively. This STITCH excess is expected, based on its concentrating all of its shear adjacent to the PIL, which we discussed above. The more broadly distributed shear in the FHC case distorts the overlying coronal null point into a current patch, facilitating the onset of breakout reconnection and removal of the high-lying magnetic tethers, earlier in the energy injection process than is true for the STITCH case. We emphasize that the elapsed time is so much longer for FHC only because its imposed vortical flows are so slow. Were we to inject magnetic free energy in the STITCH case at a comparably slow rate, we would expect the STITCH energy buildup phase to require more elapsed time than FHC, in order to accumulate the greater energy needed to initiate the STITCH eruption.

The strong concentration of the STITCH shear adjacent to the PIL can account for two more contrasting features of its energy profile compared to the FHC case. First, the STITCH magnetic energy release is only about 40% of that for FHC, ΔE_M /E_M0 ≈ 0.04 and 0.10, respectively, culminating in post-eruption energies E_M /E_M0 ≈ 1.20 and 1.10. The much narrower shear distribution achieved by STITCH reduces the volume of free energy containing magnetic flux, relative to that of FHC, and throttles back the amount of flux that is processed and energy that is released by the eruption. Second, the STITCH energy release timescale is only Δt ≈ 1000 s, a factor of four less than the Δt ≈ 4000 s timescale for the FHC energy release. (Notice that the STITCH energy release phase only appears to be much more gradual than that for the FHC phase, due to the difference in the absolute time axes for the two cases.) We would expect the smaller volume processed in the STITCH case to correspond to a shorter timescale, and this is exactly what we find. The greater concentration of the STITCH shear flux leads to a rather more explosive energy release, however, by a factor of 0.4 × 4.0 = 1.6, or a 60% increase relative to the FHC energy release rate.

For clarity, the energy storage and release rates dE_M /dt, normalized to the initial magnetic energy, are displayed separately in Figure 5(b). The storage rate for STITCH is just over 15× faster than that for FHC, whereas, as mentioned, its release rate is only about 60% faster. The energy release appears to subside about four times slower for STITCH versus FHC on this plot, due to the disparate time axes used, but in fact it subsides about four times faster, as we also noted above.

The kinetic energy in the simulation, normalized to the initial magnetic energy, is shown in Figure 5(c). This ratio is ${E}_{k}/{E}_{M0}=\langle {V}^{2}\rangle /\langle {V}_{{\rm{A}}}^{2}\rangle =\langle {M}_{{\rm{A}}}^{2}\rangle$, where V_A and M_A are the Alfvén speed and Alfvén Mach number, respectively. In both cases, the ratio during the long phase of slow driving is on the order of 10⁻⁴, or 〈M_A〉 ≈ 1%. The value jumps steeply during the explosive eruption phase to a few times 10⁻², or 〈M_A〉 ≈ 15%–20%. This rather high range highlights the pervasive Alfvénic outflows of the fast eruption, whose plasma occupies only a small fraction of the total volume. The FHC case has a larger volume of ejected flux and a higher peak kinetic energy content, some 50% more, than the STITCH case.

The shear field proxy, $| {B}_{\phi }{| }_{\max }$, is shown in Figure 5(d). This evolution, discussed previously in our work on the FHC simulation (Dahlin et al. 2019), exhibits a strong initial increase as the shear flux begins to accumulate above the PIL in the low corona. Eventually, although the total amount of shear flux continues to rise, the peak shear strength begins to decline, as the accumulating shear magnetic pressure inflates the overlying arcade loops and expands the volume occupied by the flux. This steady decrease continues until eruption onset. Although much of the shear flux is ejected upward as core field of the CME flux rope, a lower-lying portion is left behind in the reconnected flare loops. These loops relax downward, and their shear field strength increases as the occupied volume contracts. This generates a second, post-eruption peak in the shear proxy. The FHC and STITCH cases exhibit entirely similar behaviors, but due to the greater concentration of the STITCH filament channel, both of its peaks—during the early expansion and after the eruption (110/80 G, respectively)—are noticeably higher than those of the FHC filament channel (70/70 G). The local-minimum shear field strength at eruption onset, on the other hand, is very nearly the same in the two cases (about 40 G).

These results demonstrate that driving the filament channel formation process via the STITCH model produces an evolution leading to eruption that is very similar to our previous findings using FHC (Dahlin et al. 2019). The major difference between the two is the highly compressed characteristic timescale for the STITCH energization: whereas FHC relies on slow, subsonic vortical flows that gradually build up the sheared filament channel, STITCH relies on the much faster, but sub-Alfvénic, direct injection of horizontal magnetic flux into the corona. A more minor, but physically very significant, difference between them is in the width of the resulting filament channel: the minimum FHC channel width is determined by the characteristic diameter d of the vortical flows, whereas the STITCH channel width is determined entirely by the characteristic scale w of the vertical magnetic field profile. The limited numerical resolution available to us for our FHC calculation produced results for a case in which d is comparable to, but slightly larger than, w. In the true subgrid-scale regime where d ≪ w, in contrast, we would expect the resulting FHC channel width to be determined by the magnetic field length scale w. The FHC and STITCH results then should come into even better alignment than for the cases presented here, and they should show an even greater disparity in their respective timescales, further enhancing the advantages of STITCH for studies of filament channel formation.

The STITCH case presented above assumed uniform ζ throughout the minority-polarity region, best matching the uniform vortical flows in the FHC simulation. The resulting eruptions displayed in Figures 3 and 4 (best seen in the online animations) encompass the entire elliptical filament channel, proceeding first on the southern segment of the PIL and subsequently on the northern segment. This scenario is not universal, but it is observed routinely on the Sun (e.g., Wang et al. 2012). The sequence implies that the energy buildup is relatively uniform all along the PIL; hence, the entire sheared structure approaches its critical point and transitions to eruption nearly instantaneously. For the configuration discussed in Section 4, for example, the gradient in B_r across the PIL varies relatively little (by less than a factor of two) along the length. Therefore, the rate of STITCH flux injection likewise is quite uniform (Figures 2(b) and (e)), and the resulting eruption is a two-sided ejection of the entire channel.

More typically, however, solar eruptions occur along restricted segments of the filament channels and present one-sided ejections, even along quasi-circular PILs. For example, such an event is described by Mason et al. (2021). The implied local concentration of energy and shear buildup could arise as a result of several mechanisms, e.g., variable vortical granular/supergranular and/or large-scale shear flows, disordered flux cancellation, flux emergence, etc. The STITCH model allows us to abstract any and all such mechanisms into our simulations by simply adopting a spatially varying ζ profile, thereby localizing the helicity and energy injection to a selected region. We modified our uniform-ζ simulation in this way by centering a simple cosine variation on the southern PIL,

$$\begin{eqnarray}&&\zeta ={\zeta }_{0}\cos \left(\pi \displaystyle \frac{\phi -{\phi }_{0}}{2{\phi }_{w}}\right)\cos \left(\pi \displaystyle \frac{\psi -{\psi }_{0}}{2{\psi }_{w}}\right).\end{eqnarray} \tag{ 20 }$$

Here ϕ is longitude, ψ is latitude, ϕ₀ = 0°, ψ₀ = 225, ϕ_w = 225, and ψ_w = 1125. The amplitude ζ₀ and the restriction to the minority-polarity region are the same as in the previous simulation. Figures 2(c) and (f) show the resulting injection rates for B_ϕ and B_θ . The principal feature is that the injection is weakened substantially along the eastern, northern, and western segments of the PIL relative to the uniform-ζ case. A secondary feature is that the distribution of flux injection spreads from those segments of the PIL into the interior of the minority-polarity region, due to contributions to the injection rates made by the ζ gradients. Based on the results of our simulation shown previously, we ramped down the tangential field injection between elapsed times t = 6000 s and t = 7000 s, reducing the injection rate during the eruption. We also reran the uniform-ζ case with this change in the time profile, for consistency.

Snapshots of the evolving configurations at time t = 6200 s are shown in Figure 6 for the uniform-ζ (left; panels (a)–(c)) and localized-ζ (right; panels (d)–(f)) cases. Several significant differences are clearly evident from this comparison. In the uniform-ζ case, strongly sheared field lines populate the filament channel all along the PIL (Figure 6(a)); in the localized-ζ case, on the other hand, such field lines are observed only along the southern segment of the PIL, with weak shear present along the other segments by design (Figure 6(b)). The red arcade field lines overlying the southern filament channel are inflated to a greater extent in the localized-ζ case, due in part to less competition from the much more weakly sheared, northern segment of the PIL, and in part to the reversed shear that has been imparted to these field lines. This reversed shear arises from the relatively small, but finite, injection of oppositely signed B_ϕ flux near the center of the minority-polarity region due to the ζ gradients in the STITCH term (noted in the description of Figure 2(f)).

Side views of the azimuthal magnetic field B_ϕ show how localizing the STITCH injection (Figure 6(d)) greatly reduces the shear along the northern PIL segment and introduces the reversed shear in the center of the minority polarity (both shaded light red), relative to the uniform STITCH injection (Figure 6(c)). Obvious features of the uniform-ζ case (Figure 6(c)) are the roughly equal strengths of its northern and southern filament channels and the extents to which they bulge upward and outward as they inflate toward nearly simultaneous eruptions. In contrast, only the southern filament channel is clearly primed for eruption in the localized-ζ case (Figure 6(d)). These views also highlight a slight northward tilt of the southern filament channel field and its overlying arcade (shaded blue) in this case.

Side views of magnetic field lines in the four flux systems of the configuration further illustrate these features in the uniform-ζ (Figure 6(e)) and localized-ζ (Figure 6(f)) cases. The arcade of blue field lines in the north, along with its enclosed orange field lines, resides much lower in height in the latter versus the former. The color shading of radial velocity V_r in these panels captures the outward motion of the central arcade of red field lines in both cases, and it also illustrates breakout reconnection downflows along the flanks of this arcade in the localized-ζ case (Figure 6(f)).

Global diagnostics for these two simulations, uniform-ζ (black curves) and localized-ζ (green curves), are compared in Figure 7. (In these plots, unlike the previous similar Figure 5, the time axes for the two cases [both STITCH] are identical.) Figures 7(a), (b), and (c) show the energies for the uniform-ζ results at 50% amplitude, drawn with dashed curves. Coincidentally, these reduced curves track very well with those for the localized-ζ case, in the buildup to eruption for the magnetic energy (Figures 7(a) and (b)) and beyond into the eruption phase for the kinetic energy (Figure 7(c)). This agreement indicates that, during the buildup, the northern and southern segments of the PIL store roughly equal amounts of energy for uniform ζ. The magnetic energy release rate during the eruption is essentially the same in the two cases, although the release subsides somewhat faster for localized ζ (Figure 7(b)). The kinetic energy for that case peaks slightly above 50% of the uniform-ζ value (Figure 7(c)). Peak strengths of the shear field attained early in the energization phase are essentially identical for the two cases (Figure 7(d)). Early in the evolution for localized ζ, a transient expansion of the field overlying the channel briefly reduces the shear strength, producing a double peak. Late in the evolution, the shear field does not drop quite as low prior to eruption, nor does it increase quite as high after eruption, compared to the uniform-ζ case. These two features indicate that the filament channel and its overlying arcade neither expand quite as freely nor contract quite as much, when the shear is localized to the southern PIL.

Zoom In Zoom Out Reset image size

The idealized flux distribution used in the preceding examples has sufficed to show some of the capabilities of the STITCH model. Actual solar flux distributions can be much more complex, even fragmented, posing far greater challenges to eruptive event modeling of such regions and requiring far greater computational resources. To explore the potential usefulness of STITCH for investigations of this kind, we created a synthetic active region with a corrugated PIL by superposing 30 subsurface magnetic dipoles, all centered nominally at latitude ψ₀ = 225 and longitude ϕ₀ = 0°. We randomized their displacements away from this center, ∣ψ − ψ₀∣ ≤ 75 and ∣ϕ∣ ≤ 225, and their horizontal orientations away from due north, within ±45°. Each dipole was placed at a depth of 80 Mm below the surface and was given a tangential field strength of 6.7 G at the surface. Figures 8(a) and (b) show views of the resulting initial magnetic field with its corrugated PIL.

Zoom In Zoom Out Reset image size

The active region was energized using the STITCH profiles of tangential magnetic field injection shown in Figure 9. We adopted ζ as expressed in Equation (20), but here we used ϕ_w = 45° and ψ_w = 27°. As can be seen in the figure, STITCH flux injection was applied in both polarities of the active region, as well as in the surrounding background magnetic field. In this case, we used a smaller value of ζ₀ = 5 × 10⁵ km² s⁻¹, and the gradients in the radial magnetic field B_r were gentler owing to our constructing the active region from a set of well-submerged dipoles. Consequently, the resulting values of d B _s /dt were about an order of magnitude smaller than in our idealized cases (see Figure 2). We allowed the usual ramp-up interval of 1000 s duration in this simulation and then held the tangential field injection rates fixed thereafter.

Zoom In Zoom Out Reset image size

The elapsed time required to build a strongly sheared filament channel and induce it to erupt was about 2.5× larger for this case, i.e., some 15,000 s. Snapshots of the formation and eruption stages in the evolution are shown in Figures 8(c) and (d). The development is very similar qualitatively to that shown in our previous examples, with the filament channel highly localized to the PIL and generating a fast CME when it erupts. The key point here is that essentially no change in the STITCH driving procedure from the simple case was required for handling this more "realistic" PIL, whereas driving by small-scale twists, shear flows, or flux rope insertion would have required the tedious construction of drivers that followed the detailed geometry of the meandering PIL.

In this section, we briefly describe the utility of STITCH to event modeling of solar eruptions (for reviews see, e.g., Toriumi & Wang 2019; Jiang et al. 2022). A full discussion of the topic is beyond the scope of our paper and ambition. To frame the presentation, we note that event modeling attempts to replicate an observed solar eruption by taking one of three general approaches (for details see, e.g., Chintzoglou et al. 2019; Toriumi & Wang 2019): (1), (2) Energize some (arbitrarily) chosen initial configuration with the model, generating a single-snapshot end state that resembles the observed pre-eruptive solar structure; this modeled state may correspond only qualitatively to the observed state ("data-inspired"), or it may be based quantitatively on measurements of the observed state ("data-constrained"). (3) Evolve the observed solar structure within the model over an extended series of snapshots, constraining the model evolution so that it replicates the multistep evolution observed in the data at each step ("data-driven").

The usefulness of STITCH for energizing filament channels is demonstrated by the results presented in the preceding sections. It has a number of key advantages in comparison to some other methods commonly used in data-inspired or data-constrained models to generate the coronal structure preceding a solar eruptive event (for a review, see Patsourakos et al. 2020). These methods include flux emergence from below the surface (e.g., Manchester et al. 2004; Leake et al. 2013), surface shear flows (e.g., Lynch et al. 2008; Wyper et al. 2017), surface flux cancellation (e.g., Amari et al. 2010; Aulanier et al. 2010; Hassanin et al. 2022), and procedures for direct insertion of coronal flux ropes (e.g., Titov & Démoulin 1999; van Ballegooijen 2004; Fan 2005; Borovikov et al. 2017; Titov et al. 2021). One problem for these methods is matching the normal flux at the boundary of the numerical system to photospheric observations. This is particularly difficult for emergence, shear flows, and cancellation, because these processes generally change the flux distribution. Another problem is dealing with complex-geometry PILs, which can be especially challenging for insertion methods that require the specification of multiple different flux ropes carefully fitted to match the observed filament channel (e.g., Torok et al. 2018). Such intricate iterative procedures would be challenging to use in real time for space-weather predictive modeling. A third problem is that after an equilibrium has been calculated the eruption must be triggered without disturbing the original normal flux distribution. By construction, STITCH avoids all of these difficulties. Unlike flux rope insertion, the shear injection by STITCH is continuous; hence, it can be halted readily at any point to examine the pre-eruptive structure and possible trigger mechanisms near marginal stability. Note that STITCH could be run automatically as part of any 3D MHD heliospheric model. Furthermore, STITCH is so flexible and easily implemented that it can be used in combination with the other filament channel formation and energization mechanisms. For example, one could insert one or more equilibrium flux ropes where desired and then use STITCH to drive the system to eruption while preserving the matched surface boundary conditions. Beyond its foundation in a compelling physical model for filament channel formation and energization, helicity condensation, STITCH clearly has substantial practical utility for use in single-snapshot, data-inspired or data-constrained eruptive event modeling.

Multiple-snapshot, data-driven modeling of solar events is a much more challenging undertaking, requiring the model to be constrained at each step over a time sequence of observational measurements. Numerous approaches have been developed and investigations pursued with this aim in recent years (e.g., Mackay et al. 2011; Cheung & DeRosa 2012; Fisher et al. 2012; Kazachenko et al. 2014; Cheung et al. 2015; Lumme et al. 2017; Hayashi et al. 2018, 2022; Pomoell et al. 2019; Hoeksema et al. 2020; Toriumi et al. 2020). Unlike any single-snapshot method, including STITCH, data-driven models must incorporate a horizontal electric field E _s that replicates the evolution of the normal magnetic field component B _n at the Sun's surface,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{B}}}_{{\boldsymbol{n}}}}{\partial t}=-{{\boldsymbol{\nabla }}}_{s}{\boldsymbol{\times }}{{\boldsymbol{E}}}_{{\boldsymbol{s}}}.\end{eqnarray} \tag{ 21 }$$

In order to fully constrain the evolution of the surface magnetic field, observational data must be available and reliable for the rate of change of both B _n and B _s . Due to observational limitations, these boundary values are frequently underspecified, and a variety of methods have been developed to complement the measurements (see discussion in Jiang et al. 2022). We note that it is precisely the horizontal components of the magnetic field that are critical to the formation and energization of the filament channels that culminate in solar eruptions.

Among those models that do not attempt to constrain the evolution of the horizontal field B _s , one class (Cheung & DeRosa 2012; Cheung et al. 2015; Lumme et al. 2017; Pomoell et al. 2019) imposes a mathematical closure on the electric field E _s by supplementing the constraint on its curl (Equation (21)) with a specification for its divergence, via one of a hierarchy of prescribed forms:

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{s}\cdot {{\boldsymbol{E}}}_{{\boldsymbol{s}}}=\left\{\begin{array}{l}{C}_{0};\\ {C}_{1}{\boldsymbol{n}}\cdot {\boldsymbol{B}};\\ {C}_{2}{\boldsymbol{n}}\cdot \left({{\boldsymbol{\nabla }}}_{s}{\boldsymbol{\times }}{\boldsymbol{B}}\right).\end{array}\right.\end{eqnarray} \tag{ 22 }$$

The constant C_i is determined empirically by comparing the evolution of the observed and modeled systems. The STITCH term in our induction equation (Equation (5)) corresponds to the potential (curl-free) electric field

$$\begin{eqnarray}&&{{\boldsymbol{E}}}_{{\boldsymbol{s}}}^{{\boldsymbol{S}}}={{\boldsymbol{\nabla }}}_{s}\left(\zeta {B}_{n}\right),\end{eqnarray} \tag{ 23 }$$

whose divergence is

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{s}\cdot {{\boldsymbol{E}}}_{{\boldsymbol{s}}}^{{\boldsymbol{S}}}={{\boldsymbol{\nabla }}}_{s}^{2}\left(\zeta {B}_{n}\right).\end{eqnarray} \tag{ 24 }$$

In regions where ζ is uniform, this expression may be recast into the simple form

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}_{s}\cdot {{\boldsymbol{E}}}_{{\boldsymbol{s}}}^{{\boldsymbol{S}}}={C}_{3}{\boldsymbol{n}}\cdot \left({{\boldsymbol{\nabla }}}_{s}^{2}{\boldsymbol{B}}\right).\end{eqnarray} \tag{ 25 }$$

Mathematically, STITCH is similar to the next element in the Cheung hierarchy (Equation (22)) at one higher level in derivatives of B (Equation (25)). Mackay et al. (2014) estimated the solar value for ζ to be (6 ± 4) × 10² km² s⁻¹. It would be highly informative to see the results of an event modeling experiment that assumed Equation (25) with a value for C₃ in this ζ range. Notice that the STITCH mathematical form is motivated by the physics of helicity condensation, and the value of its parameter is constrained a priori by observations. Hence, a successful experiment could lead to a substantial advance in our broad physical understanding of filament channel formation and evolution on the Sun.

A second class of models in which the horizontal field B _s is unconstrained is surface flux transport simulation (e.g., Mackay et al. 2011, 2014, 2018; Hoeksema et al. 2020), in which prescribed physical processes of differential rotation, meridional flow, flux cancellation due to supergranular diffusion, and emergence of fresh magnetic flux from below are all included. The evolution of the normal field B _n in such models is more loosely constrained, on a rather long (typically one Carrington rotation) timescale, by periodically assimilating up-to-date maps of the normal surface field into the simulation to correct any accumulating errors. STITCH was introduced initially in collaborative studies with Mackay (2014, 2018; referred to therein as Statistically Averaged Helicity Condensation) that showed how adding our subgrid-scale model to his simulations improved the agreement between the observed and simulated magnetic maps. In particular, sufficiently fast helicity condensation was able to counter the shear induced by differential rotation along the high-latitude PILs of the polar crown.

Following on these studies, the STITCH formalism was used by Mackay and coworkers to predict the state of the corona during the 2017 August 21 solar eclipse (Mikic et al. 2018; Yeates et al. 2018). Those simulations included an electric field of the STITCH form (Equation (23)), in which ζ was replaced by a masking coefficient M employed to localize the induced B _s to the neighborhood of the PIL (Mikic et al. 2018, p. 919), similar to the procedure that we have described in detail here in Section 5 on localized helicity injection. The eclipse simulations (and others) described the masked electric field as corresponding to "emerging" transverse (horizontal) magnetic field locally within the filament channels (Yeates et al. 2018, p. 99). All of the simulations shown in this paper assume a reflecting condition (i.e., a zero value) on the normal component of the plasma velocity V at the inner radial boundary; consequently, flux emergence plays no role in the formation of our filament channels. Furthermore, we are not convinced that the highly complex photosphere−corona interaction accompanying flux emergence can be captured by a simple STITCH-like formalism. Rather, the STITCH term that we originally derived in Mackay et al. (2014) and employ here seeks to capture effectively the cumulative effects of surface vortical flows, induced horizontal shear field, and reconnection-driven transport of the shear field across the surface.

STITCH is our acronym for STatiscal InjecTion of Condensed Helicity, a subgrid-scale physical model for the formation and evolution of filament channels on the Sun. In this article, we have used STITCH to inject magnetic helicity and free energy into simulated coronae that evolve self-consistently according to the otherwise-familiar equations of magnetohydrodynamics. The only modification to standard MHD that is required to implement STITCH is the addition of a new, mathematically simple term to the induction equation (Faraday's law). Including the STITCH term amounts to adjusting the electric field, i.e., modifying Ohm's law in the coronal plasma, at the base of the atmosphere.

The mathematical simplicity of the STITCH term obscures, but simultaneously encapsulates, a great deal of physical complexity. As derived originally (Mackay et al. 2014, 2018), it represents the injection of helicity by vortical cells characterized by rotation rate ω_ℓ and radius ℓ, the formation of current sheets at the boundaries between adjacent cells, and efficient reconnection of the induced horizontal magnetic fields across the current sheets. The STITCH model distills these detailed processes of helicity condensation described by Antiochos (2013) into a single term with one parameter, ζ ≡ 〈ℓ² ω_ℓ 〉/2, plus an assumption about the height h over which the helicity is injected into the magnetic field. In practice, the simplest approach is to perform the injection into the bottom-most grid cell adjacent to the lower boundary, as has been done here and by others (Mackay et al. 2014, 2018; Lynch et al. 2021).

Implementing and using the STITCH term in MHD models are straightforward:

• 1.  
  Assign a value to the helicity injection parameter ζ.
• 2.  
  Evaluate the STITCH product ζ B_n at the bottom boundary, n = n₀, where B_n is the normal component of the magnetic field and n = x (Cartesian) or n = r (spherical).
• 3.  
  Difference this product along the two horizontal coordinates, (y, z) or (θ, ϕ), setting the sign according to the rules for the curl (see Equations (11) and (12)), to compute the tangential flux change rates.
• 4.  
  Divide the flux change rates by the appropriate vertical cell-face areas to calculate the corresponding tangential field change rates (see Equations (13) and (14)).
• 5.  
  Multiply these rates by the numerical time increment, Δt, finishing the calculation of the tangential field changes.

To our knowledge, the procedure above would be trivial to implement in any of the MHD codes currently used by the space science community. Indeed, it already has been implemented in one such model (Mikic et al. 2018; Yeates et al. 2018). Moreover, the helicity injection via STITCH may be tailored spatially to occur only within one polarity of the magnetic field, or only within one or more restricted region(s) of the bottom-boundary surface, or any combination of both. This is accomplished by introducing suitable spatial variations into the ζ parameter adopted in step 1, as we have done in this paper. The results for the tangential field changes will be mathematically smooth so long as the product ζ B_n goes to zero smoothly at the boundary of the tailored region, i.e., either ζ vanishes (by construction) or B_n vanishes (at a PIL).

Although the STITCH model is derived from the intricate physical processes inherent to the FHC model, we emphasize that its essence is to accumulate helicity at the large-scale boundaries of magnetic flux systems, in particular at PILs. Therefore, in general it represents the outcome of an inverse cascade of magnetic helicity, from small scales to large, due to any mechanism that injects helicity into the corona: vortical motions, shearing motions, flux emergence, flux cancellation, and so forth. Our results in this article demonstrate that STITCH can be used as a generic tool for coronal MHD modeling that forms filament channels and energizes magnetic field configurations. Here we have used it to evolve both highly idealized and more elaborately corrugated active region structures, and we have compared it successfully, both qualitatively and quantitatively, with a fully detailed helicity condensation calculation. In other work, we have used STITCH to perform full-Sun, sunspot-cycle scale modeling of filaments in a flux transport model (Mackay et al. 2014, 2018) and to energize a high-latitude filament channel leading to the eruption of a "stealthy" CME (Lynch et al. 2021). Others have employed the method to form filament channels in a global corona model used to predict solar eclipse observations (Mikic et al. 2018; Yeates et al. 2018).

The STITCH model possesses the following efficiencies and advantages compared to fully detailed MHD descriptions of filament channel formation and eruption:

• 1.  
  Straightforward implementation within existing MHD models.
• 2.  
  Easy application to realistic, complex PILs, as well as to localized regions along extended PILs, as is likely to be required for modeling observed events.
• 3.  
  Energizes the system while leaving the normal component of the magnetic field untouched, thereby greatly facilitating event modeling.
• 4.  
  Reduced computational expense (15× in our detailed comparison) relative to FHC because STITCH is amenable to sub-Alfvénic (flux injection) driving rather than requiring subsonic (vortical-flow) driving, allowing much faster formation and eruption of filament channels.
• 5.  
  Reduced computational expense (not exploited in our detailed comparison) relative to FHC because the small-scale vortical flows and current structures need not be resolved, allowing the use of coarser grids.
• 6.  
  Elimination of complex small-scale structure and dynamics, so that the large-scale topology and dynamics of the filament channel evolution can be ascertained more readily.

In summary, STITCH demonstrates substantial utility for physics-based modeling of filament channel formation and evolution culminating in CMEs and, potentially, great promise as an operational tool for future space-weather prediction capabilities. The work required to demonstrate and validate this promise is already underway.

Our work was sponsored by NASA's H-LWS, H-SR, and H-ISFM research programs. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Center for Climate Simulation (NCCS) at Goddard Space Flight Center. J.T.D. was supported by a NASA Postdoctoral Program fellowship administered by Universities Space Research Association at NASA Goddard Space Flight Center. J.T.D. acknowledges support from Grant No. 80NSSC21K1313. S.K.A. acknowledges support from Grant No. 80NSSC22K0892 to the University of Michigan.