SOLAR FLARE CHROMOSPHERIC LINE EMISSION: COMPARISON BETWEEN IBIS HIGH-RESOLUTION OBSERVATIONS AND RADIATIVE HYDRODYNAMIC SIMULATIONS

IBIS is a dual Fabry–Perot system, capable of full-Stokes dual beam polarimetry of different spectral lines. The spectral lines and six polarization states are scanned sequentially and reconstructed into images of the Stokes vector during the data reduction process. During our observations, we scanned the chromospheric Hα 6563 Å and Ca ii 8542 Å lines and the photospheric Fe i 6302 Å line. Each line scan (with all polarization states and wavelength points) took about 25 s, giving an overall cadence of 95 s (including overhead such as filter wheel motions).

The IBIS data reduction was performed using our graphical user interface and included corrections for the dark and gain, the alignment of all channels, the wavelength shift due to the collimated Fabry–Perot mount, and the prefilter transmission profile. Speckle-reconstructed images from a simultaneous broadband (white light) camera were used to destretch all spectral channels to remove variations due to seeing. A polarimetric calibration was also performed but is not relevant for our study because we only analyze the intensity profiles. The flare occurred relatively late during the observing day, by which time the atmosphere around the telescope is usually heated up and the seeing is not ideal. Depending on the spectral line, only 3–4 of the 10 scans had stable and favorable seeing to be used for this paper.

The line profiles were averaged along the ribbon. We estimated the area of the footpoints with pixels whose values are above 2500 DN, which is equal to a flux of 31,250 DN s⁻¹. This corresponds to ≈50% of the maximum intensity. The contours of this intensity level are illustrated in Figure 2.

Zoom In Zoom Out Reset image size

2.1.1. Hα Observations

2.1.2. Ca ii 8542 Å Observations

RHESSI had full coverage of the impulsive phase of the flare and part of the decay phase until its sunset, which started at 19:23:52 UT (see Figure 1).

Using the CLEAN algorithm, we reconstructed images at 25–50 keV integrated over 52 s intervals from 19:08:52 to 19:35:44. These images revealed two HXR footpoint sources and a loop top source at later times (see Figure 3). The western footpoint is brighter and moves faster than the eastern one, which is located in a stronger magnetic field. Such asymmetric footpoint emission is understood to be a result of asymmetric magnetic mirroring (e.g., Wang et al. 1995; Jin & Ding 2007; Liu et al. 2009; Yang et al. 2012).

Zoom In Zoom Out Reset image size

2.2.1. Inferring the Non-thermal Electron Distribution

We can infer the non-thermal electron distribution, and thus its energy flux, by analyzing RHESSI X-ray spectra. The moderate count rates of this M3.0 flare allow imaging spectroscopy of spatially resolved sources only during the HXR peak of the impulsive phase. To cover the temporal evolution of the flare, we thus chose to fit the spatially integrated spectra which are dominated by the western HXR footpoint (see Figure 3) that coincides with the flare kernel position observed by IBIS. This provides reasonable diagnostics for the electron energy flux as an input to our RADYN simulation of that kernel.

Because each of RHESSI's nine germanium detectors has slightly different characteristics and makes independent measurements, we analyzed the spectra of individual detectors separately. We then used the means and standard deviations of the fitting parameters to obtain the best-fit parameter and uncertainties. This has several advantages (e.g., avoiding energy smearing) over the conventional approach of directly fitting the average count spectra of all detectors (for details, see, Liu et al. 2008; Milligan & Dennis 2009). We excluded detectors 2 and 7 because of their abnormally high threshold and/or low energy resolution (Smith et al. 2002). We also excluded detector 6 due to its consistently higher χ² values of the fitting results than those of the other detectors. For each of the remaining six detectors, we obtained photon spectra by integrating 30 s intervals from 19:09:30 to 19:13:00 UT and 60 s intervals from 19:13:00 to 19:23:52 UT until the spacecraft night.

Using the standard Object Spectral Executive software package (Brown et al. 2006), we applied corrections for albedo, instrumental emission lines, pulse pileup, and the detector response matrix. We fitted the spectra with a thermal component plus a thick-target, non-thermal component consisting of a broken power law,

$$\begin{eqnarray}&&F(E)=(\delta -1)\frac{{{\mathcal{F}}_{c}}}{{{E}_{c}}}{{\left( \frac{E}{{{E}_{c}}} \right)}^{-\delta }},\end{eqnarray} \tag{ 1 }$$

