Back-Tracing and Flux Reconstruction for Solar Events with PAMELA
Abstract
The PAMELA satellite-borne experiment is providing first direct measurements of Solar Energetic Particles (SEPs) with energies from 80 MeV to several GeV in near-Earth space. Its unique observational capabilities include the possibility of measuring the flux angular distribution and thus investigating possible anisotropies related to SEP events. This paper focuses on the analysis methods developed to estimate SEP energy spectra as a function of the particle pitch angle with respect to the Interplanetary Magnetic Field (IMF). The crucial ingredient is provided by an accurate simulation of the asymptotic exposition of the PAMELA apparatus, based on a realistic reconstruction of particle trajectories in the Earth’s magnetosphere. As case study, the results of the calculation for the May 17, 2012 event are reported.
1 Introduction
Solar Energetic Particles (SEPs) are high energy particles associated with explosive phenomena occurring in the solar atmosphere, such as solar flares and Coronal Mass Ejections (CMEs). They can significantly perturb the Earth’s magnetosphere producing a sudden increase in particle fluxes and, consequently, in the radiation exposure experienced by spacecrafts and their possible crew. SEPs constitute a sample of solar material and provide important information about the sources of the particle populations, and their angular distribution can be used to investigate details of the particle transport in the interplanetary medium.
SEP measurements are performed both by in-situ detectors on spacecrafts and by ground-based Neutron Monitors (NMs). While the former are able to measure SEPs with energies below some tens of MeV, the latter can only register the highest energy SEPs ( 1 GeV) during the so-called Ground Level Enhancements (GLEs). A large energy gap exists between the two groups of observations.
New accurate measurements of SEPs are being provided by the PAMELA experiment (Adriani et al., 2011a, 2015b). In particular, the instrument is able to detect SEPs in a wide energy interval ranging from 80 MeV up to several GeV, hence bridging the low energy data by space-based instruments and the GLE data by the worldwide network of NMs. In addition, PAMELA is sensitive to the particle composition and it is able to reconstruct the flux angular distribution, thus enabling a clearer and more complete view of the SEP events.
This paper reports the analysis methods developed for the estimate of SEP energy spectra with the PAMELA apparatus, as a function of the particle asymptotic direction of arrival. As case study, the results used in the analysis of the May 17, 2012 solar event (Adriani et al., 2015b) are presented.
2 The PAMELA Experiment
PAMELA is a space-based experiment designed for a precise measurement of the charged cosmic radiation in the kinetic energy range from some tens of MeV up to several hundreds of GeV (Picozza et al., 2007; Adriani et al., 2014). The Resurs-DK1 satellite, which hosts the apparatus, was launched into a semi-polar (70 deg inclination) and elliptical (350–610 km altitude) orbit on June the 15th 2006. In 2010 it was changed to an approximately circular orbit at an altitude of about 580 km. The spacecraft is 3-axis stabilized; its orientation is calculated by an onboard processor with an accuracy better than 1 deg. Particle directions are measured with a high angular resolution ( 2 deg). Details about apparatus performance, proton selection, detector efficiencies and measurement uncertainties can be found elsewhere (see e.g. Adriani et al. (2011b, 2013)).
3 Geomagnetic Field Models
The analysis of SEP events described in this work relies on the IGRF-11 (Finlay et al., 2010) and the TS07D (Tsyganenko & Sitnov, 2007; Sitnov et al., 2008) models for the description of the internal and external geomagnetic field, respectively: the former employs a global spherical harmonic implementation of the main magnetic field; the latter is a high resolution dynamical model of the storm-time geomagnetic field in the inner magnetosphere, based on recent satellite measurements. Consistent with the dataset coverage, the model is valid within the region delimited by the magnetopause (based on Shue et al. 1998) and by a spherical surface with a radius of 30 Earth’s radii (Re). Solar Wind (SW) and IMF parameters are obtained from the high resolution (5-min) Omniweb database (King & Papitashvili, 2004). The TS07D model is more flexible and accurate with respect to all past empirical models in reconstructing the distribution of storm-scale currents, so it is particularly adequate for the study of SEP events.
4 Back-tracing Techniques
Cosmic Ray (CR) cutoff rigidities and asymptotic arrival directions (i.e. the directions of approach before encountering the Earth’s magnetosphere) are commonly evaluated by simulations, accounting for the effect of the geomagnetic field on the particle transport (see e.g. Smart et al. 2000 and references therein). Using spacecraft ephemeris data (position, orientation, time), and the particle rigidity ( = momentum/charge) and direction provided by the tracking system, trajectories of all detected protons are reconstructed by means of a tracing program based on numerical integration methods (Smart & Shea, 2000, 2005), and implementing the afore-mentioned geomagnetic field models.
To reduce the computational time, geomagnetically trapped (Adriani et al., 2015a) and most albedo (Adriani et al., 2015c) particles (originated from the interactions of CRs with the Earth’s atmosphere) are discarded by selecting only protons with rigidities GV, where is the McIlwain’s parameter (McIlwain, 1966). Then, each trajectory is back propagated from the measurement location and traced, with no constraint limiting the total path-length or tracing time, until one of the two following conditions is satisfied:
-
1.
it reaches the model magnetosphere boundaries (see Section 3);
-
2.
or it reaches an altitude111Such a value approximately corresponds to the mean production altitude for albedo protons. of 40 km.
The two categories correspond to “allowed” and “forbidden” trajectories, respectively: the former includes contributions from Solar and Galactic CRs (hereafter SCRs and GCRs), while events satisfying the latter condition, including albedo particles with rigidities greater than , are excluded from the analysis.

