EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 654 (October 2021) A&A, 654 (2021) A82 Full HTML
Free Access
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/0004-6361/202141183
Published online 15 October 2021

© 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 X-ray, and gamma-ray observations combined with modelling of transport and secondary production. In the case of energetic ions nuclear reactions give rise to gamma-ray lines in the 0.5–8 MeV range through nuclear de-excitation 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, well-validated 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 gamma-ray 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 gamma-radiation 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 space-dependence 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 rL = γ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 ∼104 cm. This is to be compared with a typical active region lengthscale of 109 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 gyro-motion 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 rL ≪ B/|∇B|, where B is the magnitude of the magnetic field and ∇B is its gradient. With this condition the magnetic moment μ= p ⊥ 2 /(2mB) $ \mu=p_\perp^2/(2mB) $, 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)

d ( γ m v ) d t = q v × B ( r ) , $$ \begin{aligned} \frac{\mathrm{d}(\gamma m \mathbf v )}{\mathrm{d}t}= q\mathbf v \times \mathbf B(r) , \end{aligned} $$(1)

where v and B(r) are the particle velocity and magnetic field vectors, respectively. As usual γ = (1−v2/c2)−1/2. In a homogeneous magnetic field B(r) = B0, and the motion of the particle is well-known: rectilinear motion parallel to B0 and oscillatory circular motion with radius rL.

We write B ( r ) = B ( r ) b ̂ $ \mathbf{B(r)} = B(\mathbf{r})\hat{\mathbf{b}} $ 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

v = v ‖ b ̂ + v ⊥ , $$ \begin{aligned} \mathbf v = { v}_\parallel \hat{\mathbf{b }} + \mathbf v _\perp , \end{aligned} $$(2)

where v ⊥ . b ̂ = 0 $ \mathbf{v}_\perp.\hat{\mathbf{b}} = 0 $. The GC approach (Northrop 1963; Ozturk 2012) applies as long as

r L ≪ B | ∇ B | . $$ \begin{aligned} r_L \ll \frac{B}{|\nabla B|} . \end{aligned} $$(3)

Then the magnetic moment μ= p ⊥ 2 /(2mB) $ \mu=p_\perp^2/(2mB) $ 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

r = R + ρ , $$ \begin{aligned} \mathbf r = \mathbf R + \mathbf {\rho} , \end{aligned} $$(4)

where R is the GC position vector (centre of the gyro-orbit) and ρ is the gyration radius vector, lying in the plane perpendicular to R and having magnitude rL. The GC description concentrates on the motion of R and gives the position to within a gyroradius.

Averaging Eq. (1) over one gyro-period we obtain the GC equations of motion, which are often written as follows:

d R d t = γ m v 2 2 q B 2 ( 1 + v ‖ 2 v 2 ) b ̂ × ∇ B ( R ) + v ‖ b ̂ , $$ \frac{\mathrm{d}\mathbf R }{\mathrm{d}t}=\frac{\gamma m{ v}^2}{2qB^2}\left(1+\frac{{ v}_{\parallel }^2}{{ v}^2}\right)\mathbf {\hat{b}} \times \mathbf {\nabla} B(\mathbf R ) +{ v}_{\parallel }\mathbf {\hat{b}} , $$(5)

d v ‖ d t = − μ γ 2 m b ̂ · ∇ B ( R ) . $$ \frac{\mathrm{d}{v}_{\parallel }}{\mathrm{d}t}=-\frac{\mu }{\gamma ^2 m}\mathbf {\hat{b}} \cdot \mathbf {\nabla} B(\mathbf R ). $$(6)

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 b ̂ $ \mathbf{\hat{b}} $ 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:

d r d t = v , $$ \begin{aligned}&\frac{\mathrm{d} \mathbf r }{\mathrm{d} t} \, = \, \mathbf v ,\end{aligned} $$(7)

d p d t = F . $$ \begin{aligned}&\frac{\mathrm{d} \mathbf p }{\mathrm{d} t} \, = \, \mathbf F . \end{aligned} $$(8)

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):

