STATISTICALLY DERIVED FLARING CHROMOSPHERIC–CORONAL DENSITY STRUCTURE FROM NON-THERMAL X-RAY OBSERVATIONS OF THE SUN

We have selected all flares in the RHESSI flare list,⁵ between the start of the RHESSI mission on 2002 February 5 and 2010 January 1. Accompanying the flare list are "quicklook images" in different energy bands (3–6, 6–12, 12–25, 25–50, 50–100, and 100–300 keV) typically accumulated over two minutes. Quicklook images in a certain energy band are made only if they were deemed to "reliably image" the X-ray source, that is, that an X-ray source was observed to be roughly at the same location in a majority of RHESSI sub-collimators (SCs). We examined all events for which such images above 25 keV existed. There were 838 of them.

We assume that non-thermal electrons precipitate along magnetic field lines that radially extend from the photosphere to the corona. They propagate from an acceleration region somewhere in the corona, and lose energy and emit bremsstrahlung HXR as per the thick-target model (e.g., Brown et al. 2002). Under these assumptions, the difference in source altitude, Δh_i, between the energy ε_i and a reference energy ε_ref can be used to approximate the average density between these two points (see, e.g., Brown et al. 2002, and Appendix B): the idea is that for an accelerated electron distribution with a negative power-law index propagating toward region of higher densities, non-thermal emission at energy ε spatially peaks where electrons have crossed a column density $N \approx \frac{\varepsilon ^2}{2K}$, with K being a constant. Hence, knowing the distance s between peaks at emissions ε₁ and ε₂, one can get an average density between the two peaks of emission: $n=\frac{N_2-N_1}{s}=\frac{\varepsilon _2^2-\varepsilon _1^2}{2Ks}$. In actuality, $N = \frac{\varepsilon ^2}{2K}$ is only a handy approximation: as discussed in Brown et al. (2002; see, also Xu et al. 2008), there is a dependence on the spectral distribution of injected electrons, and also on the density profile of the medium through which the electron beam is propagating. While both dependences can easily be taken into account when using forward-fitting techniques on singular events, both vary from flare to flare, and we have thus decided to use the $N = \frac{\varepsilon ^2}{2K}$ relationship throughout our statistical study.

We have computed the RHESSI visibilities (Hurford et al. 2002), accumulated over three minutes around the peak HXR flux, in the following energy bands: 6–10, 10–15, 15–20, 20–25, 25–30, 30–35, 35–40, 40–50, 50–60, 60–70, 70–80, 80–90, and 90–100. We have used the hsi_vis_fwdfit.pro routine to find positions of the centroid at different energies, using SCs 3–9. This method basically yields the position of the centroid of the flux (i.e., its center of gravity).

The software also yields error bars on centroid positions. Typically, higher energies, having less count statistics, yield larger error bars. We omitted cases where the software did not converge to a solution.

