FLUKA Simulations of Pion Decay Gamma-radiation from Energetic Flare Ions

Gamma-ray continuum at>10 MeV photon energy yields information on>0.2 - 0.3 GeV/nucleon ions at the Sun. We use the general-purpose Monte Carlo code FLUKA (FLUktuierende KAskade) to model the transport of ions injected into thick and thin target sources, the nuclear processes that give rise to pions and other secondaries and the escape of the resulting photons from the atmosphere. We give examples of photon spectra calculated with a range of different assumptions about the primary ion velocity distribution and the source region. We show that FLUKA gives results for pion decay photon emissivity in agreement with previous treatments. Through the directionality of secondary products, as well as Compton scattering and pair production of photons prior to escaping the Sun, the predicted spectrum depends significantly on the viewing angle. Details of the photon spectrum in the 100 MeV range may constrain the angular distribution of primary ions and the depths at which they interact. We display a set of thick-target spectra produced making various assumptions about the incident ion energy and angular distribution and the viewing angle. If ions are very strongly beamed downward, or ion energies do not extend much above 1 GeV/nucleon, the photon spectrum is highly insensitive to details of the ion distribution. Under the simplest assumptions, flares observed near disc centre should not display significant radiation above 1 GeV photon energy. We give an example application to Fermi Large Area Telescope data from the flare of 12 June 2010.


Introduction
Gamma-rays from solar flares carry information on energetic ion populations at the Sun (e.g. Vilmer, MacKinnon, and Hurford, 2011). In particular, the 100 MeV pion decay spectral component results from the highest energy flare ions and was first observed from solar flares with the Gamma-Ray Spectrometer (GRS) instrument on the Solar Maximum Mission (SMM) (Forrest et al., 1985). This radiation often shows a distinctive spectral flattening around m π 0 c 2 /2 = 67 MeV (here m π 0 is the neutral pion rest mass) that is inconsistent with bremsstrahlung. It thus testifies unambigously to the presence at the Sun of accelerated ions whose energies exceed 200 -300 MeV/nucleon. Observations with the Fermi Large Area Telescope (LAT - Atwood et al., 2009) have shown that it occurs more frequently than previously appreciated, even in smaller Mclass flares, that it can extend well into the GeV photon energy range and that it sometimes persists for many hours after all the other features of a flare have concluded (Ackermann et al., 2012(Ackermann et al., , 2014Share et al., 2018).
Ions whose energies exceed the relevant thresholds can produce both charged and neutral pions. The resulting (angle-averaged) gamma-ray spectrum was discussed in detail by Murphy, Dermer, and Ramaty (1987) and its depth-and angle-dependence by Mandzhavidze and Ramaty (1992) and Tang and Smith (2010). The neutral pions decay immediately to two photons giving a spectrum with a maximum at 67 MeV, symmetric in log ǫ (ǫ is the photon energy), with a width determined by the energies of the primary ions (Stecker, 1970;Dermer, 1986b). Charged pions decay to muons and then to electrons and positrons which emit gamma-rays by bremsstrahlung as they slow down (with a few % additional contribution from positron annihilation in flight). The relative proportions of π 0 and π ± are determined once the properties of the primary ions (chemical composition, energy distribution) and of the source region (chemical composition, thick or thin) are specified but the bremsstrahlung spectral component may be reduced relative to that from π 0 decay if synchrotron losses of the e ± are significant. Because the primary particles are positively charged, e − are significantly less numerous than e + .
The highest energy particles accelerated in association with solar flares pose particular challenges to theory. These are compounded for the long-lasting events observed by LAT Omodei et al., 2018;Share et al., 2018). Some role for a shock wave in the interplanetary medium seems likely. Radio observations offer some support for this idea (Gopalswamy et al., 2018) but questions persist, particularly around the feasibility of shock-accelerated particles making their way back to the Sun (e.g. Hudson, 2018). The absence of a strong correlation between the number of particles measured in space and the number deduced at the Sun (de Nolfo et al., 2019) is a further complication. By modelling the formation of the gamma-ray continuum in detail we may hope to constrain the angular and energy distribution of the primary ions and the properties of the medium where they interact, and thus to constrain particle acceleration models.
The FLUktuierende KAskade (FLUKA) code (Ferrari et al., 2005;Böhlen et al., 2014) is a general purpose Monte Carlo code for simulating the passage of energetic particles in matter, including the production and propagation of secondaries. It includes detailed, physically based, well validated models for the directions and energies of secondaries produced in hadronic interactions. One defines one or more geometrical shapes and their chemical composition, and the spatial form, species and velocity distribution of the primary particles. FLUKA follows test particles and their secondaries throughout the whole of the defined volume(s). One may set "detectors" to carry out various tasks, e.g. determine the energy deposited in a particular volume, or the flux of particles of a particular type crossing the boundary between two regions. It was introduced as a tool for modelling solar flare gamma-radiation by Tusnski et al. (2019). Here we will use FLUKA to predict the gamma-ray continuum radiation 10 MeV resulting from hypothesised energetic ion distributions interacting with the solar atmosphere and illustrate how these simulated spectra may be used to interpret observed gamma-ray spectra..
FLUKA is "not a toolkit". Its authors aim in the first instance to provide treatments of all the processes taking place, to the same degree of precision, rather than to offer the user a set of disconnected modules each providing its own treatment of a particular physical process. As a consequence it can be used to simulate a single, self-consistent gamma-ray spectrum extending from the positron annihilation line at 0.511 MeV through the nuclear deexcitation line component at ≈ 1−7 MeV, to the pion decay continuum in the GeV energy range (Tusnski et al., 2019). Nonetheless mechanisms are provided to, e.g., suppress some of the full set of processes taking place, even if this is strictly unphysical, and these prove very useful in developing physical understanding.
In the next section we recall some basic kinematics of pion production in p-p collisions, highlighting the likely importance of the angular distribution for gamma-ray spectra. In Section 3 we make some comments on the use of FLUKA for astrophysical modelling and describe how we set up a model that can be used to simulate solar flare gamma-ray spectra. The core of Section 4 is a set of FLUKA-simulated thick target gamma-ray spectra for a range of assumed ion energy and angular distributions. Before this we further illustrate some of the main physical processes. In Section 5 we compare some of the simulated spectra with Fermi LAT data of a solar flare. Section 6 gives some concluding remarks.