d r d t = γ m v 2 2 q B 2 ( 1 + ( v . b ̂ ) 2 v 2 ) b ̂ × ∇ B ( r ) + ( v . b ̂ ) b ̂ . $$ \begin{aligned} \frac{\mathrm{d} \mathbf r }{\mathrm{d} t}=\frac{\gamma m{ v}^2}{2qB^2}\left(1+\frac{\left(\mathbf v.\hat{b} \right)^2}{{ v}^2}\right)\mathbf {\hat{b}} \times \mathbf {\nabla} B(\mathbf r ) + (\mathbf v.\hat{b} )\mathbf {\hat{b}} \, . \\ \end{aligned} $$(9)

To evaluate the force F we introduce two unit vectors ξ ̂ $ \mathbf{\hat{\xi}} $ and ζ ̂ $ \mathbf{\hat{\zeta}} $ in the plane perpendicular to B. These satisfy the following,

ξ ̂ . B = ζ ̂ . B = ξ ̂ . ζ ̂ = 0 , $$ \begin{aligned} \mathbf {\hat{\xi }} .\mathbf B \, = \, \mathbf {\hat{\zeta }} .\mathbf B \, = \, \mathbf {\hat{\xi }} .\mathbf {\hat{\zeta }} \, = \, 0 , \end{aligned} $$(10)

and we discuss their explicit construction below. Then we can write

p = p ‖ b ̂ + p ⊥ ( δ ξ ̂ + ϵ ζ ̂ ) , $$ \begin{aligned} \mathbf p \, = \, p_\parallel \hat{\mathbf{b }} + p_\perp \left(\delta \mathbf {\hat{\xi }} + \epsilon \mathbf {\hat{\zeta }} \right) , \end{aligned} $$(11)

and so:

d p d t = d p ‖ d t b ̂ + p ‖ d b ̂ d t + d p ⊥ d t ( δ ξ ̂ + ϵ ζ ̂ ) + p ⊥ d d t ( δ ξ ̂ + ϵ ζ ̂ ) . $$ \begin{aligned} \frac{\mathrm{d} \mathbf p }{\mathrm{d} t} = \frac{\mathrm{d} p_\parallel }{\mathrm{d} t}\hat{\mathbf{b }} + p_\parallel \frac{\mathrm{d} \hat{\mathbf{b }}}{\mathrm{d} t} + \frac{\mathrm{d} p_\perp }{\mathrm{d} t}\left(\delta \mathbf {\hat{\xi }} + \epsilon \mathbf {\hat{\zeta }} \right)+ p_\perp \frac{\mathrm{d} }{\mathrm{d} t}\left(\delta \mathbf {\hat{\xi }} + \epsilon \mathbf {\hat{\zeta }} \right) . \end{aligned} $$(12)

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 b ̂ $ \hat{\mathbf{b}} $, ξ ̂ $ \hat{\mathbf{\xi}} $, and ζ ̂ $ \hat{\mathbf{\zeta}} $ 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, dp2/dt = 0, and so

d p ‖ d t = − p ⊥ p ‖ d p ⊥ d t . $$ \begin{aligned} \frac{\mathrm{d} p_\parallel }{\mathrm{d} t} = -\frac{p_\perp }{p_\parallel } \frac{\mathrm{d} p_\perp }{\mathrm{d} t} . \end{aligned} $$(13)

From the constancy of magnetic moment we have

d p ⊥ d t = p ⊥ 2 B ( v GC . ∇ ) B , $$ \begin{aligned} \frac{\mathrm{d} p_\perp }{\mathrm{d} t} \, = \, \frac{p_\perp }{2B}\left(\mathbf v _{\rm GC}.\nabla \right) B , \end{aligned} $$(14)

since B only changes for the particle because it is moving. We write vGC here for the velocity of the guiding centre, given by Eq. (5). In the same way,

