© J. Gao et al., Published by EDP Sciences 2022

1 Introduction

The geomagnetic field acts as a natural shield against charged energetic particles, such as galactic cosmic rays, that impact the Earth’s atmosphere. The primary cosmic ray particles are protons that are deflected when moving in the geomagnetic field (Beer et al., 2012; Bütikofer, 2018). Since energetic particles with equal rigidity R = pc/Ze (momentum per unit charge of the particle) have the same moving trajectory in a constant magnetic field, rigidity is frequently used by the cosmic ray community to study geomagnetic field shielding. Thus, the geomagnetic field shielding effect can be quantitatively estimated as the geomagnetic cutoff rigidity (Smart et al., 2000; Smart & Shea, 2009). When the time is fixed, at a given location and incidence angle, the geomagnetic cutoff rigidity R_c is a threshold value, and only particles with rigidity higher than this threshold value can reach the Earth’s atmosphere. The effective cutoff rigidity is an average of the shielding effect, considering incident particles in all incident angles, while the effective vertical cutoff rigidity, R_vc, only considers incident particles in the vertical direction (Cooke et al., 1991). The calculation of the cutoff rigidity started with the pioneering work by Størmer (1955) through theoretical analysis, which only considered the dipole component of the geomagnetic field. The cutoff rigidity is determined more accurately by tracing the trajectory of the particles using a geomagnetic field model and considering both the dipole and non-dipole components (Shea et al., 1965). Since the geomagnetic field changed constantly in the geological past, the geomagnetic shielding was also constantly changing. The long-term modulation of the galactic cosmic rays by the geomagnetic field is referred to as “geomagnetic modulation” in the cosmic ray community (e.g., Beer et al., 2012).

Understanding the “geomagnetic modulation” is fundamental to studying the long-term variation of the cosmic ray flux, which controls the cosmogenic radionuclides (e.g., ¹⁰Be, ¹⁴C, ³⁶Cl, ³⁹Ar, half-lives in the range of ~300–1.5 × 10⁶ years) production rate, cosmic ray induced ionization rate, and cosmic ray radiation intensity (Shea & Smart, 2000; Lifton et al., 2008; Usoskin et al., 2008, 2010). When using the cosmogenic radionuclide production rate to reconstruct the solar activity or for cosmogenic nuclides dating, the long-term trend caused by “geomagnetic modulation” must be removed (Desilets & Zreda, 2003; Solanki et al., 2004; Usoskin, 2017). Cutoff rigidity, as a quantitative estimation of the “geomagnetic modulation”, can be inferred from geomagnetic field models. For the current geomagnetic field, the cutoff rigidity has been extensively studied using the standard International Geomagnetic Reference Field (IGRF) model (Shea & Smart, 1990; Gerontidou et al., 2021); however, the variation of cutoff rigidity on multi-millennial timescales is much less studied. Most of the earlier studies which inferred the long-term variation of the cutoff rigidity treated the geomagnetic field as an ideal dipole field (see e.g., Lal, 1991; Gosse & Phillips, 2001; Usoskin et al., 2006). Lifton (2016) calculated the cutoff rigidity variation covering the Holocene (14–0 ka, with “ka” standing for thousand years before present, where “present” is 1950 AD as used in radiocarbon dating) derived from global paleomagnetic field models that consider non-dipole field contributions: CALS10k.1 and SHA.DIF.14k (Korte et al., 2011; Pavón-Carrasco et al., 2014). Due to the previous lack of paleomagnetic field models before 14 ka, further back in time, Lifton (2016) treated the geomagnetic field as a geocentric axial dipole (GAD) field and estimated the cutoff rigidity using virtual axial dipole moment (VADM) records that are obtained from global averages of regional field intensity records (see e.g., review by Panovska et al., 2019 for more details about the VADM). Lifton (2016) found clear regional differences when using paleomagnetic field models or a simple GAD to estimate the cutoff rigidity.

Knowledge about the spatial and temporal variations of the geomagnetic field is essential to reconstruct the long-term geomagnetic shielding effect. In the last 100 ka, global or regional geomagnetic excursions occurred several times (e.g., Singer, 2014; Laj & Channell, 2015; Panovska et al., 2019). The Laschamps excursion is the best-documented excursion that occurred at around 41 ka, when the geomagnetic field intensity dramatically decreased and the non-dipole component of the geomagnetic field strongly influenced the geomagnetic field morphology (Korte et al., 2019b). In a previous study, we showed that the influence of the non-dipole field could not be ignored when studying the geomagnetic shielding effect during this excursion (Gao et al., 2022). The Mono Lake/Auckland excursion (~34 ka) and the Norwegian–Greenland Sea excursion (~65 ka) also show demonstrably dipole field decrease along with significant directional deviations from the axial dipole field. Information about the long-term variation of the geomagnetic field comes from paleomagnetic sediment records, archaeomagnetic, and lava flow records (e.g., Brown et al., 2015a, 2015b; Constable & Korte, 2015; Panovska et al., 2018a). The knowledge about the paleomagnetic field evolution has constantly improved from individual records to regional or global stacks of relative paleofield intensity (RPI) (e.g., Guyodo & Valet, 1999; Laj et al., 2004) or equivalent virtual (axial) dipole moments (VADM or VDM) (e.g. Ziegler et al., 2011; Ziegler & Constable, 2015), and to global continuous paleomagnetic field models that resolve large-scale non-dipole field contributions (e.g., Leonhardt et al., 2009; Constable et al., 2016; Korte et al., 2019a; Panovska et al., 2021).

In recent years, several paleomagnetic field models have been established to reveal the long-term variations of the geomagnetic field covering up to the last 100 ka (see e.g., review by Panovska et al., 2019, and references therein). The available models focus on different timescales. The GGF100k model, for example, concentrates on the long time variation of the geomagnetic field spanning the past 100 ka and has similar spatial resolution but less temporal resolution compared with Holocene models (Panovska et al., 2018b). On the contrary, the LSMOD.2 model provides a more detailed view of the Laschamps excursion but only covers a limited time interval (50–30 ka) (Korte et al., 2019b). Most recently, Panovska et al. (2021) presented a new model called GGFSS70, which spans 70–15 ka with higher temporal resolution than GGF100k and indicates the characteristics of three global excursions: the Norwegian–Greenland Sea, Laschamps, and Mono Lake excursion. Those paleomagnetic field models provide the possibility to reconstruct the long-term variation of the cutoff rigidity considering for the first time also the non-dipole field contributions.

In this study, we reconstruct the geomagnetic cutoff rigidity variation over the last 100 ka using four global paleomagnetic field models, i.e., the GGF100k, GGFSS70, LSMOD.2, and CALS10k.2. They are introduced in Section 2. Global grids of the cutoff rigidities are presented in 200-year steps for the period 100–10 ka and 50-year steps for the period 10–0 ka, using a trajectory tracing program (see Supplementary Movie S1). The detailed differences in cutoff rigidity results using different paleomagnetic field models are investigated in Section 3. The models which are considered to best resolve the cutoff rigidity structure during different time intervals are selected for calculating a combined record of the effective vertical cutoff rigidity, named Rc100k. In Section 4, we present the variation of latitudinal cutoff rigidity structure, impact area, and global cosmic ray flux over the last 100 ka based on Rc100k. Discussion and conclusion are given in Section 5.