for energies E > E_c and a power law with a fixed index δ = −2 for E < E_c. Here ${{\mathcal{F}}_{c}}=\int _{{{E}_{c}}}^{\infty }F(E)\ dE$ is the total electron flux above E_c.

Figure 4(a) shows a spectrum of detector 4, as an example, during the impulsive phase, from which we can clearly see the dominance by the thermal component at low energies and by the non-thermal component at high energies. The green and blue lines show the fitted thermal and non-thermal components, respectively, and the red line is the total fit, summing all of the components.

Zoom In Zoom Out Reset image size

The temporal evolution of the power-law index δ and energy cutoff E_c of the non-thermal component are shown in Figures 5(a) and (b). We find a soft–hard–soft spectral evolution with the hardest spectra (lowest δ) occurring at 19:12–19:15 UT.

Zoom In Zoom Out Reset image size

From RHESSI we obtain the total rate of injection of electrons. To estimate the flux of electrons (or energy flux $\mathcal{E}(E)=E\times F(E)$) within the loop, it requires the knowledge of its cross-sectional area. The spatial resolution (≥7'') of the combined RHESSI detectors 3–9 used for HXR imaging is insufficient to resolve the flare footpoints. We thus approximated this area with SDO/AIA 1700 Å (1 2 resolution) kernels within the blue contours at 75% of the maximum intensity, as shown in Figure 6. This approximation is justified by the fact that the primary 1700 Å kernel of the largest area and highest intensity is cospatial with the RHESSI HXR footpoint, as shown in Figure 6. Similar cospatiality between Hα kernels and HXR footpoints has also been reported (e.g., Liu et al. 2007). However, we also note some small 1700 Å kernels without a corresponding HXR source, which could be caused in part by RHESSI's limited dynamic range of the order of 1:10.

Zoom In Zoom Out Reset image size

As the flare evolves, the contours associated with the AIA 1700 Å easternmost footpoint are no longer cospatially aligned with the RHESSI 25–50 keV emission. The western footpoint (observed by IBIS) is the main contributor to the HXR flux and still coincides with the 1700 Å emission. In addition, we also see HXR emission at different locations than 1700 Å emission (see, e.g., Figures 6(c) and (d)). One possibility is that such HXRs are emitted from the loop top source (see Figure 3(f)). Nevertheless, this mismatch could result in an overestimate of the flare loop cross-sectional area, and thus an underestimate of the electron energy flux. We estimated that the uncertainty in the inferred area is about 16%.

Images taken between 19:11 and 19:13 UT contain several saturated pixels which we discarded by imposing an upper threshold of 16000 DN. The resulting areas were then interpolated to the time intervals used for RHESSI spectral fitting and are shown in Figure 5(d). Finally, by dividing the electron power by the estimated loop cross-sectional area, we obtained the electron energy flux, whose temporal evolution is shown in Figure 5(e).

The RADYN code simultaneously solves the equations of hydrodynamics, population conservation, and radiative transfer implicitly for a one-dimensional adaptive grid (Dorfi & Drury 1987), as described by Carlsson & Stein (1992).

For radiative transfer calculations, atoms important to the chromospheric energy balance are treated in non-LTE. These include six-level plus continuum hydrogen, six-level plus continuum, singly ionized calcium, nine-level plus continuum helium, and four-level plus continuum, singly ionized magnesium. Line transitions treated in detail are listed in Table 1 of Abbett & Hawley (1999). A complete redistribution is assumed for all of the lines, except for the Lyman transitions in which partial frequency redistribution is mimicked by truncating the profiles at 10 Doppler widths (Milkey & Mihalas 1973). Other atomic species are included in the calculation as background continua in LTE, using the Uppsalla opacity package of Gustafsson (1973).

The addition of hydrodynamic effects due to gravity, thermal conduction, and compressional viscosity to the original RADYN code was described by Abbett & Hawley (1999). Later additions by Allred et al. (2005) included photoionization heating by high-temperature, soft X-ray emitting plasma, optically thin cooling due to thermal bremsstrahlung and collisionally excited metal transitions, and conductive flux limits to avoid unphysical values in the transition region of large temperature gradients.

