Open Access
Issue
A&A
Volume 639, July 2020
Article Number A105
Number of page(s) 9
Section The Sun and the Heliosphere
DOI https://doi.org/10.1051/0004-6361/201937338
Published online 20 July 2020

© S. Dalla et al. 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. 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), which have been observed in ground-based neutron monitor data (Belov et al. 2010; Nitta et al. 2012; Gopalswamy et al. 2012; McCracken et al. 2012; Mishev et al. 2018). Protons in the energy range of ∼0.5–30 GeV are thought to be the main contributors to GLEs (eg. McCracken et al. 2012). GLEs are much less frequent than SEP events that are detected by spacecraft instrumentation, which is typically sensitive to protons up to ∼100 MeV. Only 72 GLE events have been detected by the worldwide network of neutron monitors since 1942 (eg. Belov et al. 2010).

Recent SEP observations from the Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA) detectors have allowed us to fill the particle energy gap between traditional spacecraft instrumentation and neutron monitors, as well as 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 propagation 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 the Fermi Gamma-ray Space Telescope (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, which is 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 GLE-energy 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 they 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.

A well developed theory of the propagation of galactic cosmic rays (GCRs), relativistic protons originating outside the heliosphere and propagating through the IMF, has been used to describe GCR modulation over several decades (e.g. Potgieter & Vos 2017). Within GCR models, typically dealing with protons of energies above ∼1 GeV, 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 (Marsh et al. 2013; Dalla et al. 2013, 2017a,b). 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 preliminary data from its detectors (Adriani et al. 2015; Bruno et al. 2018). This is the first comparison of SEP PAMELA data with a model.

In Sect. 2 we present our model and the results of simple monoenergetic injection simulations, including a discussion of the number of 1 AU crossings. In Sect. 3 we consider a power-law distribution of relativistic protons and discuss how transport through interplanetary space shapes the 1 AU observables over a grid of locations. In Sect. 4 a comparison between our model and PAMELA intensity profiles is presented for GLE 71. We discuss our results in Sect. 5.

2. Simulations of monoenergetic populations

We model relativistic proton propagation through space by integrating particle trajectories in 3D via a full orbit test particle code (Marsh et al. 2013; 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 simulated by means of the so-called “ad-hoc scattering” method. A sequence of Poisson-distributed scattering events for each particle is generated, compatible with a mean scattering time tscat = λ/v, where is λ the specified value of the mean free path and v the particle’s speed. At each scattering event, the direction of the particle’s velocity is reassigned randomly from a uniform spherical distribution (Kelly et al. 2012; Marsh et al. 2013). This method for describing scattering within SEP test particle simulations has been used by a number of groups over the years (e.g. Kocharov et al. 1998; Pei et al. 2006; Chollet et al. 2010; Kelly et al. 2012; Marsh et al. 2013). The value of the scattering mean free path within the ad-hoc scattering method is equivalent to that of traditional diffusion descriptions. Kocharov et al. (1998) directly compared the ad-hoc scattering approach (termed small time-step isotropisation (SSI) model in their work) with a traditional diffusion-convection description: they obtained very close agreement in SEP time intensity profiles at 1 AU when the same value is used for λ in the ad-hoc test particle approach and as parallel mean free path in the diffusion-convection model (see Fig. 4 of Kocharov et al. 1998). They also compared the results of the ad-hoc scattering description, in which the pitch-angle may change by a large angle during a scattering event, with two small angle scattering descriptions within the focussed transport equation, one isotropic and one anisotropic (indicated in their work as IAS and AAS respectively). They found that for the same value of the mean free path, 1 AU time intensity profiles for all these models are very similar, with some differences in the peak intensities and closely matching decay phases and durations (see Fig. 5 of Kocharov et al. 1998).

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 we neglect the dependence of λ on energy for the relativistic particle energy range we consider.

In this initial study we do not explicitly introduce a term describing perpendicular transport associated with turbulence in the solar wind magnetic field, for example due to magnetic field line meandering (Laitinen et al. 2016). Our scattering description does implicitly result in minor random-walk of the particle’s gyrocentre across the magnetic field, of the order of a Larmor radius, at each scattering event. This finite Larmor radius effect is small and it is negligible compared to typical cross-field diffusion due to random-walk, or meandering, of turbulent magnetic field lines (e.g. Jokipii 1966; Giacalone & Jokipii 1999; Fraschetti & Jokipii 2011). Thus turbulence-associated perpendicular transport is not included in our simulations and motion across the magnetic field seen in our results is mainly 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 Rsun. 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).

2.1. Maps of 1 AU crossings

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 tf = 61 hr. The mean free path λ is 0.1 AU. The injection region, corresponding to the dark red pixels, for example 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.

thumbnail Fig. 1.

Maps of cumulative 1 AU crossings in a heliocentric coordinate system corotating with the Sun, for monoenergetic SEP proton populations, with energy as indicated in each panel. Left column: A+ configuration of the IMF; right column: A configuration. 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 vsw = 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).

