Where does the tracing of cosmic rays in real atmosphere terminate?

Majority of the human space activity remains confined to the low-Earth orbit (LEO). Cosmic rays penetrating the magnetosphere and overlying atmosphere create a complex radiation environment [1, 2]. For high-energy instruments operating in LEO such as AMS-02, theoretical models must be set based on whether an observed event belongs to the primary cosmic-ray population, i.e., if its rigidity exceeds the directional geomagnetic cut-off [3, 4] (the same technique is used for radiation-dose estimates). This classification is performed numerically by backtracing: a Runge-Kutta integrator propagates the charged-particle trajectory backward from the measured arrival position and direction under the Lorentz force [5]. If the trajectory escapes the magnetosphere and recedes to infinity, the particle can be flagged as primary (usually of galactic origin). The magnetic force is perpendicular to the velocity, and therefore, it does no work and the particle energy remains constant along the backward path.

The possible termination of a backtraced trajectory occurs near its perigee. Once the particle enters the dense lower atmosphere, it must interact with air nuclei and the backward propagation must stop. The key challenge is modeling such a termination. The present practice imposes a sharp boundary in altitude: 20 (radiation-exposure studies [6]), 40 (space-based cosmic-ray experiments [3]), or 100 km (Kármán line, recently adopted by AMS-02 for a conservative galactic-cosmic-ray selection [2, 7, 8]). However, these altitudes are ad hoc values. The actual cessation of a particle trajectory is dictated by physical processes, namely, microscopic interactions between the cosmic-ray particle and atmospheric nuclei.

Two microscopic mechanisms can terminate the backward extension of a trajectory. A single hard scattering that deflects the incident cosmic-ray particle out of the path elastically or inelastically, and a Bethe-Bloch energy loss, wherein the continuous slowdown violates the energy-conservation premise required for backtracing. Hard scattering is routinely studied at high energies through the air-shower signature recorded by ground-based indirect detectors such as the large high altitude air shower observatory [9], telescope array [10], and Pierre Auger observatory [11]. However, to the best of our knowledge, a quantitative treatment for backward tracking from LEO has not been presented; in such a treatment, the giga electronvolt range relevant for backtracing the necessary cross-section data must be extracted from accelerator measurements.

In this study, we illustrate the concept and estimate the order-of-magnitude effect by performing a controlled parametric study where the two physical processes are evaluated along the designed trajectories that neglect geomagnetic bending. Further, we derive a simple, factorized scaling result that can replace the ad-hoc sharp boundary when a binary stop/no-stop criterion is satisfied. Although the same processes can be integrated along the true curved trajectories produced by numerical backtracing and compared to the same criterion, we reserve this fully consistent treatment for a future publication.

Majority of the human space activity remains confined to the low-Earth orbit (LEO). Cosmic rays penetrating the magnetosphere and overlying atmosphere create a complex radiation environment [1, 2]. For high-energy instruments operating in LEO such as AMS-02, theoretical models must be set based on whether an observed event belongs to the primary cosmic-ray population, i.e., if its rigidity exceeds the directional geomagnetic cut-off [3, 4] (the same technique is used for radiation-dose estimates). This classification is performed numerically by backtracing: a Runge-Kutta integrator propagates the charged-particle trajectory backward from the measured arrival position and direction under the Lorentz force [5]. If the trajectory escapes the magnetosphere and recedes to infinity, the particle can be flagged as primary (usually of galactic origin). The magnetic force is perpendicular to the velocity, and therefore, it does no work and the particle energy remains constant along the backward path.

The possible termination of a backtraced trajectory occurs near its perigee. Once the particle enters the dense lower atmosphere, it must interact with air nuclei and the backward propagation must stop. The key challenge is modeling such a termination. The present practice imposes a sharp boundary in altitude: 20 (radiation-exposure studies [6]), 40 (space-based cosmic-ray experiments [3]), or 100 km (Kármán line, recently adopted by AMS-02 for a conservative galactic-cosmic-ray selection [2, 7, 8]). However, these altitudes are ad hoc values. The actual cessation of a particle trajectory is dictated by physical processes, namely, microscopic interactions between the cosmic-ray particle and atmospheric nuclei.

Two microscopic mechanisms can terminate the backward extension of a trajectory. A single hard scattering that deflects the incident cosmic-ray particle out of the path elastically or inelastically, and a Bethe-Bloch energy loss, wherein the continuous slowdown violates the energy-conservation premise required for backtracing. Hard scattering is routinely studied at high energies through the air-shower signature recorded by ground-based indirect detectors such as the large high altitude air shower observatory [9], telescope array [10], and Pierre Auger observatory [11]. However, to the best of our knowledge, a quantitative treatment for backward tracking from LEO has not been presented; in such a treatment, the giga electronvolt range relevant for backtracing the necessary cross-section data must be extracted from accelerator measurements.

In this study, we illustrate the concept and estimate the order-of-magnitude effect by performing a controlled parametric study where the two physical processes are evaluated along the designed trajectories that neglect geomagnetic bending. Further, we derive a simple, factorized scaling result that can replace the ad-hoc sharp boundary when a binary stop/no-stop criterion is satisfied. Although the same processes can be integrated along the true curved trajectories produced by numerical backtracing and compared to the same criterion, we reserve this fully consistent treatment for a future publication.

Majority of the human space activity remains confined to the low-Earth orbit (LEO). Cosmic rays penetrating the magnetosphere and overlying atmosphere create a complex radiation environment [1, 2]. For high-energy instruments operating in LEO such as AMS-02, theoretical models must be set based on whether an observed event belongs to the primary cosmic-ray population, i.e., if its rigidity exceeds the directional geomagnetic cut-off [3, 4] (the same technique is used for radiation-dose estimates). This classification is performed numerically by backtracing: a Runge-Kutta integrator propagates the charged-particle trajectory backward from the measured arrival position and direction under the Lorentz force [5]. If the trajectory escapes the magnetosphere and recedes to infinity, the particle can be flagged as primary (usually of galactic origin). The magnetic force is perpendicular to the velocity, and therefore, it does no work and the particle energy remains constant along the backward path.

The possible termination of a backtraced trajectory occurs near its perigee. Once the particle enters the dense lower atmosphere, it must interact with air nuclei and the backward propagation must stop. The key challenge is modeling such a termination. The present practice imposes a sharp boundary in altitude: 20 (radiation-exposure studies [6]), 40 (space-based cosmic-ray experiments [3]), or 100 km (Kármán line, recently adopted by AMS-02 for a conservative galactic-cosmic-ray selection [2, 7, 8]). However, these altitudes are ad hoc values. The actual cessation of a particle trajectory is dictated by physical processes, namely, microscopic interactions between the cosmic-ray particle and atmospheric nuclei.

Two microscopic mechanisms can terminate the backward extension of a trajectory. A single hard scattering that deflects the incident cosmic-ray particle out of the path elastically or inelastically, and a Bethe-Bloch energy loss, wherein the continuous slowdown violates the energy-conservation premise required for backtracing. Hard scattering is routinely studied at high energies through the air-shower signature recorded by ground-based indirect detectors such as the large high altitude air shower observatory [9], telescope array [10], and Pierre Auger observatory [11]. However, to the best of our knowledge, a quantitative treatment for backward tracking from LEO has not been presented; in such a treatment, the giga electronvolt range relevant for backtracing the necessary cross-section data must be extracted from accelerator measurements.

In this study, we illustrate the concept and estimate the order-of-magnitude effect by performing a controlled parametric study where the two physical processes are evaluated along the designed trajectories that neglect geomagnetic bending. Further, we derive a simple, factorized scaling result that can replace the ad-hoc sharp boundary when a binary stop/no-stop criterion is satisfied. Although the same processes can be integrated along the true curved trajectories produced by numerical backtracing and compared to the same criterion, we reserve this fully consistent treatment for a future publication.