We assumed a single quarter circle loop geometry in a plane-parallel model atmosphere, discretized in 191 grid points. The model loop is 10 Mm in height. We assume a symmetric boundary condition at the loop apex (z = 10 Mm). We note from the observations (Figure 6) that the footpoints appear to move during the flare but not more than its diameter so that our assumption of a single loop is a reasonable approximation.

The initial atmosphere (Figure 7) was adopted from the FP2 model of Abbett & Hawley (1999), which is generated by adding a transition region and corona to the model atmosphere of Carlsson & Stein (1997). The temperature was fixed at 10⁶ K at the loop top and no external heating was provided. This allowed the the atmosphere to relax to a hydrodynamic equilibrium state.

Zoom In Zoom Out Reset image size

Initially, the bottom boundary is located in the upper photosphere; the chromosphere is at 0.9 Mm from the bottom of the loop and the transition region is at a distance of 1.56 Mm from the bottom. During the evolution of the atmosphere, we assume open boundaries at the bottom and apex of the loop, extrapolating if necessary.

The non-thermal electron heating was calculated from the power-law spectrum provided by RHESSI spectral fits (see Figure 5). It was included in RADYN as a source of external heating in the equation of internal energy conservation (Q term of Equation (3) in Abbett & Hawley 1999) and has been updated at every integration time interval based on RHESSI data. Intermediate times have been interpolated.

In general, the hydrodynamic evolution of the atmosphere is qualitatively similar to the F09 case reported by Abbett & Hawley (1999). Here, we focus on a narrow region within a distance range of z = 0.62–2.7 Mm, as shown in Figure 8, which covers the upper chromosphere and transition region. This is the region where the dynamic evolution relevant to the Hα and Ca ii 8542 Å line emission occurs.

Zoom In Zoom Out Reset image size

The electrons injected at the loop apex travel downward, losing energy by Coulomb collision, and they heat the chromospheric plasma. This leads to an overpressure which drives chromospheric evaporation with a plasma upflow starting at about 5 s (19:09:33). This also causes the transition region to move upward (see, e.g., the violet dotted line at 19:11:43 in Figure 8).

As time proceeds and the injected non-thermal electron flux decreases, the atmosphere relaxes (see Figure 5). The temperature in the corona and the position of the transition region are directly related to the flux of the injected electrons. This explains the small displacements of the transition region which can be seen over time in Figure 8 in the lines illustrating the spatial distribution at 19:18:28 (green), 19:20:43 (yellow), and 19:23:52 (red).

In order to compare the synthetic Hα and Ca ii 8542 Å profiles with the observations, we integrate the model intensities at each wavelength for 25 and 27 s, respectively, which is the cadence of the observed profiles. A constant microturbulent velocity of 4.5 km s⁻¹ has been applied to compensate for the lack of small-scale random motions in the model (de la Cruz Rodríguez et al. 2012).

Since the atmosphere at the beginning of the run is in a condition of equilibrium, we treat the line profile at this time as a quiet Sun profile, which is subtracted from the intensity profiles at other times, obtaining the so-called excess line profile (Henoux et al. 1998; Matsumoto et al. 2008), which will allow us to better interpret how the flux of electrons affects the flare emission. The observed Hα and Ca ii 8542 Å profiles have a wavelength range of 3.8 and 4 Å which is also used to display the simulated profiles.

In order to better understand how atmospheric evolution affects the chromospheric emission, we write the formal solution of the transfer equation for emerging intensity (Carlsson & Stein 1997):

$$\begin{eqnarray}&&I_{\nu }^{0}=\frac{1}{\mu }{{\int }_{{{\tau }_{\nu }}}}{{S}_{\nu }}{{e}^{-\frac{{{\tau }_{\nu }}}{\mu }}}d{{\tau }_{\nu }}=\frac{1}{\mu }{{\int }_{z}}{{S}_{\nu }}{{\chi }_{\nu }}{{e}^{-\frac{{{\tau }_{\nu }}}{\mu }}}dz=\frac{1}{\mu }{{\int }_{z}}{{C}_{i}}dz,\end{eqnarray} \tag{ 2 }$$

where χ_ν is the monochromatic opacity per unit volume; S_ν is the source function, which is defined as the ratio between the emissivity and the opacity of the atmosphere; τ_ν is the monochromatic optical depth; and the integrand C_i is the so-called intensity contribution function, which represents the emergent intensity emanating from height z.

