Self-consistent Modelling of Gamma-Ray Spectra from Solar Flares with the Monte Carlo Simulation Package FLUKA

TL;DR

This paper employs the FLUKA Monte Carlo code to generate comprehensive gamma-ray spectra from solar flares, integrating various physical processes to improve modeling accuracy and fit observational data.

Contribution

It introduces a self-consistent modeling approach using FLUKA for gamma-ray spectra from solar flares, covering a broad energy range and incorporating multiple spectral components.

Findings

Good agreement with previous nuclear de-excitation models

Successful fit of the 2010 June 12 flare spectrum using FLUKA templates

First use of a single code for self-consistent broad-spectrum modeling

Abstract

We use the Monte Carlo particle physics code FLUKA (Fluktuierende Kaskade) to calculate $\gamma$-ray spectra expected from solar flare energetic ion distributions. The FLUKA code includes robust physics-based models for electromagnetic, hadronic and nuclear interactions, sufficiently detailed for it to be a useful tool for calculating nuclear de-excitation, positron annihilation and neutron capture line fluxes and shapes, as well as $\approx \; {\rm GeV}$ continuum radiation from pion decay products. We show nuclear de-excitation $\gamma$-ray line model spectra from a range of assumed primary accelerated ion distributions and find them to be in good agreement with those found using the code of Murphy et al. (2009). We also show full $\gamma$-ray model spectra which exhibit all the typical structures of $\gamma$-ray spectra observed in solar flares. From these model spectra we build templates which are incorporated into the software package Objective Spectral Executive (OSPEX) and used to fit the combined Fermi Gamma-ray Burst Monitor (GBM)/Large Area Telescope (LAT) spectrum of the 2010 June 12 solar flare, providing a statistically acceptable result. To the best of our knowledge, the fit carried out with the FLUKA templates for the full $\gamma$-ray spectrum can be regarded as the first attempt to use a single code to implement a self-consistent treatment of the several spectral components in the photon energy range from $\approx 100$s ${\rm keV}$ to $\approx 100$s ${\rm MeV}$.

Tables

Table 1. Table 1. : Compositions for the ambient atmosphere and the accelerated ions.

Element	Ambient atmosphere	Acc. ions (photospheric)	Acc. ions (impulsive flare)
H	$1.0$	$1.0$	$1.0$
He	$8.50\times {10}^{-2}$	$0.1;0.5$	$0.1;0.5$
C	$2.69\times {10}^{-4}$	$2.69\times {10}^{-4}$	$4.65\times {10}^{-3}$
N	$6.76\times {10}^{-5}$	$6.76\times {10}^{-5}$	$1.24\times {10}^{-3}$
O	$4.90\times {10}^{-4}$	$4.90\times {10}^{-4}$	$1.00\times {10}^{-2}$
Ne	$8.51\times {10}^{-5}$	$8.51\times {10}^{-5}$	$4.55\times {10}^{-3}$
Mg	$3.98\times {10}^{-5}$	$3.98\times {10}^{-5}$	$5.89\times {10}^{-3}$
Al	$2.82\times {10}^{-6}$	$2.82\times {10}^{-6}$	$1.57\times {10}^{-4}$
Si	$3.24\times {10}^{-5}$	$3.24\times {10}^{-5}$	$4.55\times {10}^{-3}$
S	$1.32\times {10}^{-5}$	$1.32\times {10}^{-5}$	$9.56\times {10}^{-4}$
Ca	$2.19\times {10}^{-6}$	$2.19\times {10}^{-6}$	$1.06\times {10}^{-4}$
Fe	$3.16\times {10}^{-5}$	$3.16\times {10}^{-5}$	$1.34\times {10}^{-2}$

Table 2. Table 2. : Parameters of the best-fit to the combined GBM/LAT spectrum of the 2010 June 12 solar flare obtained with the FLUKA templates for full γ 𝛾 \gamma -ray spectrum components.

Best-fit parameter	Value
normalization of the function 1pow at $300\mathrm{keV}$	$0.0169\pm 0.0003{\mathrm{cm}}^{-2}{\mathrm{keV}}^{-1}{s}^{-1}$
spectral index ${\delta }_e^{1pow}$ of the function 1pow	$-3.01\pm 0.03$
normalization of the function 1pow-exp at $300\mathrm{keV}$	$0.01{\mathrm{cm}}^{-2}{\mathrm{keV}}^{-1}{s}^{-1}$
spectral index ${\delta }_e^{\exp }$ of the function 1pow-exp	$-0.600$
normalization of the exp. of function 1pow-exp	$0.05{\mathrm{cm}}^{-2}{\mathrm{keV}}^{-1}{s}^{-1}$
pseudo-temp. $T_e^{\exp }$ of the function 1pow-exp	$3800\mathrm{keV}$
normalization $a{\left[0\right]}_{\mathrm{dir}}$ of the direct comp. template	$0.397\pm 0.018{\mathrm{cm}}^{-2}{s}^{-1}$
normalization $a{\left[0\right]}_{\mathrm{inv}}$ of the inverse comp. template	$0.441\pm 0.018{\mathrm{cm}}^{-2}{s}^{-1}$
Calculated parameter	Value
$a{\left[0\right]}_{\mathrm{dir}}/a{\left[0\right]}_{\mathrm{inv}}$	$0.90$
$w_{\mathrm{fac}}$	$0.12$
$N_p\left(>30\mathrm{MeV}\right)$	$9.22\times {10}^{31}$

Equations

$$\frac{\mathrm{d}{n}_i(E)}{\mathrm{d}E}=w_{i}N{E}^{-\delta }H(E_{\max }-E)H(E-E_{\min }),\text{ilabel}power-law:ions$$

$$N{\int }_{E_{\min }}^{E_{\max }}{E}^{-\delta }\mathrm{d}E=1.$$

$$\frac{\mathrm{d}{\varphi }_i(E)}{\mathrm{d}E}=N_{\mathrm{low}}\frac{\mathrm{d}{\varphi }_i^{\mathrm{low}}(E)}{\mathrm{d}E}+N_{\mathrm{high}}\frac{\mathrm{d}{\varphi }_i^{\mathrm{high}}(E)}{\mathrm{d}E}$$

$$N_{\mathrm{low}}=N{\int }_{E_{\min }^{\mathrm{low}}}^{E_{\max }^{\mathrm{low}}}{E}^{-\delta }\mathrm{d}E,$$

$$N_{\mathrm{high}}=N{\int }_{E_{\min }^{\mathrm{high}}}^{E_{\max }^{\mathrm{high}}}{E}^{-\delta }\mathrm{d}E.$$

$$\frac{\mathrm{d}\varphi (E)}{\mathrm{d}E}=\sum _i{w}_i\frac{\mathrm{d}{\varphi }_i(E)}{\mathrm{d}E}.$$

$$\frac{\mathrm{d}{\varphi }_{\mathrm{dir}}(E)}{\mathrm{d}E}=w_p\frac{\mathrm{d}{\varphi }_p(E)}{\mathrm{d}E}+w_{\alpha }\frac{\mathrm{d}{\varphi }_{\alpha }(E)}{\mathrm{d}E},\text{ilabel}dir-spec$$

$$\frac{\mathrm{d}{\varphi }_{\mathrm{inv}}(E)}{\mathrm{d}E}=\sum _{i\ne p,\alpha }{w}_i\frac{\mathrm{d}{\varphi }_i(E)}{\mathrm{d}E}.\text{ilabel}inv-spec$$

$$\frac{\mathrm{d}{\varphi }_{\mathrm{dir}}(E)}{\mathrm{d}E}$$

$$\frac{\mathrm{d}{\varphi }_{\mathrm{inv}}(E)}{\mathrm{d}E}$$

$$N_{30}={\left[{\int }_{30\mathrm{MeV}}^{1\mathrm{GeV}}{E}^{-\delta }\mathrm{d}E\right]}^{-1}.$$

$$r_{i}\,=\,\frac{{\rm d}\phi(E)}{{\rm d}E}\bigg{|}_{i,{\rm RMK}}-\frac{{\rm d}\phi(E)}{{\rm d}E}\bigg{|}_{i,{\rm FLUKA}}\;,$$

$$\sigma =\sqrt{\frac{1}{N-1}\sum _i{r}_i^2},$$

$$\frac{\mathrm{d}{\varphi }_{\mathrm{dir}}(E)}{\mathrm{d}E}\to \frac{\mathrm{d}{\tilde{\varphi }}_{\mathrm{dir}}(E)}{\mathrm{d}E}$$

$$\frac{\mathrm{d}{\varphi }_{\mathrm{inv}}(E)}{\mathrm{d}E}\to \frac{\mathrm{d}{\tilde{\varphi }}_{\mathrm{inv}}(E)}{\mathrm{d}E}$$

$$I_{\mathrm{dir}}={\int }_{E_{\min }^{\mathrm{phot}}}^{E_{\max }^{\mathrm{phot}}}\frac{\mathrm{d}{\varphi }_{\mathrm{dir}}(E)}{\mathrm{d}E}\mathrm{d}E,$$

$$I_{\mathrm{inv}}={\int }_{E_{\min }^{\mathrm{phot}}}^{E_{\max }^{\mathrm{phot}}}\frac{\mathrm{d}{\varphi }_{\mathrm{inv}}(E)}{\mathrm{d}E}\mathrm{d}E.$$

