CALCULATING ENERGY STORAGE DUE TO TOPOLOGICAL CHANGES IN EMERGING ACTIVE REGION NOAA AR 11112

rmv_flick attempts to smooth the temporal structure of our mask arrays by better associating regions that last for just a single timestep. This happens when the saddle points in the tessellation algorithm slosh back and forth over the threshold value. Figure 3 shows an example of this problem. Region N2 pops into existence in the middle timestep, occupying some of the territory ascribed to region N1 in the first and third timestep. We see that this is not a "real" new region, but should instead be considered part of N1 for the entire time. Figure 3, Column 2, shows the result of running rmv_flick on these data.

Zoom In Zoom Out Reset image size

The rmv_flick algorithm works on a sliding, three timestep window (initial, middle, and final timestep), and in two steps. In the first step, we find all regions Θ_i that have non-zero flux only in the middle timestep, and find all mask pixels in the initial and final timestep which overlap Θ_i. We take each overlapping region Θ_A and relabel those pixels of Θ_i which are overlapped by Θ_A in the initial and final timesteps as Θ_A.

Because our source regions change shape and size (pixel count) over time, we will usually have leftover pixels after performing the above bulk relabeling. While leftover pixels remain, we relabel the pixel of Θ_i that currently borders the greatest number of pixels of a single neighbor to that neighbor's label. We have found that, when multiple regions overlap a flickering region, this method divides those pixels between the overlapping regions in decent proportion to their past and future "control" of flickering territory.

With rmv_flick we dealt with regions that existed for just a single timestep. We now deal with the conceptual inverse problem, where a region exists for some time and then suddenly changes names in the next timestep. Figure 4 demonstrates the problem, as the region labeled N5 in the top left panel (and previous, undepicted times) cycles through a series of names and partitions within a few timesteps.

Zoom In Zoom Out Reset image size

We again begin by considering the mask array in a sliding, three timestep window. We find all non-zero flux regions within the window and determine which regions exist in which timesteps. For every region that disappears in the middle timestep, we project its area from the first timestep into the middle and last timesteps. We then find all new regions, in both the middle and final times, that overlap with that projected area. Finally, we relabel any new region whose centroid lays within this projected area. Figure 4, right column, shows the effect of running rmv_vanish on the data.

rmv_flick and rmv_vanish each use a sliding, three timestep window, so that it takes many repetitions of each for updated information propagate from the beginning to the end of the mask array. We must therefore repeat each of the algorithms many times. At the same time, rmv_flick acts upon labels that last for just a single timestep; eventually, all such problems will be fixed. Taking these processes into account, we find that the mask arrays best approximate the photospheric field evolution when we switch back and forth between the two algorithms at the beginning of the process, and run just rmv_vanish many times at the end of the process. For our NOAA AR 11112 data, we ran rmv_flick a total of 35 repetitions, and rmv_vanish a total of 1058 repetitions.

Together, these two algorithms accomplish about nine-tenths of the work in creating a consistent set of tessellated masks. Part of the remaining tenth may be accounted for by a boundary-shift algorithm, described below, between abutting regions of the same polarity. Even this does not accurately represent to data, and the final changes to the mask array required for a consistent time series must be done by hand. These changes again mostly involve the placement of the boundary between abutting regions of the same polarity, caused by a failure to distinguish between old flux and actively emerging flux. These boundaries are adjusted manually until each region's flux evolves smoothly.

Figure 5(a) shows the vertical flux in each (high-flux) region. This includes editing of the masks by both the automatic algorithms and by hand. Note that, even after a set of consistent masks have been created, the flux within each region is still rather noisy. To counteract this, we smooth our data with a 7 hr (13 timestep) boxcar function. The result is shown in Figure 5(b).

Zoom In Zoom Out Reset image size

In order to isolate the different varieties of flux change, we have developed an algorithm which estimates the change in each region's flux due only to boundary displacements between adjacent, like-signed regions. These changes necessarily come in pairs: what one region loses in a boundary shift, another gains.

Our algorithm works as follows. Consider a set of like-signed poles {P}, and split {P} into submerging and emerging sets, {P_↓} and {P_↑}, respectively. Each pole has an associated flux change, Δⁱψ ≡ ψ^{i + 1} − ψⁱ, with ψⁱ the flux at time i given by Equation (2). We iteratively find the pole with smallest unsigned flux change, P_s, and the pole with closest centroid $\mathbf {x}$ and opposite sense flux change, P_c. The flux change between these two poles is canceled, Δψ_c → Δψ_c − Δψ_s and Δψ_s → 0, and the cancellation is recorded in a change-in-connectivity matrix: $\Delta ^i\mathbb {M}_{c,s}=-\Delta ^i\mathbb {M}_{s,c} = \Delta \psi _s$. This process is repeated until no more connections can be made.

