STRUCTURE, STABILITY, AND EVOLUTION OF MAGNETIC FLUX ROPES FROM THE PERSPECTIVE OF MAGNETIC TWIST

For this study we used vector magnetograms obtained by the Helioseismic and Magnetic Imager (HMI; Hoeksema et al. 2014) onboard the Solar Dynamics Observatory (SDO; Pesnell et al. 2012). The data for HMI Active Region Patches are disambiguated and deprojected to the heliographic coordinates with a Lambert (cylindrical equal area) projection method, resulting in a pixel scale of 003 or 0.36 Mm (Bobra et al. 2014). Liu et al. (2014b) have demonstrated that different projection methods make little impact on the deprojected map of a normal AR and on the calculation of the helicity flux.

The flares produced by AR 11817 were observed with the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) onboard SDO. The six EUV passbands of AIA, i.e., 131 Å (peak response temperature $\mathrm{log}T=7.0$), 94 Å ($\mathrm{log}T=6.8$), 335 Å ($\mathrm{log}T=6.4$), 211 Å ($\mathrm{log}T=6.3$), 193 Å ($\mathrm{log}T=6.1$) and 171 Å ($\mathrm{log}T=5.8$), can be used to calculate the differential emission measure (DEM) in the logarithmic temperature range $[5.5,7.5]$. We utilized the regularized inversion code developed by Hannah & Kontar (2012) to recover DEMs, using the most recent available AIA temperature response functions. The AIA level-1 data are further processed by applying the routines AIA_DECONVOLVE_RICHARDSONLUCY and AIA_PREP from the SolarSoft package before being fed into the DEM code.

One of the flares produced by AR 11817 is captured by the 1.6 m New Solar Telescope (NST) at Big Bear Solar Observatory (BBSO) with unprecedented high spatiotemporal resolution. The telescope adopts a state-of-art high-order adaptive optics system of 308 elements to correct atmospheric disturbances. Diffraction-limited imaging is achieved with the aid of speckle image reconstruction. The spatial resolution of the Hα images used in this study is 0085 (60 km) at a cadence of 15 s (see Wang et al. 2015, for more detail).

To monitor the evolution and complexity of the AR, we calculated the relative helicity flux across its photospheric boundary S with the following formula (Berger 1984; Liu & Schuck 2012; Liu et al. 2014b):

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dH}}{{dt}}\right|}_{S}=2{\displaystyle \int }_{S}({{\boldsymbol{A}}}_{p}\cdot {{\boldsymbol{B}}}_{t}){V}_{\perp n}\;{dS}-2{\displaystyle \int }_{S}({{\boldsymbol{A}}}_{p}\cdot {{\boldsymbol{V}}}_{\perp t}){B}_{n}\;{dS},\end{eqnarray} \tag{ 1 }$$

where ${{\boldsymbol{A}}}_{p}$ is the vector potential of the reference potential field ${{\boldsymbol{B}}}_{p}$ that has the same vertical component on the photospheric boundary; t and n refer to the tangential and normal directions, respectively. ${{\boldsymbol{V}}}_{\perp }$ is the photospheric velocity that is perpendicular to magnetic field lines. To obtain ${{\boldsymbol{V}}}_{\perp }$, we applied the Differential Affine Velocity Estimator for Vector Magnetograms (DAVE4VM; Schuck 2008) to the time-series of deprojected, registered vector magnetograms. We then subtracted the field-aligned plasma flow $({\boldsymbol{V}}\cdot {\boldsymbol{B}}){\boldsymbol{B}}/{B}^{2}$ from the velocity derived by DAVE4VM to yield ${{\boldsymbol{V}}}_{\perp }$ (Liu & Schuck 2012). The window size used in DAVE4VM is chosen to be 19 pixels, following Liu et al. (2014b). The first term in the above equation represents helicity injection due to the emergence of twisted flux tubes into the corona; the second is due to photospheric motions that shear and braid field lines. Similarly, both flux emergence (${V}_{\perp n}$) and tangential motions (${{\boldsymbol{V}}}_{\perp t}$) contribute to the Poynting flux across the photospheric boundary S, which is given by (Kusano et al. 2002)

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dP}}{{dt}}\right|}_{S}=\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{S}{B}_{t}^{2}{V}_{\perp n}\;{dS}-\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{S}({{\boldsymbol{B}}}_{t}\cdot {{\boldsymbol{V}}}_{\perp t}){B}_{n}\;{dS}.\end{eqnarray} \tag{ 2 }$$

Readers are referred to Appendix B for an error analysis of helicity and Poynting flux.

To understand the magnetic connectivities within the active region and their evolution, we used the code package developed by T. Wiegelmann, which utilizes the "weighted optimization" method (Wiegelmann 2004; Wiegelmann et al. 2012) to build an NLFFF model that approximates the coronal field (see Appendix A for the quality of NLFFF). To best suit the force-free condition, the vector magnetograms are "pre-processed" (Wiegelmann et al. 2006b) before being taken as the photospheric boundary. Our calculation is performed within a box of $512\times 256\times 256$ uniformly spaced grid points, whose photospheric FOV is shown in Figure 1. Further, the free magnetic energy accumulated in the AR can be derived by subtracting the magnetic energy of the corresponding potential field from that of the NLFFF (see Appendix B for an error analysis of magnetic free energy). The potential field is calculated by the Green function method. We calculated the free magnetic energy regularly on an hourly cadence, but at the highest cadence available (12 minutes) around flares. It is clear that the NLFFF is not a good model of the coronal AR field during the impulsive phase of the flares, as plasma is accelerated, i.e., forces are significant, primarily in this phase (Zhang et al. 2001). However, the flare-related changes of the coronal field can be inferred from a comparison of the NLFFF before and after the impulsive flare phase.

Zoom In Zoom Out Reset image size

After obtaining the NLFFF, we refined the photospheric computational grid by 16 times and traced magnetic field lines pointwise with a fourth-order Runge–Kutta method to ensure high precision. The mapped footpoints of field lines were used to calculate the squashing factor Q of elemental magnetic flux tubes (Titov et al. 2002; Titov 2007). Basically, for a mapping defined by the two footpoints of a field line ${{\rm{\Pi }}}_{12}:{{\boldsymbol{r}}}_{1}({x}_{1},{y}_{1})\mapsto {{\boldsymbol{r}}}_{2}({x}_{2},{y}_{2})$, the Jacobian matrix associated with the mapping is

$$\begin{eqnarray}{D}_{12}=\left[\displaystyle \frac{\partial {{\boldsymbol{r}}}_{2}}{\partial {{\boldsymbol{r}}}_{1}}\right]=\left(\begin{array}{ll}\partial {x}_{2}/\partial {x}_{1} & \partial {x}_{2}/\partial {y}_{1}\\ \partial {y}_{2}/\partial {x}_{1} & \partial {y}_{2}/\partial {y}_{1}\end{array}\right)\equiv \left(\begin{array}{ll}a & b\\ c & d\end{array}\right),\end{eqnarray} \tag{ 3 }$$

and then the squashing factor associated with the field line is defined as (Titov et al. 2002)

$$\begin{eqnarray}&&Q\equiv \displaystyle \frac{{a}^{2}+{b}^{2}+{c}^{2}+{d}^{2}}{| {B}_{n,1}({x}_{1},{y}_{1})/{B}_{n,2}({x}_{2},{y}_{2})| },\end{eqnarray} \tag{ 4 }$$

where ${B}_{n,1}({x}_{1},{y}_{1})$ and ${B}_{n,2}({x}_{2},{y}_{2})$ are the components normal to the plane of the footpoints, in our case, the photosphere, and their ratio is equivalent to the determinant of D₁₂. Quasi-separatrix layers (QSLs) as defined by high-Q values are often complex three-dimensional structures, hence their visualization can be facilitated by Q-maps cutting across the QSL of interest. Pariat & Démoulin (2012) discussed three relevant methods and found that the method (their Method 3) utilizing the field-line mappings between the cutting plane and the footpoint planes gives optimal results, i.e., the chain rule of the Jacobian is used to calculate D₁₂, from a set of field lines threading the cutting plane in a neighborhood of ${{\boldsymbol{r}}}_{c}({x}_{c},{y}_{c})$:

$$\begin{eqnarray}&&{D}_{12}=\left[\displaystyle \frac{\partial {{\boldsymbol{r}}}_{2}}{\partial {{\boldsymbol{r}}}_{1}}\right]=\left[\displaystyle \frac{\partial {{\boldsymbol{r}}}_{2}}{\partial {{\boldsymbol{r}}}_{c}}\right]\times \left[\displaystyle \frac{\partial {{\boldsymbol{r}}}_{c}}{\partial {{\boldsymbol{r}}}_{1}}\right]\end{eqnarray} \tag{ 5 }$$

where $({x}_{1},{y}_{1})$ and (x₂, y₂) are again the corresponding photospheric footpoints and the Jacobian associated with the mapping ${{\rm{\Pi }}}_{1c}:{{\boldsymbol{r}}}_{1}({x}_{1},{y}_{1})\mapsto {{\boldsymbol{r}}}_{c}({x}_{c},{y}_{c})$ is given by its inverse, i.e.,

$$\begin{eqnarray}\left[\displaystyle \frac{\partial {{\boldsymbol{r}}}_{c}}{\partial {{\boldsymbol{r}}}_{1}}\right]=\displaystyle \frac{1}{| {B}_{n,c}({x}_{c},{y}_{c})/{B}_{n,1}({x}_{1},{y}_{1})| }\left(\begin{array}{ll}\partial {y}_{1}/\partial {y}_{c} & -\partial {x}_{1}/\partial {y}_{c}\\ -\partial {y}_{1}/\partial {x}_{c} & \partial {x}_{1}/\partial {x}_{c}\end{array}\right).\end{eqnarray} \tag{ 6 }$$

We found, however, that field lines touching the cutting plane, i.e., ${B}_{n,c}({x}_{c},{y}_{c})\to 0$, introduce spurious high-Q structures. These can be effectively eliminated via two different approaches: for a field line making a small angle to the cutting plane, one may calculate its Q-value from a set of field lines originating from the neighborhood of its photospheric footpoint, which is equivalent to Method 2 in Pariat & Démoulin (2012). One problem is that the new Q-values in replacement of the original spurious high-Q structures may be much smaller than those in the neighborhood, calculated with Equation (5), therefore giving rise to artificial abrupt changes. The latter are not intrinsic to the magnetic structure under study but rather are due to numerical errors that are generally different for the two methods. Alternatively, one may switch locally to a new plane that is perpendicular¹⁰ to this particular field line to calculate its Q-value again with Equation (5). This new value is often consistent with those in the surroundings.

We further calculated the twist number ${{ \mathcal T }}_{w}$ to measure how many turns two infinitesimally close field lines wind about each other (Berger & Prior 2006, Equation (16)):

$$\begin{eqnarray}&&{{ \mathcal T }}_{w}={\displaystyle \int }_{L}\displaystyle \frac{{\mu }_{0}{J}_{\parallel }}{4\pi B}\;{dl}={\displaystyle \int }_{L}\displaystyle \frac{{\rm{\nabla }}\times {\boldsymbol{B}}\cdot {\boldsymbol{B}}}{4\pi {B}^{2}}\;{dl}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{L}\alpha \;{dl},\quad {\rm{if}}\;{\rm{}}{\rm{\nabla }}\times {\boldsymbol{B}}=\alpha {\boldsymbol{B}},\end{eqnarray} \tag{ 8 }$$

where α is the force-free parameter, and ${\rm{\nabla }}\times {\boldsymbol{B}}\cdot {\boldsymbol{B}}/4\pi {B}^{2}$ can be regarded as a local density of twist along the field line. In Appendix C we address the relation among ${{ \mathcal T }}_{w}$, twist number of an individual magnetic field line, ${{ \mathcal T }}_{g}$, twist number of any curve about an axis in a general geometry (see Equation (12) in Berger & Prior 2006), and N, twist number of a field line about the axis of a cylindrical flux tube. Simply put, ${{ \mathcal T }}_{g}$ is the generalization of N, and ${{ \mathcal T }}_{w}$ approaches ${{ \mathcal T }}_{g}$ and N in the vicinity of the axis of a nearly cylindrically symmetric flux tube, but deviates otherwise (see Appendix C). While ${{ \mathcal T }}_{g}$ is the quantity usually considered in stability analyses, its computation depends upon the correct determination of the axis, which is non-trivial for numerically given fields and demanding for a large number of NLFFFs. Therefore, we employ ${{ \mathcal T }}_{w}$ in the analysis of this paper.

For a perfect NLFFF, α is a constant for each individual field line. In practice, we linearly interpolate ${\boldsymbol{B}}$ and ${\rm{\nabla }}\times {\boldsymbol{B}}$ on the line position and carry out the integration with a 5-point Newton–Cotes formula. We refrained from linearly interpolating α on the non-grid points because of its inherent nonlinearity. Consequently, a map of twist numbers ("twist map" hereafter) is yielded at the photosphere or in a cutting plane by assigning the twist number of each field line to the position where the line threads the plane. Alternatively, Inoue et al. (2011) obtained the photospheric twist map by multiplying the photospheric maps of α and field line length. To address the issue that α is often not the same at the conjugate footpoints, they took the mean α at the conjugate footpoints for each individual field line (see also Chintzoglou et al. 2015). A comparison between the different approaches to calculate a twist number is given in Appendix C.

The codes that are developed by R. Liu and J. Chen to calculate Q and ${{ \mathcal T }}_{w}$ are available online at http://staff.ustc.edu.cn/~rliu/qfactor.html.

AR 11817 emerged on August 10 as a simple β configuration (P1–N1; Figure 1(a)) with an almost east–west orientation obeying Hale's law. As time progressed, P1–N1 departed from each other and grew in size. Meanwhile, a new bipole P2–N2 emerged on the neutral line of P1–N1, defying Hale's law (Figures 1(c)–(e)). There were significant shearing and converging motions, with P2 moving westward and N2 moving northeastward (Figure 1(f)). As a result the pair exhibited a slow clockwise rotation, which could be explained by the emergence of a magnetic tube with left-handed writhe. The tube might be downward-kinked because of the converging motion, which mainly manifests as N2's moving toward P2. In contrast, no clear sign of writhe is seen for P1–N1. With the emergence of P2–N2, AR 11817 developed into a complex field configuration and produced a series of flares in the two days during August 11–12 (Table 1). After the M1.5 flare on August 12, however, AR 11817 became dormant and it would wait until August 14 to produce C-class flares again, none of which, however, exceeded C5.

Table 1 lists all the C-and-above flares taking place during the three days 2013 August 10–12, during which the longitudinal center of AR 11817 changed from E48° to E10°. The flares studied in this paper are numbered and boldfaced. Note that the C6.7 and C8.4 flares on August 11 are very close in time and practically can be regarded as a single flare with two peaks. However, the C6.7 flare is apparently confined, resulting in no CME, while for the subsequent C8.4 flare an ejection of a bubble-like structure leading up to a CME is seen in AIA images (not shown), which is hence conventionally classified as an eruptive flare. We ignored two C1 flares in this study: one occurred before any significant emergence of P2–N2; the other occurred during the decay phase of Flare No. 5 (see also the top panel of Figure 2); both last for a duration of only several minutes, less than the cadence of the HMI vector magnetograms (12 minutes).

Zoom In Zoom Out Reset image size