$$\frac{\mathrm{d}{\tilde{\varphi }}_{\mathrm{tot}}(E)}{\mathrm{d}E}=a[0{]}_{\mathrm{dir}}\frac{\mathrm{d}{\tilde{\varphi }}_{\mathrm{dir}}(E)}{\mathrm{d}E}+a[0{]}_{\mathrm{inv}}\frac{\mathrm{d}{\tilde{\varphi }}_{\mathrm{inv}}(E)}{\mathrm{d}E},$$

$$\frac{\mathrm{d}{\tilde{\varphi }}_{\mathrm{tot}}(E)}{\mathrm{d}E}=\frac{a[0{]}_{\mathrm{dir}}}{I_{\mathrm{dir}}}\left[\frac{\mathrm{d}{\varphi }_{\mathrm{dir}}(E)}{\mathrm{d}E}+\frac{\mathrm{d}{\varphi }_{\mathrm{inv}}^{\prime }(E)}{\mathrm{d}E}\right],$$

$$\frac{\mathrm{d}{\varphi }_{\mathrm{inv}}^{\prime }(E)}{\mathrm{d}E}=\left(\frac{a[0{]}_{\mathrm{inv}}{I}_{\mathrm{dir}}}{a[0{]}_{\mathrm{dir}}{I}_{\mathrm{inv}}}\right)\frac{\mathrm{d}{\varphi }_{\mathrm{inv}}(E)}{\mathrm{d}E}.$$

$$w_{\mathrm{fac}}=\left(\frac{a[0{]}_{\mathrm{inv}}{I}_{\mathrm{dir}}}{a[0{]}_{\mathrm{dir}}{I}_{\mathrm{inv}}}\right).$$

$$I_{\mathrm{tot}}^{\prime }$$

$$N_p(>30\mathrm{MeV})=a[0{]}_{\mathrm{dir}}\frac{4\pi R^2}{I_{\mathrm{tot}}^{\prime }}\Delta t=a[0{]}_{\mathrm{dir}}\frac{4\pi R^2}{\left(1+\frac{a[0{]}_{\mathrm{inv}}}{a[0{]}_{\mathrm{dir}}}\right)I_{\mathrm{dir}}}\Delta t,$$

Full text

\setlastpage\inarticletrue{opening}

1 Introduction

\ilabel

S-Introduction

Particle acceleration plays a key role in solar flares. Accelerated particles may constitute a substantial fraction of all the energy released in a flare Emslie et al. (2004, 2005). Their radiative signatures are characteristic of the flare impulsive phase and clearly key to the energy release process Kane (1974). Moreover, empirical knowledge of flares serves as “ground truth” for other astrophysical energy release events in more remote settings. Energetic particles are detected in situ in interplanetary space but we are reliant on electromagnetic radiation to learn about particle acceleration at the flare site and it is important that we develop reliable methods to interpret spectral, spatial and timing information.

Observations of (electron bremsstrahlung) flare hard X-rays with the Ramaty High-Energy Spectroscopic Imager (RHESSI) satellite have yielded a refined characterisation of flare electron acceleration and its relationship to other flare processes Holman et al. (2011); Kontar et al. (2011). The picture for flare ions remains less complete, partly because flares producing detectable $\gamma$-rays are much rarer but also because of the greater complexity of the nuclear radiation mechanisms. Ion acceleration in flares is studied via a variety of signatures at $\gamma$-ray wavelengths, consequences of nuclear reactions of primary accelerated ions with ambient nuclei Share and Murphy (2006); Murphy et al. (2007); Vilmer, MacKinnon, and Hurford (2011): nuclear de-excitation lines in the $0.4-7\;{\rm MeV}$ range, excited by ions of $1-100\;{\rm MeV/nucleon}$ Ramaty, Kozlovsky, and Lingenfelter (1979); Kozlovsky, Murphy, and Ramaty (2002); $2.223\;{\rm MeV}$ neutron capture and $511\;{\rm keV}$ positron annihilation lines, each resulting from ions spanning a wide energy range, from a few to $100$s of ${\rm MeV/nucleon}$ Lockwood, Debrunner, and Ryan (1997); Hua et al. (2002); Kozlovsky, Lingenfelter, and Ramaty (1987); Kozlovsky, Murphy, and Ramaty (2002); Murphy et al. (2007); continuum radiation in the $\gtrsim 0.1\;{\rm MeV}$ range, resulting from pion decay products of $\gtrsim 0.2-0.3\;{\rm GeV/nucleon}$ ions Dermer (1986b); Murphy, Dermer, and Ramaty (1987); Mandzhavidze and Ramaty (1992); Vilmer et al. (2003). Energetic neutrons detected in space or with ground-based instruments also give information on accelerated flare ions (Chupp et al., 1987; Kocharov et al., 1998; Watanabe et al., 2006, e.g.).

Interpretation of observed $\gamma$-ray spectra must proceed via best-fit from among a set of templates that span a range of assumptions about the ion energy distributions and source chemical abundances, and may also involve details of magnetic geometry, etc. Murphy et al. (2007); Vilmer, MacKinnon, and Hurford (2011). The generation and use of such templates has been developed over several decades, starting with the seminal work of Lingenfelter and Ramaty (1967). A detailed study of the production of nuclear de-excitation $\gamma$-ray lines was carried out by Ramaty, Kozlovsky, and Lingenfelter (1979), resulting in the development of a Monte Carlo code for calculating complete spectra of de-excitation $\gamma$-ray lines. The code includes both the “direct” reactions of energetic protons, 3He and $\alpha$-particles with ambient nuclei and the “inverse” reactions of energetic nuclei heavier than $\alpha$-particles with ambient H and 4He. This code was updated by Kozlovsky, Murphy, and Ramaty (2002) with new cross section data from laboratory measurements. An improved version of the code (hereafter referred to as RMK) was further developed by Murphy et al. (2009) incorporating cross sections evaluated with the nuclear reaction code TALYS Koning, Hilaire, and Duijvestijn (2005). Templates for the nuclear de-excitation $\gamma$-ray lines are available as standard fit functions in the software package Objective Spectral Executive (OSPEX) Freeland and Handy (1998), which were built from calculations with the RMK code assuming a downward isotropic angular distribution of primary accelerated ions having a power-law energy distribution with spectral index of $4.0$, a coronal abundance of accelerated ions with $\alpha/{\rm proton}$ ratio of $0.22$, a coronal abundance of ambient nuclei with ${\rm He/H}$ ratio of $0.1$, and a heliocentric angle of $60^{\circ}$. Templates for the continuum radiation from pion decay are also available in OSPEX, which were built from calculations of pion decay $\gamma$-ray spectra carried out by Murphy, Dermer, and Ramaty (1987) assuming an isotropic angular distribution of primary accelerated protons having a power-law energy distribution with several spectral indexes ($3.0,\;3.5,\;4.0,\;4.5,\;{\rm and}\;5.0$), and $\alpha/{\rm proton}$ ratio of $0.1$. OSPEX provides a specific standard fit function for the component due to the $511\;{\rm keV}$ positron annihilation line and its associated positronium continuum, but there is no specific standard fit function available for the component due to the $2.223\;{\rm MeV}$ neutron capture line. The neutron capture line component is usually included as a gaussian function. A more accurate neutron capture line component also incorporating the associated Compton-scattered continuum can be included if desired through user-supplied templates built from calculations with codes such as the neutron-production code developed by Hua et al. (2002) and the Monte Carlo N-Particle Transport Code (MCNP6) developed by Pelowitz (2013) (see e.g. Murphy and Share (2018)). OSPEX also provides a standard fit function for the component due to the $\alpha-\alpha$ lines, which was built from calculations assuming a downward isotropic distribution of primary accelerated ions having a power-law energy distribution with spectral index of $3.5$ and a heliocentric angle of $60^{\circ}$.

The Monte Carlo code FLUKA (Fluktuierende Kaskade) Ferrari et al. (2011); Battistoni et al. (2015) is a general purpose package of integrated routines for simulation of particle transport and interactions in arbitrary materials. It is “not a toolkit”: it aims to provide a unified description of all relevant processes, to a consistent degree of approximation, relevant to a particular situation. It has seen a wide variety of applications in areas such as high-energy experimental physics, nuclear physics and engineering, detector and telescope design, and medical physics, as well in applications particularly close to our interests here: calculating consequences of cosmic-ray interaction in the terrestrial atmosphere Battistoni et al. (2008); Claret et al. (2014); interpretation of Fermi Large Area Telescope (LAT) measurements of $\gamma$-rays produced by cosmic-ray impact on the Moon Ackermann et al. (2016).

