Changes over time in the accuracy of analytical approximations for geomagnetic cutoff rigidity Open Access

Rc is a quantitative estimation of the Earth's geomagnetic field shielding effect (Smart and Shea, 2005; Gao , 2022a, 2022b) since it measures the minimum energy needed by an impinging charged particle of outer space to reach a given location at the Earth. Its value can be obtained through two main approaches. One involves the precise calculation of Rc by tracing the trajectory of the charged particle in Earth's magnetic field. The other relies on energetic considerations, simplifying the Earth's field to a dipolar field, which leads to an analytical expression for Rc.

The first method is the most rigorous and is widely used for cosmic ray research (Shea , 1965; Smart and Shea, 2005, 2009). The calculation of particle trajectories in the Earth's magnetic field has been a subject of interest since its first mathematical formulation by Störmer (Störmer, 1907, 1917, 1930). The fundamental problem is that the trajectory-tracing process involves solving differential equations that do not have solutions in closed form, nor analytical. This was handled using numerical integration methods over many individual trajectories. As the power of computers has improved over time, this procedure has also improved and while the problem is still formidable, thousands of trajectories can be computed without excessive resources (Smart , 2000). Instead of calculating the trajectory of an inward-directed charged particle through the geomagnetic field, this trajectory method is based on tracing an outward-directed particle. The starting point is 20 km above the Earth's surface taking into consideration that the interactions of primary cosmic ray particles with the atmosphere become important below this altitude. The particle trajectory is traced until the particle reaches the outer boundary (defined as 25 Earth radii) or returns to the atmosphere (20 km above the surface) using numerical integration techniques. If the trajectory of a particle reaches the outer boundary, the particle is called an “allowed particle,” which means that a particle at that rigidity can reach the Earth's atmosphere. Otherwise, the particle is called a “forbidden particle,” and a particle at that rigidity cannot penetrate the geomagnetic field. If after a specified integration time neither of these options occur, the particle is considered captured corresponding to a “forbidden particle.” Thus, Rc is obtained tracing the particle trajectories at discrete steps from high rigidity to low rigidity. The process is further complicated by the intricate behavior of cosmic-ray trajectories within the cosmic ray penumbral region, and there is still the issue of the accuracy of the magnetic field models.

The determination of Rc by the trajectory-tracing method in higher order geomagnetic field models has generated reliable values, which are available in world grids from many groups for different epochs (Shea and Smart, 1975, 1983; Smart and Shea, 2008, 2009), and not only the period since 1900, that is the initial year of the International Reference Geomagnetic Field model (IGRF) (Alken , 2021), but also for the geological past and even for hypothetical geomagnetic reversal periods (Stadelmann , 2010; Gao , 2022a, 2022b).

A key aspect in this method is the magnetic field model used in Eq. (1). The limitations on the accuracy of these calculations are the precision of the geomagnetic field model used and the rigidity intervals for which the calculations are performed. The models used are in general a combination of the IGRF and Tsyganenko magnetospheric field model (Tsyganenko, 1995; Tsyganenko and Sitnov, 2007). This means that the greatest uncertainty in the application of charged particle trajectories occurs at low energies, since for low-energy particles, their motion is more easily affected by variations in magnetic fields compared to high-energy particles that have higher momenta.

Considering that the secular variation of the Earth's magnetic field over the last century consists mainly of a decrease in the dipolar component, and the expected transition toward a polarity reversal in some distant future, these analytical Rc expressions are expected to worsen. Even though concerning the reversal, Buffett and Davies (2018), applying a stochastic model based on paleomagnetic observations for the past 2 Myr, obtained less than a 2% chance of a reversed state after 20 kyr, and 11% after 50 kyr. They suggest that a geomagnetic reversal is unlikely within the next 10–20 kyr.

In the present work, the time variation of the agreement between analytical vs trajectory-based Rc estimation is analyzed. It is clear that estimating Rc using an analytical expression is far simpler than the trajectory-tracing method, making it convenient despite errors are introduced. However, how much error are we willing to accept? Or, in other words, how much accuracy in Rc are we willing to sacrifice in exchange for avoiding the computational cost of obtaining its exact value? This reasoning can be extended to the benefit of using simpler models that can shed some light and help to obtain a fundamental understanding of the system being studied, which is often not accessible with complex simulations.

