# Numerical Study on In-Situ Prominence Formation by Radiative Condensation in the Solar Corona

###### Abstract

We propose an in-situ formation model for inverse-polarity solar prominences and demonstrate it using self-consistent 2.5-dimensional magnetohydrodynamics simulations, including thermal conduction along magnetic fields and optically thin radiative cooling. The model enables us to form cool dense plasma clouds inside a flux rope by radiative condensation, which is regarded as an inverse-polarity prominence. Radiative condensation is triggered by changes in the magnetic topology, i.e., formation of the flux rope from the sheared arcade field, and by thermal imbalance due to the dense plasma trapped inside the flux rope. The flux rope is created by imposing converging and shearing motion on the arcade field. Either when the footpoint motion is in the anti-shearing direction or when heating is proportional to local density, the thermal state inside the flux rope becomes cooling-dominant, leading to radiative condensation. By controlling the temperature of condensation, we investigate the relationship between the temperature and density of prominences and derive a scaling formula for this relationship. This formula suggests that the proposed model reproduces the observed density of prominences, which is 10–100 times larger than the coronal density. Moreover, the time evolution of the extreme ultraviolet emission synthesized by combining our simulation results with the response function of the Solar Dynamics Observatory Atmospheric Imaging Assembly filters agrees with the observed temporal and spatial intensity shift among multi-wavelength during in-situ condensation.

^{†}

^{†}slugcomment: Not to appear in Nonlearned J., 45.

## 1 Introduction

Solar prominences are cool dense plasma clouds that form in the hot tenuous corona. Dense plasmas are sustained at the dipped (concave-up) structures of coronal arcade fields (Kippenhahn & Schlüter, 1957; Kuperus & Raadu, 1974). Prominences can be classified into two groups — normal- and inverse-polarity prominences — depending on their interior magnetic polarity with respect to the overlying arcade field. In normal-polarity prominences, the magnetic polarity matches that of the overlying arcade field. Normal-polarity prominences are typically modeled by the Kippenhahn–Schlüter (KS) model. In inverse-polarity prominences, the interior magnetic polarity opposes that of the overlying arcade field. This prominence group is described by the Kuperus–Raadu (KR) model, in which the lower half of the flux rope indicates the inverse-polarity. Both normal- and inverse-polarity prominences are observed in the solar atmosphere (Leroy et al., 1984).

An unsolved problem regarding prominences is the origin of their cool dense plasmas. Though several models have been proposed, the actual mechanism by which prominences form remains debatable (Mackay et al., 2010). One candidate mechanism for prominence formation is radiative condensation (Antiochos et al., 1999). Several models for triggering radiative condensation have been proposed and demonstrated by numerical simulations. One is the evaporation–condensation process (Mackay et al., 2010), in which evaporated flows from the chromosphere lead to radiative condensation in the coronal loops. The mechanism and constraints of this model have been investigated using one-dimensional hydrodynamics simulations (Karpen et al., 2006; Luna et al., 2012) and two-dimensional self-consistent magnetohydrodynamics (MHD) simulations (Xia et al., 2012). These studies have shown the formation of normal-polarity prominences categorized by the KS model, whereas Keppens & Xia (2014) showed that the evaporation–condensation model eventually evolved into the funnel prominence in which the flux rope structure is embedded.

Choe & Lee (1992) demonstrated the formation of a normal-polarity prominence by a different process. In their two-dimensional MHD simulations, the shearing motion at the footpoints of the potential arcade fields leads to radiative condensation due to the expansion of the arcade fields. Linker et al. (2001) demonstrated the formation of inverse-polarity prominences. In their simulations, the chromospheric plasmas are levitated as a result of flux rope formation by the magnetic cancellation and condensation occurs.

Condensation models triggered by evaporation or levitation rely on the injection of chromospheric plasmas. On the other hand, a recent observation of a prominence formation event showed no clear evidence of chromospheric plasma injection, and claimed the presence of in-situ condensation (Berger et al., 2012). Xia et al. (2014a) modeled in-situ condensation inside a flux rope by three-dimensional simulation, and showed that the features of prominences are well reproduced. Because they used the mechanically stable flux rope model obtained from the isothermal simulation (Xia et al., 2014b) and started their simulation with an ad hoc thermal imbalance, the process required to trigger the radiative condensation is still unclear.