p-p Collisions
To see the potential importance of the ion angular distribution for the gammaray spectrum, we recall the threshold energy calculation for π 0 production in p-p collisions. Suppose an energetic proton with mass m p and velocity v = βc collides with another proton at rest. The centre of momentum (CM) frame has where γ = (1−β 2 ) −1/2 and β is the proton speed in units of the speed of light, as usual. The total CM energy is E CM = 2 (γ + 1)m p c 2 . The proton energy at π 0 production threshold is determined by the requirement E CM = (2m p + m π 0 ) c 2 and the π 0 is produced with γ = 1.072 in the laboratory (lab) frame, moving in the same direction as the incident proton -not at rest, even although we consider the case of threshold. π 0 decay, isotropic in the CM frame, will take on some modest beaming along the incident proton direction when viewed in the lab frame. An observer in the backward direction will be more likely to see red-shifted photons and the gamma-ray spectrum will exhibit an asymmetry to lower photon energies. When it decays, the two emitted photons will in general be Doppler shifted, dependent on the direction of emission. In the extreme case that they are emitted parallel or anti-parallel to the pion lab frame velocity vector, their lab frame energies will be 43 MeV and 91 MeV. Above threshold more energy is available in the CM frame for sharing among the products of the reaction but there will still be a preference for the pions to be produced at a small angle to the incident proton. When primary ion distributions are isotropic, or if we average over the viewing angle, the π 0 decay spectrum is expected to be symmetric in log ǫ about m π 0 c 2 /2 = 67 MeV. In contrast, an ion distribution beamed away from the observer will produce a spectrum with an excess of red-shifted photons.
Secondary electrons and positrons too will exhibit angular distributions that reflect the ion angular distribution. Their energy distribution displays a maximum at an energy of ≈ 35 MeV (Dermer, 1986a;Murphy, Dermer, and Ramaty, 1987), i.e. a Lorentz γ ≃ 70. Their bremsstrahlung radiation has a characteristic angular width of ≈ 1/γ 1 • so again we can expect anisotropy of the emitting particles to have an important influence on the detected radiation.
To further illustrate the effects of directionality, Figure 1 shows the angleaveraged photon spectrum emitted into the backward hemisphere when 350 MeV protons are injected into a hydrogen target (more simulation details are discussed below, in Section 3). Two different initial angular distributions are assumed: beamed, unidirectional downwards, and isotropic in the downward hemisphere. The resulting spectra have some notable features. While some e ± bremsstrahlung continuum is produced, such low energy ions do not produce these secondaries prolifically and the π 0 component of the spectrum may be clearly distinguished, particularly in view of its narrowness. The peak of the π 0 spectrum is Doppler  Figure 1. Gamma-ray spectrum produced in the backward hemisphere from protons of initial energy 350 MeV directed downwards into a thick target source. The initial angular distribution is: unidirectional, parallel to the local vertical (yellow curve), or isotropic in the downward hemisphere (blue). The spectra are normalised to one primary proton.
shifted from 67 MeV to more like 55 MeV, a value consistent with lab frame π 0 energies as above. The peak is significantly sharper in the downward beamed case, with a factor of ≈ 5 fewer photons above 60 MeV. Effects of anisotropy are maximised just above threshold but we see in this example the potential of gamma-ray spectra to diagnose the primary ion angular distribution.