Figure 2 shows the helicity and energy injected into AR 11817 during August 10–12. The panels from top to bottom show the GOES 1–8 Å light curve, the temporal profiles of magnetic fluxes, helicity and energy fluxes, and the free magnetic energy derived from the NLFFF. Significant flux emergence begins from the end of August 10. Correspondingly, the helicity injection rate is quickly enhanced, but the shear term is dominant over the emergence term. In contrast, the Poynting flux from the shear term is often slightly lower than that from the emergence term. This appears to be common among ARs (Liu et al. 2014b). The accumulated helicity reaches −2 × 10⁴² Mx² when the first eruptive flare (No. 5) took place on August 11, and −4 × 10⁴² Mx² when the second eruptive flare (No. 6) occurred on August 12. Both the helicity and Poynting fluxes decreased from about 20 UT on August 12 onward, which might explain that AR 11817 became flare-quiet after Flare No. 6.

By the definition of MFRs (Section 1), the field lines winding around the rope axis must have similar magnetic connectivities, with the magnitude of their twist number exceeding unity. Hence, in the photospheric twist map (Figure 3(b)), taking Flare No. 3 for example, an MFR manifests itself straightforwardly as two conjugate compact regions with enhanced twist numbers of the same sign, which host the footpoints of the rope. The cross section of the MFR can be displayed by the twist number of field lines threading the cutting plane (Figures 3(e) and (h)). Comparing the photospheric (Figure 3(b)) and vertical (Figures 3(e) and (h)) twist maps with the corresponding Q-maps (Figures 3(a), (d), and (g)), one can see that the flux rope is bounded by high-Q lines, i.e., a QSL that separates the twisted field from untwisted field (Titov et al. 2002). Several short bald-patch sections of the PIL are typically identified under the MFR in the magnetograms, indicating a complex topology that includes bald patch separatrix surfaces (BPSSs; Titov & Démoulin 1999). In a bald patch the field lines graze the photosphere; such a representative field line is plotted in orange in Figure 3. Also shown are a representative twisted field line (green) and a QSL line (magenta).

Figure 4 shows perspective views of the MFR. The key to identifying and precisely locating the axis of the MFR is the vertical twist map. We have computed twist maps in many vertical cut planes, all oriented in the y direction, in a wide range of x values and traced a field line from the peak-$| {{ \mathcal T }}_{w}|$ point in each map. It was found that all such field lines traced in the range x = [74.5, 107.3] Mm coincide within the limits of numerical accuracy. This includes nearly the whole axis of the MFR, except for a very short section (${\rm{\Delta }}x\lt 2.5$ Mm) at its east end. Next we have placed a vertical cut plane at a point where the MFR axis runs horizontally (e.g., at the apex point) such that it intersects the cut plane perpendicularly and checked whether the in-plane field vectors display a rotational pattern centered at the intersection point, as expected for an MFR. This is indeed the case, and also the current density normal to the cut plane is enhanced in this area and thus consistent with the existence of an MFR (Figure 4(e)). Similarly to the Q map in Figure 3(d), the MFR displays a rather compact and vertically elongated cross section. Choosing field line start points in the cut plane such that they follow this shape, the displays of Figures 4(a)–(d) are obtained. These confirm that a largely coherent MFR has formed in the range of strongest flux cancellation (75 ≲ x ≲ 95 Mm). The MFR extends further westward but is less twisted and less coherent in this area, where flux from various photospheric flux patches joins the forming rope.

Zoom In Zoom Out Reset image size

This analysis was also performed for the NLFFF before the Flares No. 4 (August 11, 21:22 UT) and No. 6 (August 12, 10:22 UT), confirming that the peak-$| {{ \mathcal T }}_{w}|$ point in vertical twist maps locates the axis of the MFR in the considered AR. Some care must be executed, however, as ${{ \mathcal T }}_{w}$ may peak away from the MFR axis in other configurations, where the current has a broad spatial distribution (uniform across the rope or even a hollow profile). In such cases, ${{ \mathcal T }}_{w}$ has a minimum at the MFR axis (which also facilitates locating the axis; see Appendix C). Moreover, in areas where the MFR is not fully coherent it may not be possible to define an axis and ${{ \mathcal T }}_{w}$ may simply peak at a long field line that runs through a volume of high shear.

The well-defined boundary of the main flux rope gives us an opportunity to test its stability with respect to the helical kink mode against the classical theory, namely, the Kruskal–Shafranov limit of one field line turn about the MFR axis, modified for the presence of photospheric line tying, which raises it to $N=1.25$ or ${\rm{\Phi }}=2.5\pi$ for the twist angle ${\rm{\Phi }}(r)={B}_{\phi }{L}_{z}/({{rB}}_{z})$. Using ${B}_{\phi }={\mu }_{0}{I}_{z}/(2\pi r)$, a limiting axial current I_ζ for an idealized cylindrical flux rope with radius r₀ and length L_z is obtained (Schindler 2006, p. 522),

$$\begin{eqnarray}&&{I}_{\zeta }=1.25\displaystyle \frac{{(2\pi {r}_{0})}^{2}{B}_{z}({r}_{0})}{{\mu }_{0}{L}_{z}}.\end{eqnarray} \tag{ 9 }$$

Based on this, we introduced a dimensionless parameter,

$$\begin{eqnarray}&&{\rm{KS}}\equiv \displaystyle \frac{{F}_{p}}{{F}_{t}}=\displaystyle \frac{{\displaystyle \int }_{0}^{{r}_{0}}{\displaystyle \int }_{0}^{{L}_{z}}{B}_{\phi }(r)\;{dz}\;{dr}}{{\displaystyle \int }_{0}^{{r}_{0}}{\displaystyle \int }_{0}^{2\pi }{B}_{z}(r)r\;d\phi \;{dr}},\end{eqnarray} \tag{ 10 }$$

where F_p is the poloidal flux and F_t the toroidal flux. One can immediately see that KS is related to N with two assumptions (simplicities): (1) B_ϕ is distributed linearly over the radius from the center of the rope, or equivalently, I_z has a uniform distribution in the cross-section of the rope, and (2) B_z is uniformly distributed across the cross section, i.e.,

$$\begin{eqnarray*}{\rm{KS}} & \approx & \displaystyle \frac{{B}_{\phi 0}{r}_{0}{L}_{z}/2}{{B}_{z0}\pi {r}_{0}^{2}}=N({r}_{0})\\ & = & \displaystyle \frac{{\mu }_{0}{I}_{z}{L}_{z}}{{(2\pi {r}_{0})}^{2}{B}_{z0}},\end{eqnarray*}$$

where ${B}_{\phi 0}$ and ${B}_{z0}$ are given at $r={r}_{0}$. With Equation (9), one can see that the threshold of $| {\rm{KS}}| =1.25$. In our case, we integrated j_x over the rope's cross section in the cutting plane ${{\rm{C}}}_{{\rm{A}}}$ (Figure 3) to estimate F_p, i.e., ${F}_{p}={\mu }_{0}L/4\pi \int {j}_{x}\;{dS}$, where L ≈ 37.83 Mm is given by the length of the rope axis in Figure 4. Similarly, F_t is estimated by integrating B_x over the same cross section. The result gives ${\rm{KS}}=-1.46$, exceeding the KS threshold. In comparison, for a group of selected field lines near the rope axis (rainbow colored in Figure 4), we found that $\langle {T}_{w}\rangle =-1.88\pm 0.06$ with $\mathrm{max}({T}_{w})=-1.96$ and $\langle {T}_{g}\rangle =-1.28\pm 0.20$ with $\mathrm{max}({T}_{g})=-1.49$. Hence, both the integral property of the rope as characterized by KS and the peak "local" twist number (${{ \mathcal T }}_{w}$ or ${{ \mathcal T }}_{g}$) of individual field lines (near the flux rope axis) suggest that the MFR may be marginally kink-unstable. However, the exact threshold of instability depends on further parameters, in particular on the strength of the external toroidal (shear) field component (Kliem et al. 2014), and it is also higher for very slim flux ropes (Baty 2001) and if the rope is not fully coherent, as in the present case. Moreover, MFRs do not always have a well-defined high-Q boundary (see, e.g., Liu et al. 2014a); for the current MFR, the eastern footpoint region is fully enclosed by high-Q lines, but the western footpoint region is only partially enclosed (see Figure 3(b)). Hence, the above analysis cannot be easily applied to all occasions.

