A UNIFIED COMPUTATIONAL MODEL FOR SOLAR AND STELLAR FLARES

Stellar flares are the result of large explosions in the atmospheres of stars. They are produced when magnetic fields, which have been stressed by convective motion in the stellar photosphere, reconnect rapidly, releasing their stored energy. In addition to directly heating the reconnection site in the corona, much of this energy goes into accelerating charged particles to high velocity. These travel down magnetic field lines, colliding with the increasingly dense plasma, depositing their energy and momentum, and quickly heating the plasma in the reconnecting flux tubes to temperatures $\gt 20$ MK. The high temperature causes emission to dramatically brighten in virtually all regions of the electromagnetic spectrum. Flares are ubiquitous and have been observed on stars of nearly all spectral types.

Because of the Sun's close proximity, understanding solar flares is especially important. Together with coronal mass ejections, these directly affect the Earth's environment. They have significant impact on spaced-based communications, the power grid, and the manned space program. Also because of the Sun's proximity, it is possible to spatially resolve many of the features that drive solar flares. Such observations are not possible on other stars. However, despite the lack of spatial resolution, there is still much to be learned by studying stellar flares. (In this paper, we will refer to solar flares as those exclusively on the Sun. Stellar flares will refer to flares on all stars except the Sun.) Active dMe stars are known to have flares much larger than those on the Sun. These stars spin faster and generate much larger magnetic fields (Hawley et al. 2014, and references therein) resulting in flares some 10³ times more energetic than typical solar flares (e.g., Hawley & Pettersen 1991; Hawley et al. 2003; Kowalski et al. 2010). In addition, the background intensity of dMe stars is much smaller than that of the Sun. Thus, when these stars flare, the signal can be much higher. For example, Kowalski et al. (2010) observed the optical luminosity of the dMe star, YZ CMi, to brighten by ∼200 times. In comparison, in the largest solar flare an increase in irradiance of just ∼100 ppm was observed (Woods et al. 2006; Moore et al. 2014).

Accelerated electrons are known to play an important role in transporting energy during flares. Their presence can be detected from the bremsstrahlung radiation they produce as they collide with the ambient plasma. Since these electrons are a major source of flare heating, it is crucial that we accurately model them. Fortunately, RHESSI (Lin et al. 2002) observes this bremsstrahlung radiation from which it is possible to deduce the injected electron spectrum (e.g., Holman et al. 2003). Thus, to simulate how the lower atmosphere responds to this heating, we model the precipitation of these electrons from the acceleration site in the corona to the footpoints in the chromosphere and below.

In addition to electrons, ions are also likely accelerated by reconnecting magnetic fields. Emslie et al. (2012) estimated the energy in accelerated ions to be comparable to that of accelerated electrons in many flares. However, since the bremsstrahlung cross-section is inversely proportional to the square of the mass of the colliding particle (Haug 1997), their presence is much harder to directly detect. Even with this limitation, RHESSI has detected the presence of accelerated ions in several large flares (e.g., Hurford et al. 2006; Emslie et al. 2012). Even though direct evidence of their presence is scarce, models of particle acceleration predict that lower energy ions ($\lt 1\;\mathrm{MeV}$) will be accelerated—albeit on longer timescales than electrons (Petrosian & Liu 2004). Therefore, to accurately model the flaring atmosphere the effects of flare-accelerated ions must also be included.

A beam of charged particles propagating down a magnetic tube induces an electric field which drives a return current (van den Oord 1990; Zharkova & Gordovskyy 2006, and references therein). The magnitude of the current depends upon the particle beam spectrum. In fact, the return current alters the particle spectrum causing a flattening at low energies (Holman 2012). This current additionally heats the ambient plasma through Joule heating. In the flaring corona, this can be a major source of energy and could explain superhot coronal temperatures ($\gt 30$ MK) observed in several large flares (e.g., Caspi & Lin 2010; Caspi et al. 2014).

It is the purpose of this paper to present a unified computational model that can describe the atmosphere of both solar and stellar flares including the processes most important to flare dynamics. To do this we simulate the transport of a beam of non-thermal particles injected at the top of a magnetic flux tube and follow the subsequent heating of the stellar atmosphere. Our model flux tube extends from the sub-photosphere into the corona. Pressure, temperature, and density vary by many orders of magnitude across these regions, and our models must be able to accurately represent these very different conditions. For example, in the corona, radiative transfer is dominated by numerous high-temperature, non-LTE, optically thin atomic transitions. In the photosphere, however, radiative transfer is optically thick and in LTE. In the chromosphere, neither of those approximations hold and the full radiative transfer equation must be solved in detail for several important atomic transitions. Flares produce high-speed shocks (i.e., $\gt 600$ km s⁻¹). We have found that to resolve these requires a computational grid with spacing $\lt 100$ m. Currently, models that include detailed radiative transfer at such high resolution are only tractable in one-dimension.

Decades of observations have revealed that flares certainly have a three-dimensional geometry. However, the strong magnetic forces present in flaring active regions confine charged particles to flow along field lines. Thus, to a good approximation flare dynamics can be modeled using a one-dimensional geometry with that dimension being the axis of a model flux tube. To partially account for the three-dimensional nature of particular flares, the emission from individual flux tube models can be combined using timing and spatial information provided by observations (e.g., images from the Atmospheric Imaging Assembly on the Solar Dynamics Observatory (SDO/AIA) and RHESSI). In this way, important processes are resolved in ways that are currently not tractable in fully 3D simulations. However, this method does not account for the interaction between plasma on differing flux tubes. This could affect heating rates since those depend on the plasma density which can be sensitive to mixing of plasma on differing tubes. We speculate that the effect of mixing on plasma density will be small compared to that produced by chromospheric evaporation, which is captured by our 1D simulations. This is because neighboring flux tubes are likely to have been heated by similar fluxes of non-thermal particles, so will have similarly elevated densities. In this paper, we present results from parameter studies varying 1D model flux tubes and injected non-thermal spectra. In subsequent papers, we will combine these 1D models to form a more complete representation of particular flares.

In Section 2, we describe our method for solving the equations of radiation hydrodynamics and the modeling framework that we use to simulate solar and stellar flares. We discuss how thermal conduction and radiative transfer are included in the models. The chromospheric radiative transfer is of particular importance, and we describe our method to model emission from numerous optically thick, non-LTE atomic transitions which dominate that region. We also describe how radiative backwarming from coronally produced X-rays and extreme ultraviolet (XEUV) radiation is included. In Section 3, we present how our models simulate the precipitation of flare-accelerated particle beams. We present our method for modeling the direct collisional excitation and ionization of the ambient plasma by these particles in Section 4. In the region of beam impact, these dominate over the thermal rates, significantly altering the radiative transfer. Section 5 describes our method for simulating how return currents additionally heat flaring flux tubes. In Section 6, we present results of a parameter study which we have conducted to understand the range in dynamics predicted from various injected particle beams. In Section 7, we summarize our results and present conclusions.

