Updated model of cosmic-ray-induced ionization in the atmosphere (CRAC:CRII_v3): Improved yield function and lookup tables

– Cosmic rays, including galactic cosmic rays and solar energetic particles, form the main source of ionization of the low and middle atmosphere, which is important for various chemical and physical effects in the atmosphere. Realistic models able to compute the cosmic-ray-induced ionization (CRII) are used as inputs for chemistry-climate models. One of the most commonly used atmospheric ionization models is CRAC:CRII (Cosmic-Ray Atmospheric Cascade: application to CRII) initially developed in 2004 – 2006 (version 1) and signi ﬁ cantly improved in 2010 – 2011 (version 2). Here, a new updated version 3 of the CRAC:CRII model is presented which offers a higher accuracy for the middle-upper atmosphere and lower-energy cosmic rays. This is particularly important for studies of the atmospheric effects of solar particle storms. Detailed lookup tables of the ionization yield function are provided for the primary cosmic ray protons and a -particles (the latter representing also heavier cosmic-ray species) along with a practical recipe for their numerical use.


Introduction
The Earth's atmosphere is constantly bombarded by highenergy particles originating from outer space that are called cosmic rays (CRs).The most important are Galactic CRs (GCRs), which originate from the Galaxythey have omnipresent and nearly isotropic flux near Earth with a slight variability at different timescales.The energy of GCRs is typically in the range of several GeV/nucleon and can reach up to 10 20 eV.GCRs near Earth mostly consist of protons (about 90%), helium (about 9%), and 1% of heavier element nuclei, produced mostly by supernovae in the Galaxy (see, e.g., Gaisser et al., 2016, and references therein).High-energy particles can be also accelerated near the Sun in solar eruptive events such as flares or coronal mass ejectionsthey are called solar energetic particles (SEPs) and have sporadic occurrence and very strong variability (e.g., Desai and Giacalone, 2016).The energy of SEPs is typically below 100 MeV but in rare cases of very energetic events, called ground-level enhancements (GLEs), it can reach several GeV.The third type of CRs, the anomalous cosmic rays originating from interstellar neutral atoms, is not considered here as having too low energy and fluxes to affect the atmosphere.
Because of the geomagnetic field and its shielding effect, the flux of GCRs impinging on Earth is not uniform being higher in polar regions and lowest in the equatorial region.SEPs can enter the Earth's atmosphere in (sub)polar regions with only very rare events being detectable at mid and low latitudes.
When entering the Earth's atmosphere, energetic particles interact with matter losing energy for ionization of the ambient air (e.g., Bazilevskaya et al., 2008).If the particle's energy is sufficiently high, inelastic nuclear collisions become also important leading to the formation of a complicated hadron-electromagnetic-muon cascade of secondary particles, known also as extensive air shower (for details see, e.g., Dorman, 2004;Grieder, 2011;Gaisser et al., 2016).CRs form the main source of atmospheric ionization, called cosmic-ray-induced ionization (CRII) below %100 km height (e.g., Vainio et al., 2009) with the maximum of ion production lying at the altitude of about 12-15 km in the atmosphere (Regener and Pfotzer, 1935).Other effects, such as solar UV and X-rays as well as magnetospheric particles, can drive ionization at higher altitudes (e.g., Mironova et al., 2015).Besides, natural radionuclides can contribute to atmospheric ionization in the near-ground layer (e.g., Eisenbud and Gesell, 1997).CRII plays a crucial role in the physical and chemical atmospheric processes (e.g.Crutzen et al., 1975;Randall et al., 2007;Calisto et al., 2011;Semeniuk et al., 2011;Rozanov et al., 2012;Sinnhuber et al., 2012;Jackman et al., 2016;Golubenko et al., 2020).
While the physics of atmospheric ionization is well understood (e.g., Mironova et al., 2015), the exact modelling of CRII remains slightly uncertain, particularly in the upper-middle atmosphere.The ion production in the upper atmosphere can be calculated using analytical (parameterization) and/or semiempirical models (e.g., Vitt and Jackman, 1996;Turunen et al., 2009;McGranaghan et al., 2015).Such analytical models have limited applicability to only the upper polar atmospheric region and lower-energy particles (below 100 MeV) because they neglect inelastic scattering and the development of the atmospheric cascade.Their applicability at heights below 40 km and outside the polar cusp region is not validated (Desorgher et al., 2009;Velinov et al., 2013;Mironova et al., 2015).CRII in the low and middle atmosphere is well-modelled by full Monte Carlo models which simulate the atmospheric cascade using the modern knowledge of nuclear reactions.The most widely used full Monte-Carlo models are PLANETO-COSMIC (Desorgher et al., 2005(Desorgher et al., , 2009) ) and CRAC:CRII (Cosmic-Ray Atmospheric Casacde: application to CRII - Usoskin et al., 2004;Usoskin and Kovaltsov, 2006;Usoskin et al., 2010).The latter is a reference model, e.g., for the CMIP (Climate Model Intercomparison Projectsee Matthes et al., 2017;Jungclaus et al., 2017;Funke et al., 2024) phases 6 and 7. Full Monte-Carlo models precisely, within the 10% accuracy, compute CRII for the low and middle atmosphere with the residual atmospheric depth >10 g/cm 2 (e.g., Usoskin et al., 2009) but may be less certain in the upper atmosphere, where the cascade is not yet induced, and the Monte Carlo approach is not optimal.Thus, a combination of the Monte Carlo models for the lower and middle atmosphere and the analytical approach for the upper polar atmosphere is required for consistent modelling of ionization of the whole atmosphere.Such work was done firstly by Usoskin et al. (2010), and here we extend and upgrade that work by updating the computations and providing detailed lookup tables so that a user can compute cosmic-ray-induced ionization at any location and initial conditions.

