Effects of the Laschamps Excursion on Geomagnetic Cutoff Rigidities

Today's geomagnetic field can prevent energetic particles, including solar energetic particles and galactic cosmic rays, from directly hitting the Earth's atmosphere. However, when the geomagnetic field strength is significantly decreased during geomagnetic field excursions or reversals, the geomagnetic field shielding effect becomes less prominent. Geomagnetic cutoff rigidity, as a quantitative estimation of the shielding effect, can be calculated using trajectory tracing or theoretical equations. We use a recent high‐resolution continuous geomagnetic field model (LSMOD.2) to study the geomagnetic cutoff rigidity during the Laschamps excursion. Global grids of the geomagnetic cutoff rigidities are presented, in particular for the excursion midpoint when the geomagnetic field is weak and not dipole‐dominated anymore at Earth's surface. We compare the cutoff rigidity calculation results between a trajectory tracing program and theoretical equations and we find that the influence of the non‐dipole component of the geomagnetic field cannot be ignored during the excursion. Our results indicate that the exposure of Earth's atmosphere to energetic particles of cosmic and solar origin is high and nearly independent of latitude in the middle of the Laschamps excursion. Our results will be useful for future studies associated with cosmic radiation dose rate and cosmogenic isotope production rate during the Laschamps excursion.