Flares produce high-speed shocks and increase density and radiation throughout the stellar atmosphere. To model them, the radiative transport equation must be coupled with the standard equations of hydrodynamics. During flares, chromospheric plasma evaporates into the corona and waves can propagate into the photosphere. Modeling the atmospheric response to flare heating requires a model that can extend from the sub-photosphere through the chromosphere, transition region, and into the corona. The transition region is extremely narrow, and we have found that accurately resolving it requires spatial scales as small as ∼100 m. A model that includes all of these elements represents a major computational undertaking. We use the RADYN code developed by Carlsson & Stein (1992, 1995, 1997) to solve the equations of charge and population conservation coupled to the equations of radiation hydrodynamics in 1D. We briefly summarize the method of solution here. The equations of radiation hydrodynamics are,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{\partial \rho v}{\partial z}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho v}{\partial t}+\displaystyle \frac{\partial \rho {v}^{2}}{\partial z}+\displaystyle \frac{\partial (p+{q}_{{\rm{v}}})}{\partial z}+\rho g-{A}_{\mathrm{beam}}=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho e}{\partial t}+\displaystyle \frac{\partial \rho {ve}}{\partial z}+(p+{q}_{{\rm{v}}})\displaystyle \frac{\partial v}{\partial z}\\ &&\qquad +\displaystyle \frac{\partial ({F}_{{\rm{c}}}+{F}_{{\rm{r}}})}{\partial z}-{Q}_{\mathrm{cor}}-{Q}_{\mathrm{beam}}-{Q}_{\mathrm{rc}}=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {n}_{i}}{\partial t}+\displaystyle \frac{\partial {n}_{i}v}{\partial z}-(\displaystyle \sum _{j\ne i}^{{N}^{\prime }}{n}_{j}{P}_{{ji}}-{n}_{i}\displaystyle \sum _{j\ne i}^{{N}^{\prime }}{P}_{{ij}})=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mu \displaystyle \frac{\partial {I}_{\nu \mu }}{\partial z}={\eta }_{\nu \mu }-{\chi }_{\nu \mu }{I}_{\nu \mu }\end{eqnarray} \tag{ 5 }$$

where z, ρ, e, v, and p are the height, density, internal energy density, velocity, and pressure, respectively. g is the gravitational acceleration, and q_v is a viscous stress term added to aid in achieving numerical stability. F_c and F_r are the conductive and radiative fluxes, respectively. The conductive flux has the classical Spitzer form but is limited so that it does not exceed the saturation limit of Smith & Auer (1980). In the atomic level population equation, n_i is the number density in a given atomic state, and N' is the total number of atomic states considered in these simulations. We have derived a parameter, m₀, which represents the average mass of the plasma per hydrogen atom. This parameter assumes constant abundance ratios using the abundances derived by Asplund et al. (2009). It has a value of $2.26\times {10}^{-24}$ g. Thus, $\rho ={m}_{0}{n}_{{\rm{H}}}$, where n_H is the hydrogen number density. P_ij is the transition rate from state i to state j and is given by ${P}_{{ij}}={C}_{{ij}}+{R}_{{ij}}$, where C_ij and R_ij are the collisional and radiative rates, respectively. These are discussed in much more detail in Mihalas (1978, Chapter 5). In the radiative transfer equation, ${I}_{\nu \mu }$ is the frequency (ν) and angle (μ) dependent specific intensity, and ${\eta }_{\nu \mu }$ and ${\chi }_{\nu \mu }$ are the emission and absorption coefficients, respectively. The radiative flux divergence, $\partial {F}_{{\rm{r}}}/\partial z$, is obtained by integrating Equation (5) over frequency and angle. Q_cor is a coronal heating term which is necessary to maintain a hot corona. Q_beam and A_beam are terms describing flare-accelerated particle beam heating and momentum deposition, respectively. They are described in more detail in Section 3. Q_rc is a heating term due to return currents and is described in Section 5.

These coupled nonlinear equations are solved implicitly using a Newton–Raphson iteration scheme. Advected quantities are treated using the second-order upwind technique of van Leer (1977). RADYN employs an adaptive grid (Dorfi & Drury 1987) which is designed to resolve shocks and steep gradients that can form in flaring stellar atmospheres. The grid cell concentration is chosen to be proportional to the desired resolution. Dorfi & Drury (1987) define a resolution operator which is strongly dependent on the absolute value of the gradients of the radiative hydrodynamic variables. Because gradients of these variables can be large at different heights, weighting parameters have been chosen to give preference to those variables that require the most grid sensitivity. In these flare simulations, the largest weights were required on temperature, velocity, and atomic level populations to properly resolve the transition region, strong shocks that form in the loops, and non-LTE population densities that affect the convergence of the radiative transfer equation. Due to the complexity of the physical problem, we parameterize certain terms in our model using previous work. Thus, some uncertainties in our results may be due only to the uncertainties in these external models.

2.1. Optically Thick Radiative Transfer

The RADYN code solves the radiative transfer equation for the non-LTE conditions that dominate in the chromosphere. RADYN solves the atomic level population equations (Equation (4)) for a six-level with continuum hydrogen atom, a nine-level with continuum helium atom, a six-level with continuum Ca ii ion, and a four-level with continuum Mg ii ion. This allows the calculation of numerous transitions that are important to the chromospheric energy balance. Table 1 lists the line transitions and Table 2 lists the continuum transitions that we model in detail. For these transitions, the radiative transfer equation is solved for up to 100 frequency points and 5 angular points providing us with detailed line profiles.

Table 1.  Bound–Bound Transitions

Atom	${\lambda }_{{ij}}$ (Å)a	Transition	Atom	${\lambda }_{{ij}}$ (Å)	Transition
H i	1215.67	Lyα	Ca ii	8662.14	$3d{\;}^{2}{D}_{3/2}$ $\leftrightarrow$ $4p{\;}^{2}{P}_{1/2}^{o}$
	1025.73	Lyβ	⋯	8498.02	$3d{\;}^{2}{D}_{3/2}$ $\leftrightarrow$ $4p{\;}^{2}{P}_{3/2}^{o}$
	972.52	Lyγ	⋯	8542.09	$3d{\;}^{2}{D}_{5/2}$ $\leftrightarrow$ $4p{\;}^{2}{P}_{3/2}^{o}$
	949.74	Lyδ	He i	625.56	$1{s}^{2}$ ${}^{1}{S}_{0}$ $\leftrightarrow$ $1s$ $2s$ ${}^{3}{S}_{1}$
	6562.79	Hα	⋯	601.42	$1{s}^{2}$ ${}^{1}{S}_{0}$ $\leftrightarrow$ $1s$ $2s$ ${}^{1}{S}_{0}$
	4861.35	Hβ	⋯	10830.29	$1s$ $2s$ ${}^{3}{S}_{1}$ $\leftrightarrow$ $1s$ $2p$ ${}^{3}{P}_{4}^{0}$
	4340.47	Hγ	⋯	584.33	$1{s}^{2}$ ${}^{1}{S}_{0}$ $\leftrightarrow$ $1s$ $2p$ ${}^{1}{P}_{1}^{0}$
	18751.3	Paα	⋯	20581.29	$1s$ $2s$ ${}^{1}{S}_{0}$ $\leftrightarrow$ $1s$ $2p$ ${}^{1}{P}_{1}^{0}$
	12818.1	Paβ	He ii	303.79	$1s$ ${}^{2}{S}_{1/2}$ $\leftrightarrow$ $1s$ $2p$ ${}^{2}{S}_{1/2}$
	40522.8	Brα	⋯	303.78	$1s$ ${}^{2}{S}_{1/2}$ $\leftrightarrow$ $2p$ ${}^{2}{P}_{2}^{0}$
Ca ii	3968.47	H	Mg ii	2802.70	h
	3933.66	K	⋯	2795.53	k