2 Paleomagnetic field models

Over the last decades, a growing number of long-term global continuous paleomagnetic field models have been published. Those models are based on the compilation of paleomagnetic data and provide a global view of geomagnetic field evolution within the model time interval. In this section, we briefly introduce the modeling method and compare the spatial and temporal resolution of five paleomagnetic models utilized in this study, namely IMOLEe, CALS10k.2, GGF100k, LSMOD.2, and GGFSS70.

Spherical harmonic (SH) functions are commonly used to obtain a global picture of geomagnetic field evolution. In a source-free region, the time-depending internal magnetic field B can be expressed as the negative gradient of a potential V:

$$V\left(r,\theta ,\phi ,t\right)=a\sum _{l=1}^{L_{\max }}\mathrm{}\sum _{m=0}^l\mathrm{}{\left(\frac{a}{r}\right)}^{l+1}\left(g_l^m\left(t\right)\cos \left(\mathrm{m\phi }\right)+h_l^m\left(t\right)\sin \left(\mathrm{m\phi }\right)\right)\times P_l^m\left(\cos \theta \right),$$(1)

where r, θ, ϕ are the geocentric spherical coordinates, and a is the Earth’s mean radius, 6371.2 km. L_max is the maximum SH degree, and $P_l^m$ is the Schmidt quasi-normalized associated Legendre functions. Based on archaeomagnetic, volcanic, and paleomagnetic sediment data, paleomagnetic field models utilize various inversion techniques to solve for the time-varying spherical harmonic coefficients ($g_l^m\left(t\right)$, $h_l^m\left(t\right)$). A major advantage of the SH functions is that dipole and non-dipole components of the geomagnetic field can be easily separated (see e.g., review by Panovska et al., 2019 for more details).

Here, we use five long-term continuous SH paleomagnetic field models to obtain the best currently available view of geomagnetic field evolution over the last 100 ka. The IMOLEe model (Leonhardt et al., 2009), spanning 44–36 ka, was the first continuous global paleomagnetic field model covering the Laschamps excursion. The IMOLEe model is constrained by only six palaeomagnetic records, mostly located in the North Atlantic and the Gulf of Mexico. Being the oldest of the models, it might now be the least reliable given the sparse data distribution compared to more recently developed models. However, to get a good idea about differences in cutoff rigidity depending on differences in underlying models, we include results based on this model in the comparison in Section 3.

The CALS10k.2 model (Constable et al., 2016), spanning 10–0 ka, is the latest version of the CALS model series model (Korte et al., 2009, 2011; Korte & Muscheler, 2012) and supersedes all previous CALS series models (Constable et al., 2016). For the historical period of the past 400 years, CALS10k.2 is constrained by and agrees closely with the gufm1 model (Jackson et al., 2000), with the spatial and temporal resolution strongly increased compared to the earlier millennia. As our focus here lies on longer time scales, we do not consider any of the numerous other recently developed Holocene field models (see e.g., Constable & Korte, 2015 for a review). The large-scale surface features of the geomagnetic field in the last 10 ka reconstructed by the CALS10k.2 model are consistent with other Holocene models, e.g., HFM.OL1.A1 and ARCH10k.1 (Constable et al., 2016).

The GGF100k model (Panovska et al., 2018b) provides information about geomagnetic field evolution covering the full past 100 ka. It is constrained by more than 100 continuous sediment records and utilizes temporal smoothing kernels to fit both high and low-resolution records. Nevertheless, it has a lower temporal resolution than all the shorter models. It is the only available global model for some intervals within the past 100 ka.

The LSMOD.2 model (Korte et al., 2019b), spanning 50–30 ka, provides a more detailed view of the Laschamps (~41 ka) and Mono Lake/Auckland (~34 ka) excursions. The LSMOD.2 model is reconstructed from volcanic and high-resolution sediment data and the paleomagnetic data included in the LSMOD.2 model is notably increased compared to the IMOLEe model. A series of models based on different data sets were developed to test the influence of data selection. The preferred model shows similar excursion characteristics compared with other models and provides robust characteristics of the Laschamps and Mono Lake/Auckland excursions.

The GGFSS70 model (Panovska et al., 2021), spanning 70–15 ka, is a model built from nine globally distributed sediment records, considered to have good chronological information with the most up-to-date age calibrations. The GGFSS70 model provides a clear view of three excursions: the Norwegian–Greenland Sea (~65 ka), Laschamps (~41 ka), and Mono Lake/Auckland (~34 ka). It is the first continuous global paleomagnetic field model representing the Norwegian–Greenland Sea excursion in detail.

Figures 1a and 1b show the dipole power and non-dipole power of the five paleomagnetic field models evaluated at Earth’s surface. Apart from the clearly global Laschamps excursion (La, ~41 ka) with the lowest dipole power, two further excursions are manifested by clear lows in dipole power, in particular in GGFSS70, i.e, the Norwegian–Greenland Sea (NGS, ~65 ka), and Mono Lake/Auckland (ML, ~34 ka) excursions. Further excursions were observed only regionally, probably due to a less strong dipole decrease, e.g., the Post Blake (PB, ~98 ka) (Thouveny et al., 2004) and Hilina Pali (HP, ~19 ka) (Ahn et al., 2018) excursions. Those only show a moderate decrease in dipole power. Since the SH models are smoothed representations of the paleomagnetic data, some details of the excursions may not be well represented by the global SH models, in particular by the most strongly smoothed GGF100k (Panovska et al., 2019, 2021).

Fig. 1 Dipole power (a) and non-dipole power (b) evaluated at Earth’s surface over the last 100 ka, as predicted by five continuous spherical harmonic geomagnetic field models: IMOLEe (Leonhardt et al., 2009), CALS10k.2 (Constable et al., 2016), GGF100k (Panovska et al., 2018b), LSMOD.2 (Korte et al., 2019b), and GGFSS70 (Panovska et al., 2021). Five global and regional excursions are labeled on the top panel a: PB: Post Blake, NGS: Norwegian–Greenland Sea, La: Laschamps, ML: Mono Lake/Auckland, and HP: Hilina Pali. Gray areas indicate the three excursions studied in more detail here: NGS, La, and ML. The time intervals of the model, and models used to calculate Rc100k are given by the color bars at the bottom.

As described above, the SH models cover various time intervals and have different spatial and temporal resolutions. In order to get the best currently available view of the geomagnetic field variation over the last 100 ka, it is necessary to combine the models. To get an idea of their similarities and differences, Supplementary Movies S2 and S3 show the radial component B_r and geomagnetic field intensity B_total of the paleomagnetic models over the entire time interval (100–0 ka). We will discuss the influence of the differences for estimating cutoff rigidities in Section 3.

3 Cutoff rigidity during excursions from different paleomagnetic models

