Hydrogen Balmer Line Broadening in Solar and Stellar Flares

In Kowalski et al. (2015b) and Kowalski et al. (2016), we presented model flare atmospheres heated by a high flux of nonthermal electrons (10¹³ erg cm⁻² s⁻¹; F13) which produce a T ∼ 10⁴ K blackbody-like continuum flux distribution as observed in many flares from active M dwarf stars. The hot blackbody-like continuum distribution results from large continuum optical depths from a dense, heated chromospheric condensation. The charge density in the continuum-emitting layers achieves a maximum value of $\sim 5\times {10}^{15}$ cm⁻³. The opacity effects from level dissolution produces Balmer continuum flux¹⁵ at wavelengths longer than the Balmer limit (3646 Å) and higher order Balmer emission lines that fade into this continuum flux. These properties have been observed in high-resolution spectra of a large dMe flare (Fuhrmeister et al. 2008). With the S78 approximations as Voigt profiles in the RH code, the highest order Balmer lines that are predicted at this charge density are H10/H11, which is reasonable compared to the observations of some dMe flares (Kowalski et al. 2013, 2016).

Because of the computation time required for the convolution (Equation (1)), it is unfeasible to incorporate the TB09+HM88 method directly into RADYN. Instead, we input snapshots of the atmospheric state from RADYN into RH. It is well-known that the flaring atmosphere is very dynamic. The dynamics and the proper time-dependent ionization (electron density) at each snapshot are included, but statistical equilibrium for the excitation is assumed in solving the radiative-transfer equation in the RH calculation. To understand the extent of this limitation, we plot Balmer lines modeled using RADYN and RH in Figure 5. For this comparison, RH has been configured to most closely match the radiative-transfer method in RADYN, (i.e., using a 6-level hydrogen atom, complete redistribution, and S78 broadening). In this configuration, the only difference between the radiative-transfer methods employed by these codes is the assumption of statistical equilibrium in the excited levels employed in RH. Figure 5 indicates that, for the $t=2.2\,{\rm{s}}$ snapshot in the F13 model from Kowalski et al. (2015b), the statistical equilibrium limitation is a small effect, justifying our use of RH. The same is true for $t=2.2\,{\rm{s}}$ of the F12 model presented in Kowalski et al. (2015b).

Zoom In Zoom Out Reset image size

In RH, the occupational probability formalism is not included in the rate equations for statistical equilibrium as described in detail in Hubeny et al. (1994). In our implementation, we include the occupational probability prescription only in the non-LTE bound–bound and bound–free opacity and emissivity (Section 2, see also Kowalski et al. 2015b). Including occupational probabilities in the rate equations for statistical equilibrium in RH would likely affect the electron density and the populations of the upper levels of hydrogen. At this expense, we use the non-equilibrium electron density from a snapshot from RADYN. The agreement with the observed spectrum of Vega (Section 3) justifies that the method currently employed in RH is sufficiently accurate for atmospheres that are near LTE conditions in their continuum-emitting layers, as is the case for the F13 flare models at $t=2.2\,{\rm{s}}$ (see the discussion in Kowalski et al. 2015b).

We repeat the RH calculation with the TB09+HM88 profiles convolved with the Voigt function for the $\delta =3$, F13 model atmosphere at t = 2.2 s from Kowalski et al. (2016) and for the double-power-law ($\delta =3$ at $E\lt 105\,\mathrm{keV}$, $\delta =4$ at $E\gt 105\,\mathrm{keV}$) F13 model at $t=2.2\,{\rm{s}}$ (hereafter, F13 "dpl") from Kowalski et al. (2015b). The new prediction in the Balmer jump wavelength region is shown in Figure 6 for the F13 $\delta =3$ at 2.2 s, which exhibits broader Balmer lines, while H9 is the highest order Balmer line that is not completely transformed into Balmer continuum flux longward of the Balmer edge. With the S78 method, H11 is the bluest detectable (undissolved) Balmer line in emission. We use the TB09+HM88 profiles, which have been calculated with ${\beta }_{\mathrm{crit}}$ but find that these increase the peak of the H8 Balmer line and decrease the trough between H8 and H9 by ∼2% compared to a calculation without ${\beta }_{\mathrm{crit}}$ (VCS, no psuedo-continuum). Using the TB09+HM88 profiles, the H10 and H11 lines are broadened much more than with the S78 prescription (Figure 1); the occupational probability for the upper levels of these transitions are small, which causes the flux at these wavelengths to be dominated by dissolved level continuum flux in the emergent spectrum. Higher order Balmer lines in the impulsive phase that are as faint (relative to the nearby continuum) as the new F13 prediction are very rare (e.g., García-Alvarez et al. 2002).

Zoom In Zoom Out Reset image size

We convolve the F13 model spectra in Figure 6 with a representative spectral resolution ($R\sim 450$) from the flare spectral atlas of Kowalski et al. (2013) and calculate the width at 10% of maximum line flux (0.1 width) for the Hγ line to be 78 Å for the TB09+HM88 prescription and 32 Å for the S78 prescription. For the F13 dpl at $t=2.2\,{\rm{s}}$, the 0.1 width of Hγ is larger, ∼100 Å with the TB09+HM88 broadening. Typical line widths in dMe flares are 15–20 Å for large flares at moderate spectral resolution (Hawley & Pettersen 1991); the maximum observed 0.1 widths for the Hγ line at low spectral resolution are $30\mbox{--}40\,\mathring{\rm{A}}$ (Figure 4.13 of Kowalski 2012) at times of brightest line emission, but the line broadening has been observed to be larger (45–50 Å) in the mid-rise phase of large flares (Figure 4.14 of Kowalski 2012). The contribution function-weighted charge density over which the Hγ line forms is $\sim 3.5-5\times {10}^{15}$ cm⁻³ for the F13 $\delta =3$ (and $4\mbox{--}6\times {10}^{15}$ cm⁻³ for the F13 dpl) at $t=2.2\,{\rm{s}}$. The F13 model with TB09+HM88 profiles demonstrates that the F13 instantaneous spectrum at 2.2 s, and therefore this range of charge density, is not consistent with even the largest observed values of the broadening in flares. In Section 4.2, we show that adding flux spectra with lower broadening (at earlier and later times in the heating simulation) to emulate simultaneously heating and cooling loops produces a broadening that is closer to the observations. Because of the larger broadening with the TB09+HM88 profiles, the emergent line intensity in the F13 model originates from a larger physical depth range in the atmosphere, and the wavelength-integrated emergent line flux of the Balmer lines is therefore also several times (∼2.5×) larger than in the case of the S78 profiles.