The left panels show maps for an A+ configuration and the right panels for A, so that in the former case particles tend to move towards the HCS and in the latter away from it, due to gradient and curvature drift in the Parker spiral IMF (Dalla et al. 2013; Marsh et al. 2013; Battarbee et al. 2018a). This motion towards or 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 Fig. 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 Fig. 1 present some differences and similarities to the maps of 1 AU crossings presented by Battarbee et al. (2018a): in the latter study, 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 Fig. 1 do not include the effect of corotation, that is 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 Fig. 1 of Battarbee et al. 2018a), however at the energies considered Fig. 1 the effects of corotation are less evident than at lower energy.

2.2. Average number of 1 AU crossings per particle

In addition to the spatial patterns of crossings of the 1 AU sphere, it is interesting to consider , the average number of 1 AU crossings per particle, for a specified SEP 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). is needed to estimate the total number of SEPs at 1 AU from spacecraft detections of fluxes (Mewaldt et al. 2008). Therefore knowledge of , for example from transport simulations, allows one to compare 1 AU SEP numbers with the number of interacting particles at the Sun, deduced for example from γ-ray observations (de Nolfo et al. 2019).

We derive from our model, where E is particle energy, by obtaining the number of 1 AU crossings per particle for each integrated trajectory in our monoenergetic population simulation, with crossings collected over the entire 1 AU sphere and for the duration of the simulation, and calculating its average over the population. 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 values for the λ = 0.1 AU simulations displayed in Fig. 1 and, for comparison, a case with λ = 0.5 AU (see also de Nolfo et al. 2019). A strong dependence of 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 the number of crossings per particle is generally quite broad, so that the standard deviation for the averages in Table 1 is almost as large as the values themselves.

Table 1.

Average number of 1 AU crossings per particle, , as a function of SEP proton kinetic energy E, for A+ and A configurations, and mean free path λ.

Figure 2 shows the time evolution of the count rate I (number of 1 AU detections divided by accumulation time), using the whole 1 AU sphere as collection area, 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 Fig. 1, while the left hand panel has λ = 0.5 AU. There is a large difference in the time evolution of the count rate depending on the polarity of the IMF, with the A+ polarity decay being much faster than for A.

thumbnail Fig. 2.

1 AU crossing count rates versus time summed over all heliographic longitudes and latitudes, for a monoenergetic proton population of initial kinetic energy E = 1 GeV and mean free path λ = 0.5 AU (left) and λ = 0.1 AU (right), for A+ and A configurations of the IMF. Other parameters of the simulations are as in Fig. 1.

The reason for the differences between A+ and A in Fig. 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 Fig. 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 per particle 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 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.

3. 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 Sect. 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 Fig. 1 and with the same parameters. Therefore also in this analysis, we use a small 8 × 8° injection region.

