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

The solar corona routinely exhibits explosive activity, in particular coronal mass ejections and their accompanying eruptive flares, that have global-scale consequences. These events and their smaller counterparts, coronal jets, originate in narrow, sinuous filament channels. The key processes that form and evolve the channels operate on still smaller spatial scales and much longer time scales, culminating in a vast separation of characteristic lengths and times that govern these explosive phenomena. In this article, we describe implementation and tests of an efficient subgrid-scale model for generating eruptive structures in magnetohydrodynamics (MHD) coronal simulations. STITCH -- STatistical InjecTion of Condensed Helicity -- is a physics-based, reduced representation of helicity condensation: a process wherein small-scale vortical surface convection forms ubiquitous current sheets, and pervasive reconnection across the sheets mediates an inverse cascade of magnetic helicity and free energy, thereby forming the filament channels. STITCH abstracts these complex processes into a single new term, in the MHD Ohm's law and induction equation, which directly injects tangential magnetic flux into the low corona. We show that this approach is in very good agreement with a full helicity-condensation calculation that treats all of the dynamics explicitly, while enabling substantial reductions in temporal duration especially, but also in spatial resolution. In addition, we illustrate the flexibility of STITCH at forming localized filament channels and at energizing complex surface flux distributions that have sinuous boundaries. STITCH is simple to implement and computationally efficient, making it a powerful new technique for event-based, data-driven modeling of solar eruptions.


