A SIMPLE MODEL OF CHROMOSPHERIC EVAPORATION AND CONDENSATION DRIVEN CONDUCTIVELY IN A SOLAR FLARE

D. W. Longcope

doi:10.1088/0004-637X/795/1/10

1. INTRODUCTION

Solar flares are events in which large amounts of magnetic energy, stored in the coronal field, are rapidly converted to other forms. Many of the dramatic consequences of flares result less directly from the energy release than from the large mass of chromospheric material that is heated and driven upward in a process known as chromospheric evaporation (Canfield et al. 1980; Antonucci et al. 1999). Direct signatures of this process are observed in Doppler shifts of hot spectral lines (Antonucci & Dennis 1983; Zarro & Lemen 1988; Brosius & Phillips 2004; Milligan & Dennis 2009), but the tremendous increase in the emission measure of high-temperature plasma provides equally compelling, albeit indirect, evidence (Neupert 1968).

Since it was first proposed chromospheric evaporation has been investigated by numerical simulations solving one-dimensional gas-dynamic equations (Nagai 1980; Somov et al. 1981; Peres et al. 1982; McClymont & Canfield 1983). Some investigations introduce the flare energy as a beam of non-thermal electrons impacting the chromosphere (MacNeice et al. 1984; Fisher et al. 1985a, 1985b, 1985c), while others use an ad hoc heat source situated at the loop top, whose energy is carried to the chromosphere by thermal conduction (Cheng et al. 1983; MacNeice 1986). In both scenarios, beam-heated and conductive, chromospheric material is heated and driven upward at speeds comparable to, or exceeding, the local sound speed. In many simulations there is a downward flow, called chromospheric condensation, for which some spectroscopic observations provide separate evidence (Ichimoto & Kurokawa 1984; Brosius & Phillips 2004; Brosius 2009; Milligan & Dennis 2009).

It is generally believed that magnetic reconnection is ultimately responsible for releasing the coronal magnetic energy and thereby initiating a solar flare. The ad hoc energization invoked in one-dimensional gas-dynamic simulations mentioned above, either thermal or non-thermal, is thus intended to model reconnection. Emission of hard X-rays or microwaves from loop footpoints provide evidence of non-thermal electrons energizing the chromosphere. There are, however, flares in which substantial evaporation occurs in the absence of footpoint emission (Zarro & Lemen 1988; Longcope et al. 2010). In such cases the downward transport of reconnection energy can be attributed to thermal conduction. Moreover, even in those flares where it does occur, footpoint emission is observed to abate well before the flare's energy release ends. Evaporation in the later stages of these flares is presumed to be driven by thermal conduction from the reconnection site at the loop tops. It has been proposed that the increased density resulting from chromospheric evaporation comes to inhibit the propagation of non-thermal beams, thereby causing the transition to conduction-driven evaporation (Liu et al. 2006).

Currently, we can hope to conduct theoretical studies of reconnection and evaporation together only for cases of conduction-driven evaporation. While theoretical models exists for non-thermal electron energization by turbulence (Benz & Smith 1987; Hamilton & Petrosian 1992; Miller et al. 1996; Park et al. 1997) and by steady electric fields (Litvinenko 1996), no such model includes the process of magnetic reconnection self-consistently. Instead, most large-scale models of flare-related magnetic reconnection have been formulated using resistive MHD: fluid equations lacking a non-thermal electron population (Forbes & Priest 1983; Mikic et al. 1988; Amari et al. 1996; Magara et al. 1996; Nishida et al. 2009; Birn et al. 2009). In fluid models of fast reconnection, magnetic energy is ultimately converted to kinetic energy of flows whose shocks produce heat (Petschek 1964; Soward 1982; Longcope et al. 2009), which is then is carried by thermal conduction to the chromosphere where it drives evaporation (Forbes et al. 1989; Tsuneta 1996). Several investigations have succeeded in accommodating all these effects in a single, self-consistent simulation (Yokoyama & Shibata 1997, 1998; Chen et al. 1999a, 1999b).

Due to the numerical difficulties of resolving both the coronal reconnection, the solar transition region (TR) and the field-align thermal conduction between them, few investigations have followed those of Yokoyama & Shibata (1997) or Chen et al. (1999b) to simulate reconnection and evaporation together. Instead, the vast majority of theoretical investigations study magnetic reconnection in isolation, omitting the process of evaporation. This is primarily done to achieve better resolution of the coronal physics of reconnection. There is some evidence, however, that electron acceleration, and thus presumably reconnection, occurs at densities significantly enhanced by evaporation (Veronig & Brown 2004; Jiang et al. 2006; Guo et al. 2012) and even continuing to increase by ongoing evaporation (Liu et al. 2006). Thus it seems unrealistic to consider either process, reconnection or evaporation, in isolation. Moreover, many of the observational signatures with which reconnection models must ultimately make contact, such as X-ray and EUV light-curves or spectra, are direct effects of evaporation and only indirect effects of reconnection.

In place of a full, direct numerical simulation of both processes together, the effects of evaporation on coronal reconnection could be investigated using flows imposed in an otherwise coronal simulation. Doing so would, however, require a simplified relation between reconnection energy release and characteristics of the evaporation, such as density and flow velocity. Previous numerical studies of evaporation were generally not aimed at developing such a relation; fewer still sought one for conductively driven evaporation. Among the few investigations of this kind, Fisher (1987) and Brown & Emslie (1989) each offer different analytic models for evaporation from intense beam deposition. Fisher (1989) developed an analytic model for chromospheric condensation driven by either beam or conductive energy deposition. Fisher et al. (1984) presents a quantitative upper bound for the evaporation velocity in terms of evaporation temperature. The data provided in that work showed, however, actual evaporation speeds falling below this upper bound by as much as an order of magnitude. In fact, the data presented by Fisher et al. (1984) suggest a more useful relation might exist between evaporation velocity and flare energy flux for which they offer no explanation.

The present work takes up the challenge of finding a simple relationship between flare energy release and the evaporation and condensation flows it generates. It adopts a purely fluid model, with energy transported by thermal conduction rather than energy beams. This is done firstly because, as explained above, fluid models remain the only kind that can provide a complete picture of flares from energy storage, to release, and to evaporation. A second advantage is to reduce the size of available parameter space. Non-thermal beams are characterized by an energy flux, a spectral index, and a lower-energy cut-off, and evaporation has been found to depend on all of these parameters (Fisher 1989). We show here that thermal energy transport is effectively characterized by a single parameter, the energy flux F. This would in turn be determined by the amplitude (i.e., Mach number) of the shocks generating the heat, which depends in turn on the current sheet at which reconnection occurs (Longcope et al. 2009, 2010). The result of the present investigation is a set of simple relationships between F and the density and velocity of both chromospheric evaporation and chromospheric condensation. These relationships provide a very simple means of using spectroscopic observations of emission from the TR and chromosphere, to infer properties of magnetic reconnection in solar flares.

In order to facilitate the variation of parameters, this investigation uses a simplified model of the pre-flare TR. This is described, along with the governing dynamical equations, in the next section. Section 3 describes a typical numerical solution to these equations. This motivates a simplified analytical model, presented in the same section, and then compared to the numerical solution in detail. Section 4 present a large number of different numerical solutions spanning parameters. The analytical model is used to propose scaling relations which can be fit to the runs. Section 5 then applies these same scaling laws to simulations with more realistic TR and chromosphere treatments, including results from the literature. In nearly every case these results are found to conform reasonably well to the simple scaling laws.

2. THE MODEL