In this study, we model the process from flux rope formation to radiative condensation using self-consistent 2.5-dimensional MHD simulations, including thermal conduction and radiative cooling. By comparing the results of multiple simulations with different settings regarding the formation model of flux rope and the coronal heating model, we discuss the condition for triggering radiative condensation and the possible relationship between the temperature and density of the prominence condensations.

Section 2 introduces our prominence formation model, and Section 3 describes the numerical simulation setup. The simulation results are presented in Section 4. The condition for triggering radiative condensation and a scaling formula from the possible temperature–density relationships are discussed in Section 5.

## 2 Possible processes for triggering in-situ condensation

The schematic picture of our model for triggering in-situ radiative condensation is presented as follows. Initially, the coronal arcade exists in thermal equilibrium (Fig. 1 (a)). The arcade field is transformed into a flux rope structure by converging and shearing motion (van Ballegooijen & Martens, 1989; Martens & Zwaan, 2001). The relatively dense plasmas in the lower corona due to stratification are trapped inside closed loops of the flux rope, where they are elevated into the upper corona (Fig. 1 (b)). The dense plasmas increase the radiative cooling in the flux rope, leading to thermal imbalance. Because the closed field lines are disconnected from the exterior region, the heat loss by radiative cooling inside the flux rope is not conducted outward, and radiative condensation can be triggered. (Fig. 1 (c)).

## 3 Numerical Settings

To demonstrate the proposed model, we perform 2.5-dimensional MHD simulations, including thermal conduction along the magnetic field lines, radiative cooling, and gravity in a Cartesian (,) domain. To investigate the necessary condition for triggering radiative condensation, we test different direction of the converging and shearing motion for the formation of a flux rope and different coronal heating models.

### 3.1 Initial and boundary conditions

The corona model is set up in a rectangular box, in which the - and -axes are horizontal and vertical, respectively, and the -axis is orthogonal to the – plane. The corona is stratified under uniform temperature () and gravity () conditions,

(1) | |||||

(2) |

where , is Boltzmann’s constant, and is the mean molecular mass. We set with being the proton mass. The force-free arcade field iss described as

(3) | |||||

(4) | |||||

(5) |

where is the field strength at the footpoint, is the width, and is the magnetic scale height of the arcade field. Initially, the system exists in mechanical equilibrium.

The left and right boundaries are subjected to symmetric (for ) or anti-symmetric (for ) boundary conditions. A free boundary condition is applied to the top. In the region below , the converging and shearing motions are introduced. We test three types of the footpoint motions. One is the converging motion without shearing; the others are converging motion with shearing that increases the magnetic shear of the arcade field or with anti-shearing that decreases the shear of the arcade field. In all the cases, the velocity components and within this region are set as follows,

(6) | |||||

(7) |

(8) | |||||

(9) | |||||

(10) |

where , , and . In the case of no shearing motion, . The shearing and the anti-shearing motions are set as follows,

(11) | |||||

(12) |

where the plus and minus signs represent the shearing and anti-shearing motions, respectively. The magnetic fields are computed with the induction equation by coupling with the given converging and shearing motions. The magnetic fields at the bottom boundary are fixed to the initial state. The gas pressure and density are fixed by assuming hydrostatic equilibrium at a constant temperature of .

### 3.2 Basic equations

The basic equations are

(13) |

(14) |

(15) |

(16) |

(17) |

(18) |

(19) |

where is the coefficient of thermal conduction, is a unit vector along the magnetic field, and is the gravitational acceleration and is the magnetic diffusion rate. The temperature is computed by the following equation of state,

(20) |

For fast reconnection, we adopt the following form of the anomalous resistivity (e.g. Yokoyama & Shibata, 1994),

(21) | |||||

(22) |

where and . We restrict to . in Eq. (14) represents the net cooling term, which is expressed as

(23) |

where is the number density, is the radiative loss function of optically thin plasma, is the background heating rate, and represents ohmic heating. We adopt a simplified radiative loss function in Fig. 2 (a) (Hildner, 1974) and simulate two different coronal heating models. In one model, the heating rate depends on the local magnetic energy density (magnetic pressure) at temperatures above the cutoff temperature and is balanced out by the cooling rate when it falls below ,

(24) | |||||

(25) |

Initial thermal equilibrium required that

(26) |

