3D propagation of relativistic solar protons through interplanetary space

Context. Solar Energetic Particles (SEPs) with energy in the GeV range can propagate to Earth from their acceleration region near the Sun and produce Ground Level Enhancements (GLEs). The traditional approach to interpreting and modelling GLE observations assumes particle propagation only parallel to the magnetic field lines of interplanetary space, i.e. it is spatially 1D. Recent measurements by PAMELA have characterised SEP properties at 1 AU for the ~100 MeV-1 GeV range at high spectral resolution. Aims. We model the transport of GLE-energy solar protons through the Interplanetary Magnetic Field (IMF) using a 3D approach, to assess the effect of the Heliospheric Current Sheet (HCS) and drifts associated to the gradient and curvature of the Parker spiral. The latter are influenced by the IMF polarity. We derive 1 AU observables and compare the simulation results with data from PAMELA. Methods. We use a 3D test particle model including a HCS. Monoenergetic populations are studied first to obtain a qualitative picture of propagation patterns and numbers of crossings of the 1 AU sphere. Simulations for power law injection are used to derive intensity profiles and fluence spectra at 1 AU. A simulation for a specific event, GLE 71, is used to compare with PAMELA data. Results. Spatial patterns of 1 AU crossings and the average number of crossings are strongly influenced by 3D effects, with significant differences between periods of A+ and A- polarities. The decay time constant of 1 AU intensity profiles varies depending on the polarity and position of the observer, and it is not a simple function of the mean free path as in 1D models. Energy dependent leakage from the injection flux tube is particularly important for GLE energy particles, in many cases resulting in a roll-over in the fluence spectrum.


