Analytical Model for Estimating the Zenith Angle Dependence of Terrestrial Cosmic Ray Fluxes

A new model called “PHITS-based Analytical Radiation Model in the Atmosphere (PARMA) version 4.0” was developed to facilitate instantaneous estimation of not only omnidirectional but also angular differential energy spectra of cosmic ray fluxes anywhere in Earth’s atmosphere at nearly any given time. It consists of its previous version, PARMA3.0, for calculating the omnidirectional fluxes and several mathematical functions proposed in this study for expressing their zenith-angle dependences. The numerical values of the parameters used in these functions were fitted to reproduce the results of the extensive air shower simulation performed by Particle and Heavy Ion Transport code System (PHITS). The angular distributions of ground-level muons at large zenith angles were specially determined by introducing an optional function developed on the basis of experimental data. The accuracy of PARMA4.0 was closely verified using multiple sets of experimental data obtained under various global conditions. This extension enlarges the model’s applicability to more areas of research, including design of cosmic-ray detectors, muon radiography, soil moisture monitoring, and cosmic-ray shielding calculation. PARMA4.0 is available freely and is easy to use, as implemented in the open-access EXcel-based Program for Calculating Atmospheric Cosmic-ray Spectrum (EXPACS).


Introduction
High-energy galactic cosmic rays (GCR) can penetrate the Earth's magnetosphere and produce a variety of secondary particles in the atmosphere by inducing extensive air shower (EAS). Estimation of terrestrial cosmic-ray fluxes is of great importance not only for particle physics and astrophysics but also for the geosciences and radiation research. Thus, a number of studies were devoted for their evaluation on the basis of analytical approaches and Monte Carlo methods [1][2][3][4][5][6][7][8][9][10][11][12][13][14].
Previously, we developed an analytical model for estimating terrestrial cosmic-ray fluxes anywhere in Earth's atmosphere at nearly any given time [15][16][17] by modeling the results of an EAS simulation performed using the Particle and Heavy Ion Transport Code System (PHITS) [18]. The model comprised several theoretical or empirical equations with free parameters, whose numerical values were determined from the least square (LSq) fitting of EAS data. Omnidirectional fluxes of neutrons, protons, ions with a charge of up to 28 (Ni), muons, electrons, positrons, and photons could be calculated using the model over an energy range of 10 keV-1 TeV (per nucleon for ions), with the exception of neutrons, for which fluxes could be calculated down to 0.01 eV. The model was named PHITS-based Analytical Radiation Model in the Atmosphere (PARMA) and was implemented through the open-access EXcel-based Program for calculating Atmospheric Cosmic-ray Spectrum (EXPACS) [19]. PARMA and EXPACS have been extensively used in various research fields, including radiation protection [20][21][22], semiconductor design [23,24], and the geosciences [25][26][27][28][29].
In addition to the omnidirectional fluxes, their angular distributions are also requested to be evaluated in some applications. For example, the estimation of muon fluxes, particularly in the horizontal direction, is useful in the planning of muon radiography of large objects such as volcanoes [30,31], nuclear reactors [32,33], and pyramids [34]. Accurate modeling of the generation of terrestrial cosmic rays considering their angular distributions is the key issue in the Monte Carlo simulation used in designing cosmic-ray detectors [35], shielding calculation for electric devices [36], and soil moisture monitoring [37]. However, there is no model that can instantaneously calculate angular differential cosmic-ray fluxes in the atmosphere for all conditions. The cosmic-ray shower library, CRY [38], can generate cosmic-ray particle shower distribution for use as input to such Monte Carlo simulations, but it is applicable only to three discrete elevations (sea level, 2,100 m, and 11,300 m).
With such situations in mind, we extended our previously established model, PARMA3.0, to make it capable of calculating not only omnidirectional but also angular differential terrestrial cosmic-ray fluxes. This extended version was designated PARMA4.0. The results of EAS simulation performed in our previous study [17] were employed to improve the model, and we carried out additional EAS simulations by changing the ground conditions in order to investigate the effect of Earth's albedo on cosmic-ray angular distributions. A brief outline of the procedures of the EAS simulation together with analysis of the obtained angular distributions are given in the section titled "EAS Simulation", details on the extended capability of PARMA are described in the section titled "Development of PARMA4.0", and comparisons of the results obtained from PARMA4.0 with those from the EAS simulation and several experiments are presented in the section titled "Verification of PARMA4.0". The final section presents our concluding remarks.