where . Selecting , we obtain (i.e., a constant). Equation (25) prevents the temperature from decreasing below . In the other model, the heating rate depends on the local density and height at temperatures above and balances with the cooling rate below the cutoff temperature .

(27) | |||||

(28) |

For the initial thermal equilibrium, we have where . Note that ohmic heating is always much smaller than background heating in these settings. is one of the model parameters. For each heating and shearing model, we vary as MK, MK, and MK. The investigated cases are summarized in Table 1.

### 3.3 Numerical scheme and grid size

The numerical scheme is based on the CIP-MOCCT method (Kudoh et al., 1999). To prevent numerical oscillations in the advection term (left-hand side of Eqs. (13)–(16)), we apply the rational CIP scheme (Xiao et al., 1996), in which the rational function replaces the cubic function as an interpolation function of CIP. The thermal conduction, radiative cooling, and heating terms are included as source terms in non-advection phases of CIP method. Thermal conduction is explicitly solved by a second-order accurate slope-limiting method (Sharma & Hammett, 2007).

The required grid size for preventing numerical thermal instability can be estimated by the formula presented in Koyama & Inutsuka (2004). When applying this formula to our simulation, we assume that the density of the condensations is 10 times the coronal density,

(29) |

where is the grid size, is the Field length, is the cutoff temperature (also assume to be the temperature of condensations), and is the assumed density of condensations. Figure 2 (b) shows the required grid size against temperature. Under the constraint of Eq. (29), the grid size is determined to be 120 km, 60 km, and 30 km in the cases of 0.5 MK, 0.4 MK, and 0.3 MK, respectively. The horizontal grid size is uniform. The vertical grid size is uniform up to and increases to at to ensure a sufficiently distant upper boundary.

## 4 Results

The right-most column of Table 1 shows our results in the presence of the radiative condensation. We find that the necessary condition for the radiative condensation is either anti-shearing motion or the heating proportional to the local density. Panels (a)–(d) of Fig. 3 are snapshots of the time evolution and Fig. 4 is the final state in case M1, which is a the typical case for the radiative condensation. The initial state exists in mechanical and thermal equilibrium (Fig. 3 (a)). Converging motion triggers reconnection above the PIL at , and a flux rope is formed (Fig. 3 (b)). As the reconnection proceeds, the relatively dense plasmas in the lower corona are trapped inside the flux rope and lifted to the upper corona. The dense plasmas enhances radiative cooling in the flux rope, which overwhelms the background heating and leads to thermal imbalance (Fig. 3 (c)). Because thermal conduction works only along closed magnetic field lines of the flux rope, it can not compensate for the thermal imbalance. Consequently, radiative condensation is triggered, and the cool dense plasmas accumulates in the lower part of the flux rope (Fig. 3 (d)). Figure 4 (b) shows profiles of temperature, number density, and horizontal component of the magnetic field along the -axis in the quasi–static state (Fig. 4 (a)). The temperature of the prominence corresponding to is approximately in this case. The heating settings (Eqs. (25) and (28)) maintains the temperature of the cool dense plasmas at around . According to Fig. 4 (b), the magnetic polarity of the cool dense region opposes that of the initial arcade field. Figures 4 (a) and (b) clearly show that radiative condensation eventually establishes inverse-polarity prominence in the flux rope. Also, the high-density region (warm color region in the density plot of Fig. 4 (a)) resembles the tower prominence (Berger et al., 2012), and the low-density region surrounding the prominence (dark region in the same plot) can be regarded as the coronal cavity. Figure 5 shows each force at a time of 4250s. The mechanical quasi-static state is realized in the vertical direction. The upward magnetic tension force is balanced with the downward pressure gradient force and gravity, i.e., the prominence mass is supported mainly by the magnetic tension. This force balance system is also found in the parametric study by adiabatic MHD simulations of Hillier & van Ballegooijen (2013). Our result mostly corresponds to their case (e) of model 1 regarding the force balance, whereas the prominence density and plasma in our simulation are lower and higher than those in Hillier & van Ballegooijen (2013), respectively. They revealed that the magnetic tension force is given by the vertical stretching of the flux rope, and the curvature of the field lines inside the prominence can be smaller because the sufficiently large magnetic tension comes from the larger field strength by compression of the prominence mass. Our simulation realizes these features of prominence support as shown in Fig. 4 (a). Figures 6 shows the time evolution of case M5, which is a typical case with no radiative condensation. The flux rope is heated-up rather than cooled-down in this case.