In this work, as in many previous investigations (Nagai 1980; Somov et al. 1981; Peres et al. 1982; MacNeice et al. 1984), evaporation is studied using one-dimensional gas-dynamic equations for evolution of the plasma along a static flux tube, parameterized by length ℓ. Equations for mass, momentum, and energy conservation in the flux tube of uniform cross section are

$\begin{eqnarray} {\partial \rho \over \partial t} + {\partial \over \partial \ell }(\, \rho v\,) &=& 0, \end{eqnarray} \tag{ 1 }$

$\begin{eqnarray} \rho \left({\partial v\over \partial t} + v{\partial v\over \partial \ell }\right) &=& -{\partial p\over \partial \ell } + {\partial \over \partial \ell }\left({{4\over 3}}\mu \, {\partial v\over \partial \ell } \right) + \rho g_{\parallel }, \end{eqnarray} \tag{ 2 }$

$\begin{eqnarray} {c}_{_{\rm V}}\rho \left({\partial T\over \partial t} +v{\partial T\over \partial \ell } \right) &=& \!-\! p\,{\partial v\over \partial \ell } \!+\! {{4\over 3}}\mu \left({\partial v\over \partial \ell }\right)^2 \!+\! {\partial \over \partial \ell }\left(\kappa \, {\partial T\over \partial \ell } \right) \!+\! \dot{Q}, \nonumber \\ \end{eqnarray} \tag{ 3 }$

where ρ and T are mass density and temperature, respectively, g_∥ is the component of gravitational acceleration along the tube, and ${c}_{_{\rm V}}$ is the specific heat, described below. The fluid moves along the flux tube at velocity v. Heating and cooling due to the flare and radiative losses are included in the volumetric power $\dot{Q}$ which is described in detail in the next sections. The total plasma pressure p is found using the ideal gas equation,

$\begin{equation} p = {k_{_{\rm B}}\over \bar{m}}\, \rho \, T, \end{equation} \tag{ 4 }$

where $k_{_{\rm B}}$ is Boltzmann's constant and $\bar{m}$ is the mean mass per particle.

Equations (1)–(3) are solved within a region spanning the corona and the chromosphere, over which varying states of ionization will occur. The focus is, however, on flare-driven flows which occur in plasmas heated to temperatures around or above 10⁶ K, and for which the plasma is expected to be fully ionized. We therefore adopt the expedient measure of solving Equations (1)–(3) for a fully ionized plasma with coronal abundance, setting $\bar{m} = 0.593\, m_p$ , where m_p is the proton mass, and taking the specific heat

$\begin{equation} {c}_{_{\rm V}}= {3\over 2}{k_{_{\rm B}}\over \bar{m}}. \end{equation} \tag{ 5 }$

The main way that ionization would affect our results would be by removing energy from the condensation shock (CS) through an additional contribution to ${c}_{_{\rm V}}$ . We discuss this in an appendix and show that such a contribution would have relatively minor effects on our results.

The thermal conductivity κ depends on temperature according to the classical Spitzer form (Braginskii 1965)

$\begin{equation} \kappa = \kappa _0\, T^{5/2},\quad \kappa _0 = 10^{-6} \, {\rm erg\,cm^{-1}\,s^{-1}\,K^{-7/2}}. \end{equation} \tag{ 6 }$

While this classical form can become inapplicable for extremely large heat fluxes, we have found that it applies to cases where heating occurs self-consistently through shocks (Guidoni & Longcope 2010; Longcope & Bradshaw 2010). In order to simplify our analysis we use it for all of our computations. The parallel component of dynamic viscosity μ has a temperature dependence identical to classical conductivity, but is lower by a factor proportional to the Prandtl number Pr,

$\begin{equation} \mu = {\rm{Pr}}\, {\kappa \over {c}_{_{\rm V}}}. \end{equation} \tag{ 7 }$

The Prandtl number in fully ionized plasma is Pr = 0.012, reflecting the very small electron-to-ion mass ratio. While it is, in some sense, much smaller than thermal conductivity, its effects are very important, especially in flares where T has increased by an order of magnitude (Peres & Reale 1993). Since they thermalize bulk kinetic energy, shocks cannot be modeled without viscosity. Nevertheless, some flare investigations in the literature do solve gas dynamic equations without viscosity (Cheng et al. 1983). To make closer contact with that work, and to sharpen the shocks for easier identification, we perform most runs using an artificially small value Pr = 10⁻⁴. To prevent under-resolution and maintain energy conservation, we place a lower bound on the viscosity

$\begin{equation} \mu > 0.05\, \rho \, c_{\rm s}\Delta \ell, \end{equation} \tag{ 8 }$

where Δℓ is the local grid spacing and c_s is the local sound speed. This bound is typically assumed in the regions of lowest T, such as the pre-flare chromosphere.

Equations (1)–(3) are solved using a Lagrangian code similar to that described in Guidoni & Longcope (2010). In this new code the thermal conductivity term in Equation (3) is differenced implicitly to permit larger time steps even in the face of extremely high temperatures. In contrast to Guidoni & Longcope (2010), magnetic tension has been omitted from Equation (2) under the assumption that the flux tube is already in magnetostatic equilibrium. This omission removes magnetic reconnection as an explicit energy source. We follow, instead, the standard practice of replacing the omitted term with an ad hoc heating, without a corresponding force. In future work we will return to explore the effect of making this assumption compared to allowing reconnection to energize the flare loop legitimately.

2.1. The Initial Condition—a Simplified Transition Region

The initial condition for the model is a hydrostatic loop of total length L, consisting of a coronal portion, TRs, and chromospheres at each footpoint. We add to the heating term, $\dot{Q}$ , an ad hoc contribution, localized to the center of the loop (i.e., the loop top), to drive the evaporation and condensation via thermal conduction. The aim of this study is to characterize the evaporation and condensation as responses to the loop-top energy input. This characterization will depend, to some degree, on the pre-flare structure of the TR and chromosphere.

In equilibrium models, the structure of the TR and chromosphere is determined by the interplay between heating, thermal conduction, and radiation. The relative magnitude of these competing effects can be estimated from the minimum value of the logarithmic differential emission measure (DEM), $\xi =n_e^2dz/d\ln T$ , observed to be about ξ_min ∼ 10²⁷ cm⁻⁵ at the TR temperature T_tr ∼ 10⁵ K (Dere 1982). The minimum DEM corresponds to conductive flux of

$\begin{equation} |F_c| \simeq {\kappa _0\over 4k_{_{\rm B}}^2}{p^2T_{\rm tr}^{3/2}\over \,\xi _{\rm min}} = (4\times 10^5\,{\rm cm^4\,erg^{-1}})\,p^2, \end{equation} \tag{ 9 }$

where p is the equilibrium pressure at T_tr. Since active region loops have equilibrium pressures generally less than 10 erg cm⁻³, energy fluxes are of the order of 10⁷ erg cm⁻² s⁻¹ or less. Energy flux typically attributed to flares, above 10⁹ erg cm⁻² s⁻¹, completely overwhelms this. The most important factor determining the chromospheric response is therefore the distribution of mass, rather than the processes which created that mass distribution in the first place.