Compared with theoretical equations, the trajectory tracing method provides a more accurate estimation of the geomagnetic cutoff rigidity (Smart & Shea, 2009). In this section, we briefly introduce the methodology of trajectory tracing and compare the cutoff rigidity results using different geomagnetic field models during the Norwegian–Greenland Sea, Laschamps, and Mono Lake/Auckland excursions. More details about the trajectory tracing program can be found in Smart et al. (2000).

The common method to determine the trajectory of a cosmic ray particle is by tracing a negatively charged particle, as the trajectory of an outward-directed negatively charged particle is identical to the trajectory of an inward-directed positively charged particle with the same charge and the same rigidity. We traced the particle starting at 20 km above the Earth’s surface (geodetic coordinates are used) as an approximation for the top of the atmosphere. The trajectory tracing is terminated when the particle reaches the outer boundary (taken as 25 earth radii away from the Earth’s center) or returns to the atmosphere. The influence of the magnetosphere current on geomagnetic cutoff rigidity calculation is not considered in this study. If the trajectory of the particle reaches the outer boundary, it is called an “allowed particle”, otherwise, it is a “forbidden particle”, respectively. Note that the cutoff rigidity is not sharp. In the transition from the first forbidden rigidity down to the last allowed rigidity, there is a penumbra region, in which both allowed and forbidden rigidities exist. At each given location, after scanning particle rigidity at discrete steps from low to high through the energetic particles spectrum, the effective vertical cutoff rigidity R_vc is calculated, which is a weighted average over the rigidity penumbra and describes the transparency of the penumbra (Cooke et al., 1991; Bütikofer, 2018). The R_vc (in units of GV, i.e. 10⁹ Volt), tracing particles initiated in the vertical direction, is the most widely used approximation for the effective cutoff rigidity (Smart et al., 2000). Energetic particles with different incidence directions can also be considered in cutoff rigidity calculation, and other approximations for the effective cutoff rigidity, e.g., the apparent cutoff rigidity can also be determined (apparent cutoff rigidity is an average over cutoff rigidity in nine different directions Clem et al., 1997). However, the difference between the vertical and apparent cutoff rigidity is about 1 GV at the equator region for today’s geomagnetic field and is negligible globally during the Laschamps excursion (Gao et al., 2022). Details of the R_vc calculation can be found in Shea et al. (1965); Bütikofer (2018) and Gao et al. (2022).

In this study, R_vc was calculated on global grids of 5° in latitude by 5° in longitude, using paleomagnetic field models spanning up to the last 100 ka in 200-year steps for the GGF100k, GGFSS70, LSMOD.2, and IMOLEe models, and in 50-year steps for the CALS10k.2 model. Calculating the high-resolution cutoff rigidity is time-consuming even for modern computers. To reduce the computation time, the rigidity was determined in two steps, refining the scan over the penumbra region (Gao et al., 2022). The first rigidity scan is from 30 GV to 0 GV with an 0.5 GV rigidity interval (covering the spectrum of the energetic particles Mironova et al., 2015). The second rigidity scan is from the highest forbidden rigidity value found in the first step plus 0.5 GV to the lowest allowed rigidity value minus 0.5 GV using a 0.01 GV rigidity interval. This hybrid rigidity scan can significantly decrease the computation time without reducing the accuracy (Gao et al., 2022). Also, the influence of the crustal field on cutoff rigidity determination is ignored because the crustal field is insignificant in calculating the cutoff rigidity even during geomagnetic field transitions (Gao et al., 2022).

Paleomagnetic field models have uncertainties from limitations in data distribution, data reliability, and dating accuracy, but these uncertainties are not well quantified. However, differences among their predictions, e.g., among geomagnetic field evolution maps give indications of these uncertainties. We expect the effect of the uncertainties on cutoff rigidities to be strongest during the excursions when the dipole is weak. Therefore, we compare the cutoff rigidity results from different models during the three most prominent excursions in the time interval: the Laschamps, Norwegian–Greenland Sea, and Mono Lake/Auckland excursions. We then selected the most suitable model to reconstruct the cutoff rigidity variation for different time intervals, and we present here the cutoff rigidity results covering the last 100 ka, named Rc100k (see Supplementary Movie S1). The time series of paleomagnetic field model selection for calculation Rc100k can be found at the bottom of Figures 1 and 5, and the transition times at which we switched from one model to another for the final record are given in Table 1.

3.1 Laschamps excursion

The Laschamps excursion, first discovered in the Laschamps lava flows, is the best-documented excursion (see e.g., review by Laj & Channell, 2015) and is documented by records all around the world. Recent paleomagnetic field models demonstrated that the dipole power is close to zero, and the axial dipole coefficient $g_1^0$ is not only reduced to zero but also reverses its sign for a few hundred years in the middle of the excursion (Korte et al., 2019b; Panovska et al., 2021). Based on the LSMOD.2 and GGFSS70 model, the dipole moment reached its minimum of 0.36 × 10²² Am² and 0.12 × 10²² Am² at 40.98 ka and 40.97 ka, respectively, which is clearly lower than the present-day value of about 8 × 10²² Am².

There is no clear definition for the start and end times of geomagnetic excursion events. Several parameters have been proposed to indicate the excursion duration time, one of them is the paleosecular variation (PSV) index (Panovska & Constable, 2017), which considers both field directions and intensity to quantitatively estimate field complexity and variation activity and provides a threshold value for transitional field as found during excursions and reversals. According to the PSV index, the Laschamps excursion started at 41.9 ka and ended at 40.1 ka, with a duration of 1.8 kyr (Korte et al., 2019b). A detailed presentation of cutoff rigidity maps during the Laschamps excursion based on the LSMOD.2 model is given by Gao et al. (2022). We now compared four SH models covering the Laschamps excursion, i.e., GGF100k, GGFSS70, LSMOD.2, and IMOLEe model.

Figure 2 illustrates the different cutoff rigidity results covering the Laschamps excursion by four snapshots in time. The geomagnetic field variations and the cutoff rigidity changes over the full time are provided in Supplementary Movies S1 and S2, respectively. All models show a clear dipole field decrease during the excursion (Fig. 1a), leading to less latitudinal dependence of the cutoff rigidities, which are globally lower than 4 GV at around 40.95 ka. The GGF100k model has a higher dipole moment at the middle of the excursion due to its limited temporal resolution, and thus is least well suited to provide cutoff rigidities for this time. The IMOLEe model gives generally similar results as GGFSS70 and LSMOD.2, but also with some clear differences. For example, although the general latitudinal distribution of cutoff rigidities is similar in amplitudes at 42.45 ka, there are clear longitudinal differences. This is due to a local maximum of the geomagnetic field over Africa inferred from the IMOLEe model (see Movie S3), while the GGFSS70 and LSMOD.2 models indicated a local maximum over the Southern Pacific Ocean instead. We attribute this difference to the sparse data distribution of the IMOLEe model, and consider the GGFSS70 and LSMOD.2 models that are constructed by using more globally distributed records as more reliable. The global geomagnetic field patterns and resulting cutoff rigidities from the GGFSS70 and LSMOD.2 models are generally in good agreement. As the LSMOD.2 model has a slightly higher temporal resolution than the GGFSS70 model, the LSMOD.2 model is selected to reconstruct the long-term cutoff rigidity Rc100k during and around the Laschamps excursion.