Recent contributions have begun to explore FLUKA’s use as a tool for interpretation of solar flare radiations at $\gamma$-ray and other wavelengths. Since FLUKA combines physics-based models for nuclear processes with transport of primary and secondary particles and photons, MacKinnon (2014) and MacKinnon et al. (2019) were able to use it to study the effects of primary particle directionality on emergent spectra. They concentrated on the $\gamma$-ray continuum at $\approx$ GeV energies resulting from neutral pion decay and from bremsstrahlung of charged pion decay positrons and electrons. It was also used to calculate the distributions of charged secondaries emerging from a photospheric thick target, with a view to evaluating their synchrotron radiation at ${\rm mm}$ and shorter wavelengths Tuneu et al. (2017). Here we concentrate particularly on FLUKA’s capabilities for predicting de-excitation line fluxes and present, for the first time, full, self-consistent simulations of the $\gamma$-ray spectrum, including the nuclear de-excitation lines, the $2.223\;{\rm MeV}$ neutron capture line and the pion decay continuum. Our spectra also include the $511\;{\rm keV}$ positron annihilation line from pion decay products and pair creation in the solar atmosphere by $\gamma$-ray photons but exclude the normally dominant contribution from secondary radioactive positron-emitting nuclei. FLUKA is capable of modelling positron production via radioactive nuclei but we defer a detailed study of all components of the annihilation line for future work, noting that a full treatment of the shape of the line, the relative intensities of line and three-photon continuum, etc., need a description of solar atmospheric properties, particularly degree of ionization, that we do not attempt to implement in FLUKA (see Murphy et al. (2005)). With this qualification, the resulting spectra offer the possibility of an independent verification and validation of the work of Ramaty, Kozlovsky, and Lingenfelter (1979), Kozlovsky, Murphy, and Ramaty (2002), etc. De-excitation line production in FLUKA occurs at the end of a sequence of nuclear processes that model hadron-nucleus and nucleus-nucleus interaction. Thus the fluxes and line shapes resulting from many sorts of reactions (inelastic excitation as well as spallation and from both direct and inverse reactions) are modeled in detail. FLUKA’s several physics models are detailed in the Appendix, along with several references which provide their validation by comparison to laboratory measurements.

The Monte Carlo code Geometry and Tracking (GEANT4) Agostinelli et al. (2003) is also rooted in laboratory expertise and is employed in a similarly wide range of contexts. Tang and Smith (2010) used GEANT4 in a similar spirit to our work here, concentrating on calculating the pion decay continuum radiation from postulated primary accelerated ion distributions, as well as the $511\;{\rm keV}$ positron annihilation line. They did not use it to calculate de-excitation $\gamma$-ray line spectra. Since then, the Fermi Gamma-ray Burst Monitor (GBM) and LAT instruments have obtained good $\gamma$-ray spectra which can be compared with such calculations, as we do here. It is important to emphasize that the main purpose of the paper is to present the results of a first investigation on FLUKA’s potential as an effective tool for the modelling of $\gamma$-ray spectra from solar flares. A thorough analysis to validate FLUKA for such modelling would require calculations which are beyond the scope of the paper.

In Section \irefs-model we describe in detail the geometrical model we used to simulate $\gamma$-rays from flare ions with FLUKA. In Section \irefS-simulations we show illustrative nuclear de-excitation $\gamma$-ray line spectra and full $\gamma$-ray spectra obtained from this model. In Section \irefS-fitting we show a detailed fit to Fermi data. Since FLUKA provides an integrated description of all relevant processes we can attempt to provide a single fit to a full $\gamma$-ray spectrum stretching from $<1\;{\rm MeV}$ to ${\rm GeV}$ photon energies, all emerging from the same physical situation. To this end we combine Fermi GBM and LAT data to give a single spectrum including both lines and pion decay continuum (cf. Vilmer et al. (2003)). Section \irefS-summary gives a brief summary of our approach and results and discusses their usefulness and implications.

2 Solar Flare Model

\ilabel

s-model

In our simulations of nuclear processes in solar flares using the FLUKA package, we consider a model in which beams of primary accelerated ions are injected into a target with characteristics similar to those of the ambient solar atmosphere. FLUKA’s combinatorial geometry package Ferrari et al. (2011), which is an improved version of the package developed for the neutron and $\gamma$-ray transport code MORSE Emmett (1975), allows one to set up simulations with arbitrarily complex geometries. A huge variety of geometries can be built through combinations of solid bodies and surfaces obtained by boolean operations. We adopt a simple plane-parallel geometry for the vertical structure of the ambient solar atmosphere, since in general the dimensions of the solar flare emission region are much smaller than the solar radius111The assumption of plane-parallel geometry is adequate over most of the disk. For flares very near the limb it over-estimates the column density met by escaping photons and must be replaced by a spherical geometry, straightforwardly implemented in FLUKA if necessary.. Using the combinatorial geometry tools provided by FLUKA, we build a cubic box centered at the origin of a cartesian coordinate system (${\rm O}x,{\rm O}y,{\rm O}z$) with edges of length $L=2\times 10^{9}\;{\rm cm}$ and faces perpendicular to the coordinate axes. The $z$-coordinate corresponds to the vertical depth in the ambient solar atmosphere. A $xy$-plane at $z=0$ divides the cubic box into two half-spaces. The half-space at $z<0$ represents the coronal region and, for simplicity, is just filled with vacuum (such that the particles are transported but no longer interact). The half-space at $z>0$ represents the chromospheric/photospheric region and is filled with a dense, neutral material for which we assume a typical solar atmosphere composition with the abundances of 4He, C, N, O, Ne, Mg, Al, Si, S, Ca and Fe nuclei relative to H given by Asplund et al. (2009), as indicated in Table \ireftab:compAtSolarModIII. Our assumption of a neutral medium is appropriate to the photosphere and consistent with previous treatments of $\gamma$-ray line production (e.g. Murphy et al. (2007)). The beams of primary accelerated ions are injected into the chromospheric/photospheric region from a point ($0,0,z_{0}$) located in the coronal region at a distance $z_{0}$ above and very close to the $xy$-plane at $z=0$, avoiding any artefacts or errors that might result if the source were placed exactly on the boundary between two regions.

The half-space at $z>0$ is further divided into $52$ layers with a vertical density profile given by the semi-empirical VAL-C model of the chromosphere Vernazza, Avrett, and Loeser (1981) plus an additional layer, corresponding to the photosphere, with a density of $3.19\times 10^{-7}\;{\rm g/cm^{3}}$. The densities of the chromospheric layers given in the VAL-C model range from $2.35\times 10^{-15}\;{\rm g/cm^{3}}$ at the outermost layer to $3.19\times 10^{-7}\;{\rm g/cm^{3}}$ at the innermost one. However, materials with densities $<10^{-10}\;{\rm g/cm^{3}}$ are treated as vacuum by FLUKA. In order to circumvent this feature, we assign a density of $2\times 10^{-10}\;{\rm g/cm^{3}}$ to all layers with densities $<10^{-10}\;{\rm g/cm^{3}}$ and rescale their thickness while keeping the values for the column-densities as given in the VAL-C model. The total column-density obtained by summing up the contributions of all $53$ layers is $292\;{\rm g/cm^{2}}$, corresponding to the stopping-depth (or range) for a proton with energy of $\approx 1.6\;{\rm GeV}$ in a hydrogen target evaluated in the continuous slowing down approximation (CSDA) Berger et al. (1999). In this way, the half-space at $z>0$ fulfils the thick-target condition for primary accelerated ions with kinetic energies up to $1.6\;{\rm GeV/nucleon}$.

It is important to note that the modification we have implemented in the density structure (with respect to the original VAL-C model) does not introduce significant changes in the results of our present simulations, since the thick target photon yields produced by all relevant processes, with the exception of particle decay processes, depend only on the total column-density. The decay processes depend on time, and thus on the distance travelled by the decaying particles within the target. However, the decay processes which are most relevant to our simulations are the decay of neutral pions ($\pi^{0}$), charged pions ($\pi^{\pm}$) and muons ($\mu^{\pm}$), whose lifetimes are respectively $\tau_{\pi^{0}}=8.5\times 10^{-17}\;{\rm s}$, $\tau_{\pi^{\pm}}=2.6\times 10^{-8}\;{\rm s}$ and $\tau_{\mu^{\pm}}=2.2\times 10^{-6}\;{\rm s}$. A muon or pion with speed very close to $c$ thus travels $<1$ km in a rest-frame lifetime. Primary protons with energies up to $10\;{\rm GeV}$ can produce muons with Lorentz gamma factor up to $\approx 70$ (using results of Dermer (1986a)) but this still means they travel $\lesssim 50-100\;{\rm km}$ before decaying. Compared to the system scale they decay more or less where they are produced and their daughter particles will produce their full thick target yield of bremsstrahlung photons. Radioactive positron-emitter nuclei can have lifetimes of several minutes and longer but, as noted earlier, are not included in our simulations here.

In the simulations carried out in this work we consider primary accelerated ions both with “photospheric” composition and with “impulsive-flare” composition as defined in Mandzhavidze and Ramaty (1993), Ramaty, Mandzhavidze, and Kozlovsky (1996) and Murphy et al. (1997). For the photospheric composition we assume that the abundances of the primary accelerated ions heavier than $\alpha$-particles relative to protons are the same as for the ambient nuclei relative to H, i.e. the abundances of Asplund et al. (2009). For the impulsive-flare composition we assume the enhanced relative abundances for the primary accelerated heavy ions given in Murphy et al. (2007), as indicated in Table \ireftab:compAtSolarModIII, and adopt values of $0.1$ and $0.5$ for the $\alpha$/proton ratio. We do not include ${}^{3}{\rm He}$ in the compositions of the ambient nuclei and the primary accelerated heavy ions, although it may be important in various ways for a full interpretation of the $\gamma$-ray spectrum Murphy, Kozlovsky, and Share (2016). In any case FLUKA’s treatment of reactions involving 3He is incomplete at present. One should also note that we assume the same primary accelerated heavy ion abundances for both values of the $\alpha$/proton ratio, although according to Murphy et al. (2007) in the case where the $\alpha$/proton ratio is $0.1$ the abundances should be reduced by a factor of $5$. Nevertheless, one should note that our fits to observed $\gamma$-ray spectra from solar flares have been carried out by using two separate templates for the components due to the direct and the inverse reactions, such that the primary accelerated heavy ion abundances are actually multiplied by a common factor which depends on the yields of the direct and the inverse components and their normalizations determined by the fitting procedure (see Equation \ireffactor).

