ON THE SENSITIVITY OF THE GOES FLARE CLASSIFICATION TO PROPERTIES OF THE ELECTRON BEAM IN THE THICK-TARGET MODEL

One of the outstanding problems in solar physics is to identify the physical mechanism that gives rise to solar flares. The electron beam heating scenario has been a popular mechanism to explain X-ray emissions in solar flares. Under this mechanism, electrons accelerated high in the corona stream down a coronal loop toward the chromosphere, depositing their energy through collisions, and emitting high-energy X-rays via electron–ion bremsstrahlung with the ambient ions (Brown 1971; Emslie 1978). As they deposit their energy, the footpoints heat up and chromospheric ablation drives material into the corona, filling the loop and producing lower energy X-rays via thermal bremsstrahlung (Fisher et al. 1985a, 1985b, 1985c).

Solar flares are commonly classified in terms of their peak soft X-ray (SXR) flux in the 1–8 Å band, as observed by the Geostationary Operational Environmental Satellite (GOES) X-ray flux monitor at terrestrial distance (Garcia 1994). The largest flares, X-class, have peak fluxes greater than 10⁻⁴ W m⁻², with smaller flares classified by a decrease in peak flux by factors of 10 as M-, C-, B-, and A-class. The flux in this GOES passband is determined primarily by a combination of thermal and non-thermal bremsstrahlung and is therefore connected intimately with the properties of the beam of accelerated electrons.

A large number of flares (but not all) show the Neupert effect, which states that the fluence of the hard X-rays (HXRs) is proportional to the flux of the SXRs, or equivalently, that the HXR flux is proportional to the time derivative of the SXR flux (Neupert 1968; Dennis & Zarro 1993; Veronig et al. 2005). In the thick-target model, the deposition of energy by the electron beam induces the observed HXR bursts, subsequently ablating material into the corona and heating the plasma to produce the more gradual SXR light curves. Lee et al. (1995) demonstrate that the Neupert effect, combined with this model of chromospheric ablation, should lead to a linear proportionality between the non-thermal energy and the SXR flux; in contradiction to this, Warren & Antiochos (2004) show that the SXR flux (as seen by GOES) should be proportional to non-thermal energy E^1.75 (1–8 Å) and E^2.24 (0.5–4 Å).

In the past decade the Reuven Ramaty High-Energy Solar Spectroscopic Imager (RHESSI; Lin et al. 2002) and more recently the Fermi Gamma-Ray Burst Monitor (Meegan et al. 2009) have provided excellent coverage of SXR and HXR spectra in flares. The spectra obtained with these instruments can be used to derive the properties of electron beams through various inversion methods (Brown et al. 2003; Kontar et al. 2011). The basic method involves inverting the equation relating an observed intensity to parameters of the electron beam (see Equation (5), below) to solve for the mean electron distribution function (see, e.g., Holman et al. 2003; Kontar et al. 2003; Piana et al. 2003). Brown et al. (2006) evaluate four different methods of inversion, concluding that all of the methods recover the general magnitude of a given distribution but have trouble recovering sharp features.

Electron beams are primarily characterized by three parameters: the energy flux of the electrons, the spectral index of the electron distribution, and the low-energy cut-off of that distribution. The primary hindrance in determining the electron distribution accurately is the low-energy cut-off. At low energies, thermal bremsstrahlung masks the non-thermal component, rendering accurate determination of the cut-off difficult (Holman et al. 2011; Kontar et al. 2011). Further, most studies assume a simple power law with a sharp cut-off for the electron distribution, although it is not clear that this model is correct (Saint-Hilaire & Benz 2005; Holman et al. 2011).

Numerous models have been developed to study heating processes in solar flares. Nagai (1980) developed a one-dimensional hydrodynamic model of loops heated by a thermal conduction front to study the formation of SXR emission in flares. Later models, using the heating function derived by Emslie (1978), explored the hydrodynamic response of the solar atmosphere to a non-thermal beam of electrons (Nagai & Emslie 1984; MacNeice et al. 1984; Mariska et al. 1989). Most of these studies assumed that heating lasted for less than a minute, with fixed beam parameters. However, using results from RHESSI, we can now develop a model combining time-dependent beam properties with an advanced hydrodynamic model to forward model spectra and study the effects of an electron beam on observed spectra directly.

In this paper, we examine the sensitivity of the GOES classification to the beam properties. We use a numerical model to explore the effect of varying each of the beam parameters on light curves as measured by GOES. In Section 2.1 we briefly describe the numerical code used to perform the simulations and several improvements that make it suitable for the studies undertaken here. In Section 2.2 we describe the beam heating function and the specific assumptions we have made in implementing it, and in Section 2.3 we describe the calculation of bremsstrahlung emissions and the synthesis of GOES light curves from the simulations. In Section 3 we present the results of 60 numerical experiments to explore the correlation of GOES classification to each beam parameter derived for two flares observed by RHESSI. Finally, we summarize the results and discuss future work in Section 4.

2.1. The HYDRAD Code

The work presented here has been performed using the HYDRAD code, which solves the hydrodynamic equations for an isolated magnetic flux tube and a multi-fluid plasma (Bradshaw & Cargill 2013). The equations (conservation of mass, momentum, and energy) and assumptions are detailed in the Appendix of Bradshaw & Cargill (2013). However, there have been a number of important improvements to the code.

First, the pre-flare atmosphere is now based on the VAL model C of the photosphere and chromosphere (Vernazza et al. 1981), allowing for more realistic temperature and density profiles in the lower solar atmosphere. We recalculated the density distribution for hydrostatic equilibrium to be consistent with the average particle mass chosen for our model (m_i = 2.171 × 10⁻²⁴ g), which accounts for the relative abundances of hydrogen, helium, and heavier elements (roughly 90% hydrogen and 10% helium). The initial transition region and corona are derived by integrating the hydrostatic equations from the top of the VAL C atmosphere to the apex of the coronal loop. Second, the code has been modified to allow for the presence of neutrals in the lower atmosphere. The two fluids in the code are now electrons and hydrogen atoms (which includes ions and neutrals), with trace amounts of electrons due to ionization of heavier elements. The energy equations have been modified in a manner similar to MacNeice et al. (1984) to include the ionization of hydrogen to account for the effects of neutrals on energy balance between the two fluids (e.g., via collisions and thermal conduction of neutrals; Orrall & Zirker 1961). Finally, the chromospheric radiative energy balance is now based on the recipe derived by Carlsson & Leenaarts (2012), which accounts for cooling from optically thick lines and continua, heating due to the same processes, and heating due to coronal radiation (back-warming from the hot corona). We also used their ionization balance to calculate the electron number density (n_e), for consistency with this radiation calculation. Figure 1 shows the lower atmosphere density and temperature profiles at the start of Run 1.

