About the Altitude Profile of the Atmospheric Cut-Off of Cosmic Rays: New Revised Assessment

Cosmic rays, high-energy subatomic particles of extraterrestrial origin, are systematically measured by space-borne and ground-based instruments. A specific interest is paid to high-energy ions accelerated during solar eruptions, so-called solar energetic particles. In order to build a comprehensive picture of their nature, it is important to fill the gap and inter-calibrate ground-based and space-borne instruments. Here, we focus on ground-based detectors, specifically neutron monitors, which form a global network and provide continuous recording of cosmic ray intensity and its variability, used also to register relativistic solar energetic particles. The count rate of each neutron monitor is determined by the geomagnetic and atmospheric cut-offs, both being functions of the location. Here, on the basis of Monte Carlo simulations with the PLANETOCOSMICS code and by the employment of a new verified neutron monitor yield function, we assessed the atmospheric cut-off as a function of the altitude, as well as for specific stations located in the polar region. The assessed in this study altitude profile of the atmospheric cut-off for primary cosmic rays builds the basis for the joint analysis of strong solar proton events with different instruments and allows one to clarify recent definitions and related discussions about the new sub-class of events, so-called sub-ground-level enhancements (sub-GLEs).


Introduction
The Earth is permanently bombarded by high-energy subatomic particles of extraterrestrial origin known as cosmic rays (CRs), viz. mostly protons, α-particles, and small amounts of heavier nuclei. Primary CRs are distributed in a wide energy range from about 10 6 eV S. Poluianov stepan.poluianov@oulu.fi A. Mishev alexander.mishev@oulu.fi Table 1 Selected neutron monitors in the polar and sub-polar regions. Columns represent the station name, location, geomagnetic cut-off rigidity, altitude above sea level and type of the monitor. The table includes several high-altitude polar monitors and a new station (SUMT) proposed to extend the network. counting rate for the bulk of the stations, in the polar region, where the geomagnetic shielding is marginal, the atmospheric cut-off plays a major role. This is specifically important for the registration and study of SEPs. High-altitude polar stations (e.g., South Pole, Dome C, Vostok) possess notably lower atmospheric cut-offs than sea-level polar NMs (several such NMs are shown in Table 1 as examples). Therefore, high-altitude polar NMs are considerably more sensitive to SEPs, which allows one to register a recently introduced sub-class of events so-called sub-GLEs (for details see Poluianov et al., 2017). A systematic study of such events can provide the basis to fill the gap between space-borne and ground-based measurements of strong SEPs (e.g., Mishev, Poluianov, and Usoskin, 2017). So far, the atmospheric cut-off was estimated mainly at the sea level by, e.g., latitude surveys, when NMs, usually onboard ships, scan a range of magnetic latitudes, so that the NM counting rate can be measured as a function of the geomagnetic cut-off rigidity (for details see e.g., Dorman et al., 2000;Villoresi et al., 2000;Dorman et al., 2008;Nuntiyakul et al., 2014). One can see that while the NM count rate increases with the diminishing geomagnetic cut-off in low and mid-latitudes, in the polar region, where the geomagnetic cut-off rigidity is about 1 GV or below, the NM count rate reaches a plateau (see Figures 2-4 and 5-7 in Villoresi et al., 2000;Nuntiyakul et al., 2014, respectively). This is due to the fact that in the polar region, the geomagnetic cut-off is lower than the atmospheric one, therefore the atmospheric shielding is more important. Thus, the atmospheric cut-off can be roughly estimated as about 1 GV, which corresponds to 433 MeV nucleon −1 for CR protons (Dorman, 2004;Grieder, 2010). Yet, this experimental assessment is applicable only for sea level detectors, while the atmospheric cut-off as a function of the altitude is poorly studied.
Here, we quantitatively assessed the atmospheric cut-off energy for cosmic-ray protons at different altitudes. We employed Monte Carlo simulations of atmospheric cascades and a recent verified NM yield function. We estimated both: (a) the physical atmospheric cutoff, i.e., the minimum energy of primary CR necessary to produce at least one secondary particle at a given altitude (observation level); (b) the instrumental cut-off, which is related to a particular instrument, namely NM64 at a given altitude, i.e. explicitly considering the registration efficiency of a specific device. We intentionally focused on cosmic-ray protons leaving α-particles and heavier species out of scope because their contribution to the total amount of SEPs is minor (see, e.g., Desai and Giacalone, 2016;Reames, 2019a,b) This work follows a recent publication, where the definitions of GLE and sub-GLEs were revised  and is motivated by an active discussion related to registration and analysis of strong SEPs by ground-based and space-borne instruments (for details see Mishev, Poluianov, and Usoskin, 2017;Raukunen et al., 2018). It answers the question on the minimal energy of SEPs needed to get the event registered by neutron monitors at different altitudes in the polar regions.