Fig. 2 Effective vertical cutoff rigidity contour maps covering the period of Laschamps excursion (43.95–39.45 ka) calculated from the (a) GGF100k, (b) GGFSS70, (c) LSMOD.2, and (d) IMOLEe models.

3.2 Norwegian–Greenland Sea excursion

The Norwegian–Greenland Sea (NGS) excursion was first reported using Norwegian Greenland Sea sediments (Bleil & Gard, 1989; Nowaczyk & Frederichs, 1999). Although it has been documented in much fewer locations than the Laschamps, the two models covering the NGS excursion (GGF100k and GGFSS70) indicate a clear drop in dipole field strength. Using the PSV index and the GGFSS70 model, Panovska et al. (2021) indicated that the transitional epoch is over 66–56 ka, and two PSV peaks appeared at 65–63 ka and 60–58 ka. As shown in Figure 1a, there are also two clear dipole power drops at ~65 ka and ~59 ka.

Figures 3a and 3b show four snapshots of the cutoff rigidity results during the NGS excursion (66–63 ka) using the GGF100k and GGFSS70 model, respectively. The snapshots cover the first peak of the PSV value, which corresponds to the dipole power local minimum at ~65 ka. The GGFSS70 model shows a clear dipole field decrease at 65–64 ka, the lowest dipole moment is 1.87 × 10²² Am² at 64.9 ka (approximately 23% of today’s dipole moment, 8 × 10²² Am²). The GGF100k model indicates only a slight decrease in the dipole field around 65 ka. Although the dipole power of the GGF100k model is higher than that of the GGFSS70 at 63 ka, the dipole tilt angles from both models are similar and result in similar global geomagnetic field and cutoff rigidity patterns. Since the GGF100k model has a lower temporal resolution compared with the GGFSS70 model, the NGS excursion features are less well resolved by the GGF100k model. We in general, prefer the higher resolution GGFSS70 model over GGF100k and use it to calculate the Rc100k over and around the NGS excursion.

Fig. 3 Effective vertical cutoff rigidity contour maps covering the period of the Norwegian–Greenland Sea excursion (66–63 ka) calculated from the (a) GGF100k and (b) GGFSS70 models.

3.3 Mono Lake/Auckland excursion

The excursion previously and widely known as Mono Lake (ML) excursion follows soon after the Laschamps excursion at around 34 ka (e.g., Benson et al., 2003). This event has been termed the Mono Lake or Auckland excursion, depending on the locality where the excursion was first identified, and there is a long-standing discussion and new recent evidence that the event found at the location of Mono Lake had large dating errors and might be the Laschamps excursion (see Marcaida et al., 2019). Although there are strong suggestions to prefer the name Auckland excursion, we use Mono Lake/Auckland here to avoid confusion. Korte et al. (2019b) argued that the ML excursion is not a single excursion but rather a series of excursions that occurred between 36 and 30 ka. Using the GGFSS70 model, Panovska et al. (2021) suggested that the ML excursion is a double event, which has two dipole field intensity drops at 34 and 29 ka. The minimum dipole moment during the ML excursion according to this model is 3.3 × 10²² Am² at 34.4 ka, approximately 41% of today’s dipole moment.

The cutoff rigidity results during the first dipole field drop of the ML excursion calculated from the GGF100k, GGFSS70, and LSMOD.2 models are shown in Figure 4. Again, the dipole moments of both the GGFSS70 and LSMOD.2 models are lower than that of the GGF100k model during the ML excursion, and the ML excursion seems suppressed by strong temporal smoothing in the GGF100k model. The dipole field variation during the ML excursion appears consistent (see Movie S3) among the GGFSS70 and LSMOD.2 models; however, the small-scale structure of those two geomagnetic field models show regional differences. At 33.45 ka, for example, two local cutoff rigidity maxima are found over South America and the Western Pacific from the GGFSS70 model, while the intensity maximum is found over central Africa from the LSMOD.2 model. Over central Africa, the GGFSS70 model results in cutoff rigidity close to 9 GV, while the LSMOD.2 model results in 13 GV, giving differences in cutoff rigidity over central Africa up to 4 GV during the ML excursion. Since the LSMOD.2 model has a slightly higher temporal resolution; we use this model for its full-time interval, including the ML excursion. However, the dipole field drop around 29 ka is not covered by the time interval of LSMOD.2, and we then use the GGFSS70 model with higher spatial and temporal resolution than the GGF100k model.

Fig. 4 Effective vertical cutoff rigidity contour maps covering the period of the Mono Lake/Auckland excursion (36.45–33.45 ka) calculated from the (a) GGF100k, (b) GGFSS70, and (c) LSMOD.2 models.

3.4 Combining the paleomagnetic field models

For the three most prominent excursions, the most suitable model is selected to reconstruct the cutoff rigidity variation described above. The CALS10k.2 model is used to study cutoff rigidity evolution during the Holocene (10–0 ka). Furthermore, within the time range 100–70 ka and 18–10 ka, the GGF100k model is selected to calculate the Rc100k because no other paleomagnetic field model covers these periods. However, to obtain smooth transitions from one model to the other, we chose transition times when the dipole moment of the adjoining models agree closely. The chosen transition times are given in Table 1, as well as the dipole power, non-dipole power, and the dipole tilt angle of the different models at those epochs. The geomagnetic field is dipole dominated at all the transition times, and the geomagnetic field patterns from the different models are similar (see Figs. S1 and S2). The cutoff rigidity differences between different models are generally small but exceed 6 GV (see Fig. S3) at a few locations and some of the transition times.

Figure 5 shows the axial dipole moment (ADM) over the past 100,000 years obtained from the combined paleomagnetic field models (GGF100k, GGFSS70, LSMOD.2, and CALS10k.2) (Panovska et al., 2019). The comparison of the axial dipole moment from the combined model over the last 100 ka to several virtual axial dipole moment (VADM) stacks (GLOPIS-75 (Laj et al., 2004), C2018-Overall stack (Channell et al., 2018), sint-800 (Guyodo & Valet, 1999), and sint-2000 (Valet et al., 2005)) overall gives a reasonable agreement (see Fig. 5). However, the NGS excursion appears longer in our model, with a second ADM low at ~59 ka (from the GGFSS70 model), which is not found in any of the VADM stacks, and the ADM is in general somewhat lower between 59 and 50 ka. The high values of the C2018-Overall stack around 30 to 12 ka, not found in any of the other curves, likely result from regional bias due to the dominance of Atlantic records in this stack. At that time, our model agrees well with other stacks.

Fig. 5 Axial dipole moment (ADM) evolution with time as combination from four paleomagnetic field models: GGF100k (100–67.85 ka and 17.55–9.95 ka), GGFSS70 (67.85–49.95 ka and 32.66–17.55 ka), LSMOD.2 (49.95–32.65 ka), and CALS10k.2 (9.95–0 ka). Virtual axial dipole moment (VADM) stacks obtained from globally distributed relative paleointensity (RPI) records are also shown for comparison: GLOPIS-75 (Laj et al., 2004), C2018-Overall stack (Channell et al., 2018), sint-800 (Guyodo & Valet, 1999), and sint-2000 (Channell et al., 2018).