Introduction
Ions of relativistic energies can be accelerated at or near the Sun during flare and Coronal Mass Ejection (CME) events. When detected in the interplanetary medium, for example near Earth, they constitute the high energy portion of the spectrum of Solar Energetic Particles (SEPs) (Mewaldt et al. 2012;Cohen & Mewaldt 2018), whose properties are an important tracer of the acceleration processes and of the propagation through the Interplanetary Magnetic Field (IMF) .
Relativistic solar ions may produce secondary particles when they interact with Earth's atmosphere, causing so-called Ground Level Enhancements (GLEs), detected by ground-based instruments such as neutron monitors (Poluianov et al. 2017;McCracken et al. 2012;Mishev et al. 2018). Protons in the energy range ∼0.5-30 GeV are thought to be the main contributors to GLEs (eg. McCracken et al. (2012)). The exact range is different for each neutron monitor, as the minimum energy depends on the location via the geomagnetic cutoff and the maximum on the instrument sensitivity and spectrum hardness for a specific event. GLEs are much less frequent than SEP events detected by spacecraft instrumentation, which is typically sensitive to protons up to ∼100 MeV. Only 72 GLE events have been de-tected by the worldwide network of neutron monitors from 1942 to the present time (eg. Belov et al. (2010); Nitta et al. (2012); Gopalswamy et al. (2012)).
Recent SEP observations from the PAMELA detector have allowed us to fill the particle energy gap between traditional spacecraft instrumentation and neutron monitors, and routinely detect relativistic solar protons in the range from ∼100 MeV to a few GeV (Adriani et al. 2015;Bruno et al. 2018). The new observations call for modelling tools that describe the acceleration and propagation of particles at these energies. In addition, simulations of transport through the IMF are necessary to relate the detections of high energy SEPs at 1 AU to the numbers of interacting particles at the Sun which produce solar γ-ray events detected by FERMI (de Nolfo et al. 2019;Share et al. 2018;Klein et al. 2018).
A number of studies have modelled the propagation of relativistic solar protons through the IMF using spatially 1D descriptions, to interpret neutron monitor observations. The effect of magnetic field turbulence on particle propagation is typically described as pitch-angle scattering, characterised by a mean free path λ. Bieber et al. (2004) and Sáiz et al. (2005) used a model based on the focused transport equation to fit data for two GLE events. Strauss et al. (2017) used a focused transport model to calculate rise and decay times of GLEs. Heber et al. (2018) combined 1D propagation within interplanetary space of GLEenergy particles with trajectory integration through magnetospheric configurations. Li & Lee (2019) found analytical expressions for the flux profile and anisotropy of relativistic protons using a focused transport approach within specific scattering conditions, and used them to fit the 2005 January 20 GLE. The 1D approximation, which assumes that particles remain tied to the magnetic field line on which they were injected, is therefore the standard approach used to model the interplanetary propagation of solar relativistic protons, and to analyse GLE observations (e.g. Nitta et al. (2012)). Within this approximation, the effects of IMF polarity and of the heliospheric current sheet on the propagation of relativistic protons are neglected.
However when energetic protons above ∼1 GeV are modelled within the heliosphere in the context of Galactic Cosmic Ray (GCR) modulation studies (e.g. Potgieter & Vos (2017)), a spatially 3D description of particle propagation is thought to be necessary, due to effects such as IMF gradient and curvature drifts, diffusion in the direction perpendicular to the average field, and the influence of the Heliospheric Current Sheet (HCS) (e.g. Parker (1965); Kota & Jokipii (1983); Burger (2012)).
It is the aim of this paper to model the interplanetary propagation of relativistic protons by means of a fully 3D approach, allowing us to discuss the effects of the HCS and IMF polarity on 1 AU observables. Our earlier work has pointed out that drifts due to the gradient and curvature of the Parker spiral IMF do affect the propagation of SEPs, with their importance increasing with particle energy and being particularly significant for heavy ions Dalla et al. 2013Dalla et al. , 2017a. Analysis of the role played by a flat HCS (Battarbee et al. 2017) and by a wavy HCS (Battarbee et al. 2018a) on SEPs injected with power law distributions in the range 10-800 MeV demonstrated the role of injection region location and IMF polarity, and elucidated how the HCS provides an efficient means for particle transport in longitude.
In this paper, we focus on relativistic protons in the energy range from a few hundred MeV to 10 GeV and demonstrate the need for an approach that describes propagation as fully 3D, unlike the traditional approaches to GLE modelling. In particular we show that once a 3D approach is adopted and a HCS is introduced in the model, significant dependencies of 1 AU observables on the magnetic polarity of the IMF are observed. We point out how the latter affects time-intensity profiles and spectra, analysed at multiple locations defined with respect to the magnetic flux tube with nominal connection to the centre of the injection region. We also focus on a specific relativistic particle event for which PAMELA detected protons over a wide energy range, GLE 71, occurring on May 17, 2012, and compare our modelled observables with the data (Adriani et al. 2015;Bruno et al. 2018). This is the first comparison of SEP PAMELA observations with a model.
In Section 2 we present our model and the results of simple monoenergetic injection simulations, including a discussion of the number of 1 AU crossings. In Section 3 we consider a powerlaw distribution of relativistic protons and discuss how transport through interplanetary space shapes the 1 AU observables over a grid of locations. In Section 4 a comparison between our model and PAMELA intensity profiles is presented for GLE 71. We discuss our results in Section 5.