Depth of Production
Ion energy distributions certainly extend on occasion into the GeV energy range. A 1 GeV proton has a vertical stopping depth of 159 g cm −2 of hydrogen (Berger et al., 2005), far beyond the ≈ 4 g cm −2 depth of the photosphere (e.g. Vernazza, Avrett, and Loeser, 1981). Magnetic mirroring will prevent some ions from reaching such depths but any ions that can precipitate freely will emit pion decay secondary photons throughout most of this range. Considering 100 MeV photons for illustrative purposes, we find that the mean free paths against Compton scattering and pair production are 233 and 265 g cm −2 , respectively (Berger et al., 2010). These are comparable to the depths from which many secondary photons originate so radiation will escape, but with a spectrum modified to some extent by these processes. Both processes imply the production of secondary, energetic electrons throughout some region of the atmosphere surrounding the field lines on which primary ions precipitate, whose bremsstrahlung will further modify the total photon spectrum escaping from the atmosphere.

FLUKA Simulations
To set up a FLUKA simulation one prescribes the geometry, as a set of regions of various shapes, specifies the chemical composition of each region, fixes the parameters of the incident fast particle distributions and decides on a total number N tot of fast particles to be simulated. All of this information is encoded in an input file.
For example, in most of Section 4 below we used a very simple, plane-parallel, two-region model of the solar atmosphere, divided in two by the plane z = 0. z measures distance downwards, into the Sun. In the region z < 0 we have a vacuum "corona". In z > 0 there is a uniform density "atmosphere" of solar (Asplund et al., 2009) composition. While a stratified density structure was implemented in Tusnski et al. (2019), total thick target yields are insensitive to the detailed density structure so this simple model is adequate for our purposes here. The total depth of the atmosphere in z and its lateral extent in x and y were made big enough to ensure all particles stop and produce their thick target yields within it.
As described in Tusnski et al. (2019) we wrote a custom subroutine to simulate particles coming from a power-law energy distribution between two fixed energies E min and E max : Here H is the Heaviside step function. The quantity N 0 takes care of normalisation; with the functional form fixed, FLUKA automatically normalises results to one primary particle. We generally take a value of E min determined by the threshold energy for the reactions of interest (e.g. 270 MeV for pion production in p-p collisions). E max and δ are parameters of the acceleration mechanism whose influence on the emergent photon spectra will be studied. It seems likely that E max varies between events. At the upper extreme, the 26 February 1956 event produced signals at locations on Earth where the geomagnetic cutoff rigidity was ≈ 17 GV (Rishbeth, Shea, and Smart, 2009). The ground level enhancement (GLE) event of 29 September 1989 evidently involved particles at Earth of up to 25 GV rigidity (Swinson and Shea, 1990) or possibly even higher (Miroshnichenko, De Koning, and Perez-Enriquez, 2000). For the 24 May 1990 event, Debrunner et al. (1997) showed from energetic neutron timing that E max ≈ 2 GeV. Alexander, Dunphy, and MacKinnon (1994) argued that SMM/GRS spectra in the > 10 MeV range were best interpreted in terms of E max in the 0.5-0.8 GeV range. Thus the E max values found at the Sun could lie in a fairly wide range, from 1 to > 10 GeV. To study its consequences here we parametrise this maximum energy via a sharp cutoff at E max , rather than, e.g., an exponential roll-off so that any curvature in the simulated spectra is clearly due to radiation processes and transport, rather than the departure of the ion distribution from a power law.
We define µ = cos θ, where θ is the angle between the ion velocity vector and the local (downward) vertical. Thus µ = 1 corresponds to an ion directed vertically downwards and µ = 0 to an ion that starts out parallel to the solar surface. The angular distribution of primary ions is described by the function M (µ). We implemented four possible forms: i.e. unidirectional vertically downwards (we denote the Dirac delta function by δ D to avoid confusion with the energy spectral index), downward isotropic, moderately beamed downward, or a moderate "pancake" distribution, concentrated at large angles to the vertical. Once started, FLUKA generates N tot incident particles with energies drawn from the distribution given in Equation 2 and directions µ drawn from whichever form of M (µ) and follows them and their secondaries until they stop or, if they escape the dense atmosphere into the region z < 0, reach the boundary of the box.
Like all Monte Carlo codes, FLUKA can be made to output very large quantities of data describing details of the evolution of all particles. In this work, however, we followed the FLUKA authors' recommendation and implemented standard "detectors", devices in the code that monitor all particles crossing the boundary between two regions, or evaluate the energy deposited in a particular volume by simulation particles. We set such a detector to detect photons passing through the boundary plane z = 0 from the dense atmosphere to the corona. This gives us the spectrum of photons emitted into the backward hemisphere, rapidly and easily at the expense of some more detailed information about each photon. In the next few sections we integrate over all directions but it is also possible to study directionality by binning photons in solid angle, count statistics allowing.
FLUKA carries out several repetitions of the specified Monte Carlo simulation to enable error estimation on derived quantities. We started the simulations described here by taking five repetitions (FLUKA default) with N tot = 6.4 × 10 5 each time, if necessary then increasing N tot and/or carrying out further repetitions so that the error in almost all photon energy bins was less than 3% (except possibly where fluxes get very low at the highest energies).