In Figure 7, we plot the widths of the Hγ line for the evolution of the F13 dpl simulation compared to the width of the line from the F12 dpl simulation from Kowalski et al. (2015b) in order to bracket the range of electron densities relevant for flares (${10}^{13}-5\times {10}^{15}$ cm⁻³). For the F13 dpl model evolution in Figure 7, the slope of the relation is approximately 2/3, which is the expected scaling of the wavelength shift from electric pressure broadening (see, e.g., Johns-Krull et al. 1997):

$$\begin{eqnarray}&&{\rm{\Delta }}\lambda \propto {n}_{e}^{2/3}{\left\{\displaystyle \frac{{n}^{2}}{{n}^{2}-4}\right\}}^{2}[n(n-1)+2].\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

It is interesting to understand the difference between density estimates using TB09+HM88 and S78 for the $t=2.2\,{\rm{s}}$ of the F13 $\delta =3$ simulations in Figure 6. Using the S78 theory gives a 0.1 Hγ width of 16 Å, which is a factor of nearly five lower than with the TB09+HM88 profiles. Using the relation from the TB09+HM88 broadening models in Figure 7, a broadening of 16 Å gives an electron density that is a factor of 5–10 lower than the range in the RHD (RADYN) calculation. The results of Johns-Krull et al. (1997) demonstrated a similar discrepancy between the S78 treatment and the modified impact theory of Kepple & Griem (1968) and Bengtson et al. (1970).

According to Equation (2), the widths of higher order lines should increase significantly (e.g., Hα compared to Hγ), in the regime where thermal broadening is comparatively small. However, the interline optical depth variation results in much less extreme width differences (Svestka 1962, 1963). In our models (and generally in the observations too), the broadening is not strongly dependent on the Balmer transition because the much larger optical depth for Hα near line center leads to a much smaller (∼100×) physical depth range over which the line is produced compared to Hγ, thus giving a larger value of the 0.1 width for Hα than indicated by Equation (2). We also note that directly comparing to the width of the S(α) profiles results in largely erroneous values of the electron density; due to the optical depth over the line, the values of S(α) in the far wings at S($\alpha )\sim {10}^{-4}$ are the important determinant in the 0.1 width of the line when the optical depth near line center and the emissivity in the line wings are both large, as in the Balmer line profiles of Vega.

In Figure 7, the width of Hγ is shown from a new electron beam heating simulation with a flux of 5 × 10¹² erg cm⁻² s⁻¹ (5F12) that is calculated with the RADYN flare code, using a high, low-energy cutoff of ${E}_{c}=150\,\mathrm{keV}$ so that the heating occurs deep in the atmosphere (see Appendix A). The emergent continuum spectrum has a small Balmer jump ratio and the low-order Balmer lines are broad and in emission with a prominent central reversal. Like for Vega the Hγ line is formed over a large range of electron densities, given by the error bar range in Figure 7, which demonstrates that a line width value does not always correspond to a unique electron density value; this is also true for the F13 model atmospheres but for a narrower range of electron density. Interestingly, the broadening (75 Å) in the Hγ emission line in the F13 model is comparable to the broadening (nearly 70 Å) in the Hγ absorption line in Vega, while there is an order of magnitude difference in the maximum charge density over the formation heights of the line profile. The Hγ line in the 5F12 E_c = 150 keV model forms over an electron density that is comparable to and larger than the electron density that produces the line in Vega, but the broadening is much larger in the Vega spectrum.

The model predictions of the highest order Balmer line in absorption (at least H16) in the Vega spectrum (Figure 3 top) differ significantly from the highest order Balmer line that is in emission (H9) in the F13 spectrum. Therefore, the amount by which the higher order Balmer lines fade into the dissolved level continuum flux at wavelengths longer than the Balmer edge can provide an additional constraint to break such a degeneracy and ambiguity in the charge density-line width relation in flare atmospheres. In the 5F12 model, the highest order Balmer line that is not completely dissolved is ≈H13, and the far wing of Hγ (${\lambda }_{\mathrm{rest}}+20$ Å) forms over an electron density that is three times as great as that in Vega (at the same wavelength) but a factor of three lower than the F13 model (at the same wavelength). This ordering by electron density is consistent with the highest order Balmer line in each model spectrum. Kowalski et al. (2015b) discussed a dissolved line decrement¹⁶ using H11/Hγ that can be used to quantify the degree by which the higher order lines fade into continuum flux; we intend to pursue this diagnostic in a future work with NLTE flare models and the TB09+HM88 broadening.

F13 beam flux model atmospheres produce dense, downward-directed, heated compressions ("chromospheric condensations") and have been found to adequately explain the hot (10⁴) blackbody-like, optical and NUV continuum flux distribution of flares on active M dwarf stars (Kowalski et al. 2015b, 2016). However, the new broadening predictions show that these white-light emitting compressions result in a charge density (∼5 × 10¹⁵ cm⁻³) that over-broadens the Balmer lines as in Figure 6 compared to observations in the literature. If such dense compressions are formed in dMe flares, then they would have to account for a large fraction of the continuum flux and a small fraction of the line flux, such that the spatially integrated flare spectrum exhibits narrower lines than the TB09+HM88, F13 $t=2.2\,{\rm{s}}$ prediction.