Simulation Procedure
The procedure for our EAS simulation is described in detail in our previous paper [17], and therefore only a brief description is given in this paper. The atmosphere is divided into 28 concentric spherical shells, and its maximum altitude is assumed to be 86 km. The densities of each shell are determined by referring to US Standard Atmosphere 1976. The Earth is represented as a sphere with a radius of 6,378.14 km, and its composition is presumed to be the same as that of the air at sea level in order to obtain terrestrial cosmic-ray fluxes under the ideal condition, i.e., without disturbance from the ground. This ground condition is referred to "ideal atmosphere" in this paper.
In the EAS simulation, cosmic rays are incident from the top of the atmosphere in the isotropic irradiation geometry. GCR protons and heavy ions with energies and charges of up to 1 TeV/n and 28 (Ni), respectively, were considered as source particles. The energy spectra of GCR were determined by the model proposed by Matthiä et al. [39]. EAS simulations were conducted for near solar minimum and maximum conditions, i.e., W = 0 and 150, respectively, and for 21 geomagnetic field locations with vertical cut-off rigidities r c ranging from 0 to 20 GV. Note that the solar activity index, W, roughly represents the sun spot number, but its numerical values were determined from the count rates of several neutron monitors in our study. Note that the trapped particles as well as reentrant albedo particles were not taken into account as the source particles in the simulation.
The atmospheric propagation of incident cosmic rays and their associated EAS were simulated using PHITS version 2.73. Default nuclear reaction models and data libraries adopted in PHITS2.73, such as INCL4.6 [40], were employed in the simulation, except in the total reaction cross section model, which was particularly adjusted for high-energy particle transport simulations [41]. In the simulation, all particles were traced down to 10 keV, with the exception of neutrons, which were transported down to 0.01 eV. A variance reduction technique was adopted for transporting low-energy electrons, positrons, and photons in order to reduce the computational time. The angular distribution of particles crossing 18 surfaces at altitudes of between 0 and 52 km were scored as a function of cos(θ), where θ is the zenith angle to the downward direction. As shown in our previous paper [17], the accuracy of the EAS simulation was carefully confirmed using experimental omnidirectional flux data for neutrons and vertical downward flux data for the other particles.
In order to investigate the effect of Earth's albedo on cosmic-ray angular distributions, we conducted two additional EAS simulations that involved changing the composition of Earth: one simulation employed a realistic ground surface, while the other featured a virtual black hole that absorbed all radiation, i.e., it had no albedo effect. The composition of the realistic ground surface was assumed to be 60% SiO 2 , 20% Al 2 O 3 , and 20% H 2 O by mass. These simulations were performed for only the polar region at the solar minimum condition, i.e., at W = 0 and r c = 0 GV, and only the angular distributions at the ground surface were scored. The results obtained from the virtual black hole condition can be used in source generation for underground cosmic-ray transport simulation, which is useful in applications such as soil moisture monitoring with neutrons [37].