To quantitatively assess the influence of the model selection on cutoff rigidity estimation, the global average of the cutoff rigidity differences obtained from different models are shown in Supplemental Figure S4 (see Gao et al., 2022, Eq. (11)). For comparison, global averages of the value differences between the Rc100k and the cutoff rigidity calculated using the model with limited maximum SH degree are also shown (see Supplemental Fig. S5). Differential maps of R_vc over the full time are provided in Movies S4 and S5. Due to the systematic difference in dipole moment level between GGF100k and GGFSS70 (see Fig. 1), the globally averaged cutoff rigidity difference between these two models is rather large, and similarly, the difference between IMOLEe and LSMOD.2 during the Laschamps excursion is large due to the clear difference in non-dipole power (see also Fig. 1) at a time when the field is not dipole dominated. Otherwise, the differences among cutoff rigidity obtained from different models are mostly close to the difference between the Rc100k and the estimates based on GAD. Apart from the influence of the non-dipole field, the model selection is also crucial to constrain the cutoff rigidity variation over the last 100 ka.

3.5 Uncertainty estimates

The paleomagnetic field data clearly contain uncertainties both in magnetic data and their ages, but it is difficult to get reasonable uncertainty estimates for the paleomagnetic field models, as standard methods tend to underestimate the model uncertainties for regions where data are sparse or non-existent. Here, we used a bootstrap technique to obtain the formal uncertainties of the GGF100k model and propagated these to cutoff rigidity uncertainties. The bootstrap technique to obtain model uncertainties has previously been applied to Holocene geomagnetic field models (Korte & Constable, 2008; Korte et al., 2011) (see Supplementary Text S1 for more details about the bootstrap technique applied to the GGF100k model).

In general, paleomagnetic field uncertainties increase with time, i.e., smaller at the recent end and larger close to 100 ka. A direct error propagation from the magnetic field components to the cutoff rigidity are not straightforward. Therefore, we estimated the upper and lower bound for the cutoff rigidity for two epochs by considering the minimum and maximum magnetic field values at selected epochs. The cutoff rigidity uncertainties are up to 2 GV at the equator when the geomagnetic field is dipole dominated. During the Laschamps excursion, the geomagnetic field is less dipole dominated, and the uncertainty of the cutoff rigidity is also smaller, about 1.5 GV at the equator. The cutoff rigidity differences between the different models (e.g., the LSMOD.2 and the GGFSS70 model) are larger than the cutoff rigidity uncertainties, and the GGF100k uncertainty estimates probably underestimate the true model uncertainties. The actual cutoff rigidity uncertainties are likely in the order of the cutoff rigidity differences between using the two high-resolution models, LSMOD.2 and GGFSS70, which are based on somewhat different data sets.

4 Results

We reconstructed the cutoff rigidity variation over the last 100 ka from the combined model using trajectory tracing. Global grids (5° in latitude, 5° in longitude) of the effective vertical cutoff rigidities R_vc are presented covering the last 100 ka, and we name this product Rc100k (see Movie S1). Based on the Rc100k, the latitudinal cutoff rigidity profile, impact area, and the global cosmic ray flux during the last 100 ka are presented in this section.

4.1 Latitudinal cutoff rigidity structure

With today’s strongly dipole-dominated geomagnetic field, the cutoff rigidity is inversely correlated with the geomagnetic latitude. The cosmic ray intensity observed on the ground increases with geomagnetic latitude from the equator to about 45° geomagnetic latitude because the geomagnetic shielding is weaker than the atmosphere shielding above this latitude. This latitude is called the “knee” of the cosmic ray latitude curve, at about 45° geomagnetic latitude. The cutoff rigidity value at the latitude “knee” is dependent on both the calculation method and the geomagnetic field model, and is ~6.3 GV using equations considering only the dipole field, or ~4.0 GV using a trajectory tracing program that considered both dipole and non-dipole field (Shea & Smart, 1990; Smart & Shea, 2009).

However, the inverse correlation is no more valid during excursions, when the geomagnetic field is no more dipole dominated. For example, at the midpoint of the Laschamps excursion, the cutoff rigidities are globally lower than 4 GV and nearly independent of latitude (Fig. 2), and the latitudinal dependence also largely breaks down during the NGS excursion (Fig. 3b). During geomagnetic polarity transitions, there is no clear definition of the “geomagnetic latitude”. We show the mean value of R_vc averaged over longitude at geographic latitudes 45° N, 0° N, and 45° S in Figure 6a. Furthermore, the latitudinal and longitudinal structure of the cutoff rigidity is provided in Movie S1g. The latitudinal (longitudinal) structure of the cutoff rigidity is an average of cutoff rigidity over longitude (latitude).

Fig. 6 (a) Mean value of Rvc over longitude at latitudes 45° N, 0° N, and 45° S; (b) impact area of energetic particles at Rvc equal to 2.5 GV, 5 GV, and 7.5 GV; (c) impact area of energetic particles at Rvc equal to 2.5 GV and 5 GV, with solid lines (dashed lines) indicating the impact area in the northern hemisphere (southern hemisphere), respectively; (d) global cosmic ray proton flux accessing the Earth’s atmosphere; (e) cosmic ray proton flux reaching the Earth’s atmosphere due to its rigidities during the last 100 ka. For comparison, today’s values calculated using the International Geomagnetic Reference Field (IGRF) at epoch 2020 (Alken et al., 2021) are included in panels a and b as dashed lines.

For today’s geomagnetic field, the R_vc is ~15 GV at the equator and ~5 GV at latitude ±45°. During the excursions, as the geomagnetic field strength decreased, the global cutoff rigidity is accordingly decreased. The R_vc at the equator regionally reduces to ~3 GV during the NGS and the Laschamps excursions and reduces to ~7 GV during the ML excursion. The mean value of R_vc at Earth’s equator versus the dipole moment of the geomagnetic field is provided in Figure S6. The R_vc at the equator depends nearly linearly on the dipole moment of the geomagnetic field, but this relation changes when the dipole moment gets very low, around 2 × 10²² Am², and the non-dipole field becomes important. A hemispheric asymmetry of R_vc is observed during the NGS excursion. The R_vc at latitude 45° N is decreased to about 1 GV during the first dipole field drop (65–63 ka), and the R_vc at latitude 45° S is decreased to about 1.5 GV during the second dipole field drop (60–58 ka).

4.2 Impact area

Compared to the latitudinal cutoff rigidity structure, the impact area is more convenient for describing the global cutoff rigidity, especially during excursions when the inverse correlation between the cutoff rigidity and geomagnetic latitude no longer exists (Vogt et al., 2007; Stadelmann et al., 2010). The impact area, expressed as a percentage, is defined as the area where the energetic particles at a given cutoff rigidity can reach the atmosphere normalized to the total area of Earth. The impact area covering all rigidities derived from Rc100k is shown in Movie S1i, and the impact area at 2.5 GV, 5 GV, and 7.5 GV over the last 100 ka is shown in Figure 6b.