3.4.1. Evolution of the Hα Line Profile

The top row of Figure 9 shows the synthesized Hα excess profiles, which are asymmetric during the early stage of the flare and become almost symmetric later when the atmosphere relaxes. The line profile including the quiet Sun emission (red solid line in the right column of Figure 9) shows a dip in the line core.

Zoom In Zoom Out Reset image size

By integrating the intensity along the line profile, and subtracting the quiet Sun emission, we estimated the evolution of the Hα excess flux over time. Figure 5(f) shows the Hα light curve where flux has been averaged during the integration time of each RHESSI spectrum. As the plasma is pushed upward, chromospheric evaporation takes place and the Hα flux decreases. After 19:14:24, the atmosphere is more stable and the Hα flux varies with the flux of the injected non-thermal electrons.

Between 19:10 and 19:13, the density population at the energy level n₃ of the hydrogen atom decreases by almost a factor of two with respect to the density population at the energy level n₂ (see Figure 10). Therefore, the ratio n₃/n₂ at this time range decreases. The fact that an Hα photon is emitted by the transition from n₃ to n₂ explains the decrease of the Hα flux at these times. Calcium atoms have a similar behavior, which explains the similar decrease shown in Figure 5(g).

Zoom In Zoom Out Reset image size

The top row of Figure 11 shows the intensity contribution function, C_i (increasing from bright to dark), for Hα at the same times as in Figure 9. The line frequencies are in velocity units, where positive velocities represent plasma moving upward, toward the corona, and negative velocities denote material moving downwards. The blue line represents the atmospheric velocity stratification and the black line shows the line profile (including the quiet Sun emission, as the red profiles of Figure 9). The green line represents the height at which τ_ν = 1, showing us that the height formation of Hα wings is coming from the lower chromosphere (≈0.2 Mm) and is constant over time, while the height formation of the core varies over time during the two minutes, moving toward lower heights.

Zoom In Zoom Out Reset image size

Studying how the contribution function changes in time (Figure 11 top) with height and wavelength, we found that the Hα emission profile becomes broad and centrally reversed not only due to non-thermal effects when the atmosphere is bombarded by energetic electrons (Canfield et al. 1984), but also due to the temperature spatial distribution and the sudden behavior change of the source function in a very thin atmospheric layer.

3.4.2. Evolution of the Ca ii 8542 Å Line Profile

To compare the observed line profiles with the synthetic profiles, we first have to calibrate the spectral lines of IBIS and RADYN to the same reference system. We use the continuum emission of the Fourier Transform Spectrometer (FTS) atlas taken at the McMath–Pierce Telescope (Brault & Neckel 1999) as reference. By doing so, both the observed and synthetic profiles can be calibrated and normalized to the continuum.

In order to calibrate the lines to the continuum of the FTS atlas, we multiplied the quiet Sun line intensity by a factor such that the distance between the line and the FTS atlas is minimal at the continuum (see Figure 12).

Zoom In Zoom Out Reset image size

As mentioned in Section 2.1, the Hα and Ca ii 8542 Å lines observed by IBIS had a wavelength width of 2 and 2.3 Å with respect to the line center, not reaching the continuum. Therefore, we fit the available spectral range to the FTS atlas. The synthetic profiles were properly adjusted to the FTS continuum.

Our synthesized profiles for the quiet Sun are narrower and steeper than the observations, in agreement with Leenaarts et al. (2009). As several authors (Cauzzi et al. 2008; de la Cruz Rodríguez et al. 2011, 2012) explained previously, there are three principal reasons for this discrepancy.

• 1.  
  In general, 1D simulation (even if they take into account the dynamics of the atmosphere) cannot catch all of the structuring and small scales that are present in the chromosphere.
• 2.  
  If the spatial resolution of the simulation is not high enough, the width of the average (spatio-temporal) profile is lower because the hydrodynamic simulations do not contain the necessary small-scale turbulence and the small-scale motions are missing in the model. As Leenaarts (2010) mentions, the increase in grid resolution causes an increase in amplitude of the velocity variations in the middle and upper chromosphere, and hence also causes a widening of the average profile.
• 3.  
  Because of the higher opacities, synthetic lines usually show a much darker line core than the FTS atlas. This may be attributed to a low heating rate in simulations.