Motivated by the foregoing argument we wish to explore how different TR mass distributions affect the evaporative response to flare heating. Toward this end we perform a series of experiments initialized with an artificial TR, of set thickness Δ_tr, separating perfectly uniform coronal and chromospheric plasmas (Brannon & Longcope 2014). In these experiments gravitational stratification is omitted (g_∥ = 0) so the equilibrium has uniform pressure p₀. Radiative losses are also omitted. Coronal plasma at uniform temperature T_{co, 0} is maintained by heating functions at ℓ = Δ_tr and ℓ = L − Δ_tr. Temperature drops to a pre-set chromospheric level, T_{ch, 0}, across the TRs. This is maintained by cooling functions, at ℓ = 0 and ℓ = L, equal and opposite to the heating functions. This construction is implemented by setting $\dot{Q}$ in equilibrium to

$\begin{equation} H_{\rm eq}(\ell) = A[S(\ell -\Delta _{\rm tr}) - S(\ell) + S(L-\Delta _{\rm tr}-\ell) - S(L-\ell) ], \end{equation} \tag{ 10 }$

where the heating and cooling are distributed over a distance 2w according to the shape function

$\begin{equation} S(x) = \left\lbrace \begin{array}{lcl}1 - (x/w-1)^2, & & 0 < x < 2w \\ 0, & & \hbox{otherwise}\end{array} \right.\!\!\!\!, \end{equation} \tag{ 11 }$

centered at x = w and vanishing at x = 0. The coefficient A is chosen to obtain the desired temperature ratio, T_{co, 0}/T_{ch, 0} = R_tr, across each TR.

The initial condition is a static, isobaric equilibrium with $\dot{Q}=H_{\rm eq}(\ell)$ . From Equation (3), we find the temperature distribution of the equilibrium

$\begin{equation} T^{7/2}(\ell) = T_{\rm{ch},0}^{7/2} - \displaystyle \frac{7}{2}\kappa _0^{-1}\, \int ^{\ell }\,d\ell ^{\prime }\int ^{\ell ^{\prime }}\,d\ell ^{\prime \prime } H_{\rm eq}(\ell ^{\prime \prime }). \end{equation} \tag{ 12 }$

An example, shown in Figure 1, has T_{co, 0} = 2 × 10⁶ K, T_{ch, 0} = 2 × 10⁴ K, Δ_tr = 3 Mm, and w = 0.75 Mm. The temperature changes only within the range, 0 < ℓ < Δ_tr + w = 4.5 Mm, but mostly at the outside edge of the cooling function (i.e., the left side or the left TR);¹ this point is used to define the total length of the loop. The flux between the heating and cooling sections is |F_c| ≃ 10⁷ erg cm⁻² s⁻¹, which is slightly higher than in actual TRs, but still far below the flaring energy flux. This heat flux is lower than the free-streaming heat flux limit, $F_c^{\rm (fs)}=({3/2})m_en_ev_{{\rm th},e}^3$ (Campbell 1984) by more than an order of magnitude throughout the TR.

The Lagrangian grid points are arranged to keep roughly constant mass between them. This has the effect of concentrating grid points in the chromosphere. The initial locations are shown along the bottom of the top panel of Figure 1.

2.2. Flare Heating

Flare energy release is modeled using a flat-topped ad hoc heating function centered at the loop top and extending a total distance of Δ_fl,

$\begin{equation} H_{\rm fl}(\ell) = \left\lbrace \begin{array}{@{}lcl}\displaystyle {F\over \Delta _{\rm fl}/2}, && |\ell - L/2 | < \Delta _{\rm tr}/2 \\ 0, & & \hbox{otherwise}\end{array} \right.\!\!\!\!. \end{equation} \tag{ 13 }$

This is added to the equilibrium heating function, $\dot{Q}=H_{\rm eq}+H_{\rm fl}$ , to provide the term in Equation (3). That contribution is held fixed for the entire run.

Since the solution begins with a state which is in equilibrium with heating H_eq, the addition of H_fl has the effect of ramping the flare energy input instantaneously. In models where flare energy is released by fast magnetic reconnection the heating occurs at a slow shock which is generated by the retraction of the reconnected flux tube (Petschek 1964; Soward & Priest 1982). A particular flux tube passes the reconnection point in an instant and thereafter begins retracting (Longcope et al. 2009). It is therefore appropriate to replace this kind of self-consistent heating with an ad hoc heating term with an instantaneous turn-on. Since evaporation is the focus of this work, the heating term is never turned off.

The parameter F defines the total energy flux delivered to one side of the loop, and ultimately to a single footpoint. This is one of the key parameters expected to dictate the chromospheric response. We perform runs for a range of parameters. The most significant variations occur for variations in energy flux F, loop length L, equilibrium values of pressure p₀ and TR temperature ratio R_tr = T_{co, 0}/T_{ch, 0}.

3. THE SIMULATIONS

Figure 2 shows the solution for a flux tube of length L = 53 Mm, pressure p₀ = 1.0 erg cm⁻³, subjected to a flare energy flux F = 3.5 × 10¹⁰ erg cm⁻² s⁻¹, distributed over Δ_fl = 10 Mm. Its initial TR was the same one shown in Figure 1. The heating produces an evolution resembling those reported in many previous simulations of this kind (Cheng et al. 1983; Emslie & Nagai 1985). The temperature within the heated region (between the vertical dashed lines) rises rapidly. At the same time, thermal conduction creates fronts moving rapidly in both directions. By t = 1.5 s these fronts have reached the TRs and the loop top temperature slows its rise. After t = 2.0 s the apex temperature has achieved a steady state of T_fl = 3.7 × 10⁷ K. After that the temperature profile changes very little over the bulk of the loop. The heat flux F is being delivered to the chromosphere where it generates evaporation flow.

The loop's DEM undergoes a characteristic evolution during as the conduction front travels to the TR and initiates evaporation. Figure 3 shows the logarithmic DEM, $\xi (T)=n_e^2d\ell /d\ln T$ , throughout the evolution depicted in Figure 2. The initial DEM has a sharp peak at the coronal temperature (T = 2 × 10⁶ K) and a sloping TR below. (The simplified TR creates a minimum min (ξ) ≃ 3 × 10²⁷ cm⁻⁵, slightly higher, and at higher temperature, than a real TR.) The conduction front shifts the emission from the coronal peak into an increasing slope up to the flare peak at T ≃ 3 × 10⁷, making no change to the TR which has not yet been disturbed. As the conduction front crosses the TR, over the times 1.3 ⩽ t ⩽ 1.5 s, the DEM below the corona is rapidly eroded away: a steep drop reaches T = 7 × 10⁵ K and T = 2 × 10⁵ K at t = 1.3 and 1.4 s, respectively. This emission is piled into a peak around 3 × 10⁶ K. At later times the TR has been thinned enough that ξ ≃ 3 × 10²⁷ cm⁻⁵, and the peak at 3 × 10⁶ K grows by evaporation.

The evaporation phase begins once the conduction front has reached the TR and the peak coronal temperature has plateaued: by about t = 2.0 s in the present case. Figure 4 shows details of this phase from t = 3.5 s when the conduction front has penetrated 1.0 Mm into the chromosphere. This penetration occurs at the head of the conduction front which takes the form of a downward propagating shock of extremely high Mach number: the pressure jumps by a factor 480 and the density jump, a factor of 3.85, is nearly 4, the maximum value permitted by Rankine Hugoniot relations. The post-shock material is downflowing at v_c ≃ −370 km s⁻¹, in a chromospheric condensation (Fisher et al. 1985a); we henceforth refer to this shock as the CS. An upward propagating shock, the evaporation shock (ES), has reached ℓ = 1.5 Mm by this time. The density there jumps by a factor slightly greater than four and the peak velocity is v_e = +750 km s⁻¹. The shock is well defined because we have used the anomalous Prandtl number Pr = 10⁻⁴. Between these two shocks lies a rarefaction wave (RW) over which the velocity changes continuously between downflow and upflow, and the density increases by a factor close to 100.