Most human space activity is still confined to low-Earth orbit (LEO). Cosmic rays that penetrate the magnetosphere and the overlying atmosphere create a complex radiation environment [1, 2]. For high-energy instruments operating in LEO, such as AMS-02, theoretical models must decide whether an observed event belongs to the primary cosmic-ray population, i.e. whether its rigidity exceeds the directional geomagnetic cutoff [3, 4] (the same technique is used for radiation-dose estimates). This classification is performed numerically by backtracing: a Runge-Kutta integrator propagates the charged-particle trajectory backward from the measured arrival position and direction under the Lorentz force alone [5]. If the trajectory ultimately escapes the magnetosphere and recedes to infinity, the particle is flagged as primary (usually of galactic origin). Because the magnetic force is always perpendicular to the velocity, it does no work and the particle energy remains constant along the backward path.

Possible termination of a backtraced trajectory occurs near its perigee. Once the particle enters the dense lower atmosphere it interacts with air nuclei and the backward propagation must stop. The key question is how to exactly model such termination. Present practice imposes a sharp boundary in altitude: 20 km (radiation-exposure studies [6]), 40 km (space-based cosmic-ray experiments [3]), or 100 km (the Kármán line, recently adopted by AMS-02 for a conservative galactic-cosmic-ray selection [2, 7, 8]). These altitudes are, however, ad hoc. The actual cessation of a particle's trajectory should be dictated by physical processes, namely the microscopic interactions between the cosmic-ray particle and atmospheric nuclei.

Two microscopic mechanisms can terminate the backward extension of a trajectory. The first is a single hard scattering that deflects the incident cosmic-ray particle out of the path elastically or inelastically. The second is Bethe-Bloch energy loss: the continuous slowdown violates the energy-conservation premise on which backtracing relies. Hard scattering is routinely studied at high energies through the air-shower signature recorded by ground-based indirect detectors such as the Large High Altitude Air Shower Observatory [9], the Telescope Array [10] and the Pierre Auger Observatory [11]. However, to the best of our knowledge, no quantitative treatment suited to backward tracking from LEO has yet been presented, in which at the GeV range relevant for backtracing the necessary cross-section data must be extracted from accelerator measurements.

To illustrate the concept and estimate the order-of-magnitude effect, we perform a controlled parametric study in which the two physical processes are evaluated along designed trajectories that neglect geomagnetic bending. A simple, factorized scaling result is derived that can readily replace the ad-hoc sharp boundary whenever a binary stop/no-stop criterion is sufficient. In parallel, the same processes can be integrated along the true curved trajectories produced by numerical backtracing and compared to the same criterion; we reserve this fully consistent treatment for a future publication.

Models. Bethe-Bloch energy loss [12–14] arises from the cumulative effect of a continuous series of binary Coulomb scatterings. The stopping power is

where e is the elementary charge, $ m_{e} $ the electron mass, and u the unified atomic mass unit. The projectile velocity in units of c is $ \beta=v/c $ (we adopt natural units $ \hbar=1 $ and $ c=1 $), and $ Z_{\text{P}} $ is its atomic number. For each target element T, $ \rho_{\text{T}}(h) $ is the mass density contributed by that element at altitude h, and $ Z_{\text{T}} $, $ A_{\text{T}} $ are its atomic and mass numbers. $ I_{\text{T,min}} $ is the mean excitation energy that sets the maximum impact parameter for individual collisions; we use the values given in Table 1.1 of Ref. [13]. The maximum transferable energy, corresponding to the minimum impact parameter, is

with $ m_{\text{P}} $ the projectile mass.

In a magnetic field the local gyroradius is $ r=\mathfrak{R}/B $, where the rigidity $ \mathfrak{R} = p/q = p/(Z_\text{P}e) $ encodes both particle species and kinematics and therefore completely determines the trajectory in backtracing. Using $ E^2=p^2+m^2 $ ($ EdE=pdp $) and $ \beta=p/E $, the integrated rigidity loss is

Throughout the numerical work we treat β as constant rather than decreasing along the path ¹, i.e. the slowing-down is included only perturbatively. Then the fractional rigidity loss reads

This crude estimate clearly fails when $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim 1 $. However, in the following we will see that the eventual interested $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $ can be obtained by multiplying such results with additional small factors from perigee altitude and curvature radius, in such cases this perturbative calculation holds conveniently. Note that the integrand depends on the projectile only through $ m_{\text{P}} $ inside the logarithm and via $ \beta=\mathfrak{R}Ze/((\mathfrak{R}Ze)^2+m_{\text{P}}^2)^{1/2} $. If we evaluate the integral once for protons and neglect the weak $ m_{\text{P}} $ dependence, the result for any ion scales as $ Z_\text{P}^2/(E\beta^4) $. In this approximation a single backtracing performed for protons suffices; the fractional rigidity shift for other species is obtained by scaling.

For single hard scattering we use the Glauber-Gribov formalism to obtain proton- and nucleus-nucleus cross sections from the underlying nucleon-nucleon data, adopting the simplifications of Ref. [15]. The total (subscript "tot" omitted) hadron-nucleus and nucleus-nucleus cross sections are [15]

with the nuclear radii

where $ r_{0}=1.1 $ fm. The nucleon-nucleon level cross section is weighted from the proton-proton, neutron-neutron and proton-neutron ² cross section

and for proton-nucleon scattering we set $ Z_{\rm{P}}=A_{\rm{P}}=1 $.

The nucleon-nucleon cross sections used in this work are shown in Fig. 1. Proton-proton and proton-neutron data are taken directly from the 2022 PDG compilation [16]; when several measurements exist at the same beam momentum we average them, and we apply a Savitzky-Golay filter [17] (window length 5, polynomial order 1) to reduce residual scatter. The lowest momentum measured is $ \sim0.2 $ GeV. We do not extrapolate below this value because, as shown later, hard scattering is already negligible compared with Bethe-Bloch energy loss in that regime.

Figure 1.  (color online) Left: Proton-proton total cross section from the 2022 PDG compilation (red points; averaged when multiple measurements exist at the same beam momentum) together with the Savitzky-Golay smoothed curve (blue) adopted in this work. Middle: Same for the proton-neutron total cross section. Right: Neutron-neutron total cross section (blue curve) obtained from a Savitzky-Golay fit to three data sets: low-momentum points (orange dashed) extracted as $ \sigma_\text{nd}-\sigma_\text{np} $, high-momentum points (purple) assumed equal to $ \sigma_\text{pp} $, and intermediate-momentum data (red) taken directly from Ref. [18].

Free neutron-neutron scattering cannot be measured directly. Following standard practice [18, 19],$ \sigma_\text{nn} $ is extracted by subtracting the neutron-proton cross section from neutron-deuteron data. For momenta $ \lesssim0.8 $ GeV multiple scattering inside the weakly-bound deuteron is small and the simple subtraction is accurate [20]. At high momenta ($ \gtrsim3 $ GeV) isospin symmetry implies $ \sigma_\text{nn}\approx\sigma_\text{pp} $ up to tiny electromagnetic corrections, so we adopt the proton-proton values there. The intermediate region (0.8-3 GeV) is interpolated from the dedicated analysis of Ref. [18], with corrections from Refs. [21–23]. This piece-wise $ \sigma_\text{nn} $ is displayed in the right-hand panel of Fig. 1.

Eventually the expected scattering number along the path is

We adopt the US Standard Atmosphere 1976 [24] for the atmospheric density profile. Being independent of geographic location and season, it yields universal results; unlike more recent MSISE00 model [25]. The carbon dioxide level has risen from 0.0314% in 1976 to about 0.0425% in 2025 [26]; we have updated the atmospheric composition accordingly. Because the microscopic interactions treated below are atomic, molecular species are first decomposed into their constituent atoms.

Figure 2 shows the partial densities of the updated US Standard Atmosphere 1976. Below 86 km the atmosphere is taken to be homogeneously mixed, so the relative abundances of nitrogen, oxygen, argon, carbon, helium and hydrogen are constant. Throughout this region each component follows an exponential profile, $ \rho_ {\rm{T}}\propto e^{-0.14h/ {\rm{km}}} $, to excellent accuracy.

Figure 2.  (color online) Partial densities $ \rho_ {\rm{T}} $ of the six dominant species-hydrogen (red), helium (orange), carbon (green), nitrogen (cyan), oxygen (blue) and argon (violet)-as functions of altitude in the updated US Standard Atmosphere 1976 model. The simple exponential fit $ \rho_ {\rm{T}}\propto e^{-0.14h/ {\rm{km}}} $ (grey dashed) is shown for comparison.