Simulations of monoenergetic populations
We model relativistic proton propagation through space by integrating particle trajectories in 3D via a full orbit test particle code Dalla & Browning 2005). Particle acceleration is not modelled and injection characteristics of the accelerated population are specified as input. The IMF is characterised by two polarities separated by a model wavy HCS (Battarbee et al. 2018a). Using standard terminology from GCR studies, the configuration with magnetic field pointing outwards in the northern hemisphere and inwards in the south is referred to as A + and that with opposite polarity as A − .
Scattering due to turbulence in the interplanetary magnetic field is modelled by generating a sequence of scattering events for each particle, compatible with a specified value of the mean free path λ, as described in Marsh et al. (2013). At each scattering event the direction of the particle's velocity is reassigned randomly from a uniform spherical distribution. There is no consensus within the literature about the degree of scattering experienced by GLE energy protons in their travel to 1 AU. In the simulations of Bieber et al. (2004) and Sáiz et al. (2005), fitting to GLE data, within their 1D model, yielded λ ∼ 0.1 AU. Li & Lee (2019) were able to reproduce observations only by assuming different scattering conditions for different phases of a GLE event: at the beginning of the event they used λ=4 AU, meaning near scatter-free conditions, while later in the event strong scattering, with λ an order of magnitude smaller, was required to fit the data. In our simulations we consider a variety of mean free paths, kept constant over time and location in the heliosphere, and we neglect the dependence of λ on energy for the relativistic particle energy range we consider. In this initial study we do not introduce any perpendicular scattering, therefore any transport across the magnetic field seen in the simulations is due to drift and HCS effects.
It is instructive to analyse the propagation of monoenergetic populations of relativistic protons, to visualise how 1 AU observables vary with particle energy. Each monoenergetic population consists of N=10,000 particles, injected instantaneously from a small region at the Sun of angular extent 8×8 • , located at r=2 R sun . While in actual SEP events the source region may in fact be much more extended, the key properties of the propagation are revealed more clearly if the injection is localised within the model.
The magnetic field in the simulation is given by a Parker spiral field. We use the method described by Battarbee et al. (2018a) to include a HCS. When the presence of a HCS is taken into account, parameters of the HCS such as the tilt angle α nl , the polarity of the IMF and the position of the injection region with respect to the HCS, have a strong influence on the particle propagation (Battarbee et al. 2018a) . Figure 1 shows longitude-latitude maps of crossings of the 1 AU sphere summed over the entire duration of the simulations, for monoenergetic populations at 500 MeV, 1 GeV and 10 GeV, where these populations were followed up to a time t f =61 hr. The mean free path λ is 0.1 AU. The injection region, corresponding to the dark red pixels, e.g. in the top left plot, is located above the HCS, at longitude 76 • and latitude 11 • . The tilt of the HCS is α nl =37 • . The maps show 1 AU crossings in a heliographic coordinate system that is corotating with the Sun.

Maps of 1 AU crossings
The left panels show maps for an A + configuration and the right panels for A − , so that in the former case particles tend 180 150 120 90  60  30 0  30  60  90 120 150  Longitude   60  50  40  30  20  10  0  10  20  30  40  50   Latitude   500 MeV, A +   180 150 120 90  60  30 0  30  60  90 120 150  Longitude   60  50  40  30  20  10  0  10  20  30  40  The 8×8 • injection region at the Sun is located at longitude 76 • and latitude 11 • , above the HCS, and the zero of the coordinate system on the 1 AU map has been shifted so that the flux tube through the injection region appears at N11W76 on this map. The tilt of the HCS is α nl =37 • . All simulations used N=10,000 protons, solar wind speed v sw =400 km s −1 and mean free path λ=0.1 AU. Contour lines are plotted for the following values of the number of crossings: 1000 (green), 316 (blue), 100 (lilac), 31 (red) and 10 (black).
to move towards the HCS and in the latter away from it, due to gradient and curvature drift in the Parker spiral IMF Marsh et al. 2013;Battarbee et al. 2018a). This motion towards/away from the HCS follows standard GCR patterns (Jokipii et al. 1977). Since gradient and curvature drift effects increase with energy, the 10 GeV particles show the largest transport across the field, and for the latter population the peak counts location is southwards of the injection region for the A + configuration and northwards for A − . In addition to gradient and curvature drift, HCS drift also affects the spatial patterns in Figure  1. In the A + case, as they reach the HCS by drifting southwards, particles experience a strong westward HCS drift that spreads them efficiently in longitude. In the A − case, a drift along the HCS is also observed, though it is less pronounced compared to the A + situation, because gradient and curvature drift tend to move particles away from the HCS, and it is in the opposite direction (eastwards).
Looking at the bottom panels, for 10 GeV, one can see that although the injection region is only 8 • ×8 • in extent, the entire 1 AU sphere is accessible to particles, despite the fact that the injection was localised. It is interesting to note that at these energies, although rapid transport across the field allows particle access to regions far away from the injection region, it also quickly dilutes the population, making it more difficult for it to be detected above the GCR background. Looking at the two bottom rows, it is clear that over the energy range of interest for GLEs, interplanetary propagation is fully 3D.
The patterns seen in Figure 1 present some differences and similarities to the maps of 1 AU crossings presented by Battarbee et al. (2018a), where a power law proton population in the energy range 10-800 MeV was considered. Their maps were therefore dominated by ∼10 MeV particles, which experience much smaller drift compared to relativistic protons, resulting in a less pronounced drift along the HCS for starting locations that were not directly located on the HCS itself. The overall qualitative dependence of patterns on A + versus A − is the same as in Battarbee et al. (2018a).
The panels of Figure 1 do not include the effect of corotation, i.e. the fact that, in the spacecraft frame, magnetic flux tubes filled with particles cross a number of heliospheric longitudes over time. Corotation increases the spatial extent of the event in longitude (for an example of maps with and without the inclusion of corotation see Figure 1 of (Battarbee et al. 2018a)), however at the energies considered Figure 1 the effects of corotation are less evident than at lower energy.