**Figure 4.** Structure of the evaporation flow at the last time shown in Figure 2. The right column shows the entire left side of the flux tube and the left column shows the region around the TR. The rows shows, from top to bottom, pressure, velocity, temperature, and electron density. The dotted curve shows the initial profiles of temperature and density (initial pressure and velocity are both uniform).
Download figure:
Standard image High-resolution image

3.1. An Analytic Model of the Evaporation

The foregoing simulation suggests a simple model for the evaporation. The conduction front overtakes the TR rapidly enough that the entire region is raised to a uniform temperature, T_*, before any dynamical response can occur. The TR then functions as an initial pressure discontinuity whose subsequent evolution is an isothermal Riemann problem. Since the TR is now at uniform temperature, the pressure ratio across the jump matches the density ratio which is inversely related to the temperature ratio across the pre-flare TR,

$\begin{equation} {p_{\rm{ch}}\over p_{\rm{co}}} = {\rho _{\rm{ch},0}\over \rho _{\rm{co},0}} = {T_{\rm{co},0}\over T_{\rm{ch},0}} = R_{\rm tr}. \end{equation} \tag{ 14 }$

Under isothermal dynamics, the initial pressure jump decomposes into a shock and a RW (Fabbro et al. 1985). If the initial jump is thin enough the RW will be self-similar (Landau & Lifshitz 1959; Fabbro et al. 1985)

$\begin{eqnarray} v(\ell,t) &=& {\ell \over t} + a, \end{eqnarray} \tag{ 15 }$

$\begin{eqnarray} \rho (\ell,t) &=& \rho _0\, \exp \left(-{\ell \over at} \right) = \rho _0\, \exp \left(1 - {v\over a} \right), \end{eqnarray} \tag{ 16 }$

where $a=\sqrt{p/\rho }=\sqrt{k_{_{\rm B}}T_*/\bar{m}}$ is the isothermal sound speed and ρ₀ is a constant determined from the surrounding solution.

In the complete solution, illustrated in Figure 5, the left edge is the CS, which heats the plasma to its uniform temperature of T_*, with (adiabatic) sound speed $c_{s,*}=\sqrt{5/3}\,a$ . Because the CS has an extremely high Mach number it increases the density to 4ρ_{ch, 0} = 4R_tr ρ_{co, 0}. The post-shock flow speed for a hypersonic shock propagating into stationary plasma is

$\begin{equation} v_{c} = -{3\over \sqrt{5}}\,c_{s,*} = -\sqrt{3}\, a, \end{equation} \tag{ 17 }$

directed downward (v_c < 0). This is the condensation velocity, lying at the left-hand edge of the RW and labeled A in Figure 5. Applying these two conditions to point A allows ρ₀ to be eliminated from Equation (16)

$\begin{equation} {\rho (\ell,t)\over \rho _{\rm{co},0}} = 4R_{\rm tr}\, \exp \left(-\sqrt{3} - {v\over a} \right). \end{equation} \tag{ 18 }$

**Figure 5.** Schematic plot of the Riemann problem solution. The left column shows v(ℓ)/a (top) and ρ(ℓ) (bottom, on a logarithmic scale) vs. position ℓ. The right column shows the same quantities plotted against the other (i.e., v/a vs. ρ, in the upper right). The dash dotted line beginning at A and passing though B shows the equation for an isothermal, self-similar RW, namely, Equation (18). The dotted curves show the relation for an isothermal shock. These intersect at point B, marking the rightmost point in the RW.
Download figure:
Standard image High-resolution image

The ES propagates into the stationary coronal plasma to the right. In the reference frame of the ES, the coronal material is flowing leftward with isothermal Mach number, $M^{(it)}_{\rm{es}}$ . Because it is an isothermal shock the post-shock flow, in the ES reference frame, has isothermal Mach number $1/M^{(it)}_{\rm{es}}$ . In the non-moving reference frame, the evaporation flows at velocity v_e to the right, while the pre-shock material is stationary. The velocity difference between pre-shock and post-shock velocities is independent of reference frame so

$\begin{equation} {v_e\over a} = M^{(it)}_{\rm{es}} - {1\over M^{(it)}_{\rm{es}}}. \end{equation} \tag{ 19 }$

The density jump across the isothermal shock is $\rho _{e}/\rho _{\rm{co},0}=[M^{(it)}_{\rm{es}}]^2$ .

The right-hand edge of the RW, labeled B in Figure 5, must match the conditions of the ES described above. Applying these conditions to Equation (18) yields the relation

$\begin{equation} \left [M^{(it)}_{\rm{es}}\right ]^2 = 4R_{\rm tr}\,\exp \left(-\sqrt{3} - M^{(it)}_{\rm{es}} + {1\over M^{(it)}_{\rm{es}} } \right), \end{equation} \tag{ 20 }$

which must be satisfied by $M^{(it)}_{\rm{es}}$ . Figure 6 shows the solution for a range of TR density ratios R_tr. Diamonds follow an empirical fit

$\begin{equation} M^{(it)}_{\rm{es}} \simeq 2.670 + 1.209\, \log (R_{\rm tr}/100)\, [ 1 + 0.126\, \log (R_{\rm tr}/100) ]. \end{equation} \tag{ 21 }$

Due to its logarithmic dependence on R_tr the isothermal Mach number of the ES falls within the narrow range between 2.5 and 3.5 for most reasonable assumptions above the pre-flare TR.

**Figure 6.** Isothermal Mach number of the evaporation shock, $M^{(it)}_{\rm{es}}$ , as a function of R_tr (solid curve), satisfying Equation (20). Diamonds follow the empirical fit, Equation (21). The dashed curve shows the upper bound, Equation (22), derived by Fisher et al. (1984).
Download figure:
Standard image High-resolution image

**Figure 6.** Isothermal Mach number of the evaporation shock, $M^{(it)}_{\rm{es}}$ , as a function of R_tr (solid curve), satisfying Equation (20). Diamonds follow the empirical fit, Equation (21). The dashed curve shows the upper bound, Equation (22), derived by Fisher et al. (1984).
Download figure:
Standard image High-resolution image

This model can be compared to an upper bound obtained by Fisher et al. (1984). They used the isothermal momentum equation under the assumption that the entire pressure drop matched the pre-flare values, and the kinetic energy in the condensation was ignorable. By doing so, they obtained the upper bound on the evaporation velocity,

$\begin{equation} v_{_{\rm FCM}} = \sqrt{(6/5)\,\ln R_{\rm tr}}\,c_{s,*} = \sqrt{2\ln R_{\rm tr}} \, a. \end{equation} \tag{ 22 }$

This is plotted as a dashed line in Figure 6. The discrepancies are due to the different assumptions used in the two derivations.

3.2. Modeling the Simulation

Figure 7 shows the numerical solution from Figure 4, at a somewhat later time (t = 6.5 s), plotted in the same manner as the sketch in Figure 5. Because the solution includes viscosity, albeit artificially low, the ES is somewhat broadened. The CS satisfies the structure assumed in the simplified model: ρ_c = 4R_trρ_{co, 0} and $v_c=-\sqrt{3}\, a$ . Its actual location, plotted with a square, lies almost on top of this theoretical one, marked with a triangle in the right column of Figure 7. The ES also conforms to the assumption of an isothermal shock, marked by a dotted curve in the right column.