Results of the EAS Simulation
In this subsection, the angular distributions of terrestrial cosmic-ray fluxes for various conditions are presented in units of sr -1 , which are to be determined by the analytical model proposed in the next section. Note that cos(θ) = 1 indicates the vertical downward direction. The integral of the angular differential fluxes is normalized to 1.0. Fig 1 shows examples of angular distributions for neutron and proton fluxes for the solar minimum and maximum conditions. It is evident from the figure that the angular distributions for the two solar conditions are nearly identical to each other. Similar tendencies are also seen for other particles and condition cases, indicating that the angular distributions of terrestrial cosmic-ray fluxes are almost independent of the solar activity. We therefore adopt the mean angular distributions between the solar minimum and maximum conditions in the further analysis presented in this paper.  show examples of angular distributions for neutrons, protons, He ions, μ ± , e ± , and photons, respectively, for various altitudes at r c = 1 and 15 GV. The results obtained from the analytical model proposed in the next section, PARMA4.0, are also shown in the figures. As the angular distributions of electron and positron fluxes are very similar to each other, they are not distinguished in the analysis presented in this paper. The statistical uncertainties of the EAS data are less than 20% in most cases, but they become larger in some situations, particularly around cos(θ) = 0, because particle fluxes were calculated based on the number of particles passing through a surface divided by cos(θ) in the EAS simulation. The statistical uncertainties of the He ion data are also very large at locations where the primary ions cannot penetrate, e.g., no high-energy He ions are observed at sea level.
It is evident from the figures that the angular distributions vary significantly with altitude, but they are not as sensitive to the vertical cut-off rigidity, except in the cases of high-energy protons and He ions at high altitudes. It is also seen that most EAS data cannot be expressed on the basis of simple equations such as cos n (θ) that are frequently used for representing the angular distributions of high-energy cosmic-rays at ground level. At lower altitudes, the    angular distributions for high-energy particles exhibit strong downward directivity and become closer to isotropic with decreasing energy. With an increase in altitude, the downward directivity becomes less significant, and the angles of peak flux gradually shift toward the horizontal direction. This tendency can be explained in terms of the Pfotzer maximum, which arises because primary cosmic rays must travel a certain distance through the atmosphere in order to fully form a cascade of secondary particles; as the mean distance that cosmic rays travel at a given altitude becomes longer with decreasing cos(θ), the peak flux correspondingly shifts toward the horizontal at high altitudes. Conversely, the angular distributions for particles dominantly composed of primary cosmic rays, such as protons and He ions over 1 GeV/n at 52 km and 1 GV, are nearly constant between 0 < cos(θ) < 1, as shown in panel (D) of Figs 3 and 4. The primary and secondary cosmic rays have a dominant contribution only at lower and higher r c , respectively, and this is the reason why the angular distributions for the proton and He ion data at high altitudes significantly depend on the vertical cut-off rigidity. Fig 8 shows examples of the angular distributions of low-energy neutron fluxes at sea level obtained from EAS simulation employing a realistic ground surface. The peak observed around cos(θ) = 0 is probably caused by large statistical uncertainties of the EAS data in the horizontal direction. The results calculated using the equations proposed in the next section are also shown in the figure. It is seen that the upward fluxes are slightly larger than the downward fluxes at the thermal energy, i.e., 0.025 eV, while the angular distributions are almost isotropic for energies above 1 eV; this is because the neutrons tend to be thermalized in ground as suggested in our previous study [15]. Note that the ground surface has little influence on the angular distributions for higher energy neutrons as well as for other particles. Fig 9 shows the ratios of particle fluxes at sea level between the EAS simulation results obtained employing the virtual black hole and those obtained using the ideal atmosphere as the Earth. The ratio is equal to 0.0 in the upward direction, i.e., cos(θ) < 0, because there are no albedo particles from the virtual black hole. In the downward direction, the ratios become larger (and asymptotic to 1.0) with increasing energy and cos(θ). As indicated in the figure, these data were also reproduced by the equations proposed in the next section.  unit of the output flux is either cm −2 s −1 MeV −1 or cm −2 s −1 MeV −1 sr -1 , with input parameters including atmospheric depth d in g/cm 2 , vertical cut-off rigidity r c in GV, solar modulation index W, kinetic energy E in MeV, and the zenith angle θ, with the last factor required only for calculating angular differential fluxes. Note that the ion fluxes are output in cm −2 s −1 (MeV/ n) −1 sr -1 or cm −2 s −1 (MeV/n) −1 and their kinetic energies should be given in MeV/n.
In PARMA4.0, the angular differential flux of terrestrial cosmic rays, φ, is expressed as the product of their omnidirectional flux in the ideal atmosphere, ϕ omni , the correction factor for local geometry effects, f L , and the normalized angular distribution, F, as given by where x is cos(θ) and g is the local geometry parameter, e.g., water density in the ground or aircraft mass, as described in [15]. The values of ϕ omni and f L were obtained from a slightly revised version of PARMA3.0, while F was calculated using mathematical functions proposed later in this section. We assumed that F was independent of the solar activity, W, based on the discussion given in the previous section. Note that the azimuth angle dependence is not modeled in PARMA4.0 because the east-west effect was not considered in our EAS simulation. The particle energies covered by PARMA4.0 are 10 keV-1 TeV (per nucleon for ions), except for neutrons and muons, whose fluxes can be calculated down to 0.01 eV and up to 100 TeV, respectively. However, the calculated fluxes for neutrons, e ± , and photons with energies above approximately 5 GeV are not precise enough for the proposed method, as will be discussed later. In terms of global conditions, the applicable altitude, geomagnetic, and temporal ranges of PARMA4.0 are from sea level to the top of the atmosphere, from the polar region to the equatorial region, and from minimum to maximum solar activities recorded in the last 400 years, respectively. Note that the normalized angular distributions are assumed to remain constant above 52 km, the altitude up to which angular differential fluxes were scored in our EAS simulation.