Zoom In Zoom Out Reset image size

The loop geometry is semi-circular, along the field-aligned direction. The initial temperature and density profiles were found by integrating the hydrostatic equations from the chromosphere (VAL model C; Vernazza et al. 1981) to the apex of the coronal loop (as done, for example, in Aschwanden et al. 2001). The electron and ion populations were assumed to be in thermal equilibrium at the beginning of simulations. The loops are cool and tenuous, before significant heating drives chromospheric ablation (compare, e.g., MacNeice et al. 1984; Nagai & Emslie 1984; Mariska et al. 1989).

2.2. Beam Heating Function

We have implemented a modified form of the beam heating function derived by Emslie (1978). Following Mariska et al. (1989) and Li (1991), we add a low-energy knee to the beam electron distribution while keeping the total beam energy constant (also compare the simulations of Warren 2006). As those authors note, this form of the electron distribution produces smoother temperature and density variations than a sharp cut-off, without significantly altering the X-ray spectrum (see Sections 3.3 and 3.4 of Holman et al. 2011 for comparisons of different shapes of low-energy cut-offs). The distribution function becomes

$\begin{equation} \mathfrak {F}(E_{0}, t) = \frac{4 (\delta - 2) F_{0}(t)}{(\delta + 2) E_{c}^{2} } \left\lbrace \begin{array}{@{}l@{\quad }l@{}} \left(\frac{E_{0}}{E_{c}} \right)^{2} & {\rm if }\; E_{0} \le E_{c} \\ \left(\frac{E_{0}}{E_{c}} \right)^{- \delta } & {\rm if }\; E_{0} \ge E_{c} \end{array}\right., \end{equation} \tag{ 1 }$

where δ is the spectral index above the cut-off energy, E_c, and F₀(t) is the energy flux of the beam as a function of time (erg s⁻¹ cm⁻²). The usual caveats apply to this distribution: the number of electrons accelerated is extremely high, there may be instabilities due to Langmuir wave generation, and there may be a return current generated (see Brown & Melrose 1977 for more information). Following Li (1991), the heat input as a function of position and time is

$\begin{eqnarray} && H_{e}(s, t) = \frac{4 \pi e^{4} n_{{\rm H}} \gamma (\delta - 2)}{\mu _{0} (\delta + 2)} \frac{F_{0}(t)}{E_{c}^{2}} \nonumber\\ &&\!\quad\! \times \left\lbrace \begin{array}{@{}l@{\quad }l@{}}z^{- \delta / 2} B \left[\frac{\delta }{2}, \frac{2}{4 + \beta } \right] &\!\! {\rm if }\; z > 1 \\ z \int _{z}^{1} \frac{ dy }{y^{2} (1 - y)^{2/3}} + z^{- \delta / 2} \int _{0}^{z} y^{\delta / 2 - 1} (1 - y)^{- 2 / 3} dy &\!\! {\rm if }\; z \le 1 \end{array}\right., \nonumber\\ \end{eqnarray} \tag{ 2 }$