case | Heating | Shearing | (MK) | Condensation |
---|---|---|---|---|

M1 | 0.3 | Yes | ||

M2 | 0.4 | Yes | ||

M3 | 0.5 | Yes | ||

M4 | 0 | 0.5 | No | |

M5 | 0.5 | No | ||

D1 | 0.3 | Yes | ||

D2 | 0.4 | Yes | ||

D3 | 0.5 | Yes | ||

D4 | 0 | 0.5 | Yes | |

D5 | 0.5 | Yes |

## 5 Discussion

As shown in Table 1, the radiative condensation occurs either when the anti-shearing motion is imposed or when the heating model is proportional to the density, i.e., . The radiative condensation is triggered by the cooling-dominant thermal imbalance in the flux rope composed of the thermally isolated closed loops. The interior of the flux rope is relatively dense because of the converging motion as well as the initial stratification. Because the radiative cooling inside the flux rope is enhanced by the dense plasma, the cooling-dominant thermal imbalance is reasonable if the heating changes only slightly. Because we set the heating depending on the magnetic pressure (cases M1–M5, hereafter called M-cases) and local density (cases D1–D5, hereafter called D-cases), the change in heating rate is not negligible. In some cases (M5 and M6), a heating-dominant thermal imbalance is created and no condensations occur. In the following subsections, we discuss the mechanism for creating the cooling/heating-dominant thermal imbalance in the M- and D-cases. Then, our model is compared with previous theoretical studies and observations.

### 5.1 Thermal imbalance in M-cases

In M-cases where the heating rate is proportional to the magnetic pressure, the anti-shearing and shearing motion causes a larger discrepancy between the cooling and heating rates, establishing the cooling- and heating-dominant thermal imbalances, respectively. Figure 7 shows the cooling and heating rates along the -axis at time= in case M3–M5. In case M3 (anti-shearing case), the cooling rate is the highest of the three cases, whereas the heating rate is the lowest. Case M5 (shearing case) has the opposite tendency, where the cooling rate is the lowest and the heating rate is the highest. We discuss how these tendencies are established from a perspective of the evolution of the magnetic field and its influence on the cooling and heating rates.

In the anti-shearing case, the flux rope is pushed upward by the reconnection and confined by the downward magnetic tension of the overlying arcade field. Figure 8 shows the vertical tension force along the -axis at time= in case M3 (anti-shearing), M4 (no shearing), and M5 (shearing). The downward magnetic tension force in the anti-shearing case is larger than those in the shearing and no shearing cases. The reason the anti-shearing case has the larger downward tension force is as follows: the shear angle is related to the height of the arcade field as

(30) |

where is the arcade height, is the shear angle against the positive -axis in the -plane, is the arcade width, and we assume a linear force-free arcade field for simplicity. Equation (30) indicates that the reduction of magnetic shear by the anti-shearing motion makes the arcade field shorter. Due to the shrink of the arcade field, the flux rope in the anti-shearing case experiences a larger downward magnetic tension force.

The flux rope in the anti-shearing case is pinched by the overlying arcade field and pushed by the upward magnetic tension of the reconnected loops, leading to a density enhancement by compression. Figure 9 shows the time evolution of area and average density inside the dash-dotted line in Fig. 6. The average density is computed by dividing the total mass inside the dash-dotted line by the corresponding area. The total mass is almost the same between these three cases and is conserved after the flux rope is detached from the bottom boundary. In Fig. 9 (b), the anti-shearing case shows the enhancement of the average density. This result indicates that the density enhancement inside the flux rope in the anti-shearing case is due to the compression by the magnetic tension as well as the converging motion.

The anti-shearing motion lead to the higher cooling rate as well as the lower heating rate. The anti-shearing motion reduces the magnetic component parallel to the PIL, corresponding to in our simulations, and decreases the magnetic pressure, leading to the lowest heating rate. Thus, the interior of the flux rope in case M3 suffers from unstoppable cooling most efficiently. In case M5 (shearing case), on the other hand, the heating rate increases as the magnetic pressure increases, and the enhancement of the cooling rate is too small to overwhelm the heating rate. Eventually, the shearing motion is likely to establish the heating-dominant thermal imbalance, resulting only in the formation of hot flux rope. These results suggest that when the coronal heating depends on the magnetic energy density, the continuous shearing motion can not trigger the radiative condensation, whereas the temporary anti-shearing motion can indeed trigger it, leading to the prominence formation.