Normalized Angular Distributions in the Ideal Atmosphere
As shown in Figs 2-7, the shapes of the normalized angular distributions are so complicated that they cannot be expressed by a simple function. For example, sudden increases or decreases of the fluxes are observed at cos(θ) = 0, particularly in the data for higher altitudes; the highenergy fluxes reach their maximum level at an intermediate angle around cos(θ) = 0.5 at the altitude of 21 km. We therefore divided the angular distributions into four components, namely upward, horizontal, downward, and extra-downward fluxes, and expressed them as follows: Fðx; EÞ ¼ a 4 ðEÞ þ a 5 ðEÞx a 6 ðEÞ for where F(x,E) is the normalized angular distribution in the ideal atmosphere at cos(θ) = x for a particle with energy E, the parameters x UH , x HD , and x DE represent the switching angles from upward to horizontal, horizontal to downward, and downward to extra-downward fluxes, respectively, and a 1 to a 7 are free parameters. Both F and a i depend not only on E, but also on d and r c , although their dependencies are not explicitly denoted in these equations because they are not continuous functions of d and r c . The numerical values of a i were determined from LSq fitting to the EAS data for each energy, altitude, and vertical cut-off rigidity. However, the values for He ions below an altitude of 8 km were assumed to be the same as those at 8 km, owing to significantly high statistical uncertainties of the EAS data, as shown in Fig 4. For lower altitudes at which F smoothly and continuously increases from the upward to the downward direction, Eq (3) and Eq (5) are not necessary, and the numerical values of a 1 and a 4 should be identical. At around cos(θ) = 0, however, F dramatically changes for higher altitudes and thus is simply determined from the linear interpolation of F(x UH ,E) and F(x HD ,E) obtained from Eqs (2) and (4), respectively, in our model. The numerical values of x UH and x HD were fixed at -0.05 and 0.05, respectively, for all conditions. For the case in which a peak is observed at an intermediate angle of approximately cos(θ) = 0.5, F above the peak angle is also determined from the linear interpolation of F(x DE ,E) and a 7 , where the values of x DE and a 7 are fitted to reproduce the peak angle and F(1.0,E), respectively, obtained from the EAS simulation. Thus, x DE is regarded as a free parameter a 8 , unlike the fixed x UH and x HD . Note that the value of a 8 is equal to 1.0 for most cases, and, therefore, Eq (5) is not used for calculating F for those conditions.
As examples of the results from the LSq fitting, Fig 10 shows the evaluated a 1 and a 4 for neutrons at r c = 1 GV for various altitudes as a function of energy. In order to obtain smooth curves, the EAS data were averaged over nearly one decade of energy in the LSq fitting, e.g., the a 1 value for 1 MeV was determined from the mean of the EAS data between approximately 0.3 and 3 MeV. The evaluated a 1 and a 4 are asymptotic to 1/4π with decreasing energy, indicating the isotropic distribution for low-energy neutrons. The values of a 1 and a 4 are close to each other except at higher altitudes, where sudden increases or decreases of the fluxes are observed in the horizontal direction. Fig 11 shows the evaluated a 6 parameters for neutrons, protons, μ ± , e ± , and photons above 10 MeV at sea level and r c = 1 GV. The angular distributions for these conditions can be approximated by a frequently used function of cos n (θ), and thus a 6 is a rough indication of n. It is seen from the graph that the a 6 parameters for neutrons and protons increase with energy, while the relation is reversed for other particles. The reason for decreasing a 6 for high-energy leptons is that they are predominantly generated through extremely high-energy muons, whose angular distributions are nearly isotropic or even increasing in the horizontal direction, as discussed later. Note that the sudden increase in the a 6 parameter for high-energy photons is probably caused by the poor statistics of the EAS data owing to the very small fluxes of such high-energy photons.
One of the fundamental development strategies of PARMA is to reduce the number of free parameters as much as possible. We therefore express the energy dependences of the parameters a 1 to a 8 as the sum of three sigmoid functions: where b i,1 to b i,10 are free parameters determined for each altitude and vertical cut-off rigidity where EAS data are available. This equation is introduced only to reproduce the complicated energy dependences of a i , and its form has little physical meaning. In the cases where a i has a relatively simple energy dependence, the second and/or third sigmoid functions are not necessary, i.e., b i,5 and/or b i,8 = 0. The numerical values of b i,j were determined from LSq fitting to the evaluated a i parameters, with the results of the fitting shown in Figs 10 and 11. It is evident from the figures that Eq (6) can reproduce the data very well. There are, however, two problems in the calculation of a i using Eq (6). When upward fluxes are much smaller than downward fluxes, as is the case for high-energy particles at sea level, the value of a 2 calculated by Eq (6) is occasionally smaller than (-a 1 ), and, consequently, F(-1.0,E) calculated by Eq (2) becomes negative. In such cases, the a 2 value is adjusted to produce a positive F(-1.0,E), though this adjustment sometimes induces a discontinuity in the energy spectrum. Similar adjustments are also necessary for a 5 when the downward fluxes are much smaller than the upward fluxes. The other problem is the necessity of renormalization of F. The integral of F(x,E) with respect to x between -1 and 1 should be equal to 1/2π, but this is not true when the a i calculated by Eq (6) are directly substituted into Eqs (2) to (5). Thus, the angular distributions must be renormalized by multiplying by the ratio of 1/2π and the original integral value of F.
In contrast to energy, the altitude and vertical cut-off rigidity dependences of b i,j are so complex that they cannot be expressed in simple form. Thus, we evaluated b i,j for 18 altitudes and 7 vertical cut-off rigidities and determined F for intermediate conditions by simply interpolating F where the evaluated b i,j were available.
The calculated angular distributions are shown in Figs 2-7. The agreements between the EAS and PARMA4.0 results are generally satisfactory, except for the cases in which the EAS data are significantly fluctuated owing to large statistical uncertainties such as those in the He ion data at sea level. Disagreements are also observed when the EAS data have a sharp peak in the horizontal direction, as is true for the high-energy neutron data at 52 km. This occurs because the angular distributions between -0.05 < cos(θ) < 0.05 are simply interpolated using Eq (3) in PARMA4.0. A more comprehensive comparison between the EAS and PARMA4.0 data will be provided in the next section.