Average number of 1 AU crossings
In addition to the spatial patterns of 1 AU crossings, it is interesting to consider N cross , the average number of crossings of the 1 AU sphere by SEP protons of a specified kinetic energy. Particles may cross 1 AU more than once as they scatter back and forth in their propagation, so that this parameter is a strong function of the mean free path (Chollet et al. 2010). N cross is needed to  estimate the total number of SEPs at 1 AU from spacecraft detections of fluxes (Mewaldt et al. 2008). Therefore knowledge of N cross , for example from transport simulations, allows one to compare 1 AU SEP numbers with the number of interacting particles at the Sun, deduced e.g. from γ-ray observations (de Nolfo et al. 2019). We estimate N cross (E) from our simulations, where E is particle energy, by taking the average of the number of 1 AU crossings over each monoenergetic population considered, with crossings collected over the entire 1 AU sphere. It should be noted that particles do decelerate as they propagate through interplanetary space (see e.g. Dalla et al. (2015)), however the effect is less prominent at the energies considered here, so that it is a reasonable assumption to take the initial energy as E. Table 1 displays N cross values for the λ=0.1 AU simulations displayed in Figure 1 (see also de Nolfo et al. (2019)). A strong dependence of N cross on the IMF polarity is therefore deduced from our simulations, with the number of crossings being much larger for A − polarity than for A + . This behaviour is equivalent to the polarity dependence of fluence that was discussed by Battarbee et al. (2018a). It should be noted that the distribution of N cross values is generally quite broad, so that the standard deviation for the values in Table 1 is almost as large as the values themselves. Figure 2 shows how the number of crossings evolves with time. Here the count rate I (number of crossings divided by accumulation time) is plotted, using a collection area equal to the whole 1 AU sphere, for injection energy E=1 GeV and for the two polarities. The right hand panel (λ=0.1 AU) corresponds to the same simulations displayed in the central row of Figure 1, while the left hand panel has λ=0.5 AU. There is a large difference in the time evolution of number of crossings depending on the polarity of the IMF, with the A + polarity decay being much faster than for A − .
The reason for the differences between A + and A − in Figure  2 and Table 1 is that in the former configuration, drift along the HCS is more prominent, so that protons move towards the outer heliosphere faster than for A − and a significantly lower number of 1 AU crossings occur. The reason why the two curves in Figure 2 are very similar at early times is that it takes a finite amount of time for particles to drift down to the HCS in the A + case, at which point HCS drifts set in. Our findings on the influence of IMF polarity on number of crossings is confirmed by a completely independent test particle simulation code with flat HCS (Chollet et al. 2010;de Nolfo et al. 2019). We note that changing the parameters of the HCS (for example the tilt angle) does not affect N cross strongly and that its energy dependence (fewer crossings at higher energies) is a result of the particle populations at high energies propagating faster towards the outer heliosphere.