### 5.2 Thermal imbalance in D-cases

From the results of the D-cases in which the heating is proportional to the density, it is found that the cooling-dominant thermal imbalance can be achieved only by the difference between the heating and cooling rates in the dependence on the density. Figure 10 shows the cooling and the heating rates in the cases D3–D5. The cooling rate in case D3 (anti-shearing case) is highest because the density enhancement is largest in the anti-shearing case, as discussed in the previous subsection. Unlike case M3, however, the heating rate in case D3 is also highest. The combination of the higher cooling rate and the higher heating rate can be a disadvantage to enhancing the thermal imbalance. What actually sets the cooling-dominant thermal imbalance in D-cases is the different density dependence between the heating () and cooling terms (). The increment of the heating rate against the density increase is always smaller than that of the cooling rate. Consequently, the cooling overwhelms the heating, and the radiative condensation is triggered in all D-cases. These results suggest that if the coronal heating is only proportional to the density, the flux ropes in the corona are likely to cool down. The exact modeling of the coronal heating is required to confirm this mechanism.

### 5.3 Comparison with other theoretical studies

In evaporation–condensation models (Antiochos et al., 1999; Karpen et al., 2006), radiative condensation is triggered by the transport of dense plasma from lower altitudes. Our model is consistent with this mechanism, but it adopts a different approach to suppress the thermal conduction effects. In the evaporation–condensation model, thermal conduction effects are suppressed by preparing long magnetic field lines, whereas in our model, they are suppressed by closed loop formations. In other words, our model altered the magnetic field topology. Both our study and that of Choe & Lee (1992) assume that photospheric motion triggers the prominence formation. However, the direction of motion differs between the studies. In Choe & Lee (1992), radiative condensation arises by shearing motion along the PIL and subsequent expansion of the arcade field. In our study, converging motion towards the PIL and subsequent formation of closed magnetic fields (flux rope) are essential precursors of radiative condensation. Moreover, whereas their models reproduce normal-polarity prominences, our models form inverse-polarity prominences. The formation of the inverse-polarity prominence is found both in the results of our model and that of Linker et al. (2001). The difference is that our model does not require chromospheric plasma injection to trigger the radiative condensation; hence, our model reproduce the in-situ condensation. Xia et al. (2014a) performed three-dimensional simulation of the in-situ condensation inside the flux rope. The flux rope system in Xia et al. (2014a) was created by the converging and shearing motion in the isothermal simulation of Xia et al. (2014b). Because the strategy used to establish cooling-dominant thermal imbalance in Xia et al. (2014a) was the parameterized heating, it was still unclear why the thermal imbalance was created inside the flux rope. We ensure that the cooling-dominant thermal imbalance can be created inside the flux rope by two different mechanisms: the anti-shearing motion or the different dependence on the local density between the heating and cooling rates. We find that the direction of the shear motion strongly affected the thermal processes when the heating depends on the magnetic energy density, and it controls the presence of in-situ condensation. When the heating is proportional to the density, on the other hand, the cooling-dominant thermal imbalance is always established by the difference between the cooling and heating terms in the dependence on density, leading to radiative condensation.

The formation mechanism through the anti-shearing motion smoothly connects to some of the eruptive models. These include instability models such as kink and torus instability (Kliem et al., 2004; Török & Kliem, 2005; Kliem & Török, 2006; Fan & Gibson, 2007). The kink instability is likely to be triggered because the anti-shearing motion increases the twist of the flux rope. The torus instability also fit because the reduction of shear can alleviate the critical height of torus instability. The other model that fit with our mechanism is the reversed shear model of Kusano et al. (2004). We stopped the anti-shearing motion before the shear reversal. If we continued to impose the anti-shearing motion, the flux rope formation, radiative condensation, and eruption would be successively occurred even without the converging motion. Both the shearing and anti-shearing motions have been considered as the possible process of eruptions (Kusano, 2002; Amari et al., 2003). Our results suggest that continuous shearing motion results in the eruption of the hot flux rope corresponding to flare eruption, whereas the anti-shearing motion causes the radiative condensation of the prominence formation. The eruptive mechanism coupling with flux emergence is also the subject to investigate (Chen & Shibata, 2000; Shiota et al., 2005; Kaneko & Yokoyama, 2014). The results of our simulation with the good mechanical balance can be directly adopted as the initial condition for these eruptive studies. The model including the dense materials of the prominence may reveal the more detailed structure of coronal mass ejection such as the dense core and its evolution.