The actual RW does not, however, conform to the self-similar, isothermal structure. The velocity is not linear with distance, and the density is not exponential. The path between squares, in the right column of Figure 7, falls off the isothermal model shown by dash-dotted lines. It is temperature variation along the RW that causes its track to fall beneath the isothermal (dash-dotted) one in the lower right. This curve intersects the shock curve at a point of lower density and lower velocity. Thus the prediction of the isothermal model, i.e., Equation (20), provides an upper bound on the isothermal Mach number of the actual shock. For this case with R_tr = 100, the isothermal model predicts $M^{(it)}_{\rm{es}}=2.67$ , and thus ρ_e/ρ_{co, 0} = (2.67)² = 7.13 and v_e/a = 2.67 − (2.67)⁻¹ = 2.30. The actual values, 4.69 and 1.50, respectively, fall below those values. The isothermal Mach number of this shock, $M^{(it)}_{\rm{es}}=\sqrt{4.69}=2.16$ , lower than the prediction from Equation (20). This is the effect of temperature variation along the RW.

At higher values of heat flux F the evaporation/condensation structure becomes more isothermal owing to the larger conductivity at higher temperatures. A solution of this kind, shown in Figure 8, has a RW falling closer to the isothermal model (dash-dotted curve). At these higher fluxes, however, the RW begins to affect the CS, and the latter departs from the assumption of a hypersonic shock, as evident in the separation between square and triangle in the right column of plots. The density behind the CS is lower than 4ρ_{ch, 0}, and the velocity greater than $-\sqrt{3}a$ . We discuss in the next section the level of energy flux required to produce this departure from the analytic model.

**Figure 8.** Solutions from a run like Figure 7 but with F = 5.5 × 10¹¹ erg cm⁻² s⁻¹, plotted at t = 2.3 s in the same manner as Figure 7.
Download figure:
Standard image High-resolution image

The analytic model breaks down for different reasons in the opposite limit of very small energy fluxes. In this case, the conduction fronts move very slowly from the loop top. In one case (not shown) with F = 10⁹ erg cm⁻² s⁻¹, the steady-evaporation phase is achieved only after 30 s compared to 2 s, in Figure 2. The loop-top heating also drives flows downward, as evident on the right column of Figure 4. When the conduction front arrives late, this flow impinges on the evaporation flow. The result is a departure from the analytic model of Figure 5, and a systematically low evaporation flow.

Using the physical viscosity, characterized by Pr = 0.012, provides a more realistic, but less clear picture of the evaporation process. Figure 9 shows the results of making this change to the simulation discussed above. The larger viscosity spreads the ES considerably, causing it to overlap with the RW. As a result, the phase space curve of the ES deviates from the isothermal curve (dotted) by veering toward the RW curve (dash-dotted). The point of maximum velocity is adopted at the point separating these two structures, now largely merged. This occurs at a lower velocity (v_e/a = 1.20) and higher density (ρ_e/ρ_{co, 0} = 6.30) than in the case with a well-defined ES. In addition, the ES continues expanding, causing these values to change with time. We therefore continue using the anomalously low value, Pr = 10⁻⁴ to obtain clear values. We recognize that these values will differ from those with physical viscosity in the manner just described.

4. SCALING OF THE EVAPORATION

4.1. Variation in Flare and Loop Parameters

We next explore how the structure of the evaporation and condensation varies with parameters. Twenty-eight different simulations are performed with different values of L, F, p₀, and R_tr. Loop lengths range from 9 to 86 Mm, fluxes from 10⁹ to 10¹² erg cm⁻² s⁻¹, initial pressures from 0.3 to 3.0 erg cm⁻³, and R_tr from 40 to 100. For simplicity we hold fixed the TR structure by keeping the values Δ_tr = 3 Mm, w = 750 km; the chromosphere is kept at T_{ch, 0} = 20, 000 K, but the corona is at T_{co, 0} = R_trT_{ch, 0} which varies as R_tr is varied. The heating profile was also kept fixed with Δ_fl = 10 Mm. Each simulation was run past the time the apex temperature plateaued at T_fl. At that point, the shocks, ES and CS, were identified and characterized.

The peak apex temperature, T_fl, is achieved when the flare energy flux balances the power driving the evaporation. The heat flux reaching the evaporating plasma, at lower temperature, will be ${\sim }\kappa _0T_{\rm fl}^{7/2}/L$ . Equating this with the input flux yields an expected scaling

$\begin{equation} T_{\rm fl} = C_T\, ({FL}/\kappa _0)^{2/7}, \end{equation} \tag{ 23 }$

where C_T is a dimensionless constant. The same form is given in the Appendix of Fisher (1989). The left panel of Figure 10 shows T_fl for the 28 runs (plusses) versus the product FL, where L is the full loop length. The dashed line shows Equation (23) with C_T = 1.46, found from a fit to all 28 runs.

**Figure 10.** Flare temperature T_fl (left) and ES density enhancement ratio R_es (right) of flare simulations. The plusses are the runs with artificial TR described in Section 4. Other symbols show runs with more realistic TR treatments from the literature and this work, as described Section 5. The left panel shows T_fl, in Kelvin, vs. FL (in units of 10²⁰ erg cm⁻¹ s⁻¹), with a dashed line showing Equation (23). The right panel shows R_es vs. $F\,p_0^{-5/2}$ (in units of 10⁸ erg^−3/2 cm^−9/2 s⁻¹), and the dashed line shows relation (25).
Download figure:
Standard image High-resolution image

According to the analytic model of Section 3.1, the isothermal Mach number of the ES, $M^{(it)}_{\rm{es}}$ , depends only on the pre-flare chromospheric/coronal density ratio R_tr (and on that only logarithmically). Application to an actual solution showed, however, that temperature gradient within the RW led to a slightly lower value of $M^{(it)}_{\rm{es}}$ . For the cases of small viscosity, however, the ES is a simple isothermal shock and therefore its density enhancement is

$\begin{equation} R_{\rm{es}} = {\rho _e\over \rho _{\rm{co},0}} = [M^{(it)}_{\rm{es}} ]^2. \end{equation} \tag{ 24 }$

This will have an upper bound dependent on R_tr, and a value depending on other factors responsible for the temperature variation within the RW. Figure 10 shows that the observed ES density enhancement, R_es, is ordered by the quantity, $F/p_0^{5/2}$ , identified using multivariate regression to the 28 plusses. An empirical relation,

$\begin{equation} R_{\rm{es}} = \displaystyle\frac{1}{2} \ln \left(F\,p_0^{-5/2}/C_M\right) + 7.13, \end{equation} \tag{ 25 }$

lies near, and slightly above, much of the data, where C_M = 2.8 × 10¹² erg^−3/2 cm^−9/2 s⁻¹. Expression (21) predicts that the upper bound to the density enhancement of, R_es = (2.67)² = 7.13. The empirical fit, Equation (25), would exceed the maximum for $F\,p_0^{-5/2}>C_M$ . This value is above the point where the analytic model becomes untenable due to the reduction in Mach number of the CS, as shown in Figure 8.

The analytic solution shown in Figure 5 has a total kinetic energy (per unit area)

$\begin{equation} E_K = \int \displaystyle \frac{1}{2}\rho v^2\, d\ell = \frac{1}{2}\rho _{\rm{co},0}\, a^3\, t\, \int {\rho (\tilde{\ell })\over \rho _{\rm{co},0}}\, [M^{(it)}(\tilde{\ell }) ]^2d\tilde{\ell }, \end{equation} \tag{ 26 }$