$\Delta ^i\mathbb {M}$ is antisymmetric and the {j, k}th element is the flux change of the edge joining j and k. Each row j describes the flux change for a given pole, and each column k in a row records how much of pole j's total flux change is in partnership with k. If all of pole j's flux change is due to boundary shifts, then j's total flux change is given by summation along the jth row of $\Delta ^i\mathbb {M}$:

$$\begin{equation} \psi ^{i+1}_j = \psi ^i_j + \sum _{k}\Delta ^i\mathbb {M}_{j,k}. \end{equation} \tag{ 3 }$$

In this way, we have paired as much shrinking flux with like-signed increasing flux as possible, but have only used the total flux and centroid location of each region.

We now refer back to the photospheric mask for the last time. For each region in the mask array we find every like-signed region with which it shares a border: call this subset of poles {b}. To allow for regions separated by a few pixels to share a boundary,² we pad each region by 5 pixels and look for overlap with other regions. The value of $\Delta ^i\mathbb {M}$ for each pair of such regions is then the amount of flux they exchanged through a shift in the mask boundary.

We have chosen to distribute the flux change and then restrict to the subset {b}, rather than the reverse, in order to avoid the ill-defined problem of assigning shifts when one region shares boundaries with multiple other regions, which in turn share boundaries with multiple other regions, and so on. We see this, for instance, in the triple boundary between regions P4, P15, and P19. Qualitatively, any discrepancy is small given that we connect regions in a nearest-to-farthest order, and nearby regions tend to share boundaries.

The sum of all such changes for a given pole gives that pole's total flux change due to boundary shifts. All other change must be ascribed to submergence or emergence through the photosphere. We may describe this mathematically as

$$\begin{equation} \psi ^{i+1}_j = \psi ^i_j + \sum _{b}\Delta ^i\mathbb {M}_{j,b} + \sum _k\Delta ^i\mathbb {S}_{j,k}. \end{equation} \tag{ 4 }$$

The middle term on the right-hand side (rhs) is the pole's boundary change, determined in this section; the final term quantifies the pole's flux change due to submergence and emergence through the photosphere, which we now quantify.

Our algorithm for determining the flux change due to submergence and emergence rests on two assumptions. First, pairs of poles submerge and emerge together, and second, a single region can either consist of submerging or emerging flux, but not both. These assumptions naturally break our connectivity-change graph into two disconnected subgraphs: one consisting of submerging vertices, the other emerging. Each subgraph is composed of both positive and negative polarity poles, as flux connects only regions of opposite polarity.

In determining the flux change in each domain, we use the same algorithm as for boundary shifts, only with the roles of polarity and sense of change reversed: these domains connect opposite polarity poles with same sense flux change, and represent pairs of poles submerging or emerging together. The graph generated in this manner reflects our physical intuition about how these systems should interact. Namely, we expect poles to be more connected to those poles closest to them, rather than those farther away. We also expect that poles with little flux change should not have that small change spread between a large number of partners.

We begin by finding the flux difference for each pole between times i and i + 1, from which we subtract the change due to boundary shifts, found previously. The remainder is each pole's flux change budget, which must be paired with other poles. To accomplish this, we again iteratively find the least changing pole and cancel its flux change with its closest appropriately signed and sensed neighbor. This cancellation is the weight of each edge and is recorded in the source-change matrix, $\Delta ^i\mathbb {S}$. Just like the boundary-shift matrix, $\Delta ^i\mathbb {S}$ is defined such that, if all of a pole P_j's flux change between times i and i + 1 is due to real photospheric changes, then $\psi _j^{i+1} =\psi _j^i + \sum _k\Delta ^i\mathbb {S}_{j,k}$; that is, summation over all columns of a given row returns the total flux change for that row's pole.