### 5.4 Comparison with observation

We compare the simulated process of in-situ condensation with the observation by Berger et al. (2012) by synthetic emission through the extreme ultraviolet (EUV) filters of the Solar Dynamics Observatory Atmospheric Imaging Assembly (SDO/AIA). The emission of a certain wavelength channel is expressed as

(31) |

where , , , and represent the pixel response to the photon flux, electron number density, temperature response function of the AIA filters, and the distance along the line-of-site (LOS), respectively. The temperature response functions are obtained from aia_get_response.pro in the Solar Software Library. To compute the synthetic emission, we choose the -axis as the LOS. In our 2.5-dimensional simulation, the physical variables are assumed to be constant along -axis; hence, the integration along the LOS is altered to the multiplication of an arbitrary constant as follows,

(32) |

Figure 11 is the time evolution of the synthetic emission of each EUV filter. Note that Fig. 11 covers the time only before the temperature reaches the cutoff. Figure 11 clearly shows that the areas of strong emission progressively shift from the 193 channel, through 171 channel, to 131 channel both in time and altitude. The emission shift was claimed to be evidence of in-situ condensation in Berger et al. (2012), and our simulation has validated their argument. The shift from the 211 channel to the 193 channel is not so clear because the initial coronal temperature is set to 1MK. The shift in altitude in Fig. 11 results from the descent of the flux rope and high density plasma. It can also be realized by some other mechanisms such as magnetic reconnection in the weighted dip (Liu et al., 2012; Keppens & Xia, 2014), or cross-field slippage (Low et al., 2012a, b). These mechanisms are not realized in our simulations because we set the cutoff temperature such that the simulated prominence density can not be sufficiently large.

Figure 12 shows the relationship between the temperature and the maximum density of the condensations obtained from our results. The lower the condensation temperature, the higher the density of the condensates. The resultant density is independent of the heating model. The temperature and density of the prominence are related by the following scaling formula,

(33) |

The scaling relationship (Eq. (33)) is plotted as a solid line in Fig. 12. This formula indicates that the dense condensation results from re-stratification along the magnetic field lines with temperature . Because the mass is conserved along the field lines and confined to the length scale of the scale height , the condensate density is inversely proportional to temperature. From Eq. (33), we concluded that our model can reproduce the observed prominence density, which is 10–100 times larger than the surrounding coronal density at typical temperatures (below ).

The direct simulation without the temperature cutoff is necessary to verify the realistic prominence density and compare DEM with observation. Including the chromosphere in our model is also future work. Due to the thermal conduction from the corona to chromosphere, the background heating rate must be larger than the cooling rate. When the initial coronal density is lower, the magnitude of the initial thermal imbalance is lower, and the present model without the initial thermal imbalance will be validated. It is the next step for our model to be extended to three dimensions. Self-consistent three-dimensional MHD simulation including coronal thermal conduction is still a challenging issue. The comparison with observations including the LOS effect (Luna et al., 2012; Xia et al., 2014a) must contribute to the further understanding of the prominence formation.

## 6 Conclusion

We proposed a new in-situ formation model of inverse-polarity prominence. The model was demonstrated in 2.5-dimensional MHD simulations, involving thermal conduction and radiative cooling. Either flux rope formation by anti-shearing motion or a heating model proportional to the local density is necessary for a in-situ prominence formation. From the results of a parameter study on cutoff temperature, we derived an empirical scaling formula that related the density of the prominences to their temperature. Our study quantitatively demonstrated that our in-situ formation model can reproduce the observed density of prominences. The study also explained the observed intensity shift among multi-wavelength EUV emission of in-situ condensation.

## References