Correction for the Ground and Black Hole Effects
As shown in Fig 8, the presence of ground changes only the angular distributions of neutrons around the thermal energy. In order to express this change, we introduced the corrected a parameter, a G , for a 1 , a 2 , a 4 and a 5 as given by: where c i,1 to c i,3 are free parameters determined from LSq fitting to the EAS data at sea level for the ground condition. The evaluated values of c i,2 and c i,3 are around -7.0 and 0.5, respectively, and thus, the values of a G , i are nearly equal to a i for neutrons with energies above approximately 1 eV because of a large value of exp{[log 10 (E)−c i,2 ]/c i,3 }. This result confirms that it is not necessary to consider the ground effect in the calculation of neutron angular distributions except around the thermal energy. Replacing a i with a G , i in Eqs (2) to (5), the normalized angular distributions for neutron fluxes at the ground level can be calculated. The results of these calculations are shown in Fig 8, from which it is evident that the normalized angular distributions obtained from the EAS simulation and Eqs (2) to (5) with a G , i agree very closely. For calculating the angular differential fluxes incident to the virtual black hole, the correction factor for the local geometry effects, f L (g,E), in Eq (1) should be replaced by that for the virtual black hole effect, f B (x,E). This correction factor can be determined from the ratio of particle fluxes at sea level obtained from EAS simulations respectively employing the virtual black hole and the ideal atmosphere as the Earth, examples of which are shown in Fig 9. To express f B (x,E), we propose the following equations: f B ðx; EÞ ¼ a 9 ðEÞ þ ½a 10 ðEÞ À a 9 ðEÞx a 11 ðEÞ for 0 x 1; ð9Þ where a 9 to a 11 are free parameters determined from LSq fitting to the EAS data. The a 9 and a 10 parameters represent f B (0,E) and f B (1,E), respectively, and they are asymptotic to 1.0 with increasing E. The energy dependences of a 9 to a 11 are also expressed by Eq (6).
The calculated f B based on the evaluated a i parameters are shown in Fig 9. It is evident that Eqs (8) and (9) can reproduce the EAS data satisfactorily, although some discrepancies are observed owing to large statistical uncertainties in the EAS data, particularly for low-energy electrons and positrons. Note that the a 9 to a 11 parameters for protons, He ions, and muons are assumed to be 1.0 because most albedo-charged particles other than electrons and positrons do not return to the ground owing to their small scattering cross sections at large angles.