Results. Effect of perigee altitude. The exponential fall-off of atmospheric density guarantees that both $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ receive their dominant contributions from the segment closest to Earth. The natural control parameter is therefore the perigee altitude, h.

As a first baseline we consider a straight-line trajectory that is tangent to the atmosphere at altitude h (impact parameter $ b=r_{\oplus}+h $, with $ r_{\oplus}=6371.2 $ km); this is the red curve in Fig. 3. Although this path ignores geomagnetic bending, it captures the correct order of magnitude of the two stopping mechanisms.

Figure 3.  (color online) The two trajectory geometries considered: a straight line tangent to the atmosphere (red) and an upward-bent trajectory with radius of curvature r (blue).

Figure 4 shows the integrated fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), the expected hard-scattering count $ \langle N\rangle $ (middle) and their sum (right) along an infinite straight-line trajectory tangent to the atmosphere for protons of five rigidities. Although $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ describe distinct physical mechanisms, both are dimensionless measures of the breakdown of energy-conserving backtracing. We heuristically sum them to define a unified criterion for trajectory termination. All quantities exhibit the same $ \propto e^{-0.14h/ {\rm{km}}} $ dependence on perigee altitude h and can therefore be factorized out. The altitude range of 30-100 km is chosen to highlight the variation around the critical perigee altitude h where interactions become significant, which is much lower than the observer's altitude (e.g., $ \sim400 $ km for AMS-02).

Figure 4.  (color online) The cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), expected hard-scattering count $ \langle N\rangle $ (middle) and their sum (right) versus perigee altitude h for a straight-line tangent trajectory. Proton rigidities: 0.2 GV (red), 0.5 GV (orange), 1 GV (green), 2 GV (blue) and 5 GV (black). The altitude dependence is well described by the common exponential factor $ \exp(-0.14\,h/\text{km}) $ (grey dashed) and can therefore be factorized out.

Effect of perigee curvature radius. The straight-line tangent ignores bending by the Lorentz force; realistic trajectories must be obtained by (back-)tracing. Near perigee, however, any path is characterised by its local radius of curvature $ r=\mathfrak{R}/B $, almost independent of the remainder of the orbit. We therefore replace the tangent line by an upward-bending circular arc of radius r that reaches the same perigee altitude h (blue curve in Fig. 3). Downward-bending geometries with $ r>r_{\oplus}+h $ can also yield perigee h, but are not considered here. The limit $ r\to\infty $ recovers the straight-line case. The parameterized integrated trajectory is indeed a circle with radius r, while we have checked that the arc around the perigee almost contributes the whole integration. For simplicity the arc is taken to lie in the plane containing the Earth's centre and the perigee point.

Figure 5 compares the numerical integration of $ \langle N \rangle $ with its asymptotic and fitted forms. The behaviour of $ \Delta\mathfrak{R}/\mathfrak{R} $ is identical. The global fit

Figure 5.  (color online) The expected hard-scattering count $ \langle N \rangle $ versus local radius of curvature r for an upward-bent trajectory with perigee altitude 50 km. The black curve gives the full calculation; grey dashed curves show the low-r fit and high-r asymptote; the red curve is the global fit.

normalised to its $ r\to\infty $ limit, reproduces the calculation quite well. Combining this curvature factor with the altitude dependence gives the compact factorized forms, which works quite well if normalized with its $ r\to\infty $ asymptotic value. With this fitting for local curvature radius r around perigee and the above fitting for perigee altitude h, the two dimensionless variables can be factorized as

Particle type and rigidity dependence. Figure 6 shows the sum $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ versus rigidity for five cosmic-ray species on a straight-line trajectory tangent to the atmosphere at $ h=50 $ km. Bethe-Bloch energy loss, scaling approximately as $ \beta^{-2} $, dominates at low rigidities, while hard scattering, with its nearly energy-independent cross section, takes over at high rigidities. The crossover occurs at $ \sim0.57 $ GV for protons and rises to $ \sim1.4 $ GV for iron. These reference values are inserted into Eq. 13 to obtain the full altitude- and curvature-corrected estimates. Again, values $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim1 $ obtained from the perturbative calculation are simply rescaled by Eq. 13 to the region $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $.

Figure 6.  (color online) The cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (dotted), expected hard-scattering count $ \langle N \rangle $ (dashed) and their sum (solid) versus incident rigidity for five elements-proton (black), helium (red), oxygen (orange), silicon (green) and iron (blue)-in a straight-line trajectory with perigee altitude 50 km.

For protons at 50 km, $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle\gtrsim1 $ at all rigidities in the straight-line tangent geometry, so this altitude may serve as a reasonable lower boundary when a sharp boundary is imposed in backtracing. The chosen threshold of $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ should be smaller yet close to unity. Note, however, that the values in Fig. 6 neglect the finite-radius-of-curvature correction. For iron the sum is roughly ten times larger, so the same criterion requires the altitude-curvature factor to be $ \sim0.1 $, equivalent to raising the tangent altitude by $ \sim16 $ km.

Discussion. The physics-based termination criterion improves ad hoc backtracing boundaries. For a 1 GV proton in the straight-line tangent geometry we obtain $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle=0.93 $ at 50 km, but 3.4 at 40 km. Primary charged cosmic rays therefore cannot originate from a perigee just above the formerly adopted 40 km cutoff [3].

As a second illustration we again consider vertical incidence at $ 0^\circ $ longitude, $ 0^\circ $ latitude and 400 km altitude. Using the conventional algorithm of Ref. [2] (IGRF-13 field for epoch 2023 with the traditional 100 km sharp boundary) the geomagnetic cutoff rigidity is 12.04 GV. A very dedicated search shows that for $ \mathfrak{R}=12.03809986288667 $ GV the backtraced trajectory reaches a perigee of 98.3 km before escaping to infinity. The old sharp-boundary prescription would discard this event, whereas the physics-based criterion developed here accepts it.

Because the physics-based criterion $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ is intrinsically probabilistic, the very notion of a penumbra becomes questionable. Penumbra denotes the narrow rigidity interval - usually around the vertical cutoff - where numerical backtracing switches repeatedly between "allowed" and "forbidden" as the rigidity is scanned in fine steps. Bethe-Bloch energy loss is deterministic, but hard scattering supplies only an expected number of collisions; even when the path integral gives $ \langle N\rangle=1 $ there remains a finite survival probability. Thus the boundary between allowed and forbidden rigidities is unavoidably smeared, irrespective of whether the backtracing boundary is set at the physics-based 50 km, or any other fixed altitude.

The allowed cone (i.e., the set of permitted incident directions for cosmic rays in LEO) can be approximated at high rigidities by the solid angle above the line of sight tangent to the sharp boundary [2]. For an arrival altitude of 400 km, the use of a 20 km boundary yields an allowed cone that is 1.0% larger for protons (1.5% for iron) compared to a 50 km (66 km for iron) boundary. Conversely, a 100 km boundary yields a cone that is 1.7% smaller for protons (1.2% for iron). Consequently, space dosimetry [6] is overestimated by approximately 1.0%-1.5%. Furthermore, the AMS-02 experiment can utilize this additional 1.2%-1.7% of events vaively, which corresponds to billions of events given its 15 years of data collection.

A practical limitation of the present study is that the local radius of curvature around perigee is known only after a full backtracing calculation. In such a code the integrals for $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ can be evaluated event-by-event, without recourse to the parametric factor $ g(r) $ in Eq. 11, yielding the exact energy-loss and scattering probability for each trajectory and permitting any desired numerical threshold to be imposed a posteriori. This refinement in backtracing removes the uncertainties introduced by the sharp-boundary altitude approximation, enabling a single backtracing calculation to suffice for all cases.

Bethe-Bloch energy loss [12–14] can be attributed to the cumulative effect of a continuous series of binary Coulomb scatterings. The stopping power is expressed as