When geomagnetic field excursions occurred, the impact area increased significantly. The impact area at 2.5 GV rises to ~90% and ~50% at the midpoints of the Laschamps and NGS excursions, respectively, and shows a slight rise to ~30% during the ML excursion. Moreover, the impact area at 7.5 GV rises to 100% for all these three excursions, which indicates that the geomagnetic field during excursions cannot prevent energetic particles with rigidity higher than 7.5 GV from directly hitting the atmosphere.

The impact area also has hemispheric asymmetry during the NGS and ML excursions. Figure 6c illustrates the impact area at 2.5 GV and 7.5 GV in the northern hemisphere (solid lines) and in the southern hemisphere (dashed lines), respectively. For the first dipole field drop of the NGS excursion (65–63 ka), the impact area is higher in the northern hemisphere. On the contrary, for the second dipole drop (60–58 ka), the impact area is higher in the southern hemisphere. Additionally, the impact area in the southern hemisphere shows a steep increase during the second dipole field drop of the ML excursion (30–28 ka).

4.3 Cosmic ray flux

The cutoff rigidity controls the flux of cosmic ray particles which can access the Earth’s atmosphere. Cosmic ray particles with rigidity lower than the geomagnetic cutoff rigidity cannot reach the Earth’s atmosphere. Using the Rc100k, the variation of global cosmic ray flux over the last 100 ka is reconstructed. The global cosmic ray particle flux is the total number of cosmic ray particles which access the Earth’s atmosphere through the geomagnetic field and is calculated by integration over the cosmic ray spectrum (Shea & Smart, 2004). The galactic cosmic ray particles are mainly protons (~90%) (Gaisser et al., 2016), and we only calculated the cosmic ray proton flux for simplicity. The cosmic ray flux for alpha particles (~9%) and heavier nuclei (~1%) can be calculated in a similar way. The cosmic ray spectrum is mainly modulated by the solar activity, and can be determined by the modulation parameter Φ (Steinhilber et al., 2008). In general, the cosmic ray intensity is high when solar activity is low. We utilized the modulated cosmic ray spectrum at modulation parameter Φ = 400, 500, and 600 MeV, corresponding to the cosmic ray spectrum at solar minimum (see Fig. S7) (Vos & Potgieter, 2015). The solar activity before the Holocene is generally not known, and the solar minimum is chosen for simplicity (Solanki et al., 2004; Usoskin, 2017). The global cosmic ray flux is calculated by integrating the cosmic ray spectrum (Fig. S7) from the local cutoff rigidity to infinity over geographic latitude and longitude:

$$F_{\mathrm{global}}={\int }_{R_{\mathrm{vc}}}^{\mathrm{\infty }}\mathrm{}\pi f\left(R,\mathrm{\Phi }\right)\mathrm{ }\mathrm{d}\mathrm{\Omega }\mathrm{d}R={\int }_{\theta =0}^{\pi }\mathrm{}{\int }_{\phi =0}^{2\pi }\mathrm{}{\int }_{R_{\mathrm{vc}}}^{\mathrm{\infty }}\mathrm{}\pi r^{2}f\left(R,\mathrm{\Phi }\right)\sin \theta \mathrm{d}\theta \mathrm{d}\phi \mathrm{d}R.$$(2)

Here, f(R, Φ) is the cosmic ray intensity at rigidity R and modulation parameter Φ. Radius r = 6391.2 km is 20 km above the Earth’s mean radius. The factor π appears because it is a conversion between downward cosmic ray flux on the top of the atmosphere and the cosmic ray flux in the interplanetary space (see Poluianov et al., 2016, Eq. (3)). Note that the cosmic ray spectrum provided in Figure S7 is in units of MeV and the rigidity-energy conversion equation can be found in Vogt et al. (2007), equation (12) and Gao et al. (2022), equation (15).

Figure 6d illustrates the global cosmic ray intensity over the last 100 ka. The global cosmic ray intensity was significantly increased during the excursions. At the midpoint of the Laschamps excursion (41 ka) and the solar modulation parameter Φ = 400 MeV, the number of cosmic ray particles accessing the Earth’s atmosphere were approximately two times higher than today’s value. During the NGS and the ML excursions, the maximum global flux was ~1.8 and ~1.3 times of today’s value, respectively. Figure 6e shows the global cosmic ray intensity at different rigidities at Φ = 400 MeV over the last 100 ka. The cosmic ray flux access to the atmosphere is higher at low rigidity (~1 GV) because the cosmic ray spectrum shows a maximum value at about 1 GV (~432 MeV). During the Laschamps excursion, a large number of cosmic ray particles accessed the Earth’s atmosphere, with most of the particles in the rigidity range of 0–4 GV, and led to high cosmic ray exposure. At the NGS and ML excursions, the cosmic ray intensity also increased, but most of the particles which accessed the atmosphere had rigidity lower than 2 GV.

5 Discussion and conclusion

Using a combination of four global paleomagnetic field models, i.e., GGF100k, GGFSS70, LSMOD.2, and CALS10k.2, we reconstructed the geomagnetic cutoff rigidity over the last 100 ka in a collection of global grids named Rc100k (the effective vertical cutoff rigidity covering the last 100 ka), which is a quantitative estimation of the geomagnetic field shielding effect. The grids are calculated using a trajectory tracing program, which considers both the dipole and non-dipole components of the geomagnetic field. The most suitable model is selected to reconstruct the cutoff rigidity variation during three excursions, i.e., the Norwegian–Greenland Sea excursion (~65 ka), Laschamps excursion (~41 ka), and Mono Lake/Auckland excursion (~34 ka). After comparing the results for overlapping periods, we found that the model selection is crucial to constrain the cutoff rigidity variation. The comparison of the Rc100k and cutoff rigidity results using models with limited maximum SH degrees confirms earlier findings that the large-scale non-dipole contributions of the geomagnetic field are not negligible for estimating the long-term cutoff rigidity. The Rc100k is used to predict how strongly the impact area and global cosmic ray flux increase during the excursions. At the midpoint of the Laschamps excursion, the global cosmic ray flux reaching the Earth’s atmosphere was up to three times as much as today’s value, while it was about twice as high during the earlier NGS excursion.

The paleomagnetic field models, constructed by compilation of paleomagnetic data, provide a global view of the geomagnetic field evolution, however, with individual spatial and temporal resolution. Moreover, they contain uncertainties that are not well quantified and thus cannot be directly considered in estimating cutoff rigidities. The GGF100k, GGFSS70, LSMOD.2, and CALS10k.2 models are built by a similar modeling strategy (see e.g. Korte & Constable, 2003). Since the paleomagnetic field data always contain uncertainties, regularization in space and time are applied to determine the simplest model that explains the data to the desired level. For the four models utilized in this study, regularization limits the effective spatial resolution to about SH degree 5. When calculating the Rc100k, full resolution models are utilized to reveal the most detailed structure of the cutoff rigidity map, but details should not be over-interpreted, while the large-scale structure of the cutoff rigidity maps is considered to be reliable. The CALS10k.2 model spanning only the recent 10,000 years has the highest temporal resolution, while the GGF100k model has the lowest temporal resolution, clearly lower than 200 years. The Rc100k is a series of cutoff rigidity snapshots over the last 100 ka. Cutoff rigidity values at any given time within the last 100 ka can be derived by using linear interpolation from Rc100k; however, the temporal variation of the Rc100k should be considered no higher than 200 years for the period 100–10 ka and no higher than 50 years for 10–0 ka.