Correction for Muon Fluxes
Estimation of the angular differential energy spectra of high-energy muon fluxes at ground level, particularly in the horizontal direction, is of great importance in the planning of muon radiography. However, as shown in Fig 12, the model proposed above cannot reproduce the measured muon fluxes at ground level for large zenith angles [42,43] or even in the vertical direction [44,45] when the muon energy is over 100 GeV. This is because such muons are generally produced by primary cosmic rays with energies over 1 TeV/n, which are not considered in our EAS simulation. The predictability of such high-energy muon fluxes is crucial in the model for applying to the design of muon radiography, and thus, we introduce the correction factors to fix the discrepancy.
The discrepancy observed in high-energy muon fluxes for the vertical direction is attributable to the inaccuracy of PARMA3.0 calculating their omnidirectional fluxes instead of their angular distributions, as the model was designed to reproduce the EAS data up to around 100 GeV. The power index for expressing high-energy muon spectra at sea level used in PARMA3.0 is -3.23, which is much larger than the values obtained from experimental data [44,45]. Thus, we revised PARMA3.0 to reduce the power index of muons with energy above a certain threshold, E t , and normalized the higher energy muon fluxes by E g d t , where γ d is the difference between the power indices below and above E t . To reproduce the experimental data, the numerical values of E t and γ d were determined to be approximately 30 GeV and 0.4, respectively. The vertical muon fluxes calculated using high-energy correction are shown in panel (B) of Fig 12. Excellent agreement is seen between the measurements and the revised calculations, even for energies above 100 GeV. Owing to this correction, the applicable energy of muons covered by PARMA4.0 is extended up to 100 TeV. Note that this correction is considered in calculating the omnidirectional muon fluxes not only at ground level but also at higher altitudes.
On the other hand, the discrepancies observed in muon fluxes for large zenith angles are attributed to the inaccuracy of calculating their angular distributions, particularly at ground where φ G represents the angular differential fluxes of ground-level muons and f G is the ratio between muon fluxes at cos(θ) = x and 1.0, which is assumed to be independent of global conditions. According to [46], for extremely high-energy muons, f G should be equal to 1/cos(θ), but for lower energies, it should decrease as 1/cos(θ) increases. Fig 13 shows values of f G obtained from the measured muon fluxes at θ above 78° [43], i.e. 1/cos(θ) > 5, divided by the corresponding vertical fluxes calculated by our model. To express this dependence, we propose the following function: where a 12 to a 14 are free parameters depending on energy. The first term of the numerator can be approximated by 1/x and a 12 for small and large 1/x, respectively, while the second term expresses the decrease of f G with increasing 1/x. The denominator is introduced to normalize f G (1,E) to 1.0, and f G is asymptotic to a 12 −1 for large 1/x. The numerical values of a 12 to a 14 were determined from LSq fitting to the data shown in Fig 13 for each energy, and their energy dependences are expressed by Eq (6) for a 12 and the third-order polynomial functions of log 10 (E) for a 13 and a 14 . The values of f G obtained from Eq (11) substituting the calculated a i are shown in Fig 13. It is evident that Eq (11) can reproduce the experimentally determined values very well.
Substituting the value of f G obtained from Eq (11) into Eq (10), the ground-level muon fluxes can be precisely calculated, even for larger zenith angles. However, f G can be evaluated Fig 13. Ratio between muon fluxes at a specific zenith angle θ to those in the vertical direction, f G , as a function of 1/cos(θ). Dots are determined from the ratio between measured fluxes [43] and PARMA4.0-calculated vertical fluxes, while lines are calculated by Eq (11).
doi:10.1371/journal.pone.0160390.g013 only above 11.5 GeV, which is the minimum energy of the muons that were measured in [43]. Below this energy, the ground-level muon spectra are calculated by Eq (1) in the same manner as those in the ideal atmosphere, but their absolute values are normalized to smoothly connect with φ G at E = 11.5 GeV.
The calculated muon fluxes after considering both high-energy and ground-level corrections are shown in panel (A) of Fig 12. It is evident from the figure that these corrections significantly improve the capability of reproducing experimental data, including the data measured by Jokisch et al. [42], which were not employed in the LSq fitting for determining the a 12 to a 14 parameters. This result indicates the accuracy of our calculation for various conditions, even when experimental data are not available. Note that the ground-level correction is considered when ground instead of ideal atmosphere is selected in PARMA4.0, and, in that case, the integral of the angular differential fluxes with respect to angle is not equal to the corresponding omnidirectional flux.