5 Asymptotic Arrival Directions
The asymptotic arrival directions are evaluated with respect to the IMF direction, with polar angles and denoting the particle pitch and gyro-phase angle, respectively. Both Geographic (GEO) and Geocentric Solar Ecliptic (GSE) coordinates are used. Figure 1 reports some sample trajectories reconstructed in the GSE coordinate system. The total magnetic field intensities obtained with the IGRF-11 + TS07D models ( 30 ) are also shown (color coding) for the X-Y (top), the X-Z (middle) and the Y-Z (bottom) GSE planes. The crosses denote the estimated position of the bow shock nose (King & Papitashvili, 2004).

The trajectory analysis allows a deeper understanding of SEP events. To improve the interpretation of results, the directions of approach and the entry points at the model magnetosphere boundaries can be visualized as a function of particle rigidity and orbital position (Bruno et al., 2014). As an example, Figure 2 reports the proton results ( 3 GV) for the May 17, 2012 SEP event (Adriani et al., 2015b), associated with the first GLE of the 24rd solar cycle. Only the first PAMELA polar pass which registered the event is included, corresponding to the interval 01:57 – 02:20 UT.
Left panels show the reconstructed asymptotic directions for the selected proton sample (counts) in terms of GEO (top and middle) and GSE (bottom) coordinates; in the top panel colors denote the number of proton counts in each bin, while in middle and bottom panels they refer to the particle rigidity. Distributions are integrated over the polar pass. The spacecraft position is indicated by the grey curve. The contour curves represent values of constant pitch angle with respect to the IMF direction, denoted with crosses. As PAMELA is moving (eastward) and changing its orientation along the orbit, observed asymptotic directions rapidly vary performing a (clockwise) loop over the region above Brazil (see middle panel).
Right panels in the same figure display the distributions of asymptotic directions as a function of UT, and particle rigidity (top and middle) and pitch-angle (bottom); colors in the top panel denote the number of proton counts in each bin while, in middle and bottom panel, they correspond to particle pitch-angle and rigidity respectively. Solid curves denote the estimated Störmer vertical cutoff (Störmer, 1950; Shea et al., 1965) for the PAMELA epoch ( GV).
Since PAMELA aperture is about 20 deg, the observable pitch-angle range at a given rigidity is quite small (a few deg) except in the penumbral regions around the local geomagnetic cutoff, where particle trajectories become complex (chaotic trajectories) and both allowed and forbidden bands of CR trajectories are present (Cooke et al., 1991). Corresponding asymptotic directions rapidly change with particle rigidity and looking direction. Conservatively, these regions are excluded from the analysis.
6 Flux Evaluation
6.1 Apparatus Gathering Power
The factor of proportionality between fluxes and counting rates, corrected for selection efficiencies, is by definition the gathering power () of the apparatus:
(1) |
where is the solid angle domain limited by the instrument geometry, is the detector surface area, d is the effective element of area looking into , is the flux angular distribution (varying between 0 and 1) and is the directional response function (Sullivan, 1971).
In the case of the PAMELA instrument, is rigidity dependent due to the spectrometer bending effect on particle trajectories: it decreases with decreasing rigidity since particles with lower rigidity are more and more deflected by the magnetic field toward the lateral walls of the magnetic cavity, being absorbed before reaching the lowest plane of the Time of Flight system, which provides the event trigger (see Figure 3).