where e, $ m_{e} $, and u represent the elementary charge, electron mass, and unified atomic mass unit, respectively. The projectile velocity in units of c is $ \beta=v/c $ (we adopt natural units $ \hbar=1 $ and $ c=1 $), and $ Z_{\text{P}} $ represents its atomic number. For each target element T, $ \rho_{\text{T}}(h) $ represents the mass density contributed by that element at altitude h, and $ Z_{\text{T}} $ and $ A_{\text{T}} $ are its atomic and mass numbers. Further, $ I_{\text{T,min}} $ represents the mean excitation energy that sets the maximum impact parameter for individual collisions; we use the values given in Table 1.1 of Ref. [13]. The maximum transferable energy corresponding to the minimum impact parameter is given as

$ m_{\text{P}} $ is the projectile mass.

In a magnetic field, the local gyroradius is $ r=\mathfrak{R}/B $, where the rigidity $ \mathfrak{R} = p/q = p/(Z_\text{P}e) $ encodes both particle species and kinematics and determines the trajectory in backtracing. Using $ E^2=p^2+m^2 $ ($ E{\rm d}E=p{\rm d}p $) and $ \beta=p/E $, the integrated rigidity loss is given as

Throughout the numerical work we treat β as constant rather than decreasing along the path ¹, i.e. the slowing-down is included only perturbatively. Then the fractional rigidity loss reads

This crude estimate fails when $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim 1 $; however, we see that the eventual interested $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $ can be obtained by multiplying such results with additional small factors from perigee altitude and curvature radius. In such cases, this perturbative calculation holds conveniently. The integrand depends on the projectile only through $ m_{\text{P}} $ inside the logarithm and via $ \beta=\mathfrak{R}Ze/((\mathfrak{R}Ze)^2+m_{\text{P}}^2)^{1/2} $. If we evaluate the integral once for protons and neglect the weak $ m_{\text{P}} $ dependence, the result for any ion scales as $ Z_\text{P}^2/(E\beta^4) $. In this approximation, a single backtracing performed for protons suffices, and the fractional rigidity shift for other species is obtained by scaling.

For single hard scattering, we use the Glauber-Gribov formalism to obtain proton- and nucleus-nucleus cross-sections from the underlying nucleon-nucleon data, adopting the simplifications of Ref. [15]. The total (subscript "tot" omitted) hadron-nucleus and nucleus-nucleus cross-sections are [15]

with the nuclear radii

where $ r_{0}=1.1 $ fm. The nucleon-nucleon level cross section is weighted from the proton-proton, neutron-neutron and proton-neutron ² cross section

For proton-nucleon scattering, we set $ Z_{\rm{P}}=A_{\rm{P}}=1 $.

The nucleon-nucleon cross-sections used in this work are presented in Fig. 1. Proton-proton and proton-neutron data are obtained directly from the 2022 PDG compilation [16]. When several measurements exist at the same beam momentum, we average them and apply a Savitzky-Golay filter [17] (window length 5, polynomial order 1) to reduce residual scatter. The lowest measured momentum is $ \sim0.2 $ GeV. We do not extrapolate below this value because hard scattering is already negligible compared with the Bethe-Bloch energy loss in that regime.

Figure 1.  (color online) Left: Proton-proton total cross-section from the 2022 PDG compilation (red points; averaged when multiple measurements exist at the same beam momentum) with the Savitzky-Golay smoothed curve (blue) adopted in this work. Middle: Same for the proton-neutron total cross-section. Right: Neutron-neutron total cross-section (blue curve) obtained from a Savitzky-Golay fit to three data sets: low-momentum points (orange dashed) extracted as $ \sigma_\text{nd}-\sigma_\text{np} $, high-momentum points (purple) assumed equal to $ \sigma_\text{pp} $, and intermediate-momentum data (red) taken directly from Ref. [18].

Free neutron-neutron scattering cannot be measured directly. Following standard practice [18, 19], $ \sigma_\text{nn} $ is extracted by subtracting the neutron-proton cross-section from neutron-deuteron data. For momenta $ \lesssim0.8 $ GeV, multiple scattering inside the weakly-bound deuteron is small and the simple subtraction is accurate [20]. At high momenta ($ \gtrsim3 $ GeV), isospin symmetry implies $ \sigma_\text{nn}\approx\sigma_\text{pp} $ up to tiny electromagnetic corrections, and therefore, we adopt the proton-proton values. The intermediate region (0.8-3 GeV) is interpolated from the dedicated analysis of Ref. [18], with corrections from Refs. [21–23]. This piece-wise $ \sigma_\text{nn} $ is displayed in the right-hand panel of Fig. 1.

Eventually, the expected scattering number along the path is

We adopt the US Standard Atmosphere 1976 [24] for the atmospheric density profile. This profile is independent of the geographic location and season and yields universal results, unlike the more recent MSISE00 model [25]. The carbon dioxide level has risen from 0.0314% in 1976 to about 0.0425% in 2025 [26], and thus, we updated the atmospheric composition accordingly. The microscopic interactions treated below are atomic, and therefore, molecular species are decomposed into their constituent atoms.

Figure 2 shows the partial densities of the updated US Standard Atmosphere 1976. Below 86 km, the atmosphere is considered homogeneously mixed. Thus, the relative abundancies of nitrogen, oxygen, argon, carbon, helium, and hydrogen reamin constant. Throughout this region, each component follows an exponential profile, $ \rho_ {\rm{T}}\propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $, realizing excellent accuracy.

Figure 2.  (color online) Partial densities $ \rho_ {\rm{T}} $ of the six dominant species, i.e., hydrogen (red), helium (orange), carbon (green), nitrogen (cyan), oxygen (blue), and argon (violet), as functions of altitude in the updated US Standard Atmosphere 1976 model. The simple exponential fit $ \rho_ {\rm{T}}\propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $ (grey dashed) is shown for comparison.

Bethe-Bloch energy loss [12–14] can be attributed to the cumulative effect of a continuous series of binary Coulomb scatterings. The stopping power is expressed as

where e, $ m_{e} $, and u represent the elementary charge, electron mass, and unified atomic mass unit, respectively. The projectile velocity in units of c is $ \beta=v/c $ (we adopt natural units $ \hbar=1 $ and $ c=1 $), and $ Z_{\text{P}} $ represents its atomic number. For each target element T, $ \rho_{\text{T}}(h) $ represents the mass density contributed by that element at altitude h, and $ Z_{\text{T}} $ and $ A_{\text{T}} $ are its atomic and mass numbers. Further, $ I_{\text{T,min}} $ represents the mean excitation energy that sets the maximum impact parameter for individual collisions; we use the values given in Table 1.1 of Ref. [13]. The maximum transferable energy corresponding to the minimum impact parameter is given as

$ m_{\text{P}} $ is the projectile mass.

In a magnetic field, the local gyroradius is $ r=\mathfrak{R}/B $, where the rigidity $ \mathfrak{R} = p/q = p/(Z_\text{P}e) $ encodes both particle species and kinematics and determines the trajectory in backtracing. Using $ E^2=p^2+m^2 $ ($ E{\rm d}E=p{\rm d}p $) and $ \beta=p/E $, the integrated rigidity loss is given as

Throughout the numerical work we treat β as constant rather than decreasing along the path ¹, i.e. the slowing-down is included only perturbatively. Then the fractional rigidity loss reads

This crude estimate fails when $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim 1 $; however, we see that the eventual interested $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $ can be obtained by multiplying such results with additional small factors from perigee altitude and curvature radius. In such cases, this perturbative calculation holds conveniently. The integrand depends on the projectile only through $ m_{\text{P}} $ inside the logarithm and via $ \beta=\mathfrak{R}Ze/((\mathfrak{R}Ze)^2+m_{\text{P}}^2)^{1/2} $. If we evaluate the integral once for protons and neglect the weak $ m_{\text{P}} $ dependence, the result for any ion scales as $ Z_\text{P}^2/(E\beta^4) $. In this approximation, a single backtracing performed for protons suffices, and the fractional rigidity shift for other species is obtained by scaling.

For single hard scattering, we use the Glauber-Gribov formalism to obtain proton- and nucleus-nucleus cross-sections from the underlying nucleon-nucleon data, adopting the simplifications of Ref. [15]. The total (subscript "tot" omitted) hadron-nucleus and nucleus-nucleus cross-sections are [15]

with the nuclear radii