where $\tilde{\ell } = \ell /at$ is the similarity variable of the solution, and M^(it) = v/a is the scaled velocity plotted in the figure. The integral in the final expression depends on the scaled analytic solution and therefore on $M^{(it)}_{\rm{es}}$ only. Assuming the flare energy flux, F, exactly balances this energy requirement (i.e., F = dE_K/dt), the isothermal sound speed within the evaporation flow is

$\begin{equation} a \sim ({F/\rho _{\rm{co},0}})^{1/3}. \end{equation} \tag{ 27 }$

If we neglect the logarithmic dependence of $M^{(it)}_{\rm{es}}$ , the evaporation flow speed will scale in the same manner

$\begin{equation} v_e = C_e\,(F/\rho _{\rm{co},0})^{1/3}. \end{equation} \tag{ 28 }$

The best fit to the 28 plusses, shown in Figure 11, is for C_e = 0.38. The points at the very left and very right of the range trend below this fit. The explanation may be the departure from the analytic model at high and low fluxes, described briefly in Section 3.2.

Fisher (1989) estimates the height of the pressure peak separating evaporation from condensation to scale as p_e ∼ (F²ρ_{co, 0})^1/3 (see also Fisher et al. 2012). While he does not go on to estimate the evaporation flow speed we can take the additional step of identifying the sound speed within the evaporation as $a\sim \sqrt{p_{e}/\rho _{\rm{co},0}}$ to obtain a scaling identical to Equation (27). While Fisher's logic differs from that following Equation (26) above, both produce the same scaling, as would simple dimensional analysis (Fisher 1989).

It is noteworthy that we found CSs in every run, including those of the smallest flare energy fluxes. Thermal conduction therefore differs from non-thermal deposition, which generates condensation only when the flux exceeds a threshold F_thr. This threshold is found to be between F_thr = 3 × 10⁹ and 3 × 10¹¹ erg cm⁻² s⁻¹ depending on pre-flare pressure and the beam's spectral properties (Fisher 1989). Non-thermal energy fluxes below the threshold produce gentle evaporation, in which there is little observable condensation. Fluxes above the threshold drive explosive evaporation, accompanied by chromospheric condensation. We find that thermal conduction always leads to explosive evaporation; a fact that has been noted before (Fisher 1989).

If the isothermal model of Section 3.1 applied all the way to the CS, then the condensation velocity would scale as v_c ∼ a ∼ (F/ρ_{co, 0})^1/3, just like the evaporation velocity. The significant temperature variation within the RW, however, invalidates this conclusion. Instead, using multivariate regression to the 28 plusses, we find an empirical scaling

$\begin{equation} v_c = C_c\,(F/\rho _{\rm{ch},0})^{1/2}, \end{equation} \tag{ 29 }$

where C_c = 10⁻⁴ cm^−1/2 s^1/2, as shown by the dashed line in the right panel of Figure 11.

Fisher (1989) used a different analytical approach to find a comparable scaling for condensation velocity. For the case of thermal conduction or beam heating lower than the explosive threshold (so-called gentle evaporation) he finds condensation velocity (given in his Equation (34a))

$\begin{equation} v^{\rm (F89)}_{c} \simeq 0.4\left({F/\rho _{\rm{ch},0}}\right)^{1/3}, \end{equation} \tag{ 30 }$

although his expression identified by F only that portion of the energy flux deposited above the CS (this relation is plotted with a dashed line on the right panel of Figure 11). The scaling and pre-factor of Equation (30) match our expression for evaporation velocity, Equation (28). Fisher's derivation does not match ours precisely, although it is also possible to find the scaling from dimensional analysis alone. Curiously, the condensation velocity data, on the right panel of Figure 11, is well fit by the scaling F/ρ_{ch, 0} to a one-half power rather than the one-third power of Equation (30) shown as a dotted line. The discrepancy may result from a systematic variation in the fraction of the total flux F reaching the CS.

4.2. Variation in TR

In order to provide maximum flexibility, a simplified TR model has been used under the presupposition that the detailed TR structure does not affect the evaporation or condensation. Figure 12 shows that this presupposition is warranted. Velocity and density are shown at t = 4 s from solutions with different TR structures. The thickness of the TR is varied from w = 1 Mm to w = 3 Mm. The size of the heat sources range from Δ_tr = 0.25 Mm to Δ_tr = 0.75 Mm. It is evident that the resulting flow structure is insensitive to the structure of the initial TR.

**Figure 12.** Evaporation flows at t = 4 s from runs with different TR models. Velocity v in km s⁻¹ (left column) and the electron density *n_e* (right column) are plotted for each run, but the x axis is shifted to line up the plots. Except for w and Δ_tr, all have the same parameters as those in Figure 4. The top row (a) shows that same run, with w = 3 Mm and Δ_tr = 0.75 Mm. The next two rows have w = 3 Mm and Δ_tr = 0.25 Mm (b), and w = 1 Mm and Δ_tr = 0.25 Mm (c). The bottom row (d) is the same as (a) but with a specific heat ${c}_{_{\rm V}}$ modified to simulate ionization, as explained in the Appendix. Dotted lines in the left column repeat the velocity profile for row (a) for reference.
Download figure:
Standard image High-resolution image

**Figure 12.** Evaporation flows at t = 4 s from runs with different TR models. Velocity v in km s⁻¹ (left column) and the electron density *n_e* (right column) are plotted for each run, but the x axis is shifted to line up the plots. Except for w and Δ_tr, all have the same parameters as those in Figure 4. The top row (a) shows that same run, with w = 3 Mm and Δ_tr = 0.75 Mm. The next two rows have w = 3 Mm and Δ_tr = 0.25 Mm (b), and w = 1 Mm and Δ_tr = 0.25 Mm (c). The bottom row (d) is the same as (a) but with a specific heat ${c}_{_{\rm V}}$ modified to simulate ionization, as explained in the Appendix. Dotted lines in the left column repeat the velocity profile for row (a) for reference.
Download figure:
Standard image High-resolution image

5. APPLICATION TO MORE REALISTIC SIMULATIONS

The foregoing considered simulations with an artificial TR in order to fully explore the dependence of evaporation on parameters. In cases with more realistic TRs these parameters cannot be varied independently. The results derived above do, however, apply to these cases.

5.1. Simulation with Radiation

We use the same code to perform runs with a more realistic TR by including radiation and gravitational stratification. The heating function in this case is

$\begin{equation} \dot{Q} = -n_e^2\Lambda (T) + H_{\rm eq}(\ell) + H_{\rm fl}(\ell), \end{equation} \tag{ 31 }$

where n_e = 0.874ρ/m_p is the electron number density assuming complete ionization. The radiative loss function Λ(T) is a piece-wise power-law fit to the output of Chianti 7.1 with coronal abundances (Landi et al. 2013; Dere et al. 1997). We make this function monotonic, Λ ∼ T^1.66, over the range 19, 000 K < T < 78, 700 K. This omits a peak at T ≃ 20, 000 K, due to silicon which would make the lowest temperature regions susceptible to radiative instabilities. For T < 19, 000 K we set Λ = 0 in order to prevent any thermal instability in the isothermal chromosphere. The equilibrium heating, H_eq(ℓ) is taken to be uniform except in the chromospheric layer, where it is set to exactly balance the radiative losses in equilibrium. In order to produce a stratified chromosphere, we take g_par = ±274 m s⁻² within the chromosphere only.