We run separate simulations for the reactions involving each primary accelerated ion-species $i$ and all ambient nuclei. The primary accelerated ions are assumed to have power-law energy distributions given by:

[TABLE]

where $w_{i}$ is the relative abundance for the ion-species $i$, $E$ is the primary accelerated ion kinetic energy per nucleon in the range from $E_{\rm min}=1\;{\rm MeV/nucleon}$ to $E_{\rm max}=1\;{\rm GeV/nucleon}$, $\delta$ is the power-law spectral index, $H$ is the Heaviside step function, and $N$ is a normalization constant defined such that

[TABLE]

FLUKA allows the user to customise various aspects of the simulation by changing the code in certain subroutines. We modified the subroutine source.f so that it either draws primary ion energies from the power-law distribution Equation \irefpower-law:ions, or equivalently draws them randomly from a uniform distribution but assigns weights given by the power-law distribution MacKinnon et al. (2019). The latter strategy accounts for the contribution of the higher-energy particles to the photon spectrum without suffering large statistical fluctuations resulting from their comparative rarity. We generally found that $10^{7}-10^{8}$ primaries were necessary for a statistically reliable spectrum. Our modified version of source.f provides a number of options for primary ion angular distribution MacKinnon et al. (2019) but in this work we always assume that ion directions are drawn from a uniform distribution in the downward hemisphere (“downward isotropic”).

The algorithm implemented by FLUKA follows the evolution of primary and secondary particles individually as they propagate through the target, tracking their interactions with the ambient medium until they leave the target region, come to rest, or reach a low-energy threshold for transport. Several kinds of estimation tools are available, referred to as “detectors”, which allow one to calculate different quantities of interest. In the simulations carried out with the model described above we use a specific detector, named USRBDX, to determine the energy spectrum of photons escaping from the chromospheric/photospheric region to the coronal region. This kind of detector calculates the double differential distribution of the flux of photons crossing the $xy$-plane at $z=0$ per energy and solid angle intervals, ${\rm d}^{2}\phi/{\rm d}E{\rm d}\Omega$, in units of photons ${\rm GeV}^{-1}\;{\rm cm}^{-2}\;{\rm sr}^{-1}$ per primary accelerated ion. Integrating the flux of photons in the upward direction over the solid angle and the whole of the surface area defined by the $xy$-plane at $z=0$, we obtain the corresponding energy spectrum of escaping photons, ${\rm d}\phi/{\rm d}E$, in units of photons ${\rm GeV}^{-1}$ per primary accelerated ion – i.e. the photon spectrum emitted into the backward hemisphere (coronal region). The user determines the energy resolution of the detector by specifying the number of bins and the maximum and minimum energies of photons to be recorded. We chose these by experiment to obtain a reasonable compromise between spectral resolution and statistical reliability. We also set additional detectors with narrower, $3\;{\rm keV}$ energy bins to give better resolution around the strong $511\;{\rm keV}$ and $2.223\;{\rm MeV}$ lines, combining these with the full spectrum in the results shown below. In Section \irefS-Lines we compare spectra of nuclear de-excitation $\gamma$-ray lines obtained in simulations with FLUKA and in calculations carried out with a copy of the RMK code kindly supplied by R. J. Murphy. Because the RMK calculated spectra are obtained by integrating the flux of photons in all directions, spectra calculated with FLUKA are expected to be nearly a factor of $2$ less intense than those calculated with the RMK code. Thus, for a proper comparison of the spectra calculated with the two codes, in the figures shown we will decrease the RMK spectra by a factor two.

As pointed out before, the contributions to the solar flare $\gamma$-ray spectrum from the emission of nuclear de-excitation, neutron capture and positron annihilation lines are produced mostly by primary accelerated ions with energies in the range from $\approx 1$ to $100\;{\rm MeV/nucleon}$. The continuum emission contributions from pion decay processes, on the other hand, are produced by primary ions with energies above $\approx 200-300\;{\rm MeV/nucleon}$. Thus, in order to improve the statistics for the photon yields generated by all emission processes we run separate simulations for two energy ranges, namely $E_{\rm min}^{\rm low}=1\;{\rm MeV/nucleon}$ to $E_{\rm max}^{\rm low}=200\;{\rm MeV/nucleon}$ and $E_{\rm min}^{\rm high}=200\;{\rm MeV/nucleon}$ to $E_{\rm max}^{\rm high}=1\;{\rm GeV/nucleon}$. The energy spectrum of photons, ${\rm d}\phi_{i}(E)/{\rm d}E$, for a given ion-species $i$ is obtained by summing up the contributions from the two energy ranges with the appropriate weighting factors,

[TABLE]

where the weights $N_{\rm low}$ and $N_{\rm high}$ are given by:

[TABLE]

Thus, by summing up the contributions from each ion-species $i$ weighted by their corresponding relative abundances, $w_{i}$, we obtain the total energy spectrum of photons in units of photons ${\rm GeV}^{-1}$ per primary accelerated proton in the energy range from $E_{\rm min}$ to $E_{\rm max}$,

[TABLE]

Similarly, we can obtain the energy spectra of photons produced by the direct reactions (involving interactions of primary accelerated protons and $\alpha$-particles with all ambient nuclei) and by the inverse reactions (involving interactions of primary accelerated ions heavier than $\alpha$-particles with ambient nuclei of H and 4He), given by

[TABLE]

[TABLE]

Since the RMK spectra are normalized to one primary accelerated proton in the energy range from $30\;{\rm MeV}$ to $1\;{\rm GeV}$, we renormalize the FLUKA spectra of nuclear de-excitation $\gamma$-ray lines in order to allow for a proper comparison, i.e.,

[TABLE]

where $N_{30}$ is the normalization constant for the power-law energy distributions considered in the calculations with the RMK code, given by:

[TABLE]

The same renormalization was carried out for the full $\gamma$-ray spectra obtained in simulations with FLUKA shown in Section \irefS-Full.

3 Results of Simulations

\ilabel

S-simulations

In this section we present the results obtained in simulations carried out with FLUKA using the geometrical model for the ambient solar atmosphere described above. In Section \irefS-Lines we compare the nuclear de-excitation $\gamma$-ray line spectra obtained in simulations with FLUKA and in calculations with the RMK code. In Section \irefS-Full we show the results obtained in simulations with FLUKA for the full $\gamma$-ray spectra, i.e. with other relevant photon emission processes apart from the positron annihilation line produced via radioactive positron-emitter daughter nuclei taken into account. In all cases, both for RMK and FLUKA spectra, we consider downward isotropic beams of primary accelerated ions with power-law energy distribution of spectral index $\delta=4$ in the range from $1\;{\rm MeV/nucleon}$ to $1\;{\rm GeV/nucleon}$.

3.1 Nuclear De-excitation $\gamma$-Ray Line Spectrum

\ilabel

S-Lines

The RMK code Ramaty, Kozlovsky, and Lingenfelter (1979); Kozlovsky, Murphy, and Ramaty (2002); Murphy et al. (2009) is one of the main tools used for the analysis of solar flare $\gamma$-ray data. The code is based on a Monte Carlo algorithm which allows one to calculate the “explicit” narrow and broad nuclear de-excitation $\gamma$-ray lines and the so called “unresolved” component. By “explicit” line is meant one of the more than $200$ lines that are explicitly treated by the code. These are the strongest lines resulting from transitions of the excited nuclei produced by interactions of the most abundant nuclei in the solar atmosphere. The unresolved component is composed of the thousands of remaining lines produced by the interactions which are relatively weak and closely spaced. They are treated by the code as a “quasi-continuum”. The nuclear de-excitation $\gamma$-ray line spectra generated with the RMK code are evaluated by integrating over solid angle the flux of photons in all directions and take into account the dependence of the Doppler-shifts and the relativistic beaming of moving emitter nuclei on the heliocentric angle $\theta_{\rm obs}$. The updated version of the code Murphy et al. (2009) provides the option to include or not the effects of Compton scattering. The RMK spectra shown in this work are evaluated with Compton scattering turned off.

The production of nuclear de-excitation $\gamma$-ray lines in FLUKA occurs at the final stage of a sequence of nuclear processes implemented through a model called PEANUT (Pre-Equilibrium Approach to Nuclear Thermalization), developed by Ferrari and Sala (1998) (see the Appendix for details). In order to concentrate on the nuclear de-excitation $\gamma$-ray line spectra we use the functionalities of FLUKA which allow one to suppress the processes of bremsstrahlung, production and annihilation of electron-positron pairs, neutron capture, Compton scattering and pion decay. The FLUKA calculations include production of de-excitation-lines from excited nuclei produced by inelastic excitation reactions and from spallation reactions. We do not include fusion reactions in the FLUKA calculations. Fusion reactions are included in the RMK calculations, but apart from the $\alpha$-$\alpha$ fusion lines near $0.45\;{\rm MeV}$ and $3.6\;{\rm MeV}$, fusion reactions do not play a major role in the spectra of nuclear de-excitation $\gamma$-ray lines. It is possible to enable fusion reactions in FLUKA and model the $\alpha$-$\alpha$ feature if desired.