where $ r_{0}=1.1 $ fm. The nucleon-nucleon level cross section is weighted from the proton-proton, neutron-neutron and proton-neutron ² cross section

For proton-nucleon scattering, we set $ Z_{\rm{P}}=A_{\rm{P}}=1 $.

The nucleon-nucleon cross-sections used in this work are presented in Fig. 1. Proton-proton and proton-neutron data are obtained directly from the 2022 PDG compilation [16]. When several measurements exist at the same beam momentum, we average them and apply a Savitzky-Golay filter [17] (window length 5, polynomial order 1) to reduce residual scatter. The lowest measured momentum is $ \sim0.2 $ GeV. We do not extrapolate below this value because hard scattering is already negligible compared with the Bethe-Bloch energy loss in that regime.

Figure 1.  (color online) Left: Proton-proton total cross-section from the 2022 PDG compilation (red points; averaged when multiple measurements exist at the same beam momentum) with the Savitzky-Golay smoothed curve (blue) adopted in this work. Middle: Same for the proton-neutron total cross-section. Right: Neutron-neutron total cross-section (blue curve) obtained from a Savitzky-Golay fit to three data sets: low-momentum points (orange dashed) extracted as $ \sigma_\text{nd}-\sigma_\text{np} $, high-momentum points (purple) assumed equal to $ \sigma_\text{pp} $, and intermediate-momentum data (red) taken directly from Ref. [18].

Free neutron-neutron scattering cannot be measured directly. Following standard practice [18, 19], $ \sigma_\text{nn} $ is extracted by subtracting the neutron-proton cross-section from neutron-deuteron data. For momenta $ \lesssim0.8 $ GeV, multiple scattering inside the weakly-bound deuteron is small and the simple subtraction is accurate [20]. At high momenta ($ \gtrsim3 $ GeV), isospin symmetry implies $ \sigma_\text{nn}\approx\sigma_\text{pp} $ up to tiny electromagnetic corrections, and therefore, we adopt the proton-proton values. The intermediate region (0.8-3 GeV) is interpolated from the dedicated analysis of Ref. [18], with corrections from Refs. [21–23]. This piece-wise $ \sigma_\text{nn} $ is displayed in the right-hand panel of Fig. 1.

Eventually, the expected scattering number along the path is

We adopt the US Standard Atmosphere 1976 [24] for the atmospheric density profile. This profile is independent of the geographic location and season and yields universal results, unlike the more recent MSISE00 model [25]. The carbon dioxide level has risen from 0.0314% in 1976 to about 0.0425% in 2025 [26], and thus, we updated the atmospheric composition accordingly. The microscopic interactions treated below are atomic, and therefore, molecular species are decomposed into their constituent atoms.

Figure 2 shows the partial densities of the updated US Standard Atmosphere 1976. Below 86 km, the atmosphere is considered homogeneously mixed. Thus, the relative abundancies of nitrogen, oxygen, argon, carbon, helium, and hydrogen reamin constant. Throughout this region, each component follows an exponential profile, $ \rho_ {\rm{T}}\propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $, realizing excellent accuracy.

Figure 2.  (color online) Partial densities $ \rho_ {\rm{T}} $ of the six dominant species, i.e., hydrogen (red), helium (orange), carbon (green), nitrogen (cyan), oxygen (blue), and argon (violet), as functions of altitude in the updated US Standard Atmosphere 1976 model. The simple exponential fit $ \rho_ {\rm{T}}\propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $ (grey dashed) is shown for comparison.

Bethe-Bloch energy loss [12–14] can be attributed to the cumulative effect of a continuous series of binary Coulomb scatterings. The stopping power is expressed as

where e, $ m_{e} $, and u represent the elementary charge, electron mass, and unified atomic mass unit, respectively. The projectile velocity in units of c is $ \beta=v/c $ (we adopt natural units $ \hbar=1 $ and $ c=1 $), and $ Z_{\text{P}} $ represents its atomic number. For each target element T, $ \rho_{\text{T}}(h) $ represents the mass density contributed by that element at altitude h, and $ Z_{\text{T}} $ and $ A_{\text{T}} $ are its atomic and mass numbers. Further, $ I_{\text{T,min}} $ represents the mean excitation energy that sets the maximum impact parameter for individual collisions; we use the values given in Table 1.1 of Ref. [13]. The maximum transferable energy corresponding to the minimum impact parameter is given as

$ m_{\text{P}} $ is the projectile mass.

In a magnetic field, the local gyroradius is $ r=\mathfrak{R}/B $, where the rigidity $ \mathfrak{R} = p/q = p/(Z_\text{P}e) $ encodes both particle species and kinematics and determines the trajectory in backtracing. Using $ E^2=p^2+m^2 $ ($ E{\rm d}E=p{\rm d}p $) and $ \beta=p/E $, the integrated rigidity loss is given as

Throughout the numerical work we treat β as constant rather than decreasing along the path ¹, i.e. the slowing-down is included only perturbatively. Then the fractional rigidity loss reads

This crude estimate fails when $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim 1 $; however, we see that the eventual interested $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $ can be obtained by multiplying such results with additional small factors from perigee altitude and curvature radius. In such cases, this perturbative calculation holds conveniently. The integrand depends on the projectile only through $ m_{\text{P}} $ inside the logarithm and via $ \beta=\mathfrak{R}Ze/((\mathfrak{R}Ze)^2+m_{\text{P}}^2)^{1/2} $. If we evaluate the integral once for protons and neglect the weak $ m_{\text{P}} $ dependence, the result for any ion scales as $ Z_\text{P}^2/(E\beta^4) $. In this approximation, a single backtracing performed for protons suffices, and the fractional rigidity shift for other species is obtained by scaling.

For single hard scattering, we use the Glauber-Gribov formalism to obtain proton- and nucleus-nucleus cross-sections from the underlying nucleon-nucleon data, adopting the simplifications of Ref. [15]. The total (subscript "tot" omitted) hadron-nucleus and nucleus-nucleus cross-sections are [15]

with the nuclear radii

where $ r_{0}=1.1 $ fm. The nucleon-nucleon level cross section is weighted from the proton-proton, neutron-neutron and proton-neutron ² cross section

For proton-nucleon scattering, we set $ Z_{\rm{P}}=A_{\rm{P}}=1 $.

The nucleon-nucleon cross-sections used in this work are presented in Fig. 1. Proton-proton and proton-neutron data are obtained directly from the 2022 PDG compilation [16]. When several measurements exist at the same beam momentum, we average them and apply a Savitzky-Golay filter [17] (window length 5, polynomial order 1) to reduce residual scatter. The lowest measured momentum is $ \sim0.2 $ GeV. We do not extrapolate below this value because hard scattering is already negligible compared with the Bethe-Bloch energy loss in that regime.

Figure 1.  (color online) Left: Proton-proton total cross-section from the 2022 PDG compilation (red points; averaged when multiple measurements exist at the same beam momentum) with the Savitzky-Golay smoothed curve (blue) adopted in this work. Middle: Same for the proton-neutron total cross-section. Right: Neutron-neutron total cross-section (blue curve) obtained from a Savitzky-Golay fit to three data sets: low-momentum points (orange dashed) extracted as $ \sigma_\text{nd}-\sigma_\text{np} $, high-momentum points (purple) assumed equal to $ \sigma_\text{pp} $, and intermediate-momentum data (red) taken directly from Ref. [18].

Free neutron-neutron scattering cannot be measured directly. Following standard practice [18, 19], $ \sigma_\text{nn} $ is extracted by subtracting the neutron-proton cross-section from neutron-deuteron data. For momenta $ \lesssim0.8 $ GeV, multiple scattering inside the weakly-bound deuteron is small and the simple subtraction is accurate [20]. At high momenta ($ \gtrsim3 $ GeV), isospin symmetry implies $ \sigma_\text{nn}\approx\sigma_\text{pp} $ up to tiny electromagnetic corrections, and therefore, we adopt the proton-proton values. The intermediate region (0.8-3 GeV) is interpolated from the dedicated analysis of Ref. [18], with corrections from Refs. [21–23]. This piece-wise $ \sigma_\text{nn} $ is displayed in the right-hand panel of Fig. 1.