Thin Target
We first show some results for a thin target, a source which primary particles traverse with negligible energy degradation. We study initially isotropic ions released at the centre of a spherical target, aiming to eliminate features of the photon spectra that arise from particle anisotropy (cf. Section 2.1). Resulting photon spectra may be compared with previous calculations (Murphy, Dermer, and Ramaty, 1987;Alexander, Dunphy, and MacKinnon, 1994;Vilmer et al., 2003) and may also be relevant to ions contained in a low-density region.
We consider a sphere of hydrogen. Primary ions start from the origin of the sphere (r = 0), with an isotropic angular distribution and a specified energy distribution. The sphere has uniform density n H and a total radius R. Primary ions scatter negligibly in angle (e.g. Spitzer, 1956) so every ion in the calculation will traverse a column depth negligibly different from n H R. We integrate photons and other escaping secondaries over the surface of the sphere and include all  Figure 2. FLUKA calculation of the gamma-ray spectrum produced by 1 GeV protons released isotropically at the centre of a sphere of hydrogen whose radius corresponds to a column depth of 0.084 g cm −2 (blue curve). For comparison the yellow curve shows the photon flux expected from the same situation calculated using the cross section given in Murphy, Dermer, and Ramaty (1987).
escaping particles, so that the results are effectively integrated over angle to the sphere surface normal.
In Figure 2 we show the photon spectrum escaping from such a source with n H R = 0.084 g cm −2 . All protons have an initial energy of 1 GeV. With this value of n H R they lose ≈ 4 × 10 −4 of their energy before leaving the target. The spectrum is dominated by π 0 decay photons with the expected form: symmetric in log ǫ, with a maximum at 67 MeV (e.g. Stecker, 1970). Also shown is the thin target flux expected from this situation calculated using the cross section for photon production via π 0 decay given in Murphy, Dermer, and Ramaty (1987). The two independent calculations agree to within a factor of order unity, as we should expect in these spherically symmetric conditions. The result found using their results is slightly greater, by up to about 30%, which we attribute to the different cross-section calculations of the two treatments. Similar agreement is obtained for other proton energies.
The source is too thin for secondary e ± to reveal themselves via bremsstrahlung but these particles are present nonetheless. Figure 3 shows the electron and positron distributions that cross the boundary of the sphere. The positron distribution ( Fig. 3a) has the form expected from π + → µ + → e + decay, roughly symmetric around a maximum at ≈ 35 MeV (Murphy, Dermer, and Ramaty, 1987). The sphere becomes a thick target for e ± of around 1 MeV so the distribution at the lowest energies has been modified by energy loss en route to the surface.
In Fig. 3b we show the total escaping electron distribution as well as the contributions produced via decay of each of π 0 , π + and π − separately. The threshold energy for π − production is 870 MeV and the cross section is still very small at 1 GeV (Dermer, 1986a). Electrons are produced via π − decay only at about 10 −6 of the level of positrons. Instead the Dalitz decay π 0 → γe + e − , which occurs for 1.2% of π 0 (Olive and Particle Data Group, 2014), accounts for most of the electrons above ≈ 10 MeV. For electrons (panel b) we show the total electron distribution (thick, black curve), and the contributions via decay of: π − (blue), π 0 (red), π + (green).
Below ≈ 3.3 MeV the distribution is dominated by knock-on electrons, i.e. ambient electrons scattered to higher energies by primary protons. Consider an electron initially at rest, gaining energy at the expense of a proton moving initially at speed βc. Conserving energy and momentum (e.g. Jackson, 1999;Abraham, Brunstein, and Cline, 1966) we find the maximum energy ∆E gained by the electron to be where γ = (1 − β 2 ) −1/2 as usual. This expression is very close to the correct value as long as γ ≪ m p /m e . For a 1 GeV primary proton Equation 4 gives a maximum knock-on electron energy of 3.3 MeV, as seen in Fig. 3b. The dominance of knock-on secondaries over pion decay products at low energies has already been discussed in the cosmic ray context by Abraham, Brunstein, and Cline (1966). Implications of knock-on electrons for studying flare ions will be discussed elsewhere. Even π + decay products can contribute to the electron flux. Positrons will (Bhabha) scatter ambient electrons before leaving the sphere. Because positrons are so much more numerous this population may be comparable to electrons produced by π 0 or π − decay. We see from Fig. 3b that they are the dominant contribution to the flux leaving the sphere between 3.3 and ≈ 10 MeV.
The relative magnitudes of the various components depend on specifics of this situation: the thickness of the sphere and the primary proton energy. In particular, once we consider protons further above the threshold for π − production we would expect secondary electrons via this channel to dominate those resulting from Bhabha scattering by secondary positrons. Nonetheless this example first reassures us that results from FLUKA are comparable to those obtained using other calculational approaches, and begins to illustrate the modifications introduced by transport in the source. . Angle-integrated gamma-ray spectrum, per primary proton, produced in a spherical thick target source (blue curve). 1 GeV protons are released isotropically from the centre of a sphere of hydrogen with n H R = 160 g cm −2 . The other two curves show the photon spectrum without the contribution from secondary e ± (orange), and additionally with Compton scattering and pair production suppressed (green).