Figure 12 shows the resulting IBIS and RADYN lines fit to the FTS atlas after the calibration. The vertical error bars associated to the IBIS profile in Figure 12(a) show the difference between the observed and the synthetic emission during the flare at different wavelength positions. The intensity difference between both Hα lines in the wings is due to the poor fit of the IBIS Hα profile to the continuum, because it is a very broad line and the observations covered only a wavelength range of 4 Å.

After the quiet Sun profiles were calibrated and normalized to the continuum, we applied the normalization factor to the flaring line profiles and in Figure 13 compared them with the IBIS Hα and Ca ii 8542 Å observations Hα and Ca ii 8542 Å at 19:22:40 and 19:22:15 UT, respectively.

Zoom In Zoom Out Reset image size

For a better comparison of the shape of both line profiles, the synthetic profile has been shifted by 0.05 erg s⁻¹ cm⁻² Å⁻¹ in order to align the wings of both lines. The line has been scaled to fit at the core, as shown by the green line of Figure 13.

Figure 13(a) compares the Hα excess line profile obtained from RADYN (red) and observed by IBIS (blue) from 19:22:40 to 19:23:07 UT. As mentioned in Section 4.1, the intensity shift between both profiles at the wings of the line is due to the uncertainty of the continuum fit.

The core of the simulated line is ≈23% less bright than the observation and the wings are slightly narrow, as discussed in Section 4.1. The Hα core is sensitive to the temperature pattern due to the low mass of the hydrogen atom (Leenaarts et al. 2012), which can contribute to the difference in the core emission. The green profile in Figure 13(a) is the synthetic profile scaled to fit at the core and at the continuum of the IBIS profile and the vertical error bar represents the difference between the observed and synthetic emission at different wavelength positions, calculated for this time.

The standard NLTE line formation assumption of statistical equilibrium does not properly fit the observations because of the slow collisional and radiative transition rates when compared to the hydrodynamical timescale. In addition, the Lyα and β lines at least need to be modeled with partial frequency redistribution (PRD) because of their strong influence on the Hα line (Leenaarts 2010). The inclusion of these effects could significantly increase the Hα opacity. Thus, proper modeling of the Hα line requires fully time-dependent radiative transfer with PRD in tandem with the hydrodynamic evolution, a Herculean task that has not been done in 1D hydrodynamic simulations because all 1D simulations have thus far performed radiative transfer assuming complete redistribution (Allred et al. 2005; Kašparová et al. 2009; Varady et al. 2010).

As Leenaarts et al. (2012) explain, PRD effects can be approximated by truncating the Lyman line profiles. RADYN truncates the Lyman line profiles ±64 km s⁻¹ away from the line center frequency. This is a reasonable approximation, with the effect that we obtain higher formation heights than if we assumed PRD.

Figure 13(b) compares IBIS (blue line) and RADYN (red line) in the time range 19:22:15–19:22:42 UT, showing good agreement between both profiles. The green line profile is the synthetic profile scaled to fit the core and the continuum of the IBIS profile.

As mentioned by Smith & Drake (1987), assuming complete spectral redistribution for the scattered photons may be a poor approximation for Ca ii 8542 Å, which could explain the increased intensity in the core. Leenaarts & Carlsson (2009) and Leenaarts et al. (2009) investigated the formation of Ca ii 8542 Å in 3D MHD models that extended into the corona, finding that the 3D effects are important, especially in the core of the line. Considering these two statements and that both lines differ in the core for only ≈ 2.4%, our synthetic Ca ii 8542 Å profile fits the observations very well. As discussed in Section 4.1, the wings area is narrower in the simulated profile because of a lack of small-scale dynamics in the model.

From the study of the monochromatically optical depth of the second panel, τ_ν = 1 (green line in Figure 11), we determine that Ca ii 8542 Å is formed between 0.15 Mm in the wings and 1.05 Mm in the core, while Hα is formed at 0.2 Mm in the wings of the line and 1.15 Mm at the core. Even if both lines had a nearly similar formation height range and the same atmospheric conditions, Ca ii 8542 Å better fits the observations than does Hα. Since the 3 d ²D_3/2,5/2 levels are metastable, they can only be populated from below by collisional excitation, strengthening the sensitivity of the Ca ii infrared triplet to the local temperature. On the other hand, the lower energy level of Hα is 10 eV higher than Ca ii 8542 Å, where the population is very sensitive to the atmospheric parameters, explaining the different behavior of the two lines.