Simulations of power-law populations
We consider a proton population injected with a distribution of energies that follows a power law, and propagate it through interplanetary space using the same HCS configuration as in Section 2. We choose a spectral index at injection γ=2 for the energy range 100 MeV-1 GeV. The population is injected from the same location as the monoenergetic runs shown in Figure 1 and with the same parameters. Therefore also in this analysis, we use a small 8×8 • injection region. The mean free path is λ=0.1 AU.

Intensity profiles
To produce intensity profiles, counts are collected over 10 • ×10 • portions of the 1 AU sphere that mimic a variety of observer locations with respect to the injection region. Here the observer is not corotating with the Sun but is in the so-called spacecraft frame, in which magnetic flux tubes rooted at the Sun rotate by 14.3 • per day. Figure 3 shows intensity profiles at a variety of locations for the energy ranges 100-400 MeV (blue) and 700-1000 MeV (green). The top grid refers to an A + IMF configuration and the bottom grid to A − . Observer locations are specified using labels [∆φ 1AU , ∆δ 1AU ], where ∆φ 1AU is the heliographic longitude and ∆δ 1AU the heliographic latitude of the observer relative to the Parker spiral field line through to the centre of the particle injection region. The panel labelled [0,0] (red label) corresponds to an observer connected to the centre of the injection region at the time of injection, and the other panels to less well connected observers (black labels). In a 1D model, intensities would be zero everywhere apart from the well connected panel, [0,0].
Moving from left to right along a row in Figure 3 one can see count rate profiles for observers at the same latitude and progressively more western longitudes (i.e. source region becoming more eastern). Here one can see the important effect of corotation, in the lower energy range, resulting in a less sharp rise phase and later time of peak intensity as the source region becomes more eastern, as already noted in our previous studies Dalla et al. 2017a;Laitinen et al. 2018). Different rows correspond to different observer latitudes, becoming more southern as one moves downwards. The observer locations for A + (A − ) have been chosen to reflect the fact that, as shown in the maps of Figure 1, the spatial extent is mostly downwards (upwards). Figure 3 shows that the event duration is much shorter in the 700-1000 MeV range compared to the 100-400 MeV range. This is due to the combination of two effects: 1) the higher energy protons travel away from the inner heliosphere faster and 2) they experience stronger transport across the field due to drift effects (as shown in Figure 1), resulting in much faster dilution of the population. Therefore more efficient drift across the field does not necessarily mean a higher probability of detection at far away locations, since dilution works against detection above background at a given spacecraft. At lower energies, particles are confined inside a 'cloud' around the injection flux tube and as a result of corotation they can produce significant count rates over extended times.
Comparing the top and bottom sets of grids in Figure 3, two main differences between A + and A − are observed: the overall spatial extent of the event is larger for the A − case, in agreement with the monoenergetic 1 AU maps shown in Figure 1, and at many observers the decay phase tends to last longer in the A − configuration compared to A + , replicating the behaviour seen for the global crossings in Figure 2.
The slope of the decay phase varies significantly for different locations for a given polarity configuration, as well as between A + and A − . Thus in 3D this parameter is not simply a reflection of the value of the mean free path λ, as would be the case in a 1D model, but it is the result of a number of processes that include IMF polarity and HCS effects and dilution due to transport across the magnetic field.
Considering the peak intensities versus observer location in Figure 3, one can see that in the higher energy range they decay faster with ∆φ 1AU than for the lower one. Therefore a gaussian fit to the peak intensities versus ∆φ 1AU would yield a smaller width σ of the particle distribution at the higher energies.