Note.

^aThese are vacuum (air) wavelengths for ${\lambda }_{{ij}}$ below (above) 2000 Å.

Download table as:  ASCIITypeset image

Table 2.  Bound–Free Transitions

Atom	${\lambda }_{{ic}}$ (Å)	Initial State	Atom	${\lambda }_{{ic}}$ (Å)	Initial State
H i	911	n = 1	He i	504	$1{s}^{2}$ ${}^{1}{S}_{0}$
	3646	n = 2	⋯	2600	$1s$ $2s$ ${}^{3}{S}_{1}$
	8204	n = 3	⋯	3121	$1s$ $2s$ ${}^{1}{S}_{0}$
	14584	n = 4	⋯	3421	$1s$ $2p$ ${}^{3}{P}_{4}^{0}$
	22787	n = 5	⋯	3679	$1s$ $2p$ ${}^{1}{P}_{1}^{0}$
Ca ii	1044	$4s{\;}^{2}{S}_{1/2}$	He ii	228	$1s$ ${}^{2}{S}_{1/2}$
	1218	$3d{\;}^{2}{D}_{3/2}$	⋯	911	$2s$ ${}^{2}{S}_{1/2}$
	1219	$3d{\;}^{2}{D}_{5/2}$	⋯	911	$2p$ ${}^{2}{P}_{1/2}^{0}$
	1417	$4p{\;}^{2}{P}_{1/2}^{o}$	Mg ii	824	$2{p}^{6}$ $3s$ ${}^{2}{S}_{1/2}$
	1422	$4p{\;}^{2}{P}_{3/2}^{o}$	⋯	1168	$2{p}^{6}$ $4p$ ${}^{2}{P}_{1/2}^{0}$
	⋯	⋯	⋯	1169	$2{p}^{6}$ $4p$ ${}^{2}{P}_{3/2}^{0}$

Download table as:  ASCIITypeset image

2.2. Optically Thin Radiative Transfer

The densities in the transition region and corona are typically low enough that the "coronal approximation" applies. In this case, the radiative transfer is dominated by numerous optically thin lines. These lines are formed when ions are collisionally excited and radiatively de-excited. Since the atmosphere is optically thin in these regions, this radiation is assumed to escape, thus having a net cooling effect on the corona. We model the radiative transfer in these regions using a radiative loss function which is obtained by summing all transitions in the CHIANTI database (Dere et al. 1997; Landi et al. 2013) except those already accounted for in Tables 1 and 2. To generate this function the CHIANTI calculations were performed with the assumption of a constant electron density of 10¹⁰ cm⁻³. Figure 1 shows this radiative loss function.

Zoom In Zoom Out Reset image size

2.3. XEUV Backwarming

Half of the optically thin radiative losses described above are directed outward and leave the stellar atmosphere. The other half are directed downward and will get absorbed in deeper and denser regions. This additional source of heating and ionization in the lower atmosphere becomes especially important during flares when the coronal X-ray and extreme ultraviolet (XEUV) emission can become elevated by orders of magnitude. Previously, Allred et al. (2005, hereafter, A05) developed a method to model how heat is deposited from XEUV backwarming, but that method did not account for the increased photoionizations from this flux. Here we present a new technique that self-consistently includes heating and photoionizations. We have used the CHIANTI database to tabulate emissivities for numerous transitions as a function of temperature and wavelength. RADYN calculates the XEUV spectrum produced from a model loop by integrating the product of these emissivities with the emission measures from the transition region and coronal portions of the loop. Finally, the radiative transfer equation is solved for several optically thick chromospheric lines, as described above, assuming that this emission is incident on the chromosphere from above. In solving the radiative transfer equation, photoionization cross-sections for XEUV emission are calculated and added to the rates as described by Wahlstrom & Carlsson (1994).

To understand how XEUV backwarming affects our loop models, we have generated loops which include and exclude this heating term. Figure 2 compares the structure of the QS.SL.HT loop model (see Section 2.7) which has been generated with and without XEUV backwarming and using the technique of A05. We find that the XEUV backwarming term results in a chromosphere that is 1000–2000 K hotter than would otherwise be expected. The technique presented here and that of A05 produce similar temperature structures. However, the radiative transfer predicted by these methods is significantly different. This is illustrated in Figure 3 which compares the emission from the Ca ii H and He i 10830 Å lines for loops generated with and without XEUV backwarming and using the technique of A05. Our technique produces a He i 10830 Å line with a much deeper absorption profile. This line is formed when continuum photons from the photosphere are absorbed by neutral helium atoms in the $2s$ ${}^{3}{S}_{1}$ excited state. The XEUV flux increases the photoionization rate of helium. These ions quickly recombine into the $2s$ ${}^{3}{S}_{1}$ state (Golding et al. 2014), resulting in more absorption than would be expected without the XEUV flux. The Ca ii H line can be understood similarly. Without the inclusion of the XEUV flux, the region of the chromosphere where the Ca ii H line center forms is cooler, resulting in an overabundance of ground-state ions and an overall absorption profile. For both XEUV backwarming techniques, the Ca ii H line has a central reversal peaking in the near wings. The technique of A05 predicts more flux in the central reversal than our technique. Since their technique does not include direct photoionizations by the XEUV flux, it predicts an overabundance of Ca ii ions relative to Ca iii ions, and hence an increased flux in the Ca ii H central reversal.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.4. Opacity

2.5. Line Broadening

Observations of flares on dMe stars (Hawley & Pettersen 1991; Hawley et al. 2003; Kowalski et al. 2010) and the Sun (Johns-Krull et al. 1997) have shown that the hydrogen Balmer lines can become extremely broadened. These authors speculate the cause to be Stark broadening, and that conclusion is supported by the models of Allred et al. (2006) and Paulson et al. (2006). In order to test this against other possible sources of broadening, such as thermal or turbulent broadening, we have implemented a technique in RADYN to model Stark line broadening. The amount of Stark broadening depends upon the local electron density. Therefore, determining it can further constrain models of flaring stellar atmospheres.

That the orbital angular momentum states of hydrogen labeled by the quantum number, l, are not degenerate in the presence of an electric field perturbation is described by the well-known Stark effect (e.g., Kepple & Griem 1968; Vidal et al. 1971; Seaton 1990). In stellar atmospheres, this perturbation is due to a net electric microfield from fluctuations in the ambient electron, proton, and ion density. The first-order (linear) energy level shifts are proportional to the electric microfield strength (which has a probability distribution typically modeled as a Holtsmark or Hooper distribution; Nayfonov et al. 1999), the principal quantum number, n, and a quantum number, q, which describes the relationship between the parabolic quantum numbers q₁ and q₂ where $q={q}_{1}-{q}_{2}$ and ${q}_{1}+{q}_{2}=n-| m| -1$ (Condon & Shortley 1935). Therefore, the higher order hydrogen lines within a series experience a larger amount of Stark broadening.