Assessment of the Atmospheric Cut-Off Energy from Atmospheric Cascade Simulations
Recently, essential progress of Monte Carlo simulations of CR propagation in the atmosphere was achieved, based mostly on increased computation abilities and new data from hadron accelerators (e.g. Engel, Heck, and Pierog, 2011, and the references therein). A good example is the GEANT4 (Agostinelli et al., 2003) based PLANETOCOSMICS  code for simulations of interactions of cosmic rays with planetary atmospheres. It represents a tool for detailed studies of the evolution of CR-induced cascade, here employed specifically for the Earth atmosphere. The tool simulates interactions and decays of various nuclei, hadrons, mesons, electrons, and photons in the atmosphere up to high and very-high energies. The produced in interactions secondary particles are tracked through the atmosphere until they undergo reactions with air nuclei or decay. The result of the simulations is detailed information about the flux, spectrum, energy deposit, angle distribution of particles at given selected altitude(s) above sea level. The results from the atmospheric cascade simulations induced by primary CRs are in good agreement with experimental data, specifically the secondary neutron flux at different depths (e.g. Goldhagen, Clem, andWilson, 2003, 2004;Pioch et al., 2011;Mishev, 2016;Woolf et al., 2019). This is particularly important for the estimation of the threshold energy of a primary CR particle to produce secondary particles at given altitude (depth) in the atmosphere. In the energy range below one GeV nucleon −1 , individual contributions of secondary components in the CR-induced cascade at the ground level are dominated by neutrons because of the threshold energy required for the production of secondary leptons and mesons, and specifics of shower development (see Clem and Dorman, 2000;Engel, Heck, and Pierog, 2011, and the references therein for details).
Hence, we simulated CR-induced cascades due to primary protons around the energy range where previous studies reported an atmospheric cut-off of about 1 GV in rigidity or 433 MeV nucleon −1 in energy (Dorman, 2004;Grieder, 2010). We simulated cascades with isotropic and vertical incidence in the energy range 200 -500 MeV nucleon −1 employing an updated version of the PLANETOCOSMICS code with the GEANT4 physics list QGSP_BIC_HP (the Quark-Gluon String model for high-energy interactions; Binary Cascade model; High-Precision neutron libraries) and a spherical atmospheric model NRLMSISE-00 (Picone et al., 2002). We considered isotropic and vertical incidence of particles because they determine the two extreme cases of the possible anisotropy of SEPs entering the atmosphere. The total amount of secondary particles, specifically neutrons, at several selected depths ranging from 600 g cm −2 (4.5 km of altitude) to 1033 g cm −2 (sea level) was derived. We built a distribution of the total amount of secondary particles vs. the energy of the primary, for several depths in the atmosphere. For example, the distribution of secondary CRs at sea level is presented in Figure 1A, and at 650 g cm −2 in Figure 1B, respectively. Note that the depth of 650 g cm −2 approximately corresponds to the highest Figure 1 Secondary particles as a function of the energy for isotropic CRs. Panel A corresponds to sea level (1033 g cm −2 of the atmospheric depth), while panel B corresponds to 650 g cm −2 .