defined as follows. Only particles with rigidity higher than a certain threshold value can reach the Earth's atmosphere, and that threshold value is called geomagnetic cutoff rigidity. The geomagnetic cutoff rigidity has been first studied by Størmer (1955) and Elsasser et al. (1956) through theoretical analysis and then quantitatively determined by trajectory tracing (Shea et al., 1965). The most important limitation of cutoff rigidity calculation is the accuracy of the geomagnetic field model (Shea & Smart, 1990).
Previous studies about energetic particles during geomagnetic polarity transitions utilized simplified magnetic field models (Glassmeier et al., 2009;Stadelmann et al., 2010;Vogt & Glassmeier, 2000). Considering the configuration of the paleomagnetic field as a reduced dipole or quadrupole, the cutoff rigidity during a simulated polarity reversal process was studied by Vogt et al. (2007). Stadelmann et al. (2010) treated a mixed dipole and quadrupole as a possible configuration during geomagnetic polarity transitions and found that the region where energetic particles penetrate the geomagnetic field depends on the geomagnetic field configuration. During reversals and excursions, the configuration of the geomagnetic field can significantly deviate from dipolar or quadrupolar, the morphology of the geomagnetic field may be highly complex and higher-degree multipole components of the paleomagnetic field become important .
Recently, a spherical harmonic (SH) paleomagnetic field model based on sediment and volcanic data and covering the Laschamps and Mono Lake excursions (in the past 50-30 thousand years) was presented by Brown et al. (2018), with an updated version by . This model, called LSMOD.2, provides the possibility to study the geomagnetic field shielding against energetic particles during the Laschamps excursion based on a realistic global field structure. The model suggests that the axial dipole strength is close to zero in the middle of the Laschamps excursion.
In this study, we calculated global grids of geomagnetic cutoff rigidities using the LSMOD.2 model by both trajectory tracing and two simplified theoretical equations. Furthermore, latitudinal and longitudinal profiles of the cutoff rigidity and the impact area of the energetic particles during the Laschamps excursion are presented. The geomagnetic cutoff rigidity strongly influences the cosmic radiation dose rate (e.g.,  and cosmogenic isotope (e.g., 10 Be, 14 C, 36 Cl) production rate (Poluianov et al., 2016), which is the basis of the cosmogenic nuclides dating method (Gosse & Phillips, 2001). Moreover, it plays a role for cloud condensation nuclei (e.g., Kirkby et al., 2011;Pierce & Adams, 2009) and possible biospheric effects (e.g., Wei et al., 2014) during the excursion. Our cutoff rigidity results provide the possibility to accurately estimate the cosmogenic isotope production rate during the excursion. This paper is organized as follows. The magnetic core field model (LSMOD.2) is introduced in Section 2. The cutoff rigidity calculation using two simplified equations is presented in Section 3, and the calculation using a trajectory tracing program is described in Section 4. The cutoff rigidity results are compared in Section 5. Discussion and conclusion are given in Section 6.

Magnetic Core Field Model
Spherical harmonic functions are widely used to represent the global distribution and evolution of the geomagnetic core field. The basic methodology of constructing the LSMOD.2 model has been used for a large range of global continuous magnetic field models (see reviews by Constable and Korte [2015] and Panovska et al. [2019]). Here we briefly outline the modeling method, more details can be found in Bloxham and Jackson (1992) and Constable (2003, 2008).
In a non-magnetic source region, the time-depending geomagnetic core field vector B(t), neglecting the crustal and external fields, can be expressed as the gradient of a scalar potential V(t). In spherical polar coordinates (r, θ, φ), V(t) can be expanded as ( ) +1 ( cos ( ) + ℎ sin ( )) × (cos ) ( ) where a = 6371.2 km is the mean radius of Earth's surface. (cos ) are the Schmidt quasi-normalized associated Legendre functions of degree l and order m. The maximum SH degree L max determines the maximum spatial resolution of the model. For a given epoch t, the coefficients and are related to the SH coefficients and (also called Gauss coefficients) as M k (t) are cubic B-splines temporal basis functions, and the ith component M i is nonzero only if t lies in the interval (t i , t i + 4 ) of the knot points t k , k = 1 to K. The LSMOD.2 model utilizes equally spaced knot points every 50 years, which determines the maximum temporal resolution of the model. In practice, however, regularizations in both space and time are applied, that trade off the fit to the data against smoothness of the model, because the data contain uncertainties, and fitting them too closely would lead to artificial structure in the model. Therefore, the regularization is used to determine the simplest model that explains the data to the desired level, and it determines the effective spatial and temporal resolution. In the case of LSMOD.2, the regularization limits the spatial resolution to about SH degree 5. The model spans the interval 50 to 30 ka Before Present (BP). In the paleomagnetic community, Before Present time t BP is widely used. The conversion from Anno Domini/Before Christ (AD/ BC) time to Before Present (BP) time is t BP = 1950 − t AD/BC , that is, at 1950 AD, the BP time is 0. In this paper, unless otherwise stated, BP time in unit of thousand years (ka) is used.
The SH coefficients of degrees l = 1, 2, 3 describe the dipole, quadrupole, and octupole magnetic field contributions, respectively. Using SH functions, some field characteristics can be easily represented (Lowes, 1966). For example, the dipole power (DP) and non-dipole power (NDP), evaluated at radius r, are given by There is no clear definition for the start and end of geomagnetic transitional events, in particular in the global view. Several parameters have been suggested to indicate global transitional field behavior, such as the globally averaged paleosecular variation (PSV) index (Panovska & Constable, 2017), the time at which the power spectrum becomes flat (Brown et al., 2018), and the ratio between DP and large-scale NDP at the core-mantle boundary (CMB; . Using the LSMOD.2 model and the PSV index,  suggested that at Earth's surface the Laschamps excursion starts at 41.9 ka, reaches its midpoint at 40.6 ka, and ends at 40.1 ka, with a duration of 1.8 kyr (thousand years). Considering DP/NDP at the CMB, however, they found a duration of 4.8 kyr, and the minimum dipole moment is seen at 40.95 ka in LSMOD.2. Figure 1a shows the DP and NDP variation during the Laschamps excursion from LSMOD.2 evaluated at Earth's surface. It is clear that before and after the Laschamps excursion, the DP is significantly larger than NDP at Earth's surface, which indicates that the geomagnetic field is dipole-dominated. At the excursion midpoint, DP is smaller than NDP and reaches its minimum, so that the non-dipole field contributions dominate for a while even at Earth's surface during the excursional process.
The radial component B r and magnetic field intensity B total  of the LSMOD.2 model over the entire interval 50-30 ka are provided in Movies S1a and S1b, respectively. Panels b-d of Figure 1 show three snapshots of B r at Earth's surface at 48.35, 40.95, and 35.55 ka, respectively. These three epochs illustrate the dipole-dominance before and after the Laschamps excursion and the complex field structure in the middle of the event. LSMOD.2 indicates that the excursion process is governed by axial dipole decay and then recovery. The field intensity is very weak in the middle of the excursion. According to the LSMOD.2 model, the dipole component is almost zero at the excursion midpoint but does not reverse to opposite polarity.

Cutoff Rigidity Calculation Using Theoretical Equations
There are different types of energetic particles in the heliosphere, such as galactic cosmic rays (GCR), typically in the GeV range, solar energetic particles (SEP), typically in the MeV range, and anomalous cosmic rays (ACR), typically in the MeV range (Mironova et al., 2015;Simnett, 2017). Those energetic particles continuously bombard the Earth, and the geomagnetic field is a natural shield against them. The primary energetic particles are mainly protons and therefore get deflected by the geomagnetic field (Gaisser et al., 2016;Lifton et al., 2008). Due to the high electric conductivity in the magnetosphere (Bütikofer, 2018), the electric force can be ignored, and the electromagnetic force on a charged particle is described by the Lorentz force where m 0 is the rest mass of the particle, Z is the charge of the particle, e is the elementary charge, is the particle speed and ⃗ is the magnetic field vector. In the special relativity frame (Lambourne, 2010), γ is the Lorentz factor and = 1 . c is the speed of light.
The rigidity is defined as R = pc/Ze, where p = γm 0 v is the momentum of the particle. In a constant magnetic field, the particles with equal rigidity, equal charge sign, and equal initial conditions have the same moving trajectory. Particles that have rigidity higher than the geomagnetic cutoff rigidity can penetrate the geomagnetic field and hit the Earth's atmosphere, and particles with rigidity lower than the geomagnetic cutoff rigidity cannot reach the Earth's atmosphere.
Theoretical studies about cutoff rigidity calculation, starting from the pioneering work by Størmer (1955), found that the geomagnetic cutoff rigidity is inversely correlated with geomagnetic latitude λ. The present geomagnetic field can be treated as a dipole field as a first-order approximation. Elsasser et al. (1956) calculated the cutoff rigidity in a dipole field as where R cd (in unit of GV, 10 9 V) is the cutoff rigidity calculated in dipole configuration. μ 0 = 4π × 10 −7 N/A 2 is the permeability of free space. r (in unit of km) is the radius distance from the Earth center. M (in unit of Am 2 ) is the dipole moment of the geomagnetic field. The geomagnetic field dipole moment is not constant and has decreased ∼10% over the past 200 years (Brown et al., 2018). Using the changing dipole moment, the time variation R cd can be quantitatively calculated using Equation 6. The cutoff rigidity calculated through dipole approximation is widely used in cosmic ray and cosmogenic isotope production rate research (e.g., Gosse & Phillips, 2001;Lal, 1991).
Although the present geomagnetic field can be treated as a dipole field as a first-order approximation, the non-dipole contributions of the geomagnetic field cannot be ignored. A clear example of the regionally different non-dipole field influence is the intensity low in South Atlantic Anomaly (SAA) region. The non-dipole field contributions can reach about 20% during the last centuries and millennia  and could have a strong influence on the geomagnetic cutoff rigidity calculation (Shea et al., 1968). Following Rothwell (1958) and Dunai (2000Dunai ( , 2001 considered both the dipole and non-dipole field and used the local geomagnetic field horizontal intensity and inclination to calculate the cutoff rigidity R cnd (in units of GV) as = 0 16 × 10 11 (1 + 0.25tan 2 ) 1.5 , where H is the geomagnetic field horizontal intensity (in units of nT) and I is the geomagnetic field inclination (in units of radians). At a given location, if the local geomagnetic field vector is known, R cd and R cnd can be calculated using Equations 6 and 7. Note that the actual cutoff rigidity additionally contains daily, seasonal, and solar activity variations . The cutoff rigidity calculations from core field models thus only give an average condition about the real-time cutoff rigidity, both in the case of using these equations or trajectory tracing.

Cutoff Rigidity Calculation Using Trajectory Tracing
The more rigorous method to determine the cutoff rigidity is by tracing the trajectory of the energetic particles, and trajectory tracing has been widely used for cosmic ray research , 2004Shea et al., 1965;Smart & Shea, 1994. In this section, we briefly introduce the basic methodology of the trajectory tracing method, which is available at https://ccmc.gsfc.nasa.gov/modelweb/sun/cutoff.html (see e.g., Smart et al. [2000] for more details).
The trajectory of a charged particle moving in a magnetic field (Equation 5), is a problem that has no mathematical closed-form solution . Since it is difficult to calculate the trajectory of an inward-directed charged particle through the geomagnetic field, the common method to determine the trajectory is by tracing a charged particle with a negative charge, as the trajectory of an outward negatively charged particle is identical to the trajectory of an inward positively charged particle. The starting point of the particle tracing is from 20 km above the Earth's surface because the Earth's atmosphere is considered to extend 20 km above the surface (geodetic coordinates are used to calculate the Earth's surface elevation [Moritz, 1980]). The particle trajectory is traced until the particle reaches the outer boundary (defined as 25 earth radii) or returns to the Earth's atmosphere (20 km above the surface). The Bulirsch-Stoer numerical integration technique for the fourth-order Runge-Kutta iteration is used for tracing the particle trajectory (Stoer & Bulirsch, 2013), and the length of each trajectory step is approximately 1% of the gyro-radius of the charged particle.
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. On the contrary, if the trajectory returns to the atmosphere, the particle is called a "forbidden particle" and a particle at that rigidity cannot penetrate the geomagnetic field. After tracing the particle trajectories at discrete steps from high rigidity to low rigidity through the rigidity spectrum (the energy spectrum of the energetic particles), the geomagnetic cutoff rigidity is described as R u , R l , and R c . The upper rigidity R u is the highest forbidden rigidity value. The lower rigidity R l is the lowest allowed rigidity value. R c is the effective cutoff rigidity, which is an average over the rigidity penumbra and describes the transparency of the penumbra (Cooke et al., 1991). The rigidity penumbra is a series of allowed and forbidden bands between R u and R l . The detailed calculation method of R c is shown in Section 4.1. Since cutoff rigidity calculations are in general based on the present geomagnetic field, we performed several tests to ensure the robustness and reliability of calculating R c using the trajectory tracing program for the Laschamps excursion, including the effect of rigidity scan step size and cutoff sky map that indicates cutoff rigidity for different particle incoming angles. Moreover, we test whether the Earth's crustal field would have to be considered when the core field is weak during the excursion.
Apart from tracing an outward-directed particle starting near the Earth's surface, Stadelmann et al. (2010) found a method to trace inward-directed particles starting outside the magnetosphere. They traced the trajectories of several thousand energetic particles which were homogeneously distributed on the magnetospheric boundary and then calculated the cutoff rigidity map. They found that the cutoff latitudes resulting from their method and the results from the trajectory tracing program (Smart & Shea, 2001) are similar. Figure 2 shows the penumbra of vertical cutoff rigidity at five different latitudes at three epochs, that is, before, in the middle, and after the Laschamps excursion. Vertical cutoff rigidity means that the trajectory of the outward-directed particle was initiated in the vertical direction, that is, radial from the Earth's surface. The cutoff rigidity penumbra is not sharp. There are chaotic structures of the penumbra bands between R u and R l , which indicates that both allowed and forbidden rigidities exist in the penumbra region.

Penumbra Structure of Cutoff Rigidity
The effective vertical cutoff rigidity, R vc , is calculated as (Bütikofer, 2018) for forbidden rigidity 1, for allowed rigidity.
The LSMOD.2 model indicates that the geomagnetic field is dipole-dominated at 48.35 and 35.33 ka and is relatively weak and non-uniformly distributed at 40.95 ka. This is reflected by R vc , which is about 13 GV in the equatorial region and nearly zero in both polar regions before and after the Laschamps excursion, but close to zero in the equatorial region and no longer zero in the polar regions in the middle of the excursion. The penumbra Rigidity (GV) region is also wider at the middle latitude (latitude ±40°) when the geomagnetic field is dipole-dominated and is relatively narrow in all latitudes when the geomagnetic field is weak.

First and Second Rigidity Scan
The resolution of the rigidity scan has long depended on available computing power  and high-resolution trajectory tracing is still time-consuming even on a modern computer. To decrease the computation time, a hybrid rigidity step is used in this study.
The penumbra bands are chaotic between R u and R l and a high-resolution scan of the rigidity should be applied in this energy range. However, for rigidity higher than R u or lower than R l , the particle is fully allowed or fully forbidden and the resolution of the rigidity scan can be lower. Thus, a hybrid method using first and second rigidity scans is applied. For the first rigidity scan, a 0.5 GV rigidity interval was selected and the rigidity was scanned from 30 to 0 GV. The scanning range covers the spectrum of the energetic particles (Mironova et al., 2015). After the first rigidity scan, the first R u and R l are obtained. Then, the second rigidity scan with 0.01 GV rigidity interval from first R u plus 0.5 GV to first R l minus 0.5 GV is applied. The penumbra structure can be well determined at 0.01 GV resolution (Shea et al., 1968). The final R u and R l are determined from the second scan and R vc is calculated by Equation 8. Figure S1 in Supporting Information S1 shows the penumbra after the first rigidity scan and Figure 2 shows the penumbra after the second rigidity scan. The first rigidity scan gives the cutoff rigidity with low resolution, while the penumbra structure between R u and R l is well determined after the second rigidity scan. The cutoff rigidity R vc resulting from the first and second scans during the Laschamps excursion are shown in Movies S1c and S1d. The difference of R vc between the first and second scans is shown in Movie S1e. For most areas of the Earth's surface, the difference of R vc is lower than 0.4 GV during the Laschamps excursion. The difference is particularly small in the middle of the Laschamps excursion.

Cutoff Sky Map
The cutoff rigidity at any given location is a function of the zenith and azimuth angles of the energetic particles' incidence. A complete calculation of the effective cutoff rigidity needs to consider multiple directions and the differential rigidity spectrum of the energetic particles radiation at different directions (Rao et al., 1963;Shea et al., 1965). However, tracing the particles for different directions is demanding and the vertical cutoff rigidity is the usual approximation for the effective cutoff rigidity (Smart & Shea, 2009). To test whether this is a valid approximation in times of weak transitional geomagnetic field, we give several examples of the cutoff sky map and compare the difference between the apparent cutoff rigidity and the vertical cutoff rigidity.
At a given location, the cutoff sky map is a systematic calculation of the effective cutoff rigidity for different directions of energetic particles arrival . Figure 3 shows cutoff sky maps at longitude 0°E, latitude 40°N, 0°N, and 40°S at three epochs, that is, before, at the middle, and after the Laschamps excursion. In each cutoff sky map, the radial coordinate indicates the zenith angle and the angular coordinate indicates the azimuth angle. The cutoff sky map with arrival zenith angle larger than 75° is excluded because the energetic particle flux near the sea level is negligible Smart & Shea, 2009). At each direction, the effective cutoff rigidity is calculated by trajectory tracing.
As shown in Figure 3, at the midpoint of the Laschamps excursion, the cutoff rigidity is low at all latitudes and some irregular structure is found in the cutoff sky maps. This may result from higher degree multipole components of the geomagnetic field, which indicates that non-dipole field contributions become important for the shielding effect during the Laschamps excursion. In times of dipole-dominated field before and after the excursion, the cutoff rigidity is higher in the east direction and lower in the west direction. The higher cutoff rigidity in the east direction results from the Lorentz force (Equation 5), which is radially outward on the incident particles in the east direction.
The apparent cutoff rigidity is an approximation for the effective cutoff rigidity which considers the different incidence directions. Following Bieber et al. (1997), the apparent cutoff rigidity, R ac , can be calculated as an average over the cutoff rigidity in the vertical direction and the cutoff rigidity in eight azimuthal directions at zenith angle 30°. Those directions are marked by red dots on the cutoff sky map. R vc and R ac are both approximations for the effective cutoff rigidity, and R vc and R ac during the Laschamps excursion are shown in Movies S2c and S2d. The difference between R vc and R ac is shown in Movie S2e. Furthermore, the cutoff rigidity for nine directions, which is used to calculate the apparent cutoff rigidity, is shown in Movie S3.
When the geomagnetic field is dipole-dominated, cutoff sky maps show that the gradients of cutoff rigidity at higher values are steeper. This effect tends to make R ac higher than R vc , especially in the equatorial region where the cutoff rigidity is high. The systematic deviation between R ac and R vc is shown in Movie S2e. R ac is approximately 1 GV higher than R vc in the equatorial region before and after the Laschamps excursion. However, a systematic deviation is not obvious at middle and high latitudes. The systematic deviation becomes lower when the geomagnetic field is weaker and is negligible at the excursion midpoint.
For many purposes, R vc is sufficient for studying the shielding effect of the geomagnetic field (Smart & Shea, 2005). Movie S2e demonstrates that the difference between R vc and R ac is about 1 GV near the equator before and after the excursion, while the background value is about 13 GV. During the Laschamps excursion, the difference is particularly small. This indicates that both R vc and R ac are good approximations for the effective cutoff rigidity even considering the different incidence directions, especially during the Laschamps excursion. However, if a high accuracy result of the cutoff rigidity is required, such as accuracy higher than 1 GV in the equatorial region when the geomagnetic field is dipole-dominated, the R ac would be more appropriate.

Effect of the Crustal Field
In addition to the core magnetic field, the Earth also possesses the crustal field which is generated by magnetized rocks in the Earth's lithosphere. Since the crustal field is much smaller than the present core field, the effect of the crustal field is ignored when calculating the cutoff rigidity. However, during a reversal or excursion, the core field is rather weak. In this section, we will examine whether the crustal field can also be ignored during the Laschamps excursion when calculating the cutoff rigidity.
We use the crustal field estimate from the global Kalmag model (Baerenzung et al., 2020), an SH magnetic field model that distinguishes some of the different field sources of the Earth's magnetic field. SH coefficients up to degree 100 are used to calculate the Earth's crustal field. Figure 4a illustrates the global distribution of the crustal field radial component B r evaluated at Earth's surface calculated from this model. Compared to the core field, the crustal field is small-scale and globally inhomogeneously distributed.
The trajectory tracing becomes computationally time-consuming when adding the crustal field model in the calculations. Therefore, we only study the effect of the crustal field in a small region where the crustal field has one of the strongest anomalies on Earth, called Bangui magnetic anomaly (Ouabego et al., 2013) and shown in the black box in Figure 4a. The regional crustal field B r map and regional core field B r map are shown in Figures 4b  and 4c, respectively.
The time of the LSMOD.2 core field model shown in Figure 4c is chosen at 41.05 ka, when the core field intensity in the studied region is lowest. The geomagnetic core field is still mostly 10 times larger than the crustal field. Figures 4d and 4e show the vertical cutoff rigidity R vc maps with and without considering the crustal field, respectively. The difference between them is lower than 0.1 GV. We additionally studied a slightly different time, 40.95 ka, which is the midpoint of the Laschamps excursion (see Figure S2 in Supporting Information S1). The influence of the crustal field on cutoff rigidity is also negligible at this epoch. We conclude that the crustal field is insignificant in calculating the cutoff rigidity even during magnetic field transitions and we ignored the Earth's crustal field when calculating the cutoff rigidity during the Laschamps excursion in the following.

Results
We have calculated global grids (5° in latitude, 5° in longitude) of the cutoff rigidities over the duration of the LSMOD.2 model (50-30 ka), including the Laschamps excursion, in 200 year steps and at an altitude of 20 km above the Earth's surface, using both the simplified equations and trajectory tracing. We will compare the results in this section. Moreover, we investigate results calculated using the LSMOD.2 model with different maximum SH degree. We also present latitudinal and longitudinal structures, impact area, and particle trajectories during the Laschamps excursion. Figure 5 shows the cutoff rigidity contour maps before, during and after the Laschamps excursion. Apparent cutoff rigidity R ac , cutoff rigidity in a dipole field R cd , cutoff rigidity considering the non-dipole field R cnd , the difference of cutoff rigidity value between R ac and R cnd (D ac,cnd = R ac -R cnd ) and percentage differences between R ac and R cnd (P ac,cnd ) at three epochs are shown, respectively. P ac,cnd is given by Lifton et al. (2008)  As shown in Figure 5 and Movie S2, while the maps of R ac and R cnd look similar, R cd clearly gives an extremely simplified estimate. However, the difference D ac,cnd does not seem negligible. While the geomagnetic field is dipole-dominated, the maximum value of D ac,cnd for example, reaches about 4 GV in the equatorial region at 48.35 ka. At middle and low latitudes, D ac,cnd is relatively lower before and after the excursion. D ac,cnd is larger where the intensity of the geomagnetic field shows a larger gradient. For each epoch, apparent cutoff rigidity R ac , cutoff rigidity in a dipole field R cd , cutoff rigidity considering the non-dipole field R cnd , value differences between R ac and R cnd (D ac,cnd ) and percentage differences between R ac and R cnd (P ac,cnd ) are shown. The cutoff rigidity is calculated at the same altitude (20 km above the Earth's surface) for both trajectory tracing and using the simplified equations.

Cutoff Rigidity During Laschamps Excursion
In the middle of the Laschamps excursion (40.95 ka), both R ac and R cnd maps show similar patterns of cutoff rigidity provided mainly by non-dipole geomagnetic field, and R cd is globally close to zero. D ac,cnd becomes lower when the geomagnetic field intensity decreases, but P ac,cnd becomes large, exceeding 100% in several areas. After the Laschamps excursion, D ac,cnd increases and P ac,cnd decreases as the field gets more dipolar again. Since the post-excursion dipole moment is lower than its pre-Laschamps strength, both R ac and D ac,cnd are lower than before the excursion.
The cutoff rigidity changing rate, averaged in 1 kyr windows, is provided in Movies S1f. In the equatorial region, the cutoff rigidity changing rate is about −10 GV/kyr during the dipole field decrease (43-41 ka), and about 5 GV/kyr during the dipole field increase (41-36 ka). Mean value of the cutoff rigidity at Earth's equator versus dipole moment of the geomagnetic field from LSMOD.2 model is provided in Figure S3 of Supporting Information S1. The cutoff rigidity at equator changes nearly linear with the dipole moment of geomagnetic field.
For quantitative assessment the global averages of the differences D ac,cnd and P ac,cnd , and , respectively, are calculated as Similarly, the global averages and are determined (see Figure S4 in Supporting Information S1). Figure 6 shows the variation of , , between 50 and 32 ka. The difference varies with the magnetic (dipole) field strength, while becomes large during the excursion when the field intensity is low. During 42-40 ka, is higher than 100%, indicating that R cnd is no longer a reliable estimation of the cutoff rigidity when the surface geomagnetic field is weak.

Cutoff Rigidity Calculated Using LSMOD.2 Model With Different Maximum SH Degrees
To investigate the effect of the high degree components of the LSMOD.2 model, we calculated the cutoff rigidities using LSMOD.2 truncated at different maximum SH degrees by using the trajectory tracing program. Figure 7 shows maps of R vc using LSMOD.2 up to SH degree 5 (effective model resolution), 3, 1 (tilted dipole), and for only the first coefficient 0 1 (geocentric axial dipole, GAD) for three epochs, respectively. R vc contour is the global average of the value differences between R ac and R cnd .
is the global average of the percentage differences between R ac and R cnd . Obviously, the resolution of the cutoff maps decreases when the LSMOD.2 model with lower L max is used. The R vc maps using LSMOD.2 up to SH degree 10 (Movie S1d) and SH degree 5 are very similar because of the applied regularization, that limits the effective spatial model resolution to about SH degree 5, so that in fact the SH coefficients larger than SH degree 5 of LSMOD.2 are insignificant for the cutoff rigidity calculation. A truncation to only the first SH coefficient (geocentric axial dipole, GAD) gives similar results as the dipole equation. In the middle of the excursion, both the R vc maps calculated using the model up to L max = 1 and GAD are nearly zero. The geomagnetic field configuration during the Laschamps excursion is not dipole-dominated; thus, as Figures 7e and 7f show, non-dipole components have to be considered when studying the cutoff rigidities during the excursion.

Latitudinal and Longitudinal Profiles, and Impact Area
Based on the R vc results during the Laschamps excursion, we can further investigate the latitudinal and longitudinal structure of the cutoff rigidity and the impact area of the energetic particles during the excursion (Smart &  Shea, 2009;Stadelmann et al., 2010). Figure 8 shows the mean values of the R vc , R u , and R l over latitude (a, d, and g) and longitude (b, e, and h) for the three epochs before, during and after the excursion (See Movies S1g and S1h for those values over the whole time interval).
Before and after the excursion, the geomagnetic cutoff rigidity is clearly inversely correlated with geomagnetic latitude, nearly zero in both polar regions and about 13 GV in the equatorial region. However, the inverse correlation no longer exists at the midpoints of the excursion, geomagnetic cutoff rigidities are globally low at all latitudes, even having relatively higher values in the north polar region. A pronounced regional surface field intensity low analogous to the present-day South Atlantic Anomaly (SAA) structure was observed over 48-46 ka (Brown et al., 2018). This leads to a cutoff rigidity longitudinal minimum at around 270°E at 48.35 ka. The SAA-like structure is transient and disappeared before the Laschamps excursion started. In the middle of the Laschamps excursion, the cutoff rigidity has small longitudinal variations at about 1.5 GV globally, with some irregular R u increase at some longitudes. After the Laschamps excursion and before  The impact area illustrates the portion of the Earth where the energetic particles at a given cutoff rigidity can reach the atmosphere and is expressed as percentage. The impact area of the LSMOD.2 model and LSMOD.2 model up to SH degree 1 are shown in Figures 8c, 8f, and 8i. The impact area of the International Geomagnetic Reference Field (IGRF) model at epoch 2020 (Alken et al., 2021) is also shown for comparison. The impact areas for both the LSMOD.2 model and its dipole contribution only (SH degree 1) are mostly larger than the impact area today at most of the cutoff rigidity range, due to the lower geomagnetic dipole moment at that time compared to today's value. In the middle of the excursion, the impact area estimated from the LSMOD.2 model rises to 100% over 5 GV. The impact area for only the dipole part of the LSMOD.2 model rises to 100% to 1 GV, which again indicates that the multipole components of the paleomagnetic field become important for the shielding effect during the excursion (See also Movie S1i).

Particle Trajectory
The direct output of the trajectory tracing program is the trajectory of the traced particle starting near the Earth's surface. Here, we show several trajectories of the traced particle both before the excursion and at the excursion midpoint, which have quite different characteristics.
Here, we introduce a convenient equation for the rigidity-energy conversion (in SI units), and we also calculated the particle speed.
As the rigidity is defined as R = pc/Ze, the relativistic momentum can be represented as p = RZe/c. In Lorentz-covariant mechanics, the Lorentz factor of a particle with the relativistic momentum p and rest mass m 0 can be calculated as The particle kinetic energy is calculated as Thus, the rigidity-energy conversion equation is , the particle velocity is calculated as  (Figure 9b). As the proton kinetic energy increases (432 MeV ∼1.27 GeV), when the geomagnetic field is dipole-dominated, the particle follows complex loops before returning to the Earth. Tracing particles with 1.27-2.20 GeV kinetic energy, the particles from dipole-dominated geomagnetic field show more complex loops, some particles escape, and some particles return, which gives the penumbra with its series of allowed and forbidden bands. At higher kinetic energy (>2.2 GeV), the particles escape the geomagnetic field both at 48.35 ka and 40.95 ka (Figure 9e). Note that the trajectory depends on geomagnetic field configuration, starting location, and particle kinetic energy. The trajectory of traced protons starting at a different location (latitude 60°N, longitude 0°E) is shown in Figure  S5 of Supporting Information S1. At middle and high latitudes, the particle trajectory from a dipole-dominated field generally more easily develops complex loops (bouncing between both geomagnetic poles) than the particle trajectory during the Laschamps excursion. The shorter particle trajectory during the Laschamps excursion indicates that the geomagnetic field can trap less energetic particles, and the ring current and the radiation belt cannot stably exist (Tsareva, 2019;Tsareva et al., 2020;Vogt & Glassmeier, 2000).

Discussion and Conclusions
We presented global grids of the geomagnetic cutoff rigidities derived from the LSMOD.2 geomagnetic field model spanning 50-30 ka using trajectory tracing. For the first time, the cutoff rigidity during the Laschamps excursion is calculated from a realistic, data-based field model, in which the dipole moment is close to zero and the field is strongly non-dipole dominated at the excursion midpoint. We performed several tests to ensure the reliability of the cutoff rigidity calculation, such as analyzing the influence of the rigidity scanning step, investigating the incident particles from different directions, and a potential influence of the Earth's crustal field when the core field is weak. Both R vc (the vertical cutoff rigidity) and R ac (the apparent cutoff rigidity) are good estimations for the effective cutoff rigidity during the Laschamps excursion. More precisely, the difference between R vc and R ac is about 1 GV near the equator before and after the Laschamps excursion and is negligible during the excursion. At the midpoint of the Laschamps excursion, R ac is lower than 4 GV globally and nearly independent of latitude, which indicates that the exposure of Earth's atmosphere to energetic particles is high.
Theoretical equations have been widely used to estimate today's cutoff rigidity. We calculated both R cd (the cutoff rigidity considering only the dipole moment) and R cnd (the cutoff rigidity considering non-dipole field contributions at individual locations) with these equations using LSMOD.2. During the Laschamps excursion, the percentage difference between R ac and R cnd is high, which indicates that the simplified equations are not sufficient to estimate the cutoff rigidity during the excursion (42-40 ka). Note that R cnd gives the cutoff rigidity only considering the magnetic field at the studied location, but not the full morphology of the magnetic field through which the particle moves. Since the morphology of the geomagnetic field becomes complex during the excursion, the trajectory tracing provides a more realistic estimate of the cutoff rigidity compared to the theoretical equation.
The cutoff rigidity is a quantitative estimate of the Earth's geomagnetic field shielding effect. For today's geomagnetic field, the shielding effect is maximum at the geomagnetic equator and decreases with the increase of geomagnetic latitude. Above 60° geomagnetic latitude (called the latitude "knee"), the cosmic ray radiation intensity is nearly constant because the geomagnetic shielding is weaker than that of the atmosphere Smart & Shea, 2009). The cosmic radiation exposure above the latitude "knee" is approximately twice as much as the average cosmic radiation exposure at the geomagnetic equator. In the middle of the Laschamps excursion, the cutoff rigidity is globally lower than 4 GV (the cutoff rigidity at the latitude "knee") and nearly independent of latitude, which indicates that the cosmic radiation around 41 ka was globally approximately twice as much as the cosmic radiation at today's geomagnetic equator. During the excursion, most areas of the Earth' atmosphere were exposed to energetic particles of cosmic origin (Movie S1i).
The cosmogenic isotope (e.g., 10 Be, 14 C, or 36 Cl) production rate is controlled by the geomagnetic cutoff rigidity (Poluianov et al., 2016). The production rate of the cosmogenic isotopes is the basis of cosmogenic nuclides dating methods and sunspot number reconstruction (e.g., Solanki et al., 2004;Steinhilber et al., 2012). However, most of the cosmogenic isotope production rate calculation is based on an ideal dipole geomagnetic field with varying dipole moment (e.g., Gosse & Phillips, 2001;Lal, 1991). Lifton et al. (2014) and Lifton (2016) calculated the cutoff rigidity based on realistic geomagnetic field models, for example, CALS7k.2 and CALS3k.3 (Korte et al., , 2009), but the calculation is limited to Holocene time scale (14-0 ka). The cutoff rigidity results presented in this study provide the possibility of calculating the cosmogenic isotope production rate based on a realistic geomagnetic field model during the Laschampes excursion. Since the exposure of Earth's atmosphere to energetic particles of cosmic origin increased strongly during the excursion, the cosmogenic isotope production rate would also increase at that time (Black, 1967;Desilets & Zreda, 2003;Pigati & Lifton, 2004).
We used the internal magnetic field for calculating the cutoff rigidity. However, the external magnetic fields (e.g., magnetospheric magnetic field) would also influence the accuracy of the cutoff rigidity calculation at a low rigidity range (<2 GV). The magnetosphere is significantly altered during the excursion, and a paleomagnetosphere is formed (Saito et al., 1978;Siscoe & Chen, 1975). Although several paleomagnetosphere models have been presented (Stadelmann et al., 2010;Vogt & Glassmeier, 2001), none of them is based on a realistic internal field model like LSMOD.2, which is derived from paleomagnetic data. Moreover, when the geomagnetic field is no more dipole dominated and non-dipole components became important during an excursion, the magnetospheric current systems would change qualitatively (Vogt & Glassmeier, 2001;Wei et al., 2014). A realistic paleomagnetosphere model for the Laschamps excursion is thus not available. For today's magnetosphere, the magnetospheric field becomes important when tracing the particle with rigidity lower than 2 GV (Chu & Qin, 2016;Smart et al., 2000). Considering the Tsyganenko magnetospheric model (Tsyganenko & Usmanov, 1982) for different magnetic activity levels, Smart and Shea (2001) found that the cutoff rigidity is globally decreased during the magnetically active time. During geomagnetic polarity transitions, the influence of the magnetospheric field for different solar wind and geomagnetic conditions on the cutoff rigidity calculation needs to be further explored.
The accuracy of the cutoff rigidity calculation is limited by the accuracy of the geomagnetic field model (Shea & Smart, 1990). Although the LSMOD.2 model was developed to reconstruct the robust characteristics of the geomagnetic field during the Laschamps excursion, some limitations exist. The accuracy of the LSMOD.2, as a model based on paleomagnetic sediment data, is limited by the uncertainties in sediment chronologies, the quality of paleomagnetic sediment data, and the spatial and temporal coverage of the globe with paleomagnetic data. In principle, the SH model is a smoothed representation of the paleomagnetic data, and smaller-scale structures of the paleomagnetic field may not be well determined. Our future studies will compare cutoff rigidity results calculated from LSMOD.2 model and other paleomagnetic field models that cover the Laschamps and further geomagnetic excursions (e.g., the GGF100k model by Panovska et al. [2018]), to robustly resolve the cutoff rigidity structure in both spatial and temporal distribution during geomagnetic field transitions.