Therefore, we tested the stability of the MFR by using the NLFFF as the initial condition in an MHD code, employing the same grid. We used a simple density model which yields a slow decrease of the Alfvén velocity with distance from the flux concentrations (see, e.g., Török & Kliem 2005). No perturbation was applied. Only relatively small, nearly vertical motions of the MFR ensued spontaneously, which subsequently decreased monotonically. The resulting stable numerical MHD equilibrium shows an MFR close to the one in the corresponding initial NLFFF. Similar results are obtained for the MFR before Flares No.4 (21:22 UT on August 11) and 6 (10:22 UT on August 12). We conjecture that the stability of the MFR in the NLFFF only shortly before or just at the onset of the flares results from the minimum-current assumption in the disambiguation of the vector magnetograms.

Further, we calculated in the cutting plane C_A the decay index $n=-d\mathrm{ln}{B}_{t}/d\mathrm{ln}h$, where h is the height above the photosphere and B_t is the transverse component of the potential field. Numerical experiments have demonstrated that an MFR becomes torus unstable when the decay index of the external poloidal field reaches a critical value of $\sim 1.5$ (e.g., Török & Kliem 2007), and B_t is a good proxy of the external poloidal field at the axis of the MFR (e.g., Kliem et al. 2013). The twist and Q maps in C_A clearly show that the MFR is located below about 4 Mm above the photosphere (Figures 3(e) and (h)). In comparison with the map of the decay index (Figure 3(f)), one sees that the MFR resides in a region stable to the torus instability unless it experiences a large disturbance that pushes it up to much higher altitudes. In fact, for all the flares studied, the MFR is located below the n = 1 contour (see the insets in Figures 6–10). In cases No. 4 (confined) and 6 (eruptive), the MFR is relatively high-lying, reaching the n = 1 contour before the flare (Figures 9 and 10), but still well below the threshold of the torus instability ($n\sim 1.5$).

Hence, throughout the time series analyzed, it appears that the MFR is close to the threshold of the helical kink instability but far away from that of the torus instability. We propose a scenario in which the helical kink mode commences at the onset of all flares studied and the consequent lifting of the MFR does not reach the threshold height for the onset of the torus instability in the confined Flares No. 1–3, but, starting from a greater height, reaches it in the eruptive Flare No. 6 as well as the combined eruption seen as Flares No. 4 and 5.

The vertical twist map in Figure 3(e) reveals an unexpected layer of reversed twist (and α) on the north edge of the MFR and an area of slightly enhanced (negative) twist values under the MFR and the layer. The nearby vertical cut in Figure 4(e) shows enhanced current density and a rotational pattern of the in-plane field components in this area, i.e., an indication for the existence of a second, smaller MFR under the main MFR. Such a double-decker flux rope configuration was suggested previously as a candidate for partial eruptions (Liu et al. 2012; Cheng et al. 2014; Zhu & Alexander 2014). We have investigated the structure using our series of vertical twist maps, also rotating the cut planes about the vertical, to better display the axial current and azimuthal field structure (Figure 5(b)). A shallow peak of the twist number, ${{ \mathcal T }}_{w}\approx -1.1$, is found, which coincides with a shallow maximum of the current density and the center of a rotational structure of the azimuthal field component (Figure 5(b)). The rotational structure exists all around the center point, clearly indicating that the region is a flux rope, not simply a sheared volume, and that the peak-$| {{ \mathcal T }}_{w}|$ point lies on the MFR axis. This is substantiated by the field line plot in Figure 5(a) which displays a double-decker configuration of two left-handed flux ropes, and by the isosurface of ${{ \mathcal T }}_{w}=-1$ (Figure 5(c)).

Zoom In Zoom Out Reset image size

A shear layer is a natural consequence if both flux ropes in a double-decker configuration have the same handedness. The field lines in the bottom layer of the upper rope then have a direction different from the field lines in the top layer of the lower rope (see Figure 12 in Liu et al. 2012). If the system is driven, such that the ropes evolve in flux content and position, currents are likely induced in the shear layer, i.e., $\alpha \ne 0$ and hence ${{ \mathcal T }}_{w}\ne 0$. If the direction of the axial field in the configuration is uniform, like in the present case, then the axial current, α, and ${{ \mathcal T }}_{w}$ in the shear layer are opposite in sign to the two MFR, in agreement with Figure 3(e). Since oppositely directed currents repel each other, the shear layer should have a stabilizing influence on the configuration.

We calculated a series of photospheric twist maps based on the NLFFF modeling, focusing on a small region (black rectangle in Figure 1(f)) that encloses the major PIL, where the magnetic field is highly sheared and most shearing and converging flows are detected. Signatures of an MFR along the major PIL start to appear in the twist map from about 12:46 UT on August 11, about 1.5 hr before Flare No. 1, and become quite prominent before the occurrence of the C2.1 flare (No. 3) at 19:27 UT on August 11. One can see from Figure 3 that both the footpoint regions and the cross sections of the MFR are featured by dark blue colors, in agreement with the negative helicity dominating in this AR, and enclosed by high-Q lines (Figure 3), suggesting the presence of a QSL at the rope boundary. If a coherent MFR appears in the twist map calculated in a cutting plane by visual inspection, we recorded the peak $| {{ \mathcal T }}_{w}|$ within the MFR's cross section. The result is plotted in the bottom panel of Figure 2, which shows a common trend that the peak $| {{ \mathcal T }}_{w}|$ temporarily increases before the flare while it decreases after the flare peak, both for confined and eruptive flares. This is also seen in the environment of the peak-$| {{ \mathcal T }}_{w}|$ point. We will now examine this trend case by case.