Figure 2
Atmospheric cut-off energies for CRs with isotropic incidence at selected depths in the atmosphere. altitude with nearly null rigidity cut-off CR stations, namely DOMC/DOMB and VSTK; see Poluianov et al. (2015) and Table 1 for details.
The derived distribution represents the total amount of secondaries (normalized per unit primary particle flux) impinging the unit area on top of the atmosphere, at a selected depth, as a function of the primary particle energy. We assume one secondary particle reaching the selected depth as the threshold value. Therefore, the atmospheric cut-off is determined as the minimum energy of a primary CR proton necessary to induce a such cascade, that at least one secondary particle (in average) reaches the given atmospheric depth. Hence, we estimated the atmospheric cut-off at sea level to about 428 ± 9 MeV nucleon −1 , and 282 ± 5 MeV nucleon −1 at 650 g cm −2 , respectively. The altitude dependence of the atmospheric cut-off for protons with isotropic incidence is presented in Figure 2, the details are given in Table 2.
Accordingly, similar computations were performed for primary CR protons with vertical incidence, i.e., particles traversing minimal amount of air mass, which have a naturally lower cut-off. Note, in most cases of GLE analysis, SEPs are usually assumed with only vertical incidence (Bütikofer et al., 2009;Plainaki et al., 2014;Mishev et al., 2021). It is actually the case also for sub-GLEs. The distribution of secondary particles at sea level as a function of the energy of the primary CR proton with vertical incidence is given in Figure  3. The altitude dependence of the atmospheric cut-off for protons with vertical incidence is presented in Figure 4, the details are given in Table 2. Detailed information about the estimated atmospheric cut-off for primary CR protons with isotropic and vertical incidence at selected depths is summarized in Table 2. The derived atmospheric cut-off for particles with an isotropic incidence at sea level is comparable with previous reports as well as with some recent estimations (e.g Dorman, 2004;Grieder, 2010;Raukunen et al., 2018). However, naturally the Monte Carlo estimation of the atmospheric cut-off results in slightly lower values, because the detector efficiency is not considered here. Therefore in addition to the physical cut-off, it is necessary to perform similar estimations for a given type of detector, namely neutron monitor. This is particularly important for justification and clarification of the current definition of sub-GLEs, and GLEs concerning the energy of registered SEPs (for details see Poluianov et al., 2017;Raukunen et al., 2018).