While the combination of our algorithms dealing with boundary shifts and emergence do quite well in most circumstances, occasionally they produce unphysical connections. Consider the timestep shown in Figure 7. Here, the emerging positive regions (P22, P29, P48) each have some flux change, and our knowledge of the location of emergence says that this change is most likely a combination of boundary shifts and emergence for P48 and P29, as both N87 and N58 are also actively emerging, and primarily boundary shifts for P22. However, if the boundary shift between P22 and P29 is larger than the emergence of P48, then our boundary algorithm ascribes the boundary shift to P29–P48, leaving the emergence mostly between P22 and N87, which is clearly unphysical. For this reason, we allow for "preferred connections," which accounts for emergence between a subset of poles first, and then proceeds to calculate boundary shifts and submergence and emergence for all poles. For this region, we first pair emergence between regions P22, N55, N58, and N87 for the first 3 days (77 timesteps), and then between regions P48, N55, N58, and N87 during the remaining time. As with the generation of consistent mask structures in the previous section, we ultimately wish to automate this entire process. However, until we develop algorithms that accurately capture the physical evolution of emerging regions, we will rely on a combination of automatic procedures and ad hoc prescriptions to most accurately model these systems.

Zoom In Zoom Out Reset image size

Figure 7 shows a graphical representation of the changing domain fluxes due to submergence and emergence between 12:04 and 12:35 UT on October 16th, 7 hr before the flare. Black lines show submerging domains, white lines emerging. Table 1 shows the flux change in each domain. Regions P4, P15, P22, P29, and N55 have flux change (for these timesteps) due solely to boundary displacements with neighboring regions.

Table 1. Non-zero Changes in Domain Fluxes between 12:04 and 12:35 on 2010 October 16

Domain	Sense	$\Delta ^{12:04}\mathbb {S}$
		(1016 Mx)
P19–N1	Submerge	245.764
P19–N38	Submerge	30.6279
P19–N49	Submerge	145.569
P19–N58	Submerge	257.030
P48–N87	Emerge	849.874

Download table as:  ASCIITypeset image

As discussed above, we begin pairing P48, which increased by 1506.4 × 10¹⁶ Mx. The closest (and only) negative pole with the correct sense of flux change is N87, which increased by 849.9 × 10¹⁶ Mx. These two poles are paired, the domain flux change is set to 849.9 × 10¹⁶ Mx, and P48's flux budget is reduced to 565.5 × 10¹⁶ Mx. There are no more available connections between our preferred poles, so we next account for boundary shifts. After doing so, the least changing region is N38, which loses 30.6 × 10¹⁶ Mx. The closest positive pole with the correct sense of flux change is P19, which lost 1528.9 × 10¹⁶ Mx. These two poles are paired, the domain flux change is set to 30.6 × 10¹⁶ Mx, and P19's flux budget is reduced to 1498.3 × 10¹⁶ Mx. We again select the least changing region (N49), and the process continues until no more pairings can be made. The final results for this timestep are shown in Table 1.

The pairings represent a compromise between reasonable inference and mathematical necessity. Due to their proximity and common increase, it is reasonable to assume that P48–N87 are feet of the same emerging flux tubes. P19, on the other hand, is relatively isolated but has steadily decreasing flux. It is necessary to pair it with decreasing negative poles, which turn out to be arranged along the periphery of the old polarity. Being separated by ∼100'', these regions cannot be "submerging" in any real sense. However, we believe this decrease represents steady cancellation with the small-scale field surrounding these old flux regions. This process should be studied in detail in another investigation, but for the moment we note that it is well represented by a reduction in, for instance, the P19–N49 domain flux.

We repeat this process for each pair of consecutive timesteps. Figure 8 provides a representative example, showing the cumulative domain flux changes (designated by the symbol $\hbox{\textsf {c}}\hspace{-7.5pt}\sum$) due to the emergence, and slight amount of submergence near the end of our analysis, of P29, one element of the newly emerging flux. We see here that it emerges primarily with region N58, part of the newly emerging flux.

Zoom In Zoom Out Reset image size

The effectiveness of the combined algorithms can be appreciated in Figure 6, showing the three components of the diffuse western region. The blue line illustrates the effectiveness of our algorithms at capturing the system's behavior due to submergence and emergence. As the field shears, the boundaries between the three subregions have a tendency to discontinuously jump, resulting in large variations in the flux in each region, accounted for by the boundary algorithm. However, the combined flux of the three subregions, shown in black, instead displays a fairly steady decrease: this whole region is submerging. The blue line shows our reconstruction of this submergence using the source-change matrix at each time. Each pole has an initial flux value, to which we add the elements of $\Delta ^i\mathbb {S}$, so that each pole's flux at time i is given by $\psi _j^i = \psi _j^0 + \sum _{l=0}^{i-1}\sum _k \Delta ^{l}\mathbb {S}_{j,k}$. This measure completely discounts the (artificial) boundary shifts between photospherically adjacent regions. We then sum the flux in each pole thus reconstructed, and this is our estimate for the submergence of diffuse western region. We note that the blue and black lines follow each other reasonably well, indicating that our algorithms accurately capture the submerging trend of this region.