As can be inferred from Figure 1a, the agreement of the models is not always ideal at the ends of the time intervals, and we chose the epochs where different models are joined in some cases to minimize any jumps in particular in the dipole moment strength. At the transition times, i.e., at 67.85 ka, 49.95 ka, 32.65 ka, 17.55 ka, and 9.95 ka, the dipole moment and the dipole axis direction of the adjoining models are very similar (see Table 1 and Fig. 5). However, the combined model unavoidable contains slight jumps at the transition times.

The Rc100k provides the possibility to estimate the cosmogenic isotope (e.g., ¹⁰Be, ¹⁴C, ³⁶Cl) production rate, cosmic ray induced ionization rate, and cosmic ray radiation intensity covering the last 100 ka. Most of the previous studies estimated the cosmogenic isotope production rate using VADM reconstructions, which treated the geomagnetic field as an ideal dipole field (e.g. Gosse & Phillips, 2001), although other studies indicated that the non-dipole field is also important for cosmogenic isotope production rate estimation (Nishiizumi et al., 1989; Dunai, 2000, 2001). A growing number of studies argued that the high cosmic ray flux during excursions enhanced the cosmic radiation dose and the nucleation of clouds and may cause climatic and environmental changes (e.g., Shea & Smart, 2000; Carslaw et al., 2002; Courtillot et al., 2007; Cooper et al., 2021). High cosmic ray ionization may result in a high nucleation rate of sulfuric acid aerosols in the atmosphere, which could, in turn, act as cloud condensation nuclei (see review by Reid, 2000). The cosmic ray flux is modulated by both the solar activity and the geomagnetic field. Using the Rc100k, assuming a constant solar activity, we reconstructed the cosmic ray flux variation over the last 100 ka, providing an opportunity to advance the understanding of the cosmic radiation dose rate and clouds nucleation rate in this time interval. Furthermore, combining the cosmogenic isotope production rate and the cosmic ray flux, the Rc100k would, in turn, be useful for long-term solar activity reconstruction.

The cutoff rigidity calculation contains additional uncertainties in addition to those resulting from uncertainties in the paleomagnetic field models. The effective vertical cutoff rigidity, R_vc, is an approximation of the effective cutoff rigidity, which only considers incident particles in the vertical direction (Rao et al., 1963; Shea et al., 1965). For today’s geomagnetic field, the difference between the R_vc and the R_ac (Apparent cutoff rigidity: an average over cutoff rigidity in nine different directions (Clem et al., 1997)) can be up to 1 GV in the equator region (Clem et al., 1997; Gao et al., 2022). Moreover, the real-time cutoff rigidity contains daily, seasonal, and solar activity variations and could also be influenced by the external field created by current systems outside the Earth, e.g., magnetospheric magnetic fields (Smart et al., 2000; Smart & Shea, 2001, 2005). Those short-term variations (up to decades timescales) of the cutoff rigidity are averaged out on multi-millennial timescales. The most important limitation of cutoff rigidity estimation is the accuracy of the geomagnetic field model (Shea & Smart, 1990). One of the greatest challenges is the characteristics of the regional geomagnetic field excursions, e.g. Post Blake (~98 ka) excursion and Hilina Pali (~19 ka) excursion, which may not be well resolved by current global paleomagnetic field models. Since the amount of paleomagnetic field data increases continuously, future paleomagnetic field models will reveal more details of the geomagnetic field evolution and thus can further improve our understanding of the geomagnetic field shielding effect.

Supplementary materials

Supplementary Tables

Table S1: The effective vertical cutoff rigidity (in unit of GV) calculated using trajectory tracing from a combination of four global paleomagnetic field models, i.e., GGF100k, GGFSS70, LSMOD.2, and CALS10k.2 over the last 100 ka. Timing resolution is 200 year before 10 ka, and 50 year after 10 ka.

Table S2: The effective vertical cutoff rigidity (in unit of GV) calculated using the GGF100k model. Timing resolution is 200 year.

Table S3: The effective vertical cutoff rigidity (in unit of GV) calculated using the GGFSS70 model. Timing resolution is 200 year.

Table S4: The effective vertical cutoff rigidity (in unit of GV) calculated using the LSMOD.2 model. Timing resolution is 200 year.

Table S5: The effective vertical cutoff rigidity (in unit of GV) calculated using the CALS10k.2 model. Timing resolution is 50 year.

Movie S1. The effective vertical cutoff rigidity R_vc contour maps derived from (a) GGF100k, (b) GGFSS70, (c) LSMOD.2, (d) IMOLEe, (e) CALS10k.2, and (f) Rc100k. (g) The latitudinal structure of Rc100k. (h) The longitudinal structure of Rc100k. (i) the impact area derived from Rc100k. Access here

Movie S2. The geomagnetic field radial component Br from (a) GGF100k, (b) GGFSS70, (c) LSMOD.2, (d) IMOLEe, (e) CALS10k.2, and (f) the combined model used to calculated the Rc100k over the last 100 ka. Access here

Movie S3. The geomagnetic field intensity B_totai from (a) GGF100k, (b) GGFSS70, (c) LSMOD.2, (d) IMOLEe, (e) CALS10k.2, and (f) the combined model used to calculated the Rc100k over the last 100 ka. Access here

Movie S4. The cutoff rigidity difference between (a) GGF100k and GGFSS70, (b) GGFSS70 and LSMOD.2, (c) LSMOD.2 and IMOLEe, and (d) CALS10k.2 and GGF100k. Access here

Movie S5. The R_vc contour maps using models with maximum spherical harmonic (SH) degree (a) geocentric axial dipole (GAD), (b) 1, (c) 2, and (d) 3. The difference between Rc100k and cutoff rigidity results using model with maximum SH degrees (e) geocentric axial dipole (GAD), (f) 1, (g) 2, and (h) 3. Access here

Figure S1. Geomagnetic field radial component Br maps from different paleomagnetic field models at (a) 67.85 ka, (b) 49.95 ka, (c) 32.65 ka, (d) 17.55 ka, (e) 9.95 ka. The left and middle panels show the Br maps of the adjoining models. The differences between two models are shown in the right panels.

Figure S2. Geomagnetic field intensity Btotai maps from different paleomagnetic field models at (a) 67.85 ka, (b) 49.95 ka, (c) 32.65 ka, (d) 17.55 ka, (e) 9.95 ka. The left and middle panels show the Btotal maps of the adjoining models. The differences between two models are shown in the right panels.