- Amari et al. (2003) Amari, T., Luciani, J. F., Aly, J. J., Mikic, Z., & Linker, J. 2003, ApJ, 585, 1073
- Antiochos et al. (1999) Antiochos, S. K., MacNeice, P. J., Spicer, D. S., & Klimchuk, J. A. 1999, ApJ, 512, 985
- Berger et al. (2012) Berger, T. E., Liu, W., & Low, B. C. 2012, ApJ, 758, L37
- Chen & Shibata (2000) Chen, P. F., & Shibata, K. 2000, ApJ, 545, 524
- Choe & Lee (1992) Choe, G. S., & Lee, L. C. 1992, Sol. Phys., 138, 291
- Fan & Gibson (2007) Fan, Y., & Gibson, S. E. 2007, ApJ, 668, 1232
- Hildner (1974) Hildner, E. 1974, Sol. Phys., 35, 123
- Hillier & van Ballegooijen (2013) Hillier, A., & van Ballegooijen, A. 2013, ApJ, 766, 126
- Kaneko & Yokoyama (2014) Kaneko, T., & Yokoyama, T. 2014, ApJ, 796, 44
- Karpen et al. (2006) Karpen, J. T., Antiochos, S. K., & Klimchuk, J. A. 2006, ApJ, 637, 531
- Keppens & Xia (2014) Keppens, R., & Xia, C. 2014, ApJ, 789, 22
- Kippenhahn & Schlüter (1957) Kippenhahn, R., & Schlüter, A. 1957, ZAp, 43, 36
- Kliem et al. (2004) Kliem, B., Titov, V. S., & Török, T. 2004, A&A, 413, L23
- Kliem & Török (2006) Kliem, B., & Török, T. 2006, Physical Review Letters, 96, 255002
- Koyama & Inutsuka (2004) Koyama, H., & Inutsuka, S.-i. 2004, ApJ, 602, L25
- Kudoh et al. (1999) Kudoh, T., Shibata, K., & Matsumoto, R. 1999, in Astrophysics and Space Science Library, Vol. 240, Numerical Astrophysics, ed. S. M. Miyama, K. Tomisaka, & T. Hanawa, 203
- Kuperus & Raadu (1974) Kuperus, M., & Raadu, M. A. 1974, A&A, 31, 189
- Kusano (2002) Kusano, K. 2002, ApJ, 571, 532
- Kusano et al. (2004) Kusano, K., Maeshiro, T., Yokoyama, T., & Sakurai, T. 2004, ApJ, 610, 537
- Leroy et al. (1984) Leroy, J. L., Bommier, V., & Sahal-Brechot, S. 1984, A&A, 131, 33
- Linker et al. (2001) Linker, J. A., Lionello, R., Mikić, Z., & Amari, T. 2001, J. Geophys. Res., 106, 25165
- Liu et al. (2012) Liu, W., Berger, T. E., & Low, B. C. 2012, ApJ, 745, L21
- Low et al. (2012a) Low, B. C., Berger, T., Casini, R., & Liu, W. 2012a, ApJ, 755, 34
- Low et al. (2012b) Low, B. C., Liu, W., Berger, T., & Casini, R. 2012b, ApJ, 757, 21
- Luna et al. (2012) Luna, M., Karpen, J. T., & DeVore, C. R. 2012, ApJ, 746, 30
- Mackay et al. (2010) Mackay, D. H., Karpen, J. T., Ballester, J. L., Schmieder, B., & Aulanier, G. 2010, Space Sci. Rev., 151, 333
- Martens & Zwaan (2001) Martens, P. C., & Zwaan, C. 2001, ApJ, 558, 872
- Sharma & Hammett (2007) Sharma, P., & Hammett, G. W. 2007, Journal of Computational Physics, 227, 123
- Shiota et al. (2005) Shiota, D., Isobe, H., Chen, P. F., et al. 2005, ApJ, 634, 663
- Török & Kliem (2005) Török, T., & Kliem, B. 2005, ApJ, 630, L97
- van Ballegooijen & Martens (1989) van Ballegooijen, A. A., & Martens, P. C. H. 1989, ApJ, 343, 971
- Xia et al. (2012) Xia, C., Chen, P. F., & Keppens, R. 2012, ApJ, 748, L26
- Xia et al. (2014a) Xia, C., Keppens, R., Antolin, P., & Porth, O. 2014a, ApJ, 792, L38
- Xia et al. (2014b) Xia, C., Keppens, R., & Guo, Y. 2014b, ApJ, 780, 130
- Xiao et al. (1996) Xiao, F., Yabe, T., Nizam, G., & Ito, T. 1996, Computer Physics Communications, 94, 103
- Yokoyama & Shibata (1994) Yokoyama, T., & Shibata, K. 1994, ApJ, 436, L197