Cosmic ray cutoff prediction using magnetic field from global magnetosphere MHD simulations

. Relativistic particles entering the Earth’s magnetosphere, i.e. cosmic rays and solar energetic particles, are of prime space weather interest because they can affect satellite operations, communications, and the safety of astronauts and airline crews and passengers. In order to mitigate the hazards that originate from such particles one needs to predict the cutoff latitudes of such particles as a function of their energies and the state of the magnetosphere. We present results from a new particle tracing code that is used to determine the cutoff latitudes of 8–15 Me n − 1 alpha particles during the 23/24 April, 1998 geomagnetic storm and the preceding quiet time. The calculations are based on four different geomagnetic ﬁeld models and compared with SAMPEX observations of alpha particles in the same energy range. The geomagnetic ﬁeld models under consideration are: (i) the International Geomagnetic Reference Field (IGRF) model, (ii) the Tsyganenko “89” model (T89c), (iii) the Tsyganenko “96” model (T96), and (iv) a global magnetohydrodynamic (MHD) model of Earth’s magnetosphere. Examining 11 SAMPEX cutoff latitude observations we ﬁnd that the differences between the observed and the predicted cutoff latitudes are 2.3 ◦ ± 2.0 ◦ (mean) and 7.9 ◦ (maximum difference) for the IGRF model; 3.9 ◦ ± 2.4 ◦ (mean) and 6.9 ◦ (maximum difference) for the T89c model; 4.0 ◦ ± 1.4 ◦ (mean) and 5.5 ◦ (maximum difference) for the T96 model; and 2.5 ◦ ± 1.7 ◦ (mean) and 7.0 ◦ (maximum difference) for the MHD model. All models gen-erally predict cutoff latitudes equatorward of the SAMPEX observations. The MHD model results also show steeper cutoff energy gradients with latitude compared to the empirical models and more structure in the cutoff energy versus latitude function, presumably due to the presence of boundary layers in the MHD model.