Figures 6–10 show the twist maps across each flare studied (Table 1). The decrease in magnetic twist is often related to a decrease in the free magnetic energy of the overall AR (except Flares No. 2 and 5), but due to the relatively small volume of the MFR, the transient increase in magnetic twist does not translate to a similar increase in the free magnetic energy, although the latter possesses a gradually increasing profile with time (Figure 2). It is remarkable that, despite that Flares No. 4 and 5 are lumped together, the twist map still indicates a temporary increase in magnetic twist just prior to Flare No. 5. In comparison, although it clearly decreases across the combined eruption seen as Flares No. 4 and 5 (Figure 9), the free magnetic energy does not change appreciably across No. 5. In Flare No. 2 (Figure 7), the weakest of our sample, the flare-related energy release is too weak to stand out against the general trend of energy build-up in the AR, but the change of magnetic twist across this flare is visible. For Flare No. 1 (Figure 6), the free magnetic energy decreases already 10–20 minutes before the flare (see also Jing et al. 2009). A similar issue happens to the magnetic twist when its change is detected 10–20 minutes before Flare No. 3 (Figure 8). This might indicate a limited accuracy of the underlying magnetograms or of the NLFFF model, but it could also indicate that the flare is triggered by a photospheric evolution that starts already before the flare. We also noted that the height of the MFR gradually increases with time during the sequence studied (compare the insets in Figures 6–10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We assume that the rise of the peak twist number $| {{ \mathcal T }}_{w}|$ at the axis of the MFR just prior to the flares is correlated with the on-axis value of the classical twist number $| {{ \mathcal T }}_{g}|$. The relationship between these twist numbers is analyzed in Appendix C, where we find that, at the axis, ${{ \mathcal T }}_{w}$ approaches the sum of ${{ \mathcal T }}_{g}$ and a geometry-dependent term which measures the deviation from cylindrical symmetry (see Equation (21)). The expectation that such a deviation is of minor influence on the stability of the MFR compared to the twist underlies our assumption, which should be substantiated in future investigations.

The peak magnetic twist at the axis of the MFR can increase due to a twisting of the rope by rotational motions at its footpoints, or due to a lengthening of the rope by reconnection. The latter involves a change in the footpoints and likely also leads to a greater height of the rope. In the derived velocity maps by DAVE4VM, there are no obvious counterclockwise rotations associated with the footpoint regions as expected for the buildup of the left-handed twist, but we did find positive correlations among the peak $| {{ \mathcal T }}_{w}|$, its height, and the length of the field line threading the peak $| {{ \mathcal T }}_{w}|$ position (bottom panel of Figure 2). Thus, primarily reconnection should have caused the transient increase in $| {{ \mathcal T }}_{w}|$ prior to the flares.

For eruptive flares, it is readily understood that the decrease in magnetic twist must result from the expulsion of the flux rope, although a flux rope still survives the eruption in both cases studied (Flares No. 5 and 6), suggesting that most likely a partial expulsion of the flux rope occurs during the eruption. In Flare No. 5, the remnant MFR is significantly reduced in flux, height, and twist (Figure 9), eventually disappears, and then reforms about 2 hr before Flare No. 6 (see the bottom panel of Figure 2). In that flare, the MFR splits into two remaining branches at a lower height immediately after the flare peak (Figure 10), suggesting that a reconfiguration of the flux rope might be ongoing. It is interesting that the two branches take the form of two interlocking Js, which is the reverse of the double-J-to-S transformation that is typical of tether-cutting reconnection (e.g., Liu et al. 2010). For confined flares, the decrease in magnetic twist must indicate a conversion of magnetic energy primarily into plasma heating, but what mechanism causes this decrease is unclear until we have further analyzed Flare No. 3 using the high-resolution NST data (Section 3.5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we highlight the confined C2.1 flare taking place at 19:27 UT on August 11 (Flare No. 3). With the aid of NST observations (see the detailed analysis by Wang et al. 2015), this event sheds some light on the flaring mechanism. Figure 11 shows the flaring process observed in AIA 131 Å, which manifests as the brightening of a small sigmoid extending less than $50^{\prime\prime}$ in the east–west direction. In the high-resolution Hα images, however, one sees a filament with a similar reverse S shape. Wang et al. (2015) have demonstrated that the filament traces the MFR. The filament axis is apparently writhed after the flare onset (Figures 11(f)–(g)), with the counterclockwise rotation being consistent with the conversion of left-handed twist (Green et al. 2007). The deformation of the filament axis appears to be moderate and the filament soon settles back down in approximately six minutes to a similar configuration as before the flare, indicating that the flare is confined. The fine threads of the filament do not appear to be twisted by significantly more than one turn, which is consistent with the range of twist values of the MFR derived from the NLFFF modeling (see Appendix C). Superimposing on the AIA 131 Å image a representative field line of the MFR in the pre-eruption NLFFF (red dotted line in Figure 11), one can see that the MFR matches the sigmoid very well. This is consistent with the sigmoid being energized beneath the erupting MFR (e.g., Titov & Démoulin 1999). In the process of eruption the coherence of the filament is distorted, with some brightenings occurring in its rising part (Figures 11(f)–(g)). This suggests that reconnection also occurs with the overlying field, which halts the eruption according to the analysis of the decay index in Section 3.2. Thus, although the writhe must be back-converted into twist as the filament settles back near its original shape, the MFR loses part of its twist by reconnection with the overlying field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A diagnosis of the flaring plasma with the DEM method can shed more light on the energy release mechanism of this flare. One can see that immediately after the flare onset, the sigmoid is visible in each of the six EUV passbands of AIA (Figure 12). We obtain the average DEMs for 3 by 3 pixels centered on three selected loci along the sigmoid (red, olive and green crosses) as well as on neighboring coronal loops (blue cross) to serve as a reference (Figure 12(e)). The DEMs of the sigmoid exhibit a significant peak at about 10 MK, besides the peak at about 1 MK that can mainly be attributed to the coronal background. Both the cool and hot DEM peaks of the sigmoid are significantly elevated above those of the neighboring coronal loops, suggesting that the sigmoid is both dense and hot. This sigmoid is obviously associated with the erupting filament observed in Hα; however, its size and shape do not follow the expanding and writhing filament, but rather remain close to the initial state (Figure 11). This corresponds to the standard picture of the sigmoid being formed by current dissipation in a separatrix or quasi-separatrix layer underneath an activated flux rope. Figure 3(d) indeed shows the highest squashing factor where the QSL intersects itself under the MFR; this intersection is suggestive of the presence, or beginning formation, of a hyperbolic flux tube, where the flare reconnection is expected to develop (e.g., Savcheva et al. 2012). A location of the sigmoid under the flux rope is also supported by its high density.

Let ${\boldsymbol{x}}(s)$ be a smooth, non-self-intersecting curve parametrized by arclength s and ${\boldsymbol{y}}(s)$ a second such curve surrounding ${\boldsymbol{x}}$, jointly defining a ribbon. Points on ${\boldsymbol{y}}$ are expressed as ${\boldsymbol{y}}(s)={\boldsymbol{x}}(s)+\varepsilon \rho (s)\hat{{\boldsymbol{V}}}(s)$, where $\hat{{\boldsymbol{V}}}(s)$ is a unit vector normal to ${\boldsymbol{x}}$, so that ${\boldsymbol{y}}$ is also parametrized by the arclength on ${\boldsymbol{x}}$, which we will refer to as the axis curve. With φ denoting the rotation angle made by ${\boldsymbol{y}}$ about ${\boldsymbol{x}}$ and $\hat{{\boldsymbol{T}}}(s)$ denoting the unit tangent vector to ${\boldsymbol{x}}(s)$, the twist number of curve ${\boldsymbol{y}}$ about the axis ${\boldsymbol{x}}$ is given by (Equation (12) in Berger & Prior 2006)

$$\begin{eqnarray}&&{{ \mathcal T }}_{g}=\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{{\boldsymbol{x}}}\displaystyle \frac{d\varphi }{{ds}}{ds}=\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{{\boldsymbol{x}}}\hat{{\boldsymbol{T}}}(s)\cdot \hat{{\boldsymbol{V}}}(s)\times \displaystyle \frac{d\hat{{\boldsymbol{V}}}(s)}{{ds}}\;{ds}.\end{eqnarray} \tag{ 13 }$$

Applying this equation to magnetic field lines, we have $\hat{{\boldsymbol{T}}}=\hat{{\boldsymbol{B}}}\equiv {\boldsymbol{B}}/B$. In order to derive the relationship between the general definition of twist number, ${{ \mathcal T }}_{g}$, in a magnetic field and the alternative twist number ${{ \mathcal T }}_{w}$ (Equation (7)), suggested in the literature to approximate ${{ \mathcal T }}_{g}$ in the vicinity of ${\boldsymbol{x}}$ ($\varepsilon \ll 1$), we express the rate of direction change $d\hat{{\boldsymbol{V}}}/{ds}$ through the structure (the derivatives) of the magnetic field. Denote the infinitesimally small distance between ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$ at point s by $\delta {\boldsymbol{r}}(s)=\varepsilon \rho (s)\hat{{\boldsymbol{V}}}(s)={\boldsymbol{y}}(s)-{\boldsymbol{x}}(s)$. This quantity changes at a rate

$$\begin{eqnarray}&&\displaystyle \frac{d}{{ds}}\delta {\boldsymbol{r}}=\displaystyle \frac{d{\boldsymbol{y}}}{{ds}}-\displaystyle \frac{d{\boldsymbol{x}}}{{ds}}.\end{eqnarray} \tag{ 14 }$$

We have $d{\boldsymbol{x}}/{ds}=\hat{{\boldsymbol{T}}}(s)=\hat{{\boldsymbol{B}}}({\boldsymbol{x}}(s))$. Using arclength ${s}^{\prime }$ on ${\boldsymbol{y}}$ and denoting the unit tangent vector at ${\boldsymbol{y}}(s)$ by ${\hat{{\boldsymbol{T}}}}^{\prime }(s)$, we can write $d{\boldsymbol{y}}/{ds}=d{\boldsymbol{y}}/{{ds}}^{\prime }\cdot {{ds}}^{\prime }/{ds}$ $={\hat{{\boldsymbol{T}}}}^{\prime }(s)\cdot {{ds}}^{\prime }/{ds}$ $\;=\;\hat{{\boldsymbol{B}}}({\boldsymbol{y}}(s))\cdot {{ds}}^{\prime }/{ds}$. For sufficiently small ε the difference between the two arclengths can be neglected, ${{ds}}^{\prime }\approx {ds}$, so that

$$\begin{eqnarray}&&\displaystyle \frac{d}{{ds}}\delta {\boldsymbol{r}}\simeq \hat{{\boldsymbol{B}}}({\boldsymbol{x}}(s)+\delta {\boldsymbol{r}}(s))-\hat{{\boldsymbol{B}}}({\boldsymbol{x}}(s)).\end{eqnarray} \tag{ 15 }$$

Taylor expanding the first term yields

$$\begin{eqnarray}&&\displaystyle \frac{d}{{ds}}\delta {\boldsymbol{r}}\simeq \delta {\boldsymbol{r}}\cdot \displaystyle \frac{\partial \hat{{\boldsymbol{B}}}}{\partial {\boldsymbol{r}}},\end{eqnarray} \tag{ 16 }$$

so that

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{{\boldsymbol{V}}}}{{ds}}\simeq \hat{{\boldsymbol{V}}}\cdot \displaystyle \frac{\partial \hat{{\boldsymbol{B}}}}{\partial {\boldsymbol{r}}}-\displaystyle \frac{1}{\rho }\displaystyle \frac{d\rho }{{ds}}\hat{{\boldsymbol{V}}},\end{eqnarray} \tag{ 17 }$$

where all quantities are taken at position ${\boldsymbol{x}}(s)$. The local density of ${{ \mathcal T }}_{g}$ can now be written as

$$\begin{eqnarray}\displaystyle \frac{d{{ \mathcal T }}_{g}}{{ds}} & = & \displaystyle \frac{1}{2\pi }\hat{{\boldsymbol{T}}}\cdot \hat{{\boldsymbol{V}}}\times \displaystyle \frac{d\hat{{\boldsymbol{V}}}}{{ds}}\\ & & \simeq \displaystyle \frac{1}{2\pi }\hat{{\boldsymbol{T}}}\cdot \hat{{\boldsymbol{V}}}\times \left(\hat{{\boldsymbol{V}}}\cdot \displaystyle \frac{\partial \hat{{\boldsymbol{B}}}}{\partial {\boldsymbol{r}}}\right),\end{eqnarray} \tag{ 18 }$$

In evaluating $\hat{{\boldsymbol{V}}}\cdot \frac{\partial \hat{{\boldsymbol{B}}}}{\partial {\boldsymbol{r}}}$ we split $\partial {\boldsymbol{B}}/\partial {\boldsymbol{r}}$ into symmetric and antisymmetric parts (a similar technique is used, e.g., in the proof of the Cauchy–Helmholtz theorem for the form of the local velocity field in a fluid element; see, e.g., Kiselev et al. 1999, p. 37):

$$\begin{eqnarray*}{\left[\hat{{\boldsymbol{V}}}\cdot \displaystyle \frac{\partial \hat{{\boldsymbol{B}}}}{\partial {\boldsymbol{r}}}\right]}_{i} & = & \displaystyle \frac{1}{B}{\widehat{V}}_{j}\displaystyle \frac{\partial {B}_{i}}{\partial {x}_{j}}-\displaystyle \frac{{B}_{i}}{{B}^{2}}{\widehat{V}}_{j}\displaystyle \frac{\partial B}{\partial {x}_{j}}\\ & = & \displaystyle \frac{1}{B}{\widehat{V}}_{j}\left[\displaystyle \frac{1}{2}\left(\displaystyle \frac{\partial {B}_{i}}{\partial {x}_{j}}+\displaystyle \frac{\partial {B}_{j}}{\partial {x}_{i}}\right)+\displaystyle \frac{1}{2}\left(\displaystyle \frac{\partial {B}_{i}}{\partial {x}_{j}}-\displaystyle \frac{\partial {B}_{j}}{\partial {x}_{i}}\right)\right]\\ & & -\displaystyle \frac{{\widehat{T}}_{i}}{B}{\widehat{V}}_{j}\displaystyle \frac{\partial B}{\partial {x}_{j}}\\ & = & \displaystyle \frac{1}{B}\left({{ \mathcal S }}_{{ij}}{\widehat{V}}_{j}+\displaystyle \frac{1}{2}{\mu }_{0}{\epsilon }_{{jik}}{\widehat{V}}_{j}{J}_{k}\right)-\displaystyle \frac{{\widehat{T}}_{i}}{B}{\widehat{V}}_{j}\displaystyle \frac{\partial B}{\partial {x}_{j}}.\end{eqnarray*}$$

Here ${\mu }_{0}\;{\boldsymbol{J}}={\rm{\nabla }}\times {\boldsymbol{B}}$ and ${{ \mathcal S }}_{{ij}}\equiv {[{\mathbb{S}}]}_{{ij}}$ denotes the symmetric part of $\partial {\boldsymbol{B}}/\partial {\boldsymbol{r}}$. Generally, the orthogonal triad of eigenvectors of ${\mathbb{S}}$ has an arbitrary orientation with respect to $\hat{{\boldsymbol{B}}}=\hat{{\boldsymbol{T}}}$, so that ${\mathbb{S}}\cdot \hat{{\boldsymbol{V}}}={c}_{1}\;\hat{{\boldsymbol{T}}}+{c}_{2}\;\hat{{\boldsymbol{V}}}+{c}_{3}\;\hat{{\boldsymbol{T}}}\times \hat{{\boldsymbol{V}}}$, where the coefficients c₁, c₂, and c₃ generally depend on both ${\mathbb{S}}$ and $\hat{{\boldsymbol{V}}}$. Bearing also in mind that $\hat{{\boldsymbol{T}}}\cdot \hat{{\boldsymbol{V}}}\times \hat{{\boldsymbol{V}}}=\hat{{\boldsymbol{T}}}\cdot \hat{{\boldsymbol{V}}}\times \hat{{\boldsymbol{T}}}=0$ and $\hat{{\boldsymbol{T}}}\cdot \hat{{\boldsymbol{V}}}\times (\hat{{\boldsymbol{T}}}\times \hat{{\boldsymbol{V}}})=1$, we conclude that only the third term of ${\mathbb{S}}\cdot \hat{{\boldsymbol{V}}}$ and the antisymmetric part of $\partial {\boldsymbol{B}}/\partial {\boldsymbol{r}}$ contribute to $d{{ \mathcal T }}_{g}/{ds}$ to give

$$\begin{eqnarray}\displaystyle \frac{d{{ \mathcal T }}_{g}}{{ds}} & \simeq & \displaystyle \frac{1}{2\pi B}\left[{c}_{3}-\displaystyle \frac{1}{2}{\mu }_{0}\hat{{\boldsymbol{T}}}\cdot \hat{{\boldsymbol{V}}}\times (\hat{{\boldsymbol{V}}}\times {\boldsymbol{J}})\right]\\ & = & \displaystyle \frac{1}{2\pi B}\left[{c}_{3}+\displaystyle \frac{1}{2}{\mu }_{0}(\hat{{\boldsymbol{T}}}\cdot {\boldsymbol{J}})\right]\\ & = & \displaystyle \frac{{c}_{3}}{2\pi B}+\displaystyle \frac{{\mu }_{0}{J}_{\parallel }}{4\pi B}.\end{eqnarray} \tag{ 19 }$$

The contributions to ${{ \mathcal T }}_{g}$ from the c₃ and ${J}_{\parallel }$ terms are referred to hereafter as ${{ \mathcal T }}_{{gc}}$ and ${{ \mathcal T }}_{{gj}}$, respectively, i.e.,

$$\begin{eqnarray}&&{{ \mathcal T }}_{{gc}}={\displaystyle \int }_{{\boldsymbol{x}}}\displaystyle \frac{{c}_{3}}{2\pi B}\;{ds}\quad {\rm{and}}\quad {{ \mathcal T }}_{{gj}}={\displaystyle \int }_{{\boldsymbol{x}}}\displaystyle \frac{{\mu }_{0}{J}_{\parallel }}{4\pi B}\;{ds}.\end{eqnarray} \tag{ 20 }$$

Since all quantities in Equation (19) refer to the axis field line ${\boldsymbol{x}}(s)$, ${{ \mathcal T }}_{{gj}}$ is not equivalent to ${{ \mathcal T }}_{w}$ (Equation (7)) which is evaluated at the field line of interest, ${\boldsymbol{y}}(s)$. Given a flux rope with a well defined axis, the term ${{ \mathcal T }}_{{gj}}$ is identical for all field lines in the rope, while ${{ \mathcal T }}_{{gc}}$ differs for different field lines because $\hat{{\boldsymbol{V}}}$ differs. In the limit $\varepsilon \to 0$ the twist number ${{ \mathcal T }}_{w}$ approaches ${{ \mathcal T }}_{{gj}}$, so that

$$\begin{eqnarray}&&\underset{\varepsilon \to 0}{\mathrm{lim}}{{ \mathcal T }}_{w}(\varepsilon )={{ \mathcal T }}_{g}-{\displaystyle \int }_{{\boldsymbol{x}}}\displaystyle \frac{{c}_{3}}{2\pi B}\;{ds}.\end{eqnarray} \tag{ 21 }$$

It is now clear that ${{ \mathcal T }}_{w}$ can serve as a reliable approximation of ${{ \mathcal T }}_{g}$ only when two conditions are satisfied. First, the field line must be sufficiently close to the axis such that ${J}_{\parallel }/B$ on ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$ are approximately equal (as also required in Berger & Prior 2006). Second, the contribution from ${\mathbb{S}}$ proportional to c₃ must be negligible. Particularly, in cylindrical symmetry, ${B}_{r}=0$, ${B}_{\phi }={B}_{\phi }(r)$, and ${B}_{z}={B}_{z}(r)$, all elements of ${\mathbb{S}}$ vanish identically except

$$\begin{eqnarray*}&&{S}_{r\phi }\equiv \displaystyle \frac{1}{2}r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{B}_{\phi }}{r}\right).\end{eqnarray*}$$