In terms of the zenith and the azimuth angles describing downward-going particle directions in the PAMELA frame222The PAMELA reference system has the origin in the center of the spectrometer cavity; the Z axis is directed along the main axis of the apparatus, toward the incoming particles; the Y axis is directed opposite to the main direction of the magnetic field inside the spectrometer; the X axis completes a right-handed system.:
(2) |
where is the apparatus response function in units of area (), and the factor accounts for the trajectory inclination with respect to the instrument axis.
For non-ideal detectors, it is necessary to account for the effect of elastic and inelastic particle interactions inside the apparatus (especially with its mechanical parts), on the measured coincidence rate. Consequently, is different for each particle species (protons, electrons, ions, etc.).
6.2 Isotropic Flux Exposition
For an isotropic particle flux, the gathering power does not depend on looking direction (i.e. = 1), and it is usually called the geometrical factor .
6.2.1 Monte Carlo Methods
A technically simple but efficient solution for the calculation of the geometrical factor of the apparatus is provided by Monte Carlo methods (Sullivan, 1971). The solid angle is subdivided into a large number of () bins, with the angular domain limited to downward-going directions. For each bin:
-
•
particles are produced with random position on a plane generation surface with area , placed just above the apparatus opening aperture, so that a large number of events corresponds to the intensity incident on the instrument.
-
•
For each particle, a random (i.e. random and ) direction is chosen in the selected () bin.
-
•
Trajectories are propagated through the apparatus and tested at successively deeper layers: only particles satisfying all the detector geometrical constraints defining the apparatus Field of View (FoV) are selected.
The procedure is repeated until the desired statistical precision (see below) is achieved. Then, for each rigidity value, the geometrical factor is obtained as:
(3) |
with
(4) |
where and are the numbers of selected and generated trajectories in each angular bin, and the factor accounts for the probability of a particle trajectory with direction ():
(5) |
The statistical uncertainty on can be evaluated333A more rigorous treatment is provided by Bayesian approaches. with binomial methods:
(6) |
Advantage is taken of the angular subdivision, by varying directions over the solid angle rather than a full hemisphere. In particular, a number of incident trajectories proportional to (with taken at bin center) is chosen in order to minimize the statistical error associated to the angular bins with a small .


The dependence of the instrument response on particle rigidity is studied by performing the calculation for 30 rigidities in the range 0.39–10 GV, and the value of at intermediate is evaluated through interpolation methods. An accurate estimate of the PAMELA geometrical factor based on the Monte Carlo approach can be found in Bruno (2008).
As an example, Figure 4 reports the ratio (color codes) over the PAMELA FoV as a function of local and , for 0.39 GV (top panel) and 4.09 GV (bottom panel) protons. The four peaks structure reflects the rectangular section of the apparatus; the differences in the FoV distributions at different rigidities are due to the bending effect of the magnetic spectrometer.
6.3 Anisotropic Flux Exposition
However, in presence of an anisotropic flux exposition (), the gathering power depends also on the flux angular distribution. SCR fluxes can be conveniently expressed in terms of asymptotic polar angles (pitch angle) and (gyro-phase angle) with respect to the IMF direction: . The corresponding gathering power can be written as:
(7) |
with = and =. The flux angular distribution is unknown.
6.3.1 Asymptotic Effective Area
For simplicity, we assume that SCR fluxes depend only on particle rigidity and asymptotic pitch angle , estimating an apparatus effective area (cm2) as:
(8) |
by averaging the directional response function over the angle. The method is valid also for isotropic fluxes (independent on ): in this case, the effective area is related to the geometrical factor by:
(9) |
The approach is analogous to the one developed for the measurement of geomagnetically trapped protons (Bruno et al., 2015), with and denoting the polar angles with respect to the local geomagnetic field direction. But while local angles can be calculated from by means of basic rotation matrices of the FoV involving only the spacecraft orientation with respect to the local magnetic field, the conversion to asymptotic angles depends on the particle propagation in the magnetosphere, so trajectory tracing methods are needed.
6.3.2 Area Calculation
The effective area definition given in Equation 8 is based on the assumption of approximately isotropic fluxes within small pitch-angle bins. Consequently, can be derived from Equation 3 by integrating the directional response function over the () directions corresponding to pitch angles within the interval :
(10) |
The approach accuracy depends on the number of () bins used in the angular partitioning (see Section 6.2.1), while the width of the bins is chosen accounting for the detector angular resolution.