The nuclear de-excitation $\gamma$-ray line spectra generated in our FLUKA simulations (as well as the full $\gamma$-ray spectra) are evaluated by integrating over solid angle only the flux of photons in the upward direction (i.e., the flux of photons escaping from the chromospheric/photospheric region to the coronal region). In this way, for downward isotropic beams of primary accelerated ions we obtain photon spectra which correspond to an average of spectra for flares at heliocentric angles between $\theta_{\rm obs}=0^{\circ}$ (disk-centered) and $\theta_{\rm obs}=90^{\circ}$ (limb) and so exhibit Doppler-shifts similar to that of a spectrum for a flare at $\theta_{\rm obs}=60^{\circ}$, considering that the magnitude of the Doppler-shifts roughly changes linearly with ${\rm cos}(\theta_{\rm obs})$ from [math] for a limb flare to a minimum negative value for a disk-centered flare Harris et al. (2007). Thus, in order to perform a more consistent comparison between the nuclear de-excitation $\gamma$-ray line spectra obtained in simulations with FLUKA and in calculations with the RMK code, we generate the RMK spectra for $\theta_{\rm obs}=60^{\circ}$. Photon spectra which approximately correspond to spectra for specific heliocentric angles $\theta_{\rm obs}$ could be obtained in our FLUKA simulations by setting up a high-resolution solid angle binning and integrating over very small solid angle intervals, cf. the angle-dependent spectra at higher energies obtained using GEANT4 by Tang and Smith (2010). We would have to substantially increase run times to obtain such angle-dependent spectra with good statistics and small enough energy bins to see the effects on line shapes, however, and we leave such an effort to future work.

In Figures \irefF-RMK-Lines-D4 and \irefF-FLUKA-Lines-D4 we show nuclear de-excitation $\gamma$-ray line spectra obtained, respectively, in calculations with the RMK code and in simulations with FLUKA for downward isotropic beams of primary accelerated ions with photospheric and impulsive-flare compositions (assuming $\alpha/{\rm proton}=0.1$ and $\delta=4$ in both cases).

In Figures \irefF-RMK-FLUKA-Lines-D4-direct, \irefF-RMK-FLUKA-Lines-D4-inverse and \irefF-RMK-FLUKA-Lines-D4-total we compare FLUKA and RMK spectra, showing respectively the contributions from direct and inverse reactions and the total spectrum. For a better comparison, we have rebinned the higher-resolution RMK spectra into the same binning as that of the FLUKA spectra. As one can observe, the spectra obtained in the simulations with FLUKA can satisfactorily reproduce those obtained in calculations with the RMK code. The $0.4-0.5\;{\rm MeV}$ spectral feature from $\alpha-\alpha$ fusion reactions (the superposition of $\gamma$-ray lines produced by de-excitation of 7Be and 7Li nuclei) does not appear in the FLUKA spectra shown because, as noted above, fusion reactions are not included in the calculations.

Figures \irefF-RMK-FLUKA-Lines-D4-direct to \irefF-RMK-FLUKA-Lines-D4-total show very good agreement between the FLUKA and RMK results, particularly considering that the two codes were developed independently from different starting points. Most of the major line features are present in the FLUKA results, albeit sometimes with poorer statistical quality. In the case of the strongest lines, those produced by de-excitation of 12C at $4.438\;{\rm MeV}$ and 16O at $6.129$, $6.916$ and $7.115\;{\rm MeV}$, FLUKA also gives similar results for line widths and shapes. The normalised residuals displayed in the bottom panels of the figures show the differences between the FLUKA and RMK photon spectra. In the $i$th energy bin we calculate the difference between the photon fluxes produced by the two codes,

[TABLE]

and divide by the standard deviation $\sigma$,

[TABLE]

where $N$ is the total number of photon energy bins. The residuals focus attention on a few noteworthy discrepancies. First, as already noted, we expect disagreement around $\approx 0.4-0.5$ MeV because the $\alpha$-$\alpha$ lines are not included. This omission could be remedied for more detailed comparison with data. Although the shapes of the strongest lines are similar there are some slight differences, partly attributable to FLUKA’s treatment of primary ion transport, modelling the detailed development of the initially downward isotropic ion distribution in the atmosphere, but mostly reflecting noise in the FLUKA results. Discrepancies close to the two strongest lines between 1 and 2 MeV, at $1.37$ MeV from 24Mg and $1.63$ MeV from 20Ne, are also attributable to statistical fluctuations in the FLUKA spectrum, especially combined with the width of the photon energy bins used in our FLUKA results. Lastly, although the low absolute fluxes mean it is not highlighted by the residuals, the nuclear continuum above $\approx 7$ MeV is systematically substantially lower in the FLUKA results. This reflects the different treatments of the continuum in PEANUT (used in FLUKA) and in TALYS (used in RMK). Further detailed study will be needed to fully resolve this discrepancy. Lacking narrow line features, the contribution from the inverse reactions is less sensitive to the low resolution of the energy binning and all the main features from the inverse reactions shown in the RMK spectra are well reproduced (Figure \irefF-RMK-FLUKA-Lines-D4-inverse).

To make this comparison we ignore attenuation of photons originally moving upward away from the Sun (i.e., Compton scattering and pair production above $1.022\;{\rm MeV}$ are turned off). This is reasonable because the (hydrogen) photon absorption length is greater than $\approx 5\;{\rm g/cm^{2}}$ for photon energies $>0.5\;{\rm MeV}$ Berger et al. (2010), whereas protons with energies $<100\;{\rm MeV}$ all stop in the atmosphere higher than this, at column-densities $\lesssim 4\;{\rm g/cm^{2}}$ Berger et al. (1999). However, Compton down-scattering and escape of $>1\;{\rm MeV}$ photons originally moving downward into the Sun does contribute to the spectrum of nuclear de-excitation $\gamma$-ray lines at energies $<1\;{\rm MeV}$. In general, this Compton-scattered component has no significant impact on fits to typical $\gamma$-ray spectra from solar flares observed with instruments such as RHESSI which have a strong off-diagonal detector response Murphy and Share (2018). As expected, in the case of beams of primary accelerated ions with impulsive-flare composition the contribution from the inverse reactions significantly increases, since the relative abundance of the heavy ions is enhanced.

3.2 Full $\gamma$-Ray Spectrum

\ilabel

S-Full

In our FLUKA simulations of the full $\gamma$-ray spectrum we turn on the processes of bremsstrahlung, production and annihilation of electron-positron pairs, neutron capture and pion decay. As previously noted, the contribution from the decay of radioactive positron-emitter daughter nuclei to the $511\;{\rm keV}$ positron annihilation line is not included but can be accounted for, if desired, by enabling the radioactive decay of nuclei and the transport of decay radiation in the FLUKA simulations. We also turn on Compton scattering, so that the attenuation of the $2.223\;{\rm MeV}$ neutron capture line and its associated Compton-scattered continuum are taken into account.

In Figure \irefF-FLUKA-Full-D4 we show the full $\gamma$-ray spectra obtained for downward isotropic beams of primary accelerated ions with both photospheric and impulsive-flare compositions (assuming $\alpha/{\rm proton}=0.1$ and $\delta=4$ in both cases). As one can observe, the FLUKA spectra exhibit all the typical structures of $\gamma$-ray spectra observed in solar flares: positron annihilation line (though excluding the substantial contribution from radioactive nuclei), neutron capture line, nuclear de-excitation $\gamma$-ray lines and continuum emission components from neutral-pion decay and bremsstrahlung of secondary electrons and positrons from charged-pion decay. One can also see that the component due to the direct reactions clearly exhibits the Compton-scattered continuum from the neutron-capture line at energies below $2.223\;{\rm MeV}$. This component is intense enough that it conceals the similar Compton component from Compton scattering of the de-excitation lines described by Murphy and Share (2018). The continuum emission component from the bremsstrahlung of primary accelerated electrons is not included here. As well as the abundances of the heavier species, the intensity of the pion decay continuum from inverse reactions is influenced by the larger cross-sections and lower energy thresholds for pion production in reactions with primary accelerated heavy ions. The relative magnitudes of the inverse and direct reaction fluxes are in general agreement with those found by Kafexhiu et al. (2018). When spectral features are well-defined enough to allow a direct comparison, e.g. in the case of the $2.223\;{\rm MeV}$ neutron capture line, yields from inverse reactions are further diminished by the $Z^{2}/A$ dependence of the energy loss rate.

We note that the continuum component from Compton scattering of the $2.223\;{\rm MeV}$ line is much less evident when produced via inverse rather than direct reactions. Although we focus here on total spectra emerging from the atmosphere, FLUKA allows us also to investigate the spatial distribution of particular species or events through the atmosphere. We found that the shorter range of heavier ions, diminished as $Z^{2}/A$, results in a greater proportion of secondary neutrons capturing higher in the atmosphere, above the $\approx 12\;{\rm g/cm^{2}}$ Compton scattering length for $\approx 2\;{\rm MeV}$ photons. This accounts for the relatively-weaker Compton tail from the inverse reactions.