d b ̂ d t = ( v GC . ∇ ) b ̂ . $$ \begin{aligned} \frac{\mathrm{d} \hat{\mathbf{b }}}{\mathrm{d} t} \, = \, \left(\mathbf v _{\rm GC}.\nabla \right)\hat{\mathbf{b }} . \end{aligned} $$(15)

We still need to specify how to calculate the derivative d ( δ ξ ̂ + ϵ ζ ̂ ) / d t $ \mathrm{d}\left(\delta \mathbf{\hat{\xi}} + \epsilon\mathbf{\hat{\zeta}}\right)/\mathrm{d}t $, 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 x ̂ $ \mathbf{ \hat{x}} $, y ̂ $ \mathbf{ \hat{y}} $, and z ̂ $ \mathbf{ \hat{z}} $. Then we can construct ξ ̂ $ \mathbf{\hat{\xi}} $ and ζ ̂ $ \mathbf{\hat{\zeta}} $ as

ξ ̂ = 1 1 − b z 2 b ̂ × z ̂ = 1 1 − b z 2 ( b y − b x 0 ) , $$ \begin{aligned}&\mathbf {\hat{\xi }} = \frac{1}{\sqrt{1-b_z^2}}\mathbf {\hat{b}} \times \mathbf {\hat{z}} =\frac{1}{\sqrt{1-b_z^2}} \begin{pmatrix} b_{ y}\\ -b_x\\ 0 \end{pmatrix}, \end{aligned} $$(16)

ζ ̂ = ξ ̂ × b ̂ = 1 1 − b z 2 [ z ̂ − ( b ̂ . z ̂ ) b ̂ ] = 1 1 − b z 2 ( − b x b z − b y b z 1 − b z 2 ) . $$ \begin{aligned}&\mathbf {\hat{\zeta }} =\mathbf {\hat{\xi }} \times \mathbf {\hat{b}} = \frac{1}{\sqrt{1-b_z^2}}\left[\mathbf {\hat{z}} -\left(\mathbf {\hat{b}} .\mathbf {\hat{z}} \right)\;\mathbf {\hat{b}} \right] = \frac{1}{\sqrt{1-b_z^2}} \begin{pmatrix} -b_x b_z\\ -b_{ y} b_z\\ 1-b_z^2 \end{pmatrix} . \end{aligned} $$(17)

Equations (16) and (17) fail if at any position B becomes parallel to z ̂ $ \mathbf{ \hat{z}} $ to within working precision. In this case we substitute x ̂ $ \mathbf{ \hat{x}} $ or y ̂ $ \mathbf{ \hat{y}} $ for z ̂ $ \mathbf{ \hat{z}} $, amounting to no more than a rotation of the already random gyrophase. Then finally we have