The general prescription for line cross-section is given by a Voigt profile with a damping parameter, Γ, that includes separate terms for radiative damping, ${{\rm{\Gamma }}}_{\mathrm{rad}}$, and for collisions among neutrals (resonance and van der Waals), ${{\rm{\Gamma }}}_{\mathrm{res}/\mathrm{vdW}}$. A microturbulence of 2 km s⁻¹ is included in the Doppler width. The Voigt profile is then convolved with the Stark profile to give the total line cross-section (Mihalas 1978). However, as an approximation to the more complete treatment, the Voigt profile damping parameter, Γ, can be modified to also include a Stark damping term, ${{\rm{\Gamma }}}_{{\rm{s}}}$, such that ${\rm{\Gamma }}={{\rm{\Gamma }}}_{\mathrm{rad}}+{{\rm{\Gamma }}}_{\mathrm{res}/\mathrm{vdW}}+{{\rm{\Gamma }}}_{{\rm{s}}}$. This has been found to work well for solar hydrogen lines in the infrared (Carlsson & Rutten 1992), so we have chosen this method for modeling Stark broadening in RADYN.

The Stark damping parameter for hydrogen is taken from the approximate line shape formulae in Method #1 of Sutton (1978). In using this method, we are consistent with other NLTE radiative transfer codes such as RH (Uitenbroek 2001). The Stark damping profile is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{s}}}=0.6{g}_{q}({j}^{2}-{i}^{2}){n}_{e}^{2/3}\end{eqnarray} \tag{ 6 }$$

where ${g}_{q}=0.642$ for the α transition of each series and ${g}_{q}=1$ otherwise, and j and i are the upper and lower levels of the transition, respectively. For transitions of helium, calcium, and magnesium, ${{\rm{\Gamma }}}_{{\rm{s}}}={q}_{{\rm{s}}}{n}_{e}$, where g_s is a constant resulting in Stark broadening which is approximately three orders of magnitude less than for hydrogen.

2.6. Boundary Conditions

2.7. Initial Loops

A focus of this work is to study how particle beams, which are known to be accelerated during flares, deposit their energy in magnetic flux loops in stellar and solar atmospheres. To this end, we have generated several initial loop states with diverse lengths, temperatures, and densities which we will use in this study. These extend from the sub-photosphere through the corona. The corona is kept hot by adding a heating term, Q_cor, which is chosen to just balance the conductive and radiative losses in the upper coronal portion of the loop. Heat is also added to the sub-photosphere to balance radiative and convective flux losses there. With these heat sources, the loops are allowed to relax until a state of near equilibrium is reached.

We have generated solar-type initial loops using three free parameters. These are the photospheric temperature, loop-length, and coronal temperature. The coronal density is dependent on loop-length and coronal temperature by an RTV-type (Rosner et al. 1978) scaling law so cannot be independently varied. We have chosen photospheric temperatures (i.e., the temperature in our loops where ${\tau }_{5000}=1$) of 5000 and 5800 K which correspond with sunspot and "quiet Sun" conditions. Loop lengths were chosen to be 10 Mm for short loops and 100 Mm for long loops. So that we can model flaring flux tubes in dMe stellar atmospheres, we have also produced a model loop appropriate for M dwarf atmospheres. This loop has a higher surface gravity (562 m s⁻²), cooler photosphere (3500 K), and hotter corona (6 MK) than the solar loops. The parameters describing these loops are summarized in Table 3 and the temperature and density of the loops are plotted in Figure 4.

Zoom In Zoom Out Reset image size

Table 3.  Initial Loops

Label	Photospheric Temperature (K)	Loop-length (Mm)	Coronal Electron Density (cm−3)	Coronal Temperature (MK)
QS.LL.HT	5800	100	$3\times {10}^{8}$	3.0
QS.LL.LT	5800	100	$2\times {10}^{7}$	1.0
QS.SL.HT	5800	10	$6\times {10}^{9}$	3.0
QS.SL.LT	5800	10	$6\times {10}^{8}$	1.0
SS.LL.HT	5000	100	$3\times {10}^{8}$	3.0
SS.LL.LT	5000	100	$2\times {10}^{7}$	1.0
SS.SL.HT	5000	10	$6\times {10}^{9}$	3.0
SS.SL.LT	5000	10	$6\times {10}^{8}$	1.0
M dwarf	3500	10	$2\times {10}^{10}$	6.0

Download table as:  ASCIITypeset image

As described in Section 3, the energy deposition rate from particle beams depends upon the magnetic field in the loops. To account for this we assume a magnetic field strength in the photosphere of 1 kG for solar loops and 5 kG for the M dwarf loop. Spatially averaged observations of magnetic fields on active M dwarf stars have found that the majority of the stellar surface (60%–70%) has a magnetic field strength of $3-4$ kG (Saar & Linsky 1985; Saar 1994; Johns-Krull & Valenti 1996). Since these are spatial averages and large flares are likely to occur where the field is strongest, we have estimated a photospheric field of 5 kG in the loop. At the loop tops, the magnetic field strength is assumed to be 100 G for the shorter (10 Mm) solar loops, 10 G for the longer (100 Mm) solar loops, and 500 G for the M dwarf loop. These values were chosen to be consistent with active region coronal magnetic field extrapolations (e.g., Tadesse et al. 2012). The magnetic field in all cases is assumed to exponentially decrease between these boundary conditions.

Accelerated particle beams are the main source of heating in the lower atmosphere during flares, so it is crucial that we accurately model their propagation and energy deposition as they move through magnetic flux loops. The rate of energy lost by a particle beam depends upon the details of the beam's energy spectrum as well as the composition and ionization states of the ambient plasma. We use the Fokker–Planck method coupled with results from our radiative hydrodynamic simulation described in Section 2 to determine the distribution function, $f({\boldsymbol{x}},{\boldsymbol{p}},t)$, for flare-accelerated particles. Because the magnetic field constrains the beam particles to move along field lines, we can replace the vector position, ${\boldsymbol{x}}$, with the coordinate, z, which measures the distance along the field line. We replace the vector momentum, ${\boldsymbol{p}}$, with the kinetic energy, E, and particle pitch angle, α. In the following discussion, we will simply refer to $\mu =\mathrm{cos}(\alpha )$ as the pitch angle. The transit time for beam particles from injection source to footpoints is fast relative to changes in the beam particle energy distribution and hydrodynamic timescales, so a time-independent distribution function is assumed. Thus, we will solve the kinetic equation for the distribution, $f(E,\mu ,z)$, where $f(E,\mu ,z){dEd}\mu {dz}$ is the number density of beam particles with energy between E and ${dE}$, pitch angle between μ and $d\mu$, and position between z and ${dz}$. In solving the kinetic equation, we include the effects of particles moving at relativistic speeds, Coulomb collisions, synchrotron emission, pitch angle scattering, and magnetic field gradients. However, we neglect effects due to plasma turbulence, electric fields, and self-absorption of synchrotron emission. We follow McTiernan & Petrosian (1990) and write the Fokker–Planck equation as

$$\begin{eqnarray}&&\mu \displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}-\displaystyle \frac{d\mathrm{ln}B}{2{dz}}\displaystyle \frac{\partial }{\partial \mu }[(1-{\mu }^{2}){\rm{\Phi }}]\\ &&\quad =\displaystyle \frac{1}{{\beta }^{2}}\displaystyle \frac{\partial }{\partial E}\{[C+S{\beta }^{3}{\gamma }^{2}(1-{\mu }^{2})]{\rm{\Phi }}\}\\ &&\qquad -\displaystyle \frac{S}{\beta \gamma }\displaystyle \frac{\partial }{\partial \mu }[\mu (1-{\mu }^{2}){\rm{\Phi }}]\\ &&\qquad +\displaystyle \frac{C^{\prime} }{{\beta }^{4}{\gamma }^{2}}\displaystyle \frac{\partial }{\partial \mu }[(1-{\mu }^{2})\displaystyle \frac{\partial {\rm{\Phi }}}{\partial \mu }]+\displaystyle \frac{{\rm{\Sigma }}}{c{\beta }^{2}}\end{eqnarray} \tag{ 7 }$$