Verification of PARMA4.0 Comparison with EAS Simulation
Figs 14-16 show the angular differential energy spectra of cosmic ray fluxes for θ = 0 and 75°at the altitudes of 0, 11 and 52 km, respectively, obtained from the EAS simulation and from PARMA4.0. As expected from Figs 2-7, PARMA4.0 can reproduce most EAS data very well and predict smooth curves, even for conditions in which the EAS data are not robust because of large statistical uncertainties. However, some discrepancies are observed in the neutron, e ± , and photon data over 5 GeV at 0 and 11 km. This discrepancy is attributed to the inaccuracy in PARMA3.0's calculation of the omnidirectional fluxes in a manner similar to its overestimation of high-energy muon fluxes discussed in the previous section.
To more quantitatively verify the accuracy of PARMA4.0 in reproducing the EAS data, we calculated the coefficients of determination, R 2 , of the angular distributions for each particle and global condition using the following equation: where F x i ;E k ;EAS are the normalized angular distributions obtained from the EAS simulation for the i-th angular bin and the k-th energy bin, whose central values are x i and E k , respectively; F x i ;E k ;EAS is the mean value of F x i ;E k ;EAS ; and F(x i ,E k ) are the corresponding data calculated using Eqs (2) to (5). Fig 17 shows the calculated R 2 for r c = 1 and 15 GV as a function of atmospheric depth. Note that the R 2 values for He ions above 364 g/cm 2 , i.e., below 8 km, are not plotted because their angular distributions are assumed to be the same as those at 8 km, as previously mentioned. The calculated R 2 values above 50 g/cm 2 , i.e., below 21 km, are very close to 1.0, while they decrease at lower atmospheric depth particularly for the higher r c case. Two reasons are considered for this decreasing R 2 : the large statistical uncertainties in the EAS data and the difficulty in reproducing the peak structure around the horizontal direction by Eqs (2) to (5).  However, primary cosmic rays are the dominant contributor to the fluxes at high altitudes, and PARMA4.0 can reproduce the angular distributions of primary cosmic rays very well, as indicated by the high R 2 values for protons and He ions at r c = 1 GV. Based on these considerations, we concluded that PARMA4.0 allows for the instantaneous estimation of the differential cosmic-ray fluxes under most practical situations with an accuracy equivalent to that of the EAS simulation.