In the case of beams of primary accelerated ions with photospheric composition the total spectrum for energies below $\approx 2\;{\rm MeV}$ and above $\approx 10\;{\rm MeV}$ is dominated by the component due to the direct reactions, while in the range from $\approx 2$ to $10\;{\rm MeV}$ the contributions of the components due to direct and inverse reactions are comparable. One should note that for energies below $\approx 2\;{\rm MeV}$ the component due to the direct reactions is enhanced by the Compton-scattered continuum from the $2.223\;{\rm MeV}$ neutron capture line. In the case of beams of primary accelerated ions with impulsive-flare composition the total spectrum for energies below $\approx 10\;{\rm MeV}$ is dominated by the component due to the inverse reactions, while for energies above $\approx 10\;{\rm MeV}$ the contributions of the components due to direct and inverse reactions are comparable. For the impulsive-flare composition, the contribution from inverse reactions is approximately two orders of magnitude larger than for the photospheric composition, consistent with the substantially enhanced heavy-ion abundances. In both cases, we note that the structures corresponding to the narrow nuclear de-excitation $\gamma$-ray lines in the energy range from $\approx 2$ to $10\;{\rm MeV}$ are very well-defined in the spectrum of the component due to the direct reactions. For energies below $2\;{\rm MeV}$, on the other hand, the spectrum of the component due to the direct reactions is dominated by the continuum emission from the Compton scattering of the $2.223\;{\rm MeV}$ neutron capture line, such that the structures corresponding to the narrow nuclear de-excitation $\gamma$-ray lines become nearly invisible.

Aiming to compare our results with previous calculations, we extracted the yields of the neutron-capture line at $2.223\;{\rm MeV}$ and the 12C de-excitation line at $4.438\;{\rm MeV}$ from the FLUKA simulation for primary accelerated ions with impulsive flare composition. The yield of the $2.223\;{\rm MeV}$ line was obtained from the total, full $\gamma$-ray spectrum. For the yield of the $4.438\;{\rm MeV}$ line, however, we used the direct component of the de-excitation $\gamma$-ray line spectrum, which has better statistics and energy binning resolution. The resulting value obtained for the fluence ratio of the $2.223\;{\rm MeV}$ line to the $4.438\;{\rm MeV}$ 12C line is $8.7$, in remarkably good agreement with the value of $8.4$ calculated by Murphy et al. (2007) for similar model parameters.

Lastly, in Figure \irefF-RMK-FLUKA-Pion-D4 we show the pion decay continuum above $10\;{\rm MeV}$ resulting from direct reactions (assuming $\alpha/{\rm proton}=0.1$ and $\delta=4$). For comparison we also show the corresponding spectral template supplied in OSPEX, based on Murphy, Dermer, and Ramaty (1987). We consider the two spectra to be in reasonable agreement given: 1. the OSPEX template is calculated for primary particles extending to $10\;{\rm GeV/nucleon}$ rather than the $1\;{\rm GeV/nucleon}$ used here; 2. since the FLUKA results are integrated only over the backward hemisphere we reduced the OSPEX spectrum by a factor of 2, as discussed in Section \irefs-model, but some modest anisotropy will still reduce the backward hemisphere spectrum compared to an average over $4\pi$ steradians; 3. the FLUKA calculation includes Compton scattering of photons before they leave the solar atmosphere; 4. the FLUKA spectrum also shows the 15.1 MeV 12C line which is not included in the OSPEX pion decay template. A more detailed discussion of high-energy continuum, directionality, etc. will be provided in MacKinnon et al. (2019).

4 Fits to Observed $\gamma$-Ray Spectra

\ilabel

S-fitting

In this section we present the results for the fit to the $\gamma$-ray emission spectrum of the 2010 June 12 solar flare obtained using templates built from simulations carried out with FLUKA for the full $\gamma$-ray spectrum, which are incorporated into OSPEX. The fit is implemented by combining the FLUKA templates with standard functions available in OSPEX to account for the continuum emission from the bremsstrahlung of primary accelerated electrons. We consider downward isotropic beams of primary accelerated ions with power-law energy distribution of spectral index $\delta=4$ in the range from $1\;{\rm MeV/nucleon}$ to $1\;{\rm GeV/nucleon}$ and impulsive-flare composition (assuming $\alpha/{\rm proton}=0.1$).

Fits to $\gamma$-ray spectra from solar flares have been performed either by using two separate templates for the components of the nuclear de-excitation $\gamma$-ray line spectrum due to the direct reactions (narrow lines) and the inverse reactions (broad lines) or a single template for both components. For the fit to the spectrum of the 2010 June 10 solar flare, Ackermann et al. (2012) used a single template built with the RMK code. Here we choose to perform the fit by using two separate templates. It is important to note that in the case of fits performed with a single template the ratio between the contributions of the direct and inverse components to the total best-fit model spectrum is fixed, while in the case of fits performed with two separate templates the ratio between the contributions is determined by the fitting procedure. One should also note that although the FLUKA spectra lack the contribution to the $511\;{\rm keV}$ line of radioactive positron emitters, we do not add a gaussian line component to the spectrum at this energy. The line is barely evident above the electron bremsstrahlung continuum in the GBM data we will analyze, and inclusion of a further line component would not significantly improve the fit to data.

As pointed out in Section \irefs-model, the photon spectra used to build the FLUKA templates are renormalized to one primary accelerated proton in the energy range from $30\;{\rm MeV}$ to $1\;{\rm GeV}$. Moreover, in order to perform the fit with OSPEX222see https://hesperia.gsfc.nasa.gov/ssw/packages/xray/dbase/nuclear˙template˙normalization˙readme.txt. the FLUKA templates are further renormalized to one photon in the energy range considered, i.e.,

[TABLE]

where $I_{\rm dir}$ and $I_{\rm inv}$ are respectively the yields of the direct and the inverse components, given by:

[TABLE]

The FLUKA templates for the full $\gamma$-ray spectra cover the energy range from $E_{\rm min}^{\rm phot}=100\;{\rm keV}$ to $E_{\rm max}^{\rm phot}=1\;{\rm GeV}$.

The sum of the contributions of the direct and inverse components to the total best-fit model spectrum is given by:

[TABLE]

where $a[0]_{\rm dir}$ and $a[0]_{\rm inv}$ are respectively the normalizations of the templates for the direct and inverse components determined from the fitting procedure, in units of photons ${\rm cm}^{-2}\;s^{-1}$ at Earth. One should note that the ratio between the contributions of the direct and inverse components is given by $a[0]_{\rm dir}/a[0]_{\rm inv}$ and so is determined by the fitting procedure.

Equation \irefsum-inv-dir can be written as:

[TABLE]

where

[TABLE]

This corresponds to the contribution to the total model spectrum one would obtain from a single template built by combining the direct and inverse spectra given in Equations \irefdir-spec and \irefinv-spec with the relative abundances of the primary accelerated heavy ions multiplied by a common factor, i.e., $w_{i}\rightarrow w^{\prime}_{i}=w_{\rm fac}w_{i}$, where

[TABLE]

The total yield of this single template is given by

[TABLE]

It follows from the above discussion that the number of primary accelerated protons at the Sun with energies $>30\;{\rm MeV}$ is given by:

[TABLE]

where $R=1\;{\rm AU}$ and $\Delta t$ is the time interval over which the observed spectrum is integrated. The corresponding numbers of primary accelerated $\alpha$-particles and heavy ions of species $i$ are respectively given by $N_{\alpha}=w_{\alpha}N_{p}$ and $N_{i}=w^{\prime}_{i}N_{p}$.

The GOES M2-class solar flare of 2010 June 12 (SOL2010-06-12T00:57) was observed by the two instruments of the Fermi satellite, the Gamma-ray Burst Monitor (GBM) and the Large Area Telescope (LAT), in the NOAA active region 11081 at the heliographic coordinates N23W43 Ackermann et al. (2012). In our analysis of this event we use the Continuous Spectroscopy (CSPEC) data from the GBM bismuth germanate (BGO) detector Meegan et al. (2009) and the LAT Low Energy (LLE) data from LAT Pelassa et al. (2010).

Using OSPEX, we fit the GBM background-subtracted spectrum in the energy range from $0.3$ to $10\;{\rm MeV}$, accumulated between 00:55:40 and 00:58:50 UT333One should note that in Ackermann et al. (2012) they fit the background-subtracted GBM and LAT spectra integrated over a 50 s time interval (00:55:40 - 00:56:30 UT), but for the delayed neutron capture and annihilation lines they fit the spectra integrated over a 250 s time interval (00:55:40 - 00:59:50 UT). In order to account for the delay of the neutron capture and annihilation lines, we choose to fit the spectra integrated over a single larger time interval of 190 s (00:55:40 - 00:58:50 UT) for all spectral components. and the LAT background-subtracted spectrum in the energy range from $30$ to $300\;{\rm MeV}$ accumulated in the same time interval. We fit these spectra with four components: the two FLUKA templates for the direct and inverse components of the full $\gamma$-ray spectrum and, following the analysis by Ackermann et al. (2012), a single power-law function (1pow) and a power-law function multiplied by an exponential (1pow-exp). The last two components account for the bremsstrahlung of primary accelerated electrons and are both provided as standard fitting components by OSPEX.