Monte Carlo model
Since energetic particles initiate a complicated cascade in the atmosphere, it is important to model it in a realistic way.At each hadronic interaction, about one-third of primary CR particle energy is transferred to the electromagnetic component of the developing shower, that is electrons, positrons and c-quanta.Secondary hadrons experience subsequent interactions, leading to most of the primary particle's energy being dissipated by ionization energy losses of the electromagnetic component (see, e.g., chapter 16 in Gaisser et al., 2016).The most suitable way to model this process is a Monte Carlo approach which takes into account all known processes (Desorgher et al., 2005;Usoskin and Kovaltsov, 2006).Such models are based on a full Monte Carlo simulation of CR particle propagation and interaction in the atmosphere, employing the recent achievements of particle physics and the corresponding hadron generation simulators.The main advantage of Monte Carlo transport codes is that they consider realistically all the physical processes, including secondary and subsequent-generation particle production and propagation as well as energy losses, which contribute eventually to the ionization of air, specifically at lower and middle altitudes (e.g., Bazilevskaya et al., 2008;Usoskin et al., 2009;Mironova et al., 2015, and references therein).
There are different Monte-Carlo-based models of CRII, but some of them have limitations on the energy range or particle species (e.g.Wissing and Kallenrode, 2009;Mishev and Velinov, 2010).Two most useful models are CRAC:CRII (versions 1 and 2 - Usoskin and Kovaltsov, 2006;Usoskin et al., 2010, respectively) based on COSRIKA simulation tool (COsmic Ray SImulations for KAscade -Heck et al., 1998) and PLANE-TOCOSMICS (Desorgher et al., 2005) based on GEANT4 tool (GEometry ANd Tracking -Agostinelli et al., 2003).It is important that at the present level of knowledge, different Monte Carlo tools would produce similar averaged results (Usoskin et al., 2009) if based on the same atmospheric model and physics list, considering the intrinsic shower fluctuations (see discussion in Engel et al., 2011;Pierog, 2017).Thus, an intermodel comparison is not crucial.
We note that presently only the CRAC:CRII model provides look-up tables in the form of yield functions making fast yet accurate CRII computations possible for any given location and energetic-particle spectrum (see Usoskin and Kovaltsov, 2006, also available at https://cosmicrays.oulu.fi/CRII/CRII.html).
Herewith we update and extend the CRAC:CRII model (Usoskin and Kovaltsov, 2006;Usoskin et al., 2010), upgrading it to version 3. First, we performed the high-statistic Monte Carlo simulations of the atmospheric cascade and the deposited energy for two types of the primary CR particles, protons and aparticles, the latter effectively representing also heavier species (Koldobskiy et al., 2019), assuming their isotropic flux on the top of the atmosphere.The simulations were performed using a Fortran-based CORSIKA code (v.6.617 -Heck et al., 1998) extended with FLUKA (v.2006.3b -Fass`o et al., 2001) in the low-energy range.We applied a realistic spherical geometry of the atmosphere which is important for the isotropic flux, the atmospheric parameters were considered according to the NRLMSISE-00 model (Picone et al., 2002).The statistic of the simulated primary particles was improved to keep the statistical uncertainty generally below 1%the number of simulated showers varied between 10 3 and 10 7 for each energy point depending on the primary particle's energy.
We used an improved grid over energy -20 logarithmically spaced energy values per order of magnitude between 10 MeV and 1000 GeV; and over atmospheric depth -0.01 g/cm 2 steps between 0 and 0.1 g/cm 2 , 0.1 g/cm 2 steps between 0.1 and 1 g/cm 2 , 1 g/cm 2 between 1 and 10 g/cm 2 , and 5 g/cm 2 for the depths greater than 10 g/cm 2 .This grid is much denser than that used in CRAC:CRII versions 1 and 2 (two points over an order of magnitude in energy and 50 g/cm 2 over the atmospheric depth), which allows us to improve the accuracy of the computations.The energy deposition was converted into the ion-pair production assuming 35 eV per ion pair (Porter et al., 1976).The obtained grid values of the CRII were locally smoothed in two dimensions (energy and atmospheric depth).
However, Monte Carlo models can be imprecise in the upper atmosphere, at the atmospheric depth <10 g/cm 2 , where the cascade is not developed and CRII is dominated by the direct ionization by the primary CR particle as illustrated in Fig- ure 1.This process is not always accurately considered in the Monte Carlo simulation tools and may lead to imprecise results in the upper atmosphere (e.g., Mishev andVelinov, 2010, 2014).To overcome this problem we combined the Monte Carlo simulation with the analytical calculation of the direct ionization of the ambient air by the primary particle before its first inelastic nuclear interaction.
3 Analytical solution for the upper atmosphere For the uppermost layer of the atmosphere, with a depth <10 g/cm 2 , ionization can be calculated analytically without employing a full Monte Carlo model.
Let us consider a locally flat atmosphere as illustrated in Figure 2. The vertical dimension is the atmospheric depth (or thickness) h in units of g/cm 2 where h = 0 corresponds to the top of the atmosphere.An energetic particle (proton) with energy E enters the atmosphere at the zenith angle h (henceforth we denote l = cos(h)).If other energy losses but ionization of the ambient air can be neglected as well as the scattering of the particle, the proton's energy e(h) at the depth h can be defined as where X = h/l is the proton's pathlength to the depth h, R(e) is the full pathlength of the proton with energy e before full breaking by ionization loss, defined as where dE dR E 0 ð Þ is the ionizing loss of a particle with energy E 0 in air (Berger et al., 2017).The maximum pathlength of the particle with the initial energy E is denoted as R * = R(E).Let us denote the stopping power (energy loss per unit thickness) of the particle in the air as where energy e is defined from equation (1).Equation ( 1) is obtained for one primary particle impinging on the atmosphere with the zenith angle h.Let us denote the distribution of the particle intensity on the cosine l of the zenith angle on the top of the atmosphere as dI/dl, then one can find the energy deposited (per unit depth) at the atmospheric depth h as where l min = h/R(E) is the minimum cosine of the zenith angle so that the particle with initial energy E can still reach the depth h before losing all its energy.For the isotropic angular distribution (normalized per one particle on the top of the atmosphere) of primary particles, where dI dl ðlÞ ¼ 2l, equation (4) becomes Since for any given R * and h, the cosine l can be expressed via R as and equation (4) becomes The stopping-power range in dry air is known (e.g., Berger et al., 2017) and can be used to compute the deposited energy (per unit atmospheric depth) G for any given E and h.9)) on the atmospheric depth h for primary protons with different kinetic energies, as indicated in the legend.Two components can be observed, viz.direct ionization by the primary particles which dominate the upper atmosphere, and the ionization by secondaries of the atmospheric cascade which dominates in the lower atmosphere.
An extension of the cascade-induced ionization to the upper atmosphere is depicted by thin dashed curves.
In addition to these ionization energy losses, inelasticcollision losses of particles need to be considered (see, e.g., Section 3 in Usoskin et al., 2010).This process is usually neglected in the upper atmosphere but naturally taken care of by the Monte Carlo models in the low and middle atmosphere.
Accordingly, here we performed a smooth junction of the full Monte Carlo modelling performed in most of the atmosphere and direct analytical solutions for the upper atmospheric level.

Yield function and lookup tables
The results of the modelling are presented in the form of tabulated ionization yield functions which allow a user to compute CRII at any given location straightforwardly without the full computationally heavy modelling of the atmospheric cascade.Here, we provide a brief description of the approach, while more details can be found in our earlier works (Usoskin and Kovaltsov, 2006;Usoskin et al., 2010).
From the Monte Carlo simulations, we obtained the mean number of ion pairs per g/cm 2 produced by a single primary particle with energy E in a thin atmospheric layer at the atmospheric depth h, called the ion production function S, as where G(E,h) is the mean energy deposition per unit atmospheric thickness of all the primary and secondary particles at the atmospheric depth h, and E i = 35 eV is the average energy for the production of an ion pair in air (Porter et al., 1976).The production function was computed separately for protons, F p , and a-particles, F a .Usually, the yield-function Y (ion pairs sr cm 2 g À1 ) is considered which quantifies the ion production not for a single primary particle but for a unit intensity of the primary particles with the fixed energy E on the top of the atmosphere.For the isotropic flux of primary particles, the two quantities are related as Y ðh; EÞ ¼ p Á Sðh; EÞ; ð9Þ where p is the geometrical normalization factor.The tabulated values for both protons and a-particles are given in the Supplementary materials to this paper.The values of Y for a-particle are provided per nucleon of the primary a-particle.The ionization rate in the air at the depth h is then defined as an integral of the product of the yield function Y n (h,E) of primary particles of type n with the energy spectrum J n (E) where the summation is over the types n of primary CR particles (protons, a-particles) characterized by the charge and atomic numbers Z n and A n , respectively, and integration is over energy above the energy corresponding to the local geomagnetic rigidity cutoff P c (Cooke et al., 1991) as where E r is the proton's rest mass.
This approach works for all primary particleswhether lower-energy solar energetic particles (SEPs) with a soft spectrum or more energetic galactic cosmic rays (GCRs) with a harder spectrum.Since SEPs consist of mostly protons, the sum in equation ( 10) disappears and only the yield function for protons is considered.GCRs have richer compositions, but heavier (Z > 2) species are similar to a-particles in the sense of the Z/A % 2 ratio implying that they experience similar heliospheric and magnetospheric modulation.Accordingly, it is customary to account for Z > 2 nuclei by scaling the flux a-particle (Koldobskiy et al., 2019).
Here we provide the pre-computed and verified lookup tables of CRII as a function of the atmospheric depth h and energy E. The lookup tables provide the ionization yield functions Y(h, E) with a sufficiently dense grid.This is illustrated in Figures 3 and 4. Figure 3 depicts the tabulated CRII as a function of the primary proton's energy for two layers of h = 1 and 3 g/cm 2 (about 50 and 40 km for the standard atmosphere, respectively), for the results of CRAC:CRII_v3 presented here (red curve) and the previous version 2. As seen, the new denser grid covers the lower-energy range with much higher accuracy, specifically in the energy range of several hundred MeV which is important for SEPs.Accordingly, the new results improve the accuracy of the CRII computation due to SEPs in the upper atmosphere.Figure 4 depicts the tabulated CRII as a function of the atmospheric depth for two values of the primary proton's energy of 0.1 GeV and 0.5 GeV.As seen, the altitude profile of CRII was distorted in model version 2 for low-energy particles in the middle atmosphere h < 25 g/cm 2 (above %25 km), while the results are consistent between the versions for higher energies and deeper atmospheric layers.

Recipe for the use of lookup tables
The provided grid is sufficiently dense so that the values between the grid points can be interpolated using a power-law interpolation (linear in the double logarithmic space) as follows.A simple linear interpolation between the grid points may lead to essential inaccuracy.
Let us estimate the value of a function Y(h,E), where the values of h and E lie between the grid points h i À h i+1 and E j À E j+1 , respectively.Let y, g and e be logS, logh and logE, respectively.Then where The ion production rate Q(h) for any air depth h can be calculated using equation ( 10) where integration is replaced by summation over the grid points in energy.
where F j,n = J(E j )ÁY n (h, E j ) is the value of the integrand at the energy gridpoint j for the particle of type n.Because of the very steep integrand function, the local power-law approximation of the integrand between the grid points j and j+1 is recommended: 6 Example of CRII computations As an example of the use of the lookup tables, we present in Figure 5 the computed atmospheric ionization by GCRs as a function of the geomagnetic latitude and atmospheric depth for the conditions of 2015.The spectrum of GCR was modelled using the force-field parameterization with the value of the modulation potential / = 645 MV corresponding to the year 2015 (Usoskin et al., 2017).The geomagnetic shielding was calculated using the standard vertical geomagnetic rigidity cutoff methodology (see Nevalainen et al., 2013, for details).The geomagnetic field was considered as the International Geomagnetic Reference Field (IGRF -Thébault et al., 2015) for the epoch 2015.As one can see from Figure 5, the maximum of the ion production occurs at the atmospheric level of 50-100 g/cm 2 (10-20 km) in the polar atmosphere.This corresponds to the Regener-Pfotzer maximum where the nucleonic cascade effectively starts.At deeper layers, the intensity of the cascade decreases leading to lower ionization rates.The relatively high ionization rate in the upper layers of the polar atmosphere (h < 50 g/cm 2 ) is caused by the high flux of lower-energy cosmic rays which can directly ionize the air.Such low-energy particles are absent at lower latitudes because of the geomagnetic shielding leading to low CRII in the upper atmosphere.The ionization rate decreases towards the equator due to the geomagnetic shielding but also has a higher location at the Regener-Pfotzer maximum.Thus, the present model reproduces all the main characteristic patterns of the CRII (e.g., Bazilevskaya et al., 2008;Mironova et al., 2015).

Summary
Here we present an updated version (v.3) of the CRAC:CRII model of cosmic-ray-induced ionization of the atmosphere, which is one of the most commonly used ones in related atmospheric research.The updated version includes the core Monte-Carlo-based cascade simulations (similar to that of version 2) for the low and middle atmosphere extended by the results of an analytical integration, using the stopping-power approach, of the direct ionization in the upper atmosphere.The lookup tables of the ionization yield function are provided for protons and a-particles (the latter effectively representing heavier cosmic-ray species) in the energy range from 1 MeV to 1000 GeV and the atmospheric-depth range from the ground (1033 g/cm 2 ) to the top of the atmosphere.The grid at which the lookup tables are presented is much denser than in the previous versions: 20 logarithmically spaced points per decade in energy, 10 equispaced points per decade in the atmospheric depth h 10 g/cm 2 , and 5 g/cm 2 resolution for h > 10 g/cm 2 .A detailed recipe for the numerical use of lookup tables is provided so that it is straightforward to compute the CRII rate at any location and time (providing the cosmic-ray spectrum and geomagnetic shielding are known independently) and thus to implement the results of CRAC:CRII_v.3model to any kind of atmospheric model.The CRAC:CRII_v.3 model is suited to compute the atmospheric ionization induced by galactic cosmic rays, solar energetic particles or any other energetic protons or heavier species with a nearly isotropic spatial distribution.I.G.Usoskin et al.: J. Space Weather Space Clim. 2024, 14, 20

Figure 2 .
Figure2.Scheme of the upper atmosphere used for the analytical solution.The upper horizontal blue line denotes the top of the atmosphere, where a particle with energy E enters at the zenith angle h.Energy deposition at the atmospheric depth h is considered in a thin layer of dh thickness where the particle has energy e.The pathlength of the particle in the air between the top of the atmosphere and the depth h is denoted as X = h/cos(h).

Figure 1 .
Figure 1.Dependence of the ionization yield function Y(h) (see Eq. (9)) on the atmospheric depth h for primary protons with different kinetic energies, as indicated in the legend.Two components can be observed, viz.direct ionization by the primary particles which dominate the upper atmosphere, and the ionization by secondaries of the atmospheric cascade which dominates in the lower atmosphere.An extension of the cascade-induced ionization to the upper atmosphere is depicted by thin dashed curves.

Figure 3 .
Figure3.Comparison between the ionization yield functions Y in the CRAC:CRII versions 2 (blue) and 3 (red) as functions of the primary particles energy for two layers in the upper atmosphere of 1 g/cm 2 (panel a) and 3 g/cm 2 (panel b).

Figure 4 .
Figure 4. Comparison between the ionization yield functions Y in the CRAC:CRII versions 2 (blue) and 3 (red) as functions of the atmospheric depth h for two values of the particle's energy of 0.1 GeV (panel a) and 0.5 GeV (panel b).

Figure 5 .
Figure5.Cosmic-ray-induced ionization (CRII) as a function of the geomagnetic latitude k geom and common logarithm of the atmospheric depth h (in g/cm 2 ) for the geomagnetic field and solar modulation corresponding to the year 2015.