3.1. 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. Fig. 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 [Δϕ1 AU, Δδ1 AU], where Δϕ1 AU is the heliographic longitude and Δδ1 AU 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].

thumbnail Fig. 3.

Proton count rates versus time for A+ (top) and A (bottom) configurations of the IMF, at a variety of 1 AU locations with respect to the best connected location ([0, 0]), for a power law population, for the proton energy ranges 100–400 MeV (blue) and 700–1000 MeV (green).

Moving from left to right along a row in Fig. 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 (Marsh et al. 2015; 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 Fig. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 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.

3.2. Fluence spectra

The fluence spectra for the same locations and configurations as in Fig. 3 are presented in Fig. 4. 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.

thumbnail Fig. 4.

Fluence energy spectra for A+ (top) and A (bottom) configurations of the IMF, at a variety of 1 AU locations with respect to the best connected location ([0, 0]), for a power law population. The solid lines in the [0, 0] panels give the slope of the injection spectrum. Parameters of the simulations are as in Fig. 3.

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 Sect. 2.1 at the best connected location the spectrum is no longer a power-law but displays a roll-over. Rollover features are observed in PAMELA spectra (Bruno et al. 2018). At locations away from the well connected ones a variety of features are observed, connected to dilution at high energies and the fact that lower energy particles drift across the field less efficiently. In addition, at the lower end of the energy range shown in Fig. 4 adiabatic deceleration affects the spectra (Dalla et al. 2015). Our simulation shows that any mechanism that produces energy-dependent escape from the flux tube “processes” the injection spectrum.

In addition to the simulation for γ = 2 (as shown in Fig. 4), we also analysed spectra for a population with initial injection spectrum with γ = 3. In the latter case we found a qualitative behaviour very similar to that for the case γ = 2, apart from spectra being generally softer.

4. Comparison with PAMELA data for GLE 71

In addition to performing idealised runs, as described in Sects. 2 and 3, we applied our 3D relativistic proton simulations to a specific SEP event for which PAMELA 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, centred on the flare location, which was N11W76. The final time of the simulation was 61 hours. A solar wind speed vsw = 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), so that only the nose of the shock is an efficient accelerator at the high energies. Therefore we chose an injection region of 40 × 40° (with the size being the same at all energies considered here) and assume 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 Fig. 3 of Battarbee et al. 2018b). The HCS for this event is more “wavy” than the one seen in Fig. 1.

We carried out simulations for two values of the mean free path, assumed to be independent of energy, λ = 1.0 AU and 0.3 AU. 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 Sect. 2.1 and it intersects the HCS, a strong HCS drift (eastward because of the A IMF configuration) is observed.

thumbnail 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 centre of the injection region is at N11W76 in this plot. The mean free path is λ = 0.3 AU.

The source region for GLE 71 was magnetically well connected to Earth, so that an Earth observer was located at a position [Δϕ1 AU, Δδ1 AU] = [−2°, −13°] with respect to the centre of the injection region at the time of the flare, that is within the 40 × 40° injection region. Therefore drift along the HCS did not play a 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). The PAMELA data shown are based on omni-directional measurements taken in low Earth orbit. However, they account for a correction factor related to the pitch angle anisotropy registered during the first few hours, in particular the first polar pass (see Adriani et al. 2015). The PAMELA intensity time profiles are broadly consistent with those measured by the Energetic Proton, Electron, and Alpha Detector (EPEAD) and High Energy Proton and Alpha Detector (HEPAD) instruments on board the Geostationary Operational Environmental Satellites (GOES) 13 and 15.

thumbnail Fig. 6.

Time-intensity profiles of protons in four energy ranges, for the GLE 71 event. The left and middle panels show simulation results for an observer at the location of Earth, for λ = 1.0 AU and 0.3 AU respectively. The right panel shows preliminary data from PAMELA.