d d t ( δ ξ ̂ + ϵ ζ ̂ ) = δ ( v GC . ∇ ) ξ ̂ + ϵ ( v GC . ∇ ) ζ ̂ , $$ \begin{aligned} \frac{\mathrm{d} }{\mathrm{d} t}\left(\delta \mathbf {\hat{\xi }} + \epsilon \mathbf {\hat{\zeta }} \right) \, = \, \delta \left(\mathbf v _{\rm GC}.\nabla \right) \mathbf {\hat{\xi }} + \epsilon \left(\mathbf v _{\rm GC}.\nabla \right) \mathbf {\hat{\zeta }} , \end{aligned} $$(18)

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 Runge-Kutta (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 ∂Bx/∂x, ∂Bx/∂y,...∂Bz/∂y, ∂Bz/∂z at the vertices of the grid and have them available for calculation of the gradients of b ̂ $ \mathbf{\hat{b}} $, ξ ̂ $ \hat{\xi} $, and ζ ̂ $ \hat{\zeta} $. These are approximated by simple differencing:

d B x d x ≃ B x ( x + d , y , z ) − B x ( x − d , y , z ) 2 d . $$ \begin{aligned} \frac{\mathrm{d} B_x}{\mathrm{d} x} \, \simeq \, \frac{B_x(x+d,{ y},z)-B_x(x-d,{ y},z)}{2d} . \end{aligned} $$(19)

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,

d d s = 1 v d d t , $$ \begin{aligned} \frac{\mathrm{d} }{\mathrm{d} s} \, = \, \frac{1}{{ v}}\frac{\mathrm{d} }{\mathrm{d} t} , \end{aligned} $$(20)

so Eq. (5) becomes

d r d s = ( v ̂ . b ) b + γ m v 2 q B 2 ( 1 + v ‖ 2 v 2 ) b × ∇ B ( r ) . $$ \begin{aligned} \frac{\mathrm{d} \mathbf r }{\mathrm{d} s} \, = \, \left(\mathbf {\hat{v}} .\mathbf b \right) \mathbf b + \frac{\gamma m{ v}}{2qB^2}\left(1+\frac{{ v}_{\parallel }^2}{{ v}^2}\right)\mathbf b \times \mathbf {\nabla} B(\mathbf r ). \end{aligned} $$(21)

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 cross-field drift term will normally be negligible in any one step so we can neglect it to note

| d r d s | = | ( v ̂ . b ̂ ) | . $$ \begin{aligned} \left|\frac{\mathrm{d} \mathbf r }{\mathrm{d} s} \right|\, = \, \left|\left(\mathbf {\hat{v}} .\mathbf {\hat{b}} \right) \right|. \end{aligned} $$(22)

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 slowing-down 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 Quark-gluon 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 plane-parallel 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 × 104 km and faces perpendicular to the coordinate axes. The z-coordinate 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 4He, 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 gamma-ray line production (Murphy et al. 2007; Trottet et al. 2015). A xy-plane at z = 0 divides the box into two half-spaces. The half-space at z > 0 represents the coronal region for which we assume a density ρ = 10−13  g cm−3. The half-space 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 gamma-ray 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 power-law energy distribution over a range Emin − Emax:

d n ( E ) d E = N 0 E − δ H ( E − E min ) H ( E max − E ) . $$ \begin{aligned} \frac{\mathrm{d}n(E)}{\mathrm{d}E}=N_0 \; E^{-\delta } H(E-E_{\min }) H(E_{\max }-E) . \end{aligned} $$(23)

Here E is the kinetic energy of the primary proton, δ is the spectral index, H is the Heaviside step function, and N0 is a normalisation constant defined such that

∫ E min E max d n ( E ) d E d E = 1 . $$ \begin{aligned} \int ^{E_{\max }} _{E_{\min }} \frac{\mathrm{d}n(E)}{\mathrm{d}E}\mathrm{d}E=1 . \end{aligned} $$(24)

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 low-energy threshold for transport. Geant4 includes robust and well validated physics-based 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 z-axis

B ( r ) = B 0 z ̂ , $$ \begin{aligned} \mathbf B (\mathbf r )=B_0 \; \mathbf {\hat{z}} , \end{aligned} $$(25)

where B0 is the (constant) magnitude of the magnetic field at any point r = x x ̂ + y y ̂ + z z ̂ $ \mathbf{r}=x \mathbf{\hat{x}} + \mathit{y} \mathbf{\hat{y}} + z \mathbf{\hat{z}} $ of the model space.

The general expression for the magnetic field produced by a magnetic dipole is given by

B ( r ) = 3 [ μ · ( r − r d ) ] ( r − r d ) − | r − r d | 2 μ | r − r d | 5 , $$ \begin{aligned} \mathbf B (\mathbf r )=\frac{3[\mathbf {\mu} \cdot (\mathbf r -\mathbf r_d ) ](\mathbf r -\mathbf r_d )-|\mathbf r -\mathbf r_d |^2 \mathbf {\mu} }{|\mathbf r -\mathbf r_d |^5} , \end{aligned} $$(26)

where rd 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 x-axis, μ = μ x ̂ $ \mathbf{\mu}=\mu \mathbf{\hat{x}} $, placed at a point r d = − d z z ̂ , d z > 0 $ \mathbf{r_d}=-d_z \mathbf{\hat{z}}, d_z > 0 $. The Cartesian components of the magnetic field vector B(r) are then given by

B x = μ 2 x 2 − y 2 − ( z + d z ) 2 [ x 2 + y 2 + ( z + d z ) 2 ] 5 / 2 , B y = μ 3 x y [ x 2 + y 2 + ( z + d z ) 2 ] 5 / 2 , B z = μ 3 x ( z + d z ) [ x 2 + y 2 + ( z + d z ) 2 ] 5 / 2 , $$ \begin{aligned}&B_x=\mu \;\frac{2x^2-{ y}^2-(z+d_z)^2}{[x^2+{ y}^2+(z+d_z)^2]^{5/2}}, \nonumber \\&B_{ y}=\mu \;\frac{3x{ y}}{[x^2+{ y}^2+(z+d_z)^2]^{5/2}}, \nonumber \\&B_z=\mu \;\frac{3x(z+d_z)}{[x^2+{ y}^2+(z+d_z)^2]^{5/2}} , \end{aligned} $$(27)

thumbnail 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 xz-plane and is assumed to cross the xy-plane at the footpoints, as shown in Fig. 1. Hence,

μ = B f ( R f 2 + d z 2 ) 5 / 2 4 R f 4 + 5 R f 2 d z 2 + d z 4 , $$ \begin{aligned} \mu =B_f \;\frac{ (R_f^2+d_z^2)^{5/2}}{\sqrt{4R_f^4+5R_f^2 d_z^2 + d_z^4}} , \end{aligned} $$(28)

where Bf is the magnitude of the magnetic field at the footpoints and Rf 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 dz, Rf, and Bf.

3.2. Propagation and slowing down in the guiding centre limit

In order to verify the equivalence between the slowing-down 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™ i7-7500 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 z-axis with magnitude B0 = 500 G, for which there is no cross-field 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 Ek0 = 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: Ek0 = 0.5 GeV, Ek0 = 1.0 GeV, and Ek0 = 5.0 GeV.

In Fig. 2 we compare the primary proton trajectories projected in the xz-plane 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 stopping-depth obtained for the primary proton with Ek0 = 0.5 GeV is larger than that obtained for the primary proton with Ek0 = 0.25 GeV and the same initial pitch angle. For primary protons with Ek0 = 1.0 GeV and Ek0 = 5.0 GeV the stopping depths decrease, being even smaller for Ek0 = 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 x-axis is orders of magnitude smaller than the scale set for the z-axis, such that for the NL trajectories the change of the gyroradius rL 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.

thumbnail 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 Ek0 = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panels: initial pitch angle θ = 35° and initial kinetic energies Ek0 = 0.5 GeV (A), Ek0 = 1.0 GeV (B) and Ek0 = 5.0 GeV (C).

In Fig. 3 we show the plots for the height z as a function of the ratio Ek/Ek0 (where Ek 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 Ek0 = 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 slowing-down 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 Ek0 = 0.5 GeV (curve labelled A in the figure) loses its energy continuously, coming to rest at a stopping-depth that is larger than that obtained for a primary proton with Ek0 = 0.25 GeV and the same initial pitch angle (curves B and C). Primary protons with Ek0 = 1.0 GeV and Ek0 = 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.

thumbnail 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 Ek0 = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panel: initial pitch angle θ = 35° and initial kinetic energies Ek0 = 0.5 GeV (A), Ek0 = 1.0 GeV (B) and Ek0 = 5.0 GeV (C).

thumbnail Fig. 4.

Trajectories of a primary proton with initial kinetic energy Ek0 = 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: Bf = 1000 G, dz = 0.012 Rs, and Rf = 0.010 Rs (where Rs = 6.96 × 105 km is the solar radius). The top of the reference magnetic field line is located at a height h = 10, 623 km above the xy-plane at z = 0 and the magnitude of the magnetic field at this point is Bh = 121.5 G. The loss-cone angle (the minimum pitch angle for a particle to be trapped in the low-density coronal region) is then given by sin2(θc) = Bh/Bf, 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 Ek0 = 1.0 GeV and two different values of initial pitch angle, θ = 15° and θ = 35°, which lie inside and outside the loss-cone, 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 cross-field 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 Ek/Ek0 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: Ek0 = 0.5 GeV, Ek0 = 1.0 GeV and Ek0 = 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 Ek0 = 0.5 GeV (curve A) continuously loses all its energy by ionisation collisions until it comes to rest, while the primary protons with Ek0 = 1.0 GeV (curve B) and Ek0 = 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 slowing-down rates are nearly the same for simulations with the GC and NL approaches.

thumbnail 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 Ek0 = 0.5 GeV (A), Ek0 = 1.0 GeV (B), and Ek0 = 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 Ek0 = 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 105 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.

Table 1.

Run times for simulations of a primary proton in a uniform magnetic field with the GC and NL approaches.

Table 2.

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 H8DGT-HIBQF composed of two 64 bits CPUs AMD Magny-Cours 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.

Table 3.

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 gamma-ray spectra produced by primary protons in a dipole magnetic field obtained in simulations with the GC approach. The spectra correspond to angle-integrated energy distributions of photons crossing the xy-plane 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 power-law energy distribution. The panels show the total spectrum as well as the contributions from e+ annihilation, neutron-capture, proton inelastic processes, e± bremsstrahlung, and π0 decay. As we can see, with the exception of nuclear de-excitation gamma-ray 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 gamma-ray spectra observed in solar flares: the e+-annihilation line at 511 keV, the neutron-capture 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.

thumbnail Fig. 6.

Gamma-ray 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 power-law 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, neutron-capture, proton inelastic processes, e± bremsstrahlung, and π0 decay.

thumbnail 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 power-law photo spectrum, and a power-law 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 Gamma-ray 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 power-law 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 power-law 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 de-excitation 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 de-excitation 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 power-law 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).

thumbnail Fig. 8.

Geometry of the double-dipole magnetic field model.

thumbnail Fig. 9.

Trajectories of 1000 primary protons in a double-dipole 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 power-law 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 xy-plane at z = 0 that are produced in the chromosferic–photospheric region by 3 × 106 primary protons released from the sources depicted in panels a and b, respectively.

3.4. Double-dipole 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 mid-loop height in the middle of the two loops (Fig. 8). The dipoles have magnetic momenta parallel to the x-axis, μ = μ x ̂ $ \mathbf{\mu}=\mu \mathbf{\hat{x}} $, and are placed at r d 1 = − ( R f + d x ) x ̂ − d z z ̂ $ \mathbf{r_{d_1}}=-(R_f+d_x)\mathbf{\hat{x}}-d_z \mathbf{\hat{z}} $ and r d 2 = ( R f + d x ) x ̂ − d z z ̂ $ \mathbf{r_{d_2}}=(R_f+d_x)\mathbf{\hat{x}}-d_z\mathbf{\hat{z}} $. The two reference field lines are defined by the parameters Bf = 1000 G, dx = 500 km, dz = 0.012 Rs, and Rf = 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 half-length 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 half-length hx = hy = 345 km and hz = 100 km. In both cases the primary protons are assumed to have an isotropic angular distribution and a power-law 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 ∼105 km.

We also set up detectors of secondary photons crossing the xy-plane at z = 0 which are produced in the chromospheric–photospheric region by 3 × 106 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 fine-tuning 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 potential-field source-surface (PFSS) model of the global solar magnetic field, although without any calculation of secondaries such as Geant4 makes possible. In interplanetary space, with large-scale 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 gamma-radiation (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 gamma-ray spectra calculated in the presence of a non-uniform 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 long-lasting gamma-radiation observed after several large flares (Ackermann et al. 2014) or to model the spatial and spectral form of radio and sub-millimetre 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/24155-3, 2016/23428-4, 2017/13282-5, 2018/08180-1]; Ministry of Education Agency for the Promotion of Graduate Studies (CAPES) [grant number 88881.310386/2018-01]; National Council for Science and Development (CNPq) [grant 307722/2019-8]. 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

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 734, 116 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 745, 144 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ackermann, M., Ajello, M., Albert, A., et al. 2014, ApJ, 787, 15 [NASA ADS] [CrossRef] [Google Scholar]
  4. Agostinelli, S., Allison, J., Amako, K., et al. 2003, Nucl. Instrum. Methods Phys. Res. A, 506, 250 [Google Scholar]
  5. Allison, J., Amako, K., Apostolakis, J., et al. 2006, IEEE Trans. Nucl. Sci., 53, 270 [Google Scholar]
  6. Aschwanden, M. J., Kosugi, T., Hanaoka, Y., Nishio, M., & Melrose, D. B. 1999, ApJ, 526, 1026 [NASA ADS] [CrossRef] [Google Scholar]
  7. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  8. Berger, M. J., Coursey, J. S., Zucker, M. A., & Chang, J. 2005, ESTAR, PSTAR and ASTAR: Computer Programs for Calculating Stopping-Power and Range Tables for Electrons, Protons and Helium Ions [Google Scholar]
  9. Bertini, H. W. 1963, Phys. Rev., 131, 1801 [NASA ADS] [CrossRef] [Google Scholar]
  10. Böhlen, T. T., Cerutti, F., Chin, M. P. W., et al. 2014, Nucl. Data Sheets, 120, 211 [CrossRef] [Google Scholar]
  11. Desorgher, L., Flückiger, E. O., & Gurtner, M. 2006, in 36th COSPAR Scientific Assembly, 36, 2361 [NASA ADS] [Google Scholar]
  12. Emslie, A. G., Dennis, B. R., Shih, A. Y., et al. 2012, ApJ, 759, 71 [NASA ADS] [CrossRef] [Google Scholar]
  13. Ferrari, A., Sala, P. R., Fasso, A., & Ranft, J. 2005, FLUKA: A multi-particle transport code, Tech. Rep. CERN-2005-10, CERN [Google Scholar]
  14. Folger, G., & Wellisch, J. 2003, ArXiv e-prints [arXiv:nucl-th/0306007] [Google Scholar]
  15. Guo, J., Zeitlin, C., Wimmer-Schweingruber, R. F., et al. 2018, AJ, 155, 49 [NASA ADS] [CrossRef] [Google Scholar]
  16. Hudson, H. S., MacKinnon, A. L., De Rosa, M. L., & Frewen, S. F. N. 2009, ApJ, 698, L86 [NASA ADS] [CrossRef] [Google Scholar]
  17. Hurley, K., Boggs, S. E., Smith, D. M., et al. 2005, Nature, 434, 1098 [NASA ADS] [CrossRef] [Google Scholar]
  18. Jackson, J. D. 1999, Classical Electrodynamics, 3rd edn. (New York, NY: Wiley) [Google Scholar]
  19. Kotoku, J., Makishima, K., Matsumoto, Y., et al. 2007, PASJ, 59, 1161 [NASA ADS] [CrossRef] [Google Scholar]
  20. Li, J. P., Zhou, A. H., & Wang, X. D. 2016, Res. Astron. Astrophys., 16, 8 [Google Scholar]
  21. MacKinnon, A., Szpigel, S., Giménez de Castro, C. G., & Tuneu, J. 2020, Sol. Phys., 295, 174 [NASA ADS] [CrossRef] [Google Scholar]
  22. Matthiä, D., Ehresmann, B., Lohf, H., et al. 2016, J. Space Weather Space Clim., 6, A13 [CrossRef] [EDP Sciences] [Google Scholar]
  23. Mazziotta, M. N., Luque, P. D. L. T., Di Venere, L., et al. 2020, Phys. Rev. D, 101, 083011 [NASA ADS] [CrossRef] [Google Scholar]
  24. Merloni, A., & Fabian, A. C. 2001, MNRAS, 328, 958 [NASA ADS] [CrossRef] [Google Scholar]
  25. Miller, J. A., Cargill, P. J., Emslie, A. G., et al. 1997, J. Geophys. Res., 102, 14631 [NASA ADS] [CrossRef] [Google Scholar]
  26. Murphy, R. J., Dermer, C. D., & Ramaty, R. 1987, ApJS, 63, 721 [NASA ADS] [CrossRef] [Google Scholar]
  27. Murphy, R. J., Kozlovsky, B., Share, G. H., Hua, X.-M., & Ligenfelter, R. E. 2007, ApJS, 168, 167 [CrossRef] [Google Scholar]
  28. Najita, K., & Orrall, F. Q. 1970, Sol. Phys., 15, 176 [NASA ADS] [CrossRef] [Google Scholar]
  29. Northrop, T. G. 1963, Rev. Geophys., 1, 283 [NASA ADS] [CrossRef] [Google Scholar]
  30. Ozturk, M. K. 2012, Am. J. Phys., 80, 420 [NASA ADS] [CrossRef] [Google Scholar]
  31. Parks, G. K. 2004, Physics of space plasmas : an introduction [Google Scholar]
  32. Petkaki, P., & MacKinnon, A. L. 1997, Sol. Phys., 172, 279 [NASA ADS] [CrossRef] [Google Scholar]
  33. Tang, S., & Smith, D. M. 2010, ApJ, 721, 1174 [NASA ADS] [CrossRef] [Google Scholar]
  34. Trottet, G., Raulin, J.-P., Mackinnon, A., et al. 2015, Sol. Phys., 290, 2809 [NASA ADS] [CrossRef] [Google Scholar]
  35. Tuneu, J., Szpigel, S., de Castro, G., & MacKinnon, A. 2016, Proc. Int. Astron. Union, 12, 120 [CrossRef] [Google Scholar]
  36. 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]
  37. Vilmer, N., MacKinnon, A. L., & Hurford, G. J. 2011, Space Sci. Rev., 159, 167 [NASA ADS] [CrossRef] [Google Scholar]
  38. Wright, D., & Kelsey, M. 2015, Nucl. Instrum. Methods Phys. Res. Sect. A, 804, 175 [Google Scholar]