where n_H is the hydrogen density, μ₀ is the cosine of the pitch angle of injection, γ = xΛ + (1 − x)Λ', β = [2xΛ + (1 − x)Λ'']/[Λ' + x(Λ − Λ')] (see Appendix A for explanation of the different Coulomb logarithms Λ used here), x is the ionization fraction of hydrogen, z = N/N_c for N the column density and $N_{c} = \mu _{0} E_{c}^{2} / [ (2 + \beta / 2) \gamma 2 \pi e^{4}]$ the stopping depth of an electron with energy E_c, and finally B is the (complete) beta function.

Following the ideas of Hawley & Fisher (1994), we also generalize this expression to non-uniform ionization, which is important for recovering spectral breaks in observed spectra. As noted by those authors, the heating rate is fairly insensitive to the quantity β, but it must be constant to integrate Equation (26) of Emslie (1978). Therefore, we assume that β = 2 (corresponding to a fully ionized plasma, although we do not assume x = 1 in general) and integrate. Second, we modify the column depth to take into account the varying ionization structure (an equivalent ionized column depth, $N^{*}(N) = \int _{0}^{N} ({\gamma }/{\Lambda }) dN^{\prime }$). The heating function in this form is then evaluated numerically as a function of time and position, given the beam parameters.

2.3. Bremsstrahlung Emissions

An important element of this work is predicting the GOES class of our model flares, and so we must evaluate both thermal and non-thermal bremsstrahlung as a function of time. We note the following important properties of GOES relevant in evaluating this emission. (1) GOES observes in the bands 1–8 Å (1.55–12.41 keV) and 0.5–4 Å (3.10–24.82 keV), with flares being classified by the peak flux in the former. (2) GOES light curves are not spatially resolved; i.e., GOES does not distinguish coronal sources from footpoints. (3) GOES light curves are a summation of thermal and non-thermal bremsstrahlung, as well as background emission. We evaluate both thermal and non-thermal bremsstrahlung everywhere along the flaring loop as a function of time and sum the emission to construct a light curve in the native units of GOES (W m⁻²).

To evaluate the thermal bremsstrahlung, we use Equations (1)–(3) of Culhane & Acton (1970), derived for solar flares as observed at terrestrial distance. These equations include contributions from both free–free and free-bound emission. We approximate the emission measure by ∫_Vn²dV ≈ n_en_HAΔs, where A is an estimate of the area based on observations and Δs is the length of a given loop segment. Making these changes, we have, in units of photons s⁻¹ cm⁻² keV⁻¹,

$\begin{eqnarray} I_{{\rm thermal}} &=& 3.6 \times 10^{-39} \overline{Z^{2}} k_{B}^{0.3} T_{e}^{-0.2} \epsilon ^{-1.3} \nonumber\\ && \times \exp {\left (\frac{- \epsilon }{k_{B} T_{e}}\right)} \left [1 - \left (\frac{\epsilon }{88.0}\right)^{k_{B} T_{e} / 3}\right ]^{-1} n_{e} n_{{\rm H}} A \Delta s, \nonumber\\ \end{eqnarray} \tag{ 3 }$

where $\overline{Z^{2}} = 1.4$ is the average charge squared of ions and is the emitted photon energy (in keV). Note that this function is only valid for photon energies 1.5 ⩽ ⩽ 15 keV (Culhane & Acton 1970), which encompasses the full 1–8 Å range, and is only valid for temperatures below about 20 MK. Outside of this range, we use the following expression, derived by combining Equations (1) and (1a) of their paper:

$\begin{equation} I_{{\rm thermal}} = 3.6 \times 10^{-39} \overline{Z^{2}} T_{e}^{-0.5} \epsilon ^{-1.0} \exp {\left (\frac{- \epsilon }{k_{B} T_{e}}\right)} n_{e} n_{{\rm H}} A \Delta s. \end{equation} \tag{ 4 }$

These functions will give the emitted thermal bremsstrahlung at a given spatial location along the loop at a given time. Since GOES has no spatial resolution, the total emission is found by spatially integrating along the loop and over the time period of interest.

Non-thermal bremsstrahlung is more difficult to evaluate. It requires knowledge of the beam electron distribution function and the appropriate cross section of interactions. Following either Holman et al. (2011) or Kontar et al. (2011), thick-target bremsstrahlung observed at terrestrial distance is given in general by (in photons s⁻¹ cm⁻² keV⁻¹)

$\begin{equation} I_{{\rm thick}} = \frac{A}{4 \pi R^{2}} \int _{\epsilon }^{\infty } \mathfrak {F}(E_{0},t) \nu (\epsilon, E_{0}) dE_{0}, \end{equation} \tag{ 5 }$

where R = 1.496 × 10¹³ cm (1 AU), A is the loop area (cm²), $\mathfrak {F}(E_{0},t)$ is the electron flux spectrum of injected electrons (as before, although converted to units of electrons s⁻¹ cm⁻² keV⁻¹), E₀ is the initial electron energy, and ν (below) essentially gives the photon yield for a given electron:

$\begin{equation} \nu (\epsilon, E_{0}) = \int _{E_{0}}^{\epsilon } \frac{ n_{{\rm H}} v Q({\epsilon, E}) dE}{dE/dt} \end{equation} \tag{ 6 }$

for v the electron velocity, Q(, E) the cross section of the interaction (cm² keV⁻¹), and (dE/dt) the energy lost by the electron per unit time (keV s⁻¹). Note that ν has units of photons electron⁻¹ keV⁻¹. Following Holman et al. (2011), we take

$\begin{equation} \frac{dE}{dt} = n_{{\rm H}} v \frac{dE}{dN} = n_{{\rm H}} v \left (\frac{-2 \pi e^{4}}{E} [x \Lambda + (1-x) \Lambda ^{\prime } ] \right), \end{equation} \tag{ 7 }$

where we now evaluate over the column density, allowing for non-uniform ionization (with variables as defined in previous sections). Note that for unit consistency, e⁴ must be in units of cm² keV² with E in keV.

We can reverse the order of integration in Equation (5):

$\begin{equation} I_{{\rm thick}} = \frac{A}{4 \pi R^{2}} \int _{\epsilon }^{\infty } \frac{ n_{{\rm H}} v Q({\epsilon, E}) dE}{dE/dt} \int _{E}^{\infty } \mathfrak {F}(E_{0},t) dE_{0}, \end{equation} \tag{ 8 }$

where the inside integral is now analytic, allowing for easier numerical integration regardless of what form of cross section is chosen. Evaluation of this integral is explained in detail in Appendix B.

We must now decide what cross section is appropriate for non-thermal bremsstrahlung emissions. At the energies we are considering (1.55–12.41 keV and 3.10–24.82 keV), an appropriate choice is the Bethe–Heitler cross section (Bethe & Heitler 1934; Koch & Motz 1959; Johns & Lin 1992a, 1992b), with the Elwert correction factor (Elwert 1939), given by

$\begin{eqnarray} &&{{7}{8.5} Q(\epsilon, E)\,{=}\,\displaystyle\frac{16 \overline{Z^{2}} r_{0}^{2} \alpha }{3} \displaystyle\frac{m_{e}^{2} c^{4}}{\epsilon E (E\,{+}\,2m_{e}c^{2})} \ln {\left [ \displaystyle\frac{1\,{+}\,\left(\displaystyle\frac{(E-\epsilon)(E\,{-}\,\epsilon\,{+}\,2m_{e}c^{2})}{E (E\,{+}\,2m_{e}c^{2})} \right)^{1/2}}{1\,{-}\,\left(\displaystyle\frac{(E\,{-}\,\epsilon)(E\,{-}\,\epsilon\,{+}\,2m_{e}c^{2})}{E (E\,{+}\,2m_{e}c^{2})} \right)^{1/2}} \right]}, } \nonumber\\ &&{{7}{8.5}{\times}\,\displaystyle\frac{[E(E\,{+}\,2m_{e}c^{2})]^{1/2} [E\,{-}\,\epsilon\,{+}\,m_{e}c^{2}] \left [1\,{-}\,\exp {\left ({-}\displaystyle\frac{2 \pi \alpha [E\,{+}\,m_{e}c^{2}]}{[E(E\,{+}\,2m_{e}c^{2})]^{1/2}}\right)} \right ] }{[(E\,{-}\,\epsilon)(E\,{-}\,\epsilon\,{+}\,2m_{e}c^{2})]^{1/2} [E\,{+}\,m_{e}c^{2}] \left [1\,{-}\,\exp{\left ({-}\displaystyle\frac{2 \pi \alpha [E\,{-}\,\epsilon\,{+}\,m_{e}c^{2}]}{[(E\,{-}\,\epsilon)(E\,{-}\,\epsilon\,{+}\,2m_{e}c^{2})]^{1/2}} \right)} \right]} }\nonumber\\ \end{eqnarray} \tag{ 9 }$

with α the fine-structure constant and r₀ the classical electron radius. We can then numerically integrate the thick-target emissions with the use of Gauss–Laguerre and Gauss–Legendre quadrature, as appropriate.

To facilitate comparisons to observational data, the response of the GOES instrument must be accounted for in the bremsstrahlung calculations. Hanser & Sellers (1996) describe the response of the X-ray sensor (XRS) aboard GOES-8 in depth. Note that the sensitivity of the GOES XRS has remained nearly constant for each spacecraft. We convolve the response function (given by Equation (1) and Figure 3 of Hanser & Sellers 1996) with all bremsstrahlung calculations to improve our estimates of GOES classification for the flares studied in this paper.

We performed 60 numerical experiments exploring the correlation of the GOES classification to the beam properties. We study two flares observed with RHESSI and base our simulations on their parameters. The results are split into two sets. The first set is based on beam parameters for the 2002-04-15 M1.2 flare (Table 1) derived by Sui et al. (2005), with a cut-off energy E_c = 24 keV, a spectral index δ(t) ranging between about 6 and 10, and a beam power ranging from approximately 6 × 10²⁵ to 5 × 10²⁷ erg s⁻¹ (see Figure 6 of their paper). The second set is based on the 2002-07-23 X4.8 flare (Table 2) derived by Holman et al. (2003), with a cut-off E_c(t) ranging from about 18 to 43 keV, spectral index δ(t) ranging from about 2.5 to 8, and a beam power between 7 × 10²⁶ and 2 × 10²⁹ erg s⁻¹ (see Figure 3 of their paper). In both cases, the beam parameters were derived from RHESSI observations of the flares. The geometry of the 2002-04-15 flare loop (that is, length and cross-sectional area) is based on the estimate of Veronig & Brown (2004). The length of the 2002-07-23 flare loop is estimated from the footpoint locations in Emslie et al. (2003); the area is estimated from Holman et al. (2003), decreased slightly to 7.0 × 10¹⁸ cm² to calibrate the base run (31) to the observed GOES classification.

Table 1. The Results of Simulations for the 2002-04-15 M1.2 Flare

Run No.	Flux	Δδ	ΔEc	GOES Class	GOES	Tmax	napex, max
			(keV)	(1–8 Å)	(0.5–4 Å)	(MK)	(cm−3)
1	100%	+0	+0	M1.1	C3.1	27.3	6.3 × 1010
2	10%	+0	+0	B5.1	B3.4	5.15	2.1 × 109
3	20%	+0	+0	B9.8	B6.7	8.20	5.6 × 109
4	33%	+0	+0	C1.6	C1.1	12.2	1.4 × 1010
5	50%	+0	+0	C2.5	C1.7	17.2	2.7 × 1010
6	200%	+0	+0	M4.5	M1.1	36.6	1.1 × 1011
7	300%	+0	+0	M9.5	M2.5	42.6	1.5 × 1011
8	500%	+0	+0	X2.4	M7.1	50.1	2.1 × 1011
9	1000%	+0	+0	X10	X3.4	62.3	3.9 × 1011
10	100%	−3	+0	C5.6	C4.8	22.2	4.6 × 1010
11	100%	−2	+0	C7.4	C3.8	24.5	5.5 × 1010
12	100%	−1	+0	C9.5	C3.3	26.3	6.1 × 1010
13	100%	$-\frac{1}{2}$	+0	C9.1	C3.3	26.9	6.0 × 1010
14	100%	$+\frac{1}{2}$	+0	M1.1	C3.1	27.9	6.4 × 1010
15	100%	+1	+0	M1.2	C3.2	28.0	6.5 × 1010
16	100%	+2	+0	M1.0	C3.0	27.9	6.3 × 1010
17	100%	+3	+0	M1.1	C2.9	28.1	6.3 × 1010
18	100%	+0	−20	C7.8	C1.3	32.2	5.7 × 1010
19	100%	+0	−15	M1.2	C2.5	32.5	6.8 × 1010
20	100%	+0	−10	M1.4	C2.9	32.1	7.1 × 1010
21	100%	+0	−5	M1.3	C3.1	30.3	6.9 × 1010
22	100%	+0	+5	C5.3	C3.8	22.5	4.8 × 1010
23	100%	+0	+10	C5.7	C4.6	17.0	3.2 × 1010
24	100%	+0	+15	C6.0	C5.1	13.2	1.6 × 1010
25	100%	+0	+20	C6.3	C5.6	10.8	1.1 × 1010
26	100%	+0	+30	C6.9	C6.5	8.13	5.8 × 109
27	100%	+0	+40	C7.8	C7.6	6.51	3.4 × 109
28	100%	+0	+50	C7.8	C7.9	5.45	2.5 × 109
29	100%	+0	+75	C8.2	C8.6	3.90	1.2 × 109
30	100%	+0	+100	C8.3	C9.1	3.09	9.0 × 108

Note. Run 1 uses beam parameters taken from Sui et al. (2005), while each other simulation modifies one beam parameter, as noted.

Download table as:  ASCIITypeset image

Table 2. The Results of Simulations for the 2002-07-23 X4.8 Flare

Run No.	Flux	Δδ	ΔEc	GOES Class	GOES	Tmax	napex, max
			(keV)	(1–8 Å)	(0.5–4 Å)	(MK)	(cm−3)
31	100%	+0	+0	X3.6	X1.1	28.0	1.4 × 1011
32	10%	+0	+0	M1.2	C6.9	6.44	6.0 × 109
33	20%	+0	+0	M2.3	M1.3	11.3	1.8 × 1010
34	33%	+0	+0	M4.7	M2.1	17.3	4.7 × 1010
35	50%	+0	+0	M9.3	M3.5	21.7	7.5 × 1010
36	200%	+0	+0	X14	X5.1	35.7	2.3 × 1011
37	300%	+0	+0	X32	X8.4	40.6	3.2 × 1011
38	500%	+0	+0	X80	X21	47.3	4.7 × 1011
39	1000%	+0	+0	X300	X89	58.3	8.4 × 1011
40	100%	−3	+0	X2.0	X1.2	23.8	8.1 × 1010
41	100%	−2	+0	X2.7	X1.1	26.4	1.1 × 1011
42	100%	−1	+0	X3.2	X1.1	27.3	1.2 × 1011
43	100%	$-\frac{1}{2}$	+0	X3.4	X1.0	27.6	1.3 × 1011
44	100%	$+\frac{1}{2}$	+0	X3.6	X1.1	28.2	1.4 × 1011
45	100%	+1	+0	X3.5	X1.0	28.2	1.4 × 1011
46	100%	+2	+0	X3.8	X1.0	28.5	1.5 × 1011
47	100%	+3	+0	X3.7	X1.0	28.5	1.4 × 1011
48	100%	+0	−20	M6.1	M2.7	18.1	6.4 × 1011
49	100%	+0	−15	X2.6	M5.3	27.2	1.2 × 1011
50	100%	+0	−10	X4.0	M8.7	29.9	1.5 × 1011
51	100%	+0	−5	X3.9	M9.9	29.5	1.5 × 1011
52	100%	+0	+5	X2.8	X1.0	26.5	1.2 × 1011
53	100%	+0	+10	X1.5	X1.1	28.2	4.1 × 1010
54	100%	+0	+15	X1.6	X1.2	15.7	4.3 × 1010
55	100%	+0	+20	X1.7	X1.4	11.9	2.4 × 1010
56	100%	+0	+30	X1.8	X1.7	8.26	1.1 × 1010
57	100%	+0	+40	X2.0	X1.9	6.18	6.4 × 109
58	100%	+0	+50	X2.2	X2.1	5.04	4.3 × 109
59	100%	+0	+75	X2.3	X2.4	3.44	2.2 × 109
60	100%	+0	+100	X2.4	X2.6	2.64	1.5 × 109

Note. Run 31 uses beam parameters taken from Holman et al. (2003), while each other simulation modifies one beam parameter, as noted.

Download table as:  ASCIITypeset image

In each group of simulations, we vary one beam parameter per simulation in the following way (holding all other values constant). First, we multiply the beam flux F₀(t) by a factor of [(1/10), (1/5), (1/3), (1/2), 2, 3, 5, 10]. Second, we change the spectral index of the electron distribution δ(t) by [−3, −2, −1, −(1/2), (1/2), 1, 2, 3] (although we limit δ to 2.01 since it must be strictly greater than 2). Finally, we change the cut-off energy E_c(t) by [−20, −15, −10, − 5, 5, 10, 15, 20, 30, 40, 50, 75, 100] keV (limiting E_c to a minimum of 1 keV).

Table 1 summarizes the results for the set of simulations concerning the 2002-04-15 flare. Table 2 summarizes the results for the set of simulations concerning the 2002-07-23 flare. In both cases, we list the run number along with the changes in the beam parameters, the GOES class (in both channels), the maximum temperature (in MK), and the maximum apex density (cm⁻³).

For example, consider Figure 2, which shows the results of Run 1. The beam lasts 640 s, for a loop of length 2L = 77 Mm. The upper panels show the electron temperature (reaching a peak of 27.3 MK) and density (with a maximum apex value near 6 × 10¹⁰ cm⁻³) as functions of position, at different times as indicated. The temperature rises for around 300 s, when the beam begins to weaken. There is significant radiative cooling in the chromosphere; in the corona, thermal conduction drives energy losses initially, gradually transitioning to radiative cooling as the density rises and temperature falls. The density in the corona rises due to chromospheric ablation throughout the duration of heating, and the maximum intensity occurs when the heating ceases (a few minutes after it reaches its maximum temperature). The center left panel similarly shows the heat deposition (which includes background heating). As the loop fills due to chromospheric ablation, the heat deposition occurs higher and higher in the corona as the mean-free path of beam electrons decreases. Although the beam continues for 640 s, the beam flux begins to fall around 300 s, and thus the heat deposition begins to fall. The center right figure shows the spectrum calculated at one of the footpoints at a few selected times. Initially (t = 150 s), the emission is entirely non-thermal bremsstrahlung, but as the loop begins to heat and fill, thermal bremsstrahlung eventually becomes the dominant component of the emission at lower energies (e.g., at 600 s, the thermal emission at 1 keV is more than 100 times stronger than the non-thermal emission). Finally, the predicted GOES light curve (with 1–8 Å in red and 0.5–4 Å in blue) is shown at bottom right, with the relative amounts of thermal and non-thermal bremsstrahlung indicated. The light curve (1–8 Å) peaks at 1.1 × 10⁻⁵ W m⁻², corresponding to M1.1, while the 0.5–4 Å component peaks at 3.1 × 10⁻⁶ W m⁻² (C3.1). The actual (background subtracted) GOES light curves, starting at 2002 April 15 23:07 UT, are overlaid for comparison (peaking at M1.1 and C2.0, respectively). The light curve is initially non-thermal bremsstrahlung, but as the loop heats and fills up, thermal bremsstrahlung increases until it becomes the dominant component in the passband.

Zoom In Zoom Out Reset image size

Similarly, Figure 3 shows the evolution of Run 31. The beam lasts for 1220 s, for a loop of length of 2L = 36 Mm. The predicted light curve in this case peaks at 3.6 × 10⁻⁴ W m⁻² (X3.6) in the 1–8 Å channel and 1.1 × 10⁻⁴ W m⁻² (X1.1) in the 0.5–4 Å channel. As before, the temperature rises sharply due to the energy deposition, which then triggers chromospheric ablation. As the loop fills and heats up, the contribution of thermal bremsstrahlung to the light curve rises sharply. In this case, the loop begins to cool and drain before the beam ceases. In both cases, the GOES class is in approximate agreement with the observations. Note that the results are fairly insensitive to the initial density profile. If we had assumed an initial coronal density of 10¹⁰ cm⁻³, over a loop length of 16 Mm, we would find that only electrons with energy ≲ 12 keV would be stopped in the corona, and thus the majority of energy would still be deposited in the chromosphere (using the estimate of stopping depth from Nagai & Emslie 1984). We now turn our attention toward altering the beam parameters one at a time, to examine their effect on the GOES class.

Zoom In Zoom Out Reset image size

First, holding all other values constant, we multiply the beam flux by [(1/10), (1/5), (1/3), (1/2), 2, 3, 5, 10] for both flares (equivalent to multiplying the total non-thermal energy by the same amount). There is a clear correlation between the GOES class and the beam flux (see the first plot in Figure 4). For Runs 1–9, we find that the GOES class Ψ (in the 1–8 Å channel) is proportional to the beam energy E, such that Ψ∝E^α for α = 1.71  ±  0.084; for Runs 31–39, we find α = 1.77  ±  0.060. In the 0.5–4 Å channel, we find a similar trend: for Runs 1–9, we find Ψ∝E^α for α = 1.47  ±  0.103, and for Runs 31–39, α = 1.60  ±  0.077. In Runs 2–5 and 32–33 (and 34 in the 0.5–4 Å channel), the loops were not heated enough to produce significant thermal bremsstrahlung, and thus their GOES classes were primarily determined by non-thermal emissions. Note that non-thermal emissions are linearly proportional to the beam flux (Equation (5)); if we fit a line to just Runs 2–5, we find an exponent α = 0.98 ± .023 in the 1–8 Å channel and α = 1.00 ± 0.014 in the 0.5–4 Å channel, confirming that the non-thermal emission is proportional to the mean beam energy. In all the other runs, the peak of the light curves was a combination of thermal and non-thermal emissions, with thermal emissions being indirectly related to the beam flux. If we remove the simulations without thermal emissions and recalculate, we find exponents α = 1.94 ± 0.026 (1–8 Å) and α = 2.04  ±  0.054 (0.5–4 Å) for Runs 1 and 6–9, and 1.91 ± .014 (1–8 Å) for Runs 31 and 34–39 and 1.84 ± 0.050 (0.5–4 Å) for Runs 31 and 35–39. Thus, the maximum GOES class Ψ is related to the total beam energy E, Ψ∝E^α for α ≈ 1.7 (1–8 Å) and α ≈ 1.6 (0.5–4 Å). There are two reasons why the values differ slightly for each flare: different loop lengths and different cut-off energies, both of which affect the amount of heating and thus thermal bremsstrahlung produced.

Zoom In Zoom Out Reset image size

This result should be compared to that of Warren & Antiochos (2004), who showed that Ψ∝E^1.75 (1–8 Å) and Ψ∝E^2.24 (0.5–4 Å) using simple analytic considerations. Our results are in agreement for the 1–8 Å channel, while they are in contrast in the higher energy channel. If the emission were entirely non-thermal, then we would find a linear relationship (Ψ∝E, Equation (5)), regardless of which channel is being considered. Since the 0.5–4 Å channel is more sensitive to higher energies, we would expect it to be closer to linear than the lower energy channel as it would be more sensitive to the non-thermal emissions, and thus the exponent α should be lower. The reason for this contrast is perhaps their use of the hydrostatic scaling laws derived for the corona, whereas the non-thermal emissions will be primarily at chromospheric depths. It should also be noted that our result, as well as that of Warren & Antiochos (2004), is in contrast to what is expected from the Neupert effect. As shown by Lee et al. (1993, 1995) by using the Neupert effect relation, the maximum flux in the SXRs should be linearly proportional to the energy deposited by the electron beam. Warren & Antiochos (2004) suggest that this difference is due to the peak SXR flux being dependent on other flare parameters, to which we now turn our attention.

Next, holding all the other values constant again, we add [−3, −2, −1, −(1/2), (1/2), 1, 2, 3] to the spectral index δ(t). Varying the spectral index of the electron distribution has only a small effect on the GOES class (see the second plot in Figure 4). For Runs 1 and 10–17, with a change of δ ± 3, the GOES flux (1–8 Å) varies by a factor of two, and about a factor of three in the higher energy channel; for Runs 31 and 40–47, the flux varies by less than a factor of two in both channels. In all cases, though, an increase in the spectral index does increase the GOES class slightly in the 1–8 Å channel and decreases it slightly in the 0.5–4 Å channel. The reason is that an increase in the index reduces the proportion of high- to low-energy electrons, so that more energy is deposited higher up in the atmosphere, increasing the temperature and thus the amount of thermal bremsstrahlung, while slightly reducing the non-thermal emissions. Thus, the channel that is more sensitive to thermal emissions increases slightly, and the channel more sensitive to non-thermal emissions decreases slightly. The immediate question is whether this result is due to the assumed form of the beam electron distribution (i.e., the importance of the low-energy knee electrons). Runs 1, 10–17, 31, and 40–47 were repeated using a sharp low-energy cut-off ($\mathfrak {F}_{0}(E_{0},t) = 0$ for E₀ < E_c), but equivalent amount of total beam energy (see Figure 5). Although the slopes are larger in this case, they remain small. The conclusion is unaltered: the spectral index only marginally affects the GOES class.

Zoom In Zoom Out Reset image size

The primary reason that the spectral index minimally affects the GOES flux is because the total beam energy was held constant. Changing the spectral index is equivalent to slightly altering the proportion of high- to low-energy electrons, which in turn affects the mean location of energy deposition. Increasing the index decreases the number of high-energy electrons in the beam, which then causes deposition of energy slightly higher up the loop (and vice versa). Compare Dennis (1985), who found that the photon spectral index γ is not correlated to the peak HXR flux (and in turn, δ = γ + 1, so that the relation holds for the electron index as well). Also consider the results of Falewicz et al. (2009), who found a small positive correlation between the spectral index and the 1–8 Å flux, similar to that of the present work.

The cut-off energy of the electron distribution can alter the GOES class significantly (see the third plot in Figure 4). Slight decreases in the cut-off increase the GOES flux, but large decreases will significantly decrease the GOES flux. On the other side, as the cut-off increases more and more, the energy is deposited deeper in the dense lower atmosphere, where there is a greater heat capacity and the energy is radiated away more quickly, leading to less heating, less thermal bremsstrahlung, and thus a lower GOES class in the 1–8 Å channel. In the 0.5–4 Å channel, increases in the cut-off increase the amount of non-thermal bremsstrahlung, so the flux in this channel tends to increase with higher cut-off energies. At extremely high cut-offs, the energy is deposited too low to heat the loop, so no thermal bremsstrahlung will be emitted; however, there will still be non-thermal emissions as the beam traverses the loop, which will be determined primarily by the energy flux in the beam. Eventually, as the cut-off is increased to large enough values, the emission in both channels will be entirely non-thermal bremsstrahlung, to which the higher energy channel is more sensitive. As the cut-off decreases, more heat is deposited higher up, leading to a higher temperature and more thermal emission. If it decreases too much, though, the beam is then essentially composed of thermal electrons, leading to less non-thermal bremsstrahlung emission and a lower GOES class (in both channels). The cut-off energy is directly related to the location of energy deposition of the beam electrons. Nagai & Emslie (1984) give the mean stopping column density of an electron with energy E_c as ≈10¹⁷[E_c (keV)]² cm⁻², which can then be used to find the location of maximal energy deposition. In Figure 6, we show the predicted and actual location of maximal energy deposition as a function of time for Runs 1 and 20 in one half of the loop. In Run 20, it should be noted, the coronal density reached higher values and so the deposition location is in general higher than in Run 1, leading to more heating and more thermal bremsstrahlung, and thus a higher GOES class.

Zoom In Zoom Out Reset image size

There is also a clear relation between the maximum temperature of a flare and the maximum GOES flux. Figure 7 shows the maximum electron temperature versus the maximum GOES flux found in the 60 simulations. There is a clear correlation between the two, although there are two separate trends for each group of 30 simulations. This is due to differences in the spectral index. The index was much higher at all times in the M1.2 flare, and thus more energy was deposited higher in the atmosphere, leading to consistently higher temperatures despite lower GOES classes. The spectral index therefore determines whether a flare will be thermally driven or beam dominated. Note, however, that in each case there is a horizontal branch extending to lower temperatures with roughly constant GOES flux. These branches are the cases where there are extremely high cut-off energies, which deposit their energy too low to heat the loop, but still produce roughly constant non-thermal emissions.

Zoom In Zoom Out Reset image size

These results should be compared to Figure 6 and Equation (6) of Feldman et al. (1996), who found a similar trend, but significantly lower temperatures. These lower temperatures result from the fact that they were measuring the temperature at the time of peak emission of the Fe xxv or Ca xix channels of the Bent Crystal Spectrometer on Solar Maximum Mission, which occurs after the time of maximum temperature (the density continues to rise). For example, in Run 1, the maximum temperature was 260 s into the flare, while the maximum intensity of the Ca xix line occurred 500 s into the flare, a 4 minute discrepancy, by which time the loop had cooled by 7 MK (see Bradshaw & Klimchuk (2011) for explanation of the forward modeling of spectral lines such as the Ca xix line).

There is some observational evidence for cold flares, which have detectable GOES emission but low temperatures. For example, Fleishman et al. (2011) report a C-class flare with a maximum temperature of 6 MK. The authors suggest that the low temperature is due to a weak beam flux, which leads to little heating, so that the GOES emission would then be composed primarily of non-thermal bremsstrahlung. This explanation agrees with our results, where small beam fluxes produce low temperatures, but emission detectable by GOES (see the lower left of Figure 7). One other possibility, which seems less likely, is that the flare had a large cut-off energy so that the energy was deposited deep in the chromosphere, once again leading to minimal heating but some emission in the GOES passbands nonetheless (see the horizontal branches in Figure 7).

We have investigated the effect of individual beam heating parameters on the GOES classification of solar flares. Using two independent sets of observationally determined parameters, 60 numerical experiments have been performed, from which GOES light curves have been synthesized. We have found clear trends between the beam parameters and the GOES flux. First, the GOES classification strongly depends on the total beam energy (Ψ∝E^α with α around 1.7 in the 1–8 Å channel and around 1.6 in the 0.5–4 Å channel). This parameter dominates: the SXR spectrum is primarily determined by the amount of heat deposited in the solar atmosphere. Second, the spectral index of the beam electron distribution does not significantly alter the GOES class (whether the cut-off is sharp or a knee). There is a small positive correlation, though: an increase in δ slightly increases the GOES class. Finally, the cut-off energy, which determines the mean location of energy deposition of the beam, affects the amount of both the thermal and non-thermal bremsstrahlung produced. Thus, changes to the cut-off can either increase or decrease the GOES flux, but it is not a simple relation, unlike the other two parameters.

As noted earlier, the Neupert effect predicts a linear relation between non-thermal energy and SXR flux (Ψ∝E), which is in disagreement with the present results and those of Warren & Antiochos (2004). In this paper we have shown that the GOES flux strongly depends on the cut-off energy of the electron beam as well. A change in the cut-off energy can alter the relation between non-thermal energy and SXR flux. Consider the case of an extremely high cut-off energy where there is little heating and thus little thermal bremsstrahlung: all the emission will be non-thermal, and thus a linear relation would follow (Equation (5)). Although it seems unlikely that such extremely high cut-off energies occur, it is clear that the cut-off energy can affect the relation.

We can constrain the parameters as well. For example, Runs 38 and 39 have a GOES class larger than any flare ever observed, and thus can be discarded as unlikely to occur, setting an effective upper limit on the non-thermal energy budget in flares (the largest ever seen with GOES was around X40; Brodrick et al. 2005). The total non-thermal energy in the 2002-07-23 flare was 2.6 × 10³¹ erg (Holman et al. 2003), so there is an effective upper limit around ≈7.5 × 10³¹ erg (for similar cut-off energies). The extremely high cut-off energies used in some of the runs were also unrealistic. For example, Runs 55–60 are X-class flares with maximum temperatures beneath 10 MK and densities too low to be observed, suggesting that such high cut-offs are not realistic (at least for the entire duration of the flare, there is some evidence that they may reach high cut-offs temporarily; Warmuth et al. 2009). The SXR spectrum is only one facet of emissions produced by solar flares, though.

In future work, we will predict extreme-ultraviolet (EUV) spectra to further constrain the model. The properties of EUV lines, such as widths, Doppler shifts, and intensities, can be predicted from the model and directly compared with observations to help pin down signatures of the fundamental driving mechanism of solar flares. For example, Tsurutani et al. (2005) point out that while the 2003 November 4 flare had the largest GOES flux ever measured, the 2003 October 28 flare was more intense in the EUV. The model developed here needs to be used to explain these differences as well.

One of the authors (J.W.R.) was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program—Grant NNX11AQ54H. We thank the anonymous referee for detailed comments that substantially improved the paper. We also thank our colleagues at the SHINE and SPD conferences this summer for their helpful discussions and suggestions.

In the derivation of a heating function due to a beam of charged particles interacting with a hydrogen target, Emslie (1978) defines three effective Coulomb logarithms (the following analysis is based on Ricchiazzi 1982). In the special case that the beam particles are taken to be electrons, the first (due to interactions of beam electrons and protons) is given by (Equation (8) of Emslie 1978)

$\begin{equation} \Lambda = \ln { \frac{m_{0} v^{3}}{\nu e^{2}}} \approx \ln { \frac{m_{e} v^{3}}{\nu e^{2}}} \approx \ln { \frac{\pi ^{1/2} m_{e}^{3/2} v^{3}}{n_{e}^{1/2} e^{3}}}, \end{equation} \tag{ A1 }$

where m_e is the electron mass, m₀ = (m_em_i/m_e + m_i) ≈ m_e is the reduced mass of the system, v is the electron velocity, e is the electron charge, and ν is the plasma frequency (=(n_ee²/πm_e)^1/2 for an electron number density n_e). Taking the average electron energy to be 〈E〉 = (1/2)m_ev², we can rewrite this as

$\begin{eqnarray} \Lambda &=& \ln { \left [ 2 \langle E \rangle \left (\frac{2 \langle E \rangle }{m_{e}}\right)^{1/2} \left (\frac{\pi m_{e}}{n_{e} e^{6}}\right)^{1/2} \right ]} \nonumber \\ &=& \ln { \frac{ 2^{3/2} \pi ^{1/2}}{e^{3}}} + 1.5 \ln {\langle E \rangle } - 0.5 \ln {n_{e}} \nonumber \\ &\approx & 66.0 + 1.5 \ln {\langle E \rangle } - 0.5 \ln {n_{e}}. \end{eqnarray} \tag{ A2 }$

Note that 〈E〉 must be in units of erg here.

The second logarithm, due to interactions between beam electrons and neutral hydrogen, is given by (Equation (12) of Emslie 1978)

$\begin{eqnarray} \Lambda ^{\prime } &=& \ln { \frac{m_{e} v^{2}}{1.105 \chi _{{\rm H}}}} \nonumber \\ &=& \ln { \frac{2 \langle E \rangle }{1.105 \chi _{{\rm H}}}} \nonumber \\ &=& \ln { \frac{2}{1.105 \chi _{{\rm H}}}} + \ln {\langle E \rangle } \nonumber \\ &\approx & 25.1 + \ln {\langle E \rangle }, \end{eqnarray} \tag{ A3 }$

where χ_H = 13.6 eV is the ionization energy of hydrogen.

Finally, the third effective Coulomb logarithm is (Equation (20) of Emslie 1978)

$\begin{eqnarray} \Lambda ^{\prime \prime } &=& \ln { \frac{v}{c \alpha } } \nonumber \\ &=& \ln {\left[ \left(\frac{2 \langle E \rangle }{m_{e}} \right)^{1/2} \frac{1}{c \alpha } \right]} \nonumber \\ &=& \ln { \frac{2^{1/2}}{m_{e}^{1/2} c \alpha } } + 0.5 \ln {\langle E \rangle } \nonumber \\ &\approx & 12.3 + 0.5 \ln {\langle E \rangle }, \end{eqnarray} \tag{ A4 }$

where α is the fine-structure constant. The mean value of the energy can be found from the properties of the beam injection:

$\begin{eqnarray} \langle E \rangle &=& \frac{\int _{0}^{\infty } E \mathfrak {F}(E) dE}{\int _{0}^{\infty } \mathfrak {F}(E) dE} = \frac{\int _{0}^{E_{c}} K \frac{E^{3}}{E_{c}^{2}} dE + \int _{E_{c}}^{\infty } K \frac{E^{1 - \delta }}{E_{c}^{-\delta }} dE}{\int _{0}^{E_{c}} K \frac{E^{2}}{E_{c}^{2}} dE + \int _{E_{c}}^{\infty } K \frac{E^{- \delta }}{E_{c}^{-\delta }} dE} \nonumber \\ &=& \frac{ \frac{E_{c}^{2}}{4} - \frac{E_{c}^{2}}{2 - \delta }}{\frac{E_{c}}{3} - \frac{E_{c}}{1 - \delta }} = \frac{3}{4}\frac{\delta - 1 }{\delta - 2} E_{c}, \end{eqnarray} \tag{ A5 }$

with the condition that δ > 2.

Thick-target bremsstrahlung is evaluated with a double integral, Equation (5), where the inner integral ν(, E₀) gives the photon yield for a single electron of initial energy E₀, current energy E, and losing energy at a rate (dE/dt) through collisions with the ambient medium. The entire electron distribution is then integrated to give the photon yield for all of the electrons in the beam, as a function of time and position. However, the form of Equation (5), and in particular the form of the cross section Q(, E), does not allow for simple integration if a fully relativistic cross section is used. Many authors therefore use a simpler cross section, such as the Kramer's or non-relativistic Bethe–Heitler cross section, so that the inner integral becomes analytic, although precision will be lost in calculating bremsstrahlung emissions.

However, reversing the bremsstrahlung integral has the advantage that the form of the cross section does not matter. When reversed, the new inner integral will be analytic for almost all forms of electron distributions that are assumed in thick-target models. To this end, it is rather advantageous to reverse the integral into the form given in Equation (8) (see also Kontar et al. 2011). We now give a full explanation of the evaluation of this integral, for the case of our assumed electron distribution (other forms will be similar in their evaluation).

First, we note that the integral must be split into two parts depending on what the photon energy is relative to the cut-off energy E_c (note that if the electron distribution has a high-energy break, it must be split into three parts). We begin with the case where ⩾ E_c. Disregarding preceding constants,

$\begin{eqnarray} I &\propto & \int _{E=\epsilon }^{\infty } E\ Q(\epsilon,E)\ dE \int _{E_{0} = E}^{\infty } E_{0}^{- \delta }\ dE_{0} \nonumber\\ &\propto & \int _{E=\epsilon }^{\infty } E\ Q(\epsilon,E)\ dE\ [ E^{1-\delta } ] \nonumber\\ &\propto & \int _{E=\epsilon }^{\infty } E^{2-\delta }\ Q(\epsilon,E)\ dE. \end{eqnarray} \tag{ B1 }$

One may then make a choice of cross section, depending on what physical aspects are under consideration (see Koch & Motz 1959 for an extensive list of possible choices). For all but the simplest choices, this integral requires numerical evaluation, for which Gauss–Laguerre quadrature is highly efficient.

Next, we consider the case that < E_c. Then, the integral must be split into two parts, such that ⩽ E ⩽ E_c and E_c ⩽ E < ∞, and then summed together (since electrons of any energy greater than the photon energy contribute to the integral). In the lower range,

$\begin{eqnarray} I_{\rm {low}} &\propto & \int _{E=\epsilon }^{E_{c}} E\ Q(\epsilon,E)\ dE\ \nonumber\\ && \times \left [\int _{E_{0}=E}^{E_{c}} \frac{E_{0} ^{2}}{E_{c}^{2}}\ dE_{0}\ + \int _{E_{0}=E_{c}}^{\infty } \frac{E_{0}^{- \delta }}{E_{c}^{- \delta }}\ dE_{0} \right ] \nonumber\\ &\propto & \int _{E=\epsilon }^{E_{c}} E\ Q(\epsilon,E)\ dE\ \left [\frac{E_{c}}{3} - \frac{E^{3}}{3 E_{c}^{2}} + \frac{E_{c}}{\delta - 1} \right ] \nonumber\\ &\propto & \int _{E=\epsilon }^{E_{c}} E\ \left (E_{c}\ \left [\frac{\delta +2}{3(\delta - 1)}\right ] - \frac{E^{3}}{3 E_{c}^{2}} \right)\ Q(\epsilon,E)\ dE, \quad \quad \end{eqnarray} \tag{ B2 }$

which can be evaluated efficiently with Gauss–Legendre quadrature, for any choice of cross section. In the upper range,

$\begin{eqnarray} I_{\rm {high}} &\propto & \int _{E=E_{c}}^{\infty } E\ Q(\epsilon,E)\ dE\ \int _{E_{0}=E}^{\infty } \frac{E_{0}^{- \delta }}{E_{c}^{- \delta }}\ dE_{0} \nonumber\\ &\propto & \int _{E=E_{c}}^{\infty } E^{2 - \delta }\ Q(\epsilon,E)\ dE, \end{eqnarray} \tag{ B3 }$

which is the same as the original case, except that the lower limit here is E_c. As noted before, the latter two integrals must be summed in the case < E_c.