It is useful to comment on the qualitative differences between the two simulations and on the comparison with the data. Regarding the initial phase of the event, it is noticeable that the peak intensity is reached very quickly in the λ = 1 AU simulation, while in the λ = 0.3 AU case peak intensities are reached after a longer time, in a way that matches the observations better. Following the peak, intensities decay rapidly for λ = 1 AU, another feature that is not present in the data, suggesting a situation with stronger scattering, better reproduced by the simulation with λ = 0.3 AU. Following the initial fast decay, the slope of the intensity time profile increases in the simulations, to values closer to those of the observations.

Starting around t = 30 hrs, in the λ = 1 AU modelled intensities, an increase in the slope of the decay is seen in the lowest energy channels, which is not present in the PAMELA data. The same effect is visible, though to a lesser degree, in the λ = 0.3 AU simulation. This behaviour results from loss of magnetic connection of the observer to the flux tubes within which acceleration took place (the injection region in our simulation), as a result of corotation. In the higher energy channels the profiles are smoother and the sharp change in the decay time constant t=30 hrs is not observed, as a result of drift effects being stronger, with more efficient leaking of particles out of the injection flux tubes. It should be noted that turbulence-associated perpendicular transport, not included in our simulations, would smooth the difference in the intensity between the flux tubes connected to the acceleration region and those not connected.

Overall, 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 at higher energies is smaller than for lower energies, or may be related to an energy dependence of the mean free path. Both effects are not 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 Fig. 6. This is unlike the behaviour in 1D simulations, in which the slope is a strong function of the mean free path.

5. 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. In further work, we plan to extend the modelling to other PAMELA events.

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). In recently published work, Kocharov et al. (2020) presented an analysis of the 2017 September 10 GLE event using a 2D model of particle transport that includes perpendicular diffusion due to magnetic field turbulence. We expect that inclusion of such effects within our 3D model would smooth particle intensity profiles, eliminating discontinuities at low energies that are due to loss of connection to the flux tubes within which acceleration took place, while retaining the important effects of IMF polarity and HCS. We plan to carry out this study in future work.

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 at a given location 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 efficient outward HCS drift in the latter configuration. 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 observer dependent. Fluence spectra at 1 AU do not reflect the injection spectra. 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.

  • 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.

Therefore our analysis shows that 3D effects associated with the overall polarity of the heliospheric magnetic field and the presence of a large scale current sheet play a fundamental role in shaping relativistic solar proton propagation to Earth and the related 1 AU observables.

Acknowledgments

This work has received funding from the UK Science and Technology Facilities Council (STFC) (grant ST/R000425/1) and the Leverhulme Trust (Grant RPG-2015-094). We thank the PAMELA Collaboration for providing us with preliminary results concerning the 2012 May 17 event. G. A. de N. acknowledges support from the NASA/HSR grant NNH13ZDA001N-HSR, and the NASA/ISFM grant HISFM18. The work of J. G. was supported by NSF under Grant 1931252.