Fluence spectra
The fluence spectra for the same locations and configurations as in Figure 3 are presented in Figure 4, where counts were accumulated over the entire duration of the simulation. Although the injection spectrum is a power law with γ=2, it is evident that a variety of spectral shapes are seen at the different observers, as a result of 3D propagation effects.
The fact that drifts effects are stronger for high energies has an influence on particle spectra: as a result of the dilution effect discussed in Section 2.1 at the best connected location the spectrum is no longer a power-law but displays a roll-over. At locations away from the well connected one a variety of features are observed: in most locations the spectra display roll-over features, while at some, for example for the [20,0] and [30,0] panels in the A + simulation, the spectrum displays a hardening at high energies.
In addition to the simulation for γ=2 (as shown in Figure 4) we analysed spectra also for a case with initial injection spectrum with γ=3 and found that the qualitative behaviour in this case is very similar to that for the case γ=2, apart from spectra being generally softer.
Overall the stronger drift-associated dilution at high energies compared to lower energies is a key influence on the spectra in our 3D simulations. In addition, at the lower end of the range shown in Figure 4, adiabatic deceleration is another energydependent process affecting the observed spectral shapes . Our simulations show that any mechanism that produces energy-dependent escape from the flux tube will 'process' the injection spectrum.
Within SEP observations, rollover features are observed frequently (Bruno et al. 2018). A high-energy rollover is predicted by diffusive shock acceleration theory, accounting for particles escaping the shock region during acceleration (Ellison & Ramaty 1985;Lee & Ryan 1986;Lee 2005). The results shown in Figure 4 show that 3D transport effects are also expected to play a role contributing to the spectral features observed at 1 AU.

Comparison with PAMELA data for GLE 71
In addition to performing idealised runs, as described in Sections 2 and 3, we applied our 3D relativistic proton simulations to a specific SEP event for which preliminary data were made available to us by the PAMELA collaboration. The event is GLE 71, occurring on May 17th, 2012. Here we present the results of initial simulations in which we considered protons in the energy range 80-1300 MeV, injected instantaneously with a power law spectrum with γ=2.8 from a region at 2 solar radii, centered on 180 150 120 90  60  30 0  30  60  90 120 150  Longitude   60  50  40  30  20  10  0  10  20  30  40  50 Latitude GLE71 Fig. 5. Maps of cumulative 1 AU crossings in a heliocentric coordinate system corotating with the Sun, for protons 80-1300 MeV, for the GLE 71 event. The longitude axis has been shifted so that the center of the injection region is at longitude 76 • and latitude 11 • in this 1 AU plot. The mean free path is λ=0.3 AU. the flare location, which was N11W76. The final time of the simulation was 61 hours. A solar wind speed v sw =400 km s −1 was used. For most GLE events, detailed information on the extent of the injection of relativistic protons near the Sun is not available, but there are indications that it is much smaller compared to lower energy protons (Gopalswamy et al. 2012), i.e. that only the nose of the shock is an efficient accelerator at the high energies. We chose a relativistic proton injection region of 40×40 • (with the size being the same at all energies considered here) and assumed a constant acceleration efficiency within it. The number of particles in the simulations was N=3,000,000.
GLE 71 was studied in detail in an earlier publication which focussed on comparing simulations with multi-spacecraft SEP data for energies up to ∼100 MeV (Battarbee et al. 2018b). In that work, the tilt angle best fitting the conditions during the event was found to be α nl =57 • , within an A − IMF configuration (see Figure 3 of Battarbee et al. (2018b)). We used the same HCS parameters in the simulations described here. The HCS for the GLE 71 event is more 'wavy' than the one seen in Figure 1.
We carried out simulations for two values of the mean free path, λ=1.0 AU and 0.3 AU, assumed to be independent of energy. Figure 5 shows a map of 1 AU crossings for the simulation with λ=0.3 AU. Because the extended injection region is wider than in the simulations presented in Section 2.1 and it intersects the HCS, a strong HCS drift (eastward because of the A − IMF configuration) is observed.
The source region for GLE 71 was magnetically well connected to Earth, so that an Earth observer was located at a position [∆φ 1AU , ∆δ 1AU ]=[-13 • , -13 • ] with respect to the centre of the injection region at the time of the flare., i.e. within the 40×40 • injection region. Therefore drift along the HCS did not play a main role in determining SEP arrival in this event. To derive simulated intensity profiles near Earth we collected counts within a 10×10 • collection tile to ensure good statistics. Figure 6 shows the intensity profiles in four energy channels for simulations with mean free path λ=1.0 AU (left) and 0.3 AU (middle), compared with the PAMELA intensity profiles (right). In the simulation with λ=1 AU, an increase in the slope of the decay in the lowest energy channel is seen after ∼30 hours, at the time when connection with the magnetic flux tubes spanning the injection region is lost due to corotation. The same effect is present but to a lesser degree in the λ=0.3 AU simulation. At higher energies, drift effects smooth the intensity profile and there is no sharp change in the decay time constant. When the mean free path is shorter, the peak intensity is reached later.
The simulation with λ=0.3 AU appears to produce a better fit in the lower energy channels, although at the higher energies the simulated decay is faster than in the PAMELA data. This may  indicate that the size of the injection region at higher energies is smaller than for lower energies, or may be related to an energy dependence of the mean free path. Neither of these effects are included in the initial simulations presented here. It should be noted that although the values of mean free path are very different, the slope of decay is not dissimilar in the two simulations of Figure 6 . This is unlike the behaviour in 1D simulations, in which the slope is a strong function of the mean free path.