Assessment of the Atmospheric Cut-Off Energy with the NM Yield Function
In the previous section, the atmospheric cut-off was estimated from the point of view of an ideal detector, i.e., a detector, which registers secondary nucleonic particles with 100% efficiency. However, real instruments possess limited registration ability.
Here, we considered a standard instrument for monitoring cosmic rays and their variations, namely a neutron monitor. NM count rate records allow one to study CR flux variations at the top of the Earth's atmosphere at different time scales (e.g. diurnal, 11-year sunspot cycle, 22-year solar magnetic cycle), as well as short-term transients such as Forbush decreases and anisotropic cosmic ray enhancements (for details see Moraal, 1976;Debrunner et al., 1988;Lockwood, Debrunner, and Flueckiger, 1990;Gil et al., 2018). Neutron monitors are in practice the main detectors for continuous recording of CR intensity variations (e.g. Moraal, 1976;Debrunner et al., 1988;Gil et al., 2015;Kudela, 2016). In addition, registration of GLEs and sub-GLEs with NMs provides key information for analysis of the spectral and angular characteristics over the whole time span of those events (e.g. Shea and Smart, 1982;Cramp et al., 1997;Bombardieri et al., 2006;Vashenyuk et al., 2006;Mishev, Kocharov, and Usoskin, 2014;Mishev et al., 2018).
NM is a ground-based detector aiming the registration of secondary particles, mostly neutrons, but also protons and a small amount of muons, produced in a cosmic-ray cascade (Simpson, 1957;Clem and Dorman, 2000). The instrument consists of one or several proportional counters sensitive to thermal neutrons and filled with enriched to 10 B boron-trifluoride or 3 He gas, surrounded by a moderator, usually paraffin wax or polyethylene, a lead "producer" and an outer reflector made of the same material as the moderator (for details see Clem and Dorman, 2000;Simpson, 2000;Bütikofer, 2018, and the references therein). It was introduced as a continuous recorder of the CR intensity during the International Geophysical Year (IGY) 1957-1958 (Simpson, 1957). In 1964, the design was optimized. The improved instrument is known as the supermonitor or NM64 (see Hatton and Carmichael, 1964;Carmichael, 1968;Simpson, 2000;Stoker, Dorman, and Clem, 2000, and the references therein for details). Nowadays, the NM64 supermonitors are the main detectors in the global NM network, though other designs are also in use (e.g. . For the computation of the instrumental atmospheric cut-off, we employed the NM yield function, which represents the response of the NM to the unit flux of CRs with given energy impinging the top of the atmosphere (e.g. Clem and Dorman, 2000;Flückiger et al., 2008;Mishev, Usoskin, and Kovaltsov, 2013). The yield function incorporates the full complexity of the atmospheric cascade development, interaction(s) of the secondary particles with the detector, taking into account the registration efficiency itself. Here, we considered the standard 6NM64 detector for the computations.
Using the yield function, the NM count rate N(h) [counts s −1 ] can be computed as where Y i (P , h) is the yield function for incident CRs of ith type (protons, α-particles and heavier ions, which are scaled to and effectively accounted as α-particles (e.g. Mishev, Usoskin, and Kovaltsov, 2013)) in units of [counts m 2 sr], J i (P ) is the differential flux of GCRs of ith type in [s m 2 sr GV/nucleon] −1 , h is the atmospheric depth of the instrument in [g cm −2 ], P is the CR particle's rigidity [GV], P gm.cut is the local geomagnetic cut-off rigidity of the instrument in [GV] (Cooke et al., 1991;Smart et al., 2006). The NM yield function is continuous and only asymptotically approaches zero, i.e., there is not apparent cut-off (see Figure 1 in . Therefore, it is necessary to exclude the geomagnetic effects and elaborate a realistic signal-to-noise ratio in an attempt to assess the atmospheric cut-off.
Here, we considered a NM located in the polar region with the geomagnetic cut-off rigidity P gm.cut = 0 GV and effectively excluded the influence of the geomagnetic field. Equation 1 for ith type of CR particles can be rewritten as an integral over energy E where the differential flux J i (E) is in [s m 2 sr GeV/nucleon] −1 . Note that the energy-rigidity relation for a particle of type i is given by where E 0 = 0.938 GeV is the proton's rest energy, A and Z are the mass and charge numbers, respectively. In turn, equation (2) can be split to a sum of two integrals: where E c is the integration parameter, the details are given below.
Here, we assumed that the first term I 1,i (h, E c ) accounts the major fraction of counts, while the second term I 2,i (h, E c ) accounts the remainder due to low-energy particles, that is, the contribution from zero to E c .
Henceforth, we focus on primary protons (i = p), whilst the contribution of α-particles and heavier ions is insignificant in the particular application to the solar particle events, as discussed above. Naturally, the atmospheric cut-off energy E atm.cut is defined as the maximum energy E c , where the contribution of I 2,p (h, E c ) to the total count rate N(h) is indistinguishable.
The contribution of the term I 2,p (h, E c ) can be indistinguishable to the total count rate because of the random variability of the count rate N(h) around certain mean value N(h) with the standard deviation σ (h). The NM count rate reveals the Poisson distribution. Accordingly, its standard deviation is σ (h) = √ N(h) . The mean count rate N(h) can be computed using Equation 2. For these computations of the NM count rate at several depths in the atmosphere, we employed the newly computed NM yield function by . The new yield function is in very good agreement with the experimental latitude surveys and was recently validated by measurements (Gil et al., 2015;Nuntiyakul et al., 2018). Moreover, it was additionally verified using space-borne data (for details see Koldobskiy et al., 2019). The function is presented in tabular and parametrized forms for convenience.
Here, we considered galactic cosmic rays employing the force-field model (Gleeson and Axford, 1968), which is sufficiently accurate for the aims of this work . The differential energy CR spectrum J i (E) is expressed as where J LIS,i is the differential flux of GCRs of ith type in the local interstellar medium [s m 2 sr GeV/nuc] −1 , E is the GCR kinetic energy [GeV nucleon −1 ], φ is the potential [GV] representing the heliospheric modulation of GCRs. The local interstellar spectrum J LIS,i was taken according to Vos and Potgieter (2015) with the ratio of α-particle and heavier ion nuclei to protons as 0.353 (Koldobskiy et al., 2019).
We computed the integral normalized to the standard deviation I 2,p (h, E c )/σ (h) over a range of E c (Equation 4, Figure 5). There is a threshold T = I th 2,p (h)/σ (h), which is invariant to the total count rate N(h), where the contribution of I 2,p (h, E c ) becomes indistinguishable to the total count: Therefore, the atmospheric cut-off energy E atm.cut is a point, where I 2,p (h, E c )/σ (h) intersects T , as depicted in Figure 5. Note that the experimental latitude surveys, as discussed above, imply the atmospheric cut-off energy E atm.cut of about 433 MeV nucleon −1 at sea level (h = 1033 g/cm 2 ) and T is constant, as depicted with horizontal dashed line in Figure  5.
The computations of I 2,p (h, E c )/σ (h) were performed for the mean modulation potential φ = 652 MV (1951( , http://cosmicrays.oulu.fi/phi/phi.html, Usoskin et al. (2017), as well as for φ = 300 MV and φ = 1200 MV, which correspond to low and high solar activity conditions, respectively. Thereby, the influence of cosmic ray modulation on atmospheric cut-off was estimated. Details for those computations are given in Table 3. One can see that the solar variability causes only marginal effect of about ±1% on the derived atmospheric cut-off.
The altitude profiles of the atmospheric cut-off energies obtained with the cascade simulations and with the verified NM yield function are presented in Figure 6. Naturally, a smooth, monotonic increase of the atmospheric cut-off from low depths towards sea level, is observed. As expected, the atmospheric cut-off computed specifically for NM is greater compared to that obtained from the atmospheric cascade simulations. As an example, it is  greater by 26 MeV (7%) at 800 g/cm 2 and by 36 MeV (12%) at 600 g/cm 2 , respectively. This is due to specifics of the altitudinal profile of NM yield function leading to increase of the absolute standard deviation σ (h) of the count rate as a function of the altitude, in spite of its reduction in relative units σ (h)/ N(h) , as well as intrinsic cascade fluctuations considered in shower simulations (for details see Engel, Heck, and Pierog, 2011;.

Discussion and Conclusions
During the last decades, SEPs have been extensively studied using space probes, such as GOES (Geo-stationary Operational Environmental Satellites). However, most of the spaceborne experiments were/are focused on measurements in the low-energy range of solar ions, therefore possess apparent constraints for the precise study of the most energetic ones, e.g. leading to GLEs. Yet, SoHO/EPHIN (Solar-Heliospheric Observatory/Electron Proton and Helium Instrument) and GOES/HEPAD (High Energy Proton and Alpha Detector) have enhanced energy channels up to 500 MeV nucleon −1 and up to about 700 MeV nucleon −1 , respectively (for details see Kühl et al., 2017), as well as PAMELA (Payload for Antimatter Figure 6 Comparison of the atmospheric cut-off energy values computed with the cascade simulations and NM yield function (the latter with J LIS,i according to Vos and Potgieter (2015), φ = 652 MV, which is the mean for 1951-2019, Usoskin et al. (2017)). Blue horizontal dash lines correspond to the atmospheric depths of DOMC/DOMB and SOPO/SOPB NMs, respectively.
Matter Exploration and Light-nuclei Astrophysics) and AMS-02 (Alpha Magnetic Spectrometer) designed to measure directly more energetic particles (e.g., Adriani et al., 2016). Therefore, assessment of the atmospheric cut-off for a specific device such as NM, which is widely used for GLE studies, gives the basis to fill the gap and inter-calibrate the groundbased and space-borne instruments for SEP measurements.
Here, the most sensitive high-altitude polar NMs revealed the following atmospheric cut-off energies: 283 MeV (atmospheric cascade simulation) and 314 MeV (NM yield function employment) for DOMC/DOMB and 299 MeV (atmospheric cascade simulation) and 328 MeV (NM yield function employment) for SOPO/SOPB ( Figure 6). One can see that the atmospheric cut-off assessed with cascade simulations is slightly lower than one estimated with the NM yield function. While the atmospheric cascade simulations are related to an ideal detector with 100% registration efficiency and reveal the physical signal, the more conservative approach for a specific detector, namely NM, is more realistic. This is particularly important for clarification of recent definitions and related discussions about sub-GLEs, a specific sub-class of SEPs, observed by SOPO/SOPB and DOMC/DOMB NM and at space, but not at sea-level instruments Raukunen et al., 2018). We can conclude that the definitions of sub-GLE given by Raukunen et al. (2018) (based on the lower energy limit 300 MeV) and by Poluianov et al. (2017) (based on the simultaneous observation by SOPO and DOMC), are in very good agreement.
The assessed in this study atmospheric cut-off for primary CRs, specifically employed for SEPs, gives the basis for the joint analysis of GLEs, filling the gap between ground-based and space-borne instruments and the development of forefront methods for their studies.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.