Modeling properties of chromospheric evaporation driven by thermal conduction fronts from reconnection shocks
Magnetic reconnection in the corona results in contracting flare loops, releasing energy into plasma heating and shocks. The hydrodynamic shocks so produced drive thermal conduction fronts (TCFs) which transport energy into the chromosphere and drive upflows (evaporation) and downflows (condensation) in the cooler, denser footpoint plasma. Observations have revealed that certain properties of the transition point between evaporation and condensation (the “flow reversal point” or FRP), such as temperature and velocity-temperature derivative at the FRP, vary between different flares. These properties may provide a diagnostic tool to determine parameters of the coronal energy release mechanism and the loop atmosphere. In this study, we develop a 1-D hydrodynamical flare loop model with a simplified three-region atmosphere (chromosphere/transition region/corona), with TCFs initiated by shocks introduced in the corona. We investigate the effect of two different flare loop parameters (post-shock temperature and transition region temperature ratio) on the FRP properties. We find that both of the evaporation characteristics have scaling-law relationships to the varied flare parameters, and we report the scaling exponents for our model. This provides a means of using spectroscopic observations of the chromosphere as quantitative diagnostics of flare energy release in the corona.
The generally accepted picture of the solar flare process begins with the reconnection of magnetic field lines in the solar corona. The freshly reconnected flare loop is then free to retract under magnetic tension, which heats and compresses the loop-top plasma, forming hydrodynamic shocks (Longcope et al., 2009) and accelerating electrons near the loop apex. In the case of shocks, the steep temperature gradient between the ambient coronal plasma and the hotter post-shock plasma results in thermal conduction fronts (TCFs) that rapidly propagate down each leg of the loop (Craig & McClymont, 1976; Forbes et al., 1989; Tsuneta, 1996). In the case of accelerated electrons, the result is a large flux of non-thermal particles (NTPs) that precipitate down the loop towards the footpoints (Brown, 1973).
Although the question of which of these two models constitutes the dominant energy transport mechanism has not been resolved, the end result is similar; namely, the transport of energy down the loop which is subsequently deposited in the cooler and denser plasma in the transition region (TR) and chromosphere that lie at the loop footpoints. This deposition of energy creates a significant overpressure in the TR and upper chromosphere, and drives flows of heated plasma both up and down the loop (Fisher, 1987). These upflows and downflows are historically referred to as chromospheric evaporation (Sturrock, 1973) and condensation (Fisher, 1989), respectively, and should be distinguished from unrelated flows in the corona (such as coronal rain) that are driven by thermal instabilities due to radiative losses. Finally, the evaporation of heated dense plasma from the chromosphere fills the loop and forms the bright coronal flare loops that are visible at temperatures of several million kelvin (MK).
A critical component to understanding the subsequent flare loop development is a detailed knowledge of the characteristics of chromospheric flows during a flare. One tool for determining these characteristics are observations of Doppler spectral line shifts, which give the flow velocities within the loop at different plasma temperatures. Ideally, the Doppler line shifts would be observed with sufficient resolution (in temperature) to give a velocity profile near the point separating upflows from downflows, which would be useful in constraining the mechanism driving the flows. Unfortunately, it has proven difficult to obtain this data. Most studies instead prefer to investigate evaporation using only a small set of spectral lines, sometimes only a single one (Czaykovska et al., 1999; Wülser et al., 1994; Zarro & Canfield, 1989). Other studies concentrate instead on using observations to calculate other properties of the flare loop, such as the total quantity of evaporated plasma (Acton et al., 1982).
In the last few years, however, high-resolution spectral observations of flare footpoints have become possible thanks to the Extreme-ultraviolet Imaging Spectrometer (EIS) located onboard the Hinode spacecraft. The large number of spectral lines available with this instrument has allowed for Doppler-shifts to be derived for plasma across a broad temperature range for several different flares (Milligan et al., 2006; Milligan & Dennis, 2009; Milligan, 2011; Li & Ding, 2011). Inspection of the velocity-temperature data in these papers reveals three interesting features. First, the temperature of the point separating upflows from downflows (which we dub the “flow reversal point” or FRP) is at or above 1 MK but varies widely among observed flares (in Li & Ding (2011) it was well over 6 MK), hinting at a possible connection to properties of the flare loop. Second, the flows are broadly distributed in temperature, from K for downflows to over K for upflows. Finally, these studies used a sufficient number of spectral lines to allow a rough calculation of the velocity derivative with respect to plasma temperature, which like the flow conversion temperature varies between different flares.
Significant effort has also been devoted to modeling chromospheric evaporation using computer simulations. These simulations have generally invoked one of two candidates for energy transport: non-thermal particle (NTP) precipitation (MacNeice, 1984; Nagai & Emslie, 1984; Fisher et al., 1985a,b,c), and thermal conduction front (TCF) heating (Nagai, 1980; Cheng et al., 1984; MacNeice, 1986; Fisher, 1986). NTP models have had some success explaining flow observations. However, no model yet exists which is capable of self-consistently tracking the conversion of magnetic energy, released by reconnection, into a population of NTPs. Lacking this feature, simulations must resort to introducing the non-thermal electrons ad hoc, with a user-specified energy flux and spectrum. Properties of the evaporation flows, such as flow conversion temperature, naturally depend on this ad hoc choice.
The case of TCFs is notably different owing to the existence of comprehensive models of reconnection energy release. Large scale models of reconnection, such as the early model of Petschek (1964), have used hydrodynamic equations and thus omitted energetically significant non-thermal populations. In these models, kinetic energy is converted to thermal energy at MHD shocks, raising the post-shock loop top plasma temperature and originating TCFs due to steep temperature gradients. It is therefore possible to use these models to study how the properties of chromospheric evaporation depend on the magnetic reconnection providing the energy. As yet, however, there are no generally accepted relationships that predict evaporation velocities from flare energy. Now that observations of footpoint velocities during a flare have been made with sufficient detail to determine characteristic properties of the flows, we wish to use this characterization to infer the properties of coronal energy release.
In this paper, we use a numerical simulation code to investigate the relationship between observable properties of chromospheric evaporation during a flare and the initial properties of the flare loop. Our goal is to systematically cover a parameter space of simulation inputs in order to extract a scaling-law relationship with the output observed quantities. First, in Section 2, we develop a simplified model of a flare loop, reducing the more complex 2-D dynamics to a 1-D shocktube. In Section 3 we describe the details of our numerical simulation code, including our simplified loop atmosphere model and shock initialization. Then, in Section 4, we detail the evolution of one particular simulation including the basic hydrodynamics and the differential emission measure, develop a consistent method of extracting synthetic Doppler velocities similar to observations, and compare the results to one particular set of observed flow velocities. Finally, in Section 5 we describe the parameter survey we use to extract scaling law relationships between inputs and synthetic observations, and determine the best-fit parameters.
2 Flare Loop Model
The idealized flare loop model we use in this paper is an extension of the thin-flux-tube model developed in Longcope et al. (2009). In that model a brief, localized reconnection event is assumed to have occurred between two adjacent magnetic flux tubes previously separated by a current sheet. The sheet exists between field lines whose directions differ by less that 180 (i.e. field which is not perfectly anti-parallel). The result is a “”-shaped loop such as the one shown in the upper schematic in Figure 1 (adapted from Figure 2 in Longcope et al. (2009)) by the long dashed lines. In this picture the current sheet is located above the dashed line in the plane of the diagram, and the reconnection angle between the field lines is defined as indicated. This angle, apparent when viewing the current sheet from the side, differs from the narrow opening angle between shocks when the sheet is viewed edge-on. The latter angle has been a focus of steady-state modeling such as the seminal work of Petschek (1964); this is not the angle . Also note that a similar “V”-shaped field line would also have resulted from the reconnection, however we omit that portion in our schematic. This model also assumes that the ratio of gas pressure to magnetic pressure, known as the plasma- parameter and defined by
where is the magnetic field strength, is much less than unity. This is generally true of the pre-flare corona and transition region (TR) (Gary, 2001). The plasma- will have increased in the retracting flux tubes (i.e. outflow jets), but provided the reconnecting field was sufficiently far from anti-parallel it will still be less than unity (Longcope et al., 2009). Within the TCFs, which will be our primary concern, will lie between the initial value and that of the compressed, heated loop-top. Under the assumption of small , both the plasma and the thermal conductive flux are constrained to move only along the field line. It is also assumed in order to justify a one-dimensional treatment that the tube of reconnected flux is “thin” in the sense that the scale of variations along the loops are generally much greater than their widths. Finally, in addition to the background model developed in Longcope et al. (2009), we introduce a cool, dense chromosphere at the feet of the loop (shown as the blue portion of the tube in the upper schematic in Figure 1), which acts as a mass reservoir.
After the initial reconnection event the subsequent contraction of the field line under magnetic tension results in a shorter loop shown by the colored portion in the schematic, with the ambient coronal plasma indicated in yellow. As the loop contracts, free magnetic energy is released into accelerating the plasma downward and inward; the inward motion corresponds to motion parallel to the axis of the flux tube (Longcope et al., 2009). Starting from the initial configuration, as the contracting loop passes each of the angled solid arrows the plasma at that location is accelerated by the rotational discontinuity down and inward toward the loop center; the subsequent trajectory of the plasma is the solid arrow itself. Eventually this accelerated plasma piles up at the loop top, as shown by the red region in the schematic, resulting in heating and compression of the plasma and the formation of two slow magnetosonic shocks resembling simple gas dynamic shocks. These shocks propagate out along the loop at a hydrodynamic Mach number , determined in terms of the reconnection angle by Longcope et al. (2009) as
and they follow the trajectories given by the short dashed lines. Note that the effect of the shocks is to alter the flow of the loop plasma from downward and inward to purely downward motion (Longcope et al., 2009). At the same time, strong temperature gradients across the shock fronts, going from multi-MK post-shock plasma to 1 MK in the pre-shock coronal plasma, give rise to fast-moving thermal conduction fronts (TCFs) which move out along the loop ahead of the shocks.
Since our interest is in the effects of a shock-initiated TCF on the chromosphere, and not in the overall dynamics of the loop evolution, we narrow our focus to only that section of the flare loop indicated by the gray box in the upper schematic of Figure 1. We also adopt a reference frame that is co-moving with the contracting loop, so that the ambient coronal plasma (yellow) is stationary and the post-shock plasma (red) is being driven down the loop. We further simplify the model by neglecting gravitational stratification of the plasma. The resulting horizontal “shocktube” model of the flare loop is shown in the lower schematic in Figure 1, with the region of interest again indicated in gray. In this model, the post-shock plasma behaves as though driven by a piston, located to the right of the region of interest and moving leftward at a Mach number (referring to the pre-shock coronal plasma), and the shock front moves leftward down the tube at . Finally, we include a simplified model of the TR and chromosphere (blue in the schematic), the details of which will be discussed in Section 3.4.
3 Simulation Setup
3.1 1-D Fluid Equations
Following the above discussion we consider a one-dimensional shocktube of plasma with uniform cross-section and total length , parameterized by a coordinate , as shown in the lower schematic of Figure 1. We wish to numerically simulate the plasma hydrodynamics within the tube, beginning at an initial time forward to some later time . We begin by assuming that the plasma is everywhere of sufficient collisionality to be adequately described as a single-fluid with pressure , proton number density , average flow velocity , and temperature . In this case, we recall the 1-D hydrodynamic equations for an ideal fluid, given by
where is the proton mass, Boltzmann’s constant, is the parallel dynamic viscosity, and the thermal conductivity (discussed in Section 3.3). We adopt gas constant for a fully ionized monatomic plasma. Note that we do not include gravity in Equation (4), and hence we neglect gravitation stratification. We also do not treat explicit coronal heating or plasma radiation in Equation (5), and instead have included a single heating/cooling source term (discussed in Section 3.4) that is responsible for the equilibrium loop atmosphere. Finally, we close the system with the ideal gas law,
For this study, we define a system of dimensionless variables, where the coronal number density , temperature , sound speed , and proton mass are scaled to unity. From the equation for sound speed, , we see that the coronal pressure is rescaled to . Length is rescaled by the coronal ion mean free path, given by
after using the classical Spitzer viscosity (Spitzer & Härm, 1953), and time is rescaled to the sound transit time . Note that these new variables do not alter the form of Equations (3)-(6), except that is formally replaced by via Equation (6). Throughout the remainder of this section, we shall assume the use of the dimensionless variables.
3.2 Numerical Integration
To numerically integrate the hydrodynamic Equations (3)-(5), we first construct a staggered grid of total length and uniform cell size which defines the simulation region. The total size of the grid defined in this way is 2000 cells, to which we add two additional sets of static cells on either end to enforce the boundary conditions. These static cells are reset to their initial values after each time step. The lower boundary is completely closed (, ), and the treatment of the upper boundary will be discussed in Section 3.5. The values for the hydrodynamic variables are defined at each point on the staggered grid: bulk quantities such as and are defined at cell centers, and flux quantities such as and are defined at cell edges. We have tested our code using both 2000 and 4000 cells and found that the results do not substantially differ. We have also tested that the staggered scheme conserves mass, momentum, and energy over the simulation region, which it generally does to within 0.1% during the simulation.
With the grid and fluid variables defined, we numerically integrate Equations (3)-(5) using an explicit midpoint-stabilized stepping-algorithm for all terms except for thermal conductivity in Equation (5). Were a fully explicit scheme used the timestep size would be chosen to satisfy the Courant conditions (Courant et al., 1967),
where the minimum is taken over the full set of grid points . The first two conditions are the sound wave and flow velocity timescales, and the third is the viscous timescale. The final condition is the conductive timescale, which is in general significantly smaller than any of the other three. This is because the Prandtl number, which defines the ratio of viscosity to the thermal conductivity, is typically of order for a plasma. This results in a conductive timescale that is at least 100 times smaller than any of the other timescales, and also results in prohibitive runtimes for an explicit numerical code.
We circumvent this issue by first expanding the thermal conductive term in Equation 5 as
and then implementing an implicit Crank-Nicolson integration method (Crank & Nicolson, 1947) for the second-derivative term (the first term is folded into the normal explicit solver). This semi-implicit scheme permits us to effectively ignore the conductive Courant condition and use only the minimum of the first three terms in Equation (8). As the numerical integration proceeds, the values for , , , and for the entire grid are saved every . The entire simulation is allowed to run until the thermal conduction front, which begins at the top of the tube and propagates down, reaches the lower boundary, at which point the closed boundary condition would begin reflecting waves back up the tube. As this represents undesired (and possibly unphysical) behavior, the simulation is ended at that time.
3.3 Viscosity and Conductivity
A major obstacle to keeping the hydrodynamics well-resolved in any flare loop simulation that includes both the corona and the chromosphere is the fact that the ion mean free path given in Equation 7, which governs the length scale over which hydrodynamic quantities may vary significantly, becomes decidedly smaller as we move down from the corona into the chromosphere. In general, the chromospheric temperature is of order 100 times lower than in the corona and the density 100 times higher. Using the standard Spitzer formula for viscosity (Spitzer & Härm, 1953), and noting that for a plasma, then we see from Equation (7) that
which results in a mean free path that is six orders of magnitude smaller in the chromosphere than in the corona. For our grid spacing of this implies that there would be 50,000 mean free paths per grid cell in the chromosphere, which is inadequate to resolve fine structure hydrodynamics such as shocks.
One popular method to circumvent this issue is to use a non-uniform adaptive grid that can add or subtract grid points of varying size during the simulation to increase resolution where needed. Several established methods exist for running hydrodynamic simulations with adaptive grids, e.g. PLUTO (Mignone et al., 2007), although it is not entirely clear that such methods are able to adequately resolve shock structures in the chromosphere. Moreover, given the 1-D low plasma- nature of our hydrodynamic model, there are no additional benefits to using an adaptive grid scheme. We therefore adopt a different approach, modifying the standard Spitzer formula for viscosity by adding an additional term of the form
We see from Equation (7) that the effective mean free path then becomes
and in the corona, where we demand , , and and now are all scaled to unity, we can solve for as
To determine , we consider the mean free path in the chromosphere. In this region, due to the low temperature and high density, the first term in Equation (12) essentially vanishes, leaving
If we now impose the condition that for all points in the tube, where is a lower bound that artificially boosts the mean free path in the chromosphere, we can then solve for using Equations (13) and (14) to obtain
We find through experimentation that seems to result in adequate resolution of the hydrodynamics in the chromosphere, which thus sets and . We also performed test runs with , with the observed result that shocks became poorly resolved in the chromosphere (manifested as a slowly growing sawtooth behind the shock) and a roughly 10-20% increase in the flow reversal properties described in Section 4.3.
To determine the thermal conductivity in this model, we first recall the definition of the Prandtl number,
we adopt for the duration of this paper. Note that the modified version of in Equation (11) would result in different Prandtl numbers for the corona and the chromosphere if we used the standard Spitzer formula for conductivity . Consequently, we modify the thermal conductivity in the same manner as the viscosity, with
Substituting this expression into the expression for the Prandtl number and solving for results in .
3.4 Initial Loop Atmosphere
The left-hand portion of the simulation region contains the TR and chromosphere, as shown in Figure 1. It is well known that the structure of a static TR and chromosphere depends critically on radiation, gravity, and even ionization states (e.g. Vernazza et al. (1981), Fontenla et al. (1990), etc.). This layer responds so rapidly to the heat flux from flare reconnection, however, that these mechanisms play little role in the evaporation dynamics. The main factor determining the dynamic response is, instead, the pre-flare distribution of mass density. We therefore use a simplified physical model tuned to produce a relatively realistic initial density distribution. The chief aspect we seek to reproduce is the very large ratio of temperatures and densities, which we quantify as
Since we omit gravity the initial pressure is uniform and the density ratio is the inverse of the temperature ratio. The initial distribution is given by the expression
where is the center of the TR and is a measure of its thickness. Throughout this study we shall set (one-quarter of the way up the tube from the left-end) and . In Figure 2 we plot the temperature profile for as the solid line, and by inspection we see that our choice of results in a TR that is 10 coronal mean free paths thick.
The steep temperature gradient across the TR naturally results in a strong thermal conductive flux, given by
which transfers thermal energy from the hot corona to the much cooler chromosphere. This thermal flux (divided by 30) for the atmosphere is plotted as the dotted line in Figure 2, and clearly shows that the majority of the thermal flux occurs in the upper portion of the TR. This feature is a result of the temperature profile being defined on a log-scale, which implies that the strongest gradients in the TR will be at the higher temperature.
This term (divided by 10) is plotted as the dashed line in Figure 2, and consists of a source () in the upper layer and sink () in the lower layer. These artificial elements stand in for coronal heating and chromospheric radiation, respectively. We keep Equation (21) constant throughout the run, although we find that if we turn it off during the flare simulation there is no discernible difference in the chromospheric response. This indicates that the heating mechanism is not an essential property in determining evaporation evolution, but rather the initial mass-temperature distribution of plasma in the TR.
3.5 Initial Piston Shock
The right-hand portion of the simulation region contains the downward-propagating piston shock, which will serve to initiate and drive the thermal conduction front (TCF) into the TR and chromosphere. The classical picture of a piston shock, as shown in Figure 1, is of a plug of compressed fluid being driven at velocity down into an ambient fluid at rest with . is the Mach number of the driving piston speed as measured in the rest fluid. The shock itself is the interface between these two regions, and it propagates ahead of the piston at speed given by
The plasma compression across the shock results in an increased post-shock pressure and density, as given by the Rankine-Hugoniot conditions:
The post-shock temperature is given by Equation (6).
In the classical piston shock, the jump in velocity from to is instantaneous; the shock is a strict discontinuity in the fluid variables from pre-shock to post-shock. For a numerical simulation, however, we need to construct a smooth transition for the fluid velocity and other variables across the shock, preferably over a length scale of a few mean free paths, to ensure adequate numerical resolution of the shock (Guidoni & Longcope, 2010). In this model we initialize the shock by superimposing the post-shock plasma conditions over the upper 20% of the tube, scaled by a transition in the velocity centered at of the form
where is the initial length scale of the shock. We adopt for the remainder of this study, which results in an initial shock that is initially 10 mean free paths thick (which is identical but unrelated to the TR thickness). The post-shock region is maintained at by a flux of plasma at speed coming across the upper boundary (). This flux is the result of the boundary condition we enforce for the static cells on the right-hand side of the grid, which are reset at each timestep to the original post-shock conditions. We have tested our code by simulating shocks in the absence of the TR and thermal conduction, and found that shocks do indeed propagate at the correct speed while remaining well-resolved due to the presence of viscosity in the shock region.
4 Simulation Results
To discuss the various features of the simulated loop dynamics, we focus first on a single simulation. The qualitative results of the evolution of this particular simulation are similar for most of the runs performed in this study. We consider a TR temperature ratio and piston Mach number , and assume an ambient coronal temperature of MK and number density cm. Restoring conventional dimensions to variables (assumed throughout this section) results in a total length for the simulation region Mm, a coronal sound speed km s, and a post-shock temperature MK. The total duration of the simulation (from the initial state to the TCF reaching the lower boundary) is seconds.
The hydrodynamic evolution of the simulation region is shown in Figure 3; seven different times are plotted and color-coded according to the legend for pressure (upper left plot), velocity (upper right plot), number density (lower left plot), and temperature (lower right plot). At the initial simulation time, sec (black line), as we move up the tube from Mm we note first the cool, dense chromosphere at K and density cm. Beginning at 7 Mm, we encounter the artificial TR, which appears only in density and temperature and continues to 11 Mm. Above the TR, we have the constant temperature and density corona, which occupies more than 50% of the length of the tube. Finally, centered between 31-34 Mm, we note the initial piston shock which accelerates the post-shock plasma to km s (negative as the shock is propagating down the tube), and which heats and compresses the plasma to 9.2 MK and cm.
By sec (purple), we see the very rapid development of the TCF in the plasma temperature. Within this time, the TCF has propagated down to 20 Mm and has closed nearly half the distance between the initial shock position and the TR. We see from the density and velocity profiles in Figure 3 that the shock itself has not moved downward very far (indeed, the density profile for the shock is scarcely different than at the initial time). This is made clearer in the pressure profile where we note that the plasma pressure rises once between 20-31 Mm due to the presence of the TCF, and subsequently rises again across the shock between 31-35 Mm. This decomposition of a shock, in the presence of thermal conduction, into a TCF and a so-called isothermal sub-shock is common and has been observed in previous models (Guidoni & Longcope, 2010; Longcope et al., 2009).
After sec of evolution (blue), the TCF has propagated downward far enough that it encounters the cooler and denser TR plasma, which results in two distinct effects. First, the TCF slows dramatically due to the reduced thermal conductivity in the TR and chromosphere; indeed, the TCF takes 1 second to descend the 20 Mm between the initial piston shock and the TR, but takes another 37 seconds to clear the TR and chromosphere and reach the lower tube end. Second, the TCF begins rapidly depositing thermal energy into the stationary TR plasma. This rapid rise in plasma temperature in the upper TR, coupled with the stationary density profile, results in the development of a large overpressure clearly visible centered at 9 Mm. This TR overpressure is responsible for the initiation of the chromospheric evaporation upflows (visible between 9-11 Mm) and associated condensation downflows (visible between 8.5-9 Mm).
As the simulation continues to sec (cyan), we see that the TCF has completely cleared the TR and has begun to directly heat the chromosphere. Meanwhile, the evaporation and condensation has continued to develop and we see clearly that the upflows have higher speeds and a broader spatial distribution than the downflows. This is a reflection of the momentum balance in the TR: the TCF deposits thermal energy to the TR but no net momentum, and as evident in the density profile the condensation is occurring in a denser region than the evaporation. Also in the density profile at this time, we note the development of an evaporation front, located at 12 Mm, which is beginning to enhance the density of the upper TR and lower coronal regions, and a barely visible rarefaction region and condensation front in the lower TR.
We now jump forward to sec (green). By this time, the evaporation region has grown to encompass nearly one-third of the tube, and has developed upflow speeds of 500 km s. Meanwhile, the condensation region is restricted to speeds less than 100 km s. In the density profile we now see the fully developed three-part structure of condensation (enhanced densities between 4.5-6.5 Mm), rarefaction (decreased densities between 6.5-9 Mm), and the evaporation front (which has strongly enhanced densities up to 16 Mm). There has also begun to be significant interaction between the upward propagating evaporation front and the downward propagating subshock (centered at 23 Mm), with mildly enhanced densities in between. Finally, at this point in the simulation, we begin to observe the direct effects of the artificially-enhanced thermal conductivity , which manifests as a “shoulder” in the TCF located between 5-6 Mm.
By sec (orange), the upflows have reached and passed the maximum speed during this simulation of 510 km s. For later times the maximum upflow speed in the tube is below this. This is the result of the downward propagating subshock finally encountering the evaporation front and passing through it, which is especially evident in the pressure and density profiles at 19 Mm. The strongly negative post-shock velocities thus begin to cancel the positive evaporation velocities, although the enhanced pressure that results from the combined compression of the shock and evaporation front means that some positive velocities will remain in the post-shock region.
At the end of this particular simulation, sec (red), the TCF has fully passed through the chromosphere and has developed a distinct two-step profile due to the enhanced thermal conductivity. However, the TCF “shoulder” remains somewhat below the lower-bound of the evaporation region, and thus is not likely influencing the development of the evaporating plasma. Meanwhile, the piston subshock and evaporation front have fully passed each other, resulting in a region of highly enhanced density (1.6 cm) between 18-24 Mm. The maximum upflow speed has been reduced to 440 km s, and a uniform upflow speed of 130 km s has developed in the region between 18-24 Mm.
Finally, to conclude our discussion of the hydrodynamic evolution of the simulation we consider the ratio of the thermal conductive flux, Equation (20), to the free-streaming saturation limit, given in non-dimensional form by Longcope & Bradshaw (2010) as
where is the ratio of proton to electron masses. The ratio is plotted in Figure 4 for the same times and with the same color scheme as in Figure 3. We observe that at sec there are two peaks in the flux ratio: one for the TR centered at Mm, and another representing the shock centered at Mm. We also note that the ratio is greater than unity for a narrow range of positions centered on the initial piston shock, indicating that the thermal flux across the shock is larger than the saturation value. This might indicate a substantial problem were the thermal flux to remain supersaturated for the duration of the simulation. However, by sec, we see that the flux ratio has been reduced to 0.4 everywhere in the tube, due to the development of the TCF discussed above which quickly smoothes out the initial steep temperature gradient. This fast TCF and thermal flux development is characteristic of all simulations performed for this study, and we do not believe that this initial violation of the free-streaming saturation limit by the piston shock is of concern for the later tube evolution.
Later, as the TCF reaches the TR (1.02 sec), we note an enhancement of the flux ratio at 11 Mm as the TCF begins depositing thermal energy into the cooler TR plasma. This peak slowly begins to subside as the TCF clears the TR (5.45 sec and on), although it is never fully eliminated, remaining as a small “bump” at 11 Mm. Curiously, we also observe that the later evolution of the flux ratio (, 23.1, & 38.0 sec) somewhat mirrors that of the density and pressure. This is especially notable for positions between 15-25 Mm, where we notice the density enhancement due to the interaction of the evaporation and subshock fronts mirrored as a suppression of the flux ratio. This behavior is not indicative of any change in the thermal flux, as the TCF has long since flattened the temperature profile to nearly isothermal. Rather, it is due to the increased density resulting in a larger saturation limit, lowering the flux ratio at those positions.
4.2 Differential Emission Measure
Although the full hydrodynamic evolution (as shown in Figure 3 and described above, for example) would be the preferred method to understand the plasma dynamics in a flare loop, we are limited by observational techniques in our ability to extract information about those quantities. One observational method of tracking the plasma evolution, which combines information about the density and temperature, is the differential emission measure (DEM), defined as
where is the electron number density (identical to the proton number density in our fully-ionized hydrogen plasma). In Figure 5 we show the evolution of the DEM as a function of temperature in the tube, for the same times and with the same color scheme as in Figures 3 & 4.
At the initial time, sec, the clearest features of the DEM are the three sharp spikes located at K, 2.5 MK, and 9.2 MK. These peaks correspond to the uniform temperature chromosphere, corona, and post-shock regions seen in Figure 3. Also notable is the DEM minimum located at 2 MK which is a somewhat higher temperature than seen in other observational and modeled DEMs, although the overall magnitude of our DEM is comparable (Emslie & Nagai, 1985; Brosius et al., 1996). We attribute this higher-temperature minimum to the fact that our model atmosphere, Equation (19), has its steepest gradient at higher (2 MK) temperatures, thus resulting in the DEM being minimized at those temperatures.
By sec, the TR portion of the DEM between K and 2.5 MK remains unchanged. Only the portion of the DEM between the uniform corona and the post-shock region has been altered as the TCF begins to smooth the temperature gradient across the piston shock, resulting in some enhancement of the DEM in the 3-8 MK range. For sec, however, we observe significant changes to the DEM in the corona and upper TR; indeed, the entire range from K to 2+ MK has been enhanced by a factor of 10 to 100. At this same time, recall from Section 4.1 that upflows and downflows in the TR are beginning to form. To differentiate between the portion of the DEM that concerns upflows from the portion that concerns downflows, we have plotted several vertical dashed lines, in the same color palette, that indicate the temperature where (i.e. the flow reversal point (FRP) separating the condensation and evaporation regions). This will be referred to henceforth as the FRP temperature , which is defined in terms of the simulation as the temperature at the first position where .
For the sec profile, the evaporation temperatures range from the dashed blue line at 2.0 MK to slightly less than 4 MK. Similarly, for the sec and sec profiles, the evaporation begins at and ranges up to 5 MK and 7 MK respectively. In these two profiles, we note the evaporation front closing with the high-temperature post-TCF subshock, represented here by the peak in the DEM at K. Note that the post-shock DEM enhancement is roughly 10-fold, and since this corresponds as expected to the three-fold density increase across the subshock (see Figure 3). We also note, in the sec profile, the DEM enhancement of the condensation front at K, which is partly due to the enhanced post-condensation density and partly to the “shoulder” which forms on the TCF as described in Section 4.1.
After the evaporation front and subshock interact ( sec and sec profiles), we observe another roughly 10-fold DEM increase in the 7-9 MK range from the combined compression of the plasma. Further, the continued compression of the post-condensation front plasma has continued to enhance the DEM, forming a large peak between K and K. Finally, to conclude our description of the DEM evolution for this simulation, we observe that the flow conversion temperature first appears at a somewhat higher temperature of 2.0 MK, subsequently descends to a lower range of 1.3 to 1.4 MK, and later rises again by the end of the simulation to 1.9 MK. Although the exact temperatures vary, this decreasing-and-increasing behavior is typical of the simulations used in this study.
4.3 Synthetic Doppler Velocities
As discussed in Section 1, observations of flare loops using the Hinode/EIS instrument have revealed temperature-dependent Doppler velocity profiles for plasma at the loop footpoints during chromospheric evaporation. We would thus like to construct a similar velocity-temperature profile for the simulation results, in order to compare to these observations. However, we cannot simply use the velocity and temperature profiles as shown in Figure 3, as this is not really what is being observed by Hinode/EIS. Instead, the temperature-dependent Doppler velocities are derived from spectral lines from an exposure recorded over approximately 5-10 seconds, and which are weighted by the amount of emission coming from the plasma at that temperature.
We thus wish to construct a “synthetic” Doppler velocity to compare with data. We begin by defining the plasma emission measure
and an emission-weighted plasma velocity
We next define a binned log-temperature scale with 100 equal bins per unit interval in , and create binned versions of the emission measure, , and EM-weighted velocity, , by summing and over each temperature bin. If there are no grid points in a given bin at that time, the values of and for that bin are set to zero. Finally, we define a time-window and construct the synthetic Doppler velocity as
The variable window size allows us to consistently accommodate different simulation durations. For the simulation discussed thus far sec, which is a typical exposure duration for Hinode/EIS. Note of course that is only defined for , which for the Section 4.1 simulation corresponds to times between 4.75 sec and 33.25 sec.
In Figure 6 we have plotted the hydrodynamic velocity as a function of temperature (solid lines) and the synthetic Doppler velocity as a function of temperature (dashed lines), for five representative times during the simulation and for temperatures above K. Note that neither the times nor the color palette are identical to those plotted in Figures 3-5; this is because is not defined for the three earliest times or the final time in those plots. However, we are able to include sec (cyan) which is close to the start of the windowing, sec (green), and sec (orange). We have also added two additional times: sec (replacing purple), and sec (replacing red) which is close to the end of the windowing.
We note that one effect of the synthetic Doppler processing is to shift the velocity profiles in the upflow region to higher temperatures, particularly notable for the three earliest times. The peak Doppler upflow speed is increased by 50 km s over the hydrodynamic velocity for the 5.45 sec profile, but is generally the same or slightly reduced for later times. Downflow speeds are slightly increased for temperatures below 1 MK, again more significantly for early times, but the observed downflow speeds remain significantly less than the upflow speeds. The FRP temperature separating upflows and downflows appears mostly unaffected by the processing, however the way in which we define the FRP needs to be modified due to the temperature binning. Recall that we previously defined as the temperature at the first position where ; we now determine the temperature bin where for each time , and perform a fit to the two bins of the form
We now define the FRP temperature as
and we note that is the slope of the velocity-temperature profile at the FRP.
In the lower plot of Figure 6 we have indicated and for the sec profile as the vertical dashed line and the dash-dotted line respectively. As an inspection of this plot will indicate, however, these two quantities do vary over the course of the simulation. To track the evolution of and we have plotted them as functions of time in the upper and middle plots in Figure 7, respectively, and as functions of each other in the lower plot in Figure 7. The notable “jitteriness” of the profile is due to movement between temperature bins when tracking . We have also plotted the values for and at the five times shown in Figure 6 as solid squares in the upper two plots. With some exceptions (discussed in Section 5) the behavior of the FRP properties for other simulations is similar to that seen in Figure 7.
Finally, we define a mean FRP temperature and slope , calculated by taking a time-average of and over the range . For this simulation MK and km s, and we have plotted these values in Figure 7 as the horizontal dashed lines in the top and middle plots for and , respectively, and as an in the lower plot. Obviously, the time-dependent values for the FRP properties differ from these mean values over the course of the simulation; in fact, the lower plot in Figure 7 shows that they do not ever assume the mean values simultaneously. However, by using the mean values we can in some sense characterize the entire evolution of and for a given simulation, which will allow us in Section 5 to directly compare these values for many different simulations.
4.4 Observational Data Fit
We conclude our discussion of this particular simulation by making a comparison of our synthetic Doppler velocities to the observed flare loop Doppler velocities published by Milligan & Dennis (2009) (with a correction published in Milligan (2011)). The event studied in that paper was a GOES-class C1.1 flare that took place in NOAA AR 10978 on 2007 December 14 at 14:12 UT. Serendipitously, the Hinode/EIS instrument was rastering over one of the flare loop footpoints during the impulsive phase of the flare, with an exposure time of 10 seconds. The authors used 15 different spectral lines, with formation temperatures ranging from K to 16 MK, to derive Doppler velocities for the footpoint plasma. We have taken these Doppler velocity data and associated error ranges from Table 1 in Milligan (2011), and replotted them in Figure 8 as the square points and error bars (note that we have changed the signs on these data to match our velocity convention).
Inspection of the observational data in Figure 8 reveals similar structure to the synthetic Doppler velocity profiles in Figure 6, with downflows and upflows separated by a flow reversal temperature that lies somewhere between 1-2 MK. To estimate the FRP parameters for comparison to the simulation, we take the six velocity measurements that fall in the 1-2 MK range, and use a linear fit of the form
where is the approximate FRP slope; the approximate FRP temperature is given by . For the Milligan (2011) data we find that MK and km s, and these values are represented by the “” in Figure 7 as well as the vertical dashed line and dash-triple-dotted line, respectively, in Figure 8.
Since the observed flow conversion temperature is the quantity in which we have the most confidence, we fit the simulation to the data by selecting the time for which provides the best match to . This is found to be at sec (orange line in Figures 3-6) with MK and km s, plotted as the solid square in Figure 7. That profile for the synthetic Doppler velocity has been plotted on top of the observation data in Figure 8, along with a dash-dotted line representing the slope . Aside from the approximate match to the FRP properties, we note that the profile from the simulation matches well to the observed velocities in the range K K. Outside of this range, however, the simulation profile begins to diverge from the observed values, especially for the highest temperatures (12.5 MK and 16 MK). Indeed, these temperatures do not even exist in this simulation, which has a maximum post-shock temperature of 9.2 MK.
We also observe from Figure 7 that the mean slope is a somewhat better estimate for than the value for at secs, but that is a slight underestimate for . Nevertheless, as an order-of-magnitude estimation, the mean value is not far removed from the observational result. Thus, we feel we are justified in using the mean FRP properties and as proxies for describing the overall evolution of the evaporation during a flare.
5 Scaling Laws
Thus far we have developed a method for reducing the complex properties of the spatially- and temporally-dependent chromospheric evaporation driven by thermal conduction in a simple model atmosphere down to two scalar quantities, namely and . Further, we have shown that these quantities provide an acceptable description of the FRP properties in observed chromospheric flows. An obvious question now presents itself: is it possible to systematically relate the observed FRP properties and back to fundamental parameters in the simulation?
Of course, with only two data inputs, it will only be possible to extract at most two parameters for the simulation. The two most useful properties are likely the TR temperature ratio , which gives some insight into the pre-flare state of the loop, and the Mach number of the piston shock , which has been shown previously to relate to the initial reconnection angle of the loop (Longcope et al., 2009). We therefore seek an invertible relationship between and ; as we shall show, such a relationship does exist, but instead of we shall use the post-shock temperature . However, the post-shock temperature is dictated uniquely by the piston Mach number, as described in Section 3.5, and thus the two are effectively equivalent.
To determine the desired parameter relationships, we require a set of simulations that adequately span the parameter space. At this point, in order to avoid confusion, we shall begin labeling dimensionless quantities explicitly with a superscript “”; variables with ordinary units will be left as normal. We choose five values for the piston Mach number , given by
these values translate to the equivalent (dimensionless) post-shock temperatures,
We also choose five values for the TR temperature ratio , given by
Simulations were performed for all 25 combinations of the two parameters. For brevity, we will label and discuss these simulations based on the parameter values selected: the five values for are labeled “A”-“E”, and the five values for are labeled “1”-“5”. Thus, the simulation discussed at length in Section 4, which had and would be labeled as “A4”.
We next calculate, as in Section 4.3, the dimensionless mean synthetic Doppler velocity FRP temperature and slope for each simulation. The values of these two quantities for all 25 simulations are tabulated in Table 1. In Figure 9 we show the results of sequences “2” (dashed), “5” (broken), “A” (dotted) and “D” (dash-dot). This parameter survey shows that increasing the chromospheric density ratio leads to decreases in both FRP temperature (9c) and slope (9a). Increasing the post-shock temperature leads to an increase in each of these (9d and 9b). The flow conversion temperature lies naturally beneath the post-shock temperature solid curve in Figure 9d but can be below the coronal temperature if the chromospheric ratio is sufficiently large ( in Figures 9c and 9d). All 25 runs are fit to a pair of power law relations of the form
where , , , , , and . These fits are plotted as curves in Figure 9.
For the foregoing power-law fits we have chosen to omit the values for the “D1”, “E1”, and “E2” simulations, as doing so results in a much stronger fit for the remaining 22 simulations (the percent error between the actual values and the fit values are also tabulated in Table 1). As shown in Figure 9, the three omitted simulations have values for that fall between 10-19% below the power-law fit, and values for that fall between 14-38% below the power-law fit. We believe the reason for the poor fit for these simulations is that they fall in a “weak-TR/strong-shock” regime, where the TCF simply moves through the TR and chromosphere too quickly, and thus the time-averaged FRP properties do not reflect the later evolution seen in Figure 7.
Since the dimensionless versions of the scaling laws (37) & (38) are not especially useful for handling observational data, we rescale them by using the definitions and , and the definitions of and given above. We also use the fact that the sound speed in the corona is given by , where km s K. After some rearrangement, and making the assumption that and are adequate proxies for and , we obtain the following:
Note that these versions of the scaling relationships require the assumption of one of the parameters , , or . Recall that throughout Section 4 we restored the simulation variables to conventional units in part by setting K. Since for that simulation, it follows that K. We now generalize this assumption for all our simulations by fixing K, and we then invert Equations (39) & (40) to yield the flare parameters explicitly from observables:
where , , , , , , and , and is in units of km s.
As a final check we determine if the run presented in Section 4, which fit both the observed Doppler velocity data and flow conversion properties quite well, is actually the run that would be suggested by the above scaling laws. We adopt K and use the observed values for MK and km s in Equations (41) & (42) to obtain a suggested coronal temperature MK and post-shock temperature MK, implying a dimensionless TR ratio and post-shock temperature . From Table 1, we see that the simulation closest to these suggested values is indeed “A4”, with and , indicating that the derived scaling laws do indeed yield reasonable results.
In this paper we have developed a numerical simulation code to investigate relationships between certain observable properties of chromospheric evaporation and more fundamental (but difficult to determine) properties of the loop atmosphere and of coronal energy release. This code is an extension of the model developed in Longcope et al. (2009), in which the coronal plasma in a post-reconnection flux tube is accelerated and compressed due to the contraction of the loop under magnetic tension. The compressed plasma at the loop top results in slow magnetosonic shocks, and the post-shock plasma is heated to flare temperatures of several megakelvin from the conversion of free magnetic energy to kinetic and then to thermal energy. We have extended this model to include a highly simplified model chromosphere and transition region (TR) at the loop footpoints to act as a mass reservoir. Thermal conduction from the post-shock plasma transports heat down toward the footpoints, resulting in impulsive heating of the chromospheric and TR plasma and a disruption of the local thermodynamic equilibrium and finally in the bulk flows of plasma known as chromospheric evaporation and condensation.
In creating the numerical simulation code, we made a variety of simplifying assumptions beyond the piston shock model. First, we opted to ignore gravitational stratification and loop geometry, effectively assuming a horizontal tube of uniform cross-section. We also left out explicit radiation and coronal heating, instead setting up our model loop atmosphere using a simple function for the temperature profile. We then calculated the necessary heating input to maintain that atmosphere at equilibrium and supplied that to the loop for the entire simulation. We chose these particular simplifications in part because they allowed us to assign definite values for certain parameters of interest (e.g. the ratio of the coronal to chromospheric temperatures), and thus more easily make comparisons between simulations with different values for those parameters. We did test the effect of the heating input on the simulations by turning it off for a test run, and we found no significant impact on the result for that simulation. Although we did not present it, we also conducted a preliminary test varying the thickness of the transition region, and found only a very weak impact on the results. However, since we used the same function for the temperature profile for all simulations, we are unable to make any claims about how a different temperature and density profile might alter the results. Finally, recall that in order to adequately resolve features in the chromosphere, we artificially enhanced the viscosity and conductively at low temperatures. This has the effect of altering short-scale physics, and also reducing velocities at low-temperatures, and decreasing this enhancement by 50% results in 20% and 10% increases in and , respectively. We do not believe, however, that this substantially alters either our data comparison or our final conclusions.
To demonstrate the results of the code, we discussed in detail the hydrodynamic and DEM evolution of a single simulation, as the results for our other simulations were qualitatively similar and differed only in quantitative details. In terms of broad-scale features, we observe that our results using this code are also qualitatively similar to simulations performed in other studies of impulsive conduction-driven evaporation. This includes the rapid development of the thermal conduction front (TCF) and development of the isothermal subshock (ISS), the overpressure in the TR, and the presence of both upflows and downflows (so-called “explosive” evaporation) with upflow velocities dominating (Nagai, 1980; Cheng et al., 1984; MacNeice, 1986; Fisher, 1986). The magnitude and shape of the DEM profile was found to be similar to other studies (Emslie & Nagai, 1985; Brosius et al., 1996). This is encouraging, since our model is in many ways simpler than other models (as discussed above) but still seems to evolve in a broadly similar manner. Our model differs from most others, however, in one particularly substantial aspect: the use of a hydrodynamic shock to drive the TCF. This shock, and the associated downflows, introduce a complex interplay between the evaporation front and the isothermal subshock. It also serves to resolve the issue of the free-streaming saturation limit, which is easily violated using ad hoc heating models. Instead we find that, after a brief initial violation in the piston shock, the decomposition into the TCF and ISS quickly restores the thermal flux to below the saturation limit.
Our goal in this study was not to investigate the hydrodynamic and DEM evolution of chromospheric evaporation, which, as we have noted above, have been extensively studied elsewhere and in greater depth. Instead, our purpose was to draw a connection to footpoint velocities obtained from spectral Doppler shifts in data from instruments such as Hinode/EIS. To do this, we constructed “synthetic” Doppler velocities from our simulation results by first weighting the plasma flow velocities by the emission measure, then binning over temperature, and finally averaging over a time-window of similar duration to typical Hinode/EIS exposure times. Although this method is not as exact as, for example, using the CHIANTI database to construct the actual spectral lines and associated Doppler shifts, we find that the results are more flexible for investigating temperature-dependent properties of observable footpoints flows.
The properties of particular note for our purpose center on the existence of a flow reversal point (FRP) near the loop footpoint that separates upflows from downflows. As seen in observational data, the FRP occurs at an identifiable temperature and with an identifiable velocity-temperature slope, and these properties can be seen to vary over time and among flares (Milligan & Dennis, 2009; Milligan, 2011; Li & Ding, 2011; Raftery et al., 2009). We have identified and tracked these two properties for our synthetic Doppler data, and we find that they do indeed vary over time, but that they are also confined to a fairly narrow range of values for each simulation. We thus calculated the mean values for the FRP properties during each simulation to serve as a convenient proxy while comparing different simulations. Finally, we compared the velocity data from Milligan (2011) to both the “best-fit” and the mean FRP for one of the simulations and determined that both may be considered reasonable given the data. This fact, along with the observation that the FRP properties evolve similarly across a range of simulation inputs, indicates that the mean FRP properties are a robust proxy for the overall FRP evolution during impulsive evaporation.
Having chosen the simulation inputs (namely the Mach number of the piston and the ratio of coronal-to-chromospheric temperatures ) and outputs (the mean FRP temperature and slope), we then investigated whether a simple relationship may be determined using this simulation code. We have shown that it is possible to extract a general scaling-law relationship between these quantities, after replacing the piston shock Mach number with the (formally identical) post-shock temperature . A few exceptions to these relationships exist, particularly among the simulations where the ratio is low and the Mach number is large (the “strong-shock/weak-TR” limit), but we note that the relationships hold well if these cases are ignored. It is possible that these cases, which we note have the shortest duration among the 25 simulations, simply do not have enough time to evolve as fully as the others, particularly the FRP temperature which tends to rise as the simulations evolve. One potential solution would be to extend the length of the simulation region, particularly the chromospheric depth, to allow for runs that are of equal duration; however, this has not been tested. Finally, we observe that the scaling exponent trends make some physical sense: a larger value for (i.e. a stronger driving shock) results in higher temperatures and flow velocities and hence in larger values for the mean FRP properties, whereas a larger value for (i.e. a hotter pre-flare corona) results in lower temperatures and slower flows and hence suppresses the mean FRP properties. However at this point we have no intuition regarding the exact values for the scaling exponents, and further work will be needed to determine if these values may be derived analytically.
When applied to the data from Milligan (2011), the scaling-law relationships we have determined suggest input parameters that a reasonably close to the “best-fit” that we selected from among our performed simulations. This again strongly suggests that the mean FRP properties are an adequate proxy for the overall evolution of the evaporation dynamics. Further, the ambient coronal temperature and post-shock temperature are both reasonable and typical for active regions and flare loops. That being said, we have also performed preliminary tests of our scaling law relationships on other observed flare footpoint flow profiles (e.g. Li & Ding (2011), Raftery et al. (2009)) and we find these results somewhat less encouraging. For example, one set of Doppler velocities from Li & Ding (2011) yield an ambient coronal temperature of 24 MK and post-shock temperature of 45 MK, and data from Raftery et al. (2009) suggests MK and MK. In both cases the ambient coronal temperature is significantly higher than is typically expected, and the post-shock temperature suggests a very weak effective piston shock (). Data from both papers suggest FRP temperatures comparable to what we have presented here, but have FRP slopes that are much shallower than is seen in any of our simulations.
It therefore seems likely that there is some effect, besides viscosity, acting to suppress the plasma velocities in at least the evaporation region, and possibly the condensation region as well. We posit, without proof, that these unreasonable values for the coronal and post-shock temperature may be resolved by incorporating a flux tube “nozzle” at or near the TR. Some constriction of the magnetic flux tube is expected near the loop footpoint due to a combination of canopy expansion and magnetic pressure. The authors believe that the presence of a nozzle near the TR will modify the properties of the resulting chromospheric flows, possibly suppressing flow velocities in such a way that the shallow FRP slope is reproduced for more reasonable parameters of the flare loop. This topic will be the subject of a future study using an extension of this numerical simulation model.
|Label||Fit error||Fit error||Fit?|
|A1||100||2.0||3.67||0.929||1.8 %||2.06||-1.0 %||Yes|
|A2||150||2.0||3.67||0.764||0.36 %||1.44||-1.8 %||Yes|
|A3||200||2.0||3.67||0.659||-1.5 %||1.14||-0.38 %||Yes|
|A4||250||2.0||3.67||0.584||-3.5 %||0.960||2.1 %||Yes|
|A5||300||2.0||3.67||0.526||-5.7 %||0.841||4.8 %||Yes|
|B1||100||2.5||4.93||1.60||1.8 %||4.13||2.7 %||Yes|
|B2||150||2.5||4.93||1.35||2.7 %||2.81||-0.63 %||Yes|
|B3||200||2.5||4.93||1.18||2.1 %||2.15||-2.5 %||Yes|
|B4||250||2.5||4.93||1.05||1.1 %||1.76||-2.9 %||Yes|
|B5||300||2.5||4.93||0.959||-0.19 %||1.51||-2.9 %||Yes|
|C1||100||3.0||6.47||2.51||-3.2 %||7.42||0.49 %||Yes|
|C2||150||3.0||6.47||2.20||1.5 %||5.34||2.8 %||Yes|
|C3||200||3.0||6.47||1.96||2.9 %||4.08||0.74 %||Yes|
|C4||250||3.0||6.47||1.77||3.0 %||3.29||-1.5 %||Yes|
|C5||300||3.0||6.47||1.63||2.7 %||2.77||-2.8 %||Yes|
|D1||100||3.5||8.28||3.59||-12 %||10.5||-18 %||No|
|D2||150||3.5||8.28||3.29||-3.7 %||8.93||-0.96 %||Yes|
|D3||200||3.5||8.28||3.01||0.38 %||7.23||3.0 %||Yes|
|D4||250||3.5||8.28||2.78||2.4 %||5.94||2.6 %||Yes|
|D5||300||3.5||8.28||2.58||3.3 %||5.00||1.1 %||Yes|
|E1||100||4.0||10.4||4.78||-23 %||11.3||-47 %||No|
|E2||150||4.0||10.4||4.53||-13 %||11.9||-21 %||No|
|E3||200||4.0||10.4||4.27||-6.5 %||11.0||-6.2 %||Yes|
|E4||250||4.0||10.4||4.03||-2.4 %||9.63||-0.07 %||Yes|
|E5||300||4.0||10.4||3.81||0.13 %||8.41||2.2 %||Yes|
- slugcomment: To appear in Ap. J.
- Acton, L.W., Canfield, R.C., Gunkler, T.A., et al. 1982, ApJ 263, 409.
- Brosius, J.W., Davila, J.M., Thomas, R.J., et al. 1996, ApJS 106, 143.
- Brown, J.C. 1973, SoPh 31, 143.
- Cheng, C.C., Karpen, J.T., & Doschek, G.A. 1984, ApJ 286, 787.
- Courant, R., Friedrichs, K., & Lewy, H. 1967, IBM J. 11, 215.
- Craig, I.J.D. & McClymont, A.N. 1976, SoPh 50, 133.
- Crank, J. & Nicolson, P. 1947, Proc. Camb. Phil. Soc. 43 (1): 50-67.
- Czaykovska, A., DePontieu, B., Alexander, D., et al. 1999, ApJ 521, L75.
- Emslie, A.G. & Nagai, F. 1985, ApJ 288, 779.
- Fisher, G.H., Canfield, R.C., & McClymont, S. 1985(abc), ApJ 289, 414.
- Fisher, G.H. 1986, Lec. Notes in Phys. 255, 53.
- Fisher, G.H. 1987, ApJ 317, 502.
- Fisher, G.H. 1989, ApJ 346, 1019.
- Fontenla, J.M., Avrett, E.H., & Loeser, R. 1990, ApJ 355, 700.
- Forbes, T.G., Malherbe, J.M., Priest, E.R. 1989, SoPh 120, 285.
- Gary, G.A. 2001, SoPh 203, 71.
- Guidoni, S.E. & Longcope, D.W. 2010 ApJ 718, 1476.
- Li, Y. & Ding, M.D. 2011, ApJ 727, 98.
- Longcope, D.W., Guidoni, S.E., & Linton, M.G. 2009 ApJ 690, L18.
- Longcope, D.W. & Bradshaw, S.J. 2010 ApJ 718, 1491.
- MacNeice, P., McWhirter, R.W.P., Spicer, D.S, et al. 1984 SoPh 90, 357.
- MacNeice, P. 1986, SoPh 103, 47.
- Milligan, R.O., Gallagher, P.T., Mathioudakis, D., et al. 2006, ApJ 638, L117.
- Milligan, R.O. & Dennis, B.R. 2009, ApJ 699, 968.
- Milligan, R.O. 2011, ApJ 740, 70.
- Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228.
- Nagai, F. 1980, SoPh 68, 351.
- Nagai, F. & Emslie, A.G. 1984, ApJ 279, 896.
- Petschek, H.E. 1964, Proc-1964-Hess, 425.
- Raftery, C.L., Gallagher, P.T., Milligan, R.O., et al. 2009, A&A 494, 1127.
- Spitzer, L.Jr. & Härm, R. 1953, Phys. Rev. 89(5), 977.
- Sturrock, P.A. 1973, Proc-1973-Ramaty, 3.
- Tsuneta, S. 1996, ApJ 456, 840.
- Vernazza, J.E., Avrett, E.H., & Loeser, R. 1981, ApJS 45, 635.
- Wülser, J.P., Canfield, R.C., Acton, L.W., et al. 1994, ApJ 424, 459.
- Zarro, D.M. & Canfield, R.C. 1989, ApJ 338, L33.