In the second section of our analysis, we employ the MCC model (Longcope 2001, 1996) to use the amount of flux change in each domain to calculate the amount of free magnetic energy stored in our system as it evolves away from a potential configuration. To that end, we must define a set of topological elements. These definitions, briefly summarized below, are described more fully in Longcope & Klapper (2002), Longcope (2005), and references therein. The general idea is to model distributed photospheric sources as point sources. The potential field generated by these sources will have null points, $\mathbf {x}_\alpha$, where the magnetic vector field vanishes: $\mathbf {B}(\mathbf {x}_\alpha) = 0$. We may perform a linear expansion of the field about these points (Parnell et al. 1996; Longcope 2005, Section 2.4)

$$\begin{equation} \mathbf {B}(\mathbf {x}_\alpha +\delta \mathbf {x})\simeq \mathbb {J}^\alpha \cdot \delta \mathbf {x}, \end{equation} \tag{ 5 }$$

where the Jacobian matrix $\mathbb {J}^\alpha _{ij}\equiv \partial B_i/\partial x_j$ is both symmetric and traceless because $\nabla \times \mathbf {B}=\nabla \cdot \mathbf {B}=0$. As such, $\mathbb {J}^\alpha$ has three orthogonal eigenvectors with three real eigenvalues: the eigenvectors of the two like-signed eigenvalues define a plane (the fan); the eigenvector of the opposite-signed eigenvalue defines a line orthogonal to this plane (the spine). From this foundation the rest of the topological description follows. We define a pole, the photospheric point source, located at the flux-weighted centroid of a mask. Nulls are the zeros of the potential magnetic field, and spines, the field lines connecting a pole to a null along the single eigenvector. A fan (synonym separatrix) is the set of field lines ending in a null's like-signed eigenvector plane. A domain is a simply connected space filled with field lines connecting a given pair of poles. Finally, a separator is the field line connecting two nulls, formed at the intersection of their respective fan surfaces. In a nonpotential field, the separator may broaden into a two-dimensional (2D) ribbon. An example illustrating some of these elements is shown in Figure 9. Coronal current sheets form along the separators of the field.

Zoom In Zoom Out Reset image size

We determine the locations of all nulls in two steps. First, we find those lying in the photospheric plane via the Newton–Raphson method of Barnes et al. (2005). Second, we check the Euler characteristics given in Longcope & Klapper (2002), which provide both 2D and 3D relations between the number of nulls of each type and the number of sources of each polarity. When the 3D characteristic is not satisfied, at least one coronal null is missing. When then supply a "by-eye" list of initial (x, y, z > 0) guesses to the same Newton–Raphson root finding method until coronal nulls have been found and the Euler characteristics are satisfied. Other more automated coronal null detection methods have been developed (Barnes 2007), but the ad hoc method described here has worked well in this case.

After determining the location of all nulls, we located separators via the method of Barnes et al. (2005), with one modification. Those authors only treat photospheric nulls, whose separators lay within separatrices above the photosphere. For coronal nulls, we search for separators in all 2π directions within the null's fan surface.

Our model photosphere consists of a plane (z = 0) of isolated sources surround by regions where the normal component of the magnetic field is zero. This boundary condition may be satisfied using the method of images, as in Longcope & Klapper (2002, Section 6). With the normal field reflectionally symmetric in z, we have B_z(x, y, −z) = −B_z(x, y, z). Whenever we introduce current in the corona along a separator, we must introduce its reflection in the mirror corona to maintain this symmetry. In this way, our separators form closed current loops.

Once we have created the set of matrices describing domain flux changes due to submergence or emergence, we can create similar matrices describing the flux change due to a changing potential field. We calculate the potential field connectivity at each timestep via the Monte Carlo method of Barnes et al. (2005, Section 3.1). This produces a set of connectivity matrices $\lbrace \mathbb {P}\rbrace$, where $\mathbb {P}^i_{j,k}$ is the flux connecting sources j and k at time i. From these matrices we calculate the change in connectivity due to the changing potential field:

$$\begin{equation} \Delta ^i\mathbb {P} \equiv \mathbb {P}^{i+1} - \mathbb {P}^i. \end{equation} \tag{ 6 }$$

Summation along rows of the connectivity matrix returns the total flux of the corresponding pole; e.g., $\psi _j^i = \sum _k\mathbb {P}_{j,k}^i$. We define $\Delta ^i\mathbb {P}$ to be antisymmetric, as with the matrices $\Delta ^i\mathbb {S}$ and $\Delta ^i\mathbb {M}$ before.

The MCC model derives an energy from the discrepancy between the actual domain fluxes $\mathbb {F}^i$ at time t_i and the potential fluxes then, $\mathbb {P}^i$. In the absence of flux emergence or submergence, the actual fluxes are fixed using the potential field domain fluxes at some initial time: $\mathbb {F}^i = \mathbb {P}^0$, as in Kazachenko et al. (2010). Any flux emergence or submergence through the photosphere, quantified as $\Delta ^i\mathbb {S}$ in the previous section, modifies these actual domain fluxes, so that

$$\begin{equation} \mathbb {F}^i = \mathbb {P}^0 + \sum _{j = 0}^{i-1}\Delta ^i\mathbb {S}. \end{equation} \tag{ 7 }$$

Because we assume a field, initially potential, whose fluxes are fixed under future evolution, we determine how far removed the field is from a potential field configuration by answering the question, "What flux must be added to a potential field domain D at time 0 to get a potential field domain at time i?" This question is answered by

$$\begin{equation} \mathbb {P}_D^i = \mathbb {P}_D^0 + \sum _{j=0}^{i-1}\Delta ^j \mathbb {S}_D + \sum _{j=0}^{i-1} \Delta ^j \mathbb {R}_D. \end{equation} \tag{ 8 }$$

On the rhs, the first term is the potential field at the initial time. The second is the total flux change through the photosphere. The final term is flux change due to coronal redistribution, which must be achieved by modifying the connectivity matrix. Together, the first two terms give the actual domain flux $\mathbb {F}^i_D$ at time i. All flux changes from coronal reshuffling are then given by

$$\begin{equation} \sum _{j=0}^{i-1} \mathbb {R}_D^j = \mathbb {P}_D^i - \mathbb {F}_D^i. \end{equation} \tag{ 9 }$$

Equation (9) holds for all domains, of course, and we represent this as a matrix equation by dropping the domain subscript D. $\Delta ^i\mathbb {P}$ is antisymmetric, and we have defined $\Delta ^i\mathbb {S}$ to be antisymmetric, so $\Delta ^i\mathbb {R}$ is also antisymmetric. While the antisymmetry of these matrices carries no physical information, we will later show that it does endow $\Delta ^i\mathbb {R}$ with a nice mathematical property.

Note that the difference in Equation (9) between the actual domain fluxes and the potential field fluxes depends on the choice of the initial time t₀. When dealing with actual data, we must pick an initial time to apply the MCC flux constraint when we believe the field is in a potential configuration in order to find a meaningful difference relative to a later potential configuration. Ideally, we would track a single region from its inception up to an event which redistributes the coronal flux, as in a flare. While NOAA AR 11112 does not fit this criteria, it likely matters little in this case. As is apparent from our supporting media, and also as we will show below, NOAA AR 11112 has a very stable flux configuration for ∼40 hr prior to flux emergence. This newly emerged flux, fixed in set domains as the photospheric field continues to evolve, rapidly diverges from a potential field configuration at later times.

So far this discussion only involves the sources themselves and their interconnections. Let us now introduce a set of separators {σ}. Separators are field lines that run from a null point of one type to a null of the opposite type. In the MCC model, they are the sites of current sheets within the corona (Longcope 2005; Priest & Titov 1996). Because current may only flow in loops, the separators themselves must close along some path in the mirror corona. We can then speak of separator fluxes: the flux of domains linked by some separator σ with its closure. We will discuss closures and explain the concept of linking in detail below, but the general idea is that field lines in a linked domain cannot reconnect to another domain except by passing through the separator.

Let ψⁱ_σ be the flux linked by σ at time i in the actual field, and

$$\begin{equation} \psi _\sigma ^{(v)i} = \sum _D\mathbb {P}_D^i = \sum _{\lbrace (j,k)\rbrace }\mathbb {P}_{(j,k)\in D}^i \end{equation} \tag{ 10 }$$