Eventually, the expected scattering number along the path is

We adopt the US Standard Atmosphere 1976 [24] for the atmospheric density profile. This profile is independent of the geographic location and season and yields universal results, unlike the more recent MSISE00 model [25]. The carbon dioxide level has risen from 0.0314% in 1976 to about 0.0425% in 2025 [26], and thus, we updated the atmospheric composition accordingly. The microscopic interactions treated below are atomic, and therefore, molecular species are decomposed into their constituent atoms.

Figure 2 shows the partial densities of the updated US Standard Atmosphere 1976. Below 86 km, the atmosphere is considered homogeneously mixed. Thus, the relative abundancies of nitrogen, oxygen, argon, carbon, helium, and hydrogen reamin constant. Throughout this region, each component follows an exponential profile, $ \rho_ {\rm{T}}\propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $, realizing excellent accuracy.

Figure 2.  (color online) Partial densities $ \rho_ {\rm{T}} $ of the six dominant species, i.e., hydrogen (red), helium (orange), carbon (green), nitrogen (cyan), oxygen (blue), and argon (violet), as functions of altitude in the updated US Standard Atmosphere 1976 model. The simple exponential fit $ \rho_ {\rm{T}}\propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $ (grey dashed) is shown for comparison.

The exponential fall-off of atmospheric density guarantees that both $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ receive their dominant contributions from the segment closest to Earth. The natural control parameter is therefore the perigee altitude h.

As a first baseline, we consider a straight-line trajectory tangential to the atmosphere at altitude h (impact parameter $ b=r_{\oplus}+h $, with $ r_{\oplus}=6371.2 $ km); this is shown by the red curve in Fig. 3. Although this path ignores geomagnetic bending, it captures the correct order of magnitude of the two stopping mechanisms.

Figure 3.  (color online) Two trajectory geometries considered: a straight line tangential to the atmosphere (red) and an upward-bent trajectory with a radius of curvature r (blue).

Figure 4 shows the integrated fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), expected hard-scattering count $ \langle N\rangle $ (middle), and their sum (right) along an infinite straight-line trajectory tangential to the atmosphere for protons of five rigidities. Although $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ describe distinct physical mechanisms, both are dimensionless measures of the breakdown of energy-conserving backtracing. We heuristically add them to define a unified criterion for trajectory termination. All quantities exhibit the same $ \propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $ dependence on perigee altitude h, and therefore, they can be factorized out. The altitude range of 30-100 km is selected to highlight the variation around the critical perigee altitude h, where interactions become significant. This is considerably lower than the altitude of the observer (e.g., $ \sim400 $ km for AMS-02).

Figure 4.  (color online) Cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), expected hard-scattering count $ \langle N\rangle $ (middle), and their sum (right) versus perigee altitude h for a straight-line tangent trajectory. Proton rigidities: 0.2 GV (red), 0.5 GV (orange), 1 GV (green), 2 GV (blue), and 5 GV (black). The altitude dependence is well described by the common exponential factor $ \exp(-0.14\,{\rm h}/\text{km}) $ (grey dashed) and can therefore be factorized out.

The exponential fall-off of atmospheric density guarantees that both $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ receive their dominant contributions from the segment closest to Earth. The natural control parameter is therefore the perigee altitude h.

As a first baseline, we consider a straight-line trajectory tangential to the atmosphere at altitude h (impact parameter $ b=r_{\oplus}+h $, with $ r_{\oplus}=6371.2 $ km); this is shown by the red curve in Fig. 3. Although this path ignores geomagnetic bending, it captures the correct order of magnitude of the two stopping mechanisms.

Figure 3.  (color online) Two trajectory geometries considered: a straight line tangential to the atmosphere (red) and an upward-bent trajectory with a radius of curvature r (blue).

Figure 4 shows the integrated fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), expected hard-scattering count $ \langle N\rangle $ (middle), and their sum (right) along an infinite straight-line trajectory tangential to the atmosphere for protons of five rigidities. Although $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ describe distinct physical mechanisms, both are dimensionless measures of the breakdown of energy-conserving backtracing. We heuristically add them to define a unified criterion for trajectory termination. All quantities exhibit the same $ \propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $ dependence on perigee altitude h, and therefore, they can be factorized out. The altitude range of 30-100 km is selected to highlight the variation around the critical perigee altitude h, where interactions become significant. This is considerably lower than the altitude of the observer (e.g., $ \sim400 $ km for AMS-02).

Figure 4.  (color online) Cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), expected hard-scattering count $ \langle N\rangle $ (middle), and their sum (right) versus perigee altitude h for a straight-line tangent trajectory. Proton rigidities: 0.2 GV (red), 0.5 GV (orange), 1 GV (green), 2 GV (blue), and 5 GV (black). The altitude dependence is well described by the common exponential factor $ \exp(-0.14\,{\rm h}/\text{km}) $ (grey dashed) and can therefore be factorized out.

The exponential fall-off of atmospheric density guarantees that both $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ receive their dominant contributions from the segment closest to Earth. The natural control parameter is therefore the perigee altitude h.

As a first baseline, we consider a straight-line trajectory tangential to the atmosphere at altitude h (impact parameter $ b=r_{\oplus}+h $, with $ r_{\oplus}=6371.2 $ km); this is shown by the red curve in Fig. 3. Although this path ignores geomagnetic bending, it captures the correct order of magnitude of the two stopping mechanisms.

Figure 3.  (color online) Two trajectory geometries considered: a straight line tangential to the atmosphere (red) and an upward-bent trajectory with a radius of curvature r (blue).

Figure 4 shows the integrated fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), expected hard-scattering count $ \langle N\rangle $ (middle), and their sum (right) along an infinite straight-line trajectory tangential to the atmosphere for protons of five rigidities. Although $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ describe distinct physical mechanisms, both are dimensionless measures of the breakdown of energy-conserving backtracing. We heuristically add them to define a unified criterion for trajectory termination. All quantities exhibit the same $ \propto {\rm e}^{-0.14{\rm h}/ {\rm{km}}} $ dependence on perigee altitude h, and therefore, they can be factorized out. The altitude range of 30-100 km is selected to highlight the variation around the critical perigee altitude h, where interactions become significant. This is considerably lower than the altitude of the observer (e.g., $ \sim400 $ km for AMS-02).

Figure 4.  (color online) Cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (left), expected hard-scattering count $ \langle N\rangle $ (middle), and their sum (right) versus perigee altitude h for a straight-line tangent trajectory. Proton rigidities: 0.2 GV (red), 0.5 GV (orange), 1 GV (green), 2 GV (blue), and 5 GV (black). The altitude dependence is well described by the common exponential factor $ \exp(-0.14\,{\rm h}/\text{km}) $ (grey dashed) and can therefore be factorized out.

The straight-line tangent ignores bending by the Lorentz force. The realistic trajectories must be obtained by backtracing. However, near perigee, any path is characterised by its local radius of curvature $ r=\mathfrak{R}/B $, which is almost independent of the remainder of the orbit. Therefore, we replace the tangent line by an upward-bending circular arc of radius r that reaches the same perigee altitude h (blue curve in Fig. 3). Downward-bending geometries with $ r>r_{\oplus}+h $ can yield perigee h; however, they are not considered here. The limit $ r\to\infty $ recovers the straight-line case. The parameterized integrated trajectory is a circle with radius r, while we checked that the arc around the perigee contributes to the integration. For simplicity, the arc is considered to lie in the plane containing the centre of the Earth and its perigee point.

Figure 5 compares the numerical integration of $ \langle N \rangle $ with its asymptotic and fitted forms. The behaviour of $ \Delta\mathfrak{R}/\mathfrak{R} $ is identical. The global fit is given by

Figure 5.  (color online) Expected hard-scattering count $ \langle N \rangle $ versus the local radius of curvature r for an upward-bent trajectory with perigee altitude 50 km. The black curve indicates the full calculation; grey dashed curves represent the low-r fit and high-r asymptote; and the red curve represents the global fit.

This is normalised to its $ r\to\infty $ limit, reproducing the calculation effectively. Combining this curvature factor with the altitude dependence yields compact factorized forms, which work well if normalized with its $ r\to\infty $ asymptotic value. With this fitting for local curvature radius r around perigee and the above fitting for perigee altitude h, the two dimensionless variables can be factorized as