where ${\rm{\Phi }}=f/\beta$, $\gamma =E+1$ is the relativistic total energy with E measured in units of ${{mc}}^{2}$, $\beta c$ is the relativistic ion velocity, v, B is the magnetic field strength in the loop, and Σ is a source term for particles injected at the loop top. C and C' are coefficients which measure the beam energy loss rate and pitch angle diffusion from collisions, respectively. Similarly, S measures the loss rate and pitch angle diffusion due to synchrotron emission. McTiernan & Petrosian (1990) developed a numerical method to solve this equation for flare-accelerated electrons. We have generalized their method to solve the Fokker–Planck equation for both flare-accelerated electrons and ions. This has required altering C, C', and S to account for the physics which characterize ion beams. Here we summarize this theory.

Petrosian (1985) calculated the energy loss and pitch angle scattering rates due to synchrotron radiation. From those results, we find the synchrotron coefficient, S, for both ion and electron beams to be

$$\begin{eqnarray}&&S=\displaystyle \frac{2}{3}\displaystyle \frac{{({Ze})}^{4}{B}^{2}}{{({{mc}}^{2})}^{3}}\end{eqnarray} \tag{ 8 }$$

where e is the proton charge, and Z is the number of protons (or electrons) in the beam particle.

C is related to the energy loss due to collisions by

$$\begin{eqnarray}&&C=-\displaystyle \frac{{{dE}}_{\mathrm{col}}}{{dt}}\displaystyle \frac{\beta }{c}\end{eqnarray} \tag{ 9 }$$

McTiernan & Petrosian (1990) assumed that the ambient plasma was a "cold target" meaning that the beam particle's velocity is much greater than the thermal velocity for all constituents of the ambient plasma. Formally, this is expressed as ${x}_{i}\gg 1$ where ${x}_{i}={(v/{v}_{i})}^{2}=({m}_{i}/m)E/({{kT}}_{i})$ and m_i, v_i, and T_i are the mass, velocity, and temperature for the plasma constituent, i, and m is the mass of the beam particle. For flare-accelerated electrons even at the low-energy limit of ∼10 keV, the ambient electron temperature would have to be greater than 100 MK for the cold target approximation to fail. This far exceeds observed temperatures for even the largest flares. Therefore, beam electrons can always be treated using cold target collision theory. With ion beams, however, this is not necessarily true. Flare-accelerated protons of ∼1 MeV (which is about the low energy detection limit of current instruments) have a velocity less than that of thermal electrons at temperature $\gtrsim 6$ MK. Flares routinely produce plasma much hotter than this. Therefore when modeling how ion beams interact with the ambient plasma, we must employ a theory which can account for both hot and cold targets. Trubnikov (1965) developed this theory and found the energy loss rate of particles with energy, E, from collisions with ambient charged particles of species, i, to be

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{i}}{{dt}}=-E{\nu }_{i}\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\nu }_{i}=2{\nu }_{{0}_{i}}\displaystyle \frac{m}{{m}_{i}}[{Erf}(\sqrt{{x}_{i}})-\displaystyle \frac{2}{\sqrt{\pi }}\sqrt{{x}_{i}}{e}^{-{x}_{i}}(1+\displaystyle \frac{{m}_{i}}{m})]\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}&&{\nu }_{{0}_{i}}=\displaystyle \frac{2\pi {e}^{4}{Z}^{2}{Z}_{i}^{2}{\lambda }_{i}{n}_{i}}{{{Em}}^{2}{c}^{2}v}\end{eqnarray} \tag{ 12 }$$

Z_i is the target particle charge in units of e, and n_i is the number density for target particle, i. ${\lambda }_{i}$ is the Coulomb logarithm, and is defined by ${\lambda }_{i}=\mathrm{ln}({r}_{\mathrm{max}}/{r}_{\mathrm{min}})$. In this case r_max is the mean free path length $\eta =v/\nu$ (Emslie 1978, hereafter, referred to as E78) where ν is the plasma oscillation frequency. r_min is given by ${r}_{\mathrm{min}}={\rm{\hslash }}/2{M}_{i}v$ where M_i is the reduced mass of the beam and target particle (Huba 2011). Combining these gives

$$\begin{eqnarray}&&{\lambda }_{i}=\mathrm{ln}(\displaystyle \frac{{M}_{i}{v}^{2}}{{\rm{\hslash }}}{[\displaystyle \frac{{m}_{i}}{\pi {n}_{i}{Z}_{i}^{2}{e}^{2}}]}^{\displaystyle \frac{1}{2}}).\end{eqnarray} \tag{ 13 }$$

The collisional pitch angle diffusion coefficient, C', is obtained from the rate that beam particles are scattered out of the direction of beam propagation. It is given by

$$\begin{eqnarray}&&C^{\prime} =\displaystyle \frac{{\beta }^{2}}{{c}^{2}}\displaystyle \frac{{dv}}{{dt}}.\end{eqnarray} \tag{ 14 }$$

For collisions with neutral targets (e.g., neutral hydrogen or helium), the beam particles interact with atomically bound electrons. As long as the beam velocity is much greater than the Bohr orbital velocity, neutral atoms can be treated as cold targets. Using the fine structure constant to approximate the orbital velocity, we find that the cold target approximation for collisions with neutral hydrogen is good for protons with energy $\geqslant 25$ keV and electrons $\geqslant 0.01$ keV. The energy loss rate for a charged beam on a neutral target is given by (E78; Mott & Massey 1965)

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{n}}{{dt}}=\displaystyle \frac{-2\pi {e}^{4}}{E}\displaystyle \frac{m}{{m}_{e}}{Z}^{2}{Z}_{n}{n}_{n}{\lambda }_{n}v\end{eqnarray} \tag{ 15 }$$

where n_n is the number density for the ambient neutral atom, n. Z_n is the number of electrons in the atom, and ${\lambda }_{n}$ is an effective Coulomb logarithm for collisions with the atom. For an atom with mean ionization energy, I_n, ${\lambda }_{n}$ is given by (Evans 1955)

$$\begin{eqnarray}&&{\lambda }_{n}=\mathrm{ln}(\displaystyle \frac{2{m}_{e}{v}^{2}{\gamma }^{2}}{{I}_{n}}).\end{eqnarray} \tag{ 16 }$$

The total energy loss rate due to collisions is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{col}}}{{dt}}=\displaystyle \sum _{i}\displaystyle \frac{{{dE}}_{i}}{{dt}}+\displaystyle \sum _{n}\displaystyle \frac{{{dE}}_{n}}{{dt}}\end{eqnarray} \tag{ 17 }$$

where the sums are over charged and neutral species, respectively. In calculating the beam energy deposition rate in RADYN, we consider a plasma consisting of electrons; protons; neutral hydrogen; neutral, singly, and doubly ionized helium; and singly ionized calcium and magnesium. These are the species most important to the energetics of the chromosphere where the majority of the beam impacts. However, since the energy loss rates are inversely proportional to the target particle's mass, collisions with electrons—both ambient and those in neutral atoms—dominate the energy transfer. We use Equation (10) for calculating the energy loss rate for collisions with charged particles and ions, and Equation (15) for collisions with neutral particles. We use Equations (9) and (17) to obtain C, which is needed to solve Equation (7). The number densities for each of these species is calculated in detail in the simulations. Thus, the beam energy deposition drives the simulations which alter the beam deposition in a self-consistent way.