Comparison with Experimental Data
Figs 18-21 show a comparison between angular differential energy spectra of cosmic ray fluxes for various conditions obtained from PARMA4.0 and experiments [47][48][49][50][51][52][53][54][55][56][57][58]. For neutrons, measured angular differential fluxes are available only for very high altitudes or high energies, instead there are many experimental data of omnidirectional neutron fluxes. As shown in   Fig 19 shows that PARMA4.0 can reproduce the measured vertical proton fluxes very well, except for the high-energy data near sea level [50], which cannot be reproduced by EAS simulation neither as shown in our previous paper [17]. The measured spectrum at 105 g/cm 2 shown in panel (B) has a bump at around 3 GeV owing to the cut-off rigidity, a feature that is closely reproduced by the calculation.
As shown in Fig 20, excellent agreement can be observed between the measurements and calculations for muon fluxes. PARMA4.0 can reproduce detailed features of muon fluxes such as altitude dependence and charge ratio. Considering the agreements observed in muon fluxes for other angles and at higher energies, as shown in Fig 12, we can conclude that PARMA4.0 is well suited for the design of muon radiography.
It is seen from Fig 21 that the results obtained from PARMA4.0 agree with the corresponding experimental data fairly well, although the calculation overestimates and underestimates the high-energy electron fluxes and lower-altitude photon fluxes, respectively. The underestimation of photon fluxes may be owing to the ignorance of the primary cosmic rays above 1 TeV/n in our EAS simulation because such high-energy cosmic rays can generate a substantial number of photons at lower altitudes by inducing extensive air showers, although this hypothesis is rather contradictory to the overestimation of electron fluxes. More experimental data will therefore be necessary to validate the accuracy of PARMA4.0 in calculating electron, positron, and photon fluxes.

Conclusions
Based on the results of EAS simulation conducted using PHITS, we improved the PARMA model to version 4.0, which facilitates instantaneous estimation of not only omnidirectional but also angular differential energy spectra of cosmic ray fluxes anywhere in Earth's atmosphere at nearly any time. Both the presence of realistic ground and of the virtual black hole can be considered in the calculation of ground-level fluxes. The angular distributions of ground-level muons at larger zenith angles are specially determined using an optional function developed on the basis of experimental data. The accuracy of PARMA4.0 was carefully examined using multiple sets of experimental data obtained under various global conditions. The agreements between the calculated and measured data are generally satisfactory, particularly for muons, although some disagreements are observed in high-energy fluxes.
This extension of PARMA enlarges its applicability to wider research areas, such as the design of cosmic-ray detectors, muon radiography, soil moisture monitoring, and cosmic-ray shielding calculations. For example, it has already been used in the calculation of the attenuation length of underground muons [59]. PARMA4.0 is available freely and is easy to use, as implemented in the new version of EXPACS.