A brief explanation of the analytical Rc calculation and the trajectory-based grids considered are presented in Sec. II. The methodology of the comparative analyses is explained in Sec. III, followed by the results shown in Sec. IV and the discussion and conclusions in Sec. V.

Analytical formulas of Rc depend on the geomagnetic field. The field components needed in the calculation were estimated considering the last IGRF version, IGRF13, and also its Gauss coefficients available at https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html (Alken , 2021).

Since λ in Eq. (6) corresponds to the geomagnetic latitude, in order to obtain a grid that can be directly compared to the trajectory-based Rc grids later in this work, it has to be converted to geographic latitude. This was done considering the centered dipole approximation using the coordinate transformation equations given by Ramana (1998).

This formula is not correct, as mentioned in Comedi (2020), because its derivation assumes critical conditions that apply only to a purely dipolar configuration. This becomes evident in the step-by-step analysis provided by Tsareva (2019) when addressing the same problem of dipolar and quadrupolar axial fields' superposition, and who also does not arrive to an analytical expression for Rc. However, Jory's formula could still serve considering his assumption as a first-order approximation when the quadrupole component is much weaker than the dipole component, as observed under present conditions. In fact, the formula gives consistent results. The inclusion of the axial quadrupolar component is justified by its increase over the past 120 years, in contrast to the decrease in the dipolar component, as depicted in Fig. 1.

The trajectory-based Rc grids estimated by Gvozdevskii (2016) were used with a latitude–longitude resolution of 5° × 15°. These grids are available since 1900–2050 every 5 years, and can be downloaded from http://cosrays.izmiran.ru/dbs/MagEffect/Vertical/CutOff_IGRF/GRID_Table/. It should be noted that these grids were obtained using IGRF12, while our estimations are based on IGRF13. The field components obtained from both IGRF versions are identical from 1900 to 2010. While there are differences in 2015 and 2020, they are less than 0.5% on average in the cases of H, I, the total field intensity, and M.

Figure 2 shows as an example the different analytical approaches and the trajectory-based Rc maps for year 2020. In all the cases, as it is expected from a dominant dipolar field configuration, Rc is higher at the magnetic equator and diminishes to zero toward the magnetic poles.

The factor k_i in Eqs. (11)–(14) corresponds to the cosine of the latitude. It is introduced in order to consider that we have uniformly distributed points in a latitude–longitude grid over a spherical Earth. The reason to use cosine weighting is that, due to the meridian convergence, the width of the grid box (in kilometers) decreases toward higher latitudes proportional to the cosine of the latitude, while the number of observations is the same per degree of latitude.

MAE is a scale-dependent average measure of deviation of analytical from experimental values of Rc. The optimal value is 0 indicating that both series are identical. ME measures the average bias of analytical approximation over or underestimating Rc depending on its sign: positive or negative, respectively. While the optimal value is also 0 implying no-bias, it does not indicate that the series are identical since there can be a compensation. The relative versions of these errors are non-scale dependent, expressed in percentage according to Eqs. (13) and (14). Regarding r², even though the detection of association between spatial datasets is not a simple problem (Clifford , 1989), it measures the explained variance in Rc over the globe grid by each analytical approach. So, it quantifies how similar the spatial variabilities of analytical and experimental Rc values are. It is bounded by 1 and 0, where 1 is the optimal value indicating that the two datasets exhibit identical spatial variability.

Figure 3 shows the spatial distribution of the deviation of each analytical approach from Rc for 1900 and 2020. For both years, it can be noticed an overestimation by Rc1 and Rc2 in the region of the South Atlantic Anomaly (SAA), which becomes stronger with time and shifts following the westward secular displacement of the SAA region (Amit , 2021). This result is expected considering that the magnetic field in these cases do not present this anomaly. In the case of Rc3, the stronger deviation is an underestimation in the region around South Africa, which displaces westward toward South America, that is also almost coincident with the SAA.

Turning now to the time variation of mean errors and correlations, Fig. 4 shows MAE, ME, and r² along 1900–2020. As expected, MAE increases and r drops over time, consistent with the geomagnetic field becoming less dipolar throughout the studied period. In addition, the simplest approach, that is, Rc1, presents always the highest MAE values. Negative ME indicates an overall underestimation of Rc by all the analytical formulas. It remains stable for Rc1 and Rc3 but approaches zero for Rc2. Considering the rise in MAE in all cases, it can be inferred that the compensation between over and underestimations in the analytical formulas has grown over time, especially in the case of Rc2. Here again, the greatest underestimation is that by Rc1. Regarding spatial variability, measured by r², it can be said that there is an overall good agreement always greater than 0.9, although decreasing with time. In this case, the approach with the lowest performance is Rc3.