Kowalski et al. (2015b) used a time-average of the F13 model atmosphere over its evolution to more realistically simulate the prediction over an exposure time. Averaging over the F13 simulation is equivalent to an instantaneous measurement of a spatially unresolved observation of many F13 flare loops sequentially being initiated at regular, short time intervals corresponding to different locations on the star. This prescription for representing the flare flux from many flare "threads," or kernels, is similar to the "multithread" modeling that has been successful in reproducing spatially unresolved solar flare light curves in GOES soft X-rays Warren (2006). A method for multithread modeling with RADYN heating models lasting $\sim 15\mbox{--}20\,{\rm{s}}$ was further developed in Rubio da Costa et al. (2016) and has also been employed in Reep et al. (2016). In multithread modeling of solar flares, the flux of each thread changes over the evolution of the flare (Warren 2006) such that impulsive phase threads have a higher flux than gradual phase threads. For simplicity, we use the F13 model to represent all heated threads (kernels).

We simulate a multithread prediction for the line broadening from the F13 simulations in Kowalski et al. (2015b) and Kowalski et al. (2016). We use the snapshots from Figure 7 at t = 0.2, 1.2, 2.2, and 4.0 s (the electron beam heating is turned off at 2.3 s and allowed to relax until t = 5 s) to sample the range of atmospheric conditions and brightness values achieved in the models. The assumption of statistical equilibrium (without the occupational probability formalism of Hubeny et al. 1994) was justified in Sections 3 and 4.1 but may not be appropriate within the earliest tenths of seconds after the beam heating starts. The Balmer lines are broadest at t = 2.2 s when the chromospheric condensation has achieved the maximum electron density at $T\sim {\rm{12,000}}\mbox{--}{\rm{13,000}}\,{\rm{K}}$, thus averaging over earlier and later times decreases the line broadening in the multithread models. For the F13 dpl multithread (0–5 s average) model, the 0.1 width of Hγ is ∼50 Å, a factor of two narrower than at $t=2.2\,{\rm{s}};$ for the F13 $\delta =3$ multithread (0–5 s average) model, the 0.1 width is 40 Å (compared to 75 Å at 2.2 s; Figure 6). The 0.1 width of the F13 $\delta =3$ is near the upper range of observed values in the early impulsive phase but are still quite high. Whereas the S78 predictions from the F13 model at 2.2 s for the continuum flux ratios¹⁷ (C3615/C4170, C4170/C6010), the line-to-continuum ratio values of Hγ/C4170 (=10–20), and broadening are generally consistent with peak phase spectra of dMe flares (Kowalski et al. 2015b, 2016), the time-averaged values with TB09+HM88 broadening are more consistent with the early impulsive rise phase of some dMe flares in the literature, as concluded in Kowalski et al. (2015b) for the time-average continuum properties of the F13 dpl model. For example, the time-resolved impulsive phase data of the large IF3 flare event from Kowalski et al. (2013) exhibit values of (C3615/C4170, C4170/C6010, Hγ/C4170, 0.1 width Hγ) = (2.2, 1.8, 40, 48 Å) in the early-to-mid rise phase (S#27, cf Figure 30 of Kowalski et al. (2013)) that are very similar to these quantities (2.1, 1.8, 45, 45 Å, respectively) from the F13 $\delta =3$ multithread model (Table 1).

Table 1.  Balmer Decrements and Hγ Broadening with TB09+HM88

Model	Hα	Hβ	Hγ	Hδ	0.1 width Hγ (Å)	Hγ/C4170	C3615/C4170	C4170/C6010
pre-flarea	5.9	2.7	1.0	0.58	⋯	⋯	⋯	⋯
F13 2.2s dpl	0.3	0.7	1.0	0.8	100 (100)	45	1.8	2.1
F13 dpl multithread (0–5 s ave)	0.4	0.8	1.0	0.9	52 (60)	74	2.6	1.6
F13 2.2s $\delta =3$	0.3	0.8	1.0	0.9	75 (78)	31	1.7	2.3
F13 $\delta =3$ multithread (0–5 s ave)	0.4	0.9	1.0	1.0	40 (45)	45	2.1	1.8

Note.

^aThe pre-flare atmosphere is the most up-to-date version of the pre-flare state in our modeling, which is described in Appendix A. The continuum flux ratios C3615/C4170 and C4170/C6010, and the line-to-continuum flux ratio (Hγ/C4170) are also given. The values in parentheses for the Hγ widths are calculated with an instrumental convolution of R = 450.

Download table as:  ASCIITypeset image

Balmer decrements are defined as the ratio of the excess (i.e., background subtracted) flux in Balmer lines to that of the Hγ line and are also used to constrain the charge density variations in flares (Hawley & Pettersen 1991; Jevremovic et al. 1998; García-Alvarez et al. 2002). The observed decrements in dMe flares are typically Hβ/H $\gamma \sim 1.0\mbox{--}1.2$ and for Hδ/H $\gamma \sim 0.75\mbox{--}1.05$ (Table 4.20 of Kowalski (2012), Table 1 of Allred et al. 2006). Decrements are similarly used in solar flare studies as well, but the lack of broad wavelength coverage spectra make the measurements rare in the modern era (Johns-Krull et al. 1997; Kowalski et al. 2015a). Since Balmer decrements have been widely used, it is informative to study how they may vary in response to electric pressure broadening with the new TB09+HM88 profiles.

In Table 1, we show the decrements for the two F13 models at $t=2.2\,{\rm{s}}$ and the 0–5 s average (multithread) model (Section 4.2) calculated with the TB09+HM88 broadening. Although the line flux increases by ∼2.5× with TB09+HM88 compared to the line flux calculated with the S78 broadening, the line flux ratios do not change appreciably between the two methods. The Hβ/Hγ and Hα/Hγ line flux ratios are much smaller in the models than typical values from observations, but the Hδ/Hγ ratios are in general agreement. Such small values of Hα and Hβ to Hγ indicate a "reverse decrement" and are only observed to be as low as ∼0.8 (for Hα; Figure 4.22 of Kowalski 2012). The Hα/Hγ line flux ratios are very low in the F13 models because the optical depth necessary to produce the hot blackbody-like continuum shape results in a very large Hα optical depth and radiation thus escapes from a much smaller physical depth range of the atmosphere across this line compared to Hγ. We have also computed the line flux ratio at t = 0 s in our M dwarf model (row 1), which is in good agreement with observations (Bochanski et al. 2007, see also Table 2.8 of Kowalski et al. 2013) but with a significantly larger line flux ratio for Hβ and slightly larger ratio for Hα than in the observations. The discrepancy could result from the difference between spatially resolved and averaged observations, and 3D effects (Walkowicz 2008; Uitenbroek & Criscuoli 2011; Leenaarts et al. 2012; Wedemeyer & Ludwig 2015). As in the flare simulations, the Hδ/Hγ line flux ratio is consistent with the observations.

It is interesting to compare the values in Table 1 to the method of Drake & Ulrich (1980), which presented the Balmer line decrements for a large parameter space of electron density, temperature, and optical depth for a homogeneous slab of hydrogen. In addition, approximations were made using the escape probability from the line wing and the optical depth variation from line to line. In their Appendix B, they present the decrements for ${n}_{e}={10}^{15}$ cm⁻³, which is the highest density that they considered: the Hα/Hγ, Hβ/Hγ, and Hδ/Hγ decrements are 1.8, 1.2, and 0.85, respectively. The values for their limiting optically thick case, on the other hand, show reverse decrements for Hα and Hβ but are not as extreme as for the F13 model at $t=2.2\,{\rm{s}}$, likely due to the factor of five larger electron density in the F13 model. The results from Drake & Ulrich (1980) are qualitatively similar to ours for flare atmospheres that are so dense that they are approximately homogeneous in electron density (e.g., at $t=2.2\,{\rm{s}}$ of the F13 dpl model, the contribution function-weighted electron density varies from only $(4\mbox{--}6)\times {10}^{15}$ cm⁻³ over the Hγ line), but would not be applicable to the flare atmospheres that exhibit a large range of electron density over the line formation, such as the high, low-energy cutoff (5F12, E_c = 150 keV) flare model in Figure 7 with a range of nearly a hundred in electron density.

Following Kowalski et al. (2012), we interpret the surface flux spectra components (Figure 9) in the S#24 and S#113 spectral observations of the YZ CMi Megaflare using analogous phenomenology from high spatial resolution data of two-ribbon solar flares.

• 1.  
  F_kernels: The F13 multithread model spectrum (0–5 s average of the F13 dpl model) represents the spreading ribbons of the Megaflare. As the ribbons spread, new flare loops form, and their initiation is staggered in time so that the spatially integrated flux from this area is equivalent to the 0–5 s time-average of the F13. The area decreases over time (X_kernels decrease from 0.0009 at the time of S#24 to 0.0005 at the time of S#113); the beam flux probably also decreases over the gradual phase (Warren 2006) as lower field strengths reconnect higher in the corona. We note that newly heated flare loops are necessary to explain the temporally extended gradual phase of the soft X-ray flux in spatially unresolved GOES light curves of solar flares (Warren 2006). The spreading ribbons may correspond to the regions that produce the end of the large secondary flare (e.g., at $t=0.6\,\mathrm{hr}$ or $t=1.1\,\mathrm{hr}$ in Figure 8, before S#24) that peaked prior to the spectral observations, and/or they may be the continuation of the spreading ribbons from the BFF event with ${\rm{\Delta }}U=-5.8$ mag at $t=0.45\,\mathrm{hr}$ in Figure 8. Recent high spatial resolution images in the Hα red wing ($+0.8$ Å) with the New Solar Telescope have shown that spreading Hα ribbons are composed of many fine-scale kernels with a size of ∼100 km (Sharykin & Kosovichev 2014; Jing et al. 2016); these kernels are located at the leading edge (Isobe et al. 2007) of the ribbons and are often associated with red-wing asymmetries in Hα. Other chromospheric flare lines, such as Mg ii and Fe ii, also exhibit red-wing asymmetries. Kowalski et al. (2017) reproduces the general properties of the red-wing asymmetries in Fe ii lines in the brightest flare footpoints constituting the two flare ribbons that spread apart in the impulsive phase of the 2014 March 29 X1 solar flare. The RHD model that produces dense, heated chromospheric compressions that are consistent with the spectral observations employed a relatively high electron beam flux (5F11). The F13 model evolution consists of a dense, heated chromospheric condensation but with higher temperature, density, and continuum optical depth compared to the chromospheric condensation properties in the 5F11 model for solar flares (see Appendix C of Kowalski et al. 2017). Thus, the kernels of the spreading ribbons may consist of many ≈F13 beams in dMe flares and lower fluxes, ≈5F11 beams, in solar flares.
• 2.  
  F_decay: We use the F13 t = 4 s flux to represent the wake of the spreading flare ribbons F_kernels, resulting in long-lasting (temporally extended) decay phase flux in the Megaflare that persists for hours and is present through the spectral observations S#24 and S#113. Note that the F13 multithread model (0–5 s average) also consists of this decay phase snapshot (t = 4 s), but its contribution in the 0–5 s average represents the rapid decay phase of a kernel. The gradually decaying flux may be analogous to temporally extended decay flux in the "wake" of separating ribbons in solar flares, such as that observed after the peak of the NUV ($\lambda 2826$) continuum-emitting kernels in the 2014 March 29 X1 solar flare (Heinzel & Kleint 2014; Kowalski et al. 2017) and the temporally extended gradual phase component of optical flare kernels on the Sun (Kawate et al. 2016). We assume that the first large event in the Megaflare (or any of the secondary large events peaking before S#24) produces a wake of persistent, decaying flux behind the spreading ribbons. In our modeling, this flux component is not self-consistent in the timing relative to one of these peaks in the Megaflare; the NUV continuum exponential decay constant after $t=2.3\,{\rm{s}}$ in the F13 model is $\tau \sim 0.6\,{\rm{s}}$, meaning that the pre-flare continuum level would be reached several seconds after the beam heating ends (the F13 models are very computationally demanding, preventing us from currently following the evolution for hours after the beam heating). In Section 5, we found that this flux component is necessary in order to lower the Balmer line decrements and broadening in the model to the observed range of values in S#24. In Section 7, we speculate on directions of future work to better model and constrain the origin of the temporally extended decaying flux component. We use the high spatial resolution data in the Hα red wing from the DST/IBIS presented in Kowalski et al. (2015a) to estimate the ratio of newly formed kernel area in the spreading flare ribbons to the the area of the intensity in the wake of the ribbons. The IBIS data have a cadence of 20 s and covered several peaks in the X-ray impulsive phase of a two-ribbon, long-duration C1 solar flare, SOL2011-08-18T15:15. For each frame, we measure the area of bright flare intensity that is not bright in the previous frame. We only consider the umbral ribbon, which spreads very slowly across the umbra; most of the apparent kernel motion is parallel along the ribbon. The area of the newly formed Hα kernels is only a small fraction (0.1–0.2) of the total Hα kernel area emitting above the threshold count rate. Thus, it is reasonable to assume that the decay flux ("ribbon wake"; F13 t = 4 s) and the newly heated kernels (F13 t = 0–5 s average) could reasonably exist in the areal ratio that we employ (1/25) for the DG CVn Superflare and the YZ CMi Megaflare events.
• 3.  
  ${F}_{{hotspots}}={F}_{{new}{ribbons}}$: During the rise and peak of the secondary flare MDSF2, new flare ribbons consisting of new kernels develop. From our RHD modeling of the new flare flux in the secondary flare with the 2 × 10¹² erg cm⁻² s⁻¹ ${E}_{c}=500\,\mathrm{keV}$ beam (Appendix A), we infer that the heating scenario in the secondary flare is strikingly different than the heating in the F13 models with ${E}_{c}=37\,\mathrm{keV}$, where the electron beams lose their energy in the upper chromosphere. In our terminology, hotspots form in the lower chromosphere, which is necessary to produce the Vega-like flare spectrum, and kernels produce chromospheric condensations in the upper chromosphere. ${F}_{\mathrm{newribbons}}$ may produce F13-type heating in addition to hotspot heating. In MDSF2, we infer that the dominant heating goes into forming the hotspot. The filling factor of the hotspot, X_hotspot, is 0.001, which is remarkably similar to the filling factor inferred from a blackbody fit to the blue continuum with T = 10,000 K (Kowalski et al. 2010). The filling factor of F_decay is 12.5× greater than this, which was also inferred in Kowalski et al. (2010) by fitting an optically thin Balmer continuum (Allred et al. 2006) to the flux at $\lambda \lt 3646$ Å to represent the decaying Hα ribbons from the main peak. In the three-component phenomenological model of Kowalski et al. (2012), the F_kernels component was modeled as a second (cooler) hotspot in the photosphere.

In the superflare from DG CVn (Osten et al. 2016), F_kernels was used to represent the newly heated loops after the peak of the high-energy secondary flare event (referred to as F2). The energetic secondary flare F2 had an impulsive phase timescale that was much longer than the "big first flare" event of comparable energy an hour earlier. Thus, different heating scenarios may dominate through the BFF and F2 events in the DG CVn Superflare. A significant detection of continuum flux with a cool color temperature was obtained from the UVW2/V-band and V-band/R-band flux ratios, but the only simultaneous observations to constrain this component occurred after the peak of the F2 event. The continuum flux ratios were adequately explained by the F13 multithread model with decay phase flux without a hotspot contribution, which is the same scenario for the model for the decay phase spectrum S#24 from the Megaflare. When significant time has passed after the peaks of energetic secondary flares, ${X}_{\mathrm{hotspot}}\sim 0$, such as for the YZ CMi Megaflare at the time of S#24 and for the DG CVn Superflare F2 event at $\approx T0+11750\,{\rm{s}}$ and $\approx {\rm{T0+17,000}}\,{\rm{s}}$ (Osten et al. 2016). In future work, we will constrain X_hotspot(t) using a superposition of sequentially heated hotspot spectra, as done to produce the F13 multithread model.

In solar flares, there is also evidence that the characteristics of the beam heating can significantly change between the main flare peak and secondary flares. Warmuth et al. (2009) describe a secondary flare that produced a hard X-ray spectrum that could be explained by an electron beam with a higher, low-energy cutoff (${E}_{c}\sim 110$ keV) than in the initial impulsive phase. The parameters of the electron beam for this flare are less extreme but qualitatively similar to the parameters for the ${E}_{c}=500\,\mathrm{keV}$ electron beam that we use to produce a hotspot spectrum. Our modeling of the hotspot implies that the heating may result from ultrarelativistic electrons. Synchrotron emission from ultrarelativistic ($E\gt 500\,\mathrm{keV}$) electrons are a possible source for the sub-mm/THz radiation in solar flares (see the review in Krucker et al. 2013). Sub-THz radiation attributed to synchrotron emission has also been detected from stellar flares in active binary systems (Massi et al. 2006).

In Figures 11(b)–(e), we show representative images from the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on the Solar Dynamics Observatory (SDO) in the 1700 filter during the GOES X5.4 and X1.3 solar flares of 2012 March 7. The AIA 1700 filter shows the footpoint features (ribbons and kernels) in solar flares at high spatial resolution. There are no simultaneous spectral observations over the wavelength range of this filter, but studies of other solar flares indicate a combination of continuum, line intensity and pseudo-continuum intensity from blended lines (Cook & Brueckner 1979; Doyle & Cook 1992; Brekke et al. 1996; Qiu et al. 2013). The FUV continuum becomes bright in stellar flares (Hawley & Pettersen 1991), and it may be a good tracer of the formation of hot spots in secondary flares (Ayres 2015). The solar flare event of 2012 March 7 provides a unique comparison to the Megaflare because it produced two significant flares in the 1700 light curve (Figure 11(a)), analogous to a "big first flare" and a high-energy secondary flare, within the same active region. Notably, the thermal response detected by GOES shows a significant difference compared to the secondary flare in Warmuth et al. (2009), which does not respond in GOES, suggesting that the heating mechanism is not uniform among secondary flare events.

Zoom In Zoom Out Reset image size

In Figures 11(b)–(e), we use the spatial development in the 2012 March 7th solar flare to illustrate possible analogous development of the three continuum-emitting components (${F}_{\mathrm{kernels}},{F}_{\mathrm{decay}},{F}_{\mathrm{hotspot}}$) in the YZ CMi Megaflare. Combining the spatial development in Figure 11 with the RHD modeling of the Megaflare decay phase spectra, a possible scenario for the spatial development of the Megaflare decay phase is as follows: (1) at the start of the spectral observations, bright kernels with a filling factor of 0.1% are heated by high-flux electron beams at the leading edges of the spreading ribbons (panel d); (2) each kernel rapidly decreases in brightness as beam heating ends. The total flare area that has been swept out by these kernels and other ribbons earlier in the flare has a filling factor of $\sim 1 \% \mbox{--}2.5$% and decays in brightness on a longer timescale than the rapid brightness decay of each kernel (panel d); (3) later, the total area of the kernels at the leading edges of the separating flare ribbons has decreased in size to a filling factor of 0.05%. These kernels propagate into the locations of the active region where new ribbons with hotspots heated by ultrarelativistic electron beams are triggered and develop with a filling factor of 0.1% (panel e). Also at this time, the locations of the kernels formed in panel (d) are still decaying on a long timescale.

The Hα/Hγ and Hβ/Hγ decrements are not well-reproduced by the 0–5 s average (multithread model) using the F13 because a high optical depth over the continuum-emitting regions leads to an optical depth in these lower order lines that is very large. A larger filling factor of the F13 flux at t = 4 s (after the beam heating ends) compared to the impulsively heated kernels is required to produce a better match to the line flux ratios.

The large filling factor (1%–2% of the visible stellar hemisphere in the Megaflare) suggests that this missing component is due to the decaying flux from previous, larger events (e.g., the main event, or the secondary flare just prior to the spectral observations). However, the timescale of the flux decay in the models (∼seconds) is far too short to explain the gradual decay of flux occurring ∼70 minutes after the main peak of the Megaflare. Longer timescales after beam heating ends in any given kernel are needed to reproduce the timescales of the extended gradual decay phase white-light emission component in megaflares and in lower energy, classical flares (Davenport et al. 2014).

Several heating mechanisms have not yet been critically tested with RHD models and could have an important role in producing extended gradual decay phase flux. We speculate that 3D backwarming from X-rays and/or UV/EUV (300–3000 Å) line emission ("the metal line backheating hypothesis"; Machado et al. 1989; Hawley & Fisher 1992; Fisher et al. 2012) will help to account for the missing flux in the Ca ii K line, the flux in the Balmer lines, and the Balmer decrement for the models. The models of Hawley & Fisher (1992) have shown that irradiation from X-rays can heat the chromosphere and produce a large ratio of emergent line-to-continuum flux. The Fe ii flux in the NUV becomes bright in dMe flares (Hawley & Pettersen 1991; Hawley et al. 2007) and may also contribute significantly to the backwarming. Three-dimensional backwarming from UV/EUV lines has recently been discussed as a source of the slowly decaying component in optical flare kernels (Kawate et al. 2016). Extending the results from 1D modeling of backwarming from Balmer continuum photons (Allred et al. 2006) to 3D may also explain the observations (Metcalf et al. 1990).

After a period of fast retraction from magnetic reconnection, flare loops are expected to continue to slowly contract (Longcope & Guidoni 2011); betatron acceleration of particles in these loops may also contribute as a heating source in the extended gradual phase. Alfvén waves are expected to carry Poynting flux and heat the chromosphere with a delay from reconnection that is inversely proportional to magnetic field strength (Fletcher & Hudson 2008; Reep & Russell 2016). Our line broadening prescription will be incorporated into the publicly available version of the RH code and will be available to test the line decrement and broadening predictions from radiative-hydrodynamic multithread modeling of high-flux density electron (and proton/ion) beam heating combined with each of these additional heating scenarios. The effects of nonthermal collision rates on the Balmer line wings (Canfield & Gayley 1987; Kašparová et al. 2009) may also be revisited with the new prescription for electric pressure broadening in order to constrain the beam fluxes in the newly formed kernels. Modeling of the red-wing asymmetry in solar flares (Ichimoto & Kurokawa 1984; Sharykin & Kosovichev 2014) will also benefit by the new broadening modeling of the hydrogen wings.

We update the starting dMe pre-flare atmosphere from Kowalski et al. (2015b) using the prescription for incident XUEV radiation described in Allred et al. (2015). Charge conservation is calculated using the ionization from the detailed elements hydrogen and helium, but for calcium the electron contribution to charge conservation is calculated from LTE (as for the remaining elements). In future work, we intend to include a Ca i–ii–iii model atom in RADYN that is appropriate for the NLTE ionization and charge contribution for an M dwarf atmosphere (in the Sun, Ca ii is dominant and a calcium ion without the neutral stage is sufficient for the detailed losses in RADYN). For atoms and ions not treated in detail, we use the radiative loss function from Allred et al. (2015), which includes the most updated losses from CHIANTI at $T\gt {\rm{20,000}}\,{\rm{K}}$ but excludes several ions at low temperature at $T\lt {\rm{20,000}}\,{\rm{K}}$ that are optically thick in flares (e.g., Fe ii, Si ii, Mg ii, Al ii) as discussed by Ricchiazzi & Canfield (1983) and Hawley & Fisher (1992). We converge the atmosphere top boundary that is reflecting with a zero temperature gradient. We also apply a small amount of non-radiative heating (0.23 erg cm⁻³ s⁻¹) to all heights with $T\gt 1$ MK for a temperature gradient in the corona that produces a transition from conductive cooling to conductive heating at the transition region (Klimchuk et al. 2008). The atmosphere is relaxed with this condition and this heating is applied through the flare simulations to keep the corona from cooling below the pre-flare temperature. We also note that the chemical equilibrium of hydrogen and H₂ is included using LTE dissociation, which is important in the photosphere. Continuum wavelengths that have been added to the detailed radiative transfer are λ = 1332, 1358, 1389, 1407, 1435, 1610, 2519, 2671, 2780, 2826, 3300, 3500, 3615, 4170, and 6010 Å. These are needed for direct comparisons to data.

We apply electron beam heating for 2.3 s with ${E}_{c}=500\,\mathrm{keV}$, $\delta =7$, and an energy flux density of 2 × 10¹² erg cm⁻² s⁻¹ (2F12). The atmosphere is allowed to relax for 60 s after 2.3 s. These parameters were chosen so that the beam energy is localized to the lower atmosphere and results in a local temperature maximum near $T\sim {\rm{10,000}}\,{\rm{K}}$ within a short time to be consistent with the short pulses observed during solar and dMe flares (Aschwanden et al. 1995; Robinson et al. 1995).

In Figure 12, we show the physical parameters at $t=2.2\,{\rm{s}}$ in the lower atmosphere: the electron density variation, gas velocity, and beam heating in panel (a) the emergent flux spectra at representative times (${F}_{\mathrm{hotspot}}(t);$ see the text) in panel (b), the contribution function to the emergent continuum intensity and the continuum emissivity processes at $\lambda =4170$Å in panel (c) using the method described in Kowalski et al. (2017), and an inset of the Hγ line with the TB09+HM88 broadening calculated from RH in panel (d). The gray line in panel (c) shows the cumulative contribution function (${C}_{I}^{\prime }$) normalized from 0 to 1 on the right axis, which illustrates the depth range over which the emergent intensity is formed. The blue ($\lambda =4170$ Å) continuum forms over $z=150\mbox{--}235\,\mathrm{km}$, which corresponds to a log column mass range of log m/g cm${}^{-2}=-1.26$ to −2.05, a large range of electron density from $(5-40)\times {10}^{14}$ cm⁻³ (with a contribution-function-weighted average of $2.3\times {10}^{15}$ cm⁻³), and a temperature range from $T={\rm{10,800}}\mbox{--}{\rm{12,700}}\,{\rm{K}}$. The electron density weighted by the contribution function for the emergent intensity over the Hγ line ranges from 2 × 10¹³ cm⁻³ at line center to $1.8\times {10}^{15}$ cm⁻³ in the far wing (${\lambda }_{\mathrm{rest}}+20$ Å). The broadening is comparable to the Vega spectrum of Hγ, which forms over an electron density range that is significantly lower (${n}_{e}={10}^{12}$ cm⁻³ at line center and only 5 × 10¹⁴ cm⁻³ in the far wings). Vega has a factor of nearly 100 lower gravity compared to the M dwarf, and thus the radiation is formed over $\sim 1000\mbox{--}2000\,\mathrm{km}$ in the line wing (Figure 4) compared to 85 km in the line wing of the M dwarf ${E}_{c}=500\,\mathrm{keV}$ model. Thus, a comparable opacity is produced over a larger range of heights with lower density in Vega. In the Balmer edge region, the highest order undissolved Balmer line in absorption is ≈H13 in the ${E}_{c}=500\,\mathrm{keV}$ model (Figure 9), whereas in the Vega model (Figure 3 top) it is H16 or higher, as expected from the density dependence on the occupational probability of higher levels (see Figure 9 of Kowalski et al. 2015b).

Zoom In Zoom Out Reset image size

Using the equations in Holman (2012) for the energy lost due to the return current electric field (assuming that proton beams do not neutralize the electron beam), the energy lost over the top 8 Mm of the corona is 4 keV/electron and the Joule heating rate is ∼15 erg cm⁻³ s⁻¹. The 500 keV electrons lose a small fraction of their initial energy and thus we do not expect the beam spectrum to be modified significantly. The coronal heating rate from the return current is 65× the heating rate necessary to keep the pre-flare corona hot. This compares to nearly 8000 erg cm⁻³ s⁻¹ that results from an F13 beam with ${E}_{c}=37\,\mathrm{keV}$ and $\delta =3$. The return current drift speed is <10% of the electron thermal speed in the corona, and therefore we do not expect energy-draining double-layers to develop (Lee et al. 2008; Li et al. 2014); for the F13 beam the drift speed is three times the electron thermal speed and double layers will form.

We apply electron beam heating for 2.3 s with ${E}_{c}=150\,\mathrm{keV}$, $\delta =3$, and an energy flux density of 5 × 10¹² erg cm⁻² s⁻¹ (5F12). The atmosphere is allowed to relax for 60 s after 2.3 s. For this beam, we calculate that $\lt 20\,\mathrm{keV}$ is lost per electron in the propagation over the top 8 Mm of the loop and there is a Joule heating rate from the return current of 350 erg cm⁻³ s⁻¹ that will result in a temperature change. The drift speed of the return current is 40% of the electron thermal speed in the corona, and thus we do not expect energy-draining double-layers to develop (Lee et al. 2008; Li et al. 2014).

The model results are shown in Figure 13 for the same panels as for the E_c = 500 keV model in Appendix A.2. Compared to the phenomenological models of Cram & Woods (1982), the temperature profile, density, and emergent continuum spectrum are similar to the extreme model #5, whereas the broad Balmer wings with deep central reversals are similar to the emergent Hα profile from model #3 that exhibits a low column mass transition region. In the ${E}_{c}=150\,\mathrm{keV}$ model, the transition region is at low column mass as in the Cram & Woods (1982) model #3, but the temperature and electron density values at high column mass (log m/g cm${}^{-2}=-2$) are larger than the Cram & Woods (1982) model #5.

Zoom In Zoom Out Reset image size

Compared to the 2F12, E_c = 500, $\delta =7$ heating model (Appendix A.2), the blue ($\lambda =4170$ Å) continuum forms over comparable heights ($z\sim 150\mbox{--}250\,\mathrm{km}$; see the ${C}_{I}^{\prime }$ curves), column masses (log m/g cm${}^{-2}\gt -2.2$), temperatures ($\sim {\rm{12,000}}\mbox{--}{\rm{13,000}}\,{\rm{K}}$), and electron densities (${n}_{e}=5\mbox{--}40\times {10}^{14}$ cm⁻³ with the same contribution function-weighted electron density of $2.3\times {10}^{15}$ cm⁻³). The relative importance of each emissivity process in the formation of the continuum intensity (panel (c) in Figures 12 and 13) is very similar for the two models, while the contribution function-weighted electron density over the Hγ lines (panel d in Figures 12 and 13) is also similar. In both heating simulations, the formation of the emergent continuum intensity at $\lambda =3615$ Å is shifted 50 km higher than the formation of the $\lambda =4170$ Å emergent continuum intensity because of the larger Balmer ($n=2\to \infty$) bound–free opacity compared to the Paschen ($n=3\to \infty$) bound–free opacity. As a result, the $\lambda =3615$ Å emergent continuum intensity originates from a lower electron density than the $\lambda =4170$ Å continuum intensity. This produces much smaller Balmer jump ratios (${F}_{3615}/{F}_{4170}$) in the model spectra than for a spectrum that results from hydrogen recombination over low optical depth. In contrast, the F13 model spectra exhibit small Balmer jump ratios because of the large variation of the physical depth range as a function of wavelength, as described in Kowalski et al. (2015b, 2016) and Kowalski (2015). The F13 model also has two flaring layers with an electron density profile that increases outward, in contrast to the high E_c models, which exhibit electron density profiles that increase toward lower heights. The electron density variation over the two flaring layers combined with the optical depth variation (resulting in the physical depth range variation) as a function of wavelength both contribute to the characteristics of the emergent flux spectrum in the F13 model.

The striking differences in the emergent flux and intensity spectra of the Balmer lines and at wavelengths in the Balmer continuum for the E_c = 150 and ${E}_{c}=500\,\mathrm{keV}$ models can be explained with the panels in Figures 12 and 13. The lower energy cutoff (${E}_{c}=150$ keV) heats the mid and upper chromosphere to a much larger extent than the ${E}_{c}=500\,\mathrm{keV}$ model; the beam heating peaks at $z\sim 280\,\mathrm{km}$ in the former and at $z\sim 180\,\mathrm{km}$ in the latter (panel a of Figures 12 and 13). The larger heating rates at $z\gt 250\,\mathrm{km}$ in the ${E}_{c}=150\,\mathrm{keV}$ model produce a larger ionization fraction and excitation at these heights, which leads to a lower optical depth in the Balmer continuum and Balmer lines. For example, ${\tau }_{3615}=1$ occurs at z = 210 km in the ${E}_{c}=150\,\mathrm{keV}$ model and at z = 225 km in the ${E}_{c}=500\,\mathrm{keV}$ model. For the emergent Balmer continuum flux (e.g., $\lambda =3615$ Å) in the ${E}_{c}=150\,\mathrm{keV}$ model, the photons are formed deeper and at higher n_e: at ${\tau }_{3615}=1$ the electron density ($1.3\times {10}^{15}$ cm⁻³) is 1.8x larger than the density (7 × 10¹⁴ cm⁻³) at ${\tau }_{3615}=1$ in the ${E}_{c}=500\,\mathrm{keV}$ model, which results in a larger continuum emissivity ($\propto {n}_{e}^{2}$) at the heights where ${\tau }_{3615}=1$ and a larger emergent continuum flux. Therefore, the emergent Balmer continuum flux in the ${E}_{c}=150\,\mathrm{keV}$ model is relatively brighter than the $\lambda =4170$ Å flux (the Balmer continuum flux is "in emission") whereas the Balmer continuum flux in the ${E}_{c}=500\,\mathrm{keV}$ is fainter than the $\lambda =4170$ Å continuum flux (the Balmer continuum flux is "in absorption"). The formation of the $\lambda =4170$ Å continuum flux is similar in the two models because the ${E}_{c}=150\,\mathrm{keV}$ model has a very hard distribution ($\delta =3$) and a higher flux, which heats the deeper layers to a comparable level as the softer, lower flux, higher-energy ${E}_{c}=500\,\mathrm{keV}$ model.

At ${\lambda }_{\mathrm{rest}}+2.5$ Å, the ${E}_{c}=150\,\mathrm{keV}$ model Hγ intensity is in emission while the ${E}_{c}=500\,\mathrm{keV}$ model Hγ intensity is in absorption. Although the electron density for the line formation appears similar in panel (d) of Figures 12 and 13, there is a factor of 1.5× larger electron density for the ${E}_{c}=150\,\mathrm{keV}$ model at ${\lambda }_{\mathrm{rest}}+2.5$ Å. Because of the higher temperature at $z\gt 200\,\mathrm{km}$, the optical depth is lower in the line and Hγ is formed deeper where there is higher electron density. The higher temperature over the line formation for the ${E}_{c}=150\,\mathrm{keV}$ model also contributes to a larger emergent intensity because n_u for Hγ has a larger population density due to the scaling of the collisional rates with temperature via the Boltzmann exponential factor (n_u and n_l for Hγ are near LTE at the height of the maximum in the contribution function in these models).