For a smooth distribution of current density ${J}_{\parallel }(r)$ this term vanishes at the axis, so that ${c}_{3}=0$. For example, in a uniform-α force-free flux rope (Lundquist 1950), ${B}_{\phi }={B}_{0}{J}_{1}(\alpha r)$, hence ${S}_{r\phi }=-\alpha {B}_{0}{J}_{2}(\alpha r)/2$, where J₁ and J₂ are Bessel functions of the first kind. Also, in a uniformly twisted flux rope (Gold & Hoyle 1960), ${B}_{\phi }={br}/(1+{b}^{2}{r}^{2})$, hence ${S}_{r\phi }=-{b}^{3}{r}^{2}/{(1+{b}^{2}{r}^{2})}^{2}$. In both cases, ${S}_{r\phi }$ = 0 at the axis. Thus, one may use the ratio between the two terms in Equation (19), $2{c}_{3}/{\mu }_{0}{J}_{\parallel }$, to evaluate locally how close a flux rope is to cylindrical symmetry.

Of course, Equation (13) can be directly applied to a cylindrical flux tube. In this case, for all s, $\hat{{\boldsymbol{T}}}={\hat{{\boldsymbol{e}}}}_{z}$, $\hat{{\boldsymbol{V}}}={\hat{{\boldsymbol{e}}}}_{r}$, and