Figure S3. Effective vertical cutoff rigidity Rvc contour maps derived from different paleomagnetic field models at (a) 67.85 ka, (b) 49.95 ka, (c) 32.65 ka, (d) 17.55 ka, (e) 9.95 ka. The left and middle panels show the Rvc maps derived from the adjoining models. The differences between two Rvc maps are shown in the right panels.

Figure S4. The global average of the value differences Dvc between different paleomag- netic field models. The difference between Rc100k and cutoff rigidity calculated using only the geocentric axial dipole (GAD) from the combined model is also shown for comparison.

Figure S5. (a) The global average of the value differences $\overline{D_{\mathrm{vc}}}$ and (b) the global average of the percentage differences $\overline{P_{\mathrm{vc}}}$ between Rc100k and cutoff rigidity calculated using the truncated paleomagnetic field model, up to spherical harmonic degree 3, 2, 1, and geocentric axial dipole (GAD).

Figure S6. Mean value of the cutoff rigidity at Earth’s equator versus dipole moment of the geomagnetic field over the last 100 ka.

Figure S7. Cosmic ray spectrum for protons at solar modulation parameter $ = 300, 400, 500 and 600 MeV.

Acknowledgments

This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB 41000000), the National Natural Science Foundation of China (Grant Nos. 41922031, 41774188), the Key Research Program of the Institute of Geology and Geophysics, CAS (Grant Nos. IGGCAS-201904, IGGCAS-202102), and the Key Research Program of Chinese Academy of Sciences (Grant No. ZDBS-SSW-TLC00103). J. Gao acknowledges financial support from China Scholarship Council (202004910583). S. Panovska acknowledges the Discovery Fellowship at the GFZ Potsdam, Germany. The authors gratefully acknowledge D. Smart and M. Shea for providing a modified version of the trajectory tracing program. The geomagnetic core field models are available at https://doi.org/10.5880/GFZ.2.3.2020.005. The cutoff rigidity results, Rc100k, covering the last 100 ka are available from the EarthRef Digital Archive at https://earthref.org/ERDA/2521/ or from the supporting information of this paper. The editor thanks Alexandre Fournier and an anonymous reviewer for their assistance in evaluating this paper.

References

All Tables

All Figures

Fig. 1 Dipole power (a) and non-dipole power (b) evaluated at Earth’s surface over the last 100 ka, as predicted by five continuous spherical harmonic geomagnetic field models: IMOLEe (Leonhardt et al., 2009), CALS10k.2 (Constable et al., 2016), GGF100k (Panovska et al., 2018b), LSMOD.2 (Korte et al., 2019b), and GGFSS70 (Panovska et al., 2021). Five global and regional excursions are labeled on the top panel a: PB: Post Blake, NGS: Norwegian–Greenland Sea, La: Laschamps, ML: Mono Lake/Auckland, and HP: Hilina Pali. Gray areas indicate the three excursions studied in more detail here: NGS, La, and ML. The time intervals of the model, and models used to calculate Rc100k are given by the color bars at the bottom.

In the text

	Fig. 2 Effective vertical cutoff rigidity contour maps covering the period of Laschamps excursion (43.95–39.45 ka) calculated from the (a) GGF100k, (b) GGFSS70, (c) LSMOD.2, and (d) IMOLEe models.
In the text

	Fig. 3 Effective vertical cutoff rigidity contour maps covering the period of the Norwegian–Greenland Sea excursion (66–63 ka) calculated from the (a) GGF100k and (b) GGFSS70 models.
In the text

	Fig. 4 Effective vertical cutoff rigidity contour maps covering the period of the Mono Lake/Auckland excursion (36.45–33.45 ka) calculated from the (a) GGF100k, (b) GGFSS70, and (c) LSMOD.2 models.
In the text

Fig. 5 Axial dipole moment (ADM) evolution with time as combination from four paleomagnetic field models: GGF100k (100–67.85 ka and 17.55–9.95 ka), GGFSS70 (67.85–49.95 ka and 32.66–17.55 ka), LSMOD.2 (49.95–32.65 ka), and CALS10k.2 (9.95–0 ka). Virtual axial dipole moment (VADM) stacks obtained from globally distributed relative paleointensity (RPI) records are also shown for comparison: GLOPIS-75 (Laj et al., 2004), C2018-Overall stack (Channell et al., 2018), sint-800 (Guyodo & Valet, 1999), and sint-2000 (Channell et al., 2018).

In the text

Fig. 6 (a) Mean value of Rvc over longitude at latitudes 45° N, 0° N, and 45° S; (b) impact area of energetic particles at Rvc equal to 2.5 GV, 5 GV, and 7.5 GV; (c) impact area of energetic particles at Rvc equal to 2.5 GV and 5 GV, with solid lines (dashed lines) indicating the impact area in the northern hemisphere (southern hemisphere), respectively; (d) global cosmic ray proton flux accessing the Earth’s atmosphere; (e) cosmic ray proton flux reaching the Earth’s atmosphere due to its rigidities during the last 100 ka. For comparison, today’s values calculated using the International Geomagnetic Reference Field (IGRF) at epoch 2020 (Alken et al., 2021) are included in panels a and b as dashed lines.

In the text

	Figure S1. Geomagnetic field radial component Br maps from different paleomagnetic field models at (a) 67.85 ka, (b) 49.95 ka, (c) 32.65 ka, (d) 17.55 ka, (e) 9.95 ka. The left and middle panels show the Br maps of the adjoining models. The differences between two models are shown in the right panels.
In the text

	Figure S2. Geomagnetic field intensity Btotai maps from different paleomagnetic field models at (a) 67.85 ka, (b) 49.95 ka, (c) 32.65 ka, (d) 17.55 ka, (e) 9.95 ka. The left and middle panels show the Btotal maps of the adjoining models. The differences between two models are shown in the right panels.
In the text

	Figure S3. Effective vertical cutoff rigidity Rvc contour maps derived from different paleomagnetic field models at (a) 67.85 ka, (b) 49.95 ka, (c) 32.65 ka, (d) 17.55 ka, (e) 9.95 ka. The left and middle panels show the Rvc maps derived from the adjoining models. The differences between two Rvc maps are shown in the right panels.
In the text

	Figure S4. The global average of the value differences Dvc between different paleomag- netic field models. The difference between Rc100k and cutoff rigidity calculated using only the geocentric axial dipole (GAD) from the combined model is also shown for comparison.
In the text

	Figure S5. (a) The global average of the value differences $\overline{D_{\mathrm{vc}}}$ and (b) the global average of the percentage differences $\overline{P_{\mathrm{vc}}}$ between Rc100k and cutoff rigidity calculated using the truncated paleomagnetic field model, up to spherical harmonic degree 3, 2, 1, and geocentric axial dipole (GAD).
In the text

	Figure S6. Mean value of the cutoff rigidity at Earth’s equator versus dipole moment of the geomagnetic field over the last 100 ka.
In the text

	Figure S7. Cosmic ray spectrum for protons at solar modulation parameter $ = 300, 400, 500 and 600 MeV.
In the text