Because OSPEX does not currently have the capability to perform a joint simultaneous analysis of data obtained with more than one instrument, we make use of an iterative procedure, fixing some of the parameters using data from one instrument and then refining values of others using data from the other, repeating this to iterate to a single solution using both datasets, characterised by a minimum value of the reduced chi-squared statistic. We make no claim of uniqueness for the data fits obtained below; their main point is to show that FLUKA simulations can provide photon spectra that give a statistically acceptable fit to data across the whole of the $\gamma$-ray range, for reasonable values of fast ion numbers and energy distribution.

We start by choosing the parameters of the function 1pow-exp at values that provide a contribution to the flux similar to that obtained by Ackermann et al. (2012). This first step is necessary because the 1pow-exp component competes in the fitting procedure with the FLUKA template for the component due to the inverse reactions so, if one allows OSPEX to fit all four of the spectral components at once, it attributes too much of what should be the flux from inverse reactions to the 1pow-exp component instead444Similar issues are mentioned by Ackermann et al. (2012) although they use a single nuclear template.. The parameters of the 1pow-exp component are then held fixed throughout the rest of the iterative procedure which consists of the following steps:

i) The GBM spectrum is fitted by varying the normalization parameters of the two FLUKA templates and the parameters of the function 1pow, while keeping fixed the parameters of the function 1pow-exp.

ii) The LAT spectrum is fitted by varying the normalization parameters of the two FLUKA templates, now fixing the parameters of the functions 1pow and 1pow-exp. We allow the normalisation parameters of the two FLUKA templates to vary only within a restricted range around the values found in the previous step because the relative importance of direct and inverse components is much more strongly constrained by GBM data than by LAT.

iii) The GBM spectrum is fitted again by varying the parameters of the function 1pow, while keeping fixed the normalization parameters of the two FLUKA templates obtained in step (ii) and the parameters of the function 1pow-exp.

iv) Steps (ii) and (iii) are repeated iteratively until the best-fit to the combined GBM/LAT spectrum is obtained.

In Figures \irefF-12jun-FLUKA-Full-counts-GBM and \irefF-12jun-FLUKA-Full-counts-LAT we show the results for the fits to the GBM and LAT spectra of the 2010 June 12 solar flare obtained with the FLUKA templates for the full $\gamma$-ray spectrum components. In Figure \irefF-12jun-FLUKA-full we show the result for the best-fit to the combined GBM/LAT photon spectrum. The four best-fitting components are also shown. The model parameters obtained in this fit are listed in Table \ireftab:parametros_ajust12jun_FLUKA_Full together with their estimated uncertainties (as mentioned above the parameters of the 1pow-exp function are adopted at the outset and held fixed during the fitting procedure so we give no uncertainties for them). We note that the parameters of the 1pow component are in reasonable agreement with those found by Ackermann et al. (2012), allowing for the difference in accumulation times. Also shown in the table are the values for the ratio between the contributions of the direct and inverse components, $a[0]_{\rm dir}/a[0]_{\rm inv}$, the common factor that multiplies the relative abundances of the primary accelerated heavy ions, $w_{\rm fac}$, and the number of primary accelerated protons, $N_{p}(>30\;{\rm MeV})$, obtained from the fit.

It is useful to comment separately on the quality of the fit in the GBM and LAT energy ranges. For the GBM spectrum in the energy range from $300\;{\rm keV}$ to $10\;{\rm MeV}$ the fit is very reasonable ($\chi_{r}^{2}=2.89$), particularly for the structures corresponding to the $\gamma$-ray lines produced by the de-excitation of 12C and 16O nuclei. As shown in Figure \irefF-12jun-FLUKA-Full-counts-GBM, the normalization $a[0]_{\rm dir}$ of the FLUKA template for the component due to the direct reactions is fixed primarily by the fit to the structure corresponding to the $2.223\;{\rm MeV}$ neutron-capture line while the normalization $a[0]_{\rm inv}$ of the FLUKA template for the component due to the inverse reactions is determined mainly by the fit to the structures corresponding to the 12C and 16O lines. As discussed above (Section \irefS-Lines) our FLUKA template for the component due to the direct reactions exhibits less well-defined structures for the narrow lines between $\approx 1$ to $2\;{\rm MeV}$ range but the GBM data also do not show well-defined narrow lines in this energy range so there is little adverse impact on the quality of the fit. The predicted $2.223\;{\rm MeV}$ line appears broader than expected but this detail, also found in Ackermann et al. (2012), seems to reflect details of the provided GBM instrument response matrix.

The result obtained for the fit to the LAT spectrum in the photon energy range from $30$ to $300\;{\rm MeV}$ is only satisfactory ($\chi_{r}^{2}=6.19$). As shown in Figure \irefF-12jun-FLUKA-Full-counts-LAT, the fit is dominated by the component due to the direct reactions. For energies below $\approx 200\;{\rm MeV}$ the fit is very good. For energies above $\approx 200\;{\rm MeV}$, on the other hand, the component due to the direct reactions is weaker than that observed in the data, resulting in more significant residuals. The FLUKA templates used in these fits assume a primary ion distribution extending only to $1\;{\rm GeV/nucleon}$. Templates constructed assuming a greater maximum ion energy would certainly fit the LAT data better at higher energies.

Since the components due to the direct and inverse reactions are strongly constrained by the fit to the GBM spectrum in the range from $300\;{\rm keV}$ to $10\;{\rm MeV}$, only two iterations are necessary to obtain the best fit of the combined GBM/LAT spectrum. A full analysis of these data would search for a best fit across several templates calculated using different values of ion energy spectral index $\delta$. The exercise carried out here, with a fixed value $\delta=4$, nonetheless shows that the FLUKA simulations can give a description consistent with GBM and LAT data, across the photon energy range from $300\;{\rm keV}$ to $300\;{\rm MeV}$. Although the strength of the conclusion is limited by the GBM’s spectral resolution, it is remarkable that a single ion spectral index can achieve an acceptable fit to the data across this wide energy range.

5 Summary and Concluding Remarks

\ilabel

S-summary

We presented simulations of the $\gamma$-ray spectrum of solar flares using the Monte Carlo code FLUKA. The nuclear de-excitation $\gamma$-ray line spectra obtained in our FLUKA simulations are in reasonable agreement with the more detailed spectra generated with the RMK code, currently the main tool for analysis of the $\gamma$-ray spectrum of solar flares. Moreover, our FLUKA simulations of the full $\gamma$-ray spectrum are self-consistent, in that the main spectral components in this range (the nuclear de-excitation lines, the $511\;{\rm keV}$ positron annihilation line – excluding, for now, the contribution from the decay of radioactive positron-emitter daughter nuclei – the $2.223\;{\rm MeV}$ neutron-capture line, and the pion decay continuum) are generated from the same distribution of primary accelerated ions. The transport of accelerated ions and photons through a semi-empirical model of the ambient solar atmosphere is taken into account, along with a number of relevant nuclear processes. To the best of our knowledge, this is the first attempt to fit $\gamma$-ray spectra of solar flares using a self-consistent model of the several spectral features in the energy range from $\approx 100$s ${\rm keV}$ to $\approx 100$s ${\rm MeV}$. Previous analysis employed independent functions or templates for the nuclear de-excitation, neutron capture and annihilation lines, and the pion decay continuum. FLUKA’s usefulness for this purpose will be further enhanced when future work includes the contribution of radioactive daughter nuclei to the $511\;{\rm keV}$ positron annihilation line.

FLUKA can be a valuable tool for the study and interpretation of solar flare $\gamma$-ray spectra. Using data from Fermi (GBM and LAT) for the 2010 June 12 solar flare Ackermann et al. (2012), we show that FLUKA gives a statistically acceptable fit to the $\gamma$-ray spectrum (as measured by the reduced chi-square), implying reasonable estimates of ion numbers. It is important to emphasize that in this paper we do not intend to perform a detailed analysis of solar flare $\gamma$-ray data. However we have taken some care in Section \irefS-fitting to show how the spectral components appear, separately and in combination, when compared with data.

FLUKA’s treatment of nuclear collisions, excitation and de-excitation is sufficiently complete to capture all the main features in the photon spectrum (cf. the discussion of GEANT4’s de-excitation line capabilities in Tang and Smith (2010)). In order to completely exploit these capabilities it will be necessary to improve both the statistics and the energy resolution of the photon spectra generated in the FLUKA simulations, a task we will pursue in a future work.

Our results were obtained using the standard FLUKA USRBDX detector that monitors the flux of particles crossing a defined surface. The regularly spaced bins used by USRBDX, combined with the comparatively low cross-sections for some of the de-excitation lines resulted in a less detailed photon spectrum than we might have desired below $\approx 2\;{\rm MeV}$. Via some of the user-editable routines, it is possible to implement e.g. nonuniform bin sizes to improve the quality of the spectrum around the strongest narrow lines.

Since FLUKA evaluates the transport of both particles and photons, it is possible to investigate the directivity of the $\gamma$-ray emission, especially comparing the $2.223\;{\rm MeV}$ neutron capture line and nuclear de-excitation lines. The first is expected to be generated deep in the photosphere and susceptible to strong limb-darkening effect Hua and Lingenfelter (1987), while the latter are expected to originate in the chromosphere, without limb-darkening effects. The versatility of FLUKA opens the possibility to revisit observations of $\gamma$-ray flares from past decades to reveal more information about the elusive presence of accelerated ions during solar flares, and to predict e.g. $\gamma$-ray height structure that could be studied with future instruments with imaging capability.