The initial state is taken to be an equilibrium for heating $\dot{Q}=-n_e^2\Lambda +H_{\rm eq}$ with specified peak temperature T_{co, 0} = 2 × 10⁶ K, a chromosphere at T_{ch, 0} = 20, 000 K, and (full) loop length L = 60 Mm. Following the arguments of Rosner et al. (1978), specifying these parameters fixes the equilibrium heating, and coronal pressure, to be H_eq = 1.3 × 10⁻³ erg cm⁻³ s⁻¹ and p₀ = 0.65 erg cm⁻³. We add the flare heating profile H_fl from Equation (13), with F = 3.4 × 10⁹ and 10¹⁰ erg cm⁻² s⁻¹ in separate runs. The characteristics of the resulting evaporation and condensation, plotted as ×'s in Figures 10 and 11, fall very close to those of runs with artificial TRs. Notably, the radiative losses decreased the flare temperature (by about 25% below Equation (23) in both cases), and the ES density enhancement (by about 10% below Equation (25) in both cases); radiative losses left the evaporation velocity within 8% of relation (28). This corroborates our assumption that radiative losses do not make significant contributions to the rapid evolution of flare evaporation, at least when it is driven by conduction.

5.2. Simulations in the Literature

Numerous simulations of conductively driven flaring loops have been reported in the literature. Each of these used a different code, implementing different aspects of the relevant physics. Evaporation and condensations characteristics from some notable simulation reports are given in Table 1 and plotted as various symbols on Figures 10 and 11. These lie near enough to the dashed lines to be considered to be adequately predicted by relations described in previous section.

Table 1. Flare Simulations from the Literature

						ES^f
	L^a	F^b	$p_0^{\rm c}$	$t_{\rm shock}^{\rm d}$	$T_{\rm fl}^{\rm e}$	$n_1^{\rm g}$	$n_2^{\rm h}$	R_es	v_e	v_c
	(Mm)	(10¹⁰)		(s)	(MK)	(10⁹cm⁻³)			(km s⁻¹)	(km s⁻¹)

Cheng et al. (1983)	10	0.2	2.9	4.95	11ⁱ	1.6	4.3	2.7	200	40
Emslie & Nagai (1985)	21	0.67	0.91	10	20	0.64	2.9	4.6	370	60
	21	1.3	0.91	20	28	0.25	1.2	4.8	670	30
MacNeice (1986)	24	0.1	2.9	13.75	10	1.3	2.5	2.0	60	20
Peres & Reale (1993), r1	38	0.63	6.0	30	22	0.65	1.5	2.3	400	40
r2	38	6.3	6.0	⋅⋅⋅	42	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	750	⋅⋅⋅

Notes. ^aFull length between chromospheric footpoints. ^bHeat flux in units of 10¹⁰ erg cm⁻² s⁻¹. ^cPre-flare pressure at base of corona in units of erg cm³. ^dTime at which evaporation properties were measured. ^eFlare temperature at loop apex at t_shock. ^fProperties of the evaporative shock (ES). ^gPre-shock electron density. ^hPost-shock electron density. ⁱElectron temperature.

Download table as: ASCII Typeset image

Peres & Reale (1993) used the Palermo-Harvard numerical code (Peres et al. 1982) to perform several simulations aimed at exploring the effects of viscosity on evaporation. Two sets of these runs, named run 1 and run 2, use a steady ad hoc heating source, localized to the loop top. Theirs had a 10 Mm wide Gaussian profile in place of the flat-topped function of Equation (13). The flaring energy flux,

$\begin{equation} F = \int _0^{L/2}\, H_{\rm fl}(\ell)\, d\ell, \end{equation} \tag{ 32 }$

is F = 6.3 × 10⁹ (run 1) and 6.3 × 10¹⁰ erg cm⁻² s⁻¹ (run 2) in their two cases. Their initial condition was an equilibrium, subject to optically thin radiative losses, and gravity. The semi-circular loop had a full length of L = 38 Mm, temperature and pressure T = 3 × 10⁶ K, and p₀ = 6 erg cm⁻³, at the loop apex. Their simulation solved single-fluid hydrodynamic equations with temperature-dependent ionization of hydrogen.

Figure 13, reproduced from Peres & Reale (1993), shows the time evolution of run 1 done with zero viscosity. The ES, propagating upward until after t = 40 s, is seen clearly due to the low viscosity. The curves from t = 30 s are used to deduce properties of the shock, including pre-shock and post-shock density, n₁ and n₂, and evaporation velocity v_e. The values found are listed in Table 1 and plotted as diamonds in Figures 10 and 11. Peres & Reale (1993) report, in a table, T_fl = 4.2 × 10⁷ and v_e = 750 km s⁻¹ for run 2. Since they provide no plot analogous to our Figure 13 for that run, however, we could not deduce R_e or v_c.

One of the earliest simulations of a flaring loop was done by Cheng et al. (1983) using the NRL Dynamic Flux Tube Model (Boris et al. 1980). This code solves quasi-neutral gas dynamics equations for a fully ionized plasma with distinct electron and ion temperatures. There is no viscosity. The electrons are heated directly by an ad hoc source with a Gaussian profile with full width 2 Mm. The heating is increased linearly with time so that the energy flux, defined according to Equation (32), is F = (4 × 10⁸ erg cm⁻² s⁻²)t. The tube is quite short, with full length L = 10 Mm between the chromospheric feet of the semi-circular loop. The state at t = 4.95 s, shown in their Figures 4 and 6, clearly shows an ES structure from which we deduce shock properties. We list these values in Table 1 and plot them as squares in Figures 10 and 11.

Emslie & Nagai (1985) simulate a flaring loop energized by either a beam of non-thermal electrons or a direct ad hoc heating. They use the code reported in Nagai (1980) solving single-fluid equations with temperature-dependent ionization of hydrogen, viscosity, and a sophisticated treatment of thermal conduction. The case with ad hoc heating, applicable to our present study, has a loop-top Gaussian profile of full width 4.8 Mm, which ramps up linearly as F = (6.7 × 10⁸ erg cm⁻² s⁻²)t. With a total length L = 22 Mm, their semi-circular loop is twice as long as that of Cheng et al. (1983). From their Figures 7 and 8, showing the state of the simulation, we can deduce the properties of the shock at t = 10 and 20 s. These values are Table 1 and plotted as triangles in Figures 10 and 11.

Finally, MacNeice (1986) performed a simulation with very high spatial resolution of a flare heated by a loop-top Gaussian profile (7 Mm wide) with constant heat flux F = 10⁹ erg cm⁻² s⁻¹. His Figure 7 shows the state of the simulation at a set of times. We use the latest time shown (t = 13.75 s) to deduce properties of the evaporation which are plotted as pentagons in Figures 10 and 11. This very low evaporation velocity appears to follow the departing trend attributed above to a low energy flux.

6. DISCUSSION

The foregoing has used simplified numerical experiments to study the process of conduction-driven chromospheric evaporation and condensation. This led to an analytic model which followed from the simplified scenario of the thermal conduction front leaving the density jump of the initial TR at a uniform temperature. The analytic solution of this isothermal Riemann problem, depicted in Figure 5, predicts an ES whose Mach number, and thus density jump, depends only on the temperature jump of the pre-shock TR, R_tr—and depends on this only logarithmically. The implication of this model is that properties of the evaporation and condensation depend on very few aspects of the flaring loop. While this analytic model departs in some respects from actual solutions, it still turns out to be useful in inferring the scalings of evaporation and condensation. A set of runs are used to develop these scaling relations with more accuracy. These semi-empirical relations are given as Equations (23), (25), (28), and (29).