Once the centroid positions are determined at all energy bands and for all flares, we statistically determine the altitude difference, Δh_i, between emission at all energies ε_i and a reference energy ε_ref. The Δh_i is derived using the same method as in Sato (2006; and described in more mathematical detail in Appendix A). Figure 1 shows plots of R versus ΔR, with R being the projected (on the usual plane of the Sun, perpendicular to the observer's line of sight) distance between the HXR emission and Sun's center, and ΔR the projected distance between the centroid of emissions at different energies. The slope $\frac{\Delta {R}}{R}$ is obtained by linear fitting with a zero intercept coordinate (i.e., going through the origin). Because high energies have lesser statistics, leading to poorer positional accuracy of the centroid, and because there is no potential contamination in the determination of non-thermal centroid position by thermal emission below 35 keV, (although other effects exist and are discussed in Section 2.3), the reference energy was chosen to be 35 keV throughout. Note that there are many outliers, but that their error bars are generally larger, thus they have a lesser contribution to the determination of $\frac{\Delta R}{R}$. Finally, Δh_i = $r_s \frac{\Delta {R}}{R}$, with r_s = solar radius (see Appendix A).

Zoom In Zoom Out Reset image size

From the Δh_i values obtained in the preceding paragraph, the plot in Figure 2 (top) is constructed. It shows the altitude differences between emission at reference energy ε_ref = 35 keV and emission at other energies (we have also looked at other choices of ε_ref, and found similar results, albeit with larger error bars when ε_ref was higher than 35 keV). Also shown are curves derived from subsets of our data: blue, using only events with maximum pile-up error <20% and magenta, using only events with maximum pile-up error >20%. Pile-up issues will be discussed further on. The results presented in Figure 2 (top) show a systematic decrease in altitude as a function of energy, despite somewhat largish error bars at high energies. Note at low energies the departure from a strictly monotonic relationship. We attribute this to the contamination of our sample by thermal emission at low energies.

Zoom In Zoom Out Reset image size

We have checked the reliability of results in Figure 2 (top) by varying certain selection criteria and redoing the procedure. We have tried to select events with quick-look images above 50 keV (as opposed to 25 keV), tried other SC combinations (e.g., 1–9), tried varying accumulation time intervals, different energy bands, and different ε_ref: the curve essentially remains the same though generally with larger error bars. We have also verified that a histogram of the normalized residuals from the $\frac{\Delta {R}}{R}$ linear fitting is normally distributed, with mean 0 and standard deviation 1.

The data presented in Figure 2 (top) will be used in Section 3 to determine an average solar chromopheric/coronal density structure during flaring times. It can already be clearly seen that above ∼50 keV, flare footpoint centroids have a vertical extension of less than ∼0.5 Mm, i.e., electrons above ∼50 keV reach their stopping heights within that 0.5 Mm region. We finally discuss the potentially important sources of contamination to our data set in the next subsection.

• 1.  
  Presence of a thermal component. At low energies (e.g., usually ≲15 keV for M-flares and ≲25 keV for X-flares), the non-thermal thick-target beam model can no longer be applied, due to the presence of high-temperature, high-altitude X-ray-emitting loops.
• 2.  
  Presence of a non-thermal coronal source. (Krucker et al. 2008) These sources are often difficult to observe, due to instrumental dynamic range, but are probably present most of the time (Krucker & Lin 2008). It is likely that such high-altitude sources exist, slightly moving the center of gravity of our footpoint sources to higher projected altitudes. Assuming a coronal source altitude of 10 Mm above the footpoint(s), with 10% of the footpoint non-thermal flux, the upward shift in altitude of the emission centroid could be ∼1 Mm. As coronal source spectra are much steeper than footpoint spectra (e.g., Battaglia & Benz 2006), this effect is less pronounced at higher energies.
• 3.  
  Pulse pile-up. At high photon fluxes, detector pulse pile-up occurs (Smith et al. 2002), where two low-energy photons can be counted as a single photon with the sum of their energies. Hence, an intense thermal loop at ∼18 keV will have a trace in ∼36 keV images (in the case of attenuator state 3). A pile-up error of 10% (as given e.g., by the routine hsi_pileup_check.pro) under attenuator state 3 should roughly shift the center of gravity of 36 keV emission by ∼1 Mm to higher altitudes (assuming that the thermal loop is about 10 Mm above the footpoints). This effect will typically occur at high count rates only, and is very energy dependent, i.e., present mostly at twice the energy where most of the counts are itself dependent on the spacecraft attenuator state: ∼12 keV for attenutator state 0, ∼24 keV for attenutator state 1, and ∼36 keV for attenutator state 3.
• 4.  
  Albedo. Albedo effects (Bai & Ramaty 1978; Kontar & Jeffrey 2010) can also influence the results: up to 20%–40% of the 30–40 keV flux can be due to albedo (assuming isotropic beaming). Assuming that most of the 30–40 keV emission comes from ∼2 Mm above the photosphere, a ∼30% component reflected from the photosphere would shift the centroid of 30–40 keV emission by about ∼0.5 Mm downward. This effect is present for all flares, is expected to be greatest around 30–40 keV, and shifts the centroid to lower altitudes, contrary to the pulse pile-up and coronal source effects. Furthermore, the amplitude of the effect of albedo is maximum near the disk center, and minimum near the limb, and hence has a tendency to offset (in part) the effects of pulse pile-up and coronal sources.

To conclude, it can be said that observed heights for energies ≳50 keV are trustworthy (within their statistical limitations), but that non-thermal emission ≲40 keV may be slightly offset (typically by up to 1–2 Mm) to higher altitudes than in reality, and with the magnitude of this offset varying somewhat with energy.

With these caveats in mind, we have attempted in the next section to derive an average density–height profile from our data set.

Figure 2 (bottom) displays the density structure obtained using all the flares in our list, and using two subsets of it, with error bars. They were obtained using the method described in Brown et al. (2002) and in Appendix B.

A fully ionized corona was assumed, but high-energy electrons (producing most of the high-energy radiation) probably reach regions of low ionization in the lower chromosphere. This means that densities obtained at low altitudes (derived from the higher energies) should actually be multiplied by ∼2.8 (electron beams in a neutral plasma emit less HXR bremsstrahlung; Brown 1973; Kontar et al. 2002).

Figure 2 suggests that (1) emission at lowest energies are clearly contaminated by thermal loops, and results derived from those assuming non-thermal beam models are not to be trusted, and (2) emission above 50 keV is emitted from regions with densities >10¹³ cm⁻³, assuming fully ionized plasma. As we are probably in the chromospheric neutral region, the densities are likely ∼3 times as much, and (3) emission at intermediate energies (∼20–40 keV) appear to come from regions of different densities, depending on the amount of pile-up, itself dependent on the intensity of the non-thermal flux. From (2), one can hypothesize that >50 keV emission always comes from the same altitudes, whether we are dealing with small or large flares. This in turn tells us that emission at ∼35 keV for high pile-up events is almost 1 Mm above that of low-pile-up events (the top plot of Figure 2), which is in accordance with the discussion in Section 2.3, but does not preclude the following scenario: larger (and longer) flares, which typically have high pile-up, are expected to have had greater chromospheric ablation by the time we reach the peak in non-thermal flux (when our observations are made), making the loop denser, and the emission at intermediate energies (∼35 keV) higher. What we observe is an average of this effect over all our selected flares.

Figure 2 (bottom) suggests an atmosphere with density profile with at least two-exponential components. In the next section, we have formulated simple density models and attempted to fit it to the data in Figure 2(top).

We have attempted to fit our data with several density models, but will only discuss a two-exponential with unit step ionization (at altitude h_step: fully ionized for altitudes h>h_step and fully neutral for h < h_step) model, and briefly mention some of the others.

As in Kontar et al. (2008), we have added a data point for our fittings, which actually helps us in determining h_ref, our reference height for emission at ε_ref = 35 keV: it is the well-established hydrogen density at the photosphere n(h = 0) = n₀ = 1.16 × 10¹⁷ cm⁻³.

Zoom In Zoom Out Reset image size

We have used the routine mpfit.pro (found, e.g., in the IDL astronomy library) to make our fittings (and tried also IDL's amoeba.pro routine, with no discernible differences). We have used a Monte Carlo approach to determine error bars for our fitting parameters: for each fitting "run" (a hundred such runs were executed), a random amount was added to each data point. This amount is normally distributed, with mean 0 and standard deviation equal to the nominal error of the data value.

We show only the two-barometric component with a unit step ionization change model (Figure 3) for the low-pile-up case, because we want to minimize the influence of pulse pile-up. Data with higher pile-up and/or other more complicated models produced large χ² results. We only fitted above ∼20 keV. Events with low-pile-up errors have typically little to no thermal emission above this threshold.

Table 1 summarizes the result of the fittings. The h_step parameter is the altitude where the ionization abruptly changes from 100% (corona) to 0% (chromosphere), and ε_step is the minimal initial energy required for electrons to reach the neutral layer. h₀₁ is the altitude of the change from the low-altitude exponential component to high-altitude exponential component. The scale heights H₀ and H₁ are trustworthy, thanks to the additional n₀ data point and to the tight error bars at lower energies. The acceleration altitude h_acc and density n_acc are unsurprisingly poorly constrained. The altitude of ionization change, h_step, is found to be ≈1.3 ± 0.2 Mm above the photosphere, close to the 5%–10% ionization level in the VAL-C or FAL-P atmospheres (Figure 3). The low-altitude scale height was determined to be ∼131 ± 16 km, which corresponds to a (neutral gas) temperature of 5600 ± 700 K, consistent with chromopheric temperatures below the transition region. The high-altitude scale height of the flaring atmosphere is determined to be ∼5.4 Mm, which, assuming a fully ionized isothermal plasma, corresponds to a temperature of ∼115 kK. As this region of the lower corona is very dynamic (i.e., non-isothermal) during flares, and this value is an average over more than 800 flares, we suspect that this usually mid-transition region temperature does not have any intrinsic value.

We would like to point out that events with higher maximum pulse pile-up error required another exponential component at intermediate altitudes in order to produce a "reasonable" fit (albeit with reduced-χ²>4). This led to a fit very close to the CICM and Aschwanden et al. (2002) results at intermediate altitudes (1.5–3 Mm), and, at low and high altitudes, similar to the shown low-pile-up case. While an appealing result, we do not think it is trustworthy given the aforementioned issues.

This simple analytical model will not be used, but is shown as a starting point for other models.

Ionization level is assumed 100% throughout. The density structure is modeled using

$$\begin{equation} n(h) = n_0 e^{-h/H} = n_{\rm acc} e^{-(h-h_{\rm acc})/H} \end{equation} \tag{ C1 }$$

with n₀ = 1.16 × 10¹⁷ cm⁻³, the well-accepted photospheric value, h is the altitude above photosphere, H is the scale height, and h_acc and n_acc are the height and density of the acceleration region.

The column density is

$$\begin{eqnarray} N(h) &=& - \int _{h_{\rm acc}}^h n(h) dh = n_0 H (e^{-h/H} - e^{-h_{\rm acc}/H})\nonumber\\ & = & n_{\rm acc} H (e^{-(h-h_{\rm acc})/H} - 1), \end{eqnarray} \tag{ C2 }$$

where n_acc is the density in the acceleration region.

Using the ε² ≈ 2KN approximation, we get

$$\begin{equation} \varepsilon ^2 = 2 K n_{\rm acc} (e^{-(h-h_{\rm acc})/H} -1) \end{equation} \tag{ C3 }$$

and

$$\begin{equation} h = h(\varepsilon) = -H \ln \left(1 + \frac{\varepsilon ^2}{2KHn_{\rm acc}} \right). \end{equation} \tag{ C4 }$$

The difference in height Δh_ij between the centroid of emission at energies ε_i and ε_j:

$$\begin{equation} \Delta h_{ij} = h_i - h_j = H \ln \left(\frac{1+\frac{\varepsilon _{j}^2}{2KHn_{\rm acc}}}{1+\frac{\varepsilon _{i}^2}{2KHn_{\rm acc}}} \right) \end{equation} \tag{ C5 }$$

The data yield Δh_ij, ε_i, and ε_j, which can be used to determine the fitting parameters H and n_acc. And, finally, h_ref can be determined using Equation (C4) with ε = ε_ref. Note that only good knowledge of H is required to determine h_ref with accuracy, and is independent of what n_acc (or, alternatively, h_acc) is. As we will see, this will prove to be very useful, as n_acc or h_acc will often turn out to have quite large errors.

Ionization level is assumed 100% throughout. The density structure n(h) has h as a variable and H₀, H₁, and h₀₁ as parameters.

$$\begin{eqnarray} n(h)& = &\left\lbrace \begin{array}{@{}ll@{}} n_0 e^{-h/H_0},& h < h_{01}\\ n_{01} e^{-(h-h_{01})/H_1},& h > h_{01} \\ \end{array} \right. \end{eqnarray} \tag{ C6 }$$

with

$$\begin{eqnarray} \ \ \ n_0& = &1.16 \times 10^{17}\,\rm cm^{-3} \end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray} n_{01}& = &n_0 e^{-h_{01}/H_0}\qquad \end{eqnarray} \tag{ C8 }$$

$$\begin{eqnarray} n_{\rm acc}& = &n_{01}e^{-(h_{\rm acc}-h_{01})/H_1}. \end{eqnarray} \tag{ C9 }$$

The column density N(h) requires an additional parameter, h_acc, the height of the acceleration region:

$$\begin{equation} N(h)=- \int _{h_{\rm acc}}^h n(h) dh \end{equation} \tag{ C10 }$$

which can easily be integrated to

$$\begin{eqnarray} N(h)& = &\!\left\lbrace \begin{array}{@{}ll@{}} 0,& h > h_{\rm acc}\\ H_1 n_{\rm acc} (e^{-(h-h_{\rm acc})/H_1}-1),& h_{01} < h < h_{\rm acc}\\ H_1 n_{\rm acc} (e^{-(h_{01}-h_{\rm acc})/H_1}-1)&\\ \quad+ H_0 n_{01} (e^{-(h-h_{01})/H_0}-1),& h < h_{01}\\ \end{array} \right.\qquad \end{eqnarray} \tag{ C11 }$$

The rest of the procedure is similar to as explained in Section C.1, except for the following: (1) we have two additional fitting parameters h₀₁ and H₁, and (2) the explicit ε(h) expression that we get is no longer easily invertible into a simple h(ε) equation, as done in Equations (C3) and (C4). The inversion is hence done numerically, using interpolations.

We use the same equations as in the previous section, except we add another parameter, h_step, the altitude at which the solar atmosphere is assumed to abruptly change from fully ionized to fully neutral (Kontar et al. 2002 and reference therein). In the computation of N(h), the density in regions below h_step is weighted by a factor 1/2.818 (Kontar et al. 2002).

As in Appendices C.2 and C.3, with an additional intermediate exponential component.