Each simulated spectrum shown here typically required many days’ processing time on a high-end laptop computer. We relied on standard FLUKA detectors, particularly USRBDX, to obtain these spectra which are specific to the assumed energy and angular distribution of the primary ions and the chemical abundances of the target and accelerated ions. Many projects, for instance a comprehensive exploration of the parameter space, would be computationally challenging or even impossible with this approach. An alternative approach, however, would start from a well-defined set of primary ion energies and directions and store the details of all the photons crossing the boundary between lower atmosphere and corona, together with the details of the primary particle that gave rise to them. Emergent spectra for different injected ion distributions could then be synthesised by applying weights determined by the energies and directions of the primary ions to the photons obtained via full FLUKA simulations. Variations in source or target abundances could be accommodated similarly. Calculating a spectrum via application of these weights would involve only modest computational effort, likely light enough to enable synthesis of spectra on-the-fly for fitting, systematic exploration of the primary ion parameter space, etc. The major computational effort would be in obtaining the basic photon spectra with enough statistical reliability in the first place but this would only need to be done once. The present work has demonstrated that FLUKA simulations can be a valuable tool for interpreting flare $\gamma$-ray data; such a project, to be pursued elsewhere, offers a practical approach to making them routinely useful.

Acknowledgments

The referees’ many comments and questions helped greatly in improving the paper. We thank the Royal Society Newton Fund for supporting our UK-Brazil collaboration through project NI140209 The THz wavelength range as a window on extremes of solar flare particle acceleration. We also thank the financial support from FAPESP under grants 2009/18386-7 and 2017/13282-5. PJAS acknowledges support from the University of Glasgow’s Lord Kelvin Adam Smith Leadership Fellowship. ALM and PJAS acknowledge relevant and helpful discussions with members of the ISSI International Team on Energetic Ions: The Elusive Component of Solar Flares and with participants in the Lorentz Center Workshop on Solar Sources of GeV Gamma-rays, 26 Feb - 2 Mar 2018. DST acknowledges support from CAPES and Instituto Presbiteriano Mackenzie. CGGC research is partially supported from CNPq (Grant 305203/2016-9). We thank the FLUKA team who provide and support the code. Finally, we thank R. J. Murphy for kindly supplying the copy of the RMK code for calculating nuclear de-excitation $\gamma$-ray line spectra we used in this work.

Disclosure of Potential Conflicts of Interest

The authors declare that they have no conflicts of interest.

FLUKA simulations are set up by means of an input ASCII file built by the user which consists of a sequence of command lines, often called “cards” for historical reasons Ferrari et al. (2011). In general, a typical FLUKA input file contains cards through which the following simulation elements are defined:

• •

Particle source properties: characteristics of the beam of primary particles such as particle type, energy and angular distribution, starting position);

• •

Geometry: solid bodies and surfaces combined by boolean operations to implement the complete partition of the space of interest into regions;

• •

Materials: single-elements or compounds (either pre-defined or user-defined) which are assigned to the geometry regions;

• •

Detectors: estimators used to score physical quantities of interest;

• •

Settings: parameters, conditions and general directives which determine how calculations are performed, e.g. specification of production and transport thresholds, biasing schemes and physical effects to be included.

FLUKA can accurately simulate the transport and the interactions of about $60$ different types of particles, including electrons and muons with energies from $1\;{\rm keV}$ to $1000\;{\rm TeV}$, photons with energies from $100\;{\rm eV}$ to $10000\;{\rm TeV}$, hadrons with energies from $1\;{\rm keV}$ to $10000\;{\rm TeV}$, and all corresponding antiparticles, as well as neutrinos, low-energy neutrons (below $20\;{\rm MeV}$) down to thermal energies and heavy ions with energies up to $10000\;{\rm TeV/nucleon}$. Interactions are implemented by robust and updated physics-models, continuously benchmarked and optimized against experimental data at single interaction level. The models are fully integrated, such as to provide a consistent treatment of all physical processes to the same level of accuracy.

The models used in FLUKA to implement the transport and interactions of electromagnetic particles and muons cover a comprehensive range of processes, including bremsstrahlung, Compton and Rayleigh scattering, pair production, positron annihilation in flight and at rest, Bhabha and Moller scattering, photoelectric effect, photonuclear interactions and photomuon production. Transport and ionization energy losses are performed through approaches shared by all charged particles (electrons, positrons, muons, charged hadrons and heavy ions). Transport is implemented through an original algorithm Ferrari et al. (1992) based on the Moliere theory of multiple scattering improved by Bethe (1953), supplemented by a single scattering algorithm based on the Rutherford formula. The treatment of ionization energy losses is based on the Bethe-Bloch theory Bethe and Heitler (1934), supplemented with ionization potentials and density effect parameters taken according to the compilation of Sternheimer, Berger, and Seltzer (1984) and with shell corrections derived from the parametrized formula by Ziegler and Andersen (1977). Ionization fluctuations are implemented through an original approach (alternative to the standard ones based on the Landau and Vavilov theory) which makes use of general statistical properties of the cumulants of a distribution Fassò et al. (1997). The approach allows to combine explicit $\delta$-ray production and transport with ionization fluctuations and also to include corrections for spin-relativistic effects and distant collisions.

The treatment of hadronic interactions in FLUKA is based on a microscopic approach which makes use of several models properly tailored to the different energy ranges Ferrari and Sala (1998); Battistoni et al. (2015). At energies below $5\;{\rm GeV}$, inelastic hadron-hadron interactions are implemented through the isobar model based on the resonance production and decay of particles, while elastic and charge exchange interactions are implemented through phase-shift analysis and eikonal approximation555At energies below the pion production threshold ($\approx 290\;{\rm MeV}$ for nucleon-nucleon interactions and $\approx 170\;{\rm MeV}$ for pion-nucleon interactions), only elastic and charge exchange scattering processes take place.. At energies in the range from $5$ to $20\;{\rm TeV}$, inelastic hadron-hadron interactions are implemented through a modified version of the Dual Parton Model (DPM) Capella et al. (1994) coupled to a hadronization scheme. Inelastic hadron-nucleus interactions at energies from reaction threshold up to $20\;{\rm TeV}$ are implemented through the FLUKA model called PEANUT (Pre-Equilibrium Approach to Nuclear Thermalization) Ferrari and Sala (1994); Fassò et al. (1994); Ferrari and Sala (1998); Battistoni et al. (2006).

The PEANUT model describes the interactions of hadrons with nuclei as a sequence of the following steps:

• •

Glauber-Gribov multiple scattering cascade with formation zone;

• •

Generalized intra-nuclear cascade;

• •

Pre-equilibrium emission stage (exciton-based);

• •

Equilibrium stage: evaporation, fission, Fermi break-up, $\gamma$ de-excitation.

The $\gamma$ de-excitation process is particularly important for our simulations of $\gamma$-ray spectra from solar flares. After the evaporation process, the residual excitation energy of nuclei is dissipated through the emission of nuclear de-excitation $\gamma$-ray photons. At high excitation energies the cascade of $\gamma$-ray transitions is implemented through a statistical model by assuming a continuous density of nuclear levels, while at low excitation energies (below an arbitrary threshold) a database of tabulated discrete experimental nuclear levels is used Ferrari et al. (1996).

The extension from hadron-nucleus to nucleus-nucleus interactions is implemented by modified versions of three external event generators linked to FLUKA, namely the Boltzmann Master Equation (BME) model Cerutti et al. (2006) at energies below $0.125\;{\rm GeV/nucleon}$, the Relativistic Quantum Molecular Dynamics (RQMD-2.4) model Sorge, Stöcker, and Greiner (1989); Andersen et al. (2004) at energies in the range from $0.125$ to $5\;{\rm GeV/nucleon}$, and the Dual Parton Model and Jets (DPMJET-III) model Roesler, Engel, and Ranft (2001) at energies above $5\;{\rm GeV/nucleon}$. All three external event generators are interfaced with the PEANUT model in order to handle the equilibrium stage processes.

We end this section with a brief description of the treatment given in FLUKA to the low-energy neutrons (i.e., those with energies in the range from thermal to $20\;{\rm MeV}$), which is also relevant for our simulations of $\gamma$-ray spectra from solar flares since it provides the emission of the $2.223\;{\rm MeV}$ neutron capture line. Transport and interactions of low-energy neutrons are implemented in FLUKA through a multigroup algorithm Ferrari et al. (2011); Battistoni et al. (2015). The scattering matrices for the reaction channels are calculated from neutron cross section data provided by a dedicated library in which the energy range of interest is divided into $260$ neutron energy groups ($31$ of which are thermal) and $42$ gamma energy groups. Elastic and inelastic reactions are simulated by group-to-group transfer probabilities. In particular, the production of $\gamma$-ray photons by low-energy neutrons (including the $2.223\;{\rm MeV}$ neutron capture line) is treated by evaluating a downscattering matrix which provides the probability for a neutron in a given group to generate a photon in each of the gamma groups666One should note that the multigroup algorithm is used only for gamma generation. The transport of $\gamma$-ray photons produced by low-energy neutrons is performed through the same modules used in FLUKA for the transport of photons produced by other processes..