Discussion and conclusions
We presented 3D simulations of relativistic proton propagation from the Sun to 1 AU, which included 3D effects associated with particle drift and the presence of a HCS. We considered monoenergetic and power-law populations, injected from a small region at the Sun, to study the qualitative aspects of 3D propagation. In addition we performed initial simulations for an extended injection region aiming to reproduce PAMELA observations for the GLE 71 event. We plan to extend the modelling to other PAMELA events in future. In further work, modelling of SEP propagation through the IMF will be combined with tracing of trajectories in the magnetospheric magnetic field, allowing comparison with GLE observations. It should be stressed that our simulations have focussed on the role of IMF polarity and HCS, while other potentially important processes such as perpendicular scattering and magnetic field line meandering (Laitinen et al. 2016) have not been included. The injection of relativistic protons has been assumed to be instantaneous and located near the Sun, which is a reasonable approximation at these energies (Gopalswamy et al. 2012).
The main conclusions from our analysis are as follows: -Propagation of relativistic protons is strongly influenced by the IMF polarity (via gradient and curvature drifts) and the HCS (via HCS drift), making a 3D description necessary. Relativistic protons are not confined to the magnetic flux tube in which they were injected. They experience dilution within the interplanetary medium much faster than ∼10 MeV protons, making their detection above background at a location initially not connected to the injection region, less likely than at lower energies. Corotation is less important at relativistic energies compared to lower proton energies. In contrast to the 1D approach, leakage from the magnetic flux tube in which the particles were injected is a key new phenomenon within a 3D description. -There are significant differences in the relativistic proton propagation for A + and A − configurations, a fact not captured by 1D models, which do not include IMF polarity and a HCS. The average number of 1 AU crossings is significantly larger for A − than for A + , due to the latter configuration causing efficient outward HCS drift that helps particles leave the inner heliosphere faster. Maps of 1 AU crossings show that A + configurations are characterised by stronger HCS-induced propagation across heliolongitudes compared to A − . -In 3D, injection properties of the SEP population are processed by transport, making intensity profiles and spectra strongly dependent on observer location. Injection spectra are only weakly reflected in fluence spectra at 1 AU, which present roll-overs at many observer locations, despite a power-law injection. The decay constant of intensity profiles is not related to the mean free path in a simple way, as is the case in 1D models, and it is influenced strongly by the IMF polarity. -Comparison of our simulation results with PAMELA observations in the energy range 80 MeV -1 GeV for GLE 71 (May 17th 2012) shows that resonable agreement with data can be obtained by choosing an injection region of 40×40 • , a mean free path λ=0.3 AU and injection spectrum with γ=2.8. Varying any of these parameters, as well as modifying assumptions on the energy dependence of λ and injection region size will influence the intensity profiles. Such an analysis will be the subject of future study.