$$\begin{eqnarray*}&&\displaystyle \frac{d\hat{{\boldsymbol{V}}}}{{dz}}=\displaystyle \frac{d\phi }{{dz}}\;{\hat{{\boldsymbol{e}}}}_{\phi }.\end{eqnarray*}$$

From the field line equation in cylindrical coordinates, ${dr}/{B}_{r}={rd}\phi /{B}_{\phi }={dz}/{B}_{z}={ds}/B$, we have

$$\begin{eqnarray}&&{{ \mathcal T }}_{g}=\displaystyle \frac{1}{2\pi }\int d\phi =\displaystyle \frac{1}{2\pi }\int \displaystyle \frac{{B}_{\phi }(r)}{{{rB}}_{z}(r)}{dz}.\end{eqnarray} \tag{ 22 }$$

By cylindrical symmetry, the radial distance of the field line is independent of z, so that the classical expression for the winding number N about the z axis in a cylinder of length L_z is obtained,

$$\begin{eqnarray}&&{{ \mathcal T }}_{g}(r)=\displaystyle \frac{1}{2\pi }\displaystyle \frac{{L}_{z}{B}_{\phi }(r)}{{{rB}}_{z}(r)}=N(r).\end{eqnarray} \tag{ 23 }$$

The twist number ${{ \mathcal T }}_{w}$ can straightforwardly be computed for any field line without resorting to the geometry of an MFR, providing a convenient means to characterize and quantify magnetic configurations through maps of ${{ \mathcal T }}_{w}$ ("twist maps"). This must be done with some caution, however. The considerations of Appendix C.1 show that ${{ \mathcal T }}_{w}$ represents the classical meaning of twist—winding of field lines about an axis field line—only under certain conditions. An MFR structure possessing some degree of coherence must exist. In this case, ${{ \mathcal T }}_{w}$ is a close approximation of its twist only in the vicinity of the magnetic axis and if the MFR is approximately cylindrically symmetric. The latter condition can be quantified from Equation (19) as ${c}_{3}\ll {\mu }_{0}{J}_{\parallel }/2$.

${{ \mathcal T }}_{w}$ is also useful in locating the magnetic axis of an MFR, the nontrivial requirement for the computation of ${{ \mathcal T }}_{g}$. From Equation (7), ${{ \mathcal T }}_{w}$ may peak at any distance r from the axis, depending on the radial profile of ${J}_{\parallel }/B$. For a force-free field this is $\alpha (r)$. If the MFR possesses some degree of cylindrical symmetry, ${{ \mathcal T }}_{w}$ will have a local extremum at its axis, except in the special case that ${J}_{\parallel }/B$ is uniform in the vicinity of the axis. The uniform-α force-free rope (Lundquist 1950) has the current distribution ${J}_{z}(r)={\mu }_{0}^{-1}\alpha {B}_{0}{J}_{0}(\alpha r)$, where J₀ is the zeroth-order Bessel function. Thus, ${{ \mathcal T }}_{w}$ will peak at the axis of a force-free flux rope if ${J}_{\parallel }(r)$ is more peaked than this function and will have a minimum at the axis in the opposite case.

To demonstrate the usefulness of ${{ \mathcal T }}_{w}$ in locating the magnetic axis of an MFR as well as its limitation in approximating ${{ \mathcal T }}_{g}$, we consider again the cylindrically symmetric flux ropes, in which ${{ \mathcal T }}_{g}=N$. For the Lundquist flux rope ${{ \mathcal T }}_{w}$ approaches N at its axis,

$$\begin{eqnarray*}&&\underset{r\to 0}{\mathrm{lim}}N(r)=\underset{r\to 0}{\mathrm{lim}}\displaystyle \frac{{B}_{0}{J}_{1}(\alpha r)}{2\pi {{rB}}_{0}{J}_{0}(\alpha r)}{L}_{z}=\displaystyle \frac{\alpha }{4\pi }{L}_{z},\end{eqnarray*}$$

but in general,

$$\begin{eqnarray}\displaystyle \frac{N(r)}{{{ \mathcal T }}_{w}(r)} & = & \displaystyle \frac{{B}_{\phi }{L}_{z}/2\pi {{rB}}_{z}}{\alpha {{BL}}_{z}/4\pi {B}_{z}}=\displaystyle \frac{2{B}_{\phi }}{\alpha {rB}}\\ & = & \displaystyle \frac{2{J}_{1}(\alpha r)}{\alpha r\surd {J}_{0}^{2}(\alpha r)+{J}_{1}^{2}(\alpha r)}\lt 1.\end{eqnarray} \tag{ 24 }$$