To convert local () into asymptotic () directions and apply Equation 10, a large number of trajectories , uniformly distributed inside PAMELA field of view, has to be reconstructed in the magnetosphere, for each rigidity value and each orbital position. In order to assure a high resolution, the calculation is performed for time steps with a 1-sec width, back-tracing about 2800 trajectories for each rigidity bin, for a total of more than trajectories for each polar pass ( 23 min). At a later stage, results are extended over the full field of view of PAMELA through a bilinear interpolation.
The procedure is illustrated in Figures 5-6 for 0.39 and 4.09 GV protons respectively, at a sample orbital position (May 17, 2012, 02:07 UT). Top panels report the distribution of reconstructed directions in the PAMELA field of view, with each point associated to a given asymptotic direction (,); middle panels show the calculated (after interpolation) pitch-angle coverage444The calculation can be neglected under the gyro-tropic approximation.; bottom panels illustrate the estimated effective area as a function of the explored pitch-angle range. Final calculation results for 22 rigidity values between 0.39 and 4.09 GV (see the color code) are displayed in Figure 7: the peaks in the distributions correspond to vertically incident protons.
Figure 8 reports the asymptotic cones of acceptance of the PAMELA apparatus evaluated for the first polar pass (01:57–02:20 UT) during the May 17, 2012 SEP event (Adriani et al., 2015b). Results for sample rigidity values (see the color code) are shown as a function of GEO (top panel) and GSE (middle panel) coordinates; grey points denote the spacecraft position, while crosses indicate the IMF direction. The pitch-angle coverage as a function of orbital position is shown in the bottom panel. While moving (and rotating) along the orbit, PAMELA covers a large pitch-angle interval, approximately ranging from 0 to 145 deg. In particular, PAMELA is looking at the IMF direction between 02:14 and 02:18 UT, depending on the proton rigidity.