Figure 5 presents the relative mean errors. To ensure that the elimination of the grid points with Rc = 0 does not impact the errors, it was previously verified that MAE, ME, and r² show negligible variation with respect to Fig. 4.

Contrary to what is observed in terms of error values, in the cases of relative errors, Rc1 shows the best agreement. Another notorious fact is that MRE is mostly positive unlike ME. In order to understand the reason for this, it is useful to consider the zonal distribution of these errors, presented in Figs. 6 and 7, in addition with the decreasing weight of the zonal areas from equator to poles, as can be seen in Fig. 8.

Regarding absolute errors (MAE), Rc3 has extremely high relative values (MRAE) in the vicinity of the south magnetic pole. Even though the area weight is small, the percentage increase—arising from Rc being nearly zero in this region combined with high deviations—significantly exceeds the effect of the reduced area weight explaining the steep rise in Fig. 5(a). For the error case (ME and MRE), the same behavior is observed again in Rc3. Even though the global average of ME is negative, when relative values are estimated (MRE), the high positive deviation at the south magnetic poles compensates for the negative values at middle latitudes, as can be seen in Fig. 5(b). This pronounced disagreement is addressed in Sec. V.

In the cases of Rc1 and Rc2, the high relative positive deviations occur in the vicinity of the north magnetic pole, which is the pole with the greatest shift in position along the period considered. If we add to this what was previously mentioned about the greatest uncertainty in trajectory-based Rc values occurring at low energies, that is around the magnetic poles, we have the additional problem of a decrease in accuracy in Rc itself as an “experimental” measure.

In this study, analytical Rc approximations are compared to values derived from trajectory-based calculations provided by Gvozdevskii (2016). The accuracy of analytical approximations, derived under the assumption of a mainly dipolar magnetic field configuration, has its limitations particularly in the present context with a decline in geomagnetic dipole strength and with an increase in higher-order multipole contributions. Despite these limitations, analytical models remain advantageous due to their computational efficiency and ability to capture essential features of Rc.

Our results show that Rc analytical formulas exhibit increasing errors over time, as reflected in the rising MAE and declining r² global values in all the cases analyzed. This trend aligns with the growing complexity of the geomagnetic field and its departure from a purely dipolar configuration. According to the errors' spatial patterns, the deviations between analytical and trajectory-based Rc values are more pronounced in the SAA region, and in the vicinity of the magnetic poles. For instance, Rc1 and Rc2 consistently overestimate Rc in the SAA, while Rc3 that considers the real horizontal component of the field does not exhibit this issue but underestimates Rc. A reason for this is that the SAA is a region of weak geomagnetic field at the Earth's surface, attributed to reversed flux patches on the core-mantle boundary (Amit , 2021; de Oliveira , 2024), which is not captured if only the dipolar and axial quadrupolar terms are considered in the spherical harmonic decomposition of the magnetic potential. Multipolar terms of higher order have to be included in order to account for the SAA. Thus, there is an overestimation of Rc in the cases of Rc1 and Rc2, as can be noticed in Figs. 2(a) and 2(b) where there is no Rc reduction in the SAA region. In the case of Rc3, which considers IGRF13 to obtain the field variables, it seems to capture the SAA, as can be noticed in Fig. 2(c), with a clear decrease in Rc around this region. However, it exceeds the actual Rc decrease, resulting in negative ME values. The more pronounced decline results from the use of the magnetic field at the Earth's surface in Rc3, where the SAA is most prominent and the field decrease is strongest. At higher altitudes, the SAA weakens, that is, the magnetic field no longer exhibits such a pronounced depression relative to a dipolar configuration. Taking into account that the path of a cosmic ray particle in the Earth's field is determined not only by the field at the surface, but also by the field further out, which in general approaches the dipole field approximation with increasing distance from the Earth, the actual Rc may be expected to be between these values.