Isotropic Thick Target
With the aim of comparison with previous calculations, we next simulate a thick target source with the same spherical symmetry. Again we release an isotropic distribution of 1 GeV protons in a spherical target at r = 0. Now n H and R are chosen such that n H R = 160 g cm −2 , slightly greater than the stopping depth of a 1 GeV proton in hydrogen (158.7 g cm −2 , Berger et al., 2005). Again we show the angle-integrated energy distribution of photons crossing the boundary of the sphere, in Fig. 4.
Photons are present to ≈ 700 MeV. The spectrum that emerges from the target steepens above 100 MeV but shows no pronounced π 0 "bump". For understanding we also show the spectra that result with certain physical processes artificially suppressed. First we set a very high threshold for the creation of secondary e ± , so that only π 0 decay photons are present. Now the spectrum below 100 MeV is flat, almost energy independent. We might expect such a spectrum if the source was threaded by a very high magnetic field so that e ± lifetimes were suppressed. Lastly we also suppress Compton scattering and pair production, processes that modify the escaping spectrum. Only in this case do we reproduce the thick target π 0 spectra of previous authors (Murphy, Dermer, and Ramaty, 1987;Alexander, Dunphy, and MacKinnon, 1994;Vilmer et al., 2003), with a maximum at ≈ 67 MeV.
We contrived this spherically symmetric situation for purposes of comparison with previous calculations. The primaries start out at the maximum possible column density, a situation we do not expect to occur on the Sun. Nonetheless we can conclude that FLUKA gives results in agreement with previous calculations, and that protons in the GeV energy range may produce observable radiation in spite of their great stopping depths, but with a spectrum significantly modified by pair production and Compton scattering.

Downward Directed Primary Ions: Power Law
Here we show a set of simulated backward hemisphere gamma-ray spectra from primary protons injected with power-law energy distributions and the angular distributions specified in Equations 3 above. We concentrate on the continuum spectrum, above 10 MeV; line features at lower energies were discussed by Tusnski et al. (2019). For now spectra are averaged over the viewing angle but we discuss directionality in the next section. We use the same geometry as above and also now include He, C, N, and O in the composition of the target, with the relative abundances of Asplund et al. (2009). Figure 5 shows photon spectra resulting from primary protons injected into this source with a downward isotropic distribution. Each panel is for a fixed value of δ: 2, 3, 4 or 5 (panels a, b, c, and d respectively). Within each panel we show four spectra, calculated using different values of E max = 1, 2, 5 or 8 GeV. Higher-energy protons have a greater probability to produce pions so the net photon yield per primary increases with E max , and decreases with δ. As E max increases, the photon spectrum extends to higher photon energies and the spectrum above ≈ 100 MeV exhibits more power-law character.
With E max = 1 GeV the form of the photon spectrum barely depends on δ: proton energies do not extend far enough above threshold for the spectrum to bear an imprint of the proton energy distribution. In Figure 6 we bring  Figure 6. Gamma-ray spectrum produced in the backward hemisphere from power-law distributed protons between 0.27 MeV and 8 GeV, with δ = 2, 3, 4, 5, isotropic in the downward hemisphere. The spectra are normalised to one primary proton.  Figure 7. Gamma-ray spectra produced in the backward hemisphere from protons with δ = 3, Emax = 5 GeV and the four angular distributions of Equations 3: pancake distribution (red), downward isotropic (blue), proportional to µ (green), and unidirectional vertically down (gold).
together the four spectra calculated with E max = 8 GeV to show more clearly the steepening of the photon spectrum with increasing δ. The spectrum has an approximate power-law character ≈ ǫ −η above about 100 MeV, steepening more rapidly as the upper cutoff is approached. In the photon energy range 0.1 -2 GeV, with this value of E max an increase of 1 in δ leads to an increase of ≈ 0.5 in η.
In Figure 7 we fix δ = 3 and E max = 5 GeV and calculate emergent photon spectra for each of the four primary angular distributions in Equations 3. As expected we see that a greater proportion of downward-directed primaries leads to fewer escaping photons, with a spectrum falling off more steeply with photon energy. In the limiting case M (µ) = δ D (µ − 1) the escaping spectrum falls off very steeply and cuts off at a much lower energy than the other cases. In this extremely beamed case the escaping photon spectrum is almost completely insensitive to the details of the primary proton distribution, never extending above ≈ 700 MeV. We illustrate this further in Fig. 8, showing results for this form of M (µ) with E max = 1 or 8 GeV, and δ = 2 or 5. While the cases with δ = 2 produce a greater photon flux and the cases with E max = 8 GeV extend to slightly higher photon energies, it would be difficult observationally  Figure 8. Gamma-ray spectra produced in the backward hemisphere from protons with a unidirectional, vertically downward angular distribution, with δ = 2 and Emax = 1 GeV (blue), δ = 2 and Emax = 8 GeV (gold), δ = 5 and Emax = 1 GeV (green), δ = 5 and Emax = 8 GeV (red).
to distinguish between these spectra, despite the very different ion distributions that produce them.