References

  1. Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2015, ApJ, 801, L3 [NASA ADS] [CrossRef] [Google Scholar]
  2. Battarbee, M., Dalla, S., & Marsh, M. S. 2017, ApJ, 836, 138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Battarbee, M., Dalla, S., & Marsh, M. S. 2018a, ApJ, 854, 23 [NASA ADS] [CrossRef] [Google Scholar]
  4. Battarbee, M., Guo, J., Dalla, S., et al. 2018b, A&A, 612, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Belov, A. V., Eroshenko, E. A., Kryakunova, O. N., Kurt, V. G., & Yanke, V. G. 2010, Geomag. Aeron., 50, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bieber, J. W., Evenson, P., Dröge, W., et al. 2004, ApJ, 601, L103 [Google Scholar]
  7. Bruno, A., Bazilevskaya, G. A., Boezio, M., et al. 2018, ApJ, 862, 97 [CrossRef] [Google Scholar]
  8. Burger, R. A. 2012, ApJ, 760, 60 [NASA ADS] [CrossRef] [Google Scholar]
  9. Chollet, E. E., Giacalone, J., & Mewaldt, R. A. 2010, J. Geophys. Res. (Space Phys.), 115, A06101 [NASA ADS] [Google Scholar]
  10. Cohen, C. M. S., & Mewaldt, R. A. 2018, Space Weather, 16, 1616 [CrossRef] [Google Scholar]
  11. Dalla, S., & Browning, P. K. 2005, A&A, 436, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Dalla, S., Marsh, M. S., Kelly, J., & Laitinen, T. 2013, J. Geophys. Res. (Space Phys.), 118, 5979 [NASA ADS] [CrossRef] [Google Scholar]
  13. Dalla, S., Marsh, M. S., & Laitinen, T. 2015, Astrophys. J., 808, 62 [NASA ADS] [CrossRef] [Google Scholar]
  14. Dalla, S., Marsh, M. S., & Battarbee, M. 2017a, ApJ, 834, 167 [NASA ADS] [CrossRef] [Google Scholar]
  15. Dalla, S., Marsh, M. S., Zelina, P., & Laitinen, T. 2017b, A&A, 598, A73 [Google Scholar]
  16. de Nolfo, G. A., Bruno, A., Ryan, J. M., et al. 2019, ApJ, 879, 90 [NASA ADS] [CrossRef] [Google Scholar]
  17. Fraschetti, F., & Jokipii, J. R. 2011, ApJ, 734, 83 [Google Scholar]
  18. Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204 [Google Scholar]
  19. Gopalswamy, N., Xie, H., Yashiro, S., et al. 2012, Space Sci. Rev., 171, 23 [NASA ADS] [CrossRef] [Google Scholar]
  20. Heber, B., Agueda, N., Bütikofer, R., et al. 2018, in Solar Particle Radiation Storms Forecasting and Analysis, eds. O. E. Malandraki, N. B. Crosby, et al., Astrophys. Space Sci. Lib., 444, 179 [NASA ADS] [CrossRef] [Google Scholar]
  21. Jokipii, J. R. 1966, ApJ, 146, 480 [Google Scholar]
  22. Jokipii, J. R., Levy, E. H., & Hubbard, W. B. 1977, ApJ, 213, 861 [CrossRef] [Google Scholar]
  23. Kelly, J., Dalla, S., & Laitinen, T. 2012, ApJ, 750, 47 [NASA ADS] [CrossRef] [Google Scholar]
  24. Klein, K. L., Tziotziou, K., Zucca, P., et al. 2018, in Solar Particle Radiation Storms Forecasting and Analysis, eds. O. E. Malandraki, N. B. Crosby, et al., 444, 133 [NASA ADS] [CrossRef] [Google Scholar]
  25. Kocharov, L., Vainio, R., Kovaltsov, G. A., & Torsti, J. 1998, Sol. Phys., 182, 195 [NASA ADS] [CrossRef] [Google Scholar]
  26. Kocharov, L., Pesce-Rollins, M., Laitinen, T., et al. 2020, ApJ, 890, 13 [CrossRef] [Google Scholar]
  27. Kota, J., & Jokipii, J. R. 1983, ApJ, 265, 573 [NASA ADS] [CrossRef] [Google Scholar]
  28. Laitinen, T., Kopp, A., Effenberger, F., Dalla, S., & Marsh, M. S. 2016, A&A, 591, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Laitinen, T., Dalla, S., Battarbee, M., & Marsh, M. S. 2018, in Space Weather of the Heliosphere: Processes and Forecasts, eds. C. Foullon, & O. E. Malandraki, IAU Symp., 335, 298 [NASA ADS] [Google Scholar]
  30. Li, G., & Lee, M. A. 2019, ApJ, 875, 116 [NASA ADS] [CrossRef] [Google Scholar]
  31. Marsh, M. S., Dalla, S., Kelly, J., & Laitinen, T. 2013, ApJ, 774, 4 [Google Scholar]
  32. Marsh, M. S., Dalla, S., Dierckxsens, M., Laitinen, T., & Crosby, N. B. 2015, Space Weather, 13, 386 [NASA ADS] [CrossRef] [Google Scholar]
  33. McCracken, K. G., Moraal, H., & Shea, M. A. 2012, ApJ, 761, 101 [NASA ADS] [CrossRef] [Google Scholar]
  34. Mewaldt, R. A., Cohen, C. M. S., Giacalone, J., et al. 2008, in American Institute of Physics Conference Series, eds. G. Li, Q. Hu, O. Verkhoglyadova, et al., Am. Inst. Phys. Conf. Ser., 1039, 111 [NASA ADS] [Google Scholar]
  35. Mewaldt, R. A., Looper, M. D., Cohen, C. M. S., et al. 2012, Space Sci. Rev., 171, 97 [NASA ADS] [CrossRef] [Google Scholar]
  36. Mishev, A., Usoskin, I., Raukunen, O., et al. 2018, Sol. Phys., 293, 136 [Google Scholar]
  37. Nitta, N. V., Liu, Y., DeRosa, M. L., & Nightingale, R. W. 2012, Space Sci. Rev., 171, 61 [NASA ADS] [CrossRef] [Google Scholar]
  38. Parker, E. N. 1965, Planet. Space Sci., 13, 9 [Google Scholar]
  39. Pei, C., Jokipii, J. R., & Giacalone, J. 2006, ApJ, 641, 1222 [Google Scholar]
  40. Potgieter, M. S., & Vos, E. E. 2017, A&A, 601, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Sáiz, A., Ruffolo, D., Rujiwarodom, M., et al. 2005, Int. Cosmic Ray Conf., 1, 229 [Google Scholar]
  42. Share, G. H., Murphy, R. J., White, S. M., et al. 2018, ApJ, 869, 182 [CrossRef] [Google Scholar]
  43. Strauss, R. D., Ogunjobi, O., Moraal, H., McCracken, K. G., & Caballero-Lopez, R. A. 2017, Sol. Phys., 292, 51 [CrossRef] [Google Scholar]