Regarding the errors near the magnetic poles, particularly for Rc3, the occurrence of the strongest zonal MRAE and MRE around the south pole has a similar explanation. Charged particles movement are also ruled by the geomagnetic field at higher altitudes where the field is more dipolar, and thus, the iso-contours of Rc are expected to be centered in a position between the dipole and magnetic poles. This can be noticed in Figs. S1 and S2 in the supplementary material. The iso-contours of Rc are centered at the dipolar poles in the cases of Rc1 and Rc2, while Rc3 is centered directly at the magnetic pole taking into account that in this case H and I values of the IGRF13 surface magnetic field are used. This is specially noticed in the south where iso-contours are more circular, as is the case of auroral ovals (Tsyganenko, 2019; Zossi , 2020). For the trajectory-based Rc, the center of the lower value iso-contours lies between the two poles, with subsequent contours aligning more closely with the dipolar pole. Another point adding to this explanation is that the distance between the dipole and the magnetic poles is greatest in the southern case, as can be noticed in Fig. S3, and it has been permanently increasing along the period 1900–2020.

In terms of relative errors, Rc1 performs better, suggesting its suitability for broader, qualitative analyses despite its simplicity. However, the overestimation or underestimation in localized regions, such as the magnetic poles, is amplified when considering relative measures due to Rc values being close to zero in these areas. The zonal analysis confirms a systematic bias across latitude bands, with Rc3 showing higher deviations near the magnetic poles and compensating errors at mid-latitudes.

From the results obtained, the limits of applicability of the analytical approximations analyzed will depend on the required accuracy. If percentage errors up to approximately 30% are acceptable, the approximations become invalid at latitudes higher than ∼65° as we approach the geomagnetic and magnetic poles. If the limits of applicability are also defined in terms of absolute error, then the SAA region should also be outside the range of applicability.

A key aspect in this study is the trade-off between simplicity and accuracy. The errors introduced by the simple analytical Rc formulas are offset by their ability to provide broad insight into Rc trends and their accessibility for applications requiring rapid computations. Held (2005), in the field of climate models, emphasizes the critical role of simple models in scientific understanding, stating that while numerical simulations are invaluable for replicating phenomena, simplified models are essential for isolating and comprehending key processes. Similarly, Vallis (2016) argues that the presence of complex numerical models only increases the need for simple theoretical frameworks to elucidate basic principles. In the case of Rc, analytical approximations enable researchers to identify features, which might otherwise be obscured by the intricacies of high-resolution numerical models. For example, in Comedi (2020) using Rc1 approximation, the spatial pattern of Rc secular trends was interpreted based on the inclination of the dipole axis, which is becoming more aligned with the Earth's rotation axis.

The simplicity of analytical models ensures their relevance for applications where computational resources are limited or where a rapid understanding of Rc behavior is required. However, as the geomagnetic field continues to evolve, assuming that the axial dipolar component will continue weakening, the discrepancies between analytical and numerical approaches are expected to grow. Future work should focus on refining analytical expressions to incorporate higher-order multipoles while retaining computational efficiency. What theoretical approach might work in geomagnetic reversal scenarios? Perhaps the axial multipolar configurations analyzed analytically by Urban (1965) and Shebalin (2004) could offer some insight. Alternatively, efforts to find solutions for the superposition of axial fields proposed by Tsareva (2019) may be useful, though they remain unresolved analytically or in a simplified form. This remains a significant challenge. How long will these theoretical approximations remain valid? Currently, they already exhibit global absolute errors on the order of 20%, even in the best cases. However, on geological timescales in the future, when the field eventually returns to a predominantly dipolar state and with a dominance higher than the current of ∼80%, these formulas will regain their usefulness.

See the supplementary material for additional figures to analyze the discrepancies between Rc analytical formulas and trajectory based values in the regions around the dipolar and magnetic poles.

Authors acknowledge research projects PIP 2957 and PIUNT E756.

The authors have no conflicts to disclose.

Ana G. Elias: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Software (equal); Supervision (equal); Writing – original draft (equal). Maria A. Vega Caro: Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Facundo Abaca: Investigation (equal); Methodology (equal); Writing – review & editing (equal). Bruno S. Zossi: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Software (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available: IGRF coefficients at https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html, Alken . (2021); trajectory-based Rc grids at http://cosrays.izmiran.ru/dbs/MagEffect/Vertical/CutOff_IGRF/GRID_Table/;online calculator of cutoff rigidities at https://tools.izmiran.ru/cutoff/; and Rc dataset obtained from analytical approximations (Rc1, Rc2, and Rc3) at 10.5281/zenodo.15287681, Elias (2025).