Photon Spectrum Dependence on Viewing Angle
The FLUKA detectors measure fluxes of particular particles species at boundaries or within defined volumes. We used the detector called USRBDX to obtain photon distributions averaged over the whole of the backward hemisphere, shown in the previous section, and also broken down into intervals of solid angle referred to the normal to the boundary (in our case, between the dense atmosphere and the corona). By default we divide the backward hemisphere, 0 < θ < 90 • into ten equal intervals of solid angle. As a first example, Fig. 7 shows photon spectra calculated assuming an isotropic downward primary proton distribution with δ = 3 and E max = 5 and 8 GeV, averaged over viewing angles θ of 0 • − 26 • and 84 • − 90 • . If the magnetic field is vertical these values correspond to a flare at disc centre and at the limb, respectively and we use these words to refer to the two cases although the true situation may be more complex. We see pronounced differences between the two cases: the disc centre spectrum exhibits no photons above ≈ 1 GeV and no influence of E max , similarly to the unidirectional downward cases of Fig. 8, whereas the limb spectra exhibit power-law tails extending to ≈ 4 and 7 GeV. However the spectral flattening below ≈ 100 MeV expected from π 0 decay is more pronounced for the disc centre spectra.
Disc centre spectra are dominated by the directionality of secondaries, together with transport of photons produced deeper in the atmosphere. Figure 10 emphasises this, showing disc centre spectra with a variety of primary proton distributions. Only the normalisation of the photon spectrum varies significantly across these various cases and it would be difficult or impossible to discriminate between them observationally, even although δ, E max , and M (µ) are quite different. As a random example, note that the photon spectrum produced when δ = 2, E max = 1 Gev and a pancake angular distribution is almost identical, up to about 200 MeV photon energy, with the photon spectrum produced with δ = 4, E max = 8 Gev and a downward isotropic distribution. Figure 9. Gamma-ray spectra produced from protons with a downward isotropic angular distribution and δ = 3, for flares at disc centre with Emax = 5 GeV (blue) and 8 GeV (green), and at the limb with Emax = 5 GeV (gold) and 8 GeV (red). At disc centre the form of the photon spectrum is very insensitive to the form of the primary proton distribution. At the limb the photon spectrum extends into the GeV energy range, increasing with E max and with hardness determined by the value of δ. As the viewing angle increases from 0 • (disk centre) to 90 • (limb) there is a smooth transition between these two sorts of spectrum, as shown in Fig. 11.