6.4 Differential Directional Fluxes
Differential directional fluxes are obtained at each orbital position as:
(11) |
where is the number of proton counts in the bin , corrected by the selection efficiencies, and the denominator represents the asymptotic exposition of the apparatus integrated over the selected rigidity bin .
Averaged fluxes over the polar pass are evaluated as:
(12) |
where and the exposition is derived by weighting each effective area contribution by the corresponding lifetime spent by PAMELA at the same orbital position.
6.5 GCR Background Subtraction
Since it is not possible to discriminate between SCR and GCR signals, solar flux intensities are obtained by subtracting the GCR contribution from the total measured flux. The GCR component is evaluated by estimating proton fluxes during two days prior the arrival of SEPs. We found that the GCR flux is isotropic with respect to the IMF direction within experimental uncertainties. Consequently, the same flux is subtracted for all pitch angle bins:
(13) |
6.6 Flux Uncertainties
Flux statistical errors can be obtained by evaluating 68.27% C.L. intervals for a poissonian signal in presence of a background (Feldman & Cousins, 1998). Systematic uncertainties related to the reconstruction of asymptotic directions are evaluated by introducing a bias in the particle direction measurement from the tracking system, according to a gaussian distribution with a variance equal to the experimental angular resolution.
7 Conclusions
This paper reports the analysis methods developed for the estimate of SEP energy spectra as a function of the particle asymptotic direction of arrival. The exposition of the PAMELA apparatus is evaluated by means of accurate back-tracing simulations based on a realistic description of the Earth’s magnetosphere. As case study, the results of the calculation for the May 17, 2012 event are discussed. The developed trajectory analysis enables the investigation of flux anisotropies, providing fundamental information for the characterization of SEPs. It will prove to be a vital ingredient for the interpretation of solar events observed by PAMELA during solar cycles 23 and 24.
Acknowledgements
We acknowledge support from The Italian Space Agency (ASI), Deutsches Zentrum fr Luftund Raumfahrt (DLR), The Swedish National Space Board, The Swedish Research Council, The Russian Space Agency (Roscosmos) and The Russian Scientific Foundation. We gratefully thank N. Tsyganenko for helpful discussions, and M. I. Sitnov and G. K. Stephens for their assistance and support in the use of the TS07D model.
References
- Adriani et al. (2011a) Adriani, O., et al., 2011a, “Observations of the 2006 December 13 and 14 solar particle events in the 80 MeV n-1 – 3 GeV n-1 range from space with the PAMELA detector”, ApJ 742:102, doi:10.1088/0004-637X/742/2/102.
- Adriani et al. (2011b) Adriani, O., et al., 2011b, “PAMELA Measurements of Cosmic-Ray Proton and Helium Spectra”, Science, Vol. 332 no. 6025 pp. 69–72, doi:10.1126/science.1199172.
- Adriani et al. (2013) Adriani, O., et al., 2013, “Time dependence of the proton flux measured by PAMELA during the 2006 July – 2009 December solar minimum”, ApJ 765:91.05205, doi:10.1088/0004-637X/765/2/91.
- Adriani et al. (2014) Adriani, O., et al., 2014, Physics Reports, Volume 544, Issue 4, Pages 323–370, doi: 10.1016/j.physrep.2014.06.003.
- Adriani et al. (2015a) Adriani, O., et al., 2015a, “Trapped proton fluxes at low Earth orbits measured by the PAMELA experiment”, ApJ 799 L4, doi:10.1088/2041-8205/799/1/L4.
- Adriani et al. (2015b) Adriani, O., et al., 2015b, “PAMELA’s Measurements of Magnetospheric Effects on High Energy Solar Particles”, ApJ 801 L3, doi:10.1088/2041-8205/801/1/L3.
- Adriani et al. (2015c) Adriani, O., et al., 2015c, “Reentrant albedo proton fluxes measured by the PAMELA experiment”, J. Geophys. Res. Space Physics, 120, doi:10.1002/2015JA021019.
- Bruno (2008) Bruno, A., 2009, Ph.D. thesis, University of Bari, Bari, Italy; http://pamela.roma2.infn.it/.
- Bruno et al. (2014) Bruno, A., et al., 2014, “Solar energetic particles measured by PAMELA”, 3rd Space Radiation and Plasma Environment Monitoring Workshop, ESTEC, NL, 13–14 May.
- Bruno et al. (2015) Bruno, A., et al., 2015, “PAMELA’s measurements of geomagnetically trapped and albedo protons”, PoS(ICRC2015)288.
- Cooke et al. (1991) Cooke D. J., Humble, J. E., Smart, M. A., et al., 1991, Il Nuovo Cimento, 14C, 213.
- Feldman & Cousins (1998) Feldman, G. J., & Cousins, R. D., 1998, Phys. Rev. D57 3873–3889.
- Finlay et al. (2010) Finlay, C. C., et al., 2010, “International Geomagnetic Reference Field: the eleventh generation”, Geophysical Journal International, 183: 1216–1230.
- King & Papitashvili (2004) King, J. H., & Papitashvili, N. E., 2004, “Solar wind spatial scales in and comparisons of hourly Wind and ACE plasma and magnetic field data”, J. Geophys. Res., Vol. 110, No. A2, A02209, http://omniweb.gsfc.nasa.gov/.
- McIlwain (1966) McIlwain, C., 1966, Space Sci. Rev. 5, 585–589.
- Picozza et al. (2007) Picozza, P., et al., 2007, “PAMELA - A Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics”, Astropart. Phys., Vol 27, Pages: pp. 296–315.
- Shea et al. (1965) Shea, M. A., Smart, D. F., and McCracken, K. G., 1965, J. Geophys. Res. 70, 4117-4130.
- Shue et al. (1998) Shue, J.-H., P. Song, C. T. Russell, J. T. , et al., 1998, “Magnetopause location under extreme solar wind conditions”, J. Geophys. Res., 103, 17, 691.
- Smart et al. (2000) Smart, D.F., Shea, M.A., and Flckiger, E.O., 2000, “Magnetospheric Models and Trajectory Computations”, Space Science Reviews 93: 305–333.
- Smart & Shea (2000) Smart, D. F., & Shea, M. A., 2000, Final Report, Grant NAG5–8009, Center for Space Plasmas and Aeronomic Research, The University of Alabama in Huntsville.
- Smart & Shea (2005) Smart, D. F., & Shea, M. A., 2005, “A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft”, Adv. Space Res., 36, 2012–2020.
- Störmer (1950) Störmer, C., 1950, “The Polar Aurora”, Oxford University Press, London.
- Sullivan (1971) Sullivan, J. D., 1971, “Geometric factor and directional response of single and multi-element particle telescopes”, Nucl. Instr. and Meth. 95, 5.
- Tsyganenko & Sitnov (2007) Tsyganenko, N. A., & Sitnov, M. I., 2007, “Magnetospheric configurations from a high-resolution data-based magnetic field model”, J. Geophys. Res., 112, A06225.
- Sitnov et al. (2008) Sitnov, M. I., Tsyganenko, N. A., Ukhorskiy, A. Y., and Brandt, P. C., 2008, “Dynamical data-based modeling of the storm-time geomagnetic field with enhanced spatial resolution”, J. Geophys. Res., 113, A07218.