By solving Equation (7), we obtain the distribution function, f, from which we can obtain the particle beam heating rate (Q_beam) and momentum deposition rate (A_beam) using the following relations:

$$\begin{eqnarray}&&{Q}_{\mathrm{beam}}=\displaystyle \frac{d}{{dz}}({\displaystyle \int }_{\mu }{\displaystyle \int }_{E}\mu {vEf}\;{dEd}\mu )\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}&&{A}_{\mathrm{beam}}=\displaystyle \frac{d}{{dz}}({\displaystyle \int }_{\mu }{\displaystyle \int }_{E}\mu {vaf}\;{dEd}\mu )\end{eqnarray} \tag{ 19 }$$

where $a=\gamma {mv}$ is the relativistic momentum. Even though we have developed and presented here a method to simulate both flare-accelerated electrons and ions, in the remainder of this paper we will focus on the effects of electron beams. A detailed study of ion beams will be presented in a future paper.

During flares, the transition rates, P_ij, can be very enhanced due to direct ionizations and excitations caused by collisions of non-thermal beam particles with the ambient plasma. In this section, we describe a method for incorporating these elevated rates in the most important species (i.e., hydrogen and helium) in RADYN.

4.1. Hydrogen

The non-thermal collisional rate from a flux of charged particles is given by the general formula (Mott & Massey 1965, Chapter xvi),

$$\begin{eqnarray}&&{C}_{\mathrm{beam}}^{\mathrm{nt}}=2\pi \int \int \mu {v}_{B}f{\sigma }_{{ij}}{dEd}\mu \end{eqnarray} \tag{ 20 }$$

where f is the beam particle distribution function, μ is the particle pitch angle as described in Section 3, and ${\sigma }_{{ij}}$ is the cross section for ionization or excitation.

Since the flare-accelerated particles have much larger energy than the ionization potential of hydrogen, secondary ionizations from electrons liberated in the primary ionization need to be included. Dalgarno & Griffing (1958) have performed calculations to include these secondary ionizations, giving an average energy of 36 eV per electron–ion pair produced in the complete absorption of beam particles with energy in the range of 200–1000 eV by a cold hydrogen gas. The amount of energy produced per electron–ion pair is constant for $E\gt 200$ eV, and we assume that this extends to all electron energies above 1000 eV. Including the secondary ionizations, the non-thermal collisional ionization rate can be written in terms of the amount of energy lost by the beam to neutral hydrogen, $\frac{{{dE}}_{{\rm{H}}}}{{dt}}$ obtained from Equation (15), (Fang et al. 1993; see also Ricchiazzi & Canfield 1983),

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{{\rm{H}}}}{{dt}}=36\;\mathrm{eV}\times {n}_{1}{C}_{\mathrm{beam},1{\rm{c}}}^{\mathrm{nt}}\end{eqnarray} \tag{ 21 }$$

assuming that most of the neutral atoms are in the ground state, n₁.

A comprehensive study of non-thermal rates in solar flares was presented by Ricchiazzi (1982) and Ricchiazzi & Canfield (1983), who found that the H i rates of non-thermal collisional ionization from the ground state (denoted, 1-c) and collisional excitation from the ground to first excited state (denoted, 1–2) were important compared to their respective thermal rates, but that the non-thermal 2-c rate was negligible compared to the thermal rate for the applied range of beam fluxes. Fang et al. (1993) have presented revised H i 1-c, 1–2, 1–3, 1–4 non-thermal ionization and excitation rates, which were adopted by Kašparová et al. (2009) in their study of the hydrogen lines in response to electron beams. Notably, the 1-c rate is a factor ∼4.6 higher compared to the Ricchiazzi & Canfield (1983) 1-c rate. To facilitate comparison with the most recent studies, we use the constants (${R}^{H,\mathrm{nt}}$) from Fang et al. (1993) to derive the non-thermal collisional rates with hydrogen, according to the formula,

$$\begin{eqnarray}&&{C}_{\mathrm{beam}}^{\mathrm{nt}}={R}^{H,\mathrm{nt}}\displaystyle \frac{{{\rm{\Lambda }}}_{n}}{{n}_{e}{{\rm{\Lambda }}}_{i}+{n}_{{\rm{H}}}{{\rm{\Lambda }}}_{n}}{Q}_{\mathrm{beam}}\end{eqnarray} \tag{ 22 }$$

where n_H is the neutral hydrogen density and ${R}^{H,\mathrm{nt}}=1.73\times {10}^{10}$, $2.94\times {10}^{10}$, $5.35\times {10}^{9}$, $1.91\times {10}^{9}$ for the 1-c, 1–2, 1–3, and 1–4 transitions, respectively.

4.2. Helium

A05 found that the ionization fraction of He ii controls how a flaring flux tube transitions from gentle to explosive chromospheric evaporation. To more accurately model these dynamics, we include the non-thermal rates for the primary ionization of neutral helium (He i 1s² $\to$ He ii 1s) and singly ionized helium (He ii 1s $\to$ He iii). Younger (1981) calculated parameterized cross sections for these transitions as a function of impacting electron energy. For clarity, we reproduce his result here,

$$\begin{eqnarray}{\sigma }_{{ij}} & = & \displaystyle \frac{1}{{{uI}}^{2}}[A(1-\displaystyle \frac{1}{u})+B{(1-\displaystyle \frac{1}{u})}^{2}\\ & & +C\mathrm{ln}(u)+\displaystyle \frac{D}{u}\mathrm{ln}(u)]\end{eqnarray} \tag{ 23 }$$

where I is the ionization energy and $u=E/I$ is the normalized energy of the colliding electron. A, B, C, and D are constants provided by Younger (1981) and Arnaud & Rothenflug (1985) and have values of 17.8, $-11.0$, 7.0, $-23.2$ and 14.4, $-5.6$, 1.9, $-13.3$ in units of ${10}^{-14}$ cm² eV² for He i and He ii, respectively. With these cross sections and the beam particle distribution function, f, obtained from Equation (7), we calculate the non-thermal helium ionization rates using Equation (20).

Images of the footpoint locations of positively charged ion beams and negatively charged electron beams indicate that they are not co-spatial (Hurford et al. 2006). For an electron beam without a co-streaming positively charged beam of equal current density, the charge imbalance between the acceleration region and any given point along the magnetic loop leads to a macroscopic return current electric field (Hoyng et al. 1976; Knight & Sturrock 1977; van den Oord 1990) in a direction along the beam. The return current electric field decelerates the beam electrons while accelerating ambient electrons,⁴ which drift in the direction opposite that of the return current electric field toward the loop top with a Maxwellian distribution of speeds, assuming that the return current field is not super-Dreicer (Holman 1985). These drifting electrons form the return current, which heats the atmosphere through Ohmic (Joule) dissipation. Many aspects of the beam-return current-atmospheric system have been studied, including pitch angle modifications of the beam (Emslie 1980), return current collisional rates (Karlický et al. 2004), return current-beam instabilities and subsequent particle acceleration (Karlický & Kontar 2012; Pechhacker & Tsiklauri 2014), and turbulent effects on the beam-return current system (Xu et al. 2013; Kontar et al. 2014). The return current modifications on the classical thick target model (Brown 1971) have been used to explain the difference between looptop and footpoint hard X-ray spectral indices (Battaglia & Benz 2008; Xu et al. 2013).