The straight-line tangent ignores bending by the Lorentz force. The realistic trajectories must be obtained by backtracing. However, near perigee, any path is characterised by its local radius of curvature $ r=\mathfrak{R}/B $, which is almost independent of the remainder of the orbit. Therefore, we replace the tangent line by an upward-bending circular arc of radius r that reaches the same perigee altitude h (blue curve in Fig. 3). Downward-bending geometries with $ r>r_{\oplus}+h $ can yield perigee h; however, they are not considered here. The limit $ r\to\infty $ recovers the straight-line case. The parameterized integrated trajectory is a circle with radius r, while we checked that the arc around the perigee contributes to the integration. For simplicity, the arc is considered to lie in the plane containing the centre of the Earth and its perigee point.

Figure 5 compares the numerical integration of $ \langle N \rangle $ with its asymptotic and fitted forms. The behaviour of $ \Delta\mathfrak{R}/\mathfrak{R} $ is identical. The global fit is given by

Figure 5.  (color online) Expected hard-scattering count $ \langle N \rangle $ versus the local radius of curvature r for an upward-bent trajectory with perigee altitude 50 km. The black curve indicates the full calculation; grey dashed curves represent the low-r fit and high-r asymptote; and the red curve represents the global fit.

This is normalised to its $ r\to\infty $ limit, reproducing the calculation effectively. Combining this curvature factor with the altitude dependence yields compact factorized forms, which work well if normalized with its $ r\to\infty $ asymptotic value. With this fitting for local curvature radius r around perigee and the above fitting for perigee altitude h, the two dimensionless variables can be factorized as

The straight-line tangent ignores bending by the Lorentz force. The realistic trajectories must be obtained by backtracing. However, near perigee, any path is characterised by its local radius of curvature $ r=\mathfrak{R}/B $, which is almost independent of the remainder of the orbit. Therefore, we replace the tangent line by an upward-bending circular arc of radius r that reaches the same perigee altitude h (blue curve in Fig. 3). Downward-bending geometries with $ r>r_{\oplus}+h $ can yield perigee h; however, they are not considered here. The limit $ r\to\infty $ recovers the straight-line case. The parameterized integrated trajectory is a circle with radius r, while we checked that the arc around the perigee contributes to the integration. For simplicity, the arc is considered to lie in the plane containing the centre of the Earth and its perigee point.

Figure 5 compares the numerical integration of $ \langle N \rangle $ with its asymptotic and fitted forms. The behaviour of $ \Delta\mathfrak{R}/\mathfrak{R} $ is identical. The global fit is given by

Figure 5.  (color online) Expected hard-scattering count $ \langle N \rangle $ versus the local radius of curvature r for an upward-bent trajectory with perigee altitude 50 km. The black curve indicates the full calculation; grey dashed curves represent the low-r fit and high-r asymptote; and the red curve represents the global fit.

This is normalised to its $ r\to\infty $ limit, reproducing the calculation effectively. Combining this curvature factor with the altitude dependence yields compact factorized forms, which work well if normalized with its $ r\to\infty $ asymptotic value. With this fitting for local curvature radius r around perigee and the above fitting for perigee altitude h, the two dimensionless variables can be factorized as

Figure 6 shows the sum $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ versus rigidity curves for five cosmic-ray species on a straight-line trajectory tangent to the atmosphere at $ h=50 $ km. Bethe-Bloch energy loss, scaling approximately as $ \beta^{-2} $, dominates at low rigidities, whereas hard scattering, with its nearly energy-independent cross-section, takes over at high rigidities. The crossover occurs at $ \sim0.57 $ GV for protons and rises to $ \sim1.4 $ GV for iron. These reference values are inserted into Eq. 13 to obtain the full altitude- and curvature-corrected estimates. Values $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim1 $ obtained from the perturbative calculation are rescaled by Eq. 13 to the region $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $.

Figure 6.  (color online) Cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (dotted), expected hard-scattering count $ \langle N \rangle $ (dashed), and their sum (solid) versus the incident rigidity for proton (black), helium (red), oxygen (orange), silicon (green), and iron (blue) in a straight-line trajectory with aa perigee altitude h of 50 km.

For protons at 50 km, $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle\gtrsim1 $ at rigidities in the straight-line tangent geometry. Therefore, this altitude may serve as a reasonable lower boundary when a sharp boundary is imposed in backtracing. The selected threshold of $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ should be smaller yet close to unity. However, the values in Fig. 6 neglect the finite-radius-of-curvature correction. For iron, the sum is roughly ten times larger, and therefore, the same criterion requires the altitude-curvature factor to be $ \sim0.1 $, which is equivalent to raising the tangent altitude by $ \sim16 $ km.

Figure 6 shows the sum $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ versus rigidity curves for five cosmic-ray species on a straight-line trajectory tangent to the atmosphere at $ h=50 $ km. Bethe-Bloch energy loss, scaling approximately as $ \beta^{-2} $, dominates at low rigidities, whereas hard scattering, with its nearly energy-independent cross-section, takes over at high rigidities. The crossover occurs at $ \sim0.57 $ GV for protons and rises to $ \sim1.4 $ GV for iron. These reference values are inserted into Eq. 13 to obtain the full altitude- and curvature-corrected estimates. Values $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim1 $ obtained from the perturbative calculation are rescaled by Eq. 13 to the region $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $.

Figure 6.  (color online) Cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (dotted), expected hard-scattering count $ \langle N \rangle $ (dashed), and their sum (solid) versus the incident rigidity for proton (black), helium (red), oxygen (orange), silicon (green), and iron (blue) in a straight-line trajectory with aa perigee altitude h of 50 km.

For protons at 50 km, $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle\gtrsim1 $ at rigidities in the straight-line tangent geometry. Therefore, this altitude may serve as a reasonable lower boundary when a sharp boundary is imposed in backtracing. The selected threshold of $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ should be smaller yet close to unity. However, the values in Fig. 6 neglect the finite-radius-of-curvature correction. For iron, the sum is roughly ten times larger, and therefore, the same criterion requires the altitude-curvature factor to be $ \sim0.1 $, which is equivalent to raising the tangent altitude by $ \sim16 $ km.

Figure 6 shows the sum $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ versus rigidity curves for five cosmic-ray species on a straight-line trajectory tangent to the atmosphere at $ h=50 $ km. Bethe-Bloch energy loss, scaling approximately as $ \beta^{-2} $, dominates at low rigidities, whereas hard scattering, with its nearly energy-independent cross-section, takes over at high rigidities. The crossover occurs at $ \sim0.57 $ GV for protons and rises to $ \sim1.4 $ GV for iron. These reference values are inserted into Eq. 13 to obtain the full altitude- and curvature-corrected estimates. Values $ \Delta\mathfrak{R}/\mathfrak{R}\gtrsim1 $ obtained from the perturbative calculation are rescaled by Eq. 13 to the region $ \Delta\mathfrak{R}/\mathfrak{R}\ll1 $.

Figure 6.  (color online) Cumulative fractional rigidity loss $ \Delta\mathfrak{R}/\mathfrak{R} $ (dotted), expected hard-scattering count $ \langle N \rangle $ (dashed), and their sum (solid) versus the incident rigidity for proton (black), helium (red), oxygen (orange), silicon (green), and iron (blue) in a straight-line trajectory with aa perigee altitude h of 50 km.

For protons at 50 km, $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle\gtrsim1 $ at rigidities in the straight-line tangent geometry. Therefore, this altitude may serve as a reasonable lower boundary when a sharp boundary is imposed in backtracing. The selected threshold of $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ should be smaller yet close to unity. However, the values in Fig. 6 neglect the finite-radius-of-curvature correction. For iron, the sum is roughly ten times larger, and therefore, the same criterion requires the altitude-curvature factor to be $ \sim0.1 $, which is equivalent to raising the tangent altitude by $ \sim16 $ km.

The physics-based termination criterion improves ad hoc backtracing boundaries. For a 1 GV proton in the straight-line tangent geometry, we obtain $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle= 0.93 $ at 50 km; however, this value is 3.4 at 40 km. Thus, primary charged cosmic rays cannot originate from a perigee just above the formerly adopted 40 km cut-off [3].