Flare of 10 June 2010
The Geostationary Operational Environmental Satellite (GOES) M2 class flare of 12 June 2010 (00:50 UT) was observed by the Large Area Telescope (LAT) of the NASA Fermi satellite (Ackermann et al., 2012). High count rates in the anti-coincidence shields meant that the flare was not detected using the standard LAT data analysis method but it was possible to obtain useful data using the LAT Low Energy (LLE) procedure (detailed in the Appendix of Ajello et al., Photons consistent with a solar origin were detected in the LAT for about 60s, a much short duration than many other, highly extended events observed by LAT (e.g. Ajello et al., 2014;Share et al., 2018). In the extended events it seems likely that ions will be trapped in low density regions by either or both of magnetic inhomogeneity (mirroring) and plasma turbulence, physical processes not included in our FLUKA simulations. Even in the brief 12 June 2010 event these effects could still be important. For example, a 1 GeV proton mirroring at the 4 g cm −2 depth of the photosphere needs ≈ 40 traversals of the loop to stop, a time of about 1.5 s for a 10 9 cm loop length (neglecting any pitch-angle scattering) -much less than the event duration. However this short duration also means that radiation from ions trapped in the corona will be unimportant for the observed spectrum which is why we select this particular event for a first trial of FLUKA spectra against data. The emergent gamma-ray spectrum is determined by the energies and directions of the fast ions as they precipitate into the dense atmosphere, however complex their histories in the corona may have been. Tusnski et al. (2019) previously showed that a FLUKA angle-averaged spectrum can provide an acceptable fit to these data, giving results similar to those found using the templates included in the standard data analysis software Object SPectral EXecutive (OSPEX) (Tolbert, 2020). Here we show additionally that the LAT data can be used to discriminate between different assumptions for primary ion energy and angle distributions. We neglect the contribution from primary α-particles, straightforward to include but yielding little extra insight for present purposes.
If magnetic field lines were all oriented vertical to the local solar surface, the viewing angle would be simply determined by the heliocentric angle but the real situation is likely to be more complex. Smith et al. (2003) found that red shifts of nuclear deexcitation lines in the flare of 23 July 2002 were greater than expected on the basis of its heliocentric angle. Hudson et al. (2020), using the Bifrost MHD model, found that field lines may meander substantially in the chromosphere, presenting a wide range of angles to any particular line of sight. Thus we also allow the viewing angle to vary, rather than fixing it based on the heliocentric angle.
Following Ackermann et al. (2012) we fit the spectrum in the photon energy range 30 -300 MeV, having subtracted an instrumental background calculated by interpolation between time intervals before and after the flare. We start by assuming a single power-law photon spectrum I(ǫ) = Aǫ −s and find a best-fit spectral index s = 1.93, in good agreement with that found by Ackermann et al. (2012) (Figure 12). We also show an example of a fit to a spectrum including only a FLUKA pion decay template, calculated assuming δ = 3, E max = 1 GeV, downward isotropic primary ions and a viewing angle (i.e. angle between the line of sight and the downward vertical) in the range 45.6 • −53.1 • , appropriate to the heliocentric angle of the flare. Only the amplitude of the spectrum, equivalent to the total number of protons, is determined by fitting the LAT data. The normalised χ 2 values from these two fits are comparable, 1.63 (power-law) vs. 1.65 (pion decay). In both cases the pattern of residuals highlights systematic discrepancies between the model spectrum and the data: the observed spectrum has a curvature inconsistent with a single power-law; and the steeply falling spectrum in the ∼ 30 − 60 MeV photon energy range is inconsistent with any of the pion decay templates. Thus we next fit the LAT data using the sum of a pion decay template and a single power law. In each case the parameters of the primary proton distribution are fixed and we find the best-fit values of three parameters: the amplitude of the pion-decay template (i.e. the total number of protons) and the parameters A and s of the power-law component.
We found best fits to the LAT spectrum for a range of primary ion parameters. Results are given in Table 5. In Figure 13 we show the fit that gave the smallest χ 2 (0.71) from all those we tried, obtained for δ = 3, E max = 1 GeV, downward isotropic ions and a viewing angle of 0 • − 25.8 • . The power-law component is important in the 30-60 MeV energy range but the pion decay template dominates above this; modulating the power-law component with an exponential rollover e −ǫ/ǫ roll does not improve the fit, whether the rollover energy ǫ roll is fixed e.g. at 100 MeV or allowed to vary freely along with the other parameters. Also shown in Figure 13 is the best fit obtained using the pion decay template supplied in OSPEX for δ = 3 together with a power-law photon spectrum. The template also assumes a magnetic field strength of 300 G (determining the e ± synchrotron loss rate), and ambient density of 10 15 cm −3 . Based on a homogeneous source with an isotropic distribution of ions, the OSPEX template resembles the FLUKA angleaveraged spectra given above ( Figure 5). As measured by χ 2 the fit is poorer than the best-fit FLUKA spectrum (1.70 vs. 0.71), showing the importance of retaining particularly the angular dependence of the emission. The OSPEX template fit implies just 4.1 × 10 27 protons > 30 MeV, three orders of magnitude less than found with the FLUKA fit, as a result of the greater weight given by OSPEX to the power-law spectral component. As shown also in the distribution of residuals, the angle-averaged OSPEX template has less of the curvature needed to play a substantial role in fitting the observed spectrum; cf. the fits shown in Tusnski et al. (2019).
Many other primary proton parameters give fits that are equally acceptable, statistically, but some sets of parameters are clearly disfavoured and χ 2 values show informative systematic trends. From Table 5 we can make some general statements: • The lowest χ 2 values are obtained with δ = 3, 4. Assuming δ = 2 or 5 leads to χ 2 significantly greater than one, irrespective of other assumptions made • Assume that primary protons are injected into the source with an isotropic distribution (as might be expected for many particle acceleration mechanisms -e.g. Melrose, 1974). Keeping δ and E max fixed, the best fits (i.e. lowest χ 2 ) are then obtained for viewing angles close to 0 • , rather than the flare heliocentric angle. χ 2 increases monotonically with viewing angle. Also E max should be comparatively low, 1 or 2 GeV (as in Fig. 13). • If we now suppose that ions are injected either unidirectionally downwards, or more modestly beamed ∼ µ, the lowest χ 2 values with fixed δ and E max are generally found in the region of the flare location. χ 2 < 1 may now be obtained with E max = 5 or 8 GeV. • Viewing angles of > 70 • are always less consistent with the LAT data: χ 2 values are greater and the power-law component plays a greater role throughout the photon energy range being fit.

Discussion and Conclusions
FLUKA can be a useful tool for modelling ion transport, pion production and the observable consequences in flares. The directionality of the secondary products is modelled, as well as the transport of the escaping photons (cf. Share et al., 2018). The aim of the FLUKA authors to provide a consistent treatment of all relevant processes means that it also includes, e.g., the contribution of "knock-on" electrons, primary proton bremsstrahlung, photons, and other secondaries produced via the decay of heavier, strange mesons and baryons (e.g. K mesons, Λ's; these are produced for ion energies 2.5 GeV but make only minor contributions to the observed fluxes), etc. Our discussion of the 12 June 2010 flare shows that the spectrum can constrain the primary ion energy and angular distribution. In the case of this particular flare, too many primary ions directed towards the line of sight results in a Table 1. Parameters for various assumed primary ion distributions and viewing angles, as found from fitting the 21 June 2010 LAT spectrum in OSPEX. Np is the number of protons above 30 MeV assuming the power law continues smoothly to that energy. Ions are assumed downward isotropic except when the spectral parameters Emax and δ are marked with: an asterisk (*), in which case ions are unidirectional vertically downward; or a dagger †, in which case they are moderately beamed, ∼ cos θ. The top line gives the best fit single power law, with no contribution from ions, and a few lines give Np and χ 2 when no power-law component is included. spectrum that is harder than that observed. The best fits avoid this either by directing primary ions mostly downwards (or at least away from the observer if field lines are not vertical) or by viewing the flare from close to the effective downward vertical direction and minimising the maximum primary ion energy. The capacity to model the angular dependence of the emission, as provided by FLUKA or similar codes, is clearly essential to making such statements. In particular angle-averaged spectra resemble those seen when the emitting region is viewed from a large angle to the vertical: generally harder, extending to higher photon energies, with a less pronounced π 0 feature (cf. Fig. 9). In reality we expect a more complex situation in which the direction of the field as well as its strength and gradient vary both with height and laterally in the chromosphere (Hudson et al., 2020). The magnetic mirror force will change the velocity vectors of ions as they slow down in the source, in a way not included in our FLUKA simulations. Apart from viewing-angle effects, ions in long-duration events may be trapped in low density regions. In such a situation π 0 decay radiation will still be produced instantaneously but there will be less accompanying radiation from secondary e ± because their lifetimes are much longer in low density conditions, and/or they are suppressed by synchrotron energy losses (Murphy, Dermer, and Ramaty, 1987). Clearly our FLUKA simulations yield constraints that are useful in their own right but how to satisfy these constraints may need further consideration employing a detailed model of the chromospheric and photospheric magnetic field.
We concentrated on a short event less likely to involve radiation from particles trapped in low density regions. The forthcoming Fermi-LAT Solar Flare Catalog (M. Pesce-Rollins, private communication) emphasises "impulsive" or "gradual" as a fundamental distinction between high-energy events. It will be interesting to see if the data analysis approach here, applied across many of these events, highlights systematic differences between these two basic classes. Subsequent work will obtain best fits to observed spectra for many events, aiming to exploit the diagnostic potential discussed here, as well as considering other sorts of emission from secondary particles, e.g. (gyro)synchrotron radiation that should peak in the THz range. This has important consequences for diagnostics using high frequency observatories such as the Atacama Large Millimeter Array (ALMA) (Tuneu et al., 2017;Wedemeyer et al., 2016).
• FLUKA treats densities < 10 −10 g cm −3 as vacuum. This is the density of the chromosphere at a height of ≈ 600 km in the VAL-C (Vernazza, Avrett, and Loeser, 1981) model, so a literal model of the solar atmosphere is not feasible. This is not a major problem for present purposes since the radiative ouput from all of a thick target, or from a layer of fixed column density, is what is important, irrespective of its physical size -as long as this is much larger than the muon decay length, and as long as we also pay attention to the following. • FLUKA models the slowing down of ions assuming a neutral medium. The Coulomb logarithm may be significantly larger (typically by a factor of up to three) in an ionised medium (Hayakawa and Kitao, 1956;MacKinnon and Toner, 2003;Trottet et al., 2015). Most previous calculations of solar pion decay radiation employ neutral medium stopping rates, appropriate to the deep atmosphere (e.g. Murphy, Dermer, and Ramaty, 1987;Tang and Smith, 2010). We continue with this practice here. • However, ionisation potentials used to calculate the Bethe-Bloch stopping rate (Olive and Particle Data Group, 2014, Section 32) are by default those appropriate to molecules, e.g. H 2 , O 2 , etc., slightly different from the atomic values. Different values appropriate to single atoms may be inserted by hand. Since the ionisation potential appears in the argument of the logarithm the difference is slight in any case. Attempting to model energy loss in an ionised medium by substituting hν p for an atomic ionisation potential is too substantial a departure for FLUKA to accommodate, however. • FLUKA includes the modifications to the slowing-down rate that result from the collective response of a neutral medium (e.g. Weaver and Westphal, 2002;Olive and Particle Data Group, 2014). The consequence is that the yield of secondaries from a specified target depends not only on the total column depth, but slightly on the actual density encountered. This correction becomes most important at the highest energies but is only ever of order unity. We checked that our results did not depend on density to a degree that would be important for interpreting observed fluxes. However, the parametric (Sternheimer) form used to represent this effect breaks down, giving runtime errors, at densities slightly higher than the lower limit of 10 −10 g cm −3 mentioned above; another factor that discourages a literal model of the solar atmosphere.
• Sometimes we want information on secondaries resulting from unlikely events (e.g. Figure 1). Then statistically reliable spectra are obtained more easily by activating "biasing" in FLUKA, specifically an artificial reduction in the mean free path for hadronic interactions.