Recently, Holman (2012) derived return current heating rates in flaring coronal conditions. After the first short moments of beam propagation, a relatively steady-state is quickly attained (van den Oord 1990) in which the magnitude of the return current density is equal to the beam current density. The return current electric field can be determined from Ohm's law: ${E}_{\mathrm{rc}}=\eta {J}_{\mathrm{beam}}$, where η is the classical Spitzer resistivity and ${J}_{\mathrm{beam}}$ is the return current density. Holman (2012) determined the return-current plasma heating rate for a given electron beam flux self-consistently accounting for Joule heating and the reduction of beam electrons due to (non-collisional) thermalization caused by the return current electric field. However, this treatment does not self-consistently account for the change in flux of beam electrons due to Coulomb collisions with the ambient plasma and is therefore an overestimate, especially in the chromosphere where Coulomb collisions are large. In the corona, however, there are relatively few collisions, and this approximation works quite well. Additionally, return current heating is likely largest in the corona since it is there that beam current densities are highest. Therefore, to account for return current heating in flares, we have chosen to incorporate the formalism of Holman (2012) into RADYN. The result for the volumetric return current heating rate from a power-law electron beam spectrum with power-law index, δ, and cutoff energy, E_c, is given by

$$\begin{eqnarray}{Q}_{\mathrm{rc}}(x)=\{\begin{array}{ll}\eta {e}^{2}{{F}_{{\rm{e}}}}^{2} & x\lt {x}_{\mathrm{rc}}\\ \eta {e}^{2}{{F}_{{\rm{e}}}}^{2}(\delta \displaystyle \frac{{E}_{\mathrm{therm}}}{{E}_{{\rm{c}}}}+\displaystyle \frac{V(x)}{{E}_{{\rm{c}}}}) & \\ \times {(\displaystyle \frac{{E}_{\mathrm{therm}}}{{E}_{{\rm{c}}}}+\displaystyle \frac{V(x)}{{E}_{{\rm{c}}}})}^{1-2\delta } & x\geqslant {x}_{\mathrm{rc}}\end{array}\end{eqnarray} \tag{ 24 }$$

where x is distance from the loop top and F_e is the beam particle flux given by $\int \int \mu {vf}\;{dEd}\mu$. Beam electrons are assumed to be thermalized and removed from the beam when their energy is ${E}_{\mathrm{therm}}=2.5{kT}$. x_rc is the distance at which the lowest energy electrons in the beam are thermalized, and $V(x)$ is the electric potential energy as described in Holman (2012). A more complete treatment which includes return current losses directly in the Fokker–Planck equation (Equation (7)), similar to work done by Zharkova & Gordovskyy (2005) and Dobranskis & Zharkova (2014), is outside the scope of this paper but will be presented in a future paper.

The energy distribution of flare-accelerated particles can vary greatly between particular flares. Since the impact location strongly depends on the particle distribution, that distribution is very important in determining how stellar atmospheres respond to flare heating. To explore the range of possibilities, we have performed a parameter study varying the injected electron beam spectrum, pitch-angle distribution, and the initial loop conditions into which these particles impact.

Flare-accelerated particles are typically observed to have a power-law energy distribution with a cutoff energy, E_c, below which relatively few particles are accelerated. Cutoff energies as high as 250 keV have been observed (R. Schwartz 2015, private communication), but in many flares, the cutoff energy cannot be directly determined since thermal X-ray emission in the range below $10-20$ keV from the flare-heated plasma overwhelms the non-thermal bremsstrahlung emission. In these flares, it is only possible to determine an upper limit for the cutoff energy. To explore how the atmosphere responds to beams of particles with energies below this upper limit, we have chosen a lower limit of 5 keV. Therefore, to account for the full range in parameter space, we have varied the cutoff energy between 5 and 500 keV. The slope of the power law is given by the spectral index, δ. Typical power-law indices range from 3 to 9, but in this study we have varied it between 2.5 and 10 to ensure that we have spanned the full range of possible values. We have also varied the pitch angle distribution with which the flare-accelerated particles were injected. We modeled this distribution as a Gaussian centered around the flux loop axis. Thus, it has the form, ${e}^{{((\mu -1)/{\mu }_{0})}^{2}}$, where ${\mu }_{0}$ is a parameter that controls the width of the Gaussian distribution. In addition to the Gaussian distribution, we have included simulations in which the particles are fully beamed, i.e., they all start with a 0° pitch angle, and simulations with isotropic distributions. We have modeled the effects of beam impact varying these parameters on each of the loops in Table 3. The parameters are summarized in Table 4. To perform this study, we ran more than 11,000 simulations.

Table 4.  Parameter Study

Parameter	Range
Ec	5–500 keV
δ	2.5–10
${\mu }_{0}$	0.1–1.0, isotropic, beamed
Loop Conditions	Those listed in Table 3

Download table as:  ASCIITypeset image

6.1. Penetration Depths in Flux Loops

First, we consider how varying E_c and δ affect the location of beam impact. Figure 5 shows a few representative examples of the heating rate due to collisions with beam particles as a function of hydrogen column density in the QS.SL.HT loop. Clearly the heating extends over a range of column densities. We have determined an average value indicated by the dashed lines. In the following discussion, the term penetration depth (in units of column density) refers to this average. Figure 6 shows the beam penetration depth as a function of cutoff energy and the power-law index in the QS.SL.HT and M dwarf flux tubes. A striking result of this analysis is how strongly dependent penetration depth is on δ. Low values of δ imply a significant number of higher energy particles even for spectra with low cutoff energies. For example, an electron distribution with ${E}_{{\rm{c}}}=10$ keV and $\delta =3.0$ penetrates as deeply as one with ${E}_{{\rm{c}}}=70$ keV and $\delta =5.0$. Interestingly, for distributions that have both very low values of E_c ($\lt 8$ keV) and high δ ($\gt 6$), much of the energy is deposited directly in the corona.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Many previous studies (e.g., Abbett & Hawley 1999; Allred et al. 2005, 2006; Cheng et al. 2010) implemented the analytic expression for beam penetration depth derived by E78. That expression includes non-uniform ionization (Hawley & Fisher 1994) but does not include relativistic effects. Since the E78 expression has been so widely used, it is informative to compare it with the more complete treatment that we described in Section 3. In Figure 7, we plot the ratio of beam penetration depths obtained from E78 to that derived by solving the Fokker–Planck equation. We find E78 works quite well in the the range ${E}_{{\rm{c}}}\lt 40$ keV and $\delta \gt 4.5$ but becomes increasingly inaccurate for higher E_c and lower δ. In this regime, E78 predicts a penetration depth as much as seven times greater. As the cutoff energy increases, relativistic effects become significant. Due to the Lorentz length contraction, relativistic electrons experience a higher ambient density than would be expected from classical theory resulting in less penetration. Relativistic effects become less pronounced with increasing δ since there are fewer high energy electrons in the distribution.

Zoom In Zoom Out Reset image size

Figure 8 shows the ratio of the penetration depth for highly beamed non-thermal electrons to an isotropic pitch angle distribution. As expected, the beamed electrons penetrate more deeply. As δ increases this effect becomes more pronounced. This is because higher energy particles undergo more collisions before stopping, so their initial pitch angle is less important in determining their final stopping depth. A larger δ means fewer high energy particles are in the distribution, so the initial pitch angle has a more pronounced effect.

Zoom In Zoom Out Reset image size

Figure 9 compares the beam penetration depths as a function of cutoff energy and spectral index for several loops. The left panel shows the ratio of the penetration depths for beams propagating in QS.SL.HT to QS.SL.LT. For low values of the cutoff energy, the beam penetrates to a column mass approximately twice as great in the QS.SL.LT loop. In this loop, the coronal temperature and density are lower than in QS.SL.HT. It is the low energy electrons that are most affected by this difference. The middle panel compares beam penetration depths in QS.LL.LT with QS.SL.LT. In this case, beams with low-cutoff energy penetrate less deeply in QS.SL.LT. Since QS.SL.LT is a shorter loop than QS.LL.LT, the coronal density is greater resulting in less penetration for the lowest energy electrons. The right panel compares QS.SL.LT with SS.SL.LT. These loops are nearly identical in their coronal and upper chromospheric regions. They differ significantly in the photosphere, but very few electrons are able to penetrate to that depth so beam penetration depths are similar in these loops in the range of cutoff energies and spectral indices studied here. However, since they have very different photospheric temperatures, these loops predict significantly different levels of white light emission. This difference will be important in studying the white light emission predicted from flaring loops.

Zoom In Zoom Out Reset image size

6.2. Collision Rates

During flares, direct excitations and ionizations by beam particles can be much larger than thermal rates in the region of beam impact.⁵ To understand the relative importance of collisional ionizations by flare-accelerated particles, it is informative to compare them with the thermal rates. In Figure 10, we plot the ionization rates, ${n}_{i}{C}_{\mathrm{beam}}^{\mathrm{nt}}$, in several loops. Here n_i is the ground state number density. We find the ionization rates increase dramatically around beam impact. For example, in the solar loops the H i 1-c rates peak at $2\times {10}^{13}$ s⁻¹. To put this into perspective, the thermal rates in this region are $\sim {10}^{3}$ s⁻¹ so some 10 orders of magnitude smaller. These results are relatively independent of the initial loop conditions, with the exception of the M dwarf loop. ${n}_{i}{C}_{\mathrm{beam}}^{\mathrm{nt}}$ increases deeper in the M dwarf loop compared to the solar loops because in the former the chromosphere forms deeper.

Zoom In Zoom Out Reset image size

It is useful to model how these collision rates vary with injected electron beam spectra. Figure 11 shows the peak non-thermal collisional ionziation rates for electron beams injected into the QS.SL.HT loop as a function of E_c and δ. H i 1-c and He i 1-c rates have similar dependence on these parameters, peaking at lower E_c and higher δ. Q_beam is narrower and higher peaked in this regime (see the middle and bottom panels of Figure 5) resulting in a higher peak collision rate. However, the collision rates drop quickly for very low E_c. There much of the beam is stopped in the transition region and corona where n_i is low. The He ii 1-c rates rapidly decrease for increasing E_c. A beam with larger E_c penetrates to a deeper, cooler, and hence lower He ii density region of the atmosphere.

Zoom In Zoom Out Reset image size

6.3. Return Current Heating in the Corona

To understand how return currents are likely to affect flare dynamics, it is useful to compare their heating rates (Q_rc) to the heating rates produced directly by collisions (Q_beam). In the following discussion, it should be kept in mind that we are comparing only coronal heating rates—not those in the chromosphere where Q_beam dominates. A comparison of these rates in several loops is shown in Figure 12. In the QS.SL.HT loop, Q_rc and Q_beam are similar in size. However, in all other cases, Q_rc dominates by at least an order of magnitude, being stronger in the cooler loops. This is because the resistivity is strongly and inversely dependent on the coronal temperature. Hot loops have much less resistance so are heated by Q_rc less. Therefore, we speculate that return currents are likely to be most important in the early phases of flares before the corona has been heated to very high temperatures. This figure also compares Q_rc for the case of highly beamed electrons (solid lines) and isotropically injected electrons (dotted lines). These rates are similar in all cases indicating that the pitch-angle distribution of the injected electrons makes little difference in Q_rc.

Zoom In Zoom Out Reset image size

From the results of this parameter study, we can also compare how Q_rc varies with E_c and δ. Figure 13 shows the volumeteric heating rate due to return current in the corona of the QS.SL.HT loop as a function of these parameters. The same energy flux (10¹⁰ erg cm⁻² s⁻¹) was injected for every simulation in this study, so simulations with higher E_c have fewer total injected electrons. Return current heating scales with the particle flux (F_e) so fewer particles results in less heating. The energy per particle factors into F_e only through its velocity, so the return current heating rate is only weakly dependent on δ.

Zoom In Zoom Out Reset image size

We have developed a computational framework capable of modeling the evolution of magnetic flux loops in response to flare heating. This model can be applied to loops on the Sun and on dMe flare stars, and it includes many of the physical processes most important to flare dynamics. Our code solves the optically thick, non-LTE radiative transfer equation for many transitions important in the chromosphere where much of the flare energy is deposited allowing us to model flare emission which can be compared directly to observations. We have implemented a method for solving the Fokker–Planck equation describing how flare-accelerated particles interact with the ambient plasma and have used that to model the heating, momentum deposition, collisional excitations and ionizations, and return currents produced by these particles.

With our model, we have performed a parameter study to understand how these phenomena depend on the energy spectrum of a beam of flare-accelerated particles. We have assumed a power-law energy distribution for these particles and found that their penetration depth is strongly dependent on the beam cutoff energy and power-law index and more weakly dependent on their initial pitch-angle distribution. For distributions with low power-law index or high cutoff energy, relativistic effects become important causing a penetration to depths as much as seven times less than predicted by classical theory. Varying the conditions of the loops (i.e., the coronal temperature and density and photospheric temperature) has a small effect on the penetration depth for distributions with low cutoff energy. The collisional ionizations produced by impacts from these flare-accelerated particles also are strongly dependent on cutoff energy and power-law index. They can be 10 orders of magnitude greater than thermal collision rates. The beam particles also induce a return current which deposits heat in the corona. We have found that, depending on coronal conditions, return current heating in the corona can be more than an order of magnitude greater than coronal beam heating.

In this parameter study, we have focused on understanding the effects of flare-accelerated electrons. In future work, we will extend this to flare-accelerated ions which are also known to be important to flare energetics. This work has focused on the initial response of flaring loops. However, a key feature of our framework is its ability to model the dynamical evolution of these loops. In future work, we will present a parameter study exploring their evolution in response to flare-accelerated particles.

This work has been supported by grants through the NASA Heliophysics Supporting Research and Technology and the NASA Living With a Star programs. This research was aided by many useful discussions with the participants of the Chromospheric Flares meeting held at the International Space Science Institute (ISSI) in Bern, Switzerland. A.F.K. acknowledges the support of the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by Oak Ridge Associated Universities through a contract with NASA, and from support through UMCP GPHI Task 132. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement nr 291058 and by the Research Council of Norway through the grant "Solar Atmospheric Modelling" and through grants of computing time from the Programme for Supercomputing. We thank G. Holman for helpful discussions about return current heating.