All Tables

Table 1.

Average number of 1 AU crossings per particle, , as a function of SEP proton kinetic energy E, for A+ and A configurations, and mean free path λ.

All Figures

thumbnail Fig. 1.

Maps of cumulative 1 AU crossings in a heliocentric coordinate system corotating with the Sun, for monoenergetic SEP proton populations, with energy as indicated in each panel. Left column: A+ configuration of the IMF; right column: A configuration. 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 vsw = 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).

In the text
thumbnail Fig. 2.

1 AU crossing count rates versus time summed over all heliographic longitudes and latitudes, for a monoenergetic proton population of initial kinetic energy E = 1 GeV and mean free path λ = 0.5 AU (left) and λ = 0.1 AU (right), for A+ and A configurations of the IMF. Other parameters of the simulations are as in Fig. 1.

In the text
thumbnail Fig. 3.

Proton count rates versus time for A+ (top) and A (bottom) configurations of the IMF, at a variety of 1 AU locations with respect to the best connected location ([0, 0]), for a power law population, for the proton energy ranges 100–400 MeV (blue) and 700–1000 MeV (green).

In the text
thumbnail Fig. 4.

Fluence energy spectra for A+ (top) and A (bottom) configurations of the IMF, at a variety of 1 AU locations with respect to the best connected location ([0, 0]), for a power law population. The solid lines in the [0, 0] panels give the slope of the injection spectrum. Parameters of the simulations are as in Fig. 3.

In the text
thumbnail 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 centre of the injection region is at N11W76 in this plot. The mean free path is λ = 0.3 AU.

In the text
thumbnail Fig. 6.

Time-intensity profiles of protons in four energy ranges, for the GLE 71 event. The left and middle panels show simulation results for an observer at the location of Earth, for λ = 1.0 AU and 0.3 AU respectively. The right panel shows preliminary data from PAMELA.

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.