INTRODUCTION
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 Mm (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 corona, in particular the flare ribbons, which can be as simple as the classic, linear two-ribbon structure or as complex as circular (Masson et al. 2009) or multi-ribbon (Wang et al. 2014) configurations.
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 §2. Our magnetohydrodynamics code, ARMS, and initial configuration are described in §3. We compare STITCH to the comprehensive helicity-condensation model in §4. Thereafter, we show examples of the flexibility of STITCH in localizing the filament channel ( §5) and in being used to energize more complex and realistic magnetic configurations ( §6). To conclude, we summarize our results and discuss the implications of using STITCH in simulations of coronal energy buildup in §7.

STITCH
In his exposition of the helicity condensation model, Antiochos (2013) described how an ensemble of cyclonic motions at the photosphere builds up magnetic shear in the corona. Current sheets form at the interfaces between pairs of neighboring cyclonic cells that have the same sense (clockwise or counter-clockwise) 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 cyclonic 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 (Zhao et al. 2015;Knizhnik et al. 2015Knizhnik et al. , 2017aKnizhnik et al. ,b, 2018Dahlin et al. 2019) in which numerous small-scale cyclonic 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 (Mackay et al. 2014(Mackay et al. , 2018, a subgrid-scale representation of the model was developed and employed. The essential assumptions that were made are: (1) the cyclonic 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 to the paper by Mackay et al. (2014, Eq. A1), the resulting subgrid-scale model is expressed as an additional term in the induction equation, 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 cyclonic flows. Specifically, 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. 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 STatistical InjecTion of Condensed Helicity, or STITCH.
The STITCH term has several notable properties. First and most fundamentally, by construction it preserves the divergence-free character of the magnetic field, 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, The change in the horizontal field component B s is 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 cyclonic cells (for details, see Mackay et al. 2014), 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 1 and s 2 , whereŝ 1 ×ŝ 2 =n. Integrating Equation 5 over differential cell areas dA s1 and dA s2 , respectively, we obtain the flux change rates 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 0 , by assuming that ζ falls to zero at the top of the first layer of grid cells at height n = n 1 . Integrating Equations 7 and 8 over n from n 0 to n 1 , they simplify to where ζB n | n0 denotes the value of ζB n evaluated at the surface, n = n 0 . These expressions integrate immediately to yield the flux changes where ∆ s1 (∆ s2 ) denotes a finite difference whose argument is evaluated at cell vertices along coordinate s 1 (s 2 ). The flux changes on the four faces of any grid cell involve pairwise differences of the argument (ζB n | n0 ) evaluated at the four vertices of the (s 1 , s 2 ) 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, where A s1 (A s2 ) is the area of the cell face perpendicular to the coordinate s 1 (s 2 ). Evaluating the (ζB n | n0 ) vertex values and multiplying them by the finite-difference time increment, ∆t, yields the required changes to the field components, ∆B s1 and ∆B s2 . its PILs are colored cyan. Magnetic field lines are colored according to whether they close entirely within the active region (red), between the active region and the background magnetic field (blue and magenta), or entirely within the background field (green), or are separatrix lines bounding these flux systems (gray/silver).

NUMERICAL SIMULATIONS
Our numerical calculations were performed with the Adaptively Refined Magnetohydrodynamics Solver (ARMS; DeVore & Antiochos 2008), which has been used extensively to model both CMEs/eruptive flares (see also Lynch et al. 2008Lynch et al. , 2009Karpen et al. 2012;Masson et al. 2013;Lynch et al. 2016Lynch et al. , 2021 and the formation of filament channels (Zhao et al. 2015;Knizhnik et al. 2015Knizhnik et al. , 2017aKnizhnik et al. ,b, 2018. Our fiducial calculations adopt the magnetic field configuration shown in Figure 1. It consists of an elliptical, bipolar active region centered at 22.5 • N 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, Fig. 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 6 K and pressure P s = 4 × 10 −1 dyn cm −2 . We solved the ideal MHD equations with an adiabatic temperature equation. The domain extents were r ∈ [1R s , 30R s ], θ ∈ [π/16, 15π/16], and φ ∈ [−π, +π].
For the full helicity condensation (FHC) calculation, the system was driven by 107 tessellated, cyclonic, B r -conserving cellular flows within the black-shaded, minority-polarity region as shown in Figure 2a. The flows were subsonic and sub-Alfvénic, attaining a maximum speed |V ⊥ | ≈ 50 km s −1 after an initial sinusoidal ramp-up interval 1,000 s in duration. Elsewhere on the surface, the magnetic field was line-tied at rest (V ⊥ = 0) for the entire duration of the simulation. 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.
For our first STITCH calculation described below, we assumed that ζ was uniform at ζ 0 = 1.4 × 10 16 cm 2 s −1 within the minority-polarity region, and zero everywhere outside of that region. The resulting profiles of injected horizontal magnetic field, whose maximum rate was about 1 G s −1 , are shown in Figure 2b. 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 1,000 s of elapsed time in the simulation. The energy-buildup phase was much shorter with STITCH, in this case being only about 6,000 s in total duration and requiring approximately four days of computational time.
We note that the stated rate of field injection and time interval suggest that a tangential field strength on the order of 6 kG should build up; however, this is not the case. 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.  Figure 3 shows snapshots of the filament-channel formation and eruption phases for the FHC and STITCH cases. 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 (Fig. 3a,c). 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 −1 , indicated by dark red shading in the plane of the sky (Fig. 3b,d). The kinetic energy is near its maximum value at the times shown for both cases.

COMPARING STITCH TO FULL HELICITY CONDENSATION
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 Eq. 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 Figure 3a,c. Figure 4 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.
In order to compare the macroscopic evolution most conveniently, the respective time axes are set to align the approximate eruption times in the two cases, t ≈ 90,000 s (FHC) and t ≈ 6,000 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 M 0 , and its later reduction in conjunction with eruption, are shown in Figure 4a. We find that the accumulated magnetic free energy is some 20% higher in the  STITCH case than in the FHC case, E M /E M 0 ≈ 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 M 0 ≈ 0.04 and 0.10, respectively, culminating in post-eruption energies E M /E M 0 ≈ 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 amounts of flux that is processed and energy that is released by the eruption. Second, the STITCH energy-release time scale is only ∆t ≈ 1,000 s, a factor of four less than the ∆t ≈ 4,000 s time scale 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 time scale, 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 4b. 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 time faster, as we also noted above.
The kinetic energy in the simulation, normalized to the initial magnetic energy, is shown in Figure 4c. This ratio 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 −4 , or M A ≈ 1%. The value jumps steeply during the explosive eruption phase to a few times 10 −2 , 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 φ | max , is shown in Figure 4d. 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 full helicity condensation (Dahlin et al. 2019). The major difference between the two is the highly compressed characteristic time scale 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 subAlfvénic, direct injection of horizontal magnetic flux into the corona. A more minor, however 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 time scales, 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 cyclonic flows in the full helicity-condensation simulation. The resulting eruptions displayed in Figure 3 (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. 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 §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 (Fig. 2b), 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 onesided 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 due to 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, Here φ is longitude, ψ is latitude, φ 0 = 0 • , ψ 0 = 22.5 • , φ w = 22.5 • , and ψ w = 11.25 • . The amplitude ζ 0 and the restriction to the minority-polarity region are the same as in the previous simulation. Figure 5 shows the resulting injection rates for B φ and B θ , for comparison with those in Figure 2b. 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 = 6,000 s and t = 7,000 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 = 6,200 s are shown in Figure 6 for the uniform-ζ (left) and localized-ζ (right) 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 (Fig. 6a); 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 (Fig. 6b). 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 Fig. 5b).
Side views of the azimuthal magnetic field B φ show how localizing the STITCH injection (Fig. 6d) 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 (Fig. 6c). Obvious features of the uniform-ζ case (Fig. 6c) 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 (Fig. 6d). 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-ζ (Fig. 6e) and localized-ζ (Fig. 6f) 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 also illustrates breakout-reconnection downflows along the flanks of this arcade in the localized-ζ case (Fig. 6f).
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 4, the time axes for the two cases [both STITCH] are identical.) The energy Figures 7a,b,c also show 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 (Figs. 7a,b) and beyond into the eruption phase for the kinetic energy (Fig. 7c). 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 ζ (Fig. 7b). The kinetic energy for that case peaks slightly above 50% of the uniform-ζ value (Fig. 7c). Peak strengths of the shear field attained early in the energization phase are essentially identical for the two cases (Fig. 7d). 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.

ENERGIZING COMPLEX FLUX DISTRIBUTIONS
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 data-driven 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 ψ 0 = 22.5 • and longitude φ 0 = 0 • . We randomized their displacements away from this center, |ψ − ψ 0 | ≤ 7.5 • and |φ| ≤ 22.5 • , 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. Figure 8a,b shows views of the resulting initial magnetic field with its corrugated PIL.
The active region was energized using the STITCH profiles of tangential magnetic field injection shown in Figure 9. We adopted ζ as expressed in Equation (15), 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 ζ 0 = 5 × 10 15 cm 2 s −1 , and the gradients in the radial magnetic field B r were gentler due to our constructing the active region from a set of well-submerged dipoles. Consequently, the resulting values of dB s /dt were about an order of magnitude smaller than in our idealized cases (cf. Fig. 5). We allowed the usual ramp-up interval of 1,000 s duration in this simulation, then held the tangential-field injection rates fixed thereafter. 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 Figure  8c,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 "realistic" PIL, whereas, driving by small-scale twists or by shear flows or flux-rope insertion would have required the tedious construction of drivers that followed the detailed geometry of the PIL. 7. DISCUSSION 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 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(Mackay et al. , 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 the 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 ω /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(Mackay et al. , 2018Lynch 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 0 , 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 tangentialfield 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 presently used by the space science community. 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).
The results presented in this paper demonstrate that STITCH has a number of key advantages in comparison to other commonly used methods for filament channel formation in 3D numerical modeling of flares and solar eruptive events (e.g., Patsourakos et al. 2020). These methods include flux emergence (e.g., Manchester et al. 2004;Leake et al. 2013), shear flows (e.g., Lynch et al. 2008;Wyper et al. 2017), flux cancellation (e.g., Amari et al. 2010;Aulanier et al. 2010), and several procedures for direct insertion of a flux rope (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 challenge is dealing with PILs of complex geometry, which can pose a severe challenge for flux rope insertion, at times, requiring the specification of multiple different flux ropes carefully fitted to match an observed filament channel (e.g., Török et al. 2018). It seems unlikely that such painstaking, iterative procedures could be used for space weather predictive models. A further problem is that after an equilibrium has been calculated, the eruption must be triggered without disturbing the original normal flux distribution. STITCH inherently takes care of all these problems. In contrast to flux rope insertion, STITCH shear injection is continuous and therefore can be readily halted at any point to examine pre-eruptive structure and trigger mechanisms near marginal stability. Note also that STITCH is so flexible and easily implemented that it can be used in combination with these other filament channel formation mechanisms. For example, one could insert some equilibrium flux rope, as desired, and then use STITCH to drive the system to erupt while still matching the photospheric boundary conditions. Although the STITCH model is derived from the intricate physical processes inherent to the full helicity condensation 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, it has been used to perform full-Sun, sunspot-cycle scale modeling of filaments in a flux-transport model (Mackay et al. 2014(Mackay et al. , 2018 and to energize a high-latitude filament channel leading to the eruption of a "stealthy" CME (Lynch et al. 2021).
The STITCH model possesses the following efficiencies and advantages compared to fully detailed MHD descriptions of filament-channel formation and eruption: • Straightforward implementation within existing MHD models • Easy application to realistic, complex PILs, as well as to localized regions along extended PILs as is likely to be required for performing data-driven modeling of observed events • Energizes system while leaving the normal component of the magnetic field untouched, thereby, greatly facilitating event modeling • Reduced computational expense (15× in our detailed comparison) because STITCH is amenable to subAlfvénic (flux-injection) driving rather than requiring subsonic (vortical-flow) driving, allowing faster formation and eruption of filament channels • Reduced computational expense (not exploited in our detailed comparison) because the small-scale vortical flows and current structures need not be resolved, allowing the use of coarser grids • 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 CMEs and, potentially, great promise as an operational tool for future space-weather prediction capabilities. Our work was sponsored by NASA's H-LWS, H-SR, and H-ISFM research programs and a grant of HEC computer resources to C.R.D. at NASA's Center for Climate Simulation. J.T.D. was supported by a NASA Postdoctoral Program fellowship administered by Universities Space Research Association at NASA Goddard Space Flight Center.