In order to facilitate the exploration of parameter space we adopted a simplified model for the pre-flare TR. This was motivated by our presupposition that the TR appeared to the flare as a simple density jump. It was shown in Section 4.2 that the specifics of this simplified TR have virtually no effect on the flows. The scalings were then applied, in Section 5, to simulations with more realistic TRs done in this work, and in previous investigations from the literature. The reasonable level of conformance of these cases to the scaling laws supports our presupposition about the role of the TR. One source of complication is that density increases with depth in a realistic chromosphere. In our comparisons to realistic cases we used the local density ahead of the condensation front. Following this prescription, the condensation velocity, given by Equation (29), will decrease with depth in a chromosphere of increasing density.

In a further effort to focus our analysis, we simplified the dynamical equations by excluding aspects deemed inessential to the evaporation. The momentum equation included viscosity, often reduced below the classical value, but omitted gravity. The energy equation included classical thermal conductivity, but omitted radiative losses (at least from most runs). While radiative losses are critical in maintaining the equilibrium structure of the TR, we felt they made less critical contributions to the very rapid dynamics of evaporation. We tested this hypothesis by performing two runs, reported in Section 5.1, with optically thin radiative losses. The additional losses reduced both the peak flare temperature and ES Mach number in comparison to cases with no radiation at all. The fact that losses reduced these by less than 25% suggests that radiation was indeed too slow to significantly affect the energy budget of rapid evaporation.

Unlike the thermal conductive cases here considered, a non-thermal electron beam can deposit energy directly into layers of much higher density and lower temperature where radiative timescales will be shorter. These losses have been shown to make significant, even dominant, contributions to the evaporation dynamics (Fisher 1989). At those temperatures and densities, radiative transfer cannot be treated as optically thin, so a much more sophisticated treatment is required (McClymont & Canfield 1983; Allred et al. 2005). The most basic energetic effect of such radiative transfer is to permit losses, but prevent them from being instantaneous, as they would be under an optically thin treatment. This admittedly simplistic argument suggests that sophisticated radiative treatments of conductively driven flares should produce evaporation falling somewhere between the cases of no losses and of optically thin losses, that we report.

The relations we report will provide the means of incorporating the effects of chromospheric evaporation into simulations focused on the coronal aspects of a solar flare, including magnetic reconnection. To develop the relations we adopted the common measure of energizing the loop with an ad hoc heat source. We expect, however, that our results will be applicable to cases where the plasma is energized self-consistently, for example, by reconnection-generated shocks. Future work is aimed at verifying this as well as investigating the effect that evaporation has on the reconnection and its shocks.

The author thanks Sean Brannon, George Fisher, and Roger Scott for helpful comments and discussion, and the anonymous referee for comments which improved the manuscript. This work was funded by NASA grant NNX13AG09G.

APPENDIX: NEGLECTED EFFECT OF IONIZATION

The foregoing shows how energy input from an ad hoc source, standing in for flare reconnection, heats and accelerates plasma from the chromosphere. In order to simplify the model, complete ionization was assumed, even of the chromospheric plasma. This assumption leads to an over-estimate of the final temperature achieved, T_*, since the energy required to achieve full ionization would reduce the amount available for heating and acceleration. This under-estimate is, however, relatively small and neglecting it leads to a relatively small error.

To estimate the magnitude of the error we consider only that portion of the specific energy, $\varepsilon _*={c}^{\rm (fi)}_{_{\rm V}}T_*$ required by the condensation shock (CS) to heat the material from its initial chromospheric state. Here ${c}^{\rm (fi)}_{_{\rm V}}=(3/2)k_{_{\rm B}}/\bar{m}$ is the specific heat when fully ionized, and the initial temperature is set to zero since T_* ≫ T_{ch, 0}. Had that same specific energy, ε_*, been devoted to both heating and ionizing the plasma, the final temperature achieved would be lower by $\Delta T_{\rm ion}=\varepsilon _{\rm ion}/{c}^{\rm (fi)}_{_{\rm V}}$ , where ε_ion is the energy required to fully ionize a unit mass of chromospheric plasma. The energy required to heat the extra electrons is already accounted for by using ${c}^{\rm (fi)}_{_{\rm V}}$ for the specific heat throughout the ionization process.

We estimate the ionization energy using the fairly conservative approach of considering only hydrogen and helium, whose mass fractions we denote X = 0.74 and Y = 0.25, respectively, but counting the energy required to ionize both from purely neutral ground states. The ionization energy under this scenario is

$\begin{equation} \varepsilon _{\rm ion}= X{\chi _{_{\rm H}}\over m_p} + Y{\chi _{_{\rm He}} \over 4m_p} = {15.1\,{\rm eV}\over m_p}, \end{equation} \tag{ A1 }$

where the energy to fully ionize an atom from its neutral ground state is $\chi _{_{\rm H}} =13.6$ eV for hydrogen and $\chi _{_{\rm He}}=24.6+54.4=79$ eV for helium. The energy removed by ionization thus lowers the final temperature of the CS by

$\begin{equation} \Delta T_{\rm ion}= {\varepsilon _{\rm ion}\over {c}^{\rm (fi)}_{_{\rm V}}} = {2\over 3}{\bar{m}\over m_p}\,{15.1\,{\rm eV}\over k_{_{\rm B}}} = 69{,}000\,{\rm K}. \end{equation} \tag{ A2 }$

Since the final temperature of the CS is well over 2 × 10⁶ K in the simplified runs, the error from neglecting ionization is less than 5%.

To test the foregoing assertion we run a numerical experiment with a specific heat designed to remove energy similar to ionization. The simulation solves Equations (1)–(3) as before except with specific heat

$\begin{equation} {c}_{_{\rm V}}(T) = {c}^{\rm (fi)}_{_{\rm V}}+ \varepsilon _{\rm ion}\, f(T) = {c}^{\rm (fi)}_{_{\rm V}}\,[ 1 + \Delta T_{\rm ion}\, f(T) ], \end{equation} \tag{ A3 }$

with ΔT_ion given by Equation (A2). The function f(T) expresses the distribution of ionization temperatures and is defined to integrate to unity. It should peak at temperatures where each species is being predominately ionized. Motivated by the description above, however, we use a simple version

$\begin{equation} f(T) = {1\over 2\sigma _T}\,\hbox{sech}^2\left({T-T_0\over \sigma _T}\right), \end{equation} \tag{ A4 }$

with a single peak of width σ_T = 80, 000 K centered at T₀ = 250, 000 K. This removes all energy attributable to ionization, i.e., ε_ion, as the plasma is raised from its chromospheric temperature, T_{ch, 0} = 20, 000 K to the CS temperature >10⁶ K. The function is designed to do this gradually enough to be easily resolved by the simulation.

Aside from replacing Equation (5), where ${c}_{_{\rm V}}={c}^{\rm (fi)}_{_{\rm V}}$ , with the modified form in Equation (A3), the simulation is exactly the one performed for Figure 4, with L = 53 Mm, F = 3.5 × 10¹⁰. The result, plotted along the bottom row (d) of Figure 12 is virtually indistinguishable from the case with ${c}_{_{\rm V}}={c}^{\rm (fi)}_{_{\rm V}}$ . To produce any noticeable effect it is necessary to artificially raise the ionization energy, ε_ion, by a factor of ten. Doing so raises the density in the immediate vicinity of the CS, but otherwise leaves the solution unchanged. The conclusion is that ionization, had it been legitimately included, would not have had an appreciable effect on the evaporation or condensation flows we have studied here.

A SIMPLE MODEL OF CHROMOSPHERIC EVAPORATION AND CONDENSATION DRIVEN CONDUCTIVELY IN A SOLAR FLARE

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. THE MODEL