be flux in the potential field domains D = {(j, k)} linked by σ. For a given closure, a separator always links the same set of domains throughout time. The fluxes in those domains will generally change over time, however. When a domain's flux goes to zero, the domain no longer exists.

The FCE assumption of the MCC model states that domain fluxes are fixed as the field evolves. Minimization of the field's energy subject to these constraints shows that the coronal field will be current free except along the separator, where a current ribbon forms (Longcope 2001). This occurs whenever ψⁱ_σ ≠ ψ_σ^(v)i. We therefore define

$$\begin{equation} \psi _\sigma ^{({\rm cr})i} \equiv \Delta \psi _\sigma ^i = \psi _\sigma ^i - \psi _\sigma ^{(v)i} = \sum _D\mathbb {F}^i_D - \sum _D\mathbb {P}^i_D, \end{equation} \tag{ 11 }$$

which, from Equation (9), is

$$\begin{equation} \Delta \psi ^i_\sigma = - \sum _D\sum _{j=0}^{i-1}\Delta ^j\mathbb {R}_D. \end{equation} \tag{ 12 }$$

The double sum over elements of the redistribution matrix $\mathbb {R}$ is the difference between the flux linked by the separator in the actual and potential fields at time i. It is also the self-flux generated by current flowing along the separator, and therefore, the flux that must be added to the real field to relax the FCE constraint of MCC: $\psi _\sigma ^{({\rm cr})i} +\sum _{D,j}\Delta ^j\mathbb {R}_D = 0$.

Our method for calculating energies therefore rests on two assumptions. The first is that at some initial time 0 < i we assume that the domain fluxes are given by the potential field's domain fluxes, $\psi ^0_{(j,k)} = \mathbb {P}^0_{(j,k)}$, so that all the separator fluxes ψ^(cr)0_σ = 0. The second assumption is that the evolution of the current ribbon's self-flux depends on the evolution of the both the actual and potential field configurations:

$$\begin{eqnarray} \Delta ^i\psi _\sigma ^{({\rm cr})} = \Delta ^i\bigl (\Delta \psi _\sigma ^i \bigr) & = \Delta \psi _\sigma ^{i+1} - \Delta \psi _\sigma ^{i} \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} & = \Delta ^i\psi _\sigma - \Delta ^i\psi _\sigma ^{(v)} \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} & = -\sum _D \Delta ^i \mathbb {R}_D. \end{eqnarray} \tag{ 15 }$$

Summed over time, this is just the restatement of Equation (11) itself, namely, that the current ribbon's flux at some time i is the sum of all changes in each linked domain's flux over the life of the separator.

When ψ^(cr)i_σ is small we can estimate the properties of the current ribbon, as in Longcope & Magara (2004). For a separator of length L carrying current I, they determine the self-flux to be

$$\begin{equation} \psi ^{({\rm cr})i}_\sigma = \frac{{I L}}{4 \pi }\ln \Biggl (\frac{e I^*}{\vert I\vert }\Biggr), \end{equation} \tag{ 16 }$$

where I* is a measure of magnetic shear in the separator's vicinity. From this the excess energy of the MCC field relative to the potential field, i.e., the free energy, is (Longcope & Magara 2004; Longcope 2001)

$$\begin{equation} \Delta W_{{\rm MCC}} = \frac{{1}}{4\pi }\int _{\Psi _{{\rm potl}}}^{\Psi }I d\Psi = \frac{{L I^2}}{32 \pi ^2}\ln \left(\frac{{\sqrt{e} I^*}}{\vert I\vert }\right). \end{equation} \tag{ 17 }$$

Note that Equations (16) and (17) only take into account energy due to the self-inductance of each current loop. In general, we must include the effect of mutual inductance, though for physically separated loops we expect self-inductance to dominate the total energy.

The FCE assumption (11) requires that ψⁱ_σ = ψ^(v)i_σ + ψ^(cr)i_σ, which allows us to solve Equation (16) for the current in the current ribbon as a function of Δψⁱ_σ:

$$\begin{equation} I\big(\Delta \psi _\sigma ^i\big) = I^*\Lambda ^{-1}\big (4\pi \Delta \psi _\sigma ^i\big/LI^*\big), \end{equation} \tag{ 18 }$$

where Λ⁻¹(x) is the inverse of the function Λ(x) ≡ xln (e/|x|). Equation (18) provides a method for determining the current residing in, and hence the energy stored by, each separator: it tells us the current required to change the domain flux enclosed by a separator from that of a potential field to some other value. In the present case that value is given by the non-zero elements of $\Delta ^i\mathbb {R}$. A reverse, cumulative sum over the time index i for a given flux domain P_j–P_k gives the time history of departure from a potential field configuration.

As is implicit in Equation (12), each separator may link more than one domain: the appropriate ψ^(cr)i_σ for a separator is the sum of all changes in domain fluxes for each domain linked by that separator. Since all field lines in a given domain are topologically equivalent,³ we may establish separator linkage using a single representative field line and the separator itself. After each of these open curves is closed, their linkage is found using the Gauss linking number (Berger & Field 1984), L₁₂, for each field-line–separator pair:

$$\begin{eqnarray} &&L_{12} = L_{21} = \frac{{1}}{4\pi }\oint _{\ell _1}\oint _{\ell _2}\frac{{\mathbf {r}_1 - \mathbf {r}_2}}{\vert \mathbf {r}_1 - \mathbf {r}_2\vert ^3}\cdot (d\mathbf {r}_1\times d\mathbf {r}_2). \end{eqnarray} \tag{ 19 }$$

The linking number not only determines whether a separator links a field line (L₁₂ ≠ 0), but also the sense in which it does so (L₁₂ = ±n, with n an integer): to calculate the correct ψ^(cr)i_σ, you must sum over all linked domains, multiplied by their respective linking numbers.

In order to use Gauss' linking formula, we must have two closed curves, whereas our curves (usually⁴) begin and end at two separate sources (nulls for separators) in the photosphere, never penetrating beneath. So, in order to apply the linking formula in this case we must add to each coronal curve some curve in the mirror corona to form a complete loop. There are, in general, many ways to form the closure for each curve, and not all of them will lead to the same linking number. One requirement is that the separator be closed above the field line's closure: any closure below the field line will always give a linking number of zero. Therefore, we have chosen to close the separator with a straight line between the footpoints in the photosphere, and the field lines with a rectangular path formed by following the footpoints down and connecting them with a straight line in the z = −0.5 Mm plane. Another method would be to close the separator, say, with its projection into the z = 0 (photospheric) plane. Either method is valid, so long as the same method is used for all field lines and separators.

Another requirement for using the linking number to find linked domains is that we follow field lines and separators in a consistent direction. In general, changing the direction of one of the integrals takes $L_{12}\longrightarrow \bar{L}_{12} = -L_{12}$. For definiteness, we always trace separators from negative (A-type) nulls to positive (B-type) nulls above the photosphere, and thus close separators from B→A in the photosphere.

Because the reconnection matrix $\Delta ^i\mathbb {R}$ is antisymmetric, this same care is unnecessary when tracing the field lines, at least for the purpose of the ψ^(cr)i_σ calculation. Let us designate by Δⁱ(Δψ_{P1, N2}) the difference in flux change for the real versus potential field between times i and i + 1 (as in Equation (13)), when we trace a field line from P1 to N2 in the corona, and close the field line from N2→P1 below the photosphere. Then,

$$\begin{eqnarray} \Delta ^i (\Delta \psi _{P1,N2}) & = \Delta ^i\mathbb {R}_{1,2} \times L_{12}. \end{eqnarray} \tag{ 20 }$$

Now trace the field line from N2→P1 in the corona, with the appropriate closure for z < 0. We have

$$\begin{eqnarray} \Delta ^i (\Delta \psi _{N2,P1}) & = \Delta ^i\mathbb {R}_{2,1} \times \bar{L}_{12}\nonumber \\ & = (-\Delta ^i\mathbb {R}_{1,2} )\times (-L_{12})\nonumber \\ & = \Delta ^i (\Delta \psi _{P1,N2}). \end{eqnarray} \tag{ 21 }$$

We have therefore been justified in simply writing Δⁱψ^(cr)_σ, with no designation of tracing a field line from positive to negative or negative to positive in the corona, because Δⁱψ^(cr)_σ is just a summation of multiple such domains. This holds provided we always trace the separator in the same direction.

Once determined, the linking number for each separator/field line pair provides the correct method for calculating each separator's total ψ^(cr)i_σ, including all linked domains, and hence the energy stored by each separator.

We must account for instances where two photospheric regions merge together or split apart, as illustrated in Figure 10. Therein, we focus on four poles—P₁, P₂, N₁, and N₂—with four flux domains—ψ_{1, 1}, ψ_{1, 2}, ψ_{2, 1}, and ψ_{2, 2}—and two separators, S₁ and S₂, which link the domains in some fashion. Suppose two of the poles are each connected to the same null via its spines. If the flux of one pole is zero before (after) a split (merger), then the null between them does not exist at that time, and the flux in the domains linking the sources goes to zero. Now, if either S₁ or S₂ is attached to that null, then the separator will arise or cease to exist with that null, and the flux linked by it is appropriate. Then again, if these separators are attached to some other set of nulls and happen to link these flux domains, our calculated linked flux is still appropriate. Say flux domain ψ_{2, 2}, for instance, either folds into or breaks out from domain ψ_{2, 1}. The actual break or merger constitutes a boundary shift (with N₂ having zero flux either before or after the event) and so does not contribute to the flux linked by either separator. Any additional flux change on top of this is either appropriately picked up by S₁ or S₂, but not both.

Zoom In Zoom Out Reset image size

A second complication arises when we encounter a coronal null, as we see in the topology of NOAA AR 11112 in Figure 9 near local tangent plane coordinates (20,30) Mm. If a coronal null has any separators connected to it then it must have at least two: one carrying current up, the other down. If many separators connect to the coronal null, then some must carry current up, and others down, but all balanced according to Kirchhoff's rules. In order to calculate the flux within a separator-induced domain, we must close the separator somehow. Above, we addressed this by prescribing a photospheric closure for separators. Now, however, we must first stitch pairs of separators together in order to get back to the photosphere. If we have n separators connected to the coronal null, we require n − 1 separator pairs, or "isolating loops" in the terminology of Longcope (2001). We might assume that the currents flowing along the n − 1 loops would distribute to minimize the energy due to self- and mutual inductances. This problem remains to be solved, and for the moment we simply designate one separator linked to a coronal null as the "shared" separator. We have found that the energies we calculate do not depend heavily on the choice of common leg, laying within ≈20% of each other.

As shown in Figure 9, we have a single coronal null, B15, with spine lines leading to P22 and P19. Four separators (in red) connect to the coronal null: from null A09 between N58 and N67, A05 between N67 and N49, A08 between N49 and N38, and A07 between N38 and N1. Again, because all current flowing from the photospheric nulls up to the coronal null must be balanced by current flowing back down to the photosphere, we must combine these four separators into three, all of which share one leg; it does not matter which one. We have chosen the A05–B15 leg to be common to all three separators. As we will see below, the A05–B15–A08 separator generates the greatest amount of energy in this system. This is expected, given the high horizontal shear in the photospheric field between P22 and N49, and the greatly flux-deficient domain between these two poles, relative to the potential field configuration at the time of the flare.

We may summarize the previous several sections in the following way. In order to perform an energy calculation, we need a set of poles $\lbrace P\rbrace ^{t_f}$ at a given time final time t_f; a set of nulls at the same time $\lbrace N\rbrace ^{t_f}$; a set of separators at the same time $\lbrace \sigma \rbrace ^{t_f}$; a time history of the reconnective changes in domain fluxes from the final time, backward to a previous time t_i where we may assume the field was potential, $\lbrace \Delta ^i \mathbb {R}\rbrace$; and a set of domains $\lbrace D_\sigma \rbrace ^{t_f}$ linked by each separator.

Having all these pieces, we then calculate the total domain flux difference in each domain at time t_f relative to the initial presumed potential field. For each separator, we find all domains linked by that separator. Next, each linked separator's time history of reconnective flux changes is summed. This total flux change for each linked domain is then multiplied by the Gauss linking number for that domain/separator pair. Finally, the total domain changes are summed for each linked domain to fix the total flux difference of the separator at the final time t_f, relative to the potential field at the initial time t_i. Mathematically, this is given by Equation (12), reproduced here:

$$\begin{equation*} \psi _\sigma ^{({\rm cr})f} = -\sum _{D}\sum _{j=t_i}^{t_f-1}\Delta ^j\mathbb {R}_D. \end{equation*}$$

This flux is generated by the current along the separator, Equation (16), which may be found by the inversion (18), and plugged into the energy Equation (17) to find the energy in excess of the potential field energy at time t_f in the MCC model with the FCE constraint applied at time t_i. This energy is a lower bound for any coronal energy model.