For example, ${{ \mathcal T }}_{w}$ overestimates N by 25% at $r\approx 2.5{\alpha }^{-1}$ in this rope.

In the uniformly twisted flux rope, we have $N={{bL}}_{z}/2\pi$, where the constant $2\pi /b$ is the axial length of one field line turn, ${L}_{z}/N$, given by $b=d\phi /{dz}={B}_{\phi }/{{rB}}_{z}=2\pi N/{L}_{z}$. With ${B}_{\phi }={br}/(1+{b}^{2}{r}^{2})$ and ${B}_{z}=1/(1+{b}^{2}{r}^{2})$, we have

$$\begin{eqnarray*}{{ \mathcal T }}_{w} & = & {\displaystyle \int }_{L}\displaystyle \frac{({\rm{\nabla }}\times {\boldsymbol{B}})\cdot {\boldsymbol{B}}}{4\pi {B}^{2}}\;{ds}={\displaystyle \int }_{{L}_{z}}\displaystyle \frac{({\rm{\nabla }}\times {\boldsymbol{B}})\cdot {\boldsymbol{B}}}{4\pi {{BB}}_{z}}\;{dz}\\ & = & \displaystyle \frac{1}{4\pi }\displaystyle \frac{2b}{\surd 1+{b}^{2}{r}^{2}}{L}_{z},\end{eqnarray*}$$

and hence

$$\begin{eqnarray}&&\displaystyle \frac{N}{{{ \mathcal T }}_{w}}=\sqrt{1+{b}^{2}{r}^{2}}.\end{eqnarray} \tag{ 25 }$$

Again, ${{ \mathcal T }}_{w}$ matches N at the axis, but away from the axis ${{ \mathcal T }}_{w}$ underestimates N, contrary to the case of the Lundquist rope.

We have further checked the approach of ${{ \mathcal T }}_{w}$ and N at the flux rope axis and the location of the peak ${{ \mathcal T }}_{w}$ for an approximately force-free Titov–Démoulin flux rope equilibrium, using two different choices for the toroidal current density ${J}_{{\rm{t}}}(r)$, one roughly uniform as in the original construction of the equilibrium (Titov & Démoulin 1999), the other with a distribution of ${J}_{{\rm{t}}}(r)$ strongly peaked at r = 0. By construction, each small section of this toroidal flux rope possesses approximate cylindrical symmetry, especially for large aspect ratio. The computation was performed on a discrete grid, similar to the situation of a numerical extrapolation or an MHD simulation, in a vertical cut through the flux rope apex. For an aspect ratio of merely 3.4, ${{ \mathcal T }}_{w}$ (Equation (7)) and N (Equation (23)) are found to agree to within 5% at the axis, where they also peak in the case of the peaked ${J}_{{\rm{t}}}(r)$. For the roughly uniform ${J}_{{\rm{t}}}(r)$, both peak at the periphery of the current channel in the center of the flux rope and are minimal at the magnetic axis, with similar radial profiles.

Thus, assuming a coherent MFR, which possesses at least an approximate cylindrical symmetry, and excluding the special case of uniform ${J}_{\parallel }(r)/B(r)$, the magnetic axis is located at the local extremum of the $| {{ \mathcal T }}_{w}|$ map, which either is a peak, or a dip enclosed by a ring of high $| {{ \mathcal T }}_{w}|$ values. The flux ropes in the NLFFFs studied in this paper conform to the former case.

In the following we compare the different means of obtaining a twist number, using the flux rope shown in Figure 4. These are ${{ \mathcal T }}_{g}$ from Equation (13), its approximation broken into two terms, ${{ \mathcal T }}_{{gc}}$ and ${{ \mathcal T }}_{{gj}}$, from Equation (20), ${{ \mathcal T }}_{w}$ from Equation (7), ${{ \mathcal T }}_{w}^{*}$ from Equation (8), and the simplified estimate of the latter, ${{ \mathcal T }}_{\alpha }=({\alpha }_{\mathrm{start}}+{\alpha }_{\mathrm{end}})L/8\pi$, obtained by multiplying the average value of α at the field line footpoints with the field line length L.

Figures 14(a) and (b) examine the relationship between the exact ${{ \mathcal T }}_{g}$ (Equation (13)) and the approximate one (integration of Equation (19) along the axis), as well as the individual terms in the latter (Equation (20)). The uniform value of ${{ \mathcal T }}_{{gj}}$ is indicated by the horizontal line in (a). Both ${{ \mathcal T }}_{{gc}}$ and ${{ \mathcal T }}_{{gj}}$ contribute significantly to ${{ \mathcal T }}_{g}$, however, none of them is individually correlated with ${{ \mathcal T }}_{g}$, and near the axis they are of opposite sign. Their sum is highly correlated with ${{ \mathcal T }}_{g}$ for all field lines of the inner two sets (blue and red), indicating that the approximation of small ε underlying Equation (19) is here satisfied up to about one third of the rope's minor radius. As expected from the large radial distance and irregular appearance in the field line plot, there's no correlation for the remote (green) field lines.

Zoom In Zoom Out Reset image size

Figure 14(c) shows that ${{ \mathcal T }}_{w}$ and ${{ \mathcal T }}_{g}$ are in only moderate agreement, due to two reasons. One first sees that the blue points agree rather well with ${{ \mathcal T }}_{{gj}}$, whereas the red points here agree only moderately. This reflects the fact that the current density in our flux rope is considerably peaked at the axis (Figure 4(e)), so that ${{ \mathcal T }}_{w}$, which is based purely on the local ${J}_{\parallel }(r)/B(r)$, begins to differ from ${{ \mathcal T }}_{{gj}}$ at a smaller radial distance than the approximate ${{ \mathcal T }}_{g}$ (Equation (19)) from the exact one (Equation (13); Figure 14(b)). Second, even the blue points do not agree perfectly with ${{ \mathcal T }}_{g}$, with ratios in the range ${{ \mathcal T }}_{w}/{{ \mathcal T }}_{g}\approx 1.3{\rm{\mbox{--}}}2$. This implies that the term ${{ \mathcal T }}_{{gc}}$ is not negligible for our flux rope, i.e., the rope as a whole is not close to cylindrical symmetry, consistent with its appearance in Figure 4. We will address this in detail below.

The different forms of twist number not referring to an axis are compared in Figures 14(d) and (e). ${{ \mathcal T }}_{w}$ and ${{ \mathcal T }}_{w}^{*}$ are found to be highly correlated, indicating that the force-free condition is very well satisfied in the volume of the flux rope. On the other hand, the simplified variant, ${{ \mathcal T }}_{\alpha }$, shows significant differences. Thus, deviations from force-freeness in the extrapolated field are particularly high at the footpoints of the field lines. This problem is also illustrated by the difference of ${\alpha }_{\mathrm{start}}$ and ${\alpha }_{\mathrm{end}}$, which can be substantial (Figure 14(f)), and has also been noted by Inoue et al. (2011) and Chintzoglou et al. (2015). Further, we found that the correlation between T_α and T_g becomes much worse when α at the footpoints is obtained by ${J}_{z}/{B}_{z}$ rather than by J/B. This may be due to the large uncertainties associated with the transverse field components B_x and B_y in general and those associated with the weak B_z near the PIL in particular.

The local density of ${{ \mathcal T }}_{{gj}}$ and ${{ \mathcal T }}_{{gc}}$ is further examined in Figure 15. ${c}_{3}/2\pi B$ is minimal and much smaller than ${\mu }_{0}{J}_{\parallel }/4\pi B$ in the middle section of the flux rope, which thus is nearly cylindrically symmetric. The fact that ${c}_{3}\approx 0$ in the middle section indicates clearly that the magnetic axis of the forming flux rope was correctly located by the maximum of $| {{ \mathcal T }}_{w}|$ in our twist maps.

Zoom In Zoom Out Reset image size