Issue 
A&A
Volume 654, October 2021



Article Number  A82  
Number of page(s)  11  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202141183  
Published online  15 October 2021 
Modelling magnetised medium particle transport in the guiding centre limit with GEANT4^{⋆}
^{1}
Universidade Presbiteriana Mackenzie, Centro de Rádio Astronomia e Astrofísica Mackenzie, São Paulo, Brazil
email: jordituneu@pm.me
^{2}
Instituto de Astronomía y Física del Espacio, CONICETUBA, Buenos Aires, Argentina
^{3}
University of Glasgow, School of Physics and Astronomy, Glasgow, UK
Received:
26
April
2021
Accepted:
2
August
2021
Monte Carlo codes are a standard tool for studying energetic particle propagation, secondary production, and radiation in astrophysical settings. In magnetised plasmas such as those found in solar active regions, the enormous disparity between particle gyroradii and system scales proves to be a major computational obstacle. To address this problem we have written a new module in Geant4 using the guiding centre (GC) approach in which the particle motion is averaged over a gyrofrequency. We describe the formulation and implementation of this method in particular dealing with the uncertainty in gyrophase so that particle velocities are welldefined for input to the modules handling reactions. As far as feasible, we compare the propagation and slowing down of primary protons, secondary particle production, and run times in the GC limit with the Newton–Lorentz approach, finding very good agreement between the two methods and orders of magnitude improvement in run times in the GC case. Finally, we present an illustrative solar physics application involving two interacting dipoles, which is only achievable using the GC approach.
Key words: astroparticle physics / methods: numerical / Sun: flares / Sun: magnetic fields / Sun: particle emission
The Geant4 codes used in this work, along with running instructions, are available from a public GitHub repository (https://github.com/guigue/Geant4_Codes_for_Solar_Simulations/tree/master).
© ESO 2021
1. Introduction
Particle acceleration is a key feature of solar and stellar flares (Miller et al. 1997) and many other astrophysical phenomena (e.g., Hurley et al. 2005; Merloni & Fabian 2001). The transport of primary accelerated particles and the secondaries they produce then becomes an important topic for deducing properties and behaviour of the energy release process, understanding energy deposition in and radiation of the flaring atmosphere, among other topics. In the example of solar flares, energetic ions and electrons embody a large fraction of all the energy released (Emslie et al. 2012). They are believed to deposit this energy in the lower atmosphere, driving flare gradual phase phenomena. Their properties at the Sun are deduced from radio, hard Xray, and gammaray observations combined with modelling of transport and secondary production. In the case of energetic ions nuclear reactions give rise to gammaray lines in the 0.5–8 MeV range through nuclear deexcitation and continuum in the ≳100 MeV range from decay of secondary π^{0} (Murphy et al. 2007; Vilmer et al. 2011) as well as other observable radiation, so treatments of the production of these secondaries are particularly important for understanding observed radiation. Monte Carlo codes are an obvious tool for these problems since they can deal with arbitrary geometries and chemical abundances, and any reaction for which cross sections are known. Monte Carlo has been the modelling tool of choice for bespoke codes developed within the solar physics community (see Murphy et al. 2007). However several mature powerful codes rooted in laboratory nuclear and particle physics also exist. These include physically motivated, wellvalidated models for secondary production and have a great potential as tools for astrophysical modelling. Example applications already exist. Kotoku et al. (2007) used Geant4 (Agostinelli et al. 2003; Allison et al. 2006) to simulate electron transport in solar flares and the resulting bremsstrahlung gammaray continuum spectra. Ion transport and secondary production in solar flares have also been modelled by Tang & Smith (2010) using Geant4, and by Tusnski et al. (2019) and MacKinnon et al. (2020) using Fluka (Böhlen et al. 2014; Ferrari et al. 2005). Geant4 has also been used to study cosmic ray interaction with planetary atmospheres and surfaces (Desorgher et al. 2006; Matthiä et al. 2016; Guo et al. 2018). Fermi Large Area Telescope (LAT) observations of gammaradiation from the quiet sun (Abdo et al. 2011) have motivated studies of cosmic ray interaction with the solar atmosphere using Fluka (Mazziotta et al. 2020) and Geant4 (Li et al. 2016).
Magnetic fields play an important role in all these situations, and their influence on particle propagation needs to be included in realistic modelling. In the low field strengths found near planets or in the interplanetary medium it is feasible to do this by direct (numerical) integration of the Lorentz force, given a model for the spacedependence of the magnetic field Desorgher et al. (2006). This approach is impractical in situations where the particle gyroradius is much smaller than characteristic system scales, however, as it is in solar magnetic active regions. The particle gyroradius is given by r_{L} = γmcv_{⊥}/(qB), where m and q are the particle mass and charge, respectively, c is the speed of light, B is the magnitude of the magnetic field, v_{⊥} is the magnitude of the particle velocity perpendicular to the magnetic field, and γ is the particle total energy in units of its rest mass. For example, a 1 GeV proton in a 500 G field has a gyroradius of ∼10^{4} cm. This is to be compared with a typical active region lengthscale of 10^{9} cm. To explicitly integrate the Lorentz force the code has to take several time steps per gyroperiod to correctly describe gyration, or millions of time steps to follow a single particle as it traverses an active region loop once. If the particle is trapped by magnetic mirroring just above the photosphere, one traversal of the loop will involve the particle encountering ∼8 g cm^{−2} of material, to be compared with its stopping depth of ∼160 g cm^{−2} (Berger et al. 2005). We see that sufficiently accurate description of gyromotion involves significant computational effort in regions where very few nuclear interactions occur. Run times rapidly become impractically long if even spatially integrated photon spectra are to be obtained with adequate statistics.
It is normal in flare physics to employ a guiding centre (GC) description, in which particle motion is averaged over a gyroperiod (e.g., Northrop 1963; Parks 2004). The motion of the centre of the gyratory motion is studied instead of the precise position of the particle. This approach is valid if the gyroradius r_{L} ≪ B/∇B, where B is the magnitude of the magnetic field and ∇B is its gradient. With this condition the magnetic moment , where p_{⊥} is the magnitude of the particle momentum perpendicular to the magnetic field, is an adiabatic invariant of the motion; there is no need to model gyration explicitly; and the equations of motion assume a simplified form describing the constancy of the magnetic moment and the rates of drift across the field due to curvature and ∇B. However, GC gives a reduced description, in terms of the two components of v parallel and perpendicular to B. No information on gyrophase is retained.
Here we present a Geant4 application implementing a GC description of particle motion in the presence of a magnetic field. To do this we have to implement the GC equations of motion in the necessary coordinate system, to ascertain that a particle does indeed traverse its full gyrating path length, and to deal with the uncertainty in gyrophase so that particle velocities are well defined for input to the modules handling reactions. Section 2 addresses these matters, having first recalled the main points of a GC approach. Section 3 tests the GC approach by comparing results, as far as is computationally feasible, with those obtained using the full Lorentz force, and describes a simple illustrative application. We make concluding comments in Sect. 4.
2. Guiding centre approach in Geant4
2.1. Guiding centre basics
The motion of a relativistic charged particle in the presence of a magnetic field is described by the Newton–Lorentz (NL) equation (Jackson 1999)
where v and B(r) are the particle velocity and magnetic field vectors, respectively. As usual γ = (1−v^{2}/c^{2})^{−1/2}. In a homogeneous magnetic field B(r) = B_{0}, and the motion of the particle is wellknown: rectilinear motion parallel to B_{0} and oscillatory circular motion with radius r_{L}.
We write where hats denote unit vectors. The GC formalism is usually written in terms of the components of particle velocity parallel and perpendicular to the field
where . The GC approach (Northrop 1963; Ozturk 2012) applies as long as
Then the magnetic moment is an adiabatic invariant, giving a second constant of the motion along with the energy of the particle. The position of the particle is written as
where R is the GC position vector (centre of the gyroorbit) and ρ is the gyration radius vector, lying in the plane perpendicular to R and having magnitude r_{L}. The GC description concentrates on the motion of R and gives the position to within a gyroradius.
Averaging Eq. (1) over one gyroperiod we obtain the GC equations of motion, which are often written as follows:
Equation (5) expresses motion parallel to B combined with a slow drift perpendicular. With v_{∥} determined by Eq. (6), the magnitude of v_{⊥} is immediately obtained from conservation of particle energy. The phase defining its precise direction in the plane perpendicular to B plays no role in the GC description, however, and remains undefined. The field direction vector also changes continually so Eq. (6) does not on its own give components of dv/dt. Both of these issues must be dealt with to construct a GC description compatible with Geant4.
2.2. Guiding centre implementation in Geant4
Geant4 already includes implementations of several algorithms for numerical solution of equations of motion in the general form:
The force F is normally given by Eq. (1), but this is not a requirement, so our approach is to code the form appropriate to the GC description. First we identify position r with the guiding centre R, so from Eq. (5):
To evaluate the force F we introduce two unit vectors and in the plane perpendicular to B. These satisfy the following,
and we discuss their explicit construction below. Then we can write
and so:
The real numbers δ and ϵ satisfy δ^{2} + ϵ^{2} = 1. They may be interpreted as the sine and cosine of the gyrophase. In Eq. (12) the parallel and perpendicular components of p change because of the adiabatic invariance of the mangetic moment while the vectors , , and change because of the geometry of B.
For consistency with the rest of Geant4 we must work in a particular coordinate system. Once this system has been chosen, however, the ordinary differential equations for the components of r and p may be written down explicitly and coded from Eqs. (9) and (12).
Because particle energy is constant, dp^{2}/dt = 0, and so
From the constancy of magnetic moment we have
since B only changes for the particle because it is moving. We write v_{GC} here for the velocity of the guiding centre, given by Eq. (5). In the same way,
We still need to specify how to calculate the derivative , and in particular to choose the gyrophase (i.e. the quantities δ and ϵ). The Monte Carlo code calculates the distance to the next discrete interaction, at which point all three components of p are required. In the tenuous media of the solar atmosphere the path length to the next interaction will always correspond to very many gyrations so that the gyrophase is in practice a random variable. One possible approach would thus be to choose the gyrophase randomly just before the interaction, generating three components of p for input to the code handling the details of the reaction and the generation of secondaries. We chose not to do this, however, since it would involve modification of core elements of Geant4. Instead, we calculated δ and ϵ when the particle is initialised and fixed them throughout its subsequent evolution. As long as care is taken to initialise particles with random gyrophase, as expected from any plausible astrophysical initial state, the gyrophases at interaction will also be random. Testing establishes that no artefacts arise from this procedure.
In the testing and examples in the next sections we work in cartesian geometry with basis vectors , , and . Then we can construct and as
Equations (16) and (17) fail if at any position B becomes parallel to to within working precision. In this case we substitute or for , amounting to no more than a rotation of the already random gyrophase. Then finally we have
which can be substituted into Eq. (12).
Equations (9)–(18) now constitute equations of motion in the GC limit, which can be implemented as a method within Geant4 and used with any of the numerical routines included for ordinary differential equation solutions. This constitutes the core of our GC extension of Geant4: classes that give the components of dr/dt and dp/dt in the GC limit that can be used with one of the existing numerical integrators, plus code that gives B(r) and the derivatives of its components at any location, specifies a model of atmospheric density, and so on. We invoke the RungeKutta (RK4) ordinary differential equations integrator with default settings for precision, so that the step size is variable, determined by the need to attain maximum acceptable tolerances. The compromise between the rectilinear steps of a Monte Carlo and the curvilinear trajectory of a particle in a magnetic field involves further accuracy parameters. For the sake of comparison, we set them to the same values in both GC and NL approaches. We note, however, that in the GC approach we may have a faster solution still with adequate accuracy if we use a different set of parameters that would not be well suited to the NL approach. The results shown in the next section are all obtained using RK4.
It is possible to specify B(r) either analytically or as a 3D grid of values. In the latter case we found it expedient also to calculate all the derivatives ∂B_{x}/∂x, ∂B_{x}/∂y,...∂B_{z}/∂y, ∂B_{z}/∂z at the vertices of the grid and have them available for calculation of the gradients of , , and . These are approximated by simple differencing:
The adequacy of this approach is shown in the next section, where particles are found to follow field lines closely, as expected.
2.3. Path length
In Geant4 the equations of motion are expressed as derivatives with respect to s, the distance travelled by the particle,
so Eq. (5) becomes
Once the Monte Carlo method has calculated the distance to the next interaction, the equations of motion are integrated to determine where this occurs. The crossfield drift term will normally be negligible in any one step so we can neglect it to note
This says that the path actually traversed by the particle is greater than the distance the guiding centre has travelled along the field, by secθ where θ is the particle pitch angle. Thus, the additional distance traversed by the particle in gyration is automatically included in, for example, the slowingdown rates and mean free paths. In the next section we check whether this in fact happens as desired.
2.4. Physics list
The Geant4 physics list is the description of physics models used in simulations. The selection or construction of a physics list is an important ingredient of a Geant4 simulation. The user decides how much detail is needed for the simulation and which models are better for the user purposes. It is also important to take into account the detail required against CPU performance.
Geant4 offers several options for the physics list, but it is also possible to create a specific physics list adapting it to the user’s needs. Tang & Smith (2010) previously considered the choice of physics list for solar applications. We follow them and use a Geant4 built physics list HadronPhysics QGSP_BERT_HP, where QGSP refers to the Quarkgluon String model (Folger & Wellisch 2003); BERT refers to the Bertini cascade model (Bertini 1963; Wright & Kelsey 2015), recommended for energies ∼10 GeV; and HP refers to the High Precision neutron model which is important if results involving neutrons are expected, as they are here.
3. Model testing and simulation results
The main goal of this work is to demonstrate the accuracy of the GC approach in Monte Carlo simulations primarily by comparison with results obtained using the NL approach, and the dramatic increase in throughput with respect to the classical NL framework. Moreover, we are now in the condition to simulate more realistic flare situations to produce diagnostics for the observational data.
3.1. Solar flare model
We carried out our tests within the framework of a model in which primary accelerated protons are injected into a target with characteristics similar to those of the ambient solar atmosphere. We adopted a simple planeparallel geometry for the vertical structure of the ambient solar atmosphere. This geometry is adequate for a single active region, usually with a scale much smaller than that of the solar radius.
We consider a cubic box centred at the origin of a cartesian coordinate system (Ox, Oy, Oz) with edges of length 2L = 1.3 × 10^{4} km and faces perpendicular to the coordinate axes. The zcoordinate corresponds to the height in the ambient solar atmosphere. For scales approaching a solar radius we would use spherical polar coordinates. The box is filled with a neutral material for which we assume a typical solar atmosphere composition with ^{4}He, C, N, and O abundances relative to H given by Asplund et al. (2009). Our assumption of a neutral medium is appropriate to the photosphere and consistent with previous treatments of gammaray line production (Murphy et al. 2007; Trottet et al. 2015). A xyplane at z = 0 divides the box into two halfspaces. The halfspace at z > 0 represents the coronal region for which we assume a density ρ = 10^{−13} g cm^{−3}. The halfspace at z < 0 represents the chromospheric–photospheric region for which we assume a density ρ = 4 × 10^{−7} g cm^{−3}, obtaining a column density ρL ≃ 255 g cm^{−2}, which is greater than the stopping depth (or range) R = 158.7 g cm^{−2} of a 1 GeV proton in a hydrogen target (Berger et al. 2005). The primary accelerated protons are injected into the model solar atmosphere from a point or an extended source located in the coronal region. This highly simplified density structure is adequate for consideration of the pion decay gammaray continuum (cf. MacKinnon et al. 2020), although a more realistic density structure may be implemented if necessary, as in Tusnski et al. (2019).
We performed simulations for single primary protons released in a given direction and beams of a large number of primary protons with an isotropic angular distribution. In the case of beams we consider that the primary protons are either monoenergetic or have a powerlaw energy distribution over a range E_{min} − E_{max}:
Here E is the kinetic energy of the primary proton, δ is the spectral index, H is the Heaviside step function, and N_{0} is a normalisation constant defined such that
The algorithm implemented by the Geant4 package 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 lowenergy threshold for transport. Geant4 includes robust and well validated physicsbased models for electromagnetic, hadronic, and nuclear interactions that provide a detailed treatment of the physical processes most relevant for simulations of particle transport and interactions in solar flares.
We performed test simulations considering both uniform and dipole magnetic fields. The uniform field is assumed to be parallel to the zaxis
where B_{0} is the (constant) magnitude of the magnetic field at any point of the model space.
The general expression for the magnetic field produced by a magnetic dipole is given by
where r_{d} and μ are the dipole position and magnetic moment vectors, respectively. As shown in Fig. 1, we consider a dipole with magnetic moment parallel to the xaxis, , placed at a point . The Cartesian components of the magnetic field vector B(r) are then given by
Fig. 1. Geometry of the dipole magnetic field model. 
where μ is the magnitude of the magnetic moment.
In order to fix the value of μ we choose a reference magnetic field line that lies in the xzplane and is assumed to cross the xyplane at the footpoints, as shown in Fig. 1. Hence,
where B_{f} is the magnitude of the magnetic field at the footpoints and R_{f} is half of the distance between the footpoints. Thus, the dipole magnetic field B(r) at any point r of the model space is uniquely determined by three free parameters, namely d_{z}, R_{f}, and B_{f}.
3.2. Propagation and slowing down in the guiding centre limit
In order to verify the equivalence between the slowingdown rates of particles in simulations carried out with the GC and NL approaches we compared the evolution of energy with depth, obtained in the two cases for a single primary proton which is released from the coronal region and precipitate into the chromospheric–photospheric region. For these simulations we used 4 × Intel®core™ i77500 CPU @ 2.70 GHz machine with 16 GB RAM. The system runs over a Manjaro Linux distribution where we installed Geant4 version 10.06. We first ran simulations considering particle motion in a uniform magnetic field parallel to the zaxis with magnitude B_{0} = 500 G, for which there is no crossfield drift motion or mirroring. We considered a single primary proton released from the point x = 0, y = 0 and z = 100 km with initial kinetic energy E_{k0} = 0.25 GeV and three different values of initial pitch angle: θ = 15°, θ = 35° and θ = 80°. In a second test, we considered a single primary proton released from the same point with initial pitch angle θ = 35° and three different values of initial kinetic energy: E_{k0} = 0.5 GeV, E_{k0} = 1.0 GeV, and E_{k0} = 5.0 GeV.
In Fig. 2 we compare the primary proton trajectories projected in the xzplane obtained with the GC and NL approaches. The trajectories followed in the GC approach are consistent with the results of the exact NL approach to within a gyroradius, as expected. Differences of detail appear because the algorithm for incorporating consequences of multiple scattering has different consequences in the GC and NL cases, but these differences remain below the scale of a gyroradius. The stopping depths in both cases are consistent with the expected value of ∼400 km (Berger et al. 2005), scaling with cos θ. These figures confirm that the full column depth traversed by the particle is incorporated in the GC simulation, as anticipated in Sect. 2.3, even though the distance traversed by the guiding centre itself is smaller. In the right panels we show the results obtained for the second test, where we considered primary protons with an initial pitch angle θ = 35° and three different values of initial kinetic energy. The stoppingdepth obtained for the primary proton with E_{k0} = 0.5 GeV is larger than that obtained for the primary proton with E_{k0} = 0.25 GeV and the same initial pitch angle. For primary protons with E_{k0} = 1.0 GeV and E_{k0} = 5.0 GeV the stopping depths decrease, being even smaller for E_{k0} = 5.0 GeV. However these are not true stopping depths, but the point at which the primary particle loses its identity due to an inelastic collision, as discussed immediately below. It should be noted that the scale set for the xaxis is orders of magnitude smaller than the scale set for the zaxis, such that for the NL trajectories the change of the gyroradius r_{L} with the height z can be assessed, while the individual gyrations are not properly shown in the figures, particularly in the coronal region where the primary proton moves faster and nearly without interactions. As expected, the gyroradius decreases as the primary proton travels down into the chromospheric–photospheric region, losing energy by interactions with the ambient medium.
Fig. 2. Trajectories of a primary proton in a uniform magnetic field obtained in simulations with the GC and NL approaches. Left panels: initial kinetic energy E_{k0} = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panels: initial pitch angle θ = 35° and initial kinetic energies E_{k0} = 0.5 GeV (A), E_{k0} = 1.0 GeV (B) and E_{k0} = 5.0 GeV (C). 
In Fig. 3 we show the plots for the height z as a function of the ratio E_{k}/E_{k0} (where E_{k} is the instantaneous kinetic energy at height z), which provides a more quantitative comparison between the energy losses obtained in the simulations with the GC and NL approaches. Primary protons with initial kinetic energy E_{k0} = 0.25 GeV continuously lose all their energy, coming to rest at stopping depths which increase with the initial pitch angle (left panel). We verified from the output data that the energy losses are only due to ionisation collisions. Moreover, we can see that the slowingdown rates are nearly the same for simulations with the GC and NL approaches. In the right panel of Fig. 3 we show the results obtained for the second test. A primary proton with E_{k0} = 0.5 GeV (curve labelled A in the figure) loses its energy continuously, coming to rest at a stoppingdepth that is larger than that obtained for a primary proton with E_{k0} = 0.25 GeV and the same initial pitch angle (curves B and C). Primary protons with E_{k0} = 1.0 GeV and E_{k0} = 5.0 GeV are above the threshold for pion production and the probability of the particle undergoing an inelastic collision increases with energy. After such an event the particle has lost its identity and is no longer primary for Geant4. Thus, for higher initial energies we can only carry out the comparison of the NL and GC trajectories until the first inelastic collision. Curves B and C of the right panel of Fig. 3 show such a comparison. We verified from the output data that the trajectories terminate before losing all their energy because of an inelastic collision. Until this happens the results of the NL and GC calculations are again in good agreement.
Fig. 3. Comparison between the energy losses of a primary proton in a uniform magnetic field obtained in simulations with the GC and NL approaches. Left panel: initial kinetic energy E_{k0} = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panel: initial pitch angle θ = 35° and initial kinetic energies E_{k0} = 0.5 GeV (A), E_{k0} = 1.0 GeV (B) and E_{k0} = 5.0 GeV (C). 
Fig. 4. Trajectories of a primary proton with initial kinetic energy E_{k0} = 1.0 GeV in a dipole magnetic field obtained in simulations with the GC and NL approaches. Left panels: initial pitch angle θ = 15°. Right panels: initial pitch angle θ = 35°. 
We also compare trajectories and slowing down obtained in simulations with the GC and NL approaches for particle motion in a dipole magnetic field. For this purpose we consider a reference magnetic field line defined by the following parameters: B_{f} = 1000 G, d_{z} = 0.012 R_{s}, and R_{f} = 0.010 Rs (where R_{s} = 6.96 × 10^{5} km is the solar radius). The top of the reference magnetic field line is located at a height h = 10, 623 km above the xyplane at z = 0 and the magnitude of the magnetic field at this point is B_{h} = 121.5 G. The losscone angle (the minimum pitch angle for a particle to be trapped in the lowdensity coronal region) is then given by sin^{2}(θ_{c}) = B_{h}/B_{f}, which yields θ_{c} = 20.4°. We consider a single primary proton released from the top of the reference magnetic field line with initial kinetic energy E_{k0} = 1.0 GeV and two different values of initial pitch angle, θ = 15° and θ = 35°, which lie inside and outside the losscone, respectively.
In Fig. 4 we compare the primary proton trajectories projected in the planes xz, yz, and xy obtained with the GC and NL approaches. As expected, for an initial pitch angle θ = 15° (left panels) the primary proton precipitates into the chromospheric–photospheric region, losing energy by interactions with the ambient medium and stopping at a given height z, while for an initial pitch angle θ = 35° (right panels) it is trapped in the coronal region, bouncing back and forth between the mirroring points. In both cases the crossfield drift motion can be seen. We also observe that the stopping depths (in the case of θ = 15°) and the mirroring points (in the case of θ = 35°) are nearly the same for the simulations with the GC and the NL approaches.
In Fig. 5 we show the plots for the height z as a function of the ratio E_{k}/E_{k0} obtained in simulations with the GC and NL approaches for a primary proton with initial pitch angle θ = 15° and three different values of the initial kinetic energy: E_{k0} = 0.5 GeV, E_{k0} = 1.0 GeV and E_{k0} = 5.0 GeV. As can be seen, the results are similar to those obtained in the simulations considering a uniform magnetic field: the primary proton with E_{k0} = 0.5 GeV (curve A) continuously loses all its energy by ionisation collisions until it comes to rest, while the primary protons with E_{k0} = 1.0 GeV (curve B) and E_{k0} = 5.0 GeV (curve C) continuously lose energy only up to a given height z where they undergo an inelastic collision. Moreover, in all cases the slowingdown rates are nearly the same for simulations with the GC and NL approaches.
Fig. 5. Comparison between the energy losses of a primary proton in a dipole magnetic field obtained in simulations with the GC and NL approaches. The plots are for an initial pitch angle θ = 15° and initial kinetic energies E_{k0} = 0.5 GeV (A), E_{k0} = 1.0 GeV (B), and E_{k0} = 5.0 GeV (C). 
In Tables 1 and 2 we compare the run times obtained in the simulations with the GC and NL approaches for single primary protons in a uniform magnetic field and a dipole magnetic field, respectively. The simulations with the two approaches were carried out using the same parameters for RK4 and the same random seeds. The GC approach shows a superior performance especially at high energies. Once charged secondary and even tertiary particles are produced, the run time ratio of GC to NL increases faster than linearly. Moreover, we note the strong increase in NL run time for a pitch angle of 15° and energies in the GeV range due to the secondary electrons and positrons that have gyroradius of the order of centimetre scales. The smallest run time is that obtained for a primary proton in a dipole magnetic field with E_{k0} = 1.0 GeV and θ = 35° because we set a stopping condition; specifically, the simulation ends when the total track length of a particle trapped in the coronal region reaches 10^{5} km. This is necessary to prevent particles trapped in the coronal region from bouncing back and forth between the mirroring points indefinitely since the energy losses are negligible.
Run times for simulations of a primary proton in a uniform magnetic field with the GC and NL approaches.
Run times for simulations of a primary proton in a dipole magnetic field with the GC and NL approaches.
3.3. Secondary production in the guiding centre
A crucial test to validate our Geant4 application implementing a GC description of particle transport in magnetic fields is to verify whether the production of secondaries it provides in the simulations is statistically consistent with that obtained with the NL approach. To perform such test we run simulations with both approaches for primary protons in a dipole magnetic field defined by the same parameters used in the simulations carried out in Sect. 3.2. For these simulations we use a node of our cluster Superserver H8DGTHIBQF composed of two 64 bits CPUs AMD MagnyCours with a total of 24 cores, running at 2.2 GHz clock speed, and with 12 GB DDR3 ECC RAM. The system runs over a standard Debian 9 (Stretch) distribution where we installed Geant4 version 10.05.
In order to obtain reliable statistics for secondary production we run five simulations with different random seeds, each considering a beam of 1000 monoenergetic primary protons with initial kinetic energy of 1.0 GeV released from the top of the reference magnetic field line. We assume a random distribution of initial pitch angles lying inside the loss cone (whose angle is θ_{c} = 20.4°), such that all primary protons precipitate into the chromospheric–photospheric region. We set up ‘detectors’ (i.e. software devices, not simulated physical devices) to record the number of secondary e^{−}, e^{+}, π^{0}, π^{−}, π^{+}, photons, and neutrons produced in each of the five simulations and to evaluate the mean values n̄ and standard deviations 1σ. The average run time for simulations with the GC approach (≈15 min) is ≲500 times smaller than with the NL calculations (≈190 h). As shown in Table 3, the mean number of secondaries of each kind obtained in simulations with the GC and NL approaches agree within 1σ. The proton energy of 1 GeV is only just above the threshold for π^{−} production (cf. MacKinnon et al. 2020) so the number of these particles is very small. The long run times involved in the NL simulations limit the extent to which this can be investigated, but the NL and GC results are nonetheless in statistical agreement.
Statistics of secondary production by primary protons in a dipole magnetic field for simulations with the GC and NL approaches.
In Fig. 6 we show gammaray spectra produced by primary protons in a dipole magnetic field obtained in simulations with the GC approach. The spectra correspond to angleintegrated energy distributions of photons crossing the xyplane at z = 0 from the chromospheric–photospheric region to the coronal region. The primary protons are released from the top of the reference magnetic field line and are assumed to have an isotropic angular distribution and a powerlaw energy distribution. The panels show the total spectrum as well as the contributions from e^{+} annihilation, neutroncapture, proton inelastic processes, e^{±} bremsstrahlung, and π^{0} decay. As we can see, with the exception of nuclear deexcitation gammaray lines, whose contribution is not properly accounted for by the Geant4 model HadronPhysics QGSP_BERT_HP used for proton inelastic processes (Tang & Smith 2010), the simulated spectra exhibit all typical structures of gammaray spectra observed in solar flares: the e^{+}annihilation line at 511 keV, the neutroncapture line at 2.223 MeV, and the continuum emission components from π^{0} decay and bremsstrahlung of secondary e^{±} produced from π^{±} decay. The proton inelastic processes nevertheless produce a continuum emission component that dominates the spectra below ∼20 MeV, particularly for the simulation with spectral index δ = 3.
Fig. 6. Gammaray spectra from primary protons in a dipole magnetic field obtained in simulations with the GC approach. The primary protons released from the top of the reference magnetic field line are assumed to have an isotropic angular distribution and a powerlaw energy distribution with spectral indexes δ = 2 (left panel) and δ = 3 (right panel) in the range 1.0 MeV–1.0 GeV. The panels show the total spectrum, and the contributions from e^{+} annihilation, neutroncapture, proton inelastic processes, e^{±} bremsstrahlung, and π^{0} decay. 
Fig. 7. Measured photon fluxes as a function of energy from the flare of 12 June 2010 together with an approximate fit including the GEANT4 simulated spectrum of Fig. 6, a single powerlaw photo spectrum, and a powerlaw spectrum with exponential turnover at high energy. The parameters of the various components are shown in the bottom left corner. The fitting procedure is detailed in the text. The vertical dashed lines near 30 MeV show the precise boundaries of the energy ranges of GBM (< 30 MeV) and LAT (> 30 MeV). 
To show the usefulness of this kind of simulation for data interpretation we carried out a first illustrative comparison of the spectrum in Fig. 6 with the Fermi Gammaray Burst Monitor (GBM) and Large Area Telescope (LAT) data from the flare of 12 June 2010 (Ackermann et al. 2012) at time 00:55:40–00:58:50 UT. These data were previously compared with results of the FLUKA simulations by Tusnski et al. (2019) and MacKinnon et al. (2020). A rough fit to the data from both instruments was carried out in OSPEX and is shown in Fig. 7, combining the GEANT4 template of Fig. 6 (right) with two powerlaw components, one also with exponential turnover. The amplitude of the GEANT4 component and the parameters of the power law plus exponential are first fixed by fitting to the LAT data between 30–300 MeV. The parameters of the second powerlaw component are then fixed to give a best fit to the GBM data, which cover the energy range 300 keV–30 MeV. The GEANT4 spectrum alone gives a very good fit to the LAT data above 50 MeV (χ^{2} ≃ 1); we recall that MacKinnon et al. (2020) find that a fairly wide range of ion distributions can give spectra consistent with these data. The fit to the GBM data is much poorer, however, for several reasons. First, assuming a primary proton δ = 3 gives too few protons in the 1–100 MeV energy range to account for the general level of the deexcitation lines and the flux in the 2.223 and 0.511 MeV lines. We recall that the ratio of the fluxes at 4–7 MeV to those at ∼100 MeV has been used as a diagnostic of primary ion spectral index (Murphy et al. 1987). Evidently the data for this event require a value for δ significantly greater than 3. A fuller analysis would also have to include the contributions of αparticles and heavier accelerated nuclei (Tusnski et al. 2019) whose neglect here further decreases the line yield relative to the ∼100 MeV flux. However, the lack of a detailed deexcitation line spectrum in the GEANT4 results further worsens the fit, and the OSPEX fitting in fact attributes most of the flux at 0.3–30 MeV to the powerlaw components.
Nonetheless, this preliminary comparison with the data shows how useful such spectra can be. Future work will carry out more complete data fitting, and will compare the details of these spectra with those obtained assuming unmagnetised plasma (e.g., MacKinnon et al. 2020).
Fig. 8. Geometry of the doubledipole magnetic field model. 
Fig. 9. Trajectories of 1000 primary protons in a doubledipole magnetic field obtained in simulations with the GC approach. The primary protons are released from two different extended sources of ellipsoidal form: one centred at x = 0, y = 0 and z = 2000 km (a) and the other centred at x = 0, y = 0 and z = 5000 km (b). In both cases the primary protons are assumed to have an isotropic angular distribution and a powerlaw energy distribution of spectral index δ = 2 in the range 0.2–1.0 GeV. Panels c and d show the maps of secondary photons crossing the xyplane at z = 0 that are produced in the chromosferic–photospheric region by 3 × 10^{6} primary protons released from the sources depicted in panels a and b, respectively. 
3.4. Doubledipole magnetic field
A more representative and canonical magnetic flare configuration involves two interacting loops producing a magnetic reconnection that releases energy to accelerate particles (e.g., Najita & Orrall 1970; Aschwanden et al. 1999). To illustrate what is possible with the much faster GC simulations we show examples of secondary photon production in this magnetic field. We simulate this configuration using two dipoles and by releasing protons at midloop height in the middle of the two loops (Fig. 8). The dipoles have magnetic momenta parallel to the xaxis, , and are placed at and . The two reference field lines are defined by the parameters B_{f} = 1000 G, d_{x} = 500 km, d_{z} = 0.012 R_{s}, and R_{f} = 0.010 Rs.
We work these simulations in the GC approach alone, and for the simulations we set up detectors to obtain the trajectories of 1000 primary protons released from two different extended sources of ellipsoidal form: one centred at x = 0, y = 0 and z = 2000 km with principal axes of halflength hx = hy = 290 km and hz = 100 km, and another one centred at x = 0, y = 0 and z = 5000 km with principal axes of halflength hx = hy = 345 km and hz = 100 km. In both cases the primary protons are assumed to have an isotropic angular distribution and a powerlaw energy distribution of spectral index δ = 2 in the range from 0.2 to 1.0 GeV. As shown in the top panels of Fig. 9, the primary protons released from the sources located at heights z = 2000 km and z = 5000 km move along magnetic field lines that are directed to regions near the inner and the outer footpoints, respectively. In both cases, depending on the initial pitch angle, the primary protons can either precipitate into the chromospheric–photospheric region or become trapped in the coronal region bouncing back and forth between the footpoints. As in the simulations in a single dipole (Sect. 3.2), particles trapped in the coronal region are not tracked further when the total track length reaches ∼10^{5} km.
We also set up detectors of secondary photons crossing the xyplane at z = 0 which are produced in the chromospheric–photospheric region by 3 × 10^{6} primary protons released from the same two extended sources considered above. The results are shown in the bottom panels of Fig. 9. Most of the photons are detected around the footpoints of the magnetic loops; however, in panel d we see an intermediate site of photon accumulation. A great number of primary protons is needed to provide statistically reliable conclusions. Nevertheless we note that this kind of simulation, with highly complex magnetic topology would need an unreasonable run time in the NL framework, a large fraction of which would be expended following particles in the coronal region where they barely produce any observable radiation. We do not show examples here, but these simulations can also yield energy and angular distributions of escaping neutrons, for example.
4. Summary and conclusions
We have implemented a GC description of charged particle transport in magnetic fields within Geant4. The GC equations of motion are substituted for the exact NL equations, and are solved using essentially the same numerical methods. Single particles modelled using the GC equations of motion follow trajectories that are the same, to within a gyroradius, as those calculated using the exact NL approach: they follow field lines, drifting slowly across them at the correct rate. In a denser medium they slow down at the correct rate and stop at the correct depth. Calculating energy and angle distributions of secondaries would be computationally infeasible in the NL case, but total numbers agree to within the statistics of the simulations. Since energy and angle distributions of secondaries are matters of nuclear physics, not transport, we expect them to be consistent since the total numbers of secondaries are in agreement.
Most importantly, the GC run times are two to five orders of magnitude shorter than for the NL computation, depending on the energies of the primaries (Tables 1 and 2). This comparison was carried out adopting the same values of the RK4 precision parameters and the cord parameter (which governs the maximum discrepancy between straight line Monte Carlo steps and the exact curvilinear trajectory) for GC and NL calculations. It seems likely that further experimentation with these parameters can speed the GC calculation up further without compromising accuracy. Since we expect a dependence on detailed magnetic field model, density structure, for example, we defer such finetuning for future work applying the GC code.
The GC formalism is immediately applicable to active regions. Hudson et al. (2009), for example, used it to study energetic ions in the potentialfield sourcesurface (PFSS) model of the global solar magnetic field, although without any calculation of secondaries such as Geant4 makes possible. In interplanetary space, with largescale lengths and weak fields, a NL treatment may be equally feasible. One development of the work here would be to test the value of B/∇B and switch between GC and NL descriptions as needed. Such a development might be used in modelling aimed at interpreting the quiet sun gammaradiation (cf. Li et al. 2016; Mazziotta et al. 2020). It would also allow us to model configurations including nulls, in the vicinity of which particle adiabatic motion will break down (e.g., Petkaki & MacKinnon 1997).
The potential of the GC approach is shown in Figs. 6 and 9. Figure 6 shows gammaray spectra calculated in the presence of a nonuniform magnetic field. In MacKinnon et al. (2020) we showed examples of such spectra and compared them with observations of a solar flare. Since these were calculated neglecting any effect of magnetic field the constraints found on source parameters still needed further interpretation. The significance of fits to data will be much clearer when these have been obtained in the presence of a model of magnetic field. The primary and secondary particles in Fig. 9 have been followed through a moderately complex magnetic field, and the spatial distributions of secondaries against the solar disc calculated in consequence. This sort of detailed modelling would be computationally infeasible using the exact NL treatment, but could be used for example to test models for the longlasting gammaradiation observed after several large flares (Ackermann et al. 2014) or to model the spatial and spectral form of radio and submillimetre radiation produced by secondary positrons and electrons (Tuneu et al. 2016). Future work will explore this potential and give detailed solar physics applications.
Acknowledgments
This work was supported by the São Paulo Science State Agency (FAPESP) [grant numbers 2013/241553, 2016/234284, 2017/132825, 2018/081801]; Ministry of Education Agency for the Promotion of Graduate Studies (CAPES) [grant number 88881.310386/201801]; National Council for Science and Development (CNPq) [grant 307722/20198]. Solar physics research in the University of Glasgow is supported by an STFC consolidated grant. We thank Chris Osborne [Astronomy& Astrophysics Group, University of Glasgow] for useful discussions on C++. David Hamilton and John Annand from the Nuclear Physics Group of the University of Glasgow gave helpful advice on Geant4.
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 734, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 745, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Albert, A., et al. 2014, ApJ, 787, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Agostinelli, S., Allison, J., Amako, K., et al. 2003, Nucl. Instrum. Methods Phys. Res. A, 506, 250 [Google Scholar]
 Allison, J., Amako, K., Apostolakis, J., et al. 2006, IEEE Trans. Nucl. Sci., 53, 270 [Google Scholar]
 Aschwanden, M. J., Kosugi, T., Hanaoka, Y., Nishio, M., & Melrose, D. B. 1999, ApJ, 526, 1026 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, M. J., Coursey, J. S., Zucker, M. A., & Chang, J. 2005, ESTAR, PSTAR and ASTAR: Computer Programs for Calculating StoppingPower and Range Tables for Electrons, Protons and Helium Ions [Google Scholar]
 Bertini, H. W. 1963, Phys. Rev., 131, 1801 [NASA ADS] [CrossRef] [Google Scholar]
 Böhlen, T. T., Cerutti, F., Chin, M. P. W., et al. 2014, Nucl. Data Sheets, 120, 211 [CrossRef] [Google Scholar]
 Desorgher, L., Flückiger, E. O., & Gurtner, M. 2006, in 36th COSPAR Scientific Assembly, 36, 2361 [NASA ADS] [Google Scholar]
 Emslie, A. G., Dennis, B. R., Shih, A. Y., et al. 2012, ApJ, 759, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrari, A., Sala, P. R., Fasso, A., & Ranft, J. 2005, FLUKA: A multiparticle transport code, Tech. Rep. CERN200510, CERN [Google Scholar]
 Folger, G., & Wellisch, J. 2003, ArXiv eprints [arXiv:nuclth/0306007] [Google Scholar]
 Guo, J., Zeitlin, C., WimmerSchweingruber, R. F., et al. 2018, AJ, 155, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Hudson, H. S., MacKinnon, A. L., De Rosa, M. L., & Frewen, S. F. N. 2009, ApJ, 698, L86 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, K., Boggs, S. E., Smith, D. M., et al. 2005, Nature, 434, 1098 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, J. D. 1999, Classical Electrodynamics, 3rd edn. (New York, NY: Wiley) [Google Scholar]
 Kotoku, J., Makishima, K., Matsumoto, Y., et al. 2007, PASJ, 59, 1161 [NASA ADS] [CrossRef] [Google Scholar]
 Li, J. P., Zhou, A. H., & Wang, X. D. 2016, Res. Astron. Astrophys., 16, 8 [Google Scholar]
 MacKinnon, A., Szpigel, S., Giménez de Castro, C. G., & Tuneu, J. 2020, Sol. Phys., 295, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Matthiä, D., Ehresmann, B., Lohf, H., et al. 2016, J. Space Weather Space Clim., 6, A13 [CrossRef] [EDP Sciences] [Google Scholar]
 Mazziotta, M. N., Luque, P. D. L. T., Di Venere, L., et al. 2020, Phys. Rev. D, 101, 083011 [NASA ADS] [CrossRef] [Google Scholar]
 Merloni, A., & Fabian, A. C. 2001, MNRAS, 328, 958 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, J. A., Cargill, P. J., Emslie, A. G., et al. 1997, J. Geophys. Res., 102, 14631 [NASA ADS] [CrossRef] [Google Scholar]
 Murphy, R. J., Dermer, C. D., & Ramaty, R. 1987, ApJS, 63, 721 [NASA ADS] [CrossRef] [Google Scholar]
 Murphy, R. J., Kozlovsky, B., Share, G. H., Hua, X.M., & Ligenfelter, R. E. 2007, ApJS, 168, 167 [CrossRef] [Google Scholar]
 Najita, K., & Orrall, F. Q. 1970, Sol. Phys., 15, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Northrop, T. G. 1963, Rev. Geophys., 1, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Ozturk, M. K. 2012, Am. J. Phys., 80, 420 [NASA ADS] [CrossRef] [Google Scholar]
 Parks, G. K. 2004, Physics of space plasmas : an introduction [Google Scholar]
 Petkaki, P., & MacKinnon, A. L. 1997, Sol. Phys., 172, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Tang, S., & Smith, D. M. 2010, ApJ, 721, 1174 [NASA ADS] [CrossRef] [Google Scholar]
 Trottet, G., Raulin, J.P., Mackinnon, A., et al. 2015, Sol. Phys., 290, 2809 [NASA ADS] [CrossRef] [Google Scholar]
 Tuneu, J., Szpigel, S., de Castro, G., & MacKinnon, A. 2016, Proc. Int. Astron. Union, 12, 120 [CrossRef] [Google Scholar]
 Tusnski, D. S., Szpigel, S., Giménez de Castro, C. G., MacKinnon, A. L., & Simões, P. J. A. 2019, Sol. Phys., 294, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Vilmer, N., MacKinnon, A. L., & Hurford, G. J. 2011, Space Sci. Rev., 159, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, D., & Kelsey, M. 2015, Nucl. Instrum. Methods Phys. Res. Sect. A, 804, 175 [Google Scholar]
All Tables
Run times for simulations of a primary proton in a uniform magnetic field with the GC and NL approaches.
Run times for simulations of a primary proton in a dipole magnetic field with the GC and NL approaches.
Statistics of secondary production by primary protons in a dipole magnetic field for simulations with the GC and NL approaches.
All Figures
Fig. 1. Geometry of the dipole magnetic field model. 

In the text 
Fig. 2. Trajectories of a primary proton in a uniform magnetic field obtained in simulations with the GC and NL approaches. Left panels: initial kinetic energy E_{k0} = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panels: initial pitch angle θ = 35° and initial kinetic energies E_{k0} = 0.5 GeV (A), E_{k0} = 1.0 GeV (B) and E_{k0} = 5.0 GeV (C). 

In the text 
Fig. 3. Comparison between the energy losses of a primary proton in a uniform magnetic field obtained in simulations with the GC and NL approaches. Left panel: initial kinetic energy E_{k0} = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panel: initial pitch angle θ = 35° and initial kinetic energies E_{k0} = 0.5 GeV (A), E_{k0} = 1.0 GeV (B) and E_{k0} = 5.0 GeV (C). 

In the text 
Fig. 4. Trajectories of a primary proton with initial kinetic energy E_{k0} = 1.0 GeV in a dipole magnetic field obtained in simulations with the GC and NL approaches. Left panels: initial pitch angle θ = 15°. Right panels: initial pitch angle θ = 35°. 

In the text 
Fig. 5. Comparison between the energy losses of a primary proton in a dipole magnetic field obtained in simulations with the GC and NL approaches. The plots are for an initial pitch angle θ = 15° and initial kinetic energies E_{k0} = 0.5 GeV (A), E_{k0} = 1.0 GeV (B), and E_{k0} = 5.0 GeV (C). 

In the text 
Fig. 6. Gammaray spectra from primary protons in a dipole magnetic field obtained in simulations with the GC approach. The primary protons released from the top of the reference magnetic field line are assumed to have an isotropic angular distribution and a powerlaw energy distribution with spectral indexes δ = 2 (left panel) and δ = 3 (right panel) in the range 1.0 MeV–1.0 GeV. The panels show the total spectrum, and the contributions from e^{+} annihilation, neutroncapture, proton inelastic processes, e^{±} bremsstrahlung, and π^{0} decay. 

In the text 
Fig. 7. Measured photon fluxes as a function of energy from the flare of 12 June 2010 together with an approximate fit including the GEANT4 simulated spectrum of Fig. 6, a single powerlaw photo spectrum, and a powerlaw spectrum with exponential turnover at high energy. The parameters of the various components are shown in the bottom left corner. The fitting procedure is detailed in the text. The vertical dashed lines near 30 MeV show the precise boundaries of the energy ranges of GBM (< 30 MeV) and LAT (> 30 MeV). 

In the text 
Fig. 8. Geometry of the doubledipole magnetic field model. 

In the text 
Fig. 9. Trajectories of 1000 primary protons in a doubledipole magnetic field obtained in simulations with the GC approach. The primary protons are released from two different extended sources of ellipsoidal form: one centred at x = 0, y = 0 and z = 2000 km (a) and the other centred at x = 0, y = 0 and z = 5000 km (b). In both cases the primary protons are assumed to have an isotropic angular distribution and a powerlaw energy distribution of spectral index δ = 2 in the range 0.2–1.0 GeV. Panels c and d show the maps of secondary photons crossing the xyplane at z = 0 that are produced in the chromosferic–photospheric region by 3 × 10^{6} primary protons released from the sources depicted in panels a and b, respectively. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.