Introduction
Space Radiation is a serious hazard for satellite operations, communications, and human space flights. A type of space radiation known as solar energetic particle (SEP) event (i.e. particles with E>1 MeV at 1 AU) becomes a concern because under nominal shielding or while on extravehicular activity, individual SEP events can exceed exposure limits over a period of a few hours. Likewise, such radiation is of concern for the safety of airline passengers and crews on routes that reach high geomagnetic latitudes. Thus it is important to be able to map latitudinal cutoffs of SEP events. The latitudinal cutoff is the lowest latitude at which ions less then a specific rigidity cannot precipitate below. The cutoff energies reach a maximum near the equatorial regions and a minimum near the geomagnetic poles. However, the latitudinal cutoff is not a simple boundary, but a range of latitudes that is called the penumbra. The penumbra region consists of latitudinal bands of forbidden and allowed particle trajectories. In the direction of descending energy, the last allowed trajectory before the first forbidden trajectory is called the upper cutoff while the last allowed trajectory before the last forbidden trajectory is called the lower cutoff. The weighted average between those cutoff latitudes is the effective cutoff latitude. However, for many applications, including those presented in this paper, it is not important to consider the penumbra region in detail because other effects play a more important role. A more detailed discussion of the characteristics of geomagnetic cutoffs is given by Cooke et al. (1991).
The geomagnetic field model is the primary factor in the theoretical determination of the latitudinal cutoffs. Calculations of the cutoff latitudes in previous studies used the International Geomagnetic Reference field (IGRF) model (Peddie, 1982;Barton, 1997;Sabaka et al., 1997) and the Tsyganenko magnetic field models (Tsyganenko, 1989(Tsyganenko, , 1995(Tsyganenko, , 1996, which will be referred to as the T89c and T96 models. These models have been used to calculate the latitudinal cutoffs for the international space station during quiet and active times (Smart et al., 1999a(Smart et al., , 1999b. With the increasing importance of space weather effects and the emergence of new magnetic field models there is a renewed interest in cutoff calculations. Recent studies have begun to compare observed geomagnetic cutoff latitudes with those calculated using different magnetic field models. For example, Leske et al. (2001) determined the geomagnetic cutoff latitude as a function of K p with the 8 Me n −1 to 15 MeV n −1 channels of the SAMPEX Mass Spectrometer Telescope (MAST) and 20 MeV to 29 MeV proton channels and 29 MeV to 64 MeV proton channels of the SAMPEX Proton/Electron Telescope (PET) for six SEP events. Kahler and Ling (2000) compared the SAMPEX/PET proton count rate cutoff latitudes for energies between 20 MeV to 29 MeV and 29 MeV to 64 MeV with cutoff predictions calculated using the IGRF geomagnetic magnetic field model and the T89c geomagnetic magnetic field model. That study found that the measured results were on average different from the IGRF model by about 7 • ±3 • in latitude and were different from the T89c model by about 2.55 • ±2.45 • latitude. Such differences are significant, in particular because the differences were much larger for some passes than the averages would indicate.
The present study essentially extends the work of Kahler and Ling (2000) by including more magnetic field models in the comparison. Specifically, our study compares the 1995 IGRF, T89c, T96, and MHD geomagnetic field model cutoff energy results with one of the events in the  work. The primary goal of our study is to compare the cutoff latitudes as measured by SAMPEX/MAST with those determined from relativistic particle tracing in a magnetic field from a global magnetosphere MHD simulation during a solar energetic particle event on 23 April 1998 at ∼22:40 UT. Although global MHD models have their own limitations (Raeder, 2003), they are expected to better represent the magnetic field in the vicinity of the polar cap, as well as the polar cap locations, in particular during geomagnetic active times. In addition, we compare the SAMPEX results with the cutoff latitudes determined from the T96 magnetosphere magnetic field model, which has not previously been used to determine cutoff energies and latitudes.

Particle tracing and magnetic field models
In order to determine the cutoff energies of relativistic particles it is most efficient to calculate the reverse trajectories of particles starting from the top of Earth's atmosphere at an altitude of 20 km (Smart and Shea, 1994). Determining the cutoff energy from the magnetopause to Earth's upper atmosphere would be far more time consuming because the trajectories would need to be calculated from many positions on the magnetopause, and for a range of pitch angles and phase angles, but most of these particles would never reach the upper atmosphere.
Thus, for a given latitude and longitude we launch particles along the field lines upward, i.e. at 0 • or 180 • pitch angle and arbitrary phase angle, and determine their fate. The method of launching particles up the field lines applies at the polar region that this study is concerned with, but not at lower latitude locations. For the empirical models the choice of the phase angle did not make any significant differences and thus only one arbitrary phase angle was chosen. The MHD model magnetic field produced results that were more sensitive to the phase angle; thus we repeated each calculation for 20 randomly chosen phase angles and averaged the results. To trace the model particle's path we use the equations of a relativistic charged particle in a magnetic field. In a Cartesian coordinate system, these equations result in six simultaneous differential equations for the six unknown components of the particles' momentum p: it's location r: together with the relativistic relation between the momentum and the velocity v: where B is the magnetic field, m is the particle's mass, and c is the speed of light. We solve these equations with a fourth order Runge-Kutta integration scheme. Other investigators have used higher order schemes but we find them not beneficial. In Appendix A we explain the rationale for using this numerical scheme and how we control the numerical integration error.
In order to find the cutoff energy for a given latitude and longitude we need to calculate the trajectories of an ensemble of particles with varying energy. The trajectory calculation of the particles begins at an altitude of 20 km above the Earth's surface, which is an altitude above which particle collisions resulting in nuclear cascades would be minimal. Particles that reach a distance of 15 R E from Earth are considered on an allowed trajectory. On the other hand, particles that return to Earth to within less than 450 km altitude are considered to be on a forbidden trajectory. Such a particle would with near certainty collide with a neutral and precipitate. Likewise, particles that have traveled more than 500 R E without having met either of the previous criteria are also considered to be on a forbidden trajectory. Such particles are usually on closed drift paths around the Earth. In order to minimize the number of particle trajectories that need to be calculated we use bisection along the energy coordinate. The bisection method examines the midpoint of an energy interval where the interval is bracketed by particles with low energies (i.e. trapped particles) and high energy (i.e. particles that escape). If the particles with the energy at the midpoint of the interval escape, then this energy becomes the new upper interval bracket. If the particles of the energy at midpoint of the interval are trapped this energy becomes the new lower interval bracket. This process continues in a loop until the interval reaches some minimum size, which is less than 0.0001 of the difference between the upper and lower bracket divided by the midpoint of the interval (see Press et al. (1997) for details on the bisection method.) Using bisection has the principal downside that it does not take into account the penumbra. Although the bisection method will not fail, it will return a cutoff at an arbitrary location within the penumbra. We do not consider this uncertainty detrimental to our results for two reasons: First, because of the existence of the penumbra, it is already difficult to define a unique cutoff, and second, other errors, for example in the magnetic field interpolation, have greater effects on the cutoff energy (or the cutoff latitude for a given energy) than the penumbra effect.
The essence of this study is the use a global numerical model of Earth's space environment that is principally based on a MHD description of the plasma as the basis for Earth's magnetic field. Part of this simulation code is an ionospheric model for closure of field-aligned currents. The model solves the ideal MHD equations for the magnetosphere and a potential equation for the ionosphere. Numerical effects, such as diffusion, viscosity, and resistivity, are necessarily introduced by the numerical methods, which are discussed in detail by Raeder (2003). A pure dipole model is used inside of 5 R E , a linearly weighted superposition of the MHD and interior dipole fields is used between spherical boundaries at 5 R E and 7 R E , and the MHD results are used outside of 7 R E . This is necessary because within a few R E from Earth the dipole field is not well resolved by the numerical MHD grid. The boundaries of the blending region can be controlled, however, for this study they are set at L=5 to L=7. The cutoff latitudes for a particular energy are dependent on the values for the blending region. As the L shells for the blending region are decreased the cutoff latitudes for a specific energy decrease. The required model inputs are the solar wind magnetic field vectors, density, solar wind speed, and dynamic pressure (Raeder et al., 2001). The output relevant to our study is the MHD model's magnetic field which is then used to integrate particle trajectories for the cutoff determination.
The empirical magnetic field models developed by Tsyganenko (1989Tsyganenko ( , 1995Tsyganenko ( , 1996 provide an often used alternative to using either dipole/IGRF or the global MHD model representations of the Earth's magnetic field. The T89c (Tsyganenko, 1989) and T96 (Tsyganenko, 1996) models are a semi-empirical best-fit representation for the geomagnetic field, based on a large number of satellite observations (IMP, HEOS, ISEE, POLAR, Geotail, etc.). The models include the contributions from external magnetospheric sources: ring current, magnetotail current system, magnetopause currents and large-scale system of field-aligned currents (Tsyganenko, 1989(Tsyganenko, , 1995(Tsyganenko, , 1996. The T89c model requires both K p and dipole tilt information in order to determine the magnetic field vector. The T96 model, on the other hand, is parameterized using the values of D st , IMF B z and IMF B y to determine the magnetic field vector. Even though the Tsyganenko models are parameterized using parameters that relate to geomagnetic activity (D st , K p , IMF B z ) they remain statistical representations that are unlikely to capture the Earth's magnetic field truthfully on short time scales. We thus expect that during times of high geomagnetic activity, like the ones chosen in this study, the MHD based calculations will be closer to reality compared to those based on empirical models. This statement is based on the fact that empirical models like IGRF or Tsyganenko lack principal magnetospheric features, such as a bow shock, a magnetosheath, the cusps, and in the case of the IGRF, even a tail.
In addition to the Tsyganenko model a separate portion, the IGRF, calculates contribution from Earth's internal field. The IGRF model uses a spherical harmonics expansion of the scalar potential in geocentric coordinates. The IGRF model coefficients are based on all available data sources including geomagnetic measurements from observatories, ships, aircrafts and satellites. Other than the date, there are no other input values required for this model, because the internal magnetic field does not vary in response to solar wind or geomagnetic activity conditions.

Data
The relativistic alpha particle precipitation was measured by the Solar, Anomalous, and Magnetospheric Particle Explorer (SAMPEX) spacecraft, which was launched into a 520 by 670 km orbit with an inclination of 82 • in July 1992 (Baker et al., 1993). SAMPEX is a momentum biased, Sun pointed spacecraft that maintains the experiment-view axis in a zenith direction. The data that we use for the model comparisons were obtained by the Mass Spectrometer Telescope (MAST) (Cooke et al., 1993). This instrument measures counts per second of MeV particles within an energy range of 7 Me n −1 to 450 MeV n −1 and has an acceptance cone of approximately 100 • . Cutoff latitudes were previously determined with the MAST instrument using the Z2 rate (count rate associated with 8 Me n −1 to 15 Me n −1 ) in Leske et al. (2001). In that study the cutoff latitude was defined to be the location at which the count rate falls below half of its mean value above 70 • . Using this method the cutoff latitude is known to within ∼0.2 • . In the following, we use the Leske et al. (2001) definition of the cutoff latitude when we compare our model results with the SAMPEX data.  Figure 1a displays the cutoff energies on 23 April 1998 at 16:07 UT before a solar energetic particle event that occurs at about 22:40 UT on 23 April. Figure 1b displays the cutoff energies on 24 April at 06:48 UT after the solar energetic particle event. The color bar indicates red for high cutoff energies and dark blue for low cutoff energies. There are no contours above about 75 • latitude because only cutoff energies above 1 MeV are considered.  Leske et al. (2001) and predicted with the 1995 IGRF model with the coefficients applicable for April events in 1998. The first column gives the date and time, the second column is the K p , and the third column is the Leske et al. (2001) cutoff latitude in Solar Magnetic coordinates. The next six columns give the cutoff latitude in Solar Magnetic coordinates and the difference in cutoff latitude compared to the Leske et al. (2001) value for three different cutoff energies per nucleon: 8 MeV n −1 , 11.5 MeV n −1 , and 15 MeV n −1 , respectively. The last row of the table displays the mean of the magnitude of the latitudinal difference and the standard deviation of the mean. axis in the north and the Y axis is perpendicular to the Earth-Sun line. The X axis is not always directed toward the Sun and rocks back and forth through 11.5 • about the Earth-Sun line. Figure 1a displays the cutoff energies before a solar energetic particle event that occurs at about 22:40 UT on 23 April and Fig. 1b displays the cutoff energies after the event.
The color bar indicates red for high cutoff energies and dark blue for low cutoff energies. No contours exist above about 75 • latitude because only cutoff energies above 1 MeV are considered. The contours appear to be shifted since the solar magnetic system is not fixed with respect to the Sun but shifts with daily rotation of the magnetic dipole axis around the polar axis and with Earth's orbital motion around the Sun. The results are very similar to the results of Smart et al., (1999a,b). Furthermore, a careful examination of Figs 1a and b shows that the contour plots are the same except for the longitudinal shift because the IGRF model does not depend on geomagnetic activity or solar wind conditions and thus the cutoff energies will not change. Table 1 compares the measured cutoff latitude determined by the SAMPEX/MAST instrument in the 8 MeV n −1 to 15 MeV n −1 range with the cutoff latitude calculated using the IGRF 1995 model. The first column gives the date and time and the second column givens the K p values for each event. Three different cutoffs for the IGRF model are given for 8 MeV n −1 , for 11.5 MeV n −1 , and for 15 MeV n −1 , respectively. The three different energies per nucleon were examined because particles from 8 MeV n −1 to 15 MeV n −1 are included in the measurements, but the model has a finer resolution. Examining the three energies also demonstrates the variation in the cutoff latitude associated with the large energy range of the measurements. Of course the degree of the variability depends on the model as we will show later. The mean of the magnitude of the latitudinal difference between the measured data and the model are for 8 MeV: 2.2 • ±2.6 • , for 11.5 MeV: 2.3 • ±2.0 • , and for 15 MeV: 2.6 • ±1.8 • . The largest difference between the model and the data is 7.9 • for an energy of 11.5 MeV n −1 which occurred at 08:08 UT on 24 April. The mean differences are similar to the statistical study of Kahler and Ling (2002), which found a difference of about 6.56 • ±2.64 • for the IGRF model and the 20 MeV to 29 MeV channels of the SAMPEX/PET instrument. The differences between the model results and the observed cutoff latitudes can be explained by the independence of the IGRF model with geomagnetic activity. That is to say, the cutoff latitude found by the model does not vary with the change in K p , but only with local time. Note that in Table 1 the difference between the measure values and the model values are significant for all K p values (even low values of K p ). It is also important to note that the maximum difference, which characterizes the model with respect to space weather predictions, is much larger than the mean and reaches values up to ∼8 • . Table 2 show the cutoff latitudes for the same event as determined with the T89c model. The figures and the table have the same format as Fig. 1 and Table 1, respectively, of the previous section. Between 16:07 UT on 23 April and 06:48 UT on 24 April the K p index value changed from 1.3 to 5.7. This change is reflected in the change in the cutoff latitudes and the fact that particles of a given energy have access to a larger area in the polar regions. We also find that the structure of the contours at auroral latitudes is more complex and not as smooth as the one that was obtained with  Fig. 1 and shows the cutoff energies per nucleon that were determined with the T89c model. After the solar energetic particle event the cutoff latitudes have significantly decreased (Fig. 2b) due to the increase in K p . .5 • . The most significant difference between the model and the data is −7.0 • for an energy of 11.5 MeV n −1 and occurred at 05:10 UT on 24 April. The T89c differences are larger than those that were found with the IGRF model, but they are not statistically significant. It is also important to note that the T89c based cutoff latitude is systematically lower than the observed one, which indicates that the T89c model's polar cap area is systematically too large with a few exceptions. Table 3 show the cutoff energies as determined with the T96 model, also in the same format as Fig. 1 and Table 1, respectively. The results closely match those obtained with the T89c model. Specifically, the mean of the magnitude of the latitudinal difference between the measured cutoff latitude and the cutoff latitude derived from the model are for 8 Me n −1 : 3.2 • ±1.5 • , for 11.5 Me n −1 : 4.0 • ±1.4 • , and for 15 Me n −1 : 4.5 • ±1.6 • . The most significant difference between the model and the data is -5.5 • for an energy of 11.5 Me n −1 and occurred at 06:48 UT on 24 April. Thus, the average differences are slightly higher than those that were obtained with the T89c model. Also, as with the T89c model, the mis-prediction of the cutoff latitude is systematically too low, with one exception. Table 4 show the cutoff energies as determined with the MHD model. Clearly at the low latitudes the cutoff energy contours do not match those observed in the other three magnetic field models. This is due to the blending of the MHD model with a magnetic dipole model between L=5 to L=7. The dipole model lacks the higher order moments (quadrupole and up) that make a big difference at the lower latitudes and at the higher energies. However, the cutoff energies at auroral latitudes should not be much affected, because the magnetic field topology at those latitudes is primarily determined by the external currents in the magnetosphere. Another difference between the results of the IGRF, T89c, and T96 models and the MHD model results is the number of particles used in each model. The results from the other models are the results from a single alpha particle. The MHD results are an average of 20 different alpha particles. The difference between each of the 20 particles is the phase angle, which was randomly selected for each alpha particle. The mean of the magnitude of the latitudinal difference between the measured cutoff latitudes and the predicted latitudes are for 8 Me n −1 : 1.6 • ±1.7 • , for 11.5 Me n −1 : 2.5 • ±1.7 • , and for 15 Me n −1 : 3.3 • ±2.0 • , which is significantly less than for any of the empirical models. The largest difference between the model and the data is −7.0 • for an energy of 11.5 Me n −1 and occurred at 22:39 UT on 23 April. As in the case of the empirical models, the predicted cutoff latitude is in all cases except one systematically too low, and thus the polar cap is also too large in the MHD model.  Fig. 1 and shows the cutoff energies that were determined with the T96 model. After the solar energetic particle event the cutoff latitudes have significantly decreased due to the decrease in D st and increase in solar wind pressure similar to the results obtained with the T89c model. the 8 Me n −1 to 15 Me n −1 channels of SAMPEX/MAST. The dotted horizontal lines show the energies of 8 MeV n −1 , 11.5 MeV n −1 , and 15 MeV n −1 for reference. The panels graphically demonstrate the differences between the models. Note that the models are consistently predicting lower cutoff latitudes than observed except for at 24 April at 08:08 UT when the observed cutoff is depressed equatorward of the model predicted values. This is opposite to the findings of Smart and Shea (1994) and Ogliore (2001).

Comparison of all four models
The difference between models and data has been previously observed with HEAO-3 data (Smart and Shea, 1994;Ogliore et al., 2001), and with SAMPEX/PET (Kahler and Ling, 2002). Smart and Shea (1994) found the difference to be on the order of 5%. Kahler and Ling (2002) found differences of 6.56 • ±2.64 • between the SAMPEX/PET 20 MeV to 29 MeV proton data and the IGRF model based results, but differences of 2.47 • ±1.89 • with the T89c model. These differences were even smaller with the SAMPEX/PET 29 MeV to 64 MeV proton data 5.37 • ±2.62 • and 1.76 • ±1.54 • for the IGRF and T89c models, respectively. In both energy ranges the IGRF model consistently overestimated the cutoff latitude, but the T89c model fluctuated about a mean difference of about 0 • . The higher energy range most likely gives a better estimate of the cutoff latitude due to the larger gradient in the cutoff energy with latitude in that energy range. For the event examined in this study the T89c consistently underestimated the cutoff latitude. It is not yet clear why the model consistently underestimates the cutoff latitude.
As we have shown the agreement between the MHD model and the SAMPEX data is similar to the agreement between the other models. This agreement is influenced by the blending of the MHD model with the dipole field. Decreasing or increasing the L values of the blending region decreases or increases the cutoff energy latitude. A detailed analysis of the blending region for the MHD model and the dipole would further improve the results, but a more realistic improvement would be to combine the IGRF model and the MHD model. This has not yet been done in order to separate the MHD model results from the IGRF model and to illustrate the features of the MHD model. This combination will be performed within a later study.
In addition to the model and measured cutoff latitude differences, the panels in Fig. 5 frequently show a "band" of slightly higher cutoff energies just poleward of the sharpest decrease in cutoff energy in some of the models. For the empirical models this band may be a crude reproduction of the cosmic ray penumbra region. The penumbra is most pronounced in the T89c model. The penumbra is small, but present in the IGRF model. Both the MHD and T96 model show no penumbra at all. The MHD model most likely shows little or no penumbra because it is an average of the cutoff energies for 20 particles. An examination of the cutoff energy versus latitude for individual alpha and proton particles displays a more complicated structure. The T96 model may not exhibit any penumbra because it has smoothly varying fields. Examining the detailed cause of the model penumbra is possible by careful analysis of the particle trajectories but it beyond the scope of this paper and will be presented elsewhere.

Summary and conclusions
We used 4 different models (IGRF, T89c, T96, and MHD) to calculate the cutoff latitudes for precipitating alpha particles in the 8 MeV n −1 to 15 MeV n −1 range for the storm event of 23/24 April 1998. We compared the predicted cutoff latitudes with measured cutoff latitudes that were determined using the 8 MeV n −1 to 15 MeV n −1 channels of the SAM-PEX/MAST instrument . Table 5 summarizes our results. The table shows that mean difference for a particle energy of 11.5 MeV n −1 between the model and the data in cutoff latitude, the standard deviation, and the maximum difference in latitude is largest for the IGRF model and smallest for the T96 model. Within the standard deviation  Fig. 1 and shows the cutoff energies that were determined with the MHD magnetic field model. Note that the high latitude cutoffs for the 23 April 16:07 UT case are significantly more complex than those in the IGRF, T89c, and T96 models. Fig. 5. Cutoff energies versus solar magnetic latitude for the four different magnetic field models at the longitude of the SAMPEX spacecraft for 3 different times. The blue curve is the IGRF model, the red curve is the T89c model, the black curve is the T96 model, and the green curve is the MHD model. The vertical dashed line is the cutoff latitude determined by  for the 8 Me n −1 to 15 Me n −1 channels of SAMPEX/MAST. The dotted lines show the energies of 8 MeV n −1 , 11.5 MeV n −1 , and 15 MeV n −1 . The panel graphically demonstrates the differences between the models. the IGRF, T89c, T96, and MHD model give approximately the same results. In general, however, all models estimate the cutoff latitude equatorward of the SAMPEX observations.
Our calculations for the IGRF and T89c model confirm the earlier result of the Kahler and Ling (2002) study, which used the 20 MeV to 29 MeV proton channels (IGRF: 6.56 • ±2.64 • and T89c: 2.55 • ±1.89 • ) of the SAMPEX PET. As far as we are aware, this study is the first comparison of T96 and Thus, these models either mis-predict the polar cap boundary itself or the magnetic field topology in the vicinity of the polar cap boundary. The MHD results are slightly better. However, the MHD model also predicts systematically cutoff latitudes that are too low. Future work with the MHD model will include the use of the IGRF model at low L values, instead of the simple dipole field, and a larger number of events for a better statistical analysis.

Appendix A
The calculation of the particle trajectories requires the choice of a suitable numerical method. Since Eqs. (1-3) represent a set of autonomous, yet non-linear ordinary differential equations a variety of methods could be used. One should keep in mind, however, that the particle equations have the nice property that the momentum p, and hence the energy, is con-served. This property allows for a convenient accuracy estimation of the global error which would otherwise be difficult to achieve. Specifically, error control algorithms like extrapolation methods (Bulirsch and Stoer,1991) or Runge-Kutta-Fehlberg methods (Verne, 1978;Burden and Faires, 2001) are based on the local truncation error and thus provide no such global error measures. For that reason methods with error control are of no particular advantage for this problem and their added complexity can be avoided. With some experimentation it can easily be verified that any numerical method with local truncation errors worse than O( t 4 ) is impractical. Implicit methods offer no advantages either because the time step is limited by accuracy requirement and not by stability requirements. One is then essentially left with the explicit Runge-Kutta methods and the Adams-Bashford methods (Kincaid and Cheney, 2002;Burden and Faires, 2001). The latter are more economical because they require only one evaluation of the right hand side per time step, but their convergence order belies their formal order of accuracy because the coefficients of the leading error terms are rather large. There are also other difficulties with Adams-Bashford methods, in particular the integration startup needs to be done with a one-step method, and it is difficult to adapt the time step in the course of the integration. Runge-Kutta methods, on the other hand, are known to be robust and accurate. Because they are single-step methods they do not suffer from the startup and time step adjustment problems of the Adams-Bashford methods. Their primary disadvantage is that they require a rather large number of right hand side evaluations per step. Runge-Kutta methods of very high order have been developed, for example Verner (Verner, 1978) lists Runge-Kutta methods of up to order of 9. Runge-Kutta methods of order higher than 9 become impractical because the coefficients can no longer be represented accurately, in double precision (REAL*8) computer arithmetic. The remaining choice is then the local truncation order of the Runge-Kutta method to be employed. We found that using 6th and 7th order methods indeed offer some advantages in the cases where the right hand side, i.e. B, is given analytically. In those cases, the right hand side is differentiable to any order and the derivatives are bounded. Thus the local truncation errors are also bounded. However, when B is obtained from the gridded MHD solution the situation is entirely different. In this case one needs to distinguish two situations. First, when the entire trajectory of a particle during one time step lies in one and only one cell the right hand side of Eq. (1) is essentially a tri-linear function. Thus, the integration is exact, save for the non-linearity in the equation, with any method that is of at least of order 3 in the local truncation error. Because the time steps are very small most steps fall into this category. On the other hand, when the particle crosses a grid cell boundary during a time step the higher derivatives of the right hand side are no longer bounded. Specifically, B is continuous, the first derivative of B is non-continuous, and there is no bound for the higher derivatives. Thus, the local truncation error could become large. Because the right hand side then violates the usual as-sumptions that are made to estimate the error, namely that the right hand side be Lipshitz-continuous (Kincaid and Cheney, 2002;Burden and Faires, 2001), it is impossible to derive a firm bound on the error. This implies that there can not be an expectation that higher order methods perform better than, for example, the classical 4th order Runge-Kutta scheme. For that reason we use the latter scheme for all the calculations, and defer the error control entirely to the monitoring of the particle momentum.
Error control is then accomplished as follows: We calculate a particle's trajectory with a given time step. If at any time during the calculation the relative error in the momentum exceeds a given threshold we repeat the entire trajectory with a smaller time step. The time step is reduced as many times as necessary to remain within the error bound. We typically use a 5% relative error for this threshold. This may seem large, however, there are only very few particles that come even close to the error threshold. The vast majority of particles finishes the trajectory with a relative error in the momentum of less than a fraction of a percent. Particles with the largest errors are typically those that stay on quasitrapped orbits for a long time. Thus, these particles are likely in or near the penumbra where small disturbances might let the particle go either way. We have repeated some of the calculations with vastly different error bounds, as large as 50% and as small as 0.1%. We found little, if any, difference in the results.