All Tables

Table 1.

Run times for simulations of a primary proton in a uniform magnetic field with the GC and NL approaches.

In the text
Table 2.

Run times for simulations of a primary proton in a dipole magnetic field with the GC and NL approaches.

In the text
Table 3.

Statistics of secondary production by primary protons in a dipole magnetic field for simulations with the GC and NL approaches.

In the text

All Figures

thumbnail Fig. 1.

Geometry of the dipole magnetic field model.
In the text
thumbnail 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 Ek0 = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panels: initial pitch angle θ = 35° and initial kinetic energies Ek0 = 0.5 GeV (A), Ek0 = 1.0 GeV (B) and Ek0 = 5.0 GeV (C).
In the text
thumbnail 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 Ek0 = 0.25 GeV and initial pitch angles θ = 15° (A), θ = 35° (B) and θ = 80° (C). Right panel: initial pitch angle θ = 35° and initial kinetic energies Ek0 = 0.5 GeV (A), Ek0 = 1.0 GeV (B) and Ek0 = 5.0 GeV (C).
In the text
thumbnail Fig. 4.

Trajectories of a primary proton with initial kinetic energy Ek0 = 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
thumbnail 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 Ek0 = 0.5 GeV (A), Ek0 = 1.0 GeV (B), and Ek0 = 5.0 GeV (C).
In the text
thumbnail Fig. 6.

Gamma-ray 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 power-law 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, neutron-capture, proton inelastic processes, e± bremsstrahlung, and π0 decay.
In the text
thumbnail 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 power-law photo spectrum, and a power-law 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
thumbnail Fig. 8.

Geometry of the double-dipole magnetic field model.
In the text
thumbnail Fig. 9.

Trajectories of 1000 primary protons in a double-dipole 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 power-law 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 xy-plane at z = 0 that are produced in the chromosferic–photospheric region by 3 × 106 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.