As a second illustration, we consider a vertical incidence at $ 0^\circ $ longitude, $ 0^\circ $ latitude, and 400 km altitude. Using a conventional algorithm reported previously [2] (IGRF-13 field for epoch 2023 with the traditional 100 km sharp boundary), the geomagnetic cut-off rigidity is 12.04 GV. A dedicated search shows that for $ \mathfrak{R}= $ $ 12.03809986288667 $ GV, the backtraced trajectory reaches a perigee of 98.3 km before escaping to infinity. The old sharp-boundary prescription can discard this event, while the physics-based criterion developed accepts it.

The physics-based criterion $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ is intrinsically probabilistic, and therefore, the very notion of a penumbra becomes questionable. Penumbra represents the narrow rigidity interval, which is usually around the vertical cut-off, wherein the numerical backtracing switches repeatedly between "allowed" and "forbidden" as the rigidity is scanned in fine steps. The Bethe-Bloch energy loss is deterministic; however, hard scattering supplies only an expected number of collisions. Even when the path integral yields $ \langle N\rangle=1 $, there remains a finite survival probability. Thus, the boundary between allowed and forbidden rigidities is unavoidably smeared, irrespective of whether the backtracing boundary is set at the physics-based 50 km or any other fixed altitude.

The allowed cone (i.e., the set of permitted incident directions for cosmic rays in LEOs) can be approximated at high rigidities by the solid angle above the line of sight tangential to the sharp boundary [2]. For an arrival altitude of 400 km, the use of a 20 km boundary yields an allowed cone 1.0% larger for protons (1.5% for iron) compared to a 50 km (66 km for iron) boundary. Conversely, a 100 km boundary yields a cone 1.7% smaller for protons (1.2% for iron). Consequently, space dosimetry [6] is overestimated by ~1.0%−1.5%. The AMS-02 experiment can utilize this additional 1.2%−1.7% of events, which corresponds to billions of events given its 15 years of data collection.

One practical limitation of the present study is that the local radius of curvature around the perigee is known only after a full backtracing calculation. In such a code, the integrals for $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ can be evaluated event-by-event, without recourse to the parametric factor $ g(r) $ in Eq. 11. This yields the exact energy-loss and scattering probability for each trajectory and permit any desired numerical threshold to be imposed a posteriori. This refinement in backtracing removes uncertainties introduced by the sharp-boundary altitude approximation, which enables a single backtracing calculation to suffice for all cases.

The physics-based termination criterion improves ad hoc backtracing boundaries. For a 1 GV proton in the straight-line tangent geometry, we obtain $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle= 0.93 $ at 50 km; however, this value is 3.4 at 40 km. Thus, primary charged cosmic rays cannot originate from a perigee just above the formerly adopted 40 km cut-off [3].

As a second illustration, we consider a vertical incidence at $ 0^\circ $ longitude, $ 0^\circ $ latitude, and 400 km altitude. Using a conventional algorithm reported previously [2] (IGRF-13 field for epoch 2023 with the traditional 100 km sharp boundary), the geomagnetic cut-off rigidity is 12.04 GV. A dedicated search shows that for $ \mathfrak{R}= $ $ 12.03809986288667 $ GV, the backtraced trajectory reaches a perigee of 98.3 km before escaping to infinity. The old sharp-boundary prescription can discard this event, while the physics-based criterion developed accepts it.

The physics-based criterion $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ is intrinsically probabilistic, and therefore, the very notion of a penumbra becomes questionable. Penumbra represents the narrow rigidity interval, which is usually around the vertical cut-off, wherein the numerical backtracing switches repeatedly between "allowed" and "forbidden" as the rigidity is scanned in fine steps. The Bethe-Bloch energy loss is deterministic; however, hard scattering supplies only an expected number of collisions. Even when the path integral yields $ \langle N\rangle=1 $, there remains a finite survival probability. Thus, the boundary between allowed and forbidden rigidities is unavoidably smeared, irrespective of whether the backtracing boundary is set at the physics-based 50 km or any other fixed altitude.

The allowed cone (i.e., the set of permitted incident directions for cosmic rays in LEOs) can be approximated at high rigidities by the solid angle above the line of sight tangential to the sharp boundary [2]. For an arrival altitude of 400 km, the use of a 20 km boundary yields an allowed cone 1.0% larger for protons (1.5% for iron) compared to a 50 km (66 km for iron) boundary. Conversely, a 100 km boundary yields a cone 1.7% smaller for protons (1.2% for iron). Consequently, space dosimetry [6] is overestimated by ~1.0%−1.5%. The AMS-02 experiment can utilize this additional 1.2%−1.7% of events, which corresponds to billions of events given its 15 years of data collection.

One practical limitation of the present study is that the local radius of curvature around the perigee is known only after a full backtracing calculation. In such a code, the integrals for $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ can be evaluated event-by-event, without recourse to the parametric factor $ g(r) $ in Eq. 11. This yields the exact energy-loss and scattering probability for each trajectory and permit any desired numerical threshold to be imposed a posteriori. This refinement in backtracing removes uncertainties introduced by the sharp-boundary altitude approximation, which enables a single backtracing calculation to suffice for all cases.

The physics-based termination criterion improves ad hoc backtracing boundaries. For a 1 GV proton in the straight-line tangent geometry, we obtain $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle= 0.93 $ at 50 km; however, this value is 3.4 at 40 km. Thus, primary charged cosmic rays cannot originate from a perigee just above the formerly adopted 40 km cut-off [3].

As a second illustration, we consider a vertical incidence at $ 0^\circ $ longitude, $ 0^\circ $ latitude, and 400 km altitude. Using a conventional algorithm reported previously [2] (IGRF-13 field for epoch 2023 with the traditional 100 km sharp boundary), the geomagnetic cut-off rigidity is 12.04 GV. A dedicated search shows that for $ \mathfrak{R}= $ $ 12.03809986288667 $ GV, the backtraced trajectory reaches a perigee of 98.3 km before escaping to infinity. The old sharp-boundary prescription can discard this event, while the physics-based criterion developed accepts it.

The physics-based criterion $ \Delta\mathfrak{R}/\mathfrak{R}+\langle N\rangle $ is intrinsically probabilistic, and therefore, the very notion of a penumbra becomes questionable. Penumbra represents the narrow rigidity interval, which is usually around the vertical cut-off, wherein the numerical backtracing switches repeatedly between "allowed" and "forbidden" as the rigidity is scanned in fine steps. The Bethe-Bloch energy loss is deterministic; however, hard scattering supplies only an expected number of collisions. Even when the path integral yields $ \langle N\rangle=1 $, there remains a finite survival probability. Thus, the boundary between allowed and forbidden rigidities is unavoidably smeared, irrespective of whether the backtracing boundary is set at the physics-based 50 km or any other fixed altitude.

The allowed cone (i.e., the set of permitted incident directions for cosmic rays in LEOs) can be approximated at high rigidities by the solid angle above the line of sight tangential to the sharp boundary [2]. For an arrival altitude of 400 km, the use of a 20 km boundary yields an allowed cone 1.0% larger for protons (1.5% for iron) compared to a 50 km (66 km for iron) boundary. Conversely, a 100 km boundary yields a cone 1.7% smaller for protons (1.2% for iron). Consequently, space dosimetry [6] is overestimated by ~1.0%−1.5%. The AMS-02 experiment can utilize this additional 1.2%−1.7% of events, which corresponds to billions of events given its 15 years of data collection.

One practical limitation of the present study is that the local radius of curvature around the perigee is known only after a full backtracing calculation. In such a code, the integrals for $ \Delta\mathfrak{R}/\mathfrak{R} $ and $ \langle N\rangle $ can be evaluated event-by-event, without recourse to the parametric factor $ g(r) $ in Eq. 11. This yields the exact energy-loss and scattering probability for each trajectory and permit any desired numerical threshold to be imposed a posteriori. This refinement in backtracing removes uncertainties introduced by the sharp-boundary altitude approximation